Make C solver the default

README.md

@@ -53,6 +53,8 @@ If installing from source, [download the source code](https://github.com/mwshinn
 - Matplotlib
 - [Paranoid Scientist](<https://github.com/mwshinn/paranoidscientist>)
 - Pathos (optional, for multiprocessing support)
+- A C compiler (If you don't already have one, the easiest way to install one
+  may be by installing Cython.)
 
 
 ## Contact

ddm/csolve.c

@@ -25,7 +25,7 @@ void easy_dgtsv(int n, double *dl, double *d, double *du, double *b) {
 }
 
 double* _analytic_ddm_linbound(double a1, double b1, double a2, double b2, int nsteps, double tstep);
-int _implicit_time(int Tsteps, double *pdfcorr, double *pdferr, double* drift, double* noise, double *bound, double *ic, int Xsteps, double dt, double dx, uint drift_mode, uint noise_mode, uint bound_mode);
+int _implicit_time(int Tsteps, double *pdfcorr, double *pdferr, double *pdfcurr, double* drift, double* noise, double *bound, double *ic, int Xsteps, double dt, double dx, uint drift_mode, uint noise_mode, uint bound_mode);
 
 static PyObject* analytic_ddm_linbound(PyObject* self, PyObject* args) {
   double a1, b1, a2, b2, tstep;
@@ -47,19 +47,17 @@ static PyObject* implicit_time(PyObject* self, PyObject* args) {
   PyObject *__drift, *__noise, *__bound, *__ic;
   int nsteps, len_x0;
   int drifttype, noisetype, boundtype;
-  if (!PyArg_ParseTuple(args, "OiOiOiOddd", &__drift, &drifttype, &__noise, &noisetype, &__bound, &boundtype, &__ic, &T_dur, &dt, &dx))
+  if (!PyArg_ParseTuple(args, "OiOiOiOdddi", &__drift, &drifttype, &__noise, &noisetype, &__bound, &boundtype, &__ic, &T_dur, &dt, &dx, &nsteps))
     return NULL;
-  nsteps = (int)(T_dur/dt)+1;
-  len_x0 = (int)(1/dx*2)+1;
 
   _drift = (PyArrayObject*)PyArray_FROMANY(__drift, NPY_DOUBLE, 1, 1, NPY_ARRAY_C_CONTIGUOUS);
   _noise = (PyArrayObject*)PyArray_FROMANY(__noise, NPY_DOUBLE, 1, 1, NPY_ARRAY_C_CONTIGUOUS);
   _bound = (PyArrayObject*)PyArray_FROMANY(__bound, NPY_DOUBLE, 1, 1, NPY_ARRAY_C_CONTIGUOUS);
   _ic = (PyArrayObject*)PyArray_FROMANY(__ic, NPY_DOUBLE, 1, 1, NPY_ARRAY_C_CONTIGUOUS);
+
+  len_x0 = PyArray_SIZE(_ic);
   if (!_drift || !_noise || !_bound || !_ic)
     return NULL;
-  if (PyArray_SIZE(_ic) != len_x0)
-    return NULL;
   drift = (double*)PyArray_DATA(_drift);
   noise = (double*)PyArray_DATA(_noise);
   bound = (double*)PyArray_DATA(_bound);
@@ -68,13 +66,17 @@ static PyObject* implicit_time(PyObject* self, PyObject* args) {
 
   double *pdfcorr = (double*)malloc(nsteps*sizeof(double));
   double *pdferr = (double*)malloc(nsteps*sizeof(double));
-  _implicit_time(nsteps, pdfcorr, pdferr, drift, noise, bound, ic, len_x0, dt, dx, drifttype, noisetype, boundtype);
+  double *pdfcurr = (double*)malloc(len_x0*sizeof(double));
+  _implicit_time(nsteps, pdfcorr, pdferr, pdfcurr, drift, noise, bound, ic, len_x0, dt, dx, drifttype, noisetype, boundtype);
   npy_intp dims[1] = { nsteps };
+  npy_intp dimscurr[1] = { len_x0 };
   PyObject *corrarray = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, pdfcorr);
   PyObject *errarray = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, pdferr);
+  PyObject *currarray = PyArray_SimpleNewFromData(1, dimscurr, NPY_DOUBLE, pdfcurr);
   PyArray_UpdateFlags((PyArrayObject*)corrarray, NPY_ARRAY_OWNDATA);
   PyArray_UpdateFlags((PyArrayObject*)errarray, NPY_ARRAY_OWNDATA);
-  PyObject *ret = Py_BuildValue("(OO)", corrarray, errarray);
+  PyArray_UpdateFlags((PyArrayObject*)currarray, NPY_ARRAY_OWNDATA);
+  PyObject *ret = Py_BuildValue("(OOO)", corrarray, errarray, currarray);
   return ret;
 }
 
@@ -156,6 +158,8 @@ double* _analytic_ddm_linbound(double a1, double b1, double a2, double b2, int n
     //suminc[i] *= exp(-pow((a1+b1*i*tstep),2)/(i*tstep)/2)/sqrtpi/pow(i*tstep, 1.5);
     //itstep = i*tstep;
     suminc[i] *= oneoversqrtpi*exp(-.5*(a1+b1*i*tstep)*(a1+b1*i*tstep)/(i*tstep))*pow(i*tstep, -1.5);
+    if (suminc[i] < 0)
+      suminc[i] = 0;
   }
   suminc[0] = 0;
   free(tmp);
@@ -164,7 +168,7 @@ double* _analytic_ddm_linbound(double a1, double b1, double a2, double b2, int n
 
 
 
-int _implicit_time(int Tsteps, double *pdfcorr, double *pdferr, double* drift, double* noise, double *bound, double *ic, int Xsteps, double dt, double dx, uint drift_mode, uint noise_mode, uint bound_mode) {
+int _implicit_time(int Tsteps, double *pdfcorr, double *pdferr, double *pdfcurr, double* drift, double* noise, double *bound, double *ic, int Xsteps, double dt, double dx, uint drift_mode, uint noise_mode, uint bound_mode) {
   int dmultt=-1, dmultx=-1, nmultt=-1, nmultx=-1, bmultt=-1;
   int j;
   double dxinv = 1/dx;
@@ -177,7 +181,6 @@ int _implicit_time(int Tsteps, double *pdfcorr, double *pdferr, double* drift, d
   double *DU_copy = (double*)malloc((Xsteps-1)*sizeof(double));
   double *D_copy = (double*)malloc(Xsteps*sizeof(double));
   double *DL_copy = (double*)malloc((Xsteps-1)*sizeof(double));
-  double *pdfcurr = (double*)malloc(Xsteps*sizeof(double));
   double *pdfcurr_copy = (double*)malloc(Xsteps*sizeof(double));
   memset(pdfcorr, 0, Tsteps*sizeof(double));
   memset(pdferr, 0, Tsteps*sizeof(double));
@@ -311,7 +314,6 @@ int _implicit_time(int Tsteps, double *pdfcorr, double *pdferr, double* drift, d
     pdferr[i] *= dtinv;
   }
   // Free memory
-  free(pdfcurr);
   free(pdfcurr_copy);
   free(D);
   free(DU);

ddm/model.py

@@ -9,8 +9,6 @@
 import scipy.sparse.linalg
 from .tridiag import TriDiagMatrix
 
-import inspect
-
 from . import parameters as param
 from .analytic import analytic_ddm
 from .models.drift import DriftConstant, Drift
@@ -452,22 +450,18 @@ def simulated_solution(self, conditions={}, size=1000, rk4=True, seed=0):
     @returns(Boolean)
     def has_analytical_solution(self):
         """Is it possible to find an analytic solution for this model?"""
-        mt = self.get_model_type()
         # First check to make sure drift doesn't vary with time or
         # particle location
-        driftfuncsig = inspect.signature(mt["Drift"].get_drift)
-        if "t" in driftfuncsig.parameters or "x" in driftfuncsig.parameters:
+        if self.get_dependence("drift")._uses_t() or self.get_dependence("drift")._uses_x():
             return False
         # Check noise to make sure it doesn't vary with time or particle location
-        noisefuncsig = inspect.signature(mt["Noise"].get_noise)
-        if "t" in noisefuncsig.parameters or "x" in noisefuncsig.parameters:
+        if self.get_dependence("noise")._uses_t() or self.get_dependence("noise")._uses_x():
             return False
         # Check to make sure bound is one that we can solve for
-        boundfuncsig = inspect.signature(mt["Bound"].get_bound)
-        if "t" in boundfuncsig.parameters and mt["Bound"]!=BoundCollapsingLinear:
+        if self.get_dependence("bound")._uses_t() and self.get_dependence("bound").__class__ != BoundCollapsingLinear:
             return False
         # Make sure initial condition is a single point
-        if not issubclass(mt['IC'],(ICPointSourceCenter,ICPoint)):
+        if not issubclass(self.get_dependence("IC").__class__,(ICPointSourceCenter,ICPoint)):
             return False
         # Assuming none of these is the case, return True.
         return True
@@ -491,8 +485,7 @@ def can_solve_cn(self, conditions={}):
         # dis.get_instructions(get_bound)]" to see if it is used in the
         # function rather than looking to see if it is passed to the
         # function.
-        boundfuncsig = inspect.signature(self.get_dependence("bound").get_bound)
-        if "t" in boundfuncsig.parameters:
+        if self.get_dependence("bound")._uses_t():
             return False
         return True
 
@@ -515,10 +508,13 @@ def solve(self, conditions={}, return_evolution=False, force_python=False):
         self.check_conditions_satisfied(conditions)
         if self.has_analytical_solution() and return_evolution is False:
             return self.solve_analytical(conditions=conditions)
-        elif isinstance(self.get_dependence("bound"), BoundConstant) and return_evolution is False:
+        elif isinstance(self.get_dependence("bound"), BoundConstant) and return_evolution is False and (force_python or not HAS_CSOLVE):
             return self.solve_numerical_cn(conditions=conditions)
         else:
-            return self.solve_numerical_implicit(conditions=conditions, return_evolution=return_evolution)
+            if force_python or not HAS_CSOLVE or return_evolution:
+                return self.solve_numerical_implicit(conditions=conditions, return_evolution=return_evolution)
+            else:
+                return self.solve_numerical_c(conditions=conditions)
 
     @accepts(Self, conditions=Conditions, force_python=Boolean)
     @returns(Solution)
@@ -567,7 +563,9 @@ def solve_analytical(self, conditions={}, force_python=False):
         anal_pdf_err[0] = 0.
 
         # Fix numerical errors
-        pdfsum = (np.sum(anal_pdf_corr) + np.sum(anal_pdf_err))*self.dt
+        anal_pdf_corr *= self.dt
+        anal_pdf_err *= self.dt
+        pdfsum = np.sum(anal_pdf_corr) + np.sum(anal_pdf_err)
         if pdfsum > 1:
             if pdfsum > 1.01 and param.renorm_warnings:
                 print("Warning: renormalizing probability density from", pdfsum, "to 1.  " \
@@ -576,59 +574,73 @@ def solve_analytical(self, conditions={}, force_python=False):
             anal_pdf_corr /= pdfsum
             anal_pdf_err /= pdfsum
 
-        sol = Solution(anal_pdf_corr*self.dt, anal_pdf_err*self.dt, self, conditions=conditions)
+        sol = Solution(anal_pdf_corr, anal_pdf_err, self, conditions=conditions)
         return self.get_dependence('overlay').apply(sol)
 
     def solve_numerical_c(self, conditions={}):
-        mt = self.get_model_type()
-        driftfuncsig = inspect.signature(mt["Drift"].get_drift)
+        """Solve the DDM model using the implicit method with C extensions.
+
+        This function should give near identical results to
+        solve_numerical_implicit.  However, it uses compiled C code instead of
+        Python code to do so, which should make it much (10-100x) faster.
+
+        This does not current work with non-Gaussian diffusion matrices (a
+        currently undocumented feature).
+        """
+
         get_drift = self.get_dependence("drift").get_drift
-        if "t" not in driftfuncsig.parameters and "x" not in driftfuncsig.parameters:
+        drift_uses_t = self.get_dependence("drift")._uses_t()
+        drift_uses_x = self.get_dependence("drift")._uses_x()
+        if not drift_uses_t and not drift_uses_x:
             drifttype = 0
             drift = np.asarray([get_drift(conditions=conditions)])
-        elif "t" in driftfuncsig.parameters and "x" not in driftfuncsig.parameters:
+        elif drift_uses_t and not drift_uses_x:
             drifttype = 1
             drift = np.asarray([get_drift(t=t, conditions=conditions) for t in self.t_domain()])
-        elif "t" not in driftfuncsig.parameters and "x" in driftfuncsig.parameters:
+        elif not drift_uses_t and drift_uses_x:
             drifttype = 2
             drift = np.asarray(get_drift(x=self.x_domain(), conditions=conditions))
-        elif "t" in driftfuncsig.parameters and "x" in driftfuncsig.parameters:
+        elif drift_uses_t and drift_uses_x:
             drifttype = 3
             drift = np.concatenate([get_drift(t=t, x=self.x_domain(t=t, conditions=conditions), conditions=conditions) for t in self.t_domain()])
-        noisefuncsig = inspect.signature(mt["Noise"].get_noise)
-        get_noise = self.get_dependence("Noise").get_noise
-        if "t" not in noisefuncsig.parameters and "x" not in noisefuncsig.parameters:
+        get_noise = self.get_dependence("noise").get_noise
+        noise_uses_t = self.get_dependence("noise")._uses_t()
+        noise_uses_x = self.get_dependence("noise")._uses_x()
+        if not noise_uses_t and not noise_uses_x:
             noisetype = 0
             noise = np.asarray([get_noise(conditions=conditions)])
-        elif "t" in noisefuncsig.parameters and "x" not in noisefuncsig.parameters:
+        elif noise_uses_t and not noise_uses_x:
             noisetype = 1
             noise = np.asarray([get_noise(t=t, conditions=conditions) for t in self.t_domain()])
-        elif "t" not in noisefuncsig.parameters and "x" in noisefuncsig.parameters:
+        elif not noise_uses_t and noise_uses_x:
             noisetype = 2
             noise = np.asarray(get_noise(x=self.x_domain(), conditions=conditions))
-        elif "t" in noisefuncsig.parameters and "x" in noisefuncsig.parameters:
+        elif noise_uses_t and noise_uses_x:
             noisetype = 3
             noise = np.concatenate([get_noise(t=t, x=self.x_domain(conditions=conditions, t=t), conditions=conditions) for t in self.t_domain()])
-        boundfuncsig = inspect.signature(mt["Bound"].get_bound)
-        if "t" not in boundfuncsig.parameters:
+        bound_uses_t = self.get_dependence("bound")._uses_t()
+        if not bound_uses_t:
             boundtype = 0
             bound = np.asarray([self.get_dependence("Bound").get_bound(conditions=conditions)])
-        elif "t" in boundfuncsig.parameters:
+        elif bound_uses_t:
             boundtype = 1
             bound = np.asarray([self.get_dependence("Bound").get_bound(t=t, conditions=conditions) for t in self.t_domain()])
-        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)
+        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()))
         # TODO: Handle the pdf going below zero, returning pdfcurr, and fix numerical errors
         corr = (res[0]*self.dt)
         corr[corr<0] = 0
         err = (res[1]*self.dt)
         err[err<0] = 0
-        return self.get_dependence('overlay').apply(Solution(corr, err, self, conditions=conditions))
+        undec = res[2]
+        undec[undec<0] = 0
+        return self.get_dependence('overlay').apply(Solution(corr, err, self, conditions=conditions, pdf_undec=undec))
 
     @accepts(Self, method=Set(["explicit", "implicit", "cn"]), conditions=Conditions, return_evolution=Boolean)
     @returns(Solution)
     @requires("method == 'explicit' --> self.can_solve_explicit(conditions=conditions)")
     @requires("method == 'cn' --> self.can_solve_cn()")
     def solve_numerical(self, method="cn", conditions={}, return_evolution=False):
+
         """Solve the DDM model numerically.
 
         Use `method` to solve the DDM.  `method` can either be

ddm/models/__init__.py

@@ -13,7 +13,7 @@
                "OverlaySimplePause", "OverlayBlurredPause",
                "OverlayUniformMixture", "OverlayPoissonMixture",
                "OverlayNonDecision", "OverlayNonDecisionUniform", "OverlayNonDecisionGamma",
-           "LossFunction", "LossSquaredError", "LossLikelihood", "LossBIC"]
+           "LossFunction", "LossSquaredError", "LossLikelihood", "LossBIC", "LossRobustLikelihood", "LossRobustBIC"]
 
 from .base import *
 from .drift import *

ddm/models/base.py

@@ -7,6 +7,8 @@
 __all__ = ["Dependence"]
 
 import paranoid
+import dis
+import inspect
 
 class Dependence(object): # TODO Base this on ABC
     """An abstract class describing how one variable depends on other variables.
@@ -122,5 +124,26 @@ def __str__(self):
         return self.__repr__()
     def __hash__(self):
         return hash(repr(self))
+    def _uses(self, f, name):
+        """Check if function `f` uses a variable named `name`."""
+        # First get rid of wrappers, eh hem, paranoid
+        while "__wrapped__" in f.__dict__:
+            f = f.__wrapped__
+        # Get the names of all variables used in the function.  If it uses t,
+        # then return True.
+        vars_in_func = [inst.argrepr for inst in dis.get_instructions(f)]
+        if name in vars_in_func:
+            return True
+        # Check for the use of varargs or kwargs to be 100% safe.  Users
+        # shouldn't use these anyway, so if they do, then too bad, their
+        # function will run a bit slower.  Sucks to be them.
+        args = inspect.getargs(f.__code__).varargs
+        kwargs = inspect.getargs(f.__code__).varkw
+        if args and args in vars_in_func:
+            return True
+        if kwargs and kwargs in vars_in_func:
+            return True
+        return False
+
 
 

ddm/models/bound.py

@@ -23,6 +23,8 @@ class Bound(Dependence):
     Dependence.)
     """
     depname = "Bound"
+    def _uses_t(self):
+        return self._uses(self.get_bound, "t")
     def get_bound(self, t, conditions, **kwargs):
         """Calculate the bounds which particles cross to determine response time.
 

ddm/models/drift.py

@@ -23,6 +23,10 @@ class Drift(Dependence):
     and `required_parameters` (see documentation for Dependence.)
     """
     depname = "Drift"
+    def _uses_t(self):
+        return self._uses(self.get_drift, "t")
+    def _uses_x(self):
+        return self._uses(self.get_drift, "x")
     @accepts(Self, x=NDArray(d=1, t=Number), t=Positive0, dx=Positive, dt=Positive, conditions=Conditions)
     @returns(TriDiagMatrix)
     @ensures("return.shape == (len(x), len(x))")

ddm/models/noise.py

@@ -23,6 +23,10 @@ class Noise(Dependence):
     and `required_parameters` (see documentation for Dependence.)
     """
     depname = "Noise"
+    def _uses_t(self):
+        return self._uses(self.get_noise, "t")
+    def _uses_x(self):
+        return self._uses(self.get_noise, "x")
     @accepts(Self, x=NDArray(d=1, t=Number), t=Positive0, dx=Positive, dt=Positive, conditions=Conditions)
     @returns(TriDiagMatrix)
     @ensures("return.shape == (len(x), len(x))")