Question about： Encountered AssertionError When Adding a Migration Event

[LeiHuang66]

Hello jackkamm,

First of all, thank you for creating and maintaining this project! It has been incredibly helpful for my development work.
However, I encountered an AssertionError when trying to add a migration event, and I was hoping you could help me troubleshoot this issue:

Below is the code

#!/bin/python

-- coding: utf-8 --

import momi
import logging
logging.basicConfig(level=logging.INFO, filename="4mig1.log")

sfs = momi.Sfs.load("../momi_data/momi_data2")
model = momi.DemographicModel(
N_e=6e4,
gen_time=1,
muts_per_gen=8e-9
)
model.set_data(sfs, 784692364)

model.add_leaf("XZ", N="N_XZ")
model.add_leaf("SC", N="N_SC")
model.add_leaf("NMG", N="N_NMG")
model.add_leaf("XJ", N="N_XJ")
model.add_leaf("ELS3", N="N_ELS3")
model.add_leaf("ELS4", N="N_ELS4")
model.add_leaf("ELS2", N="N_ELS2")

model.move_lineages("NMG", "SC", t="t_split_SC_NMG", N="N_SCA")
model.move_lineages("SC", "XZ", t="t_split_XZ_SC", N="N_XZA")
model.move_lineages("ELS3", "XJ", t="t_split_XJ_ELS3", N="N_XJA")
model.move_lineages("ELS4", "ELS3", t="t_split_ELS3_ELS4", N="N_ELS3A")
model.move_lineages("ELS2", "ELS4", t="t_split_ELS4_ELS2", N="N_ELS4A")
model.move_lineages("XJ", "XZ", t="t_split_XZ_XJ", N="N_A")

model.move_lineages("XJ", "XZ", t="t_pulse_XZ_XJ", p="p_XZ_XJ")
model.move_lineages("ELS3", "XJ", t="t_pulse_XJ_ELS3", p="p_XJ_ELS3")
model.move_lineages("SC", "XZ", t="t_pulse_XZ_SC", p="p_XZ_SC")
model.move_lineages("NMG", "SC", t="t_pulse_SC_NMG", p="p_SC_NMG")
model.move_lineages("ELS2", "ELS4", t="t_pulse_ELS4_ELS2", p="p_ELS4_ELS2")

model.add_size_param("N_XZ", lower=8.0e4, upper=9.2e4)
model.add_size_param("N_SC", lower=7.5e3, upper=9.5e3)
model.add_size_param("N_NMG", lower=1.2e4, upper=1.6e4)
model.add_size_param("N_XJ", lower=4.2e4, upper=6.0e4)
model.add_size_param("N_ELS3", lower=8.5e3, upper=1.0e4)
model.add_size_param("N_ELS4", lower=1.7e4, upper=2.2e4)
model.add_size_param("N_ELS2", lower=8.5e3, upper=1.0e4)

model.add_time_param("t_split_XZ_XJ", lower=6.0e4, upper=7.7e4)
model.add_time_param("t_split_XJ_ELS3", lower=3.3e4, upper=4.5e4)
model.add_time_param("t_split_XZ_SC", lower=2.5e4, upper=3e4)
model.add_time_param("t_split_ELS3_ELS4", lower=1.7e4, upper=2.5e4)
model.add_time_param("t_split_SC_NMG", lower=4.5e3, upper=6.5e3)
model.add_time_param("t_split_ELS4_ELS2", lower=3.5e3, upper=4.8e3)

model.add_size_param("N_SCA", lower=6.5e5, upper=8.0e5)
model.add_size_param("N_XZA", lower=5.5e4, upper=7.5e4)
model.add_size_param("N_XJA", lower=1.4e6, upper=1.8e6)
model.add_size_param("N_ELS3A", lower=1.0e4, upper=1.3e4)
model.add_size_param("N_ELS4A", lower=1.7e4, upper=2.5e4)

model.add_size_param("N_A", lower=8.5e4, upper=1.0e5)

model.add_time_param("t_pulse_XZ_XJ", upper_constraints=["t_split_XZ_XJ"])
model.add_time_param("t_pulse_XJ_ELS3", upper_constraints=["t_split_XJ_ELS3"])
model.add_time_param("t_pulse_XZ_SC", upper_constraints=["t_split_XZ_SC"])
model.add_time_param("t_pulse_SC_NMG", upper_constraints=["t_split_SC_NMG"])
model.add_time_param("t_pulse_ELS4_ELS2", upper_constraints=["t_split_ELS4_ELS2"])

model.add_pulse_param("p_XZ_XJ", lower=0.05, upper=0.15)
model.add_pulse_param("p_XJ_ELS3", lower=0.05, upper=0.1)
model.add_pulse_param("p_XZ_SC", lower=0.05, upper=0.1)
model.add_pulse_param("p_SC_NMG", lower=0.05, upper=0.1)
model.add_pulse_param("p_ELS4_ELS2", lower=0.05, upper=0.1)

results = []
model_copy = model.copy()
model_copy.set_params(model.get_params(),randomize=True)

The code I provided above runs successfully without any issues. However, the error occurs after this step when I attempt to

results.append(model_copy.optimize(method="L-BFGS-B"),options={'maxiter': 20,"ftol":1e-7})

#The error message I received is as follows:

---

AssertionError Traceback (most recent call last)
/tmp/ipykernel_153943/2232076006.py in
----> 1 results.append(model_copy.optimize(method="L-BFGS-B"),options={'maxiter': 20,"ftol":1e-7})

~/anaconda3/envs/momi/lib/python3.7/site-packages/momi/demo_model.py in optimize(self, method, jac, hess, hessp, printfreq, **kwargs)
912 jac=jac, hess=hess, hessp=hessp,
913 bounds=bounds, callback=callback,
--> 914 **kwargs)
915
916 self._set_x(res.x)

~/anaconda3/envs/momi/lib/python3.7/site-packages/momi/likelihood.py in find_mle(self, x0, method, jac, hess, hessp, bounds, callback, **kwargs)
267 return _find_minimum(fun, x0, scipy.optimize.minimize,
268 bounds=bounds, callback=callback,
--> 269 opt_kwargs=opt_kwargs, gradmakers=gradmakers, replacefun=replacefun)
270
271

~/anaconda3/envs/momi/lib/python3.7/site-packages/momi/optimizers.py in _find_minimum(f, start_params, optimizer, bounds, callback, opt_kwargs, **kwargs)
87
88 ret = _find_minimum_helper(
---> 89 f, start_params, optimizer, opt_kwargs, **kwargs)
90 if fixed_params:
91 ret.x = get_x(ret.x)
...
---> 47 assert np.all(mins <= tol * maxes)
48
49 if axis is not None:

AssertionError:
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...

I verified that all dependencies and earlier steps were applied correctly.Could you help me understand what might be causing this issue? Please let me know if more details or code are needed. Thank you!

[jackkamm]

Sorry for the delay and thanks for reporting.

It appears that assertion fails in truncate0 in util.py. There's a couple places where we're dealing with a matrix that should be non-negative, but might have very small negative numbers due to numerical precision issues. This function ensures the negative numbers are very small and then sets them to 0.

The assertion is failing because the negative numbers are larger than allowed by its tolerance. It might help to print verbose output during the optimization to get an idea if there's some runaway behavior. I think calling

import logging
logging.basicConfig(level=logging.INFO)

when starting up the python session should do that.

Otherwise, invoking the interactive debugger to print variables at the time of the assertion failure may be informative. I think running a script with python -m pdb should do that. But this is rather advanced.

It's possible the tolerance is too strict and things would run fine without. But more likely the numerical precision issues would make the output unreliable. But you could try it out by disabling assertions to see what happens:
https://stackoverflow.com/questions/1273211/disable-assertions-in-python

(Note that using assertions in this way isn't best practice, and I would take a different approach if I were to implement this again)