Binding polynomial extension (#46)

* update

* util changes for gui plot figs, starting binding_polynomial c extension

* binding polynomial extension

* organization

* updates

* more changes to binding polynomial

* binding polynomial C extension

Author: hrmyd

.gitignore

@@ -3,6 +3,8 @@ __pycache__/
 *.py[cod]
 *$py.class
 
+.DS_Store
+
 # C extensions
 *.so
 

pytc/experiments/base.py

@@ -17,16 +17,16 @@
 import random, string, os
 import numpy as np
 
-AVAIL_UNITS = {"cal":1.9872036,
-               "kcal":0.0019872036,
-               "J":8.3144598,
-               "kJ":0.0083144598}
-
 class BaseITCExperiment:
     """
     Class that holds an experimental ITC measurement and a model that describes it.
     """
 
+    AVAIL_UNITS = {"cal":1.9872036,
+               "kcal":0.0019872036,
+               "J":8.3144598,
+               "kJ":0.0083144598}
+
     def __init__(self,dh_file,model,shot_start=1,units="cal",uncertainty=0.1,
                  **model_kwargs):
         """
@@ -58,10 +58,10 @@ def __init__(self,dh_file,model,shot_start=1,units="cal",uncertainty=0.1,
         # Deal with units
         self._units = units
         try:
-            self._R = AVAIL_UNITS[self._units]
+            self._R = self.AVAIL_UNITS[self._units]
         except KeyError:
             err = "units must be one of:\n"
-            for k in AVAIL_UNITS.keys():
+            for k in self.AVAIL_UNITS.keys():
                 err += "    {}\n".format(k)
             err += "\n"
 

pytc/fitters/base.py

@@ -172,6 +172,8 @@ def corner_plot(self,filter_params=(),*args,**kwargs):
 
         return fig
 
+        return fig
+
     @property
     def samples(self):
         """

pytc/indiv_models/binding_polynomial.py

@@ -12,6 +12,8 @@
 import scipy.optimize
 from .base import ITCModel
 
+from . import bp_ext
+
 class BindingPolynomial(ITCModel):
     """
     Base class for a binding polynomial fit.
@@ -67,26 +69,6 @@ def _initialize_param(self):
         self._fit_dH_list    = ["dH{}".format(i+1) for i in range(self._num_sites)]
 
         self._T_conc_free = np.zeros(len(self._S_conc),dtype=float)
-        self._numerator = np.zeros(len(self._S_conc),dtype=float)
-        self._denominator = np.ones(len(self._S_conc),dtype=float)
-        self._i_array = np.arange(len(self._fit_beta_list),dtype=float) + 1
-
-    def _dQdT(self,T_free,S_total,T_total):
-        """
-        T_total = T_free + S_total*(dln(P)/dln(T_free)), so:
-        0 = T_free + S_total*(dln(P)/dln(T_free)) - T_total
-
-        Solve this for T_free to find free titrant.
-
-        dln(P)/dln(T_free)
-
-        P = (beta1*T_free**1) *  b2*T**2
-        """
-
-        numerator =       np.sum(self._i_array*self._fit_beta_array*(T_free**self._i_array))
-        denominator = 1 + np.sum(              self._fit_beta_array*(T_free**self._i_array))
-
-        return T_free + S_total*(numerator/denominator) - T_total
 
     @property
     def dQ(self):
@@ -102,66 +84,13 @@ def dQ(self):
             self._fit_dH_array[i] = self.param_values[self._fit_dH_list[i]]
 
         S_conc_corr = self._S_conc*self.param_values["fx_competent"]
+        num_shots = len(S_conc_corr)
+        size_T = self._T_conc.size
 
-        # Find the root of the derivative of the binding polynomial, giving the
-        # free titrant concentration
-        for i in range(len(S_conc_corr)):
-
-            # If there's no titrant, nothing is free.  (avoid numerical problems)
-            if self._T_conc[i] == 0:
-                self._T_conc_free[i] = 0.0
-                continue
-
-            # Manually check that bounds for root calculation have opposite signs
-            min_value = self._dQdT(            0.0,S_conc_corr[i],self._T_conc[i])
-            max_value = self._dQdT(self._T_conc[i],S_conc_corr[i],self._T_conc[i])
-
-            # Uh oh, they have same sign (root optimizer will choke)
-            if min_value*max_value > 0:
-
-                if max_value < 0:
-                    # root is closest to min --> set to that
-                    if (max_value < min_value):
-                        self._T_conc_free[i] = 0.0
-
-                    # root is closest to max --> set to that
-                    else:
-                        self._T_conc_free[i] = self._T_conc[i]
-                else:
-                    # root is closest to max --> set to that
-                    if (max_value < min_value):
-                        self._T_conc_free[i] = self._T_conc[i]
-
-                    # root is closest to min --> set to that
-                    else:
-                        self._T_conc_free[i] = 0.0
-
-                continue
-
-            T = scipy.optimize.brentq(self._dQdT,
-                                      0,self._T_conc[-1],
-                                      args=(S_conc_corr[i],
-                                           self._T_conc[i]))
-
-            # numerical problems sometimes make T slightly bigger than the total
-            # concentration, so bring down to the correct value
-            if (T > self._T_conc[i]):
-                T = self._T_conc[i]
-
-            self._T_conc_free[i] = T
-
-
-        # calculate the average enthalpy change
-        self._numerator   = 0.0
-        self._denominator = 1.0
-        for i in range(len(self._fit_beta_array)):
-            self._numerator   += self._fit_dH_array[i]*self._fit_beta_array[i]*(self._T_conc_free**(i+1))
-            self._denominator +=                       self._fit_beta_array[i]*(self._T_conc_free**(i+1))
-
-        avg_dH = self._numerator/self._denominator
-
-        X = avg_dH[1:] - avg_dH[:-1]
+        final_array = np.zeros((num_shots-1),dtype=float)
 
-        to_return = self._cell_volume*S_conc_corr[1:]*X + self.dilution_heats
+        bp_ext.dQ(self._cell_volume, num_shots, size_T, self._num_sites, self.dilution_heats, 
+            self._fit_beta_array, self._fit_dH_array, S_conc_corr, self._T_conc, self._T_conc_free, 
+            final_array)
 
-        return to_return
+        return final_array

pytc/util/util.py

@@ -35,11 +35,12 @@ def compare_models(*models):
     **NOTE**: This comparison is only valid if the models all fit the  same
     set of observations.  The models do not need to be nested.
     """
-
     aic = []
     aic_c = []
     bic = []
 
+    plots = []
+
     num_obs = None
     for i, m in enumerate(models):
 
@@ -56,14 +57,16 @@ def compare_models(*models):
 
         aic.append(m.fit_stats["AIC"])       
         aic_c.append(m.fit_stats["AICc"])       
-        bic.append(m.fit_stats["BIC"])       
+        bic.append(m.fit_stats["BIC"])
+
+        plots.append([m.plot(), i])       
 
     out = {}
 
     out["AIC"]  = weight_stat(aic)
     out["AICc"] = weight_stat(aic_c)
     out["BIC"]  = weight_stat(bic)
 
-    return out 
+    return out, plots
 
 

setup.py

@@ -7,6 +7,13 @@
 
 # Try using setuptools first, if it's installed
 from setuptools import setup, find_packages
+from setuptools.extension import Extension
+import numpy.distutils.misc_util
+
+# set up binding polynomial C extension
+ext = Extension('pytc.indiv_models.bp_ext', 
+      ['src/binding_polynomial.c',
+       'src/_bp_ext.c'])
 
 # Need to add all dependencies to setup as we go!
 setup(name='pytc-fitter',
@@ -20,5 +27,6 @@
       download_url='https://github.com/harmslab/pytc/tarball/1.1.0',
       zip_safe=False,
       install_requires=["matplotlib","scipy","numpy","emcee","corner"],
-      classifiers=['Programming Language :: Python'])
-
+      classifiers=['Programming Language :: Python'],
+      ext_modules=[ext],
+      include_dirs=numpy.distutils.misc_util.get_numpy_include_dirs())

src/_bp_ext.c

@@ -0,0 +1,111 @@
+#include <Python.h>
+//#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+#include <numpy/arrayobject.h>
+#include "binding_polynomial.h"
+#include <stdio.h>
+
+// docstrings
+static char module_docstring[] = 
+    "calculate binding polynomial fit";
+
+static char dQ_docstring[] = 
+    "Calculate the heats that would be observed across shots for a given set of enthalpies and binding constants for each reaction. This will work for an arbitrary-order binding polynomial.";
+
+static PyObject *bp_ext_dQ(PyObject *self, PyObject *args);
+
+// methods
+static PyMethodDef module_methods[] = {
+    {"dQ", bp_ext_dQ, METH_VARARGS, dQ_docstring},
+    {NULL, NULL}
+};
+
+// init function
+PyMODINIT_FUNC PyInit_bp_ext(void)
+{
+    PyObject *module;
+    static struct PyModuleDef moduledef = {
+        PyModuleDef_HEAD_INIT,
+        "bp_ext",
+        module_docstring,
+        -1,
+        module_methods,
+        NULL,
+        NULL,
+        NULL,
+        NULL
+    };
+
+    module = PyModule_Create(&moduledef);
+    if (!module) return NULL;
+
+    // Load numpy funcionality.
+    import_array();
+
+    return module;
+}
+
+static PyObject *bp_ext_dQ(PyObject *self, PyObject *args)
+{   
+    int num_sites, num_shots, size_T_conc;
+    double cell_volume;
+    PyObject *fit_beta_obj, *fit_dH_obj, *S_conc_corr, *T_conc, *T_conc_free, *dilution_heats, *final_array;
+
+    if (!PyArg_ParseTuple(args, "diiiOOOOOOO:dQ", &cell_volume, &num_shots, &size_T_conc, &num_sites, 
+                                        &dilution_heats, &fit_beta_obj, &fit_dH_obj, &S_conc_corr, &T_conc, 
+                                        &T_conc_free, &final_array)){
+        return NULL;
+    }
+
+    // PyObjects to numpy arrays
+    PyObject *fit_beta_array = PyArray_FROM_OTF(fit_beta_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *fit_dH_array = PyArray_FROM_OTF(fit_dH_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *T_conc_array = PyArray_FROM_OTF(T_conc, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *T_conc_free_array = PyArray_FROM_OTF(T_conc_free, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *S_conc_corr_array = PyArray_FROM_OTF(S_conc_corr, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *dilution_heats_array = PyArray_FROM_OTF(dilution_heats, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *final_array_obj = PyArray_FROM_OTF(final_array, NPY_DOUBLE, NPY_IN_ARRAY);
+
+    // throw exception if cast failed
+    if (fit_beta_array == NULL || fit_dH_array == NULL || T_conc_array == NULL || T_conc_free_array == NULL || 
+        S_conc_corr_array == NULL || dilution_heats_array == NULL || final_array_obj == NULL) {
+        Py_XDECREF(fit_beta_array);
+        Py_XDECREF(fit_dH_array);
+        Py_XDECREF(T_conc_array);
+        Py_XDECREF(T_conc_free_array);
+        Py_XDECREF(S_conc_corr_array);
+        Py_XDECREF(dilution_heats_array);
+        Py_XDECREF(final_array_obj);
+        return NULL;
+    }
+
+    // pointers 
+    double *fit_beta = (double*)PyArray_DATA(fit_beta_array);
+    double *fit_dH = (double*)PyArray_DATA(fit_dH_array);
+    double *T_conc_p = (double*)PyArray_DATA(T_conc_array);
+    double *T_conc_free_p = (double*)PyArray_DATA(T_conc_free_array);
+    double *S_conc_corr_p = (double*)PyArray_DATA(S_conc_corr_array);
+    double *dilution_heats_p = (double*)PyArray_DATA(dilution_heats_array);
+    double *final_array_p = (double*)PyArray_DATA(final_array_obj);
+
+    // call function 
+    dQ(fit_beta, fit_dH, S_conc_corr_p, T_conc_p, T_conc_free_p, cell_volume, dilution_heats_p, 
+        num_sites, num_shots, size_T_conc, final_array_p);
+
+    // clean up 
+    Py_DECREF(fit_beta_array);
+    Py_DECREF(fit_dH_array);
+    Py_DECREF(T_conc_array);
+    Py_DECREF(T_conc_free_array);
+    Py_DECREF(S_conc_corr_array);
+    Py_DECREF(dilution_heats_array);
+
+    // check errors
+    /*if (value == NULL){
+        PyErr_SetString(PyExc_RuntimeError, "object return null");
+
+        return NULL;
+    }*/
+
+    // build output and return
+    return Py_BuildValue("O", final_array_obj);
+}