popgenmethods / momi2

Infer demographic history with the Moran model
GNU General Public License v3.0
47 stars 11 forks source link

momi.events.DemographyError #59

Open pengyan19 opened 1 year ago

pengyan19 commented 1 year ago

Hello, when I used the momi2, I meet a question as follows:

This the code: logging.basicConfig(level=logging.INFO,filename="momi2.log") pop1=sys.argv[1] pop2=sys.argv[2] sfs = momi.Sfs.load("_".join([pop1,pop2,"allele_sfs.gz"]))

Adding migration events

migration_model = momi.DemographicModel(N_e=7e5,gen_time=0.25, muts_per_gen=0.84e-8)

migration_model.set_data(sfs,length=470000000)

migration_model.add_leaf(pop1, N="n_north") migration_model.add_leaf(pop2, N="n_south")

migration_model.move_lineages(pop2, pop1, t="tdiv",N="N_pop1_pop2") migration_model.move_lineages(pop1, pop2, t="tmig_north_south", p="mfrac_north_south") migration_model.move_lineages(pop2, pop1, t="tmig_south_north", p="mfrac_south_north")

migration_model.add_time_param("tmig_north_south",10000) migration_model.add_time_param("tmig_south_north",10000) migration_model.add_time_param("tdiv", lower_constraints=["tmig_north_south", "tmig_south_north"]) migration_model.add_size_param("n_north",10000) migration_model.add_size_param("n_south",10000) migration_model.add_size_param("N_pop1_pop2",10000) migration_model.add_pulse_param("mfrac_north_south", upper=1) migration_model.add_pulse_param("mfrac_south_north", upper=1)

results = [] migration_model_copy = migration_model.copy() migration_model_copy.set_params(migration_model.get_params(),randomize=True) results.append(migration_model_copy.optimize(method="TNC",options={"maxiter":5000, "ftol":0.84e-8}))

migration_model_copy=migration_model.copy() n_bootstraps = 100 bootstrap_results = [] model_list=[] for i in range(n_bootstraps): print(f"migration Fitting {i+1}-th bootstrap out of {n_bootstraps}") resampled_sfs = sfs.resample() migration_model_copy.set_data(resampled_sfs,length=470000000) migration_model_copy.set_params(randomize=True) migration_model_copy.optimize() model_list.append(migration_model_copy) bootstrap_results.append(migration_model_copy.get_params()) AICs_4 = [] for model in model_list: lik = model.log_likelihood() nparams = len(model.get_params()) aic = 2nparams - 2lik print("AIC {}".format(aic)) AICs_4.append(aic) print("AIC---bootstraps") for i,j in zip(bootstrap_results,AICs_4):

print(i,j)

print("---done---")
fig = momi.DemographyPlot(migration_model_copy,[pop1, pop2],
                      linthreshy=5000, figsize=(6,8),major_yticks=yticks)
for params in bootstrap_results:
    print(params)
    fig.add_bootstrap(params,alpha=1/n_bootstraps)

import matplotlib.backends.backend_pdf pdf=matplotlib.backends.backend_pdf.PdfPages("".join([pop1,".",pop2,".","migration.bootstraps_model.pdf"])) pdf.savefig(fig.draw()) pdf.close()

ERROR: Traceback (most recent call last): File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 557, in _demo_fun return self._get_demo(None) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 540, in _get_demo G = _build_demo_graph(events, sampled_n_dict, params_dict, default_N=1.0) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 21, in _build_demo_graph e.add_to_graph(_G, sample_sizes, params_dict) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 270, in add_to_graph _check_ej_ep_pops(G, '-ep', t, i, j, pij) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 368, in _check_ej_ep_pops raise DemographyError( momi.events.DemographyError: Invalid move_lineages event at time Autograd ArrayBox with value 0.031590331937878385: pop North was already removed by previous move_lineages

The above exception was the direct cause of the following exception:

Traceback (most recent call last): File "/public/home/pengyan/bioproject/01.ACB/24.momi2/06.momi-North-South/momi2_mig.py", line 50, in migration_model_copy.optimize() File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 910, in optimize res = self._get_surface().find_mle( File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 267, in find_mle return _find_minimum(fun, x0, scipy.optimize.minimize, File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/optimizers.py", line 88, in _find_minimum ret = _find_minimum_helper( File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/optimizers.py", line 103, in _find_minimum_helper return optimizer(f, start_params, opt_kwargs) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_minimize.py", line 684, in minimize res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback, File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_tnc.py", line 416, in _minimize_tnc rc, nf, nit, x, funv, jacv = moduleTNC.tnc_minimize( File "_moduleTNC.pyx", line 201, in _moduleTNC.tnc_minimize File "_moduleTNC.pyx", line 73, in _moduleTNC.function File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 285, in fun_and_grad self._update_fun() File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 251, in _update_fun self._update_fun_impl() File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 155, in update_fun self.f = fun_wrapped(self.x) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 137, in fun_wrapped fx = fun(np.copy(x), args) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 76, in call self._compute_if_needed(x, args) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 70, in _compute_if_needed fg = self.fun(x, args) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/wrap_util.py", line 20, in nary_f return unary_operator(unary_f, x, nary_op_args, nary_op_kwargs) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/differential_operators.py", line 135, in value_and_grad vjp, ans = _make_vjp(fun, x) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/core.py", line 10, in make_vjp end_value, end_node = trace(start_node, fun, x) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/tracer.py", line 10, in trace end_box = fun(start_box) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/wrap_util.py", line 15, in unary_f return fun(*subargs, *kwargs) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 263, in fun ret = self.kl_div(x) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 175, in kl_div log_lik = self.log_lik(x) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 100, in log_lik ret = self._log_lik(x, vector=vector) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 122, in _log_lik demo = self._get_multipop_moran(x) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 134, in _get_multipop_moran demo = self.demo_func(x) File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 566, in _demo_fun raise ValueError(new_msg) from e ValueError: Exception encountered at parameter values {'tmig_north_south': 'Autograd ArrayBox with value 41.71175152831711', 'tmig_south_north': 'Autograd ArrayBox with value 22113.23235651487', 'tdiv': 'Autograd ArrayBox with value 22113.23235651487', 'n_north': 'Autograd ArrayBox with value 4037.628393936731', 'n_south': 'Autograd ArrayBox with value 14983.057878609443', 'N_pop1_pop2': 'Autograd ArrayBox with value 2079.969299777001', 'mfrac_north_south': 'Autograd ArrayBox with value 0.9948611983585026', 'mfrac_south_north': 'Autograd ArrayBox with value 0.5805977886307966'} (internal scaling x = {'tmig_north_south': 'Autograd ArrayBox with value 41.71175152831711', 'tmig_south_north': 'Autograd ArrayBox with value 22113.23235651487', 'tdiv': 'Autograd ArrayBox with value 1e-12', 'n_north': 'Autograd ArrayBox with value 8.303412767381138', 'n_south': 'Autograd ArrayBox with value 9.614675366987658', 'N_pop1_pop2': 'Autograd ArrayBox with value 7.640108412863842', 'mfrac_north_south': 'Autograd ArrayBox with value 5.265783319683204', 'mfrac_south_north': 'Autograd ArrayBox with value 0.3252278516667415'})

jackkamm commented 1 year ago

Thanks for reporting.

I would guess the error has something to do with this:

'tdiv': 'Autograd ArrayBox with value 1e-12'

In the internal scaling the tdiv parameter was running away to 0 (so tdiv happens immediately after max(t_mig_south_north, t_mig_north_south)).

Still, it shouldn't crash like this.

I will look into this more. In the meantime, if you could share your data (e.g. by putting on Google Drive and emailing me a link), that would be helpful.

Also, how reproducible was this error? Does it happen on most/all of the bootstrap resamples, or just very rarely?

jackkamm commented 1 year ago

I think the following workaround should work for this case:

import autograd.numpy as anp
...
migration_model.add_time_param("tdiv")
...
migration_model.move_lineages(
    pop2, pop1,
    t=lambda params: .001 + params.tdiv + anp.max([
        params.tmig_north_south,
        params.tmig_south_north]),
    N="N_pop1_pop2"
)

In words, rather than putting a lower constraint on tdiv, instead make the join time an explicit function of max(tmig)+tdiv+epsilon, to ensure that the join time is not the same as the pulse time due to floating precision issues.

See also: https://momi2.readthedocs.io/en/latest/tutorial.html#Using-functions-of-model-parameters-as-demographic-parameters

pengyan19 commented 1 year ago

Hello, I would like to provide the sfs file, this the code, and you can run this code, such as "python code.py South North". By the way, what the meaning of the parameter,such as "options = {"maxiter":5000, "ftol":0.84e-8})" ?

------------------ 原始邮件 ------------------ 发件人: "popgenmethods/momi2" @.>; 发送时间: 2022年10月24日(星期一) 凌晨1:38 @.>; @.**@.>; 主题: Re: [popgenmethods/momi2] momi.events.DemographyError (Issue #59)

I think the following workaround should work for this case: import autograd.numpy as anp ... migration_model.add_time_param("tdiv") ... migration_model.move_lineages( pop2, pop1, t=lambda params: .001 + params.tdiv + anp.max([ params.tmig_north_south, params.tmig_south_north]), N="N_pop1_pop2" )
In words, rather than putting a lower constraint on tdiv, instead making the join time an explicit function of max(tmig)+tdiv+epsilon, to ensure that the join time is not the same as the pulse time due to floating precision issues.

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: @.***>

jackkamm commented 1 year ago

I'm pretty sure of what the problem is now so no need to provide the SFS.

For short term, I recommend to use the workaround with lambda params I described above, and force tdiv to be strictly greater than the migration times plus some epsilon. Using lambda params is explained more in the docs here: https://momi2.readthedocs.io/en/latest/tutorial.html#Using-functions-of-model-parameters-as-demographic-parameters

Actually, I've had to use that workaround for some of my own work before when I ran into similar problems.

Give it a go and if it's still not clear, feel free to follow up with more questions here.

Fine to leave this issue open still, as I would like to eventually fix this so the workaround is unnecessary.

what the meaning of the parameter,such as "options = {"maxiter":5000, "ftol":0.84e-8})"

Those parameters are passed onto scipy.optimize.minimize, and control the convergence behavior. References here:

pengyan19 commented 1 year ago

Thank you very much. when I change the code as you description, I found that the running of the speed was very slowly. So, Shoud I change those code to run faster?

------------------ 原始邮件 ------------------ 发件人: "popgenmethods/momi2" @.>; 发送时间: 2022年10月24日(星期一) 中午11:46 @.>; @.**@.>; 主题: Re: [popgenmethods/momi2] momi.events.DemographyError (Issue #59)

I'm pretty sure of what the problem is now so no need to provide the SFS.

For short term, I recommend to use the workaround with lambda params I described above, and force tdiv to be strictly greater than the migration times plus some epsilon. Using lambda params is explained more in the docs here: https://momi2.readthedocs.io/en/latest/tutorial.html#Using-functions-of-model-parameters-as-demographic-parameters

Actually, I've had to use that workaround for some of my own work before when I ran into similar problems.

Give it a go and if it's still not clear, feel free to follow up with more questions here.

Fine to leave this issue open still, as I would like to eventually fix this so the workaround is unnecessary.

what the meaning of the parameter,such as "options = {"maxiter":5000, "ftol":0.84e-8})"

Those parameters are passed onto scipy.optimize.minimize, and control the convergence behavior. References here:

https://momi2.readthedocs.io/en/latest/api.html#momi.DemographicModel.optimize

https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.minimize-tnc.html#optimize-minimize-tnc

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: @.***>

jackkamm commented 1 year ago

Try this:

migration_model.move_lineages(
    pop2, pop1,
    t=lambda params: .001 + params.tdiv + anp.max(anp.array([
        params.tmig_north_south,
        params.tmig_south_north])),
    N="N_pop1_pop2"
)

The difference from the code above is that I wrap the list of tmig params in an anp.array before calling anp.max on it.

That's the only thing I can think of that might make it slower. Otherwise momi is essentially doing something similar internally so speed should be similar.

If it's still slow, should look into whether the likelihood evaluation is actually slower, or whether the optimizer is just taking more steps, and how reproducible it is.

pengyan19 commented 1 year ago

Hi,  when i change the code, i also meet a new question as follows:

/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/tracer.py:48: RuntimeWarning: overflow encountered in exp   return f_raw(*args, *kwargs) /public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/numpy/numpy_vjps.py:75: RuntimeWarning: invalid value encountered in double_scalars   defvjp(anp.exp,    lambda ans, x : lambda g: ans g) Traceback (most recent call last):   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 557, in _demo_fun     return self._get_demo(None)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 540, in _get_demo     G = _build_demo_graph(events, sampled_n_dict, params_dict, default_N=1.0)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 21, in _build_demo_graph     e.add_to_graph(_G, sample_sizes, params_dict)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 224, in add_to_graph     _set_sizes(G.nodes[k], t)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 411, in _set_sizes     raise DemographyError("Negative time or population size in {sizes}".format( momi.events.DemographyError: Negative time or population size in [{'t': '0.0', 'N': '1.0', 'growth_rate': 'None', 'tau': '0.0', 'N_top': '1.0'}, {'t': '0.0', 'N': 'Autograd ArrayBox with value nan', 'growth_rate': 'None', 'tau': '0.0', 'N_top': 'Autograd ArrayBox with value nan'}, {'t': '0.0', 'growth_rate': '0.0'}, {'t': 'Autograd ArrayBox with value nan'}]

The above exception was the direct cause of the following exception:

Traceback (most recent call last):   File "/public/home/pengyan/bioproject/01.ACB/24.momi2/06.momi-North-South/momi2_mig.py", line 50, in <module>     migration_model_copy.optimize()   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 910, in optimize     res = self._get_surface().find_mle(   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 267, in find_mle     return _find_minimum(fun, x0, scipy.optimize.minimize,   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/optimizers.py", line 88, in _find_minimum     ret = _find_minimum_helper(   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/optimizers.py", line 103, in _find_minimum_helper     return optimizer(f, start_params, opt_kwargs)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_minimize.py", line 684, in minimize     res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_tnc.py", line 416, in _minimize_tnc     rc, nf, nit, x, funv, jacv = moduleTNC.tnc_minimize(   File "_moduleTNC.pyx", line 201, in _moduleTNC.tnc_minimize   File "_moduleTNC.pyx", line 73, in _moduleTNC.function   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 285, in fun_and_grad     self._update_fun()   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 251, in _update_fun     self._update_fun_impl()   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 155, in update_fun     self.f = fun_wrapped(self.x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 137, in fun_wrapped     fx = fun(np.copy(x), args)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 76, in call     self._compute_if_needed(x, args)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 70, in _compute_if_needed     fg = self.fun(x, args)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/wrap_util.py", line 20, in nary_f     return unary_operator(unary_f, x, nary_op_args, nary_op_kwargs)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/differential_operators.py", line 135, in value_and_grad     vjp, ans = _make_vjp(fun, x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/core.py", line 10, in make_vjp     end_value, end_node =  trace(start_node, fun, x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/tracer.py", line 10, in trace     end_box = fun(start_box)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/wrap_util.py", line 15, in unary_f     return fun(*subargs, *kwargs)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 263, in fun     ret = self.kl_div(x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 175, in kl_div     log_lik = self.log_lik(x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 100, in log_lik     ret = self._log_lik(x, vector=vector)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 122, in _log_lik     demo = self._get_multipop_moran(x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 134, in _get_multipop_moran     demo = self.demo_func(x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 566, in _demo_fun     raise ValueError(new_msg) from e ValueError: Exception encountered at parameter values {'tmig_north_south': 'Autograd ArrayBox with value nan', 'tmig_south_north': 'Autograd ArrayBox with value 1e-12', 'tdiv': 'Autograd ArrayBox with value nan', 'n_north': 'Autograd ArrayBox with value nan', 'n_south': 'Autograd ArrayBox with value nan', 'N_pop1_pop2': 'Autograd ArrayBox with value nan', 'mfrac_north_south': 'Autograd ArrayBox with value nan', 'mfrac_south_north': 'Autograd ArrayBox with value nan'} (internal scaling x = {'tmig_north_south': 'Autograd ArrayBox with value nan', 'tmig_south_north': 'Autograd ArrayBox with value 1e-12', 'tdiv': 'Autograd ArrayBox with value nan', 'n_north': 'Autograd ArrayBox with value nan', 'n_south': 'Autograd ArrayBox with value nan', 'N_pop1_pop2': 'Autograd ArrayBox with value nan', 'mfrac_north_south': 'Autograd ArrayBox with value nan', 'mfrac_south_north': 'Autograd ArrayBox with value nan'})

------------------ 原始邮件 ------------------ 发件人: "popgenmethods/momi2" @.>; 发送时间: 2022年10月24日(星期一) 晚上9:02 @.>; @.**@.>; 主题: Re: [popgenmethods/momi2] momi.events.DemographyError (Issue #59)

Try this: migration_model.move_lineages( pop2, pop1, t=lambda params: .001 + params.tdiv + anp.max(anp.array([ params.tmig_north_south, params.tmig_south_north])), N="N_pop1_pop2" )
The difference from the code above is that I wrap the list of tmig params in an anp.array before calling anp.max on it.

That's the only thing I can think of that might make it slower. Otherwise momi is essentially doing something similar internally so speed should be similar.

If it's still slow, should look into whether the likelihood evaluation is actually slower, or whether the optimizer is just taking more steps, and how reproducible it is.

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: @.***>

pengyan19 commented 1 year ago

Hi,   Do you help deal with it?

------------------ 原始邮件 ------------------ 发件人: "空谷幽兰" @.>; 发送时间: 2022年10月25日(星期二) 上午8:44 @.>;

主题: 回复: [popgenmethods/momi2] momi.events.DemographyError (Issue #59)

Hi, when i change the code, i also meet a new question as follows:

/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/tracer.py:48: RuntimeWarning: overflow encountered in exp   return f_raw(*args, *kwargs) /public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/numpy/numpy_vjps.py:75: RuntimeWarning: invalid value encountered in double_scalars   defvjp(anp.exp,    lambda ans, x : lambda g: ans g) Traceback (most recent call last):   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 557, in _demo_fun     return self._get_demo(None)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 540, in _get_demo     G = _build_demo_graph(events, sampled_n_dict, params_dict, default_N=1.0)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 21, in _build_demo_graph     e.add_to_graph(_G, sample_sizes, params_dict)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 224, in add_to_graph     _set_sizes(G.nodes[k], t)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/events.py", line 411, in _set_sizes     raise DemographyError("Negative time or population size in {sizes}".format( momi.events.DemographyError: Negative time or population size in [{'t': '0.0', 'N': '1.0', 'growth_rate': 'None', 'tau': '0.0', 'N_top': '1.0'}, {'t': '0.0', 'N': 'Autograd ArrayBox with value nan', 'growth_rate': 'None', 'tau': '0.0', 'N_top': 'Autograd ArrayBox with value nan'}, {'t': '0.0', 'growth_rate': '0.0'}, {'t': 'Autograd ArrayBox with value nan'}]

The above exception was the direct cause of the following exception:

Traceback (most recent call last):   File "/public/home/pengyan/bioproject/01.ACB/24.momi2/06.momi-North-South/momi2_mig.py", line 50, in <module>     migration_model_copy.optimize()   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 910, in optimize     res = self._get_surface().find_mle(   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 267, in find_mle     return _find_minimum(fun, x0, scipy.optimize.minimize,   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/optimizers.py", line 88, in _find_minimum     ret = _find_minimum_helper(   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/optimizers.py", line 103, in _find_minimum_helper     return optimizer(f, start_params, opt_kwargs)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_minimize.py", line 684, in minimize     res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_tnc.py", line 416, in _minimize_tnc     rc, nf, nit, x, funv, jacv = moduleTNC.tnc_minimize(   File "_moduleTNC.pyx", line 201, in _moduleTNC.tnc_minimize   File "_moduleTNC.pyx", line 73, in _moduleTNC.function   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 285, in fun_and_grad     self._update_fun()   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 251, in _update_fun     self._update_fun_impl()   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 155, in update_fun     self.f = fun_wrapped(self.x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 137, in fun_wrapped     fx = fun(np.copy(x), args)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 76, in call     self._compute_if_needed(x, args)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 70, in _compute_if_needed     fg = self.fun(x, args)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/wrap_util.py", line 20, in nary_f     return unary_operator(unary_f, x, nary_op_args, nary_op_kwargs)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/differential_operators.py", line 135, in value_and_grad     vjp, ans = _make_vjp(fun, x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/core.py", line 10, in make_vjp     end_value, end_node =  trace(start_node, fun, x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/tracer.py", line 10, in trace     end_box = fun(start_box)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/wrap_util.py", line 15, in unary_f     return fun(*subargs, *kwargs)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 263, in fun     ret = self.kl_div(x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 175, in kl_div     log_lik = self.log_lik(x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 100, in log_lik     ret = self._log_lik(x, vector=vector)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 122, in _log_lik     demo = self._get_multipop_moran(x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/likelihood.py", line 134, in _get_multipop_moran     demo = self.demo_func(x)   File "/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py", line 566, in _demo_fun     raise ValueError(new_msg) from e ValueError: Exception encountered at parameter values {'tmig_north_south': 'Autograd ArrayBox with value nan', 'tmig_south_north': 'Autograd ArrayBox with value 1e-12', 'tdiv': 'Autograd ArrayBox with value nan', 'n_north': 'Autograd ArrayBox with value nan', 'n_south': 'Autograd ArrayBox with value nan', 'N_pop1_pop2': 'Autograd ArrayBox with value nan', 'mfrac_north_south': 'Autograd ArrayBox with value nan', 'mfrac_south_north': 'Autograd ArrayBox with value nan'} (internal scaling x = {'tmig_north_south': 'Autograd ArrayBox with value nan', 'tmig_south_north': 'Autograd ArrayBox with value 1e-12', 'tdiv': 'Autograd ArrayBox with value nan', 'n_north': 'Autograd ArrayBox with value nan', 'n_south': 'Autograd ArrayBox with value nan', 'N_pop1_pop2': 'Autograd ArrayBox with value nan', 'mfrac_north_south': 'Autograd ArrayBox with value nan', 'mfrac_south_north': 'Autograd ArrayBox with value nan'})

------------------ 原始邮件 ------------------ 发件人: "popgenmethods/momi2" @.>; 发送时间: 2022年10月24日(星期一) 晚上9:02 @.>; @.**@.>; 主题: Re: [popgenmethods/momi2] momi.events.DemographyError (Issue #59)

Try this: migration_model.move_lineages( pop2, pop1, t=lambda params: .001 + params.tdiv + anp.max(anp.array([ params.tmig_north_south, params.tmig_south_north])), N="N_pop1_pop2" )
The difference from the code above is that I wrap the list of tmig params in an anp.array before calling anp.max on it.

That's the only thing I can think of that might make it slower. Otherwise momi is essentially doing something similar internally so speed should be similar.

If it's still slow, should look into whether the likelihood evaluation is actually slower, or whether the optimizer is just taking more steps, and how reproducible it is.

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: @.***>

jackkamm commented 1 year ago

Hi,

It is hard to tell from the current information, but it seems like your model runs into numerical problems due to runaway behavior during optimization. One way to check this is to increase the logging level from INFO to DEBUG, then it will output the full optimization path.

Possibly, your model has too many parameters to reliably fit with your data, causing the runaway behavior. In that case, the best solution would be to simplify the model so that the likelihood surface is better behaved.

Alternatively, if the model usually converges but this error only occurs in a small number of bootstraps, then a practical workaround would be to wrap migration_model_copy.optimize() in a try-except clause, and then drop the failed samples, or try rerunning from different initial parameters.

Also regarding the original "Invalid move_lineages" error, I think there is a better solution than the lambda function I previously suggested. In particular, I think it would help to reverse the order of move_lineage() calls, in particular you should call them in this order:

migration_model.move_lineages(pop2, pop1, t="tmig_south_north", p="mfrac_south_north")
migration_model.move_lineages(pop1, pop2, t="tmig_north_south", p="mfrac_north_south")
migration_model.move_lineages(pop2, pop1, t="tdiv",N="N_pop1_pop2")

If the events are added in this order (backwards in time), then using the lambda function shouldn't be necessary anymore. The reason is that events are sorted by time before being evaluated, but when two events occur at the same time, then they are evaluated in the order they were added, which caused the original "Invalid move_lineages" event error when tdiv==tmig_south_north (because momi2 tried to evaluate the migration pulse after it already merged the 2 populations).

pengyan19 commented 1 year ago

Thanks, i meet some warning as follows, can it affect the results?

/public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/tracer.py:48: RuntimeWarning: overflow encountered in exp   return f_raw(*args, *kwargs) /public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/numpy/numpy_vjps.py:75: RuntimeWarning: invalid value encountered in double_scalars   defvjp(anp.exp,    lambda ans, x : lambda g: ans g) /public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py:312: RuntimeWarning: divide by zero encountered in divide   return np.log(x/np.array(1-x)) /public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_plotter.py:468: RuntimeWarning: invalid value encountered in double_scalars   -self.curr_g (nxt_t - self.curr_t)) /vol3/agis/xiaoyutao_group/pengyan/bioproject2/01.Chil/16.momi2/momi2.py:229: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.   pdf.savefig(fig.draw()) /public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/autograd/numpy/numpy_vjps.py:53: RuntimeWarning: overflow encountered in double_scalars   lambda ans, x, y : unbroadcast_f(y, lambda g: - g x / y*2)) /public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_plotter.py:468: RuntimeWarning: invalid value encountered in double_scalars   -self.curr_g (nxt_t - self.curr_t)) /public/home/pengyan/anaconda3/envs/momi2/lib/python3.10/site-packages/momi/demo_model.py:312: RuntimeWarning: divide by zero encountered in divide   return np.log(x/np.array(1-x)) /vol3/agis/xiaoyutao_group/pengyan/bioproject2/01.Chil/16.momi2/momi2.py:229: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.   pdf.savefig(fig.draw())

------------------ 原始邮件 ------------------ 发件人: "popgenmethods/momi2" @.>; 发送时间: 2022年11月7日(星期一) 上午9:13 @.>; @.**@.>; 主题: Re: [popgenmethods/momi2] momi.events.DemographyError (Issue #59)

Hi,

It is hard to tell from the current information, but it seems like your model runs into numerical problems due to runaway behavior during optimization. One way to check this is to increase the logging level from INFO to DEBUG, then it will output the full optimization path.

Possibly, your model has too many parameters to reliably fit with your data, causing the runaway behavior. In that case, the best solution would be to simplify the model so that the likelihood surface is better behaved.

If this error only occurs in a small number of bootstraps, then a practical workaround would be to wrap migration_model_copy.optimize() in a try-except, and then drop the failed samples, or try rerunning from different initial parameters.

Also regarding the original Invalid move_lineages error, I think there is a better solution than the lambda function I previously suggested. In particular, I think it would help to reverse the order of move_lineage() calls, in particular you should call them in this order: migration_model.move_lineages(pop2, pop1, t="tmig_south_north", p="mfrac_south_north") migration_model.move_lineages(pop1, pop2, t="tmig_north_south", p="mfrac_north_south") migration_model.move_lineages(pop2, pop1, t="tdiv",N="N_pop1_pop2")
If the events are added in this order (backwards in time), then using the lambda function shouldn't be necessary anymore. The reason is that events are sorted by time before being evaluated, but when two events occur at the same time, then they are evaluated in the order they were added, causing the original Invalid move_lineages event error when tdiv==tmig_north_south (because momi2 tries to evaluate the migration pulse after it already merged the 2 populations).

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: @.***>