ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/fit_module.py
Revision: 65
Committed: Fri May 20 16:33:02 2011 UTC (8 years, 5 months ago) by clausted
File size: 12653 byte(s)
Log Message:
Added much testing to all modules using docstrings and doctest.  To help this, I also added two optional parameters to mclma(), verbose and maxiter.  For example, use "mclma(rois, verbose=False, maxiter=99)" to run the regression up to 99 iterations without echoing progress to the standard output.

Line File contents
1 """
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 osprai_one as osp, numpy as np
27 >>> ba1 = osp.BiosensorArray(1,300)
28 >>> roi = ba1.roi[0]
29 >>> ## Simple 1to1 binding curve.
30 >>> ## Note KD = 2.2e-2 / 4.4e-4 = 5e-7, log10(5e-7) = -6.3
31 >>> roi.simulate(bulkRI=0.0, kon=4.4e4, koff=2.2e-2, conc=7.7e-7)
32 >>> ## Create a copy ba2 and use for fitting three parameters.
33 >>> ba2 = osp.deepcopy(ba1)
34 >>> roi = ba2.roi[0]
35 >>> roi.model = osp.simple1to1
36 >>> roi.params = osp.simple1to1_def_params()
37 >>> roi.params['rmax']['fixed'] = 'float'
38 >>> roi.params['kon']['fixed'] = 'float'
39 >>> roi.params['koff']['fixed'] = 'float'
40 >>> roi.params['conc']['value'] = 7.7e-7
41 >>> ## Do curvefitting, then check result.
42 >>> osp.mclma(ba2.roi, verbose=False)
43 2
44 >>> roi.value = osp.zeros(len(roi.value))
45 >>> roi.value = roi.model(roi.time, roi.value, roi.params)
46 >>> kD = float(roi.params['koff']['value']) / float(roi.params['kon']['value'])
47 >>> print round((np.log10(kD)), 1)
48 -6.3
49
50 """
51 __version__ = "110519"
52
53
54 ## Import libraries
55 import ba_class as ba
56 from scipy.optimize import leastsq
57 from copy import deepcopy
58 from numpy import log, exp, tanh, arctanh
59 from numpy import sum, hstack, zeros
60
61
62 def lma(roi):
63 """
64 *Normal LMA curve-fitting. (Deprecated)*
65 Least-squares curve-fitting is performed on the data in the specified ROI.
66 The ROI must contain SPR data, a model, and initial (guess) parameters.
67 The roi.params values are modified in place.
68 The maximum number of iterations is set to 500.
69
70 :param roi: ROI containing data, a model, and initial parameter guesses.
71 :type roi: RegOfInterest
72 :returns: number of iterations completed.
73 """
74 params = roi.params
75 __checkparams(params)
76 ## Create list initial values for the floating parameters.
77 p0 = [dict['value'] for key,dict in params.iteritems() if (dict['fixed'] != True)]
78 ## Error function.
79 erf = lambda p: (__model_interface(roi, p) - roi.value)
80 ## Fit.
81 #p1, success = leastsq(erf, p0, args=(), maxfev=999)
82 p1, success = leastsq(erf, p0, maxfev=500)
83 return success
84 """End of constlma() function."""
85
86
87 def __model_interface(roi, pfloat):
88 """
89 Use the model assigned to this roi to simulation data vs time.
90 Provide a list of values for the floating parameters (pfloat).
91 This function will write the values to the parameter dictionary (roi.params).
92 This allows the LMA fitter to adjust the parameters in place.
93 """
94 for i,key in enumerate(roi.params):
95 if (roi.params[key]['fixed'] != True):
96 roi.params[key]['value'] = pfloat[i]
97 return roi.model(roi.time, roi.value, roi.params)
98 """End of __model_interface() function."""
99
100
101 def __checkparams(params):
102 """
103 Check that each dictionary in the params dictionary contains the four
104 keys (value, min, max, fixed) and add them if necessary.
105 Default is 'fixed':'fixed'.
106 Also check that when parameters float, min<value<max.
107 """
108 flag = False
109 for key,dict in params.iteritems():
110 ## First check that keys are there.
111 if ('value' not in dict.keys()): dict['value'] = 0
112 if ('min' not in dict.keys()): dict['min'] = dict['value']
113 if ('max' not in dict.keys()): dict['max'] = dict['value']
114 if ('fixed' not in dict.keys()): dict['fixed'] = 'fixed'
115 ## Check that when parameters float, min<value<max.
116 if (dict['fixed'] == 'float'):
117 if (dict['min'] >= dict['max']):
118 ## Set max to be 1+min.
119 dict['max'] = float(dict['min']) + 1.0
120 flag = True
121 if (dict['value'] <= dict['min']) or (dict['value'] >= dict['max']):
122 ## Set initial value halfway between min and max.
123 dict['value'] = 0.5 * (dict['min'] + dict['max'])
124 flag = True
125 if (flag == True):
126 print "Parameter min/max errors: Modifications were made automatically."
127 return
128 """End of __checkparams() function."""
129
130
131 #### The code below is a variant allowing parameter constraints ####
132 def __transform(a, b, x0):
133 """
134 A transformation to convert x from range (-inf,+inf) to (a,b).
135 This is used to help constrain the outputs of the Levenberg-Marquart algorithm.
136 """
137 ## Considered transform: x^2 for (0,+inf)
138 ## Considered transform: a*tanh(x) for (-a,+a)
139 ## Considered transform: a+((b-a)/(1+exp(-x)) for (a,b) but it has problems.
140 #x = max(-x0, -709) ## Function math.exp(>709) overflows.
141 #x1 = a + ( (b-a) / (1+exp(-x) ))
142 ## Try transform: y = ((b-a)*tanh(x)+b+a)/2
143 x1 = ( (b-a)*tanh(x0) + b + a ) * 0.5
144 return x1
145 """End of __transform() function."""
146
147
148 def __itransform(a, b, x1):
149 """
150 The (inverse) transformation to convert x from range(a,b) to (-inf,+inf).
151 This is used to help constrain the outputs of the Levenberg-Marquart algorithm.
152 """
153 #x0 = -log( ((b-a)/(x1-a)) - 1 )
154 x0 = arctanh( (2*x1-b-a) / (b-a) )
155 return x0
156 """End of __itransform() function."""
157
158
159 def clma(roi):
160 """
161 *Constrained LMA curve-fitting. (Deprecated)*
162 Least-squares curve-fitting is performed on the data in the specified ROI.
163 The ROI must contain SPR data, a model, and initial (guess) parameters.
164 The roi.params values are modified in place.
165 The maximum number of iterations is set to 999.
166
167 :param roi: ROI containing data, a model, and initial parameter guesses.
168 :type roi: RegOfInterest
169 :returns: number of iterations completed.
170 """
171 params = roi.params
172 __checkparams(params)
173
174 ## Create list of initial values for the floating parameters.
175 p0 = [dict['value'] for key,dict in params.iteritems() if (dict['fixed'] != True)]
176 ## Adjust list based on min/max constraints.
177 for i,key in enumerate(params):
178 if (params[key]['fixed'] != True):
179 a, b = params[key]['min'], params[key]['max']
180 p0[i] = __itransform(a, b, float(p0[i]))
181
182 ## Fit.
183 p1, success = leastsq(__cerf, p0, args=(roi), maxfev=999)
184 return success
185 """End of constrained lma() function."""
186
187
188 def __cerf(pfloat, roi):
189 """
190 Error estimation function for use with constrained LMA.
191 Use the model assigned to this roi to simulation data vs time.
192 Provide a list of values for the floating parameters (pfloat).
193 This list is adjusted based on min/max constraints in the parameter dictionary.
194 This function will write the values to the parameter dictionary (roi.params).
195 This allows the LMA fitter to adjust the parameters in place.
196 """
197 params = roi.params
198 for i,key in enumerate(params):
199 if (params[key]['fixed'] != True):
200 a, b = params[key]['min'], params[key]['max']
201 params[key]['value'] = __transform(a, b, pfloat[i])
202 print "%s %.2e" % (key, params[key]['value']), # Temp.
203 erf = (roi.model(roi.time, roi.value, params) - roi.value)
204 print "sum %0.1f" % np.sum(erf) # Temp.
205 return erf
206 """End of constrained lma erf() function."""
207
208
209 def mclma(rois, maxiter=999, verbose=True):
210 """
211 *Multiple-curve Constrained Levenberg-Marquart Algorithm.*
212 Least-squares curve-fitting is performed on the data in the specified ROIs.
213 The ROIs must contain SPR data, a model, and initial (guess) parameters.
214 The roi.params values are modified in place.
215 The ``maxiter`` specifies the maximum number of iterations (default 999).
216 The ``verbose`` option is used to display progress (default True).
217
218 :param rois: ROI containing data, a model, and initial parameter guesses.
219 :type rois: List of RegOfInterest
220 :returns: number of iterations completed.
221 """
222 pval = [] ## List of values for floating parameters needed by leastsq().
223
224 for roi in rois:
225 ## Check that pd contains all 4 keys: value, min, max, fixed.
226 __checkparams(roi.params)
227 ## For each parameter with its dictionary...
228 for pkey,pd in roi.params.iteritems():
229 if (pd['fixed'] == 'fixed'):
230 ## Fixed parameter.
231 pass
232 elif (pd['fixed'] == 'float'):
233 ## Floating parameter.
234 pval.append(pd['value']) ## e.g. 1e5
235 ## Perform transform based on min/max constraints.
236 a, b = float(pd['min']), float(pd['max'])
237 pval[-1] = __itransform(a, b, float(pval[-1]))
238 ## Print intial guess and its transform, if wanted.
239 if verbose: print pkey, pd['value'], pval[-1]
240 else:
241 ## Shared, fixed to a floating parameter in another ROI.
242 pass
243
244 ## Fit.
245 p1, success = leastsq(__mcerf, pval, args=(rois,verbose), maxfev=maxiter)
246 return success
247 """End of multi-roi constrained lma() function."""
248
249
250 def __mcerf(pfloat, rois, verbose=True):
251 """
252 Error estimation function for use with constrained LMA.
253 Use the model assigned to this roi to simulation data vs time.
254 Provide a list of values for the floating parameters (pfloat).
255 This list is adjusted based on min/max constraints in the parameter dictionary.
256 This function will write the values to the parameter dictionary (roi.params).
257 This allows the LMA fitter to adjust the parameters in place.
258 """
259 i = -1 ## Index for pfloat.
260
261 ## Insert the new guesses from pfloat into ROI parameters.
262 ## First, get all of the floating parameters done.
263 for j,roi in enumerate(rois):
264 if verbose: print j,
265 ## For each parameter with its dictionary...
266 for pkey,pd in roi.params.iteritems():
267 if (pd['fixed'] == 'float'):
268 ## Floating parameter. Transform back.
269 i += 1
270 a, b = float(pd['min']), float(pd['max'])
271 pd['value'] = __transform(a, b, pfloat[i])
272 ## Print current state of floating param, if wanted.
273 if verbose: print "%s %.4f" % (pkey, pd['value']),
274
275 ## Second, do the floating parameters shared across ROIs.
276 for j,roi in enumerate(rois):
277 for pkey,pd in roi.params.iteritems():
278 if (pd['fixed'] == 'fixed') or (pd['fixed'] == 'float'):
279 pass
280 else:
281 ## Shared--fixed to a floating parameter in another ROI.
282 refroi = int(pd['fixed'])
283 pd['value'] = rois[refroi].params[pkey]['value']
284
285 ## Calculate model error function for each roi. Then concatenate them.
286 erf = zeros(0)
287 for roi in rois:
288 e = roi.model(roi.time, roi.value, roi.params) - roi.value
289 erf = hstack( (erf, e) ) ## Concatenation of nparrays.
290 ## Print sum of squares value, if wanted.
291 if verbose: print "ssq %0.1f" % sum(erf**2)
292
293 return erf
294 """End of multi-roi constrained lma error function mcerf()."""
295
296
297
298 def _test():
299 """
300 Automatic Code testing with doctest.
301 Enable use of ellipses as wildcard for returned text.
302 Doctest runs the example code in the docstrings.
303 """
304 import doctest
305 doctest.testmod(optionflags=doctest.ELLIPSIS)
306 return
307
308 if __name__ == "__main__":
309 """Use 'python fit_module.py' for code testing."""
310 _test()
311
312
313 ################################# End of module #################################