RuntimeError: probability vector not in [0, 1] #81

ZheCai · 2018-04-14T14:58:28Z

Dear Terhorst,

While using SMC++ v1.13.0, I encountered the runtime errors:

1144 smcpp.analysis.base INFO theta: 0.000100
1145 smcpp.analysis.base INFO rho: 0.000100
1146 smcpp.data_filter INFO Loading data...
2194 smcpp.data_filter INFO 0.38 Gb of data
2203 smcpp.data_filter INFO 1 population
16977 smcpp.optimize.plugins.progress_printer INFO Starting EM algorithm...
16977 smcpp.optimize.plugins.progress_printer INFO EM iteration 1 of 20...
16977 smcpp.analysis.base INFO Running E-step
18590 smcpp.analysis.base INFO E-step completed
18605 smcpp.optimize.plugins.loglikelihood_monitor INFO Loglik: -678978.832221
21244 smcpp.optimize.plugins.parameter_optimizer INFO Updating rho, bounds (0.000001, 0.010000)
21784 smcpp.optimize.plugins.parameter_optimizer INFO New rho: 0.00208413
Traceback (most recent call last):
File "/build/Cellar/python/3.6.5/bin/smc++", line 11, in
load_entry_point('smcpp==1.13.0', 'console_scripts', 'smc++')()
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/frontend/console.py", line 26, in main
cmds[args.command].main(args)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/commands/estimate.py", line 73, in main
analysis.run()
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/analysis/base.py", line 109, in run
self._optimizer.run(niter or self._niter)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/optimize/optimizers.py", line 157, in run
res = self._minimize(x0, coords)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/optimize/optimizers.py", line 113, in _minimize
method=alg)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/scipy/optimize/_minimize.py", line 477, in minimize
return _minimize_powell(fun, x0, args, callback, **options)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 2594, in _minimize_powell
fx2 = squeeze(func(x2))
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 292, in function_wrapper
return function(*(wrapper_args + args))
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 63, in call
fg = self.fun(x, *args)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/optimize/optimizers.py", line 66, in _f
q = analysis.Q()
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/analysis/base.py", line 113, in Q
qq = [self._ims[pop].Q(separate=True) for pop in self._ims]
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/analysis/base.py", line 113, in
qq = [self._ims[pop].Q(separate=True) for pop in self._ims]
File "smcpp/_smcpp.pyx", line 296, in smcpp._smcpp._PyInferenceManager.Q
File "smcpp/_smcpp.pyx", line 291, in smcpp._smcpp._PyInferenceManager.Q
RuntimeError: probability vector not in [0, 1]

My command line was:
smc++ estimate --cores 8 --timepoints 10,300000 -o aro/ 6.5e-9 ../smcpp_02/input/aro.chr*.smc.gz

The debug file:
.debug.txt

My input:
aro.zip

Is there any problem with my command line? or my input?

Thank you!

terhorst · 2018-04-14T15:35:38Z

Try running it with the default time points instead of overriding them.

ZheCai · 2018-04-15T04:29:34Z

The new error ("RuntimeError: erroneous average coalescence time
") was reported when running with default time points.

2300 smcpp.analysis.base INFO theta: 0.000100
2301 smcpp.analysis.base INFO rho: 0.000100
2304 smcpp.data_filter INFO Loading data...
3504 smcpp.data_filter INFO 0.38 Gb of data
3515 smcpp.data_filter INFO 1 population
21559 smcpp.analysis.analysis INFO calculated t1: 3415.028805 gens
22914 smcpp.optimize.plugins.progress_printer INFO Starting EM algorithm...
22915 smcpp.optimize.plugins.progress_printer INFO EM iteration 1 of 20...
22915 smcpp.analysis.base INFO Running E-step
24325 smcpp.analysis.base INFO E-step completed
24331 smcpp.optimize.plugins.loglikelihood_monitor INFO Loglik: -667287.620700
27891 smcpp.optimize.plugins.parameter_optimizer INFO Updating rho, bounds (0.000001, 0.010000)
28498 smcpp.optimize.plugins.parameter_optimizer INFO New rho: 0.00139285
451579 smcpp.optimize.plugins.progress_printer INFO Current model:

-0.52 2.13 30.27 1.99 3.01 1.49 1.64 2.64
0.60 0.49 0.32 0.19 0.13 0.12 0.17 0.49 3.28 59.80 2611.97 195216.86 17738214.66 1392122649.28 67042834348.69 1407559519781.69 9152778775215.18 13120851855391.83 3676328104266.64 272201166353.91 7553596112.07 111421800.65 1239092.69 14733.95 265.69 10.31 1.15 0.35 0.24 0.31 0.64 1.72 5.06 13.55 27.73 40.57 44.60 38.93 28.49 18.47 11.21 6.73 4.22 2.88 2.14 1.72 1.48 1.36 1.33 1.37 1.47 1.63 1.88 2.21 2.64 3.19 3.87 4.68 5.60 6.61 7.71 8.87 10.09 11.33 12.59 13.83 15.04 16.20 17.28 18.26 19.14 19.89 20.52 21.01 21.38 21.61 21.72 21.71 21.59 21.38 21.09 20.73 20.32 19.85 19.36 18.85 18.33 17.81 17.30 16.80 16.33 15.89 15.48 15.12 14.79 14.52 14.30 14.14 14.04 14.00
452379 smcpp.optimize.plugins.ascii_plotter INFO Plot of current model:
N_e
1e+18 ++--++-+++++-----+--+--+-++-++++----+---+-+-+-+++++-----+--+-+-+++
+ A* + +Pop. 1 ****** +
1e+16 ++ * * ++
+ * * +
1e+14 ++ * * ++
+ * * +
1e+12 ++ * * ++
1e+10 ++ * * ++
+ * * +
1e+08 ++ * * ++
+ * * +
1e+06 ++ * * ++
+ A* A A***** *************A
10000 + * * ** AA ++
+ + ** + + +
100 ++--++-+++++-----+--+--+-++-++++----+---+-+-+-+++++-----+--+-+-+++
1000 10000 100000
Generations

452380 smcpp.optimize.plugins.progress_printer INFO EM iteration 2 of 20...
452384 smcpp.analysis.base INFO Running E-step
454172 smcpp.analysis.base INFO E-step completed
454177 smcpp.optimize.plugins.loglikelihood_monitor INFO New loglik: -642761.904746 (old: -667287.620700 [3.675434%])
460712 smcpp.optimize.plugins.parameter_optimizer INFO Updating rho, bounds (0.000001, 0.010000)
461361 smcpp.optimize.plugins.parameter_optimizer INFO New rho: 0.0026181
Traceback (most recent call last):
File "/build/Cellar/python/3.6.5/bin/smc++", line 11, in
load_entry_point('smcpp==1.13.0', 'console_scripts', 'smc++')()
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/frontend/console.py", line 26, in main
cmds[args.command].main(args)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/commands/estimate.py", line 73, in main
analysis.run()
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/analysis/base.py", line 109, in run
self._optimizer.run(niter or self._niter)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/optimize/optimizers.py", line 157, in run
res = self._minimize(x0, coords)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/optimize/optimizers.py", line 113, in _minimize
method=alg)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/scipy/optimize/_minimize.py", line 477, in minimize
return _minimize_powell(fun, x0, args, callback, **options)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 2594, in _minimize_powell
fx2 = squeeze(func(x2))
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 292, in function_wrapper
return function(*(wrapper_args + args))
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/scipy/optimize/optimize.py", line 63, in call
fg = self.fun(x, *args)
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/optimize/optimizers.py", line 66, in _f
q = analysis.Q()
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/analysis/base.py", line 113, in Q
qq = [self._ims[pop].Q(separate=True) for pop in self._ims]
File "/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/analysis/base.py", line 113, in
qq = [self._ims[pop].Q(separate=True) for pop in self._ims]
File "smcpp/_smcpp.pyx", line 296, in smcpp._smcpp._PyInferenceManager.Q
File "smcpp/_smcpp.pyx", line 291, in smcpp._smcpp._PyInferenceManager.Q
RuntimeError: erroneous average coalescence time

Debug file:
.debug.txt

ZheCai · 2018-04-15T04:30:01Z

@terhorst

terhorst · 2018-04-16T21:02:01Z

The basic issue is that the model appears to diverge during fitting. I would try first estimating a simple piecewise model:

smc++ estimate 6.5e-9 aro/*.smc.gz

Then feed the results of that into a second estimation command where you use splines:

smc++ estimate 6.5e-9 aro/*.smc.gz --spline cubic --initial-model model.final.json

ZheCai · 2018-04-17T07:47:15Z

I will try it. Thank you for your reply @terhorst !

ZheCai · 2018-04-17T09:39:47Z

Hi, @terhorst
Following your suggestion, the estimating run seems normal for my data of aro population. But in another population, the Stacktrace infomation was reported:

stack trace:
/lib64/libc.so.6() [0x39b1632920]
/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/_smcpp.cpython-36m-x86_64-linux-gnu.so : Eigen::Matrix<double, -1, 1, 0, -1, 1>::Matrix<Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double, double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const> const, Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const> const> >(Eigen::CwiseBinaryOp<Eigen::internal::scalar_sum_op<double, double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const> const, Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<double, double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> const> const> const&)+0x8e
/build/Cellar/python/3.6.5/lib/python3.6/site-packages/smcpp/_smcpp.cpython-36m-x86_64-linux-gnu.so : ()+0x1337c7
/build/Cellar/gcc/4.8.4/lib64/libgomp.so.1 : ()+0x830a
/lib64/libpthread.so.0() [0x39b1a07851]
/lib64/libc.so.6 : clone()+0x6d

Here is the debug file:
.debug.txt

the input file:
jap.zip

ZheCai · 2018-04-17T13:06:04Z

I run SMC++ with the default time points and piecewise mode, but the error "probability vector not in [0, 1]" was still reported in some samples.

the command:
smc++ estimate --cores 8 -o ind3/ 6.5e-9 ../smcpp_02/input/ind3.chr*.smc.gz

terhorst · 2018-04-19T20:48:19Z

It's again running off to some huge or small value of Ne. You can try decreasing the -rp option. Alternatively in v1.13.1 I added options to bound the effective population size.

ZheCai · 2018-04-20T03:50:07Z

Thanks a lot for your kind help. I am running the model with v1.13.1. The running is ok so far, and I will post the results here when finished.

terhorst closed this as completed in 1858d07 Apr 18, 2018

terhorst added a commit that referenced this issue Aug 21, 2018

[estimate] add bounds. Hopefully fixes #81

cfde840

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

RuntimeError: probability vector not in [0, 1] #81

RuntimeError: probability vector not in [0, 1] #81

ZheCai commented Apr 14, 2018

terhorst commented Apr 14, 2018

ZheCai commented Apr 15, 2018

ZheCai commented Apr 15, 2018

terhorst commented Apr 16, 2018

ZheCai commented Apr 17, 2018

ZheCai commented Apr 17, 2018

ZheCai commented Apr 17, 2018

terhorst commented Apr 19, 2018

ZheCai commented Apr 20, 2018

RuntimeError: probability vector not in [0, 1] #81

RuntimeError: probability vector not in [0, 1] #81

Comments

ZheCai commented Apr 14, 2018

terhorst commented Apr 14, 2018

ZheCai commented Apr 15, 2018

ZheCai commented Apr 15, 2018

terhorst commented Apr 16, 2018

ZheCai commented Apr 17, 2018

ZheCai commented Apr 17, 2018

ZheCai commented Apr 17, 2018

terhorst commented Apr 19, 2018

ZheCai commented Apr 20, 2018