ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/fit_module.py
Revision: 54
Committed: Wed Feb 16 21:29:02 2011 UTC (8 years, 8 months ago) by clausted
File size: 11472 byte(s)
Log Message:
Added fit_module and mdl_module to sphinx documentation.
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 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__ = "110216"
37
38
39 ## Import libraries
40 import ba_class as ba
41 from scipy.optimize import leastsq
42 from copy import deepcopy
43 from numpy import log, exp, tanh, arctanh
44 from numpy import sum, hstack, zeros
45
46
47 def lma(roi):
48 """
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)
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)
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."""
70
71
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).
76 This function will write the values to the parameter dictionary (roi.params).
77 This allows the LMA fitter to adjust the parameters in place.
78 """
79 for i,key in enumerate(roi.params):
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):
87 """
88 Check that each dictionary in the params dictionary contains the four
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
94 for key,dict in params.iteritems():
95 ## First check that keys are there.
96 if ('value' not in dict.keys()): dict['value'] = 0
97 if ('min' not in dict.keys()): dict['min'] = dict['value']
98 if ('max' not in dict.keys()): dict['max'] = dict['value']
99 if ('fixed' not in dict.keys()): dict['fixed'] = 'fixed'
100 ## Check that when parameters float, min<value<max.
101 if (dict['fixed'] == 'float'):
102 if (dict['min'] >= dict['max']):
103 ## Set max to be 1+min.
104 dict['max'] = float(dict['min']) + 1.0
105 flag = True
106 if (dict['value'] <= dict['min']) or (dict['value'] >= dict['max']):
107 ## Set initial value halfway between min and max.
108 dict['value'] = 0.5 * (dict['min'] + dict['max'])
109 flag = True
110 if (flag == True):
111 print "Parameter min/max errors: Modifications were made automatically."
112 return
113 """End of __checkparams() function."""
114
115
116 #### The code below is a variant allowing parameter constraints ####
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.
121 """
122 ## Considered transform: x^2 for (0,+inf)
123 ## Considered transform: a*tanh(x) for (-a,+a)
124 ## Considered transform: a+((b-a)/(1+exp(-x)) for (a,b) but it has problems.
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 ) * 0.5
129 return x1
130 """End of __transform() function."""
131
132
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.
137 """
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."""
142
143
144 def clma(roi):
145 """
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)
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)]
161 ## Adjust list based on min/max constraints.
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 #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)
171 return success
172 """End of constrained lma() function."""
173
174
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.
179 Provide a list of values for the floating parameters (pfloat).
180 This list is adjusted based on min/max constraints in the parameter dictionary.
181 This function will write the values to the parameter dictionary (roi.params).
182 This allows the LMA fitter to adjust the parameters in place.
183 """
184 params = roi.params
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])
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."""
194
195
196 def mclma(rois):
197 """
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)
213 ## For each parameter with its dictionary...
214 for pkey,pd in roi.params.iteritems():
215 if (pd['fixed'] == 'fixed'):
216 ## Fixed parameter.
217 pass
218 elif (pd['fixed'] == 'float'):
219 ## Floating parameter.
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]))
224 ## Print out intial guess and its transform. Temporary.
225 print pkey, pd['value'], pval[-1]
226 else:
227 ## Shared, fixed to a floating parameter in another ROI.
228 pass
229
230 ## Fit.
231 p1, success = leastsq(__mcerf, pval, args=(rois), maxfev=999)
232 return success
233 """End of multi-roi constrained lma() function."""
234
235
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.
240 Provide a list of values for the floating parameters (pfloat).
241 This list is adjusted based on min/max constraints in the parameter dictionary.
242 This function will write the values to the parameter dictionary (roi.params).
243 This allows the LMA fitter to adjust the parameters in place.
244 """
245 i = -1 ## Index for pfloat.
246
247 ## Insert the new guesses from pfloat into ROI parameters.
248 ## First, get all of the floating parameters done.
249 for j,roi in enumerate(rois):
250 print j, # Temp
251 ## For each parameter with its dictionary...
252 for pkey,pd in roi.params.iteritems():
253 if (pd['fixed'] == 'float'):
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])
258 print "%s %.4f" % (pkey, pd['value']), # Temp.
259
260 ## Second, do the floating parameters shared across ROIs.
261 for j,roi in enumerate(rois):
262 for pkey,pd in roi.params.iteritems():
263 if (pd['fixed'] == 'fixed') or (pd['fixed'] == 'float'):
264 pass
265 else:
266 ## Shared--fixed to a floating parameter in another ROI.
267 refroi = int(pd['fixed'])
268 pd['value'] = rois[refroi].params[pkey]['value']
269
270 ## Calculate model error function for each roi. Then concatenate them.
271 erf = zeros(0)
272 for roi in rois:
273 e = roi.model(roi.time, roi.value, roi.params) - roi.value
274 erf = hstack( (erf, e) ) ## Concatenation of nparrays.
275 print "ssq %0.1f" % sum(erf**2) # Temp.
276
277 return erf
278 """End of multi-roi constrained lma error function mcerf()."""
279
280
281 ################################# End of module #################################