Problem simulating tree from RP model

ArtPoon commented 7 years ago

Error in sample.path(a.state, b.state, 0, z$edge.length[u], qm) : 
  The process cannot start in an absorbing state if the endpoints are different.

The corresponding line in sample.path() is:

 if (all(Q[a, ] == 0) & a != b) {
        stop("The process cannot start in an absorbing state if the endpoints are different.", Q)
    }

ArtPoon commented 7 years ago

Having trouble calculating the A matrix, current resolution (1e5) is inadequate and entries are plunging into negative values. As a temporary fix, I am zeroing out negative entries with the following line:

A.mx$mat[A.mx$mat<0] <- 0

ArtPoon commented 7 years ago

Although fixing the A matrix computation was an improvement, it didn't not address this problem :-/

ArtPoon commented 7 years ago

Turned off simulate.migrations option, now a different error is coming up:

Error in sample.int(z$m^2, size = 1, prob = as.vector(.lambdamat)) : 
  NA in probability vector

ArtPoon commented 7 years ago

Switching integration method from adams to rk4 in the A matrix function may have been part of the problem; the latter was giving nonsensical results although it complained less. Switched back and doing some more tests.

ArtPoon commented 7 years ago

Error in kids[[parent[i]]] : 
  attempt to select more than one element in get1index <real>

Tracing through the workflow, the z$edge matrix is not getting populated completely:

> z$edge
       [,1] [,2]
  [1,]  211    1
  [2,]  219    2
  [3,]  221    3
...
[188,]  383  188
[189,]   -1   -1
[190,]   -1   -1
...
[209,]   -1   -1
[210,]  213  210
...
[393,]  394  393
[394,]   -1   -1
[395,]   -1   -1
...
[416,]   -1   -1

ArtPoon commented 7 years ago

Ok, tips 189:209 are missing edges. If we look at the sorted sample heights, we get a clue as to why:

> sorted.sample.heights[188:209]
 185  180  181  184  186  187  188  189  190  191  193  195  199  201  202  203 
2768 2981 2981 2981 2981 2981 2981 2981 2981 2981 2981 2981 2981 2981 2981 2981 
 204  205  206  207  208  209 
2981 2981 2981 2981 2981 2981

ArtPoon commented 7 years ago

Output of some debugging statements:

ih h0 h1 n.extant event.type
196 2782.686 2787.183 3 coalesce
197 2787.183 2832.373 2 coalesce
198 2832.373 2836.819 1 coalesce
199 2836.819 2838.879 1 coalesce
200 2838.879 2840.724 1 coalesce
201 2840.724 2841.049 1 coalesce
202 2841.049 2842.086 1 coalesce
203 2842.086 2842.456 1 coalesce
204 2842.456 2842.459 1 coalesce
205 2842.459 2842.543 1 coalesce
206 2842.543 2843.995 1 coalesce
207 2843.995 2848.287 1 coalesce
208 2848.287 2889.107 1 coalesce
209 2889.107 2890.175 1 coalesce
210 2890.175 2914.357 1 coalesce
211 2914.357 2914.602 1 coalesce
212 2914.602 2922.509 1 coalesce
213 2922.509 2929.117 1 coalesce
214 2929.117 2932.623 1 coalesce
215 2932.623 2938.258 1 coalesce
216 2938.258 2952.577 1 coalesce
217 2952.577 2958.143 1 coalesce
218 2958.143 2961.091 1 coalesce
219 2961.091 2972.026 1 coalesce
220 2972.026 2981 22 sample

What concerns me here is that we have a bunch of coalescent events at a single lineage because the remaining 21 lineages haven't been sampled until time 2981.

haxis runs to 2990.70090
`Et$event.times$ terminates at 2981.

A.mx$mat[9900:10000,] has an abrupt transition:

[65,] 11.000000728 8.000002 2.0000000
[66,] 11.000000024 8.000002 2.0000000
[67,] 11.000000044 8.000002 2.0000000
[68,]  0.602944749 8.138382 1.8128595
[69,]  0.013494601 8.400925 1.4559388
[70,]  0.010159910 8.601264 1.1813309

This row corresponds to haxis index 9967 with height 2980.83.

ArtPoon commented 7 years ago

My reimplementation of the A matrix computation (odA) still diverges from Erik's (odA2):

> tail(odA)
            time             V             L             Ts
[9736,] 2912.030  1.100000e+01  9.000000e+00   1.000000e+00
[9737,] 2912.329  1.100000e+01  9.000000e+00   1.000000e+00
[9738,] 2912.628  1.100000e+01  9.000000e+00   1.000000e+00
[9739,] 2912.927  1.100000e+01  9.000000e+00   1.000000e+00
[9740,] 2913.226  1.100000e+01  9.000000e+00   1.000000e+00
[9741,] 2913.272 -2.130763e+95 -1.941097e+93 -4.157117e+201
> tail(odA2)
             time           V        L         Ts
 [9995,] 2989.504 0.004125392 9.411804 0.08856393
 [9996,] 2989.803 0.004292994 9.384265 0.08778073
 [9997,] 2990.103 0.004231863 9.356954 0.08718515
 [9998,] 2990.402 0.003923977 9.331999 0.08676711
 [9999,] 2990.701 0.209648932 9.077899 0.08382280
[10000,] 2990.701 0.209664931 9.077874 0.08382251

ArtPoon commented 7 years ago

Switching back to interpolation scheme seems to help. Leaving space for coalescence when tips have been sampled at time zero is also important.

ArtPoon commented 7 years ago

Closing - there are still problems but it will be more informative to write them up as specific issues.

PoonLab / Kaphi

Problem simulating tree from RP model #33