Tinkering

src/calculate/optimize/fit_coeffs.py

@@ -9,6 +9,7 @@
 from copy import deepcopy
 from scipy.optimize import curve_fit, OptimizeWarning, least_squares, approx_fprime
 from timeit import default_timer as timer
+import itertools
 
 from calculate.convert_units import OoM
 from calculate.optimize.misc_fcns import generalized_loss_fcn
@@ -268,23 +269,27 @@ def calc_s(jac):
 
         return 1/x
 
-    def obj_fcn_decorator(fit_fcn, grad_fcn, x0, idx, T, ln_k, return_sum=True):
+    def obj_fcn_decorator(fit_fcn, fit_func_jac, x0, idx, T, ln_k, bnds, return_sum=True):
         def obj_fcn(x, grad=[]):
             #x = x/s + x0[idx]
 
-            #print(x)
+            #warnings.simplefilter('ignore', OptimizeWarning)
+            #x, _ = curve_fit(fit_func, T, ln_k, p0=x, method='trf', bounds=bnds, # dogbox
+            #                 jac=fit_func_jac, x_scale='jac', max_nfev=len(p0)*1000)
 
             resid = fit_func(T, *x) - ln_k
-            if return_sum:
+            if return_sum:              
                 obj_val = generalized_loss_fcn(resid).sum()
             else:
                 obj_val = resid
 
-            #s[:] = np.abs(np.sum(loss*grad_fcn(T, *x).T, axis=1))
+            #s[:] = np.abs(np.sum(loss*fit_func_jac(T, *x).T, axis=1))
+
+            if grad.size > 0:
+                grad[:] = np.sum(fit_func_jac(T, *x).T*resid, axis=1)
 
-            #if len(grad) > 0:
-            #    grad[:] = np.sum(loss*grad_fcn(T, *x).T, axis=1)
             print(obj_val)
+
             return obj_val
         return obj_fcn
 
@@ -314,7 +319,10 @@ def obj_fcn(x, grad=[]):
     # initial guesses for falloff
     #if len(x0) != 10:
     x0 = [*x0[:6], 0.7, 200, 300, 400] # initial guesses for fitting Troe if none exist
-
+
+    #x0_falloff = list(itertools.product([0, 1], repeat=4))
+    #print(x0_falloff)
+
     x0 = np.array(x0)
 
     A_idx = [1, 4]
@@ -371,16 +379,17 @@ def obj_fcn(x, grad=[]):
             warnings.simplefilter('ignore', OptimizeWarning)
             try:
                 x_fit, _ = curve_fit(fit_func, T, ln_k, p0=p0, method='trf', bounds=bnds, # dogbox
-                                        #jac=fit_func_jac, x_scale='jac', max_nfev=len(p0)*1000)
-                                        jac='2-point', x_scale='jac', max_nfev=len(p0)*1000)
+                                        jac=fit_func_jac, x_scale='jac', max_nfev=len(p0)*1000)
+                                        #jac='2-point', x_scale='jac', max_nfev=len(p0)*1000)
             except:
                 return
 
     else:
         #s[:] = calc_s(fit_func_jac(T, *p0))
-        nlopt_fit_fcn = obj_fcn_decorator(fit_func, fit_func_jac, x0, idx, T, ln_k)
+        nlopt_fit_fcn = obj_fcn_decorator(fit_func, fit_func_jac, x0, idx, T, ln_k, bnds)
 
         opt = nlopt.opt(nlopt.GN_CRS2_LM, len(idx))
+        #opt = nlopt.opt(nlopt.GN_DIRECT, len(idx))
         #opt = nlopt.opt(nlopt.LN_SBPLX, len(idx)) # nlopt.LN_SBPLX nlopt.LN_COBYLA nlopt.LD_MMA nlopt.LD_LBFGS
 
         opt.set_min_objective(nlopt_fit_fcn)
@@ -394,7 +403,7 @@ def obj_fcn(x, grad=[]):
         opt.set_lower_bounds(bnds[0][idx])
         opt.set_upper_bounds(bnds[1][idx])
 
-        #opt.set_initial_step(calc_s(fit_func_jac(T, *p0))*1E-8)
+        opt.set_initial_step(calc_s(fit_func_jac(T, *p0))*1E-2)
         x_fit = opt.optimize(p0) # optimize!
 
     if HPL_LPL_defined: # Fit falloff

src/calculate/optimize/misc_fcns.py

@@ -10,6 +10,7 @@
 min_pos_system_value = np.finfo(float).eps*(1E2)
 max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/3)
 min_neg_system_value = -max_pos_system_value
+T_max = 6000
 
 default_arrhenius_coefNames = ['activation_energy', 'pre_exponential_factor', 'temperature_exponent']
 
@@ -141,7 +142,7 @@ def set_bnds(mech, rxnIdx, keys, coefNames):
                 if coef_x0 > 0:
                     coef_bnds['lower'].append(0)                                # Ea shouldn't change sign
                 else:
-                    coef_bnds['lower'].append(-Ru*10000*np.log(max_pos_system_value))
+                    coef_bnds['lower'].append(-Ru*T_max*np.log(max_pos_system_value))
             elif coefName == 'pre_exponential_factor':
                 coef_bnds['lower'].append(min_pos_system_value)             # A should be positive
             elif not isinstance(coefName, int):     # ints will be falloff, they will be taken care of below
@@ -151,7 +152,7 @@ def set_bnds(mech, rxnIdx, keys, coefNames):
             if coefName == 'activation_energy' and coef_x0 < 0:   # Ea shouldn't change sign
                 coef_bnds['upper'].append(0)
             elif coefName == 'temperature_exponent':
-                coef_bnds['upper'].append(np.log(max_pos_system_value)/np.log(298.15))
+                coef_bnds['upper'].append(np.log(max_pos_system_value)/np.log(T_max))
             elif not isinstance(coefName, int):
                 coef_bnds['upper'].append(max_pos_system_value)
         else: