Multiple improvements to pymcmodels

.travis.yml

@@ -30,6 +30,8 @@ script:
   - conda install --yes --quiet nose nose-timer
   # Test the package
   - cd devtools && nosetests $PACKAGENAME --nocapture --verbosity=2 --with-timer -a '!slow' && cd ..
+  # Run quickmodel
+  - pushd . && cd examples/direct-fluorescence-assay && env PYTHONPATH="./" quickmodel --inputs 'inputs_p38_singlet' && popd
   # Run IPython notebook tests
   # THIS NEEDS TO BE REPLACED WITH A PY3.x COMPATIBLE SCHEME
   #- source devtools/travis-ci/ipythontests.sh

AssayTools/analysis.py

@@ -777,6 +777,10 @@ def run_mcmc(self, dbfilename='output'):
 
         """
 
+        # DEBUG: Write model
+        print('Writing graph...')
+        pymc.graph.moral_graph(self.model, format='ps')
+
         # Sample the model with pymc
         # TODO: Allow
         mcmc = pymc.MCMC(self.model, db='sqlite', name='output', verbose=True)

AssayTools/bindingmodels.py

@@ -72,21 +72,32 @@ def equilibrium_concentrations(cls, DeltaG, Ptot, Ltot):
          Bound complex concentration, molarity.
 
       """
+      # Handle only strictly positive elements---all others are set to zero as constants
+      try:
+          nonzero_indices = np.where(Ltot > 0)[0]
+          zero_indices = np.where(Ltot <= 0)[0]
+      except:
+          nonzero_indices = range(size[0])
+          zero_indices = []
+      nnonzero = len(nonzero_indices)
+      nzeros = len(zero_indices)
 
       # Original form:
       #Kd = np.exp(DeltaG)
       #sqrt_arg = (Ptot + Ltot + Kd)**2 - 4*Ptot*Ltot
       #sqrt_arg[sqrt_arg < 0.0] = 0.0
       #PL = 0.5 * ((Ptot + Ltot + Kd) - np.sqrt(sqrt_arg));  # complex concentration (M)
 
+
       # Numerically stable variant?
-      logP = np.log(Ptot)
-      logL = np.log(Ltot)
+      PL = np.zeros(Ptot.shape)
+      logP = np.log(Ptot[nonzero_indices])
+      logL = np.log(Ltot[nonzero_indices])
       logPLK = np.logaddexp(np.logaddexp(logP, logL), DeltaG)
       PLK = np.exp(logPLK);
       sqrt_arg = 1.0 - np.exp(np.log(4.0) + logP + logL - 2*logPLK);
       sqrt_arg[sqrt_arg < 0.0] = 0.0 # ensure always positive
-      PL = 0.5 * PLK * (1.0 - np.sqrt(sqrt_arg));  # complex concentration (M)
+      PL[nonzero_indices] = 0.5 * PLK * (1.0 - np.sqrt(sqrt_arg));  # complex concentration (M)
 
       # Another variant
       #PL = 2*Ptot*Ltot / ((Ptot+Ltot+Kd) + np.sqrt((Ptot + Ltot + Kd)**2 - 4*Ptot*Ltot));  # complex concentration (M)

AssayTools/platereader.py

@@ -164,7 +164,8 @@ def select_data(file_data,section_name,selection,*args,**kwargs):
 
     # extract wavelength (only relevant for spectra)
     # if you don't include wavelength for spectra, data will be a dict of all wavelengths
-    if type(data.itervalues().next()) == dict and wavelength != None:
+    datatypes = [ type(entry) for entry in data.values() ]
+    if datatypes[0] == dict and wavelength != None:
         new_data = dict()
         for key in data.keys():
             new_data[key] = data[key][wavelength]

docs/theory.rst

@@ -199,6 +199,10 @@ In the :mod:`pymc` model, these priors are implemented via ::
   model['F_plate'] = pymc.Uniform('F_plate', lower=0.0, upper=Fmax, value=F_plate_guess)
   model['F_buffer'] = pymc.Uniform('F_buffer', lower=0.0, upper=Fmax/path_length, value=F_buffer_guess)
 
+If an estimate of :math:`F_{PL}` is known from a prior experiment, this value and its standard error can be specified via a lognormal distribution ::
+
+  model['F_PL'] = pymc.Lognormal('F_PL', mu=np.log(F_PL**2 / np.sqrt(dF_PL**2 + F_PL**2)), tau=np.sqrt(np.log(1.0 + (dF_PL/F_PL)**2))**(-2))
+
 Top/bottom detector gain.
 """""""""""""""""""""""""
 .. _detector-gain:

examples/direct-fluorescence-assay/inputs_Gef_WIP.py

@@ -8,7 +8,7 @@
     'ligand_order'  :  ['Gefitinib','Gefitinib','Gefitinib','Gefitinib'],
     'section'       :  '280_TopRead',
     'wavelength'    :  '480',
-    'Lstated'       :  np.array([20.0e-6,9.15e-6,4.18e-6,1.91e-6,0.875e-6,0.4e-6,0.183e-6,0.0837e-6,0.0383e-6,0.0175e-6,0.008e-6,0.0001e-6], np.float64), # ligand concentration
+    'Lstated'       :  np.array([20.0e-6,9.15e-6,4.18e-6,1.91e-6,0.875e-6,0.4e-6,0.183e-6,0.0837e-6,0.0383e-6,0.0175e-6,0.008e-6,0.0], np.float64), # ligand concentration
     'Pstated'       :  0.5e-6 * np.ones([12],np.float64), # protein concentration, M
     'assay_volume'  :  100e-6, # assay volume, L
     'well_area'     :  0.3969, # well area, cm^2 for 4ti-0203 [http://4ti.co.uk/files/3113/4217/2464/4ti-0201.pdf]

examples/direct-fluorescence-assay/inputs_p38_singlet.py

@@ -2,10 +2,10 @@
 import numpy as np
 
 inputs = {
-    'xml_file_path' :  "./data/singlet_copy/",
+    'xml_file_path' :  "./data/single_wavelength_copy",
     'section'       :  '280_480_TOP_120',
     'ligand_order'  :  ['Bosutinib','Bosutinib Isomer','Erlotinib','Gefitinib','Ponatinib','Lapatinib','Saracatinib','Vandetanib'],
-    'Lstated'       :  np.array([20.0e-6,14.0e-6,9.82e-6,6.88e-6,4.82e-6,3.38e-6,2.37e-6,1.66e-6,1.16e-6,0.815e-6,0.571e-6,0.4e-6,0.28e-6,0.196e-6,0.138e-6,0.0964e-6,0.0676e-6,0.0474e-6,0.0320e-6,0.0240e-6,0.0160e-6,0.0120e-6,0.008e-6,0.00001e-6], np.float64), # ligand concentration, M
+    'Lstated'       :  np.array([20.0e-6,14.0e-6,9.82e-6,6.88e-6,4.82e-6,3.38e-6,2.37e-6,1.66e-6,1.16e-6,0.815e-6,0.571e-6,0.4e-6,0.28e-6,0.196e-6,0.138e-6,0.0964e-6,0.0676e-6,0.0474e-6,0.0320e-6,0.0240e-6,0.0160e-6,0.0120e-6,0.008e-6,0.0], np.float64), # ligand concentration, M
     'protein'       :  'p38',
     'Pstated'       :  0.5e-6 * np.ones([24],np.float64), # protein concentration, M
     'assay_volume'  :  50e-6, # assay volume, L

scripts/quickmodel.py

@@ -14,6 +14,7 @@
 import string
 import json
 import numpy as np
+import traceback
 
 import seaborn as sns
 import pymbar
@@ -61,7 +62,7 @@ def quick_model(inputs):
             protein_row = ALPHABET[i]
             buffer_row = ALPHABET[i+1]
 
-            name = "%s-%s%s"%(inputs['ligand_order'][i/2],protein_row,buffer_row)
+            name = "%s-%s%s"%(inputs['ligand_order'][int(i/2)],protein_row,buffer_row)
 
             print(name)
 
@@ -71,7 +72,10 @@ def quick_model(inputs):
             try:
                 part1_data_protein = platereader.select_data(data, inputs['section'], protein_row)
                 part1_data_buffer = platereader.select_data(data, inputs['section'], buffer_row)
-            except:
+            except Exception as e:
+                # Print exception
+                print("An exception occurred while processing '%s' rows %s and %s:" % (my_file, protein_row, buffer_row))
+                print(traceback.print_exc())
                 continue
 
             reorder_protein = reorder2list(part1_data_protein,well)
@@ -111,7 +115,7 @@ def quick_model(inputs):
             [t,g,Neff_max] = pymbar.timeseries.detectEquilibration(mcmc.DeltaG.trace())
             DeltaG_equil = mcmc.DeltaG.trace()[t:].mean()
             dDeltaG_equil = mcmc.DeltaG.trace()[t:].std()
-            
+
             ## PLOT MODEL
             #from assaytools import plots
             #figure = plots.plot_measurements(Lstated, Pstated, pymc_model, mcmc=mcmc)
@@ -155,9 +159,9 @@ def quick_model(inputs):
                 clrs[idx] = (.5,.5,.5)
             for idx in gray_after:
                 clrs[idx] = (.5,.5,.5)
-           
+
             plt.subplot(312)
-            
+
             plt.bar(bin_edges[:-1],hist,binwidth,color=clrs, edgecolor = "white");
             sns.kdeplot(mcmc.DeltaG.trace()[t:],bw=.4,color=(0.39215686274509803, 0.7098039215686275, 0.803921568627451),shade=False)
             plt.axvline(x=interval[0],color=(0.5,0.5,0.5),linestyle='--')
@@ -167,8 +171,8 @@ def quick_model(inputs):
             plt.xlabel('$\Delta G$ ($k_B T$)',fontsize=16);
             plt.ylabel('$P(\Delta G)$',fontsize=16);
             plt.xlim(-35,-10)
-            hist_legend = mpatches.Patch(color=(0.7372549019607844, 0.5098039215686274, 0.7411764705882353), 
-                        label = '$\Delta G$ =  %.3g [%.3g,%.3g] $k_B T$' 
+            hist_legend = mpatches.Patch(color=(0.7372549019607844, 0.5098039215686274, 0.7411764705882353),
+                        label = '$\Delta G$ =  %.3g [%.3g,%.3g] $k_B T$'
                         %(interval[1],interval[0],interval[2]) )
             map_legend = mlines.Line2D([],[],color='black',label="MAP = %.1f $k_B T$"%DeltaG_map)
             plt.legend(handles=[hist_legend,map_legend],fontsize=14,loc=0,frameon=True);
@@ -251,7 +255,7 @@ def quick_model_spectra(inputs):
             protein_row = ALPHABET[i]
             buffer_row = ALPHABET[i+1]
 
-            name = "%s-%s-%s%s"%(protein,inputs['ligand_order'][i/2],protein_row,buffer_row)
+            name = "%s-%s-%s%s"%(protein,inputs['ligand_order'][int(i/2)],protein_row,buffer_row)
 
             print(name)