ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/osprai/osprai/trunk/fit_module.py
(Generate patch)
# Line 2 | Line 2
2   fit: Curve-fitting module for SPRI data.
3   Christopher Lausted, Institute for Systems Biology,
4   OSPRAI developers
5 < Last modified on 100425 (yymmdd)
5 > Last modified on 100518 (yymmdd)
6  
7   Example:
8   #import fit_module as fit
# Line 12 | Line 12
12   #success = lma(ba1.roi[0])
13   #for i in ba1.roi: lma(i)
14   """
15 < __version__ = "100425"
15 > __version__ = "100518"
16  
17  
18   ## Import libraries
19   import ba_class as ba
20 import numpy as np
20   from scipy.optimize import leastsq
21   from copy import deepcopy
22   from numpy import log, exp, tanh, arctanh
23 < #from math import exp
23 > from numpy import sum, hstack, zeros
24  
25  
26   def lma(roi):
# Line 152 | Line 151
151          if (params[key]['fixed'] != True):
152              a, b = params[key]['min'], params[key]['max']
153              params[key]['value'] = transform(a, b, pfloat[i])
154 <            print "%s %f" % (key, params[key]['value']), # Temp.
154 >            print "%s %.2e" % (key, params[key]['value']), # Temp.
155      erf = (roi.model(roi.time, roi.value, params) - roi.value)
156      print "sum %0.1f" % np.sum(erf) # Temp.
157      return erf
158      ## End of constrained lma erf() function.
159      
160 <    
161 < '''    
163 < #### The code below is an alternative strategy that is not very robust. ####
164 < def clma_old(roi):
160 >  
161 > def mclma(rois):
162      """
163 <    Constrained Levenberg-Marquart fitting based on SciPy.
163 >    Constrained Levenberg-Marquart Algorithm fitting based on SciPy.
164      This function takes a single ba_class RegionOfInterest object (roi).
165      It modifies roi.params in place.  It returns the number of iterations completed.
169    The error estimate is increased when a parameter moves outside the given bounds.
166      """
167 <    params = roi.params
168 <    checkparams(params)
169 <    ## Create list initial values for the floating parameters.
170 <    p0 = [dict['value'] for key,dict in params.iteritems() if (dict['fixed'] != True)]
167 >    pval = []  ## List of values for floating parameters needed by leastsq().
168 >
169 >    for roi in rois:
170 >        ## Check that pd contains all 4 keys: value, min, max, fixed.
171 >        checkparams(roi.params)
172 >        ## For each parameter with its dictionary...
173 >        for pkey,pd in roi.params.iteritems():
174 >            if (pd['fixed'] == True):  
175 >                ## Fixed parameter.
176 >                pass
177 >            elif (pd['fixed'] == False):  
178 >                ## Floating parameter.
179 >                pval.append(pd['value'])  ## e.g. 1e5
180 >                ## Perform transform based on min/max constraints.
181 >                a, b = float(pd['min']), float(pd['max'])
182 >                pval[-1] = itransform(a, b, float(pval[-1]))
183 >                ## Print out intial guess and its transform. Temporary.
184 >                print pkey, pd['value'], pval[-1]
185 >            else:  
186 >                ## Shared, fixed to a floating parameter in another ROI.
187 >                pass
188 >    
189      ## Fit.
190 <    p1, success = leastsq(cerf_old, p0, args=(roi), maxfev=10000)
190 >    p1, success = leastsq(mcerf, pval, args=(rois), maxfev=999)
191      return success
192 <    ## End of clma() function.
192 >    ## End of multi-roi constrained lma() function.
193      
194      
195 < def cerf_old(pfloat, roi):
195 > def mcerf(pfloat, rois):
196      """
197 <    Constraining Error Function.
198 <    Use the model assigned to this roi to simulation data vs time and calculate the error.  
185 <    This allows the Constrained LMA fitter to adjust the parameters in place.
197 >    Error estimation function for use with constrained LMA.
198 >    Use the model assigned to this roi to simulation data vs time.  
199      Provide a list of values for the floating parameters (pfloat).
200 <    The penalty factor is 1*(distance outside bounds).
201 <    Versus regular LMA, it's less robust--initialize with reasonable parameter guesses.
200 >    This list is adjusted based on min/max constraints in the parameter dictionary.
201 >    This function will write the values to the parameter dictionary (roi.params).
202 >    This allows the LMA fitter to adjust the parameters in place.
203      """
204 <    penalty = 1
205 <    for i,key in enumerate(roi.params):
206 <        if (roi.params[key]['fixed'] != True):
207 <            ## Substitute floating parameters in roi.params.
208 <            roi.params[key]['value'] = pfloat[i]
209 <            ## Impose penalties when outside constraints.
210 <            a = roi.params[key]['min']
211 <            b = roi.params[key]['max']
212 <            x = roi.params[key]['value']
213 <            if (x<a and a==0): penalty = penalty * (1 + (a-x))
214 <            if (x<a and a!=0): penalty = penalty * (1 + (a-x)/a)
215 <            if (x>b and b==0): penalty = penalty * (1 + (x-b))
216 <            if (x>b and b!=0): penalty = penalty * (1 + (x-b)/b)
217 <            print "%s %f" % (key, x),
218 <    print "penalty %0.2f" % (penalty)
219 <    errval = roi.model(roi.time, roi.value, roi.params) - roi.value
220 <    errval = errval * penalty
221 <    return errval
222 <    ## End of cerf() function.
223 < '''
224 <
204 >    i = -1  ## Index for pfloat.
205 >    
206 >    ## Insert the new guesses from pfloat into ROI parameters.
207 >    ## First, get all of the floating parameters done.
208 >    for j,roi in enumerate(rois):
209 >        print j, # Temp
210 >        ## For each parameter with its dictionary...
211 >        for pkey,pd in roi.params.iteritems():
212 >            if (pd['fixed'] == False): ## Floating parameter. Transform back.
213 >                i += 1
214 >                a, b = float(pd['min']), float(pd['max'])
215 >                pd['value'] = transform(a, b, pfloat[i])
216 >                print "%s %.4f" % (pkey, pd['value']), # Temp.
217 >                
218 >    ## Second, do the floating parameters shared across ROIs.
219 >    for j,roi in enumerate(rois):
220 >        for pkey,pd in roi.params.iteritems():
221 >            if (pd['fixed'] == True) or (pd['fixed'] == False):
222 >                pass
223 >            else:  ## Shared, fixed to a floating parameter in another ROI.
224 >                refroi = int(pd['fixed'])
225 >                pd['value'] = rois[refroi].params[pkey]['value']
226 >    
227 >    ## Calculate model error function for each roi. Then concatenate them.
228 >    erf = zeros(0)
229 >    for roi in rois:
230 >        e = roi.model(roi.time, roi.value, roi.params) - roi.value
231 >        erf = hstack( (erf, e) )
232 >    print "sum %0.1f" % sum(erf) # Temp.
233 >    
234 >    return erf
235 >    ## End of multi-roi constrained lma erf() function.
236 >    
237 >    
238   ################################# End of module #################################

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines