Moved hooke.plugin.fit -> polymer_fit & updated to Plugin/Command architecture.
[hooke.git] / hooke / util / fit.py
index 4ae0e18da3a42847ed9baaa1284c64b1ea19aa1e..d5e4d0364493ce757add0f6daa3716a9c8edf13e 100644 (file)
@@ -21,6 +21,7 @@
 
 from numpy import arange, ndarray
 from scipy.optimize import leastsq
+import scipy.optimize
 
 
 class PoorFit (ValueError):
@@ -29,6 +30,10 @@ class PoorFit (ValueError):
 class ModelFitter (object):
     """A convenient wrapper around :func:`scipy.optimize.leastsq`.
 
+    TODO: look into :mod:`scipy.odr` as an alternative fitting
+    algorithm (minimizing the sum of squares of orthogonal distances,
+    vs. minimizing y distances).
+
     Parameters
     ----------
     d_data : array_like
@@ -217,3 +222,54 @@ class ModelFitter (object):
                     'convergence flag': ier,
                     })
         return params
+
+# Example ORD code from the old fit.py
+#        def dist(px,py,linex,liney):
+#            distancesx=scipy.array([(px-x)**2 for x in linex])
+#            minindex=numpy.argmin(distancesx)
+#            print px, linex[0], linex[-1]
+#            return (py-liney[minindex])**2
+#
+#
+#        def f_wlc(params,x,T=T):
+#            '''
+#            wlc function for ODR fitting
+#            '''
+#            lambd,pii=params
+#            Kb=(1.38065e-23)
+#            therm=Kb*T
+#            y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
+#            return y
+#
+#        def f_wlc_plfix(params,x,pl_value=pl_value,T=T):
+#            '''
+#            wlc function for ODR fitting
+#            '''
+#            lambd=params
+#            pii=1/pl_value
+#            Kb=(1.38065e-23)
+#            therm=Kb*T
+#            y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
+#            return y
+#
+#        #make the ODR fit
+#        realdata=scipy.odr.RealData(xchunk_corr_up,ychunk_corr_up)
+#        if pl_value:
+#            model=scipy.odr.Model(f_wlc_plfix)
+#            o = scipy.odr.ODR(realdata, model, p0_plfix)
+#        else:
+#            model=scipy.odr.Model(f_wlc)
+#            o = scipy.odr.ODR(realdata, model, p0)
+#
+#        o.set_job(fit_type=2)
+#        out=o.run()
+#        fit_out=[(1/i) for i in out.beta]
+#
+#        #Calculate fit errors from output standard deviations.
+#        #We must propagate the error because we fit the *inverse* parameters!
+#        #The error = (error of the inverse)*(value**2)
+#        fit_errors=[]
+#        for sd,value in zip(out.sd_beta, fit_out):
+#            err_real=sd*(value**2)
+#            fit_errors.append(err_real)
+