Errors in using parallel processing for models with starting point bias

[zaizibai]

Dear pyddm developer,

Thanks for making this wonderful package. I am using the DDM/collapsing boundary model with starting point bias parameter. I noticed that for models with starting point parameter, the fitting speed is extremely slow which is quite different from using the other packages which utilize the Navarro likelihood. Therefore i was thinking maybe parallel computing might faster the fitting process. However, i confronted the following error which only happens in when the starting point parameter was added:

---

RemoteTraceback Traceback (most recent call last)
RemoteTraceback:
"""
Traceback (most recent call last):
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\multiprocess\pool.py", line 119, in worker
result = (True, func(*args, **kwds))
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\multiprocess\pool.py", line 44, in mapstar
return list(map(*args))
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\paranoid\decorators.py", line 114, in _decorated
return func(*args, **kwargs)
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\pyddm\model.py", line 550, in solve
return self.solve_numerical_implicit(conditions=conditions, return_evolution=return_evolution, force_python=force_python)
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\paranoid\decorators.py", line 114, in _decorated
return func(*args, **kwargs)
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\pyddm\model.py", line 913, in solve_numerical_implicit
return self.solve_numerical(method="implicit", conditions=conditions, **kwargs)
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\paranoid\decorators.py", line 114, in _decorated
return func(*args, **kwargs)
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\pyddm\model.py", line 727, in solve_numerical
return self.solve_numerical_c(conditions=conditions)
File "C:\Users\KLIPS\anaconda3\envs\pyddm\lib\site-packages\pyddm\model.py", line 676, in solve_numerical_c
res = csolve.implicit_time(drift, drifttype, noise, noisetype, bound, boundtype, self.get_dependence("IC").get_IC(self.x_domain(conditions=conditions), self.dx, conditions=conditions), self.T_dur, self.dt, self.dx, len(self.t_domain()))
File "", line 11, in get_IC
NameError: name 'np' is not defined
"""

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

NameError Traceback (most recent call last)
in
----> 1 fit_models = fit_adjust_model(sample=samp_sub_dat,model=model, lossfunction=LossRobustBIC)

~\anaconda3\envs\pyddm\lib\site-packages\pyddm\functions.py in fit_adjust_model(sample, model, fitparams, fitting_method, lossfunction, verify, method, verbose)
366 if "disp" not in fitparams.keys():
367 fitparams["disp"] = verbose
--> 368 x_fit = differential_evolution(_fit_model, constraints, **fitparams)
369 elif fitting_method == "hillclimb":
370 x_fit = evolution_strategy(_fit_model, x_0, **fitparams)

~\anaconda3\envs\pyddm\lib\site-packages\scipy\optimize_differentialevolution.py in differential_evolution(func, bounds, args, strategy, maxiter, popsize, tol, mutation, recombination, seed, callback, disp, polish, init, atol, updating, workers, constraints)
306 workers=workers,
307 constraints=constraints) as solver:
--> 308 ret = solver.solve()
309
310 return ret

~\anaconda3\envs\pyddm\lib\site-packages\scipy\optimize_differentialevolution.py in solve(self)
749 self.population_energies[self.feasible] = (
750 self._calculate_population_energies(
--> 751 self.population[self.feasible]))
752
753 self._promote_lowest_energy()

~\anaconda3\envs\pyddm\lib\site-packages\scipy\optimize_differentialevolution.py in _calculate_population_energies(self, population)
872 try:
873 calc_energies = list(self._mapwrapper(self.func,
--> 874 parameters_pop[0:nfevs]))
875 energies[0:nfevs] = np.squeeze(calc_energies)
876 except (TypeError, ValueError):

~\anaconda3\envs\pyddm\lib\site-packages\scipy\optimize_differentialevolution.py in call(self, x)
1259
1260 def call(self, x):
-> 1261 return self.f(x, *self.args)
1262
1263

~\anaconda3\envs\pyddm\lib\site-packages\pyddm\functions.py in _fit_model(xs)
348 x = p.minval
349 s(m, x)
--> 350 lossf = lf.loss(m)
351 if verbose:
352 _logger.info(repr(m) + " loss="+ str(lossf))

~\anaconda3\envs\pyddm\lib\site-packages\paranoid\decorators.py in _decorated(*args, **kwargs)
112 # Skip verification if paranoid is disabled.
113 if Settings.get("enabled", function=func) == False:
--> 114 return func(*args, **kwargs)
115 # We only run this function once for performance reasons, and
116 # then pass it as an argument to each check function.

~\anaconda3\envs\pyddm\lib\site-packages\pyddm\models\loss.py in loss(self, model)
213 @requires("model.dt == self.dt and model.T_dur == self.T_dur")
214 def loss(self, model):
--> 215 loglikelihood = -LossLikelihood.loss(self, model)
216 return np.log(self.samplesize)self.nparams - 2loglikelihood
217

~\anaconda3\envs\pyddm\lib\site-packages\paranoid\decorators.py in _decorated(*args, **kwargs)
112 # Skip verification if paranoid is disabled.
113 if Settings.get("enabled", function=func) == False:
--> 114 return func(*args, **kwargs)
115 # We only run this function once for performance reasons, and
116 # then pass it as an argument to each check function.

~\anaconda3\envs\pyddm\lib\site-packages\pyddm\models\loss.py in loss(self, model)
164 def loss(self, model):
165 assert model.dt == self.dt and model.T_dur == self.T_dur
--> 166 sols = self.cache_by_conditions(model)
167 loglikelihood = 0
168 for k in sols.keys():

~\anaconda3\envs\pyddm\lib\site-packages\pyddm\models\loss.py in cache_by_conditions(self, model)
91 """
92 from ..functions import solve_all_conditions
---> 93 return solve_all_conditions(model, sample=self.sample, method=self.method)
94 @paranoidclass
95 class LossSquaredError(LossFunction):

~\anaconda3\envs\pyddm\lib\site-packages\pyddm\functions.py in solve_all_conditions(model, sample, condition_combinations, method)
528 # since you can't broadcast commands to all processes
529 _parallel_pool.map(lambda x : paranoid_settings.set(enabled=False), [None]_parallel_pool.n_cpus2)
--> 530 sols = _parallel_pool.map(meth, conds, chunksize=1)
531 for c,s in zip(conds, sols):
532 cache[frozenset(c.items())] = s

~\anaconda3\envs\pyddm\lib\site-packages\multiprocess\pool.py in map(self, func, iterable, chunksize)
264 in a list that is returned.
265 '''
--> 266 return self._map_async(func, iterable, mapstar, chunksize).get()
267
268 def starmap(self, func, iterable, chunksize=None):

~\anaconda3\envs\pyddm\lib\site-packages\multiprocess\pool.py in get(self, timeout)
642 return self._value
643 else:
--> 644 raise self._value
645
646 def _set(self, i, obj):

NameError: name 'np' is not defined

[mwshinn]

Thanks for the bug report! Could you please post your full model? The error you're getting looks like a problem with the model rather than with the parallelisation.

Regarding the slowness with the starting points: the last released version of PyDDM only uses this (fast) analytical solution when you use ICPoint, not ICPointRelative or any other similar single point starting positions. This is fixed in the latest version which should be released in the next month or so. If you are eager, you can find it at: https://github.com/mwshinn/PyDDM/tree/dev

[mwshinn]

I should also mention: if you have parallelisation issues, one thing that usually fixes the problem is to put your custom model components in a separate python file and then import them.

[zaizibai]

Thanks for the fast respond. Here is my model:

``

  class DriftFun(Drift):
      name = "Drift depends linearly on valence, movement and their interaction"
      required_parameters = ["drift_beta0","drift_beta1","drift_beta2","drift_beta3"] # drift rate parameter
      required_conditions = ["valance","movement",] 
    
            def get_drift(self, conditions, **kwargs):
                return self.drift_beta0 + self.drift_beta1 * conditions['valance'] + self.drift_beta2 * conditions['movement'] + self.drift_beta3 * conditions['movement']*conditions['valance']

    class ICPointFun(InitialCondition):
        name = "approach bias"
        required_parameters = ["x0"]
        def get_IC(self, x, dx, conditions):
            x0 = self.x0/2 + .5 
            shift_i = int((len(x)-1)*x0)
            assert shift_i >= 0 and shift_i < len(x), "Invalid initial conditions"
            pdf = np.zeros(len(x))
            pdf[shift_i] = 1.
            return pdf
    model_1 = Model(drift=DriftFun(drift_beta0=Fittable(minval=-2, maxval=2),
                                      drift_beta1=Fittable(minval=-2, maxval=2),
                                      drift_beta2=Fittable(minval=-2, maxval=2),
                                      drift_beta3=Fittable(minval=-2, maxval=2)),
                  bound=BoundConstant(B=Fittable(minval=0.1, maxval=2.2)),
                  IC=ICPointFun(x0=Fittable(minval=-1, maxval=1)),
                  overlay=OverlayNonDecision(nondectime=Fittable(minval=0.1, maxval=3)),
                  dx=.01, dt=.01,T_dur=7.2)

``
If i delete the IC=ICPointFun(x0=Fittable(minval=-1, maxval=1)) or don't use the parallel fitting, it works fine.

[zaizibai]

BTW, thanks for mentioning the next version of pyddm will use the analytical solution if the model is standard DDM. But I feel like we should always use the same protocol for fitting and comparing different models. If models without analytical solution are included in the comparison metric, then models even with the analytical solution, like DDM should also be fitted with solving fokker-planck equation.

[mwshinn]

If you see my other comment, putting it in its own separate file and importing numpy in that file should do it. You could also just add "import numpy as np" to the first line of get_IC and that should do it too. Parallelisation is weird on some platforms but this should also be better when the new version of pyddm is released.

You can force analytical or numerical (FP) by using "solve_analytical()" or "solve_numerical()" instead of "solve()". Or in the fitting, you can pass the "method" argument of your choice to fit_adjust_model. Though the numerical one should converge to the analytical one anyway.