diff --git a/itc_experiment.py b/itc_experiment.py index b4ad7a9..873968c 100644 --- a/itc_experiment.py +++ b/itc_experiment.py @@ -96,7 +96,7 @@ def __init__(self, T, V0, injections, dQ, Cell, Syringe, skip=[], ddQ=[], Q_dil= for s in self.Cell: self.Concentrations[i][s] += self.Cell[s] * ( (1-(dV/(2.0*V0))) / (1.0+(dV/(2.0*V0))) ) - # convert raw data (in calories) to joules per mol of injectant + # convert raw data (in calories) to joules (note that this is not normalized per mol of injectant!) assert len(dQ) == self.npoints self.dQ_exp = numpy.array([J_from_cal(dQ[i]) for i in xrange(self.npoints)],dtype='d') self.dQ_fit = None @@ -156,6 +156,7 @@ def make_plot(self,hardcopy=False,hardcopydir='.',hardcopyprefix='', hardcopytyp pyplot.ylabel("%s/mol of %s"%(self.units,self.syringeRef)) pyplot.xlabel("%s / %s"%(self.syringeRef,self.cellRef)) + # For convention, normalize the heat evolved as per mol of injected reference ligand tmpx = [ self.Concentrations[i][self.syringeRef]/self.Concentrations[i][self.cellRef] for i in xrange(self.npoints) if i not in self.skip ] tmpy = [ convert_from_J(self.units,self.dQ_exp[i])/self.Syringe[self.syringeRef]/self.injections[i] for i in xrange(self.npoints) if i not in self.skip ] tmpd = [ convert_from_J(self.units,self.ddQ[i]) for i in xrange(self.npoints) if i not in self.skip ] @@ -211,7 +212,7 @@ def export_to_file(self,path,units='cal',full=False): h.write("# skip %s\n"%(",".join(map(str,self.skip)))) if full: - h.write("#\n# Ivol dQ_exp ddQ_exp dQ_fit spline skipped %s %s\n"%("\t".join(self.Cell.keys()),"\t".join(self.Syringe.keys()))) + h.write("#\n# Ivol dQ_exp ddQ_exp dQ_fit dQ_spline skipped %s %s\n"%("\t".join(self.Cell.keys()),"\t".join(self.Syringe.keys()))) else: h.write("#\n# Ivol dQ_exp\n") @@ -221,13 +222,11 @@ def export_to_file(self,path,units='cal',full=False): spline = self.spline[:] for i in xrange(self.npoints): - cell = ["%.5f"%(self.Concentrations[i][s]) for s in self.Cell] - syringe = ["%.5f"%(self.Concentrations[i][s]) for s in self.Syringe] + cell = ["%.5E"%(self.Concentrations[i][s]) for s in self.Cell] + syringe = ["%.5E"%(self.Concentrations[i][s]) for s in self.Syringe] if full: - h.write("%.5f %.5f %.5E %.5E %.5f %.5E %.5E %i %s %s\n"%( + h.write("%.5f %.5E %.5E %.5E %.5E %i %s %s\n"%( self.injections[i], - self.Concentrations[i][self.cellRef], - self.Concentrations[i][self.syringeRef], convert_from_J(units,self.dQ_exp[i]), convert_from_J(units,self.ddQ[i]), convert_from_J(units,self.dQ_fit[i]), @@ -250,16 +249,19 @@ def get_chisq(self, Q, writeback=False): Returns: (float): The goodness of the fit, as a reduced chi-square. """ - + dV = 0.0 for i in xrange(self.npoints): dV += self.injections[i] Q[i] += self.Q_dil*self.injections[i]*(1.0 - 1.0/(1.0+self.V0/dV)) # add heat of dilution Q[i] *= self.V0 * self.Concentrations[i][self.cellRef] # normalize total heat content to macromolecule concentration in the cell volume - # obtain the change in cell heat b/t each titration point + # obtain the change in cell heat between each titration point dQ = [0.0]*self.npoints for i in xrange(self.npoints): + #if self.title == "35C-01-TRAPstk-TrpA": + # print Q[i] + if i==0: dQ[i] = Q[i] + ( (self.injections[i]/self.V0)*(Q[i]/2.0) ) else: diff --git a/itc_fit.py b/itc_fit.py index 1680959..ee7c32b 100644 --- a/itc_fit.py +++ b/itc_fit.py @@ -33,7 +33,7 @@ def __init__(self, sim, method='simplex', method_args={}, verbose=False): """ self.sim = sim - self.model = self.sim.model + self.model = self.sim.get_model() self.method = method self.method_args = method_args self.verbose = verbose @@ -41,7 +41,14 @@ def __init__(self, sim, method='simplex', method_args={}, verbose=False): # obtain model-defined boundaries to enforce during fitting self.bounds = dict( (name,self.model.get_param_bounds(name)) for name in self.model.get_param_names() ) - + + def set_sim(self, sim): + self.sim = sim + self.model = self.sim.get_model() + + def get_sim(self): + return self.sim + def add_bounds(self, param, low=None, high=None): """Add low and/or high boundaries for a specified parameter. @@ -93,7 +100,7 @@ def _target(x,sim): return sim.get_chisq() * (1+m) return sim.run(writeback=False) - + # optimize parameters opt = self._fitter( _target, x0, callback ) diff --git a/itc_model.py b/itc_model.py index 9821517..09c8019 100644 --- a/itc_model.py +++ b/itc_model.py @@ -11,7 +11,7 @@ class ITCModel(): def __init__(self): self.params = OrderedDict() self.components = [] - self.units = None + self.units = "J" self._param_meta, self._component_meta = {},{} diff --git a/itc_sim.py b/itc_sim.py index 231d47d..1feb6d5 100644 --- a/itc_sim.py +++ b/itc_sim.py @@ -82,7 +82,7 @@ def set_model_params(self, *args, **kwargs ): def set_model_param(self, param, value): self.model.set_param( param, value ) - + ### def info(self): print str(self) @@ -152,10 +152,13 @@ def run( self, experiments=None, writeback=True ): """ if experiments: for E in experiments: - self.in_Queue.put( (self.model.get_params(units=_units[0]),E) ) + self.in_Queue.put( (self.model.get_params(units=self.units),E) ) else: + if self.size == 0: + print "No experiments to simulate." + return None for E in self.get_experiments(): - self.in_Queue.put( (self.model.get_params(units=_units[0]),E) ) + self.in_Queue.put( (self.model.get_params(units=self.units),E) ) queue_contents = [] while len(queue_contents) < self.size: diff --git a/ms_experiment.py b/ms_experiment.py new file mode 100644 index 0000000..5e692f4 --- /dev/null +++ b/ms_experiment.py @@ -0,0 +1,119 @@ +import os +import numpy as np + +_MATPLOTLIB_BACKEND = None #None for default + +class MSExperiment(): + def __init__(self, file): + self.title = os.path.splitext(os.path.basename(file))[0] + self.chisq = None + + fh = open(file) + self.T = float(fh.readline().split()[2]) + PConc = map(float,fh.readline().split()[2:]) + LConc = map(float,fh.readline().split()[2:]) + fh.close() + + self.npoints = len(PConc) + assert len(LConc) == self.npoints + + data = np.genfromtxt( file, usecols=xrange(1,self.npoints+1), unpack=True ) + assert len(data) == self.npoints + self.npops = len(data[0])/2 + for i in xrange(len(data)): + assert len(data[i])/2 == self.npops + + self.Concentrations,self.PopIntens,self.PopErrors,self.PopFits = [],[],[],[] + for i in xrange(self.npoints): + self.Concentrations.append({}) + self.Concentrations[i]['Lattice'] = PConc[i] / 1.0E6 + self.Concentrations[i]['Ligand'] = LConc[i] / 1.0E6 + self.PopIntens.append( data[i][:self.npops] ) + self.PopErrors.append( data[i][self.npops:] ) + self.PopFits.append( [0.0]*self.npops ) + + def make_plot(self,hardcopy=False,hardcopydir='.',hardcopyprefix='', hardcopytype='png'): + try: + if _MATPLOTLIB_BACKEND != None: + import matplotlib + matplotlib.use(_MATPLOTLIB_BACKEND) + import matplotlib.pyplot as pyplot + except: + pyplot = None + + if pyplot == None: return + if hardcopy: fig = pyplot.figure() + + pyplot.clf() + pyplot.title(self.title) + pyplot.ylabel("Experimental - Fit Abundance") + pyplot.xlabel("Mass Populations") + + ax1 = pyplot.subplot(3,1,1) + ax2 = pyplot.subplot(3,1,2) + ax3 = pyplot.subplot(3,1,3) + + xax_positions,xax_labels = [],[] + width,space,left = 0.25,0.5,0.0 + for i in xrange(self.npops): + bars1 = ax1.bar( [left + (j*width) for j in xrange(self.npoints)], [self.PopIntens[j][i] for j in xrange(self.npoints)], width=width, edgecolor='r' ) + bars2 = ax2.bar( [left + (j*width) for j in xrange(self.npoints)], [self.PopFits[j][i] for j in xrange(self.npoints)], width=width, edgecolor='r' ) + bars3 = ax3.bar( [left + (j*width) for j in xrange(self.npoints)], [self.PopIntens[j][i] - self.PopFits[j][i] for j in xrange(self.npoints)], width=width, edgecolor='r' ) + xax_positions.append( left + (self.npoints*width)/2.0 ) + xax_labels.append("%i"%i) + left += (j*width)+space + + for j in xrange(self.npoints): + #color = float(j) / self.npoints + #bars1[j].set_color( (0, 0, color) ) + #bars2[j].set_color( (0, color, 0) ) + #bars3[j].set_color( (color, 0, 0) ) + color = 1.0 - (float(j) / self.npoints) + bars1[j].set_color( (color, color, 1) ) + bars2[j].set_color( (color, 1, color) ) + bars3[j].set_color( (1, color, color) ) + + bars1[j].set_edgecolor( (0,0,0) ) + bars2[j].set_edgecolor( (0,0,0) ) + bars3[j].set_edgecolor( (0,0,0) ) + + ax1.set_ylabel("Experimental") + ax1.set_xticks( xax_positions ) + ax1.set_xticklabels( xax_labels) + ax2.set_ylabel("Fit") + ax2.set_xticks( xax_positions ) + ax2.set_xticklabels( xax_labels) + ax3.set_ylabel("Residuals") + ax3.set_xticks( xax_positions ) + ax3.set_xticklabels( xax_labels) + + pyplot.draw() + + if hardcopy: + fig.savefig( os.path.join(hardcopydir,"%s%s.%s"%(hardcopyprefix,self.title,hardcopytype)), bbox_inches='tight') + pyplot.close(fig) + else: + pyplot.show() + + def export_to_file(self,path): + fh = open(path, 'w') + fh.write("# Experimental data\n") + for i in xrange(self.npoints): + fh.write("%i %s\n"%(i,"\t".join(["%f"%f for f in self.PopIntens[i]]))) + fh.write("# Fit data\n") + for i in xrange(self.npoints): + fh.write("%i %s\n"%(i,"\t".join(["%f"%f for f in self.PopFits[i]]))) + fh.close() + + def get_chisq(self, data, writeback=False): + self.chisq = 0.0 + for i in xrange(self.npoints): + for j in xrange(self.npops): + self.chisq += (data[i][j] - self.PopIntens[i][j])**2 / self.PopErrors[i][j]**2 + + if writeback: + self.PopFits = data + + self.chisq = self.chisq / (self.npoints * self.npops) + + return self.chisq \ No newline at end of file