--- osprai/trunk/fit_module.py 2010/04/23 01:45:38 19 +++ osprai/trunk/fit_module.py 2010/04/26 21:59:45 20 @@ -2,7 +2,7 @@ fit: Curve-fitting module for SPRI data. Christopher Lausted, Institute for Systems Biology, OSPRAI developers -Last modified on 100422 (yymmdd) +Last modified on 100425 (yymmdd) Example: #import fit_module as fit @@ -12,7 +12,7 @@ #success = lma(ba1.roi[0]) #for i in ba1.roi: lma(i) """ -__version__ = "100422" +__version__ = "100425" ## Import libraries @@ -20,6 +20,8 @@ import numpy as np from scipy.optimize import leastsq from copy import deepcopy +from numpy import log, exp, tanh, arctanh +#from math import exp def lma(roi): @@ -28,22 +30,20 @@ This function takes a single ba_class RegionOfInterest object (roi). It modifies roi.params in place. It returns the number of iterations completed. """ - time = roi.time - data = roi.value params = roi.params checkparams(params) ## Create list initial values for the floating parameters. p0 = [dict['value'] for key,dict in params.iteritems() if (dict['fixed'] != True)] ## Error function. - erf = lambda p: (model_interface(roi, time, data, p) - data) + erf = lambda p: (model_interface(roi, p) - roi.value) ## Fit. - #p1, success = leastsq(erf, p0, args=(), maxfev=10000) - p1, success = leastsq(erf, p0, maxfev=10000) + #p1, success = leastsq(erf, p0, args=(), maxfev=999) + p1, success = leastsq(erf, p0, maxfev=500) return success ## End of constlma() function. -def model_interface(roi, time, data, pfloat): +def model_interface(roi, pfloat): """ Use the model assigned to this roi to simulation data vs time. Provide a list of values for the floating parameters (pfloat). @@ -53,21 +53,160 @@ for i,key in enumerate(roi.params): if (roi.params[key]['fixed'] != True): roi.params[key]['value'] = pfloat[i] - return roi.model(time, data, roi.params) + return roi.model(roi.time, roi.value, roi.params) def checkparams(params): """ Check that each dictionary in the params dictionary contains the four - keys (value, min, max, fixed) and add them if necessary. Default is 'fixed':True. + keys (value, min, max, fixed) and add them if necessary. Default is 'fixed':True. + Also check that when parameters float, min= dict['max']): + ## Set max to be 1+min. + dict['max'] = float(dict['min']) + 1.0 + flag = True + if (dict['value'] <= dict['min']) or (dict['value'] >= dict['max']): + ## Set initial value halfway between min and max. + dict['value'] = 0.5 * (dict['min'] + dict['max']) + flag = True + if (flag == True): + print "Parameter min/max errors: Modifications were made automatically." return ## End of checkparams() function. +#### The code below is an alternative, non-working strategy. #### +def transform(a, b, x0): + """ + A transformation to convert x from range (-inf,+inf) to (a,b). + This is used to help constrain the outputs of the Levenberg-Marquart algorithm. + """ + ## Considered transform: x^2 for (0,+inf) + ## Considered transform: a*tanh(x) for (-a,+a) + ## Considered transform: a+((b-a)/(1+exp(-x)) for (a,b) but it has problems. + #x = max(-x0, -709) ## Function math.exp(>709) overflows. + #x1 = a + ( (b-a) / (1+exp(-x) )) + ## Try transform: y = ((b-a)*tanh(x)+b+a)/2 + x1 = ( (b-a)*tanh(x0) + b + a ) / 2 + return x1 + ## End of transform() function. + + +def itransform(a, b, x1): + """ + The (inverse) transformation to convert x from range(a,b) to (-inf,+inf). + This is used to help constrain the outputs of the Levenberg-Marquart algorithm. + """ + #x0 = -log( ((b-a)/(x1-a)) - 1 ) + x0 = arctanh( (2*x1-b-a) / (b-a) ) + return x0 + ## End of itransform() function. + + +## Warning: The next two functions are not yet working. +def clma(roi): + """ + Constrained Levenberg-Marquart Algorithm fitting based on SciPy. + This function takes a single ba_class RegionOfInterest object (roi). + It modifies roi.params in place. It returns the number of iterations completed. + """ + params = roi.params + checkparams(params) + + ## Create list of initial values for the floating parameters. + p0 = [dict['value'] for key,dict in params.iteritems() if (dict['fixed'] != True)] + ## Adjust list based on min/max constraints. + for i,key in enumerate(params): + if (params[key]['fixed'] != True): + a, b = params[key]['min'], params[key]['max'] + p0[i] = itransform(a, b, float(p0[i])) + # Temp + print key, params[key]['value'], p0[i] + + ## Fit. + p1, success = leastsq(cerf, p0, args=(roi), maxfev=999) + return success + ## End of constrained lma() function. + + +def cerf(pfloat, roi): + """ + Error estimation function for use with constrained LMA. + Use the model assigned to this roi to simulation data vs time. + Provide a list of values for the floating parameters (pfloat). + This list is adjusted based on min/max constraints in the parameter dictionary. + This function will write the values to the parameter dictionary (roi.params). + This allows the LMA fitter to adjust the parameters in place. + """ + params = roi.params + for i,key in enumerate(params): + if (params[key]['fixed'] != True): + a, b = params[key]['min'], params[key]['max'] + params[key]['value'] = transform(a, b, pfloat[i]) + print "%s %f" % (key, params[key]['value']), # Temp. + erf = (roi.model(roi.time, roi.value, params) - roi.value) + print "sum %0.1f" % np.sum(erf) # Temp. + return erf + ## End of constrained lma erf() function. + + +''' +#### The code below is an alternative strategy that is not very robust. #### +def clma_old(roi): + """ + Constrained Levenberg-Marquart fitting based on SciPy. + This function takes a single ba_class RegionOfInterest object (roi). + It modifies roi.params in place. It returns the number of iterations completed. + The error estimate is increased when a parameter moves outside the given bounds. + """ + params = roi.params + checkparams(params) + ## Create list initial values for the floating parameters. + p0 = [dict['value'] for key,dict in params.iteritems() if (dict['fixed'] != True)] + ## Fit. + p1, success = leastsq(cerf_old, p0, args=(roi), maxfev=10000) + return success + ## End of clma() function. + + +def cerf_old(pfloat, roi): + """ + Constraining Error Function. + Use the model assigned to this roi to simulation data vs time and calculate the error. + This allows the Constrained LMA fitter to adjust the parameters in place. + Provide a list of values for the floating parameters (pfloat). + The penalty factor is 1*(distance outside bounds). + Versus regular LMA, it's less robust--initialize with reasonable parameter guesses. + """ + penalty = 1 + for i,key in enumerate(roi.params): + if (roi.params[key]['fixed'] != True): + ## Substitute floating parameters in roi.params. + roi.params[key]['value'] = pfloat[i] + ## Impose penalties when outside constraints. + a = roi.params[key]['min'] + b = roi.params[key]['max'] + x = roi.params[key]['value'] + if (xb and b==0): penalty = penalty * (1 + (x-b)) + if (x>b and b!=0): penalty = penalty * (1 + (x-b)/b) + print "%s %f" % (key, x), + print "penalty %0.2f" % (penalty) + errval = roi.model(roi.time, roi.value, roi.params) - roi.value + errval = errval * penalty + return errval + ## End of cerf() function. +''' + ################################# End of module ################################# \ No newline at end of file