diff --git a/frhodo_environment.yml b/frhodo_environment.yml index ad9723a..7012f94 100644 --- a/frhodo_environment.yml +++ b/frhodo_environment.yml @@ -4,12 +4,15 @@ channels: - conda-forge - defaults dependencies: - - cantera=2.5.0b1 - - matplotlib=3.3.1 - - nlopt=2.6.1 - - numpy=1.19.1 - - pyside2=5.13.1 + - cantera=2.6.0a1 + - matplotlib=3.4.2 + - nlopt=2.7.0 + - numpy=1.21.1 + - pip + - pyside2=5.13.2 - qtpy=1.9.0 - - requests=2.24.0 - - scipy=1.5.2 - - tabulate=0.8.7 \ No newline at end of file + - requests=2.26.0 + - scipy=1.7.0 + - tabulate=0.8.9 + - pip: + - dtcwt==0.12.0 \ No newline at end of file diff --git a/src/calculate/optimize/fit_coeffs.py b/src/calculate/optimize/fit_coeffs.py index 2cea195..b443daf 100644 --- a/src/calculate/optimize/fit_coeffs.py +++ b/src/calculate/optimize/fit_coeffs.py @@ -13,6 +13,7 @@ from calculate.convert_units import OoM from calculate.optimize.misc_fcns import generalized_loss_fcn, set_arrhenius_bnds +from calculate.smooth_data import dual_tree_complex_wavelet_filter Ru = ct.gas_constant # Ru = 1.98720425864083 @@ -27,6 +28,10 @@ 'activation_energy_inf', 'pre_exponential_factor_inf', 'temperature_exponent_inf', 'A', 'T3', 'T1', 'T2'] +#falloff_bnds = np.array([[1E-12, 1E-12, 0.1], [1E30, 1E30, 1]]) +#falloff_bnds[:,0:2] = np.log(falloff_bnds[:,0:2]) # ln(k_0), ln(k_inf), Fcent +falloff_bnds = np.log([[1E-6, 1E-6, 1E-3], [1E25, 1E25, 1]]) # ln(k_0), ln(k_inf), ln(Fcent) + troe_falloff_0 = [[0.6, 200, 600, 1200], # (0, 0, 0) # (0, 0, 1) [0.05, 1000, -2000, 3000], # (0, 1, 0) @@ -37,8 +42,6 @@ # (1, 1, 1) ] -falloff_bnds = np.log([[1E-12, 1E-12, 0.1], [1E30, 1E30, 1]]) # ln(k_0), ln(k_inf), ln(Fcent) - #troe_min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/3) #troe_max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/3) #troe_min_neg_system_value = -max_pos_system_value @@ -112,11 +115,7 @@ def ln_arrhenius_jac(T, *args): bnds[1][A_idx] = np.log(bnds[1][A_idx]) # only valid initial guesses - for n, val in enumerate(p0): - if val < bnds[0][n]: - p0[n] = bnds[0][n] - elif val > bnds[1][n]: - p0[n] = bnds[1][n] + p0 = np.clip(p0, bnds[0], bnds[1]) with warnings.catch_warnings(): warnings.simplefilter('ignore', OptimizeWarning) @@ -159,6 +158,7 @@ def x_bnds(self, x0): bnds.append(troe_all_bnds[coef]['-']) else: bnds.append(troe_all_bnds[coef]['+']) + return np.array(bnds).T def fit(self): @@ -565,6 +565,21 @@ def fit_arrhenius_rates(self): def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTANT arrhenius_coeffs, T, ln_k = self.fit_arrhenius_rates() + # set nlopt opt parameters + options = [{'algorithm': nlopt.GN_DIRECT_NOSCAL, 'xtol_rel': 1E-4, 'ftol_rel': 1E-4, # GN_DIRECT_NOSCAL GN_DIRECT + 'initial_step': 1, 'max_eval': 2000}, + #{'algorithm': nlopt.LN_COBYLA, 'xtol_rel': 1E-2, 'ftol_rel': 1E-2, + # 'initial_step': 1E-3, 'max_eval': 500}, + #{'algorithm': nlopt.LN_COBYLA, 'xtol_rel': 1E-2, 'ftol_rel': 1E-2, + # 'initial_step': 1E-4, 'max_eval': 1000} # LN_COBYLA + ] + + # set scipy curve-fit loss + if self.loss_fcn_par == [1, 1]: # huber-like + loss = 'huber' # 'soft_l1', 'huber', 'cauchy', 'arctan' + else: + loss = 'linear' + # set T, P, M as arrays P = np.unique(self.P) P[1:] = np.roll(P[1:], 1) @@ -583,96 +598,122 @@ def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTA if self.P_limit_constrained[i]: del idx_fit[i] - if self.loss_fcn_par == [1, 1]: # huber-like - loss = 'huber' # 'soft_l1', 'huber', 'cauchy', 'arctan' - else: - loss = 'linear' - res = [] - for i in range(0, np.shape(T)[1]): # keep T constant and fit log_k_0, log_k_inf, log_Fcent - # set fit function, might be a better way, but this works for now - if idx_fit == [0, 1, 2]: - fit_fcn = self.ln_Troe - M_fit = M[:,i] - ln_k_fit = ln_k[:,i] - - elif idx_fit == [0, 2]: - n = np.concatenate([[0], np.arange(2, len(M[:,i]))]) - M_fit = M[n,i] - ln_k_fit = ln_k[n,i] - - fit_fcn = lambda M, ln_k_0, Fcent: self.ln_Troe(M, [ln_k_0, ln_k[1,i], Fcent]) - - elif idx_fit == [1, 2]: - n = np.concatenate([[1], np.arange(2, len(M[:,i]))]) - M_fit = M[n,i] - ln_k_fit = ln_k[n,i] - - fit_fcn = lambda M, ln_k_inf, Fcent: self.ln_Troe(M, [ln_k[0,i], ln_k_inf, Fcent]) - - elif idx_fit == [2]: - n = np.arange(2, len(M[:,i])) - M_fit = M[n,i] - ln_k_fit = ln_k[n,i] - - fit_fcn = lambda M, Fcent: self.ln_Troe(M, [ln_k[0,i], ln_k[1,i], Fcent]) - - p_bnds = np.copy(falloff_bnds)[:, idx_fit] # ln(k_0), ln(k_inf), ln(Fcent) - - if len(res) == 0: - p0 = np.array([20, 14, -0.5])[idx_fit] # ln(k_0), ln(k_inf), ln(Fcent) - else: - p0 = np.log([res[-1][1:][i] for i in idx_fit]) - - # TODO: Change this to nlopt? - if self.scipy_curvefit: - x_fit, _ = curve_fit(fit_fcn, M_fit, ln_k_fit, p0=p0, method='trf', bounds=p_bnds, # dogbox - jac='2-point', x_scale='jac', max_nfev=len(p0)*1000, loss=loss) - else: - self.obj_vars = {'fcn': fit_fcn, 'x': M_fit, 'y': ln_k_fit} - - xtol_rel = 1E-4 - ftol_rel = 1E-4 - initial_step = 1E-8 + for n, opt_options in enumerate(options): + if n > 0: # smooth prior optimization for initial guess in next opt + res = np.log(res) + for i in [1, 2, 3]: + res[:, i] = dual_tree_complex_wavelet_filter(res[:, i], lvls=2) # smooth + + res = np.exp(res) + #print(f'smoothed res {n-1}') + #print(res) + + for i in np.arange(0, np.shape(T)[1]): # keep T constant and fit log_k_0, log_k_inf, log_Fcent + # set fit function, might be a better way, but this works for now + if idx_fit == [0, 1, 2]: + fit_fcn = self.ln_Troe + M_fit = M[:,i] + ln_k_fit = ln_k[:,i] + + elif idx_fit == [0, 2]: + j = np.concatenate([[0], np.arange(2, len(M[:,i]))]) + M_fit = M[j,i] + ln_k_fit = ln_k[j,i] + + fit_fcn = lambda M, ln_k_0, Fcent: self.ln_Troe(M, [ln_k_0, ln_k[1,i], Fcent]) + + elif idx_fit == [1, 2]: + j = np.concatenate([[1], np.arange(2, len(M[:,i]))]) + M_fit = M[j,i] + ln_k_fit = ln_k[j,i] + + fit_fcn = lambda M, ln_k_inf, Fcent: self.ln_Troe(M, [ln_k[0,i], ln_k_inf, Fcent]) + + elif idx_fit == [2]: + j = np.arange(2, len(M[:,i])) + M_fit = M[j,i] + ln_k_fit = ln_k[j,i] + + fit_fcn = lambda M, Fcent: self.ln_Troe(M, [ln_k[0,i], ln_k[1,i], Fcent]) + + # set p0 and p_bnds for optimization + if len(res) == 0: + p0 = np.array([20, 14, -0.5])[idx_fit] # ln(k_0), ln(k_inf), ln(Fcent) + #p0 = np.array([20, 14, 0.6])[idx_fit] # ln(k_0), ln(k_inf), Fcent + p_bnds = np.copy(falloff_bnds)[:, idx_fit] - self.opt = opt = nlopt.opt(nlopt.LN_SBPLX, len(idx_fit)) + else: + if n == 0: + p0 = np.log(res[i-1][1:]) + else: + p0 = np.log(res[i][1:]) + + p_bnds = np.copy(falloff_bnds) + if i > 0: + p_bnds[0][0:2] = np.log(res[i-1][1:])[0:2] # change lower rate bnds to be prior rate + + p0 = p0[idx_fit] + p_bnds = p_bnds[:, idx_fit] # ln(k_0), ln(k_inf), ln(Fcent) + + if len(p_bnds) > 0: # remove guesses outside of bounds + p0 = np.clip(p0, p_bnds[0], p_bnds[1]) + + if self.scipy_curvefit: # not setup after global/local change + x_fit, _ = curve_fit(fit_fcn, M_fit, ln_k_fit, p0=x0, method='trf', bounds=p_bnds, # dogbox + jac='2-point', x_scale='jac', max_nfev=len(x0)*1000, loss=loss) - opt.set_min_objective(self.objective) - opt.set_maxeval(5000) + else: + s = np.ones_like(p0) + self.obj_vars = {'fcn': fit_fcn, 'x': M_fit, 'y': ln_k_fit, 'coefs_0': p0, 's': s} + s = self.obj_vars['s'] = self.calc_s(np.zeros_like(p0)) - opt.set_xtol_rel(xtol_rel) - opt.set_ftol_rel(ftol_rel) + self.opt = opt = nlopt.opt(opt_options['algorithm'], len(idx_fit)) - opt.set_lower_bounds(p_bnds[0]) - opt.set_upper_bounds(p_bnds[1]) + opt.set_min_objective(self.objective) + opt.set_maxeval(opt_options['max_eval']) + opt.set_xtol_rel(opt_options['xtol_rel']) + opt.set_ftol_rel(opt_options['ftol_rel']) - opt.set_initial_step(initial_step) - #opt.set_population(int(np.rint(10*(len(idx)+1)*10))) + opt.set_lower_bounds((p_bnds[0]-p0)/s) + opt.set_upper_bounds((p_bnds[1]-p0)/s) - x_fit = opt.optimize(p0) # optimize! + initial_step = opt_options['initial_step'] + opt.set_initial_step(initial_step) + #opt.set_population(int(np.rint(10*(len(idx)+1)*10))) + x_fit = opt.optimize(np.zeros_like(p0)) # optimize! + x_fit = x_fit*s + p0 - if idx_fit == [0, 1, 2]: - res.append([T[0,i], *np.exp(x_fit)]) + #print(f'{opt.last_optimum_value()} in {opt.get_numevals()}') - elif idx_fit == [0, 2]: - res.append([T[0,i], *np.exp([x_fit[0], ln_k[1,i], x_fit[1]])]) + if idx_fit == [0, 1, 2]: + res_vector = [T[0,i], *np.exp(x_fit)] - elif idx_fit == [1, 2]: - res.append([T[0,i], *np.exp([ln_k[0,i], x_fit[0], x_fit[1]])]) + elif idx_fit == [0, 2]: + res_vector = [T[0,i], *np.exp([x_fit[0], ln_k[1,i], x_fit[1]])] - elif idx_fit == [2]: - res.append([T[0,i], *np.exp([ln_k[0,i], ln_k[1,i], x_fit[0]])]) + elif idx_fit == [1, 2]: + res_vector = [T[0,i], *np.exp([ln_k[0,i], x_fit[0], x_fit[1]])] - #cmp = np.array([T[idx], M[idx], ln_k[idx], fit_func(M[idx], *x_fit)]).T - #for entry in cmp: - # print(*entry) - #print('') + elif idx_fit == [2]: + res_vector = [T[0,i], *np.exp([ln_k[0,i], ln_k[1,i], x_fit[0]])] + + if n == 0: + res.append(res_vector) + else: + res[i] = res_vector + #cmp = np.array([T[idx], M[idx], ln_k[idx], fit_func(M[idx], *x_fit)]).T + #for entry in cmp: + # print(*entry) + #print('') + + #print(f'res {n}') + #print(np.array(res)) + res = np.array(res) - print(res) - # Fit HPL and LPL Arrhenius parameters for idx, arrhenius_type in enumerate(['low_rate', 'high_rate']): if idx in idx_fit: @@ -687,6 +728,8 @@ def ln_Troe(self, M, *x): if len(x) == 1: x = x[0] + #Fcent = x[2] + #k_0, k_inf = np.exp(x[:2]) [k_0, k_inf, Fcent] = np.exp(x) with np.errstate(all='raise'): try: @@ -695,21 +738,23 @@ def ln_Troe(self, M, *x): except: return np.ones_like(M)*max_pos_system_value - log10_Fcent = np.log10(Fcent) - C = -0.4 - 0.67*log10_Fcent - N = 0.75 - 1.27*log10_Fcent + log_Fcent = np.log10(Fcent) + C = -0.4 - 0.67*log_Fcent + N = 0.75 - 1.27*log_Fcent f1 = (log_P_r + C)/(N - 0.14*(log_P_r + C)) ln_F = np.log(Fcent)/(1 + f1**2) - ln_k_calc = np.log(k_inf*P_r/(1 + P_r)) + ln_F return ln_k_calc - def objective(self, coefs, grad=np.array([]), obj_type='obj_sum'): + def objective(self, coefs_in, grad=np.array([]), obj_type='obj_sum'): fit_fcn, x, y = self.obj_vars['fcn'], self.obj_vars['x'], self.obj_vars['y'] + coefs_0, s = self.obj_vars['coefs_0'], self.obj_vars['s'] loss_scale, loss_alpha = self.loss_fcn_par + coefs = coefs_in*s + coefs_0 + resid = fit_fcn(x, coefs) - y if obj_type == 'obj_sum': obj_val = generalized_loss_fcn(resid, a=loss_alpha, c=loss_scale).sum() @@ -725,6 +770,22 @@ def objective(self, coefs, grad=np.array([]), obj_type='obj_sum'): return obj_val + def calc_s(self, x, grad=[]): + if len(grad) == 0: + grad = approx_fprime(x, self.objective, 1E-10) + + y = np.abs(grad) + if (y < min_pos_system_value).all(): + y = np.ones_like(y)*1E-14 + else: + y[y < min_pos_system_value] = 10**(OoM(np.min(y[y>=min_pos_system_value])) - 1) # TODO: MAKE THIS BETTER running into problem when s is zero, this is a janky workaround + + s = 1/y + #s = s/np.min(s) + s = s/np.max(s) + + return s + def fit_generic(rates, T, P, X, rxnIdx, coefKeys, coefNames, is_falloff_limit, mech, bnds, mpPool=None): rxn = mech.gas.reaction(rxnIdx) diff --git a/src/calculate/optimize/mech_optimize.py b/src/calculate/optimize/mech_optimize.py index 357a6f0..e629b71 100644 --- a/src/calculate/optimize/mech_optimize.py +++ b/src/calculate/optimize/mech_optimize.py @@ -185,7 +185,7 @@ def _set_coef_opt(self): return coef_opt - def _set_rxn_coef_opt(self, min_T_range=600, min_P_range_factor=2): + def _set_rxn_coef_opt(self, min_T_range=1000, min_P_range_factor=2): mech = self.parent.mech rxn_coef_opt = [] for coef in self.coef_opt: @@ -226,7 +226,7 @@ def _set_rxn_coef_opt(self, min_T_range=600, min_P_range_factor=2): P_bnds = np.array([P_min, P_max]) if P_bnds[1]/P_bnds[0] < min_P_range_factor: - P_f_min = min_P_range_factor/2 + 1 + P_f_min = min_P_range_factor**0.5 P_bnds = np.array([P_median/P_f_min, P_median*P_f_min]) for rxn_coef in rxn_coef_opt: @@ -308,8 +308,7 @@ def _set_rxn_rate_opt(self): rxn_rate_opt['x0'] = rates(self.rxn_coef_opt, mech) # Determine rate bounds - lb = [] - ub = [] + bnds = np.array([[],[]]) for i, rxn_coef in enumerate(self.rxn_coef_opt): rxnIdx = rxn_coef['rxnIdx'] rxn = mech.gas.reaction(rxnIdx) @@ -338,18 +337,13 @@ def _set_rxn_rate_opt(self): rxn_rate_opt['x0'][i+n] = np.log(x[1]) + x[2]*np.log(T) - x[0]/(Ru*T) ln_rate = rxn_rate_opt['x0'][i:i + len(rxn_coef['T'])] - bnds = mech.rate_bnds[rxnIdx]['limits'](np.exp(ln_rate)) - bnds = np.log(bnds) # operate on ln and scale - scaled_bnds = np.sort(bnds/ln_rate, axis=0) + rxn_coef_bnds = mech.rate_bnds[rxnIdx]['limits'](np.exp(ln_rate)) + rxn_coef_bnds = np.log(rxn_coef_bnds) # operate on ln and scale + scaled_rxn_coef_bnds = np.sort(rxn_coef_bnds/ln_rate, axis=0) - if np.size(bnds) > 2: - lb.extend(list(scaled_bnds[0])) - ub.extend(list(scaled_bnds[1])) - else: - lb.append(scaled_bnds[0]) - ub.append(scaled_bnds[1]) - - rxn_rate_opt['bnds'] = {'lower': np.array(lb), 'upper': np.array(ub)} + bnds = np.concatenate((bnds, scaled_rxn_coef_bnds), axis=1) + + rxn_rate_opt['bnds'] = {'lower': bnds[0,:], 'upper': bnds[1,:]} # set mech to prior mech mech.coeffs = prior_coeffs