ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/fit_module.py
(Generate patch)
# Line 1 | Line 1
1   """
2 < fit: Curve-fitting module for SPRI data.
3 < Christopher Lausted, Institute for Systems Biology,
4 < OSPRAI developers
5 < Last modified on 100518 (yymmdd)
6 <
7 < Example:
8 < #import fit_module as fit
9 < #import mdl_module as mdl
10 < #ba1.roi[0].model = mdl.drift
11 < #ba1.roi[0].params = {'rate': {'value':1, 'min':-100.0, 'max':100.0, 'fixed':'float'} }
12 < #success = lma(ba1.roi[0])
13 < #for i in ba1.roi: lma(i)
2 > fit_module
3 > -------------------------------------------------------------------------------
4 >
5 > *Curve-fitting module for SPRI data.*
6 > These routines are based on the Levenberg-Marquart Algorithm (LMA) from SciPy.
7 > This module provides several extended features that make it
8 > particularly powerful for the analysis of high-throughput SPRI data.
9 >
10 >  * Constrained parameters.
11 >    *Invalid results such as negative kinetic rates can be avoided.*
12 >  * Simultaneous fitting of multiple curves.
13 >    *Parameters used in one ROI may be tied to those in other ROIs.*
14 >  * Multiple models may be specified.
15 >    *Different ROIs may exhibit different binding interactions.*
16 >
17 > Users should take advantage of the ``mclma()`` function.  
18 > The ``lma()`` and ``clma()`` functions have been deprecated.
19 >
20 > .. moduleauthor:: Christopher Lausted,
21 >                  Institute for Systems Biology,
22 >                  OSPRAI developers.
23 >                  
24 > Examples::
25 >
26 >  >>> import fit_module as fit
27 >  >>> import mdl_module as mdl
28 >  >>> ba.roi[0].model = mdl.drift
29 >  >>> ba.roi[1].model = mdl.drift
30 >  >>> p1 = {'rate': {'value':1.0, 'min':-100.0, 'max':100.0, 'fixed':'float'} }
31 >  >>> p2 = dict(rate=dict(value=1.0, min=-100.0, max=100.0, fixed=True))
32 >  >>> ba.roi[0].params = p1
33 >  >>> ba.roi[1].params = p2
34 >  >>> i = mclma(ba.roi)
35   """
36 < __version__ = "100518"
36 > __version__ = "110216"
37  
38  
39   ## Import libraries
# Line 25 | Line 46
46  
47   def lma(roi):
48      """
49 <    Normal Levenberg-Marquart fitting based on SciPy.
50 <    This function takes a single ba_class RegionOfInterest object (roi).
51 <    It modifies roi.params in place.  It returns the number of iterations completed.
49 >    *Normal LMA curve-fitting.  (Deprecated)*
50 >    Least-squares curve-fitting is performed on the data in the specified ROI.
51 >    The ROI must contain SPR data, a model, and initial (guess) parameters.
52 >    The roi.params values are modified in place.
53 >    The maximum number of iterations is set to 500.
54 >    
55 >    :param roi: ROI containing data, a model, and initial parameter guesses.
56 >    :type roi: RegOfInterest
57 >    :returns: number of iterations completed.
58      """
59      params = roi.params
60 <    checkparams(params)
60 >    __checkparams(params)
61      ## Create list initial values for the floating parameters.
62      p0 = [dict['value'] for key,dict in params.iteritems() if (dict['fixed'] != True)]
63      ## Error function.
64 <    erf = lambda p: (model_interface(roi, p) - roi.value)
64 >    erf = lambda p: (__model_interface(roi, p) - roi.value)
65      ## Fit.
66      #p1, success = leastsq(erf, p0, args=(), maxfev=999)
67      p1, success = leastsq(erf, p0, maxfev=500)
68      return success
69 <    ## End of constlma() function.
69 >    """End of constlma() function."""
70      
71      
72 < def model_interface(roi, pfloat):
72 > def __model_interface(roi, pfloat):
73      """
74      Use the model assigned to this roi to simulation data vs time.  
75      Provide a list of values for the floating parameters (pfloat).
# Line 53 | Line 80
80          if (roi.params[key]['fixed'] != True):
81              roi.params[key]['value'] = pfloat[i]
82      return roi.model(roi.time, roi.value, roi.params)
83 +    """End of __model_interface() function."""
84      
85      
86 < def checkparams(params):
86 > def __checkparams(params):
87      """
88      Check that each dictionary in the params dictionary contains the four
89 <    keys (value, min, max, fixed) and add them if necessary. Default is 'fixed':'fixed'.
89 >    keys (value, min, max, fixed) and add them if necessary.
90 >    Default is 'fixed':'fixed'.
91      Also check that when parameters float, min<value<max.    
92      """
93      flag = False
# Line 81 | Line 110
110      if (flag == True):
111          print "Parameter min/max errors: Modifications were made automatically."
112      return
113 <    ## End of checkparams() function.
113 >    """End of __checkparams() function."""
114  
115  
116   #### The code below is a variant allowing parameter constraints ####
117 < def transform(a, b, x0):
117 > def __transform(a, b, x0):
118      """
119      A transformation to convert x from range (-inf,+inf) to (a,b).
120      This is used to help constrain the outputs of the Levenberg-Marquart algorithm.
# Line 96 | Line 125
125      #x = max(-x0, -709)  ## Function math.exp(>709) overflows.
126      #x1 = a + ( (b-a) / (1+exp(-x) ))
127      ## Try transform: y = ((b-a)*tanh(x)+b+a)/2
128 <    x1 = ( (b-a)*tanh(x0) + b + a ) / 2
128 >    x1 = ( (b-a)*tanh(x0) + b + a ) * 0.5
129      return x1
130 <    ## End of transform() function.
130 >    """End of __transform() function."""
131      
132      
133 < def itransform(a, b, x1):
133 > def __itransform(a, b, x1):
134      """
135      The (inverse) transformation to convert x from range(a,b) to (-inf,+inf).
136      This is used to help constrain the outputs of the Levenberg-Marquart algorithm.
# Line 109 | Line 138
138      #x0 = -log( ((b-a)/(x1-a)) - 1 )
139      x0 = arctanh( (2*x1-b-a) / (b-a) )
140      return x0
141 <    ## End of itransform() function.
141 >    """End of __itransform() function."""
142  
143  
144   def clma(roi):
145      """
146 <    Constrained Levenberg-Marquart Algorithm fitting based on SciPy.
147 <    This function takes a single ba_class RegionOfInterest object (roi).
148 <    It modifies roi.params in place.  It returns the number of iterations completed.
146 >    *Constrained LMA curve-fitting.  (Deprecated)*
147 >    Least-squares curve-fitting is performed on the data in the specified ROI.
148 >    The ROI must contain SPR data, a model, and initial (guess) parameters.
149 >    The roi.params values are modified in place.
150 >    The maximum number of iterations is set to 999.
151 >    
152 >    :param roi: ROI containing data, a model, and initial parameter guesses.
153 >    :type roi: RegOfInterest
154 >    :returns: number of iterations completed.
155      """
156      params = roi.params
157 <    checkparams(params)
157 >    __checkparams(params)
158      
159      ## Create list of initial values for the floating parameters.
160      p0 = [dict['value'] for key,dict in params.iteritems() if (dict['fixed'] != True)]
# Line 127 | Line 162
162      for i,key in enumerate(params):
163          if (params[key]['fixed'] != True):
164              a, b = params[key]['min'], params[key]['max']
165 <            p0[i] = itransform(a, b, float(p0[i]))
166 <            # Temp
167 <            print key, params[key]['value'], p0[i]
165 >            p0[i] = __itransform(a, b, float(p0[i]))
166 >            #Next line for troublshooting.
167 >            #print key, params[key]['value'], p0[i]
168      
169      ## Fit.
170 <    p1, success = leastsq(cerf, p0, args=(roi), maxfev=999)
170 >    p1, success = leastsq(__cerf, p0, args=(roi), maxfev=999)
171      return success
172 <    ## End of constrained lma() function.
172 >    """End of constrained lma() function."""
173      
174      
175 < def cerf(pfloat, roi):
175 > def __cerf(pfloat, roi):
176      """
177      Error estimation function for use with constrained LMA.
178      Use the model assigned to this roi to simulation data vs time.  
# Line 150 | Line 185
185      for i,key in enumerate(params):
186          if (params[key]['fixed'] != True):
187              a, b = params[key]['min'], params[key]['max']
188 <            params[key]['value'] = transform(a, b, pfloat[i])
188 >            params[key]['value'] = __transform(a, b, pfloat[i])
189              print "%s %.2e" % (key, params[key]['value']), # Temp.
190      erf = (roi.model(roi.time, roi.value, params) - roi.value)
191      print "sum %0.1f" % np.sum(erf) # Temp.
192      return erf
193 <    ## End of constrained lma erf() function.
193 >    """End of constrained lma erf() function."""
194      
195    
196   def mclma(rois):
197      """
198 <    Constrained Levenberg-Marquart Algorithm fitting based on SciPy.
199 <    This function takes a single ba_class RegionOfInterest object (roi).
200 <    It modifies roi.params in place.  It returns the number of iterations completed.
198 >    *Multiple-curve Constrained Levenberg-Marquart Algorithm.*
199 >    Least-squares curve-fitting is performed on the data in the specified ROIs.
200 >    The ROIs must contain SPR data, a model, and initial (guess) parameters.
201 >    The roi.params values are modified in place.
202 >    The maximum number of iterations is set to 999.
203 >    
204 >    :param rois: ROI containing data, a model, and initial parameter guesses.
205 >    :type rois: List of RegOfInterest
206 >    :returns: number of iterations completed.
207      """
208      pval = []  ## List of values for floating parameters needed by leastsq().
209  
210      for roi in rois:
211          ## Check that pd contains all 4 keys: value, min, max, fixed.
212 <        checkparams(roi.params)
212 >        __checkparams(roi.params)
213          ## For each parameter with its dictionary...
214          for pkey,pd in roi.params.iteritems():
215              if (pd['fixed'] == 'fixed'):  
# Line 179 | Line 220
220                  pval.append(pd['value'])  ## e.g. 1e5
221                  ## Perform transform based on min/max constraints.
222                  a, b = float(pd['min']), float(pd['max'])
223 <                pval[-1] = itransform(a, b, float(pval[-1]))
223 >                pval[-1] = __itransform(a, b, float(pval[-1]))
224                  ## Print out intial guess and its transform. Temporary.
225                  print pkey, pd['value'], pval[-1]
226              else:  
# Line 187 | Line 228
228                  pass
229      
230      ## Fit.
231 <    p1, success = leastsq(mcerf, pval, args=(rois), maxfev=999)
231 >    p1, success = leastsq(__mcerf, pval, args=(rois), maxfev=999)
232      return success
233 <    ## End of multi-roi constrained lma() function.
233 >    """End of multi-roi constrained lma() function."""
234      
235      
236 < def mcerf(pfloat, rois):
236 > def __mcerf(pfloat, rois):
237      """
238      Error estimation function for use with constrained LMA.
239      Use the model assigned to this roi to simulation data vs time.  
# Line 213 | Line 254
254                  ## Floating parameter. Transform back.
255                  i += 1
256                  a, b = float(pd['min']), float(pd['max'])
257 <                pd['value'] = transform(a, b, pfloat[i])
257 >                pd['value'] = __transform(a, b, pfloat[i])
258                  print "%s %.4f" % (pkey, pd['value']), # Temp.
259                  
260      ## Second, do the floating parameters shared across ROIs.
# Line 234 | Line 275
275      print "ssq %0.1f" % sum(erf**2) # Temp.
276      
277      return erf
278 <    ## End of multi-roi constrained lma erf() function.
278 >    """End of multi-roi constrained lma error function mcerf()."""
279      
280      
281   ################################# End of module #################################

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines