From e46fcec6987944cad15ae8d0c59ef3a526da423f Mon Sep 17 00:00:00 2001 From: Wirawan Purwanto Date: Thu, 21 May 2015 11:32:01 -0400 Subject: [PATCH] * Imported more recent function fitting facility from Cr2 project. --- math/fitting/__init__.py | 14 +- math/fitting/funcs_pec.py | 141 +++++++++++++++++ math/fitting/funcs_physics.py | 49 ++++++ math/fitting/funcs_simple.py | 276 ++++++++++++++++++++++++++++++++++ 4 files changed, 475 insertions(+), 5 deletions(-) create mode 100644 math/fitting/funcs_pec.py create mode 100644 math/fitting/funcs_physics.py create mode 100644 math/fitting/funcs_simple.py diff --git a/math/fitting/__init__.py b/math/fitting/__init__.py index 505d0f8..19b8b6e 100644 --- a/math/fitting/__init__.py +++ b/math/fitting/__init__.py @@ -104,7 +104,8 @@ def fit_func(Funct, Data=None, Guess=None, Params=None, N is the dimensionality of the domain, while M is the number of data points, whose count must be equal to the size of y data below. - For a 2-D fitting, for example, x should be a column array. + For a 2-D curve (y = f(x)) fitting, for example, + x should be a column array. An input guess for the parameters can be specified via Guess argument. It is an ordered list of scalar values for these parameters. @@ -120,8 +121,8 @@ def fit_func(Funct, Data=None, Guess=None, Params=None, If "dy" is specified, then "w" is defined to be (1.0 / dy**2), per usual convention. - Inspect Poly_base, Poly_order2, and other similar function classes in this - module to see the example of the Funct function. + Inspect Poly_base, Poly_order2, and other similar function classes in the + funcs_poly module to see the example of the Funct function. The measurement (input) datasets, against which the function is to be fitted, can be specified in one of two ways: @@ -209,7 +210,9 @@ def fit_func(Funct, Data=None, Guess=None, Params=None, # Try to provide an initial guess # This is an older version with y-only argument Guess = Funct.Guess(y) - elif Guess == None: # VERY OLD, DO NOT USE ANYMORE! + elif Guess == None: + # VERY OLD, DO NOT USE ANYMORE! Will likely not work for anythingnonlinear + # functions. Guess = [ y.mean() ] + [0.0, 0.0] * len(x) if use_lmfit: @@ -471,7 +474,8 @@ class fit_func_base(object): - TODO: dict-like Guess should be made possible. - otherwise, the guess values will be used as the initial values. - + Refer to various function objects in wpylib.math.fitting.funcs_simple + for actual examples of how to use and create your own fit_func_base object. """ class multi_fit_opts(dict): """A class for defining default control parameters for different fit methods. diff --git a/math/fitting/funcs_pec.py b/math/fitting/funcs_pec.py new file mode 100644 index 0000000..87546e3 --- /dev/null +++ b/math/fitting/funcs_pec.py @@ -0,0 +1,141 @@ +# +# wpylib.math.fitting.funcs_pec module +# Created: 20150521 +# Wirawan Purwanto +# +# Imported 20150521 from Cr2_analysis_cbs.py +# (dated 20141017, CVS rev 1.143). +# + +""" +wpylib.math.fitting.funcs_pec module +A library of simple f(x) functions for PEC fitting + +For use with OO-style x-y curve fitting interface. +""" + +import numpy + + +class harm_fit_func(fit_func_base): + """Harmonic function object. + For use with fit_func function on a PEC. + + Functional form: + + E0 + 0.5 * k * (x - re)**2 + + Coefficients: + * C[0] = energy minimum + * C[1] = spring constant + * C[2] = equilibrium distance + """ + dim = 1 # a function with 1-D domain + param_names = ('E0', 'k', 'r0') + def __call__(self, C, x): + xdisp = (x[0] - C[2]) + y = C[0] + 0.5 * C[1] * xdisp**2 + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + fit_rslt = fit_harm(x[0], y) + self.guess_params = tuple(fit_rslt[0]) + return self.guess_params + + +class harmcube_fit_func(fit_func_base): + """Harmonic + cubic term function object. + For use with fit_func function on a PEC. + + Functional form: + + E0 + 0.5 * k * (x - re)**2 + cub * (x - re)**3; + + Coefficients: + * C[0] = energy minimum + * C[1] = spring constant + * C[2] = equilibrium distance + * C[3] = nonlinear (cubic) constant + """ + dim = 1 # a function with 1-D domain + param_names = ('E0', 'k', 'r0', 'c3') + def __call__(self, C, x): + xdisp = (x[0] - C[2]) + y = C[0] + 0.5 * C[1] * xdisp**2 + C[3] * xdisp**3 + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + fit_rslt = fit_harm(x[0], y) + self.guess_params = tuple(fit_rslt[0]) + (0,) + return self.guess_params + def Guess_xy_old(self, x, y): + imin = numpy.argmin(y) + return (y[imin], 2, x[0][imin], 0.00001) + + +class morse2_fit_func(fit_func_base): + """Morse2 function object. + For use with fit_func function. + + Functional form: + + E0 + 0.5 * k / a**2 * (1 - exp(-a * (x - re)))**2 + + Coefficients: + * C[0] = energy minimum + * C[1] = spring constant + * C[2] = equilibrium distance + * C[3] = nonlinear constant + """ + dim = 1 # a function with 1-D domain + param_names = ('E0', 'k', 'r0', 'a') + def __call__(self, C, x): + from numpy import exp + E0, k, r0, a = self.get_params(C, *(self.param_names)) + y = E0 + 0.5 * k / a**2 * (1 - exp(-a * (x[0] - r0)))**2 + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + imin = numpy.argmin(y) + harm_params = fit_harm(x[0], y) + if self.debug >= 10: + print "Initial guess by fit_harm gives: ", harm_params + self.guess_params = (y[imin], harm_params[0][1], x[0][imin], 0.01 * harm_params[0][1]) + return self.guess_params + def Guess_xy_old(self, x, y): + imin = numpy.argmin(y) + return (y[imin], 2, x[0][imin], 0.01) + + +class ext3Bmorse2_fit_func(fit_func_base): + """ext3Bmorse2 function object. + For use with fit_func function. + + Functional form: + + E0 + 0.5 * k / a**2 * (1 - exp(-a * (x - re)))**2 + + C3 * (1 - exp(-a * (x - re)))**3 + + Coefficients: + * C[0] = energy minimum + * C[1] = spring constant + * C[2] = equilibrium distance + * C[3] = nonlinear constant + * C[4] = coefficient of cubic term + """ + dim = 1 # a function with 1-D domain + def __call__(self, C, x): + from numpy import exp + E = 1 - exp(-C[3] * (x[0] - C[2])) + y = C[0] + 0.5 * C[1] / C[3]**2 * E**2 + C[4] * E**3 + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + imin = numpy.argmin(y) + harm_params = fit_harm(x[0], y) + if self.debug >= 10: + print "Initial guess by fit_harm gives: ", harm_params + self.guess_params = (y[imin], harm_params[0][1], x[0][imin], 0.01 * harm_params[0][1], 0) + return self.guess_params + + diff --git a/math/fitting/funcs_physics.py b/math/fitting/funcs_physics.py new file mode 100644 index 0000000..5d5ade0 --- /dev/null +++ b/math/fitting/funcs_physics.py @@ -0,0 +1,49 @@ +# +# wpylib.math.fitting.funcs_physics module +# Created: 20150521 +# Wirawan Purwanto +# +# Imported 20150521 from Cr2_analysis_cbs.py +# (dated 20141017, CVS rev 1.143). +# + +""" +wpylib.math.fitting.funcs_physics module +A library of simple f(x) functions for physics-related common functional fitting + +For use with OO-style x-y curve fitting interface. +""" + +import numpy + + +class FermiDirac_fit_func(fit_func_base): + """Fermi-Dirac function object. + For use with fit_func function. + + Functional form: + + C[0] * (exp((x - C[1]) / C[2]) + 1)^-1 + + Coefficients: + * C[0] = amplitude + * C[1] = transition "temperature" + * C[2] = "smearing temperature" + """ + dim = 1 # a function with 1-D domain + param_names = ('A', 'F', 'T') + # FIXME: Not good yet!!! + F_guess = 1.9 + T_guess = 0.05 + def __call__(self, C, x): + from numpy import exp + A, F, T = self.get_params(C, *(self.param_names)) + y = A * (exp((x[0] - F) / T) + 1)**(-1) + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + imin = numpy.argmin(y) + self.guess_params = (y[imin], self.F_guess, self.T_guess) + return self.guess_params + + diff --git a/math/fitting/funcs_simple.py b/math/fitting/funcs_simple.py new file mode 100644 index 0000000..d43d2e0 --- /dev/null +++ b/math/fitting/funcs_simple.py @@ -0,0 +1,276 @@ +# +# wpylib.math.fitting.funcs_simple module +# Created: 20150520 +# Wirawan Purwanto +# +# Imported 20150520 from Cr2_analysis_cbs.py +# (dated 20141017, CVS rev 1.143). +# + +""" +wpylib.math.fitting.funcs_simple module +A library of simple f(x) functions for fitting + +For use with OO-style x-y curve fitting interface. +""" + +import numpy + + +# Some simple function fitting--to aid fitting the complex ones later + +def fit_linear(x, y): + """Warning: the ansatz for fitting is + C[0] + C[1]*x + so I have to reverse the order of fit parameters. + """ + rslt = numpy.polyfit(x, y, 1, full=True) + return (rslt[0][::-1],) + rslt + + +def fit_harm(x, y): + """Do a quadratic fit using poly fit and return it in terms of coeffs + like this one: + + C0 + 0.5 * C1 * (x - C2)**2 + + => 0.5*C1*x**2 - C1*C2*x + (C0 + 0.5 * C1 * C2**2) + + Polyfit gives: + a * x**2 + b * x + c + + Equating the two, we get: + + C1 = 2 * a + C2 = -b/C1 + C0 = c - 0.5*C1*C2**2 + + This function returns the recast parameters plus the original + fit output. + """ + rslt = numpy.polyfit(x, y, 2, full=True) + + (a,b,c) = rslt[0] + C1 = 2*a + C2 = -b/C1 + C0 = c - 0.5*C1*C2**2 + + return ((C0,C1,C2),) + rslt + + + +# fit_func-style functional ansatz + +class const_fit_func(fit_func_base): + """Constant function object. + For use with fit_func function on a PEC. + + Functional form: + + C[0] + + Coefficients: + * C[0] = the constant sought + """ + dim = 1 # a function with 1-D domain + param_names = ('c') + def __call__(self, C, x): + from numpy import exp + y = C[0] + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + self.guess_params = (numpy.average(y),) + return self.guess_params + + +class linear_fit_func(fit_func_base): + """Linear function object. + For use with fit_func function. + + Functional form: + + a + b * x + + Coefficients: + * C[0] = a + * C[1] = b + """ + dim = 1 # a function with 1-D domain + param_names = ('a', 'b') + def __call__(self, C, x): + y = C[0] + C[1] * x[0] + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + fit_rslt = fit_linear(x[0], y) + self.guess_params = tuple(fit_rslt[0]) + return self.guess_params + + +class linear_leastsq_fit_func(linear_fit_func): + def fit(self, x, y, dy=None, fit_opts=None, Funct_hook=None, Guess=None): + from wpylib.math.fitting.linear import linregr2d_SZ + # Changed from: + # rslt = fit_linear_weighted(x,y,dy) + # to: + rslt = (x, y, sigma=None) + + self.last_fit = rslt[1] + # Retrofit for API compatibility: not necessarily meaningful + self.guess_params = rslt[0] + return rslt[0] + + +class exp_fit_func(fit_func_base): + """Exponential function object. + For use with fit_func function. + + Functional form: + + C[0] * (exp(C[1] * (x - C[2])) + + Coefficients: + * C[0] = amplitude + * C[1] = damping factor + * C[2] = offset + """ + dim = 1 # a function with 1-D domain + param_names = ['A', 'B', 'x0'] + A_guess = -2.62681 + B_guess = -9.05046 + x0_guess = 1.57327 + def __call__(self, C, x): + from numpy import exp + A, B, x0 = self.get_params(C, *(self.param_names)) + y = A * exp(B * (x[0] - x0)) + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + from numpy import abs + #y_abs = abs(y) + # can do linear fit to guess the params, + # but how to separate A and B*x0, I don't know. + #imin = numpy.argmin(y) + self.guess_params = (self.A_guess, self.B_guess, self.x0_guess) + return self.guess_params + + +class expm_fit_func(exp_fit_func): + """Similar to exp_fit_func but the exponent is always negative. + """ + def __call__(self, C, x): + from numpy import exp,abs + A, B, x0 = self.get_params(C, *(self.param_names)) + y = A * exp(-abs(B) * (x[0] - x0)) + self.func_call_hook(C, x, y) + return y + + +class powx_fit_func(fit_func_base): + """Power of x function object. + For use with fit_func function. + + Functional form: + + C[0] * ((x - C[2])**C[1]) + + Coefficients: + * C[0] = amplitude + * C[1] = exponent (< 0) + * C[2] = offset + """ + dim = 1 # a function with 1-D domain + param_names = ['A', 'B', 'x0'] + A_guess = -2.62681 + B_guess = -9.05046 + x0_guess = 1.57327 + def __call__(self, C, x): + from numpy import exp + A, B, x0 = self.get_params(C, *(self.param_names)) + y = A * (x[0] - x0)**B + self.func_call_hook(C, x, y) + return y + def Guess_xy(self, x, y): + from numpy import abs + #y_abs = abs(y) + # can do linear fit to guess the params, + # but how to separate A and B*x0, I don't know. + #imin = numpy.argmin(y) + self.guess_params = (self.A_guess, self.B_guess, self.x0_guess) + return self.guess_params + + +class invx_fit_func(powx_fit_func): + """Inverse of x function object that leads to 0 as x->infinity. + For use with fit_func function. + + Functional form: + + C[0] * ((x - C[2])**C[1]) + + Specialized for CBX1 extrapolation + Coefficients: + * C[0] = amplitude (< 0) + * C[1] = exponent (< 0) + * C[2] = offset (> 0) + """ + """ + /home/wirawan/Work/GAFQMC/expt/qmc/Cr2/CBS-TZ-QZ/UHF-CBS/20140128/Exp-CBX1.d/fit-invx.plt + + Iteration 154 + WSSR : 0.875715 delta(WSSR)/WSSR : -9.96404e-06 + delta(WSSR) : -8.72566e-06 limit for stopping : 1e-05 + lambda : 0.00174063 + + resultant parameter values + + A = -29.7924 + B = -13.2967 + x0 = 0.399396 + + After 154 iterations the fit converged. + final sum of squares of residuals : 0.875715 + rel. change during last iteration : -9.96404e-06 + + degrees of freedom (FIT_NDF) : 2 + rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 0.661708 + variance of residuals (reduced chisquare) = WSSR/ndf : 0.437858 + + Final set of parameters Asymptotic Standard Error + ======================= ========================== + + A = -29.7924 +/- 8027 (2.694e+04%) + B = -13.2967 +/- 196.1 (1474%) + x0 = 0.399396 +/- 21.4 (5357%) + + + correlation matrix of the fit parameters: + + A B x0 + A 1.000 + B 1.000 1.000 + x0 1.000 1.000 1.000 + + For some reason the fit code in python gives: + A,B,x0 = (-7028.1498486021028, -16.916447508009664, 2.2572321406455487e-06) + but they fit almost exactly the same in the region 1.8 <= r <= 3.0. + + """ + A_guess = -29.7924 + B_guess = -13.2967 + x0_guess = 0.399396 + def __init__(self): + from lmfit import Parameters + self.fit_method = "lmfit:leastsq" + p = Parameters() + p.add_many( + # (Name, Value, Vary, Min, Max, Expr) + ('A', -2.6, True, -1e6, -1e-9, None), + ('B', -2.0, True, None, -1e-9, None), + ('x0', 1.9, True, 1e-6, None, None), + # The values are just a placeholder. They will be set later. + ) + self.Params = p + +