search

PoonLab / Kaphi

Kernel-embedded ABC-SMC for phylodynamic inference
GNU Affero General Public License v3.0
4 stars 2 forks source link

Malformed newick error when running `run.smc` with >100 particles #82

Closed MathiasRenaud closed 7 years ago

MathiasRenaud commented 7 years ago

Comments from issue #62:

When running smc with 1000 particles I ran into this error that @0ldM4j0r wrote for issue #70:

Initializing SMC-ABC run with 1000 particles
Step  1  epsilon: 0.3844203  ESS: 950.0131 accept: elapsed: 229 s
Step  2  epsilon: 0.377911  ESS: 902.4979 accept: 0.821 elapsed: 453.6 s
Step  3  epsilon: 0.3738658  ESS: 857.8386 accept: 0.7907445 elapsed: 676.2 s
Step  4  epsilon: 0.3714554  ESS: 814.7837 accept: 0.7212614 elapsed: 888 s
Step  5  epsilon: 0.3690808  ESS: 773.2843 accept: 0.7581152 elapsed: 1091.2 s
Step  6  epsilon: 0.3667808  ESS: 734.148 accept: 0.7195254 elapsed: 1288.8 s
Step  7  epsilon: 0.3650012  ESS: 698.2046 accept: 0.6920493 elapsed: 1479.5 s
Error in value[[3L]](cond) : Malformed Newick tree string!

I am going to mask his changes to see if an error appears again and there is an issue with the simulated trees or if there is an issue with the method used to check for malformed Newick strings.

Running run.smc with 1000 particles without the malformed Newick string check allows the function to run with no errors. However, the function aborts before stopping conditions are reached, giving poor results (run.smc stopped at: Step 12 epsilon: 0.3549191 ESS: 540.891 accept: 0.1227831 elapsed: 912 s when the config contains final.epsilon: 0.1 and final.accept.rate: 0.1).

A possibility is that the same error from issue #70 (Error in kids[[parent[i]]] : attempt to select less than one element in integerOneIndex) is still being encountered, and is stopping the function without printing the error message. If this is the case, it would mean that running the function with a large number of particles leads to an error in the simulation of new trees.

MathiasRenaud commented 7 years ago

I have not been able to reproduce the error message, but the function is still ending before stopping conditions are reached. These prematurely stopped runs result in poor posterior estimations: yule-1000-bad

MathiasRenaud commented 7 years ago

Today I've had many successful runs with varying smc settings and parameter distributions when the Malformed Newick tree string! check is turned off (function does not stop prematurely). When I run the package with the check turned on under the same settings, the error still occurs.

ArtPoon commented 7 years ago

Let's build in some exception handling to generate more informative output so that when this error happens again, we'll have more to work with.

Sent from my iPhone

On Jul 6, 2017, at 2:09 PM, Mathias Renaud notifications@github.com wrote:

Today I've had many successful runs with varying smc settings and parameter distributions when the Malformed Newick tree string! check is turned off (function does not stop prematurely). When I run the package with the check turned on under the same settings, the error still occurs.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

MathiasRenaud commented 7 years ago

On Wednesday when I reported that the function was stopping before stopping conditions were met, I think I misinterpreted the output because I had set verbose=FALSE and so the final conditions were not being printed last, but rather the 2nd last iteration. I modified the output for run.smc yesterday to provide more information and it seems that there are no errors with a large number of particles (when the malformed Newick check is off).

Therefore, the Malformed Newick tree string! error may be an issue with using plot() as the check with a large number of particles. However, I am not sure why this would be an issue since the function only handles one tree at a time.

MathiasRenaud commented 7 years ago

I'll try the methods detailed in this thread to obtain a more informative output

MathiasRenaud commented 7 years ago

Replacing: 

stop("Malformed Newick tree string!")

with: 

print("It's that .to.newick error!")
traceback()

resulted in the following output during the run with apparently 12 errors in .to.newick:

Step  26  epsilon: 0.2933377  ESS: 264.5026 accept: 0.4426606 elapsed: 2544.6 s
run.smc niter:  26 
ws$epsilon:  0.2933377 
config$final.epsilon:  0.05 
result$accept.rate:  0.7963892 0.7827411 0.7494781 0.7567568 0.7293065 0.6876457 0.7028015 0.7135615 0.6859396 0.6607387 0.6571835 0.6366322 0.6067588 0.616242 0.5860927 0.5851979 0.5673759 0.5659341 0.5776515 0.524558 0.4898374 0.4968553 0.4924078 0.4146341 0.4426606 0.3791469 
config$final.accept.rate:  0.05 
Step  27  epsilon: 0.2901515  ESS: 250.7329 accept: 0.3791469 elapsed: 2598.3 s
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
[1] "It's that .to.newick error!"
No traceback available 
run.smc niter:  27 
ws$epsilon:  0.2901515 
config$final.epsilon:  0.05 
result$accept.rate:  0.7963892 0.7827411 0.7494781 0.7567568 0.7293065 0.6876457 0.7028015 0.7135615 0.6859396 0.6607387 0.6571835 0.6366322 0.6067588 0.616242 0.5860927 0.5851979 0.5673759 0.5659341 0.5776515 0.524558 0.4898374 0.4968553 0.4924078 0.4146341 0.4426606 0.3791469 0.4156479 
config$final.accept.rate:  0.05 
Step  28  epsilon: 0.286117  ESS: 238.1637 accept: 0.4156479 elapsed: 2649.3 s

After the error at iteration 27, the function continued to run without any problems (71 iterations).

I'm going to see if I can get more information from the error by replacing stop("Malformed Newick tree string!") with: 

e$message <- paste(".to.newick error:", e, sep=' ')
        stop(e)
MathiasRenaud commented 7 years ago

Finally got an actual error message!

Step  27  epsilon: 0.2901515  ESS: 250.7329 accept: 0.3791469 elapsed: 2532.6 s
Error in plot.window(...) : 
  .to.newick error: Error in plot.window(...): need finite 'xlim' values

This appears to be an issue with using plot as the function that checks the quality of the Newick string and not with run.smc. For some reason, the simulated trees passed to .to.newick in step 27 with my current settings don't result in the same xlim values for the plot window as the other steps.

After looking up the issue on stack overflow, it seems that it may be due to a value of NA in the xlim vector.

ArtPoon commented 7 years ago

Looks like a stochastic error. We need to get .to.newick to return the Newick tree string that it is calling malformed in order to diagnose this problem.

MathiasRenaud commented 7 years ago

I added the line print(write.tree(tree)) to .to.newick and the function printed the Newick string but it was truncated (below). I'm running it again but saving the Newick string to a variable instead to store the whole string.

"(((((((((((sp42:2167.869606,sp43:2167.869606)nd46:14987.01557,sp23:17154.88517)nd45:144.1673116,sp22:17299.05248)nd31:43051.32886,((sp32:9772.393965,sp33:9772.393965)nd34:50568.13258,sp5:60340.52655)nd30:9.854800695)nd25:4404.90409,sp4:64755.28544)nd20:12667.60906,sp3:77422.89449)nd11:75888.78962,((((sp18:18926.20287,sp19:18926.20287)nd43:7753.035868,sp14:26679.23874)nd41:2252.89409,(sp34:7914.73943,sp35:7914.73943)nd40:21017.3934)nd33:31415.00277,(sp16:22786.27606,sp17:22786.27606)nd32:37560.85954)nd10:92964.54851)nd8:13247.50044,((((sp20:17547.26825,sp21:17547.26825)nd35:40771.56469,sp6:58318.83295)nd28:2936.722724,((sp30:9959.671735,sp31:9959.671735)nd36:44138.63447,sp7:54098.3062)nd29:7157.249468)nd14:36736.45445,sp1:97992.01012)nd9:68567.17444)nd6:18776.63587,(sp8:40956.23079,sp9:40956.23079)nd7:144379.5896)nd4:32051.59379,((((sp49:479.6748893,sp50:479.6748893)nd42:26644.45952,sp13:27124.13441)nd23:44264.81769,((sp47:573.7714213,sp48:573.7714213)nd49:1486.435529,sp46:2060.206... <truncated>
MathiasRenaud commented 7 years ago

Full Newick strong causing an error:

"(((((((((((sp42:2167.869606,sp43:2167.869606)nd46:14987.01557,sp23:17154.88517)nd45:144.1673116,sp22:17299.05248)nd31:43051.32886,((sp32:9772.393965,sp33:9772.393965)nd34:50568.13258,sp5:60340.52655)nd30:9.854800695)nd25:4404.90409,sp4:64755.28544)nd20:12667.60906,sp3:77422.89449)nd11:75888.78962,((((sp18:18926.20287,sp19:18926.20287)nd43:7753.035868,sp14:26679.23874)nd41:2252.89409,(sp34:7914.73943,sp35:7914.73943)nd40:21017.3934)nd33:31415.00277,(sp16:22786.27606,sp17:22786.27606)nd32:37560.85954)nd10:92964.54851)nd8:13247.50044,((((sp20:17547.26825,sp21:17547.26825)nd35:40771.56469,sp6:58318.83295)nd28:2936.722724,((sp30:9959.671735,sp31:9959.671735)nd36:44138.63447,sp7:54098.3062)nd29:7157.249468)nd14:36736.45445,sp1:97992.01012)nd9:68567.17444)nd6:18776.63587,(sp8:40956.23079,sp9:40956.23079)nd7:144379.5896)nd4:32051.59379,((((sp49:479.6748893,sp50:479.6748893)nd42:26644.45952,sp13:27124.13441)nd23:44264.81769,((sp47:573.7714213,sp48:573.7714213)nd49:1486.435529,sp46:2060.20695)nd24:69328.74514)nd13:75106.94349,((sp36:7890.297862,sp37:7890.297862)nd21:67720.43795,(sp26:12347.78144,sp27:12347.78144)nd22:63262.95437)nd12:70885.15977)nd5:70891.51864)nd3:18015.58198,((((sp28:12205.9334,sp29:12205.9334)nd26:49749.36798,(sp38:3939.406508,sp39:3939.406508)nd27:58015.89488)nd18:15596.08148,(((sp40:2584.535187,sp41:2584.535187)nd48:9952.076814,sp25:12536.612)nd47:3953.680507,sp24:16490.29251)nd19:61061.09036)nd15:3537.641053,(((((sp44:2136.433198,sp45:2136.433198)nd44:20805.6963,sp15:22942.1295)nd39:11788.19381,sp10:34730.32331)nd38:930.2612248,(sp11:30581.12224,sp12:30581.12224)nd37:5079.462298)nd17:45034.56675,sp2:80695.15128)nd16:393.8726391)nd2:154313.9723)nd1;"
ArtPoon commented 7 years ago

Thanks for grabbing that. Strange I'm not able to reproduce an error with this Newick string though.

MathiasRenaud commented 7 years ago

I picked apart the string to check for anything irregular and I couldn't find any issues with it. I still got errors in R though:


err.tree <- read.tree(file='Documents/error-tree.txt')
> err.tree

Phylogenetic tree with 50 tips and 49 internal nodes.

Tip labels:
    sp42, sp43, sp23, sp22, sp32, sp33, ...
Node labels:
    nd1, nd3, nd4, nd6, nd8, nd11, ...

Rooted; includes branch lengths.
> plot(err.tree)
Error in plot.window(...) : need finite 'xlim' values
ArtPoon commented 7 years ago

A-ha. Branch lengths are too large. Try:

tr <- read.tree(text="(((((((((((sp42:2167.869606,sp43:2167.869606)nd46:14987.01557,sp23:17154.88517)nd45:144.1673116,sp22:17299.05248)nd31:43051.32886,((sp32:9772.393965,sp33:9772.393965)nd34:50568.13258,sp5:60340.52655)nd30:9.854800695)nd25:4404.90409,sp4:64755.28544)nd20:12667.60906,sp3:77422.89449)nd11:75888.78962,((((sp18:18926.20287,sp19:18926.20287)nd43:7753.035868,sp14:26679.23874)nd41:2252.89409,(sp34:7914.73943,sp35:7914.73943)nd40:21017.3934)nd33:31415.00277,(sp16:22786.27606,sp17:22786.27606)nd32:37560.85954)nd10:92964.54851)nd8:13247.50044,((((sp20:17547.26825,sp21:17547.26825)nd35:40771.56469,sp6:58318.83295)nd28:2936.722724,((sp30:9959.671735,sp31:9959.671735)nd36:44138.63447,sp7:54098.3062)nd29:7157.249468)nd14:36736.45445,sp1:97992.01012)nd9:68567.17444)nd6:18776.63587,(sp8:40956.23079,sp9:40956.23079)nd7:144379.5896)nd4:32051.59379,((((sp49:479.6748893,sp50:479.6748893)nd42:26644.45952,sp13:27124.13441)nd23:44264.81769,((sp47:573.7714213,sp48:573.7714213)nd49:1486.435529,sp46:2060.20695)nd24:69328.74514)nd13:75106.94349,((sp36:7890.297862,sp37:7890.297862)nd21:67720.43795,(sp26:12347.78144,sp27:12347.78144)nd22:63262.95437)nd12:70885.15977)nd5:70891.51864)nd3:18015.58198,((((sp28:12205.9334,sp29:12205.9334)nd26:49749.36798,(sp38:3939.406508,sp39:3939.406508)nd27:58015.89488)nd18:15596.08148,(((sp40:2584.535187,sp41:2584.535187)nd48:9952.076814,sp25:12536.612)nd47:3953.680507,sp24:16490.29251)nd19:61061.09036)nd15:3537.641053,(((((sp44:2136.433198,sp45:2136.433198)nd44:20805.6963,sp15:22942.1295)nd39:11788.19381,sp10:34730.32331)nd38:930.2612248,(sp11:30581.12224,sp12:30581.12224)nd37:5079.462298)nd17:45034.56675,sp2:80695.15128)nd16:393.8726391)nd2:154313.9723)nd1;")
tr$edge.length <- tr$edge.length / 10000
plot(tr)
MathiasRenaud commented 7 years ago

Would the solution then be to use an alternative to plot for checking the trees?

ArtPoon commented 7 years ago

Yeah - plot was a bad idea. Wasn't the original problem that .to.newick() was failing to catch malformed strings? We need some other method...

MathiasRenaud commented 7 years ago

I reopened issue #70 to brainstorm new ideas to replace plot.