traitecoevo / plant

Trait-Driven Models of Ecology and Evolution :evergreen_tree:
https://traitecoevo.github.io/plant
53 stars 20 forks source link

Find equilibrium seed rain using root finding/minimisation #85

Closed richfitz closed 10 years ago

richfitz commented 10 years ago

Should be possible to do much better than our naive iteration approach.

dfalster commented 10 years ago

Indeed

On 14 August 2014 11:56, Rich FitzJohn notifications@github.com wrote:

Should be possible to do much better than our naive iteration approach.

— Reply to this email directly or view it on GitHub https://github.com/richfitz/tree/issues/85.

richfitz commented 10 years ago

Relevant chain of discussion here: https://github.com/richfitz/tree/commit/3c2090ff7f59

richfitz commented 10 years ago

@dfalster: can you keep an eye out for a case where we're still a reasonable distance from equilibrium after 20 steps? There should be one in the logs. All we need is the parameter set, traits and seed rains. Then I can poke about with that when tuning this.

dfalster commented 10 years ago

Running this in 'copy' folder with all our results reveals a bunch:

ack --literal "*** 20:"

Seems there three possible causes

1. Churning, asking for too much accuracy (on large seed rain)

Experiment: seed_size_initial_exploration Task: simulations Id: 2 StdOut stored in seed_size_initial_exploration/jobs/seed_size_initi.e6624580 Gets stuck on first attempt:

Calculating selection gradient for s = 0.000001
1: Splitting {4} times (97)
2: Splitting {1} times (101)
*** 1: {1} -> {1305.472} (delta = {1304.472})
1: Splitting {29} times (97)
2: Splitting {1} times (126)
*** 2: {1305.472} -> {15997.8} (delta = {14692.33})
1: Splitting {32} times (97)
2: Splitting {9} times (129)
*** 3: {15997.8} -> {22489.98} (delta = {6492.177})
1: Splitting {41} times (97)
2: Splitting {9} times (138)
*** 4: {22489.98} -> {23346.82} (delta = {856.846})
1: Splitting {37} times (97)
2: Splitting {9} times (134)
3: Splitting {1} times (143)
*** 5: {23346.82} -> {23440.91} (delta = {94.0839})
1: Splitting {36} times (97)
2: Splitting {8} times (133)
3: Splitting {1} times (141)
*** 6: {23440.91} -> {23449.95} (delta = {9.047037})
1: Splitting {1} times (142)
*** 7: {23449.95} -> {23450.99} (delta = {1.031435})
*** 8: {23450.99} -> {23451.05} (delta = {0.0676576})
*** 9: {23451.05} -> {23451.1} (delta = {0.04160394})
*** 10: {23451.1} -> {23451.07} (delta = {-0.0288013})
*** 11: {23451.07} -> {23451.06} (delta = {-0.002899799})
*** 12: {23451.06} -> {23451.08} (delta = {0.01412029})
*** 13: {23451.08} -> {23451.06} (delta = {-0.01299117})
*** 14: {23451.06} -> {23451.08} (delta = {0.01306854})
*** 15: {23451.08} -> {23451.06} (delta = {-0.0130615})
*** 16: {23451.06} -> {23451.08} (delta = {0.01305788})
*** 17: {23451.08} -> {23451.06} (delta = {-0.01305564})
*** 18: {23451.06} -> {23451.08} (delta = {0.01305349})
*** 19: {23451.08} -> {23451.06} (delta = {-0.01305483})
*** 20: {23451.06} -> {23451.08} (delta = {0.0130727})

Seems we're asking for bit much in terms of accuracy for such a high seed_rain. Had stabilised by step 11, after that just churning.

2. Slow approach to equilibrium

Experiment: seed_size_initial_exploration Task: simulations Id: 2 StdOut stored in seed_size_initial_exploration/jobs/seed_size_initi.e6624580

Calculating selection gradient for s = 0.026801
1: Splitting {6} times (97)
*** 1: {1} -> {0.3509203} (delta = {-0.6490797})
*** 2: {0.3509203} -> {0.200943} (delta = {-0.1499773})
1: Splitting {2} times (103)
*** 3: {0.200943} -> {0.1387551} (delta = {-0.06218786})
*** 4: {0.1387551} -> {0.1055436} (delta = {-0.03321148})
*** 5: {0.1055436} -> {0.08509937} (delta = {-0.02044426})
*** 6: {0.08509937} -> {0.07132089} (delta = {-0.01377848})
*** 7: {0.07132089} -> {0.0614354} (delta = {-0.009885487})
*** 8: {0.0614354} -> {0.05401114} (delta = {-0.007424264})
*** 9: {0.05401114} -> {0.0482375} (delta = {-0.005773641})
*** 10: {0.0482375} -> {0.04362264} (delta = {-0.004614864})
*** 11: {0.04362264} -> {0.03985154} (delta = {-0.003771092})
*** 12: {0.03985154} -> {0.03671356} (delta = {-0.003137985})
*** 13: {0.03671356} -> {0.03406245} (delta = {-0.00265111})
*** 14: {0.03406245} -> {0.0317936} (delta = {-0.002268852})
*** 15: {0.0317936} -> {0.02983022} (delta = {-0.001963379})
*** 16: {0.02983022} -> {0.02811482} (delta = {-0.001715395})
*** 17: {0.02811482} -> {0.02660342} (delta = {-0.001511401})
*** 18: {0.02660342} -> {0.02526183} (delta = {-0.001341589})
*** 19: {0.02526183} -> {0.02406311} (delta = {-0.001198728})
*** 20: {0.02406311} -> {0.02298568} (delta = {-0.00107743})

In this case approach to equilibrium is very slow (large seed size) so takes a long time to get there.

Another in Experiment: wood_density_initial_exploration Task: mortality_height Id: 1 StdOut stored in wood_density_initial_exploration/jobs/wood_density_in.e6626174

Calculating selection gradient for rho = 4.954626
1: Splitting {13} times (92)
2: Splitting {15} times (105)
3: Splitting {9} times (120)
4: Splitting {1} times (129)
5: Splitting {1} times (130)
6: Splitting {1} times (131)
*** 1: {1} -> {0.9955209} (delta = {-0.004479118})
*** 2: {0.9955209} -> {0.9910676} (delta = {-0.004453243})
*** 3: {0.9910676} -> {0.9866401} (delta = {-0.004427569})
*** 4: {0.9866401} -> {0.982238} (delta = {-0.004402095})
*** 5: {0.982238} -> {0.9778612} (delta = {-0.004376817})
*** 6: {0.9778612} -> {0.9735094} (delta = {-0.004351735})
*** 7: {0.9735094} -> {0.9691826} (delta = {-0.004326846})
*** 8: {0.9691826} -> {0.9648804} (delta = {-0.004302149})
*** 9: {0.9648804} -> {0.9606028} (delta = {-0.00427764})
*** 10: {0.9606028} -> {0.9563495} (delta = {-0.004253319})
*** 11: {0.9563495} -> {0.9521203} (delta = {-0.004229184})
*** 12: {0.9521203} -> {0.9479151} (delta = {-0.004205232})
*** 13: {0.9479151} -> {0.9437336} (delta = {-0.004181462})
*** 14: {0.9437336} -> {0.9395757} (delta = {-0.004157873})
*** 15: {0.9395757} -> {0.9354413} (delta = {-0.004134461})
*** 16: {0.9354413} -> {0.93133} (delta = {-0.004111226})
*** 17: {0.93133} -> {0.9272419} (delta = {-0.004088165})
*** 18: {0.9272419} -> {0.9231766} (delta = {-0.004065278})
*** 19: {0.9231766} -> {0.919134} (delta = {-0.004042562})
*** 20: {0.919134} -> {0.915114} (delta = {-0.004020015})

Again, looks like slow approach to equilibrium

3. Lots of types

When there are large number of types, it's understandable they take a while to stabilise. Could be because of 1 and 2 above but hard to determine. for examples see

Experiment: compare_assembly_algorithms Task: simulations Id: 8 StdOut stored in: compare_assembly_algorithms/jobs/compare_assembl.e6616322

*** Assembler: step 15, (6 strategies), to_equilibrium
Recomputing ode times
1: Splitting {1,3,28,1,2,0} times (106,106,106,106,106,106)
2: Splitting {0,5,9,0,4,0} times (107,109,134,107,108,106)
3: Splitting {0,2,9,1,0,0} times (107,114,143,107,112,106)
4: Splitting {0,1,3,0,1,0} times (107,116,152,108,112,106)
*** 1: {28.73691,93.44078,305.4931,7.079078,23.23914,0.01094953} ->
{24.80488,90.71826,311.0706,7.231333,24.21034,0.01178602} (delta =
{-3.932029,-2.722516,5.577536,0.1522548,0.9712002,0.0008364944})
*** 2: {24.80488,90.71826,311.0706,7.231333,24.21034,0.01178602} ->
{21.3822,87.92618,316.2015,7.374076,25.17817,0.01266425} (delta =
{-3.422686,-2.792085,5.130852,0.142743,0.9678243,0.000878229})
*** 3: {21.3822,87.92618,316.2015,7.374076,25.17817,0.01266425} ->
{18.41026,85.09339,320.9278,7.50819,26.14456,0.01358706} (delta =
{-2.971941,-2.832788,4.726379,0.134114,0.9663923,0.0009228114})
*** 4: {18.41026,85.09339,320.9278,7.50819,26.14456,0.01358706} ->
{15.83501,82.24315,325.2859,7.634383,27.11104,0.01455722} (delta =
{-2.575248,-2.850243,4.358116,0.1261928,0.9664842,0.0009701569})
*** 5: {15.83501,82.24315,325.2859,7.634383,27.11104,0.01455722} ->
{13.60723,79.39311,329.2975,7.753144,28.07839,0.01557727} (delta =
{-2.227776,-2.850042,4.011602,0.1187606,0.9673501,0.001020052})
*** 6: {13.60723,79.39311,329.2975,7.753144,28.07839,0.01557727} ->
{11.68312,76.55958,332.9907,7.865071,29.04795,0.01665026} (delta =
{-1.924117,-2.83353,3.693181,0.1119271,0.9695561,0.001072989})
*** 7: {11.68312,76.55958,332.9907,7.865071,29.04795,0.01665026} ->
{10.02356,73.75581,336.3903,7.970681,30.02089,0.01777933} (delta =
{-1.659555,-2.803761,3.399527,0.10561,0.9729429,0.001129067})
*** 8: {10.02356,73.75581,336.3903,7.970681,30.02089,0.01777933} ->
{8.593975,70.99284,339.5194,8.070446,30.99828,0.01896766} (delta =
{-1.429587,-2.762973,3.129126,0.09976526,0.9773865,0.001188334})
*** 9: {8.593975,70.99284,339.5194,8.070446,30.99828,0.01896766} ->
{7.363807,68.27954,342.3991,8.164823,31.98128,0.02021885} (delta =
{-1.230167,-2.713302,2.87975,0.09437726,0.9829983,0.001251188})
*** 10: {7.363807,68.27954,342.3991,8.164823,31.98128,0.02021885} ->
{6.306253,65.62265,345.0458,8.25415,32.97062,0.02153627} (delta =
{-1.057554,-2.656891,2.646701,0.08932625,0.9893484,0.001317426})
*** 11: {6.306253,65.62265,345.0458,8.25415,32.97062,0.02153627} ->
{5.39788,63.0279,347.4778,8.338808,33.96739,0.02292386} (delta =
{-0.9083732,-2.594748,2.431943,0.08465886,0.9967634,0.001387587})
1: Splitting {0,0,1,0,0,0} times (107,117,155,108,113,106)
*** 12: {5.39788,63.0279,347.4778,8.338808,33.96739,0.02292386} ->
{4.618844,60.5051,349.6962,8.419654,34.97413,0.02438583} (delta =
{-0.7790362,-2.522797,2.218372,0.08084534,1.006747,0.001461971})
*** 13: {4.618844,60.5051,349.6962,8.419654,34.97413,0.02438583} ->
{3.950654,58.05228,351.7335,8.49657,35.99056,0.02592647} (delta =
{-0.6681897,-2.452823,2.037333,0.07691643,1.016423,0.001540643})
*** 14: {3.950654,58.05228,351.7335,8.49657,35.99056,0.02592647} ->
{3.377826,55.67048,353.5948,8.569642,37.01668,0.02754962} (delta =
{-0.5728288,-2.381797,1.861319,0.07307225,1.026118,0.001623142})
*** 15: {3.377826,55.67048,353.5948,8.569642,37.01668,0.02754962} ->
{2.88706,53.36181,355.2952,8.639195,38.05356,0.02926009} (delta =
{-0.4907658,-2.308677,1.700399,0.06955244,1.036889,0.001710469})
*** 16: {2.88706,53.36181,355.2952,8.639195,38.05356,0.02926009} ->
{2.466796,51.12672,356.8428,8.705387,39.10169,0.03106258} (delta =
{-0.420264,-2.235089,1.547643,0.06619225,1.048127,0.001802495})
1: Splitting {0,0,1,0,1,0} times (107,117,156,108,113,106)
*** 17: {2.466796,51.12672,356.8428,8.705387,39.10169,0.03106258} ->
{2.106959,48.96513,358.248,8.768699,40.16425,0.03296597} (delta =
{-0.3598363,-2.161584,1.405163,0.06331195,1.062557,0.00190339})
*** 18: {2.106959,48.96513,358.248,8.768699,40.16425,0.03296597} ->
{1.799101,46.87706,359.5151,8.828947,41.23902,0.03497173} (delta =
{-0.3078579,-2.088069,1.267131,0.0602481,1.074767,0.002005759})
*** 19: {1.799101,46.87706,359.5151,8.828947,41.23902,0.03497173} ->
{1.53582,44.86202,360.6521,8.88629,42.32654,0.03708544} (delta =
{-0.2632812,-2.015041,1.136914,0.05734233,1.087523,0.002113713})
*** 20: {1.53582,44.86202,360.6521,8.88629,42.32654,0.03708544} ->
{1.310787,42.92018,361.6699,8.940949,43.4275,0.03931223} (delta =
{-0.2250334,-1.941847,1.017849,0.05465909,1.100957,0.002226783})

On 11 September 2014 14:49, Rich FitzJohn notifications@github.com wrote:

@dfalster https://github.com/dfalster: can you keep an eye out for a case where we're still a reasonable distance from equilibrium after 20 steps? There should be one in the logs. All we need is the parameter set, traits and seed rains. Then I can poke about with that when tuning this.

— Reply to this email directly or view it on GitHub https://github.com/richfitz/tree/issues/85#issuecomment-55219071.

richfitz commented 10 years ago

Just so I know - how did you move the jobs onto the Copy volume? Just grabbed all the files manually with scp/rsync?

dfalster commented 10 years ago

from within Copy/../successional_diversity/

I ran

rsync -avz dsf576@raijin.nci.org.au:/home/576/dsf576/successional_diversity/analysis/experiments/output/
.

NB. In an attempt to "clean" unwanted files I deleted some of the *.rds.1 backups, but just remembered that some of the logs showed having loaded from that, so I may have browk the dependency chain for those.

On 12 September 2014 09:30, Rich FitzJohn notifications@github.com wrote:

Just so I know - how did you move the jobs onto the Copy volume? Just grabbed all the files manually with scp/rsync?

— Reply to this email directly or view it on GitHub https://github.com/richfitz/tree/issues/85#issuecomment-55343487.

richfitz commented 10 years ago

First case can be dealt with by implementing relative accuracy as well as absolute.

richfitz commented 10 years ago

Multidimensional secant (Broyden's method) might come in useful. Alternatively, hone in close and do a number of essentially independent root finding steps, but I suspect that's basically what Broyden's method is doing (with the initial stab at the Jacobian taking a long time).

Another solution would be to formulate as a set of ODEs and use runsteady etc to find the stable point. One of the solutions there does nonlinear equation solving nleqslv I think.

Level set method might be better longer term.

richfitz commented 10 years ago

nleqslv looks like the place to start. Code to drive this already in Revolve (repurpose that). Test that it works when initial guess is terrible.

richfitz commented 10 years ago

This one turns out to be a pain:

Working with per-capita fitness (log(out/in)) is really nice when there is a positive root, but terrible when there is not.

richfitz commented 10 years ago

See script here

richfitz commented 10 years ago

Reopening this, in light of the continued dramas. I think that starting off with a few rounds of iteration might be the way to do this. Seems like a hack, but the solver is really badly behaved around zero.

richfitz commented 10 years ago

I think that runsteady might finally work here. Once the last batch of simulations are back, if they've worked let's close this. Then the issue just becomes "can we get there faster"

richfitz commented 10 years ago

So runsteady was a dud: now trying an iterate & solve & repeat approach. Ironing out the last wrinkles in that at the moment.