import os import sys import numpy as np from numpy import sin, cos, tan, pi, exp, log, sqrt from scipy.optimize import minimize import pandas as pd from matplotlib import pyplot as plt import matplotlib.widgets as wg from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg, analyze_varstr, save_csv from tklib.tkutils import is_exist, is_file, is_dir from tklib.tkapplication import tkApplication from tklib.tkparams import tkParams from tklib.tkvariousdata import tkVariousData from tklib.tksci.tkmatrix import make_matrix1 from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me, Ea_Arrhenius, log10 from tklib.tksci.tkoptimize import mlsq_general from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA from tklib.tktransport.tkmobility import tkMobility """ 半導体の散乱モデルの最小二乗法 """ usage_str = ''' "" f"(i) usage: python {sys.argv[0]} fit mode target infile xlabel ylabel" '''[1:-1] def initialize(): #================================ # Global variables #================================ app = tkApplication(usage_str = usage_str) argv, narg = app.get_argv() app.cparams = tkParams() cparams = app.cparams cparams.debug = 0 cparams.print_level = 0 # mode: 'init', 'fit', 'sim' # cparams.mode = 'init' cparams.mode = 'fit' # cparams.target = 'Aop' # cparams.target = 'A_T0' cparams.target = 'A_T32' # cparams.target = 'VB' cparams.fix_N = False # data file to fit (.xlsx) cparams.infile = 'IO.xlsx' cparams.parameterfile = None cparams.parameterbkfile = None cparams.xlabel = 2 cparams.ylabel = 6 cparams.Nlabel = 8 # Temperature range cparams.Tmin = 50 # K cparams.Tmax = 1000 cparams.Tstep = 10 cparams.nT = int((cparams.Tmax - cparams.Tmin) / cparams.Tstep + 1) # 最適化パラメータの初期値 # cparams.varname = ("Eb", "Aop", "Eop", "a0", "p0", "a1", "p1", "a2", "p2", "a3", "p3", "a4", "p4") # cparams.varunit = ["eV", "", "eV", "", " ", "", "", "", "", "", "", "", ""] # cparams.ai0 = mu.parameter_list() # cparams.optid = [ 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0] # cparams.id_llsq = [ 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0] #============================================= # scipy.optimize.minimizeで使うアルゴリズム #============================================= #nelder-mead Downhill simplex #powell Modified Powell #cg conjugate gradient (Polak-Ribiere method) #bfgs BFGS法 #newton-cg Newton-CG #trust-ncg 信頼領域 Newton-CG 法 #dogleg 信頼領域 dog-leg 法 #L-BFGS-B’ (see here) #TNC’ (see here) #COBYLA’ (see here) #SLSQP’ (see here) #trust-constr’(see here) #dogleg’ (see here) #trust-exact’ (see here) #trust-krylov’ (see here) cparams.method = "nelder-mead" cparams.nmaxiter = 300 cparams.tol = 1.0e-3 cparams.h_diff = 1.0e-3 cparams.print_interval = 5 cparams.plot_interval = 1 #============================= # Graph configuration #============================= cparams.fig = None cparams.figsize = [6, 6] cparams.fontsize = 14 cparams.legend_fontsize = 12 cparams.graphupdateinterval = 10 return app, cparams #============================= # Treat argments #============================= def usage(app): cparams = app.cparams for s in app.usage_str.split('\n'): cmd = 'print({})'.format(s.rstrip()) eval(cmd) def update_vars(app, cparams): argv = app.argv # if len(argv) <= 1: # app.terminate(usage = usage) cparams.mode = getarg(1, cparams.mode) cparams.target = getarg(2, cparams.target) cparams.infile = getarg(3, cparams.infile) cparams.xlabel = getarg(4, cparams.xlabel) cparams.ylabel = getarg(5, cparams.ylabel) cparams.Nlabel = getarg(6, cparams.Nlabel) #Aop app.Aop_params = tkParams() Aop_params = app.Aop_params Aop_params.Eop_0 = 0.046 Aop_params.A1 = 0.014301489 Aop_params.B1 = 20.0591851 Aop_params.C1 = 0.0 Aop_params.logN01 = 19.57620089 #AT0 app.AT0_params = tkParams() AT0_params = app.AT0_params AT0_params.A2 = 0.013524202 AT0_params.B2 = 16.15118657 AT0_params.C2 = 0.006474685 AT0_params.logN02 = 18.99367178 #A_T3/2 app.AT32_params = tkParams() AT32_params = app.AT32_params AT32_params.A3 = 8.00E-36 AT32_params.B3 = 15 AT32_params.C3 = 0.00E+00 AT32_params.D3 = 0.00E+00 AT32_params.F3 = 0 AT32_params.logN03 = 18.5 #VB app.VB_params = tkParams() VB_params = app.VB_params VB_params.A6 = 4.64E-02 VB_params.B6 = 10.79576843 VB_params.C6 = 0 VB_params.D6 = 0.030576323 VB_params.logN06 = 18.36771264 class tkFit(tkParams): def __init__(self, parameter_file = None, app = None, tol = 1.0e-5, nmaxiter = 100, print_interval = 1, plot_interval = 1, **args): super(tkParams, self).__init__(**args) # super(tkParams, self).__init__(parameter_file, app, args) # self.update(**args) self.varname = [] self.unit = [] self.pk = [] self.optid = [] self.tol = tol self.print_interval = print_interval self.nmaxiter = nmaxiter self.plot_interval = plot_interval self.stop_flag = False self.x = None self.y = None self.func = None self.iter = 0 self.f = None self.pk = [] def extract_parameters(self, pk = None, optid = None): if pk is None: pk = self.pk if optid is None: optid = self.optid optpk = [] for i in range(len(pk)): if optid[i]: optpk.append(pk[i]) return optpk def recover_parameters(self, optpk, allpk = None, optid = None, set_member = True): if allpk is None: allpk = self.pk if optid is None: optid = self.optid pk = [] c = 0 for i in range(len(allpk)): if optid[i]: pk.append(optpk[c]) c += 1 else: pk.append(allpk[i]) if set_member: self.pk = pk return pk def print_parameters(self, varname = None, unit = None, pk = None, dx = None, optid = None, f = None): if varname is None: varname = self.varname if unit is None: unit = self.unit if pk is None: pk = self.pk if optid is None: optid = self.optid for i in range(len(pk)): print(f" {varname[i]:>10}={pk[i]:12.6g} {unit[i]:10} (id={optid[i]})") if f is not None: print(" f=", f) def minimize_func(self, pk, x = None, y = None): if x is None: x = self.x if y is None: y = self.y recoverpk = self.recover_parameters(pk, set_member = False) f = 0.0 n = len(x) for i in range(n): ycal = self.func(x[i], *recoverpk) d = ycal - y[i] f += d * d # print("i=", i, self.x[i], self.y[i], ycal) return f def cal_ylist(self, pk, x = None): if x is None: x = self.x return [self.func(x[i], *pk) for i in range(len(x))] def cal_mulist(self, pk, T = None): if T is None: T = slef.T return [cal_mu_from_params(T[i], *pk) for i in range(len(T))] def read_data(self, infile, xlabel, ylabel, usage = None): print("") print(f"Read [{infile}]") datafile = tkVariousData(infile) labels, datalist = datafile.Read_minimum_matrix(close_fp = True, force_numeric = False, usage = usage) xlabel, xin = datafile.FindDataArray(xlabel, flag = 'i') ylabel, yin = datafile.FindDataArray(ylabel, flag = 'i') self.labels = labels self.datafile = datafile self.xlabel = xlabel self.ylabel = ylabel self.x = xin self.y = yin def button_click(self, e): self.stop_flag = True def initial_plot(self, axis, yini, fontsize): self.plt = plt axis.tick_params(labelsize = fontsize) self.input_data = axis.plot(self.x, self.y, label = 'input', linestyle = '', marker = 'o', markersize = 5.0) self.initial_data = axis.plot(self.x, yini, label = 'initial', linestyle = '-', linewidth = 0.5, color = 'red') # self.final_data = axis.plot(self.x, yini, label = 'final', linestyle = '-', linewidth = 1.0, color = 'blue') self.fit_data = axis.plot([], [], label = 'fit', linestyle = '-', linewidth = 0.5, color = 'blue') axis.set_xlabel(self.xlabel, fontsize = fontsize) axis.set_ylabel(self.ylabel, fontsize = fontsize) axis.legend(fontsize = fontsize) self.plt.tight_layout() self.plt.subplots_adjust(top = 0.90, bottom = 0.10) self.ax_button = plt.axes([0.15, 0.95, 0.10, 0.03]) self.stop_button = wg.Button(self.ax_button, 'stop', color = '#f8e58c', hovercolor = '#38b48b') self.stop_button.on_clicked(self.button_click) self.plt.pause(0.0001) self.window = plt.get_current_fig_manager().window self.axis = axis self.yini = yini def minimize(self, method = None, tol = None, nmaxiter = None): if method is None: method = self.method if tol is None: tol = self.tol if nmaxiter is None: nmaxiter = self.nmaxiter print("") print("Optimizing parameters:") optpk = self.extract_parameters() print(" optpk=", optpk) print("Minimize:") res = minimize(self.minimize_func, optpk, method = method, jac = '3-point', callback = self.callback, tol = tol, options = {'maxiter': nmaxiter, "disp": True}) recoverpk = self.recover_parameters(res.x) return recoverpk, res.fun, res.success, res def callback(self, pk): if self.stop_flag: return False w = self.plt.get_current_fig_manager().window if w != self.window: self.stop_flag = True return False recoverpk = self.recover_parameters(pk, set_member = False) if self.iter % self.print_interval == 0: print(f"iter: {self.iter}") n = len(recoverpk) for i in range(n): print(f" {self.varname[i]:10}: {recoverpk[i]:10.4g} {self.unit[i]}") f = self.minimize_func(pk) print(f" f={f:12.6g}") if self.iter % self.plot_interval == 0: ycal = self.cal_ylist(recoverpk) self.fit_data[0].set_data(self.x, ycal) plt.tight_layout() plt.subplots_adjust(top = 0.90, bottom = 0.10) plt.pause(0.0001) self.iter += 1 return True def find_nearest_value(x, y, x0): dx_prev = 1.0e100 for i in range(len(x)): dx = abs(x[i] - x0) if dx < dx_prev: v = y[i] dx_prev = dx return v def Aop(N, A1, B1, C1, logN01): logN = log10(N) # print("A=", N, A1, B1, C1, logN01, N, logN) return A1 / (1 + exp(-B1 * (logN - logN01))) + C1 def Eop(N, Eop_0): return Eop_0 def AT0(N, A2, B2, C2, logN02): logN = log10(N) return A2 /(1.0 + exp(B2 * (logN - logN02))) + C2 def AT32(N, A3, B3, C3, D3, F3, logN03): logN = log10(N) K1 = 1.0 + exp(B3 * (logN - logN03)) a1 = A3 * N * N / K1 K2 = 1.0 + exp(-B3 * (logN - logN03)) a2 = C3 / K2 / (1.0 + D3 * N) + F3 return a1 + a2 def VB(N, A6, B6, C6, D6, logN06): logN = log10(N) return A6 * (1.0 - D6 * logN) / (1.0 + exp(B6 * (logN - logN06))) + C6 def cal_mu_from_params(T, _Eop, _Aop, _AT0, _AT32, _VB): ekBT = e / kB / T muinv = _Aop / (exp(_Eop * ekBT) - 1.0) muinv += _AT0 muinv += _AT32 * pow(T, -1.5) mu = 1.0 / muinv * exp(-_VB * ekBT) return mu def pk2params(*pk): Aop_params = tkParams() Aop_params.Eop_0 = pk[0] Aop_params.A1 = pk[1] Aop_params.B1 = pk[2] Aop_params.C1 = pk[3] Aop_params.logN01 = pk[4] AT0_params = tkParams() AT0_params.A2 = pk[5] AT0_params.B2 = pk[6] AT0_params.C2 = pk[7] AT0_params.logN02 = pk[8] AT32_params = tkParams() AT32_params.A3 = pk[9] AT32_params.B3 = pk[10] AT32_params.C3 = pk[11] AT32_params.D3 = pk[12] AT32_params.F3 = pk[13] AT32_params.logN03 = pk[14] #VB VB_params = tkParams() VB_params.A6 = pk[15] VB_params.B6 = pk[16] VB_params.C6 = pk[17] VB_params.D6 = pk[18] VB_params.logN06 = pk[19] return Aop_params, AT0_params, AT32_params, VB_params def cal_mu_from_TN(T, N, *pk): ekBT = e / kB / T Aop_params, AT0_params, AT32_params, VB_params = pk2params(*pk) # print("Aop_params=", Aop_params.__dict__) # print("AT0_params=", AT0_params.__dict__) # print("AT32_params=", AT32_params.__dict__) # print("VB_params=", VB_params.__dict__) _Aop = Aop(N, Aop_params.A1, Aop_params.B1, Aop_params.C1, Aop_params.logN01) _Eop = Eop(N, Aop_params.Eop_0) _AT0 = AT0(N, AT0_params.A2, AT0_params.B2, AT0_params.C2, AT0_params.logN02) _AT32 = AT32(N, AT32_params.A3, AT32_params.B3, AT32_params.C3, AT32_params.D3, AT32_params.F3, AT32_params.logN03) _VB = VB(N, VB_params.A6, VB_params.B6, VB_params.C6, VB_params.D6, VB_params.logN06) # print("params=", _Aop, _Eop, _AT0, _AT32, _VB) return cal_mu_from_params(T, _Eop, _Aop, _AT0, _AT32, _VB) class tkFit_TN(tkFit): def __init__(self, parameter_file = None, app = None, Tbase = None, Nbases = None, fix_N = False, tol = 1.0e-5, nmaxiter = 100, print_interval = 1, plot_interval = 1, **args): super(tkFit, self).__init__(**args) self.app = app self.varname = [] self.unit = [] self.pk = [] self.optid = [] self.fix_N = fix_N self.Tbase = Tbase self.Nbases = Nbases self.tol = tol self.print_interval = print_interval self.nmaxiter = nmaxiter self.plot_interval = plot_interval self.stop_flag = False self.x = None self.y = None self.func = None self.iter = 0 self.f = None self.pk = [] def cal_ylist(self, pk, T = None, N = None): if T is None: T = self.T if N is None: N = self.N if self.fix_N: return [self.func(T[i], self.Nbases[i], *pk) for i in range(len(T))] else: return [self.func(T[i], N[i], *pk) for i in range(len(T))] def minimize_func(self, pk, T = None, N = None, y = None): if T is None: T = self.T if N is None: N = self.N if y is None: y = self.y recoverpk = self.recover_parameters(pk, set_member = False) f = 0.0 n = len(T) for i in range(n): if self.fix_N: ycal = self.func(T[i], self.Nbases[i], *recoverpk) else: ycal = self.func(T[i], N[i], *recoverpk) d = ycal - y[i] f += d * d # print("i=", i, self.x[i], self.y[i], ycal) return f def initial_plot(self, axis, yini, fontsize): self.plt = plt axis.tick_params(labelsize = fontsize) index = range(len(self.y)) self.input_data = axis.plot(index, self.y, label = 'input', linestyle = '', marker = 'o', markersize = 5.0) if yini is not None: self.initial_data = axis.plot(index, yini, label = 'initial', linestyle = '-', linewidth = 0.5, color = 'red') self.fit_data = axis.plot([], [], label = 'fit', linestyle = '-', linewidth = 0.5, color = 'blue') axis.set_xlabel('index', fontsize = fontsize) axis.set_ylabel(self.ylabel, fontsize = fontsize) axis.legend(fontsize = fontsize) self.plt.tight_layout() self.plt.subplots_adjust(top = 0.90, bottom = 0.10) self.ax_button = plt.axes([0.15, 0.95, 0.10, 0.03]) self.stop_button = wg.Button(self.ax_button, 'stop', color = '#f8e58c', hovercolor = '#38b48b') self.stop_button.on_clicked(self.button_click) self.plt.pause(0.0001) self.window = plt.get_current_fig_manager().window self.axis = axis self.yini = yini def callback(self, pk): if self.stop_flag: return False w = self.plt.get_current_fig_manager().window if w != self.window: self.stop_flag = True return False recoverpk = self.recover_parameters(pk, set_member = False) if self.iter % self.print_interval == 0: print(f"iter: {self.iter}") n = len(recoverpk) for i in range(n): print(f" {self.varname[i]:10}: {recoverpk[i]:10.4g} {self.unit[i]}") f = self.minimize_func(pk) print(f" f={f:12.6g}") if self.iter % self.plot_interval == 0: ycal = self.cal_ylist(recoverpk) index = range(len(self.y)) self.fit_data[0].set_data(index, ycal) plt.tight_layout() plt.subplots_adjust(top = 0.90, bottom = 0.10) plt.pause(0.0001) self.iter += 1 return True def fit_mu_TN(app): cparams = app.cparams cparams.method = 'cg' cparams.sample_label = 0 cparams.formula_label = 1 cparams.Tbase = 300.0 print("") print("Fitting to mu-T Hall data by nonlinear least-squares method") print("mode : ", cparams.mode) print("target : ", cparams.target) print("infile : ", cparams.infile) print("parameter file : ", cparams.parameterfile) print("xlabel : ", cparams.xlabel) print("ylabel : ", cparams.ylabel) print("Nlabel : ", cparams.Nlabel) print("method : ", cparams.method) print("tol : ", cparams.tol) print("nmaxiter: ", cparams.nmaxiter) fit = tkFit_TN(tol = cparams.tol, nmaxiter = cparams.nmaxiter, print_interval = cparams.print_interval, plot_interval = cparams.plot_interval) print("") print(f"Read data from [{cparams.infile}]") cparams.T_label = cparams.xlabel cparams.mu_label = cparams.ylabel fit.read_data(cparams.infile, cparams.T_label, cparams.mu_label, usage = usage) N_label, N_in = fit.datafile.FindDataArray(cparams.Nlabel, flag = 'i') sample_label, sample_in = fit.datafile.FindDataArray(cparams.sample_label, flag = 'i') formula_label, formula_in = fit.datafile.FindDataArray(cparams.formula_label, flag = 'i') print("") print(f"Extract data for {cparams.Tmin} <= T <= {cparams.Tmax}") fit.sample = [] fit.T = [] fit.mu = [] fit.N = [] fit.index = [] for i in range(len(fit.x)): if cparams.Tmin <= fit.x[i] <= cparams.Tmax: fit.sample.append(sample_in[i]) fit.T.append(fit.x[i]) fit.mu.append(fit.y[i]) fit.N.append(N_in[i]) fit.index.append(sample_in[i]) fit.Tlabel = fit.xlabel fit.mulabel = fit.ylabel fit.Nlabel = N_label fit.index_label = sample_label fit.x = fit.T fit.y = fit.mu ndata = len(fit.T) print("ndata=", ndata) print(f"{'sample':10}: {fit.Tlabel:10} {fit.Nlabel:10} {fit.mulabel:10}") for i in range(ndata): print(f"{sample_in[i]:10}: {fit.T[i]:10.4g} {fit.N[i]:10.4g} {fit.mu[i]:10.4g}") # Calculate Nbase at Tbase for different samples data_set = {} for i in range(len(fit.sample)): s = fit.sample[i] if data_set.get(s, None) is None: data_set[s] = [[fit.T[i]], [fit.N[i]]] else: data_set[s][0].append(fit.T[i]) data_set[s][1].append(fit.N[i]) # Nbase = find_nearest_value(fit.T, fit.N, cparams.Tbase) fit.Nbases = [] Nbases_dict = {} for i in range(len(fit.sample)): s = fit.sample[i] Nbase = np.interp([cparams.Tbase], data_set[s][0], data_set[s][1])[0] Nbases_dict[s] = Nbase fit.Nbases.append(Nbase) # print("s=", sample_in) # print("Nbase=", fit.Nbases) print("") if cparams.fix_N: fit.fix_N = True fit.Tbase = cparams.Tbase for s in Nbases_dict.keys(): print(f"sample [{s}]: N is fixed to N={Nbases_dict[s]:10.4g} cm-3 at T={cparams.Tbase} K") else: print(f"Nbase={cparams.Nbase:10.4g} cm-3 from T={cparams.Tbase} K") print(f"N values are taken from T-dependent Hall data") #============================================= # Build pk #============================================= fit.func = cal_mu_from_TN p = app.Aop_params p.B1 = 0 p.C1 = 0.0 fit.varname = ["Eop_0", "A1", "B1", "C1", "logN01"] fit.unit = ["eV", "", "", "", ""] fit.pk = [p.Eop_0, p.A1, p.B1, p.C1, p.logN01] # fit.optid = [0, 1, 1, 1, 1] fit.optid = [0, 1, 0, 0, 0] p = app.AT0_params fit.varname += ["A2", "B2", "C2", "logN02"] fit.unit += ["", "", "", ""] fit.pk += [p.A2, p.B2, p.C2, p.logN02] fit.optid += [1, 1, 1, 1] p = app.AT32_params fit.varname += ["A3", "B3", "C3", "D3", "F3", "logN03"] fit.unit += ["", "", "", "", "", ""] fit.pk += [p.A3, p.B3, p.C3, p.D3, p.F3, p.logN03] fit.optid += [1, 1, 1, 1, 1, 1] p = app.VB_params fit.varname += ["A6", "B6", "C6", "D6", "logN06"] fit.unit += ["", "", "", "", ""] fit.pk += [p.A6, p.B6, p.C6, p.D6, p.logN06] fit.optid += [1, 1, 1, 1, 1] print("pk=", len(fit.pk), fit.pk) print("optid=", len(fit.optid), fit.optid) yini = fit.cal_ylist(fit.pk, fit.x) print("") print("Initial guess") print(f"{'i':4} {fit.Nlabel:10} {fit.Tlabel:10} {fit.mulabel:12} {fit.mulabel+'(cal)':12}") for i in range(len(fit.index)): print(f"{i:4} {fit.N[i]:10.4g} {fit.T[i]:10.4g} {fit.mu[i]:12.4g} {yini[i]:12.4g}") #======================================== # plot #======================================== fig, axes = plt.subplots(1, 1, figsize = cparams.figsize) axes = [axes] # axes = axes.flatten() fit.initial_plot(axes[0], yini, cparams.fontsize) #======================================== # Optimize #======================================== pfin, ffin, success, res = fit.minimize(cparams.method) if success: # print("") # print(res) print("") print(f"Converged at iteration: {res.nit}") else: print(f"Function did not converge") print(res) pk = fit.recover_parameters(pfin, set_member = True) print("Final parameters") n = len(pfin) for i in range(n): print(f" {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}") print(f" f={ffin:12.6g}") #======================================== # Final result #======================================== yfin = fit.cal_ylist(pfin) print("") print("Final result") print(f"{'i':4} {fit.Nlabel:10} {fit.Tlabel:10} {fit.mulabel:12} {fit.mulabel+'(cal)':12}") for i in range(len(fit.index)): print(f"{i:4} {fit.N[i]:10.4g} {fit.T[i]:10.4g} {fit.mu[i]:12.4g} {yini[i]:12.4g}") outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-mu_TN-fit.xlsx"]) df = pd.DataFrame(np.array([fit.N, fit.T, fit.mu, yini, yfin]).T, columns = ['N(cm-3)', 'T(K)', 'mu', 'mu(ini)', 'mu(fin)']) print("") print(f"Save results to [{outfile}]") df.to_excel(outfile, index = False) print("") print(f"Save parameters to {cparams.parameterfile}") cparams.save_parameters (cparams.parameterfile, section = 'Preferences', sort_by_keys = False, update_commandline = False, IsPrint = False) print("") print(f"Scattering parameters at Tbase={cparams.Tbase:10.4g} K") Aop_params, AT0_params, AT32_params, VB_params = pk2params(*pfin) for s in Nbases_dict.keys(): Nbase = Nbases_dict[s] mu_params = tkParams() mu_params.Aop = Aop(Nbase, Aop_params.A1, Aop_params.B1, Aop_params.C1, Aop_params.logN01) mu_params.Eop = Eop(Nbase, Aop_params.Eop_0) mu_params.AT0 = AT0(Nbase, AT0_params.A2, AT0_params.B2, AT0_params.C2, AT0_params.logN02) mu_params.AT32 = AT32(Nbase, AT32_params.A3, AT32_params.B3, AT32_params.C3, AT32_params.D3, AT32_params.F3, AT32_params.logN03) mu_params.VB = VB(Nbase, VB_params.A6, VB_params.B6, VB_params.C6, VB_params.D6, VB_params.logN06) print(f" sample [{s}]: N is fixed to N={Nbases_dict[s]:10.4g} cm-3 at T={cparams.Tbase} K") print(f" Aop ={mu_params.Aop:10.4g}") print(f" Eop ={mu_params.Eop:10.4g}") print(f" AT0 ={mu_params.AT0:10.4g}") print(f" AT32={mu_params.AT32:10.4g}") print(f" VB ={mu_params.VB:10.4g}") mu_params.save_parameters(cparams.parameterfile, section = s, sort_by_keys = False, update_commandline = False, IsPrint = False) index = range(len(yfin)) fit.fit_data[0].set_data(index, yfin) plt.tight_layout() plt.subplots_adjust(top = 0.90, bottom = 0.10) plt.pause(0.0001) app.terminate("", usage = usage, pause = True) def fit_mu_T(app): cparams = app.cparams cparams.method = 'cg' cparams.target_sample = '02A-01' cparams.Tbase = 300.0 # mu_params.Aop = 0.0 # mu_params.Eop_0 = 0.046 # mu_params.AT0 = 0.0671 # mu_params.AT32 = 8.912 # mu_params.VB = 0.0214 cparams.sample_label = 0 cparams.formula_label = 1 cparams.T_label = 6 cparams.mu_label = 7 cparams.N_label = 8 mu_params = tkParams() Aop_params = app.Aop_params AT0_params = app.AT0_params AT32_params = app.AT32_params VB_params = app.VB_params print("") print("Fitting to mu-T Hall data by nonlinear least-squares method") print("mode : ", cparams.mode) print("target : ", cparams.target) print("infile : ", cparams.infile) print("parameter file : ", cparams.parameterfile) print("xlabel : ", cparams.xlabel) print("ylabel : ", cparams.ylabel) print("method : ", cparams.method) print("tol : ", cparams.tol) print("nmaxiter: ", cparams.nmaxiter) fit = tkFit(tol = cparams.tol, nmaxiter = cparams.nmaxiter, print_interval = cparams.print_interval, plot_interval = cparams.plot_interval) print("") print(f"Read data from [{cparams.infile}]") fit.read_data(cparams.infile, cparams.T_label, cparams.mu_label, usage = usage) sample_label, sample_in = fit.datafile.FindDataArray(cparams.sample_label, flag = 'i') formula_label, formula_in = fit.datafile.FindDataArray(cparams.formula_label, flag = 'i') N_label, N_in = fit.datafile.FindDataArray(cparams.N_label, flag = 'i') print("") print(f"Extract data for sample [{cparams.target_sample}] and {cparams.Tmin} <= T <= {cparams.Tmax}") fit.T = [] fit.mu = [] fit.N = [] for i in range(len(fit.x)): if sample_in[i] == cparams.target_sample and cparams.Tmin <= fit.x[i] <= cparams.Tmax: fit.T.append(fit.x[i]) fit.mu.append(fit.y[i]) fit.N.append(N_in[i]) fit.Tlabel = fit.xlabel fit.mulabel = fit.ylabel fit.Nlabel = N_label fit.x = fit.T fit.y = fit.mu ndata = len(fit.T) print("ndata=", ndata) print(f"{fit.Tlabel:10} {fit.mulabel:10} {fit.Nlabel:10}") for i in range(ndata): print(f"{fit.T[i]:10.4g} {fit.mu[i]:10.4g} {fit.N[i]:10.4g}") Nbase = find_nearest_value(fit.T, fit.N, cparams.Tbase) cparams.Nbase = np.interp([cparams.Tbase], fit.T, fit.N)[0] mu_params.Aop = Aop(cparams.Nbase, Aop_params.A1, Aop_params.B1, Aop_params.C1, Aop_params.logN01) mu_params.Eop = Eop(cparams.Nbase, Aop_params.Eop_0) mu_params.AT0 = AT0(cparams.Nbase, AT0_params.A2, AT0_params.B2, AT0_params.C2, AT0_params.logN02) mu_params.AT32 = AT32(cparams.Nbase, AT32_params.A3, AT32_params.B3, AT32_params.C3, AT32_params.D3, AT32_params.F3, AT32_params.logN03) mu_params.VB = VB(cparams.Nbase, VB_params.A6, VB_params.B6, VB_params.C6, VB_params.D6, VB_params.logN06) fit.varname = [ "Eop", "Aop", "AT0", "AT32", "VB"] fit.unit = [ "", "", "", "", "eV"] fit.pk = [mu_params.Eop, mu_params.Aop, mu_params.AT0, mu_params.AT32, mu_params.VB] fit.optid = [ 0, 1, 1, 1, 1] fit.func = lambda T, *pk: cal_mu_from_params(T, *pk) print("") print(f"Scattering parameters at Nbase={cparams.Nbase:10.4g} cm-3") fit.print_parameters() optpk = fit.extract_parameters() yini = fit.cal_mulist(fit.pk, fit.T) print("") print("Initial guess") print(f"{'i':4} {'T':10} {'mu':12} {'mu(cal)':12}") for i in range(len(fit.T)): print(f"{i:4} {fit.T[i]:10.4g} {fit.mu[i]:12.4g} {yini[i]:12.4g}") #======================================== # plot #======================================== fig, axes = plt.subplots(1, 1, figsize = cparams.figsize) axes = [axes] # axes = axes.flatten() fit.initial_plot(axes[0], yini, cparams.fontsize) #======================================== # Optimize #======================================== pfin, ffin, success, res = fit.minimize(cparams.method) if success: # print("") # print(res) print("") print(f"Converged at iteration: {res.nit}") else: print(f"Function did not converge") print(res) pk = fit.recover_parameters(optpk, set_member = True) print("Final parameters") n = len(pfin) for i in range(n): print(f" {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}") print(f" f={ffin:12.6g}") #======================================== # Final result #======================================== yfin = fit.cal_ylist(pfin) print("") print("Final result") print(f"{'i':4} {'T':10} {'mu':12} {'mu(ini)':12} {'mu(fin)':12}") for i in range(len(fit.x)): print(f"{i:4} {fit.T[i]:10.4g} {fit.mu[i]:12.4g} {yini[i]:12.4g} {yfin[i]:12.4g}") outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-mu-fit.xlsx"]) df = pd.DataFrame(np.array([fit.T, fit.mu, yini, yfin]).T, columns = ['T(K)', 'mu', 'mu(ini)', 'mu(fin)']) print("") print(f"Save results to [{outfile}]") df.to_excel(outfile, index = False) print("") print(f"Save parameters to {cparams.parameterfile}") cparams.save_parameters (cparams.parameterfile, section = 'Preferences', sort_by_keys = False, update_commandline = False, IsPrint = False) mu_params.save_parameters (cparams.parameterfile, section = 'mu', sort_by_keys = False, update_commandline = False, IsPrint = False) fit.fit_data[0].set_data(fit.x, yfin) plt.tight_layout() plt.subplots_adjust(top = 0.90, bottom = 0.10) plt.pause(0.0001) app.terminate("", usage = usage, pause = True) def fit_parameters(app): cparams = app.cparams print("") print("Fitting to mu-T Hall data by nonlinear least-squares method") print("mode : ", cparams.mode) print("target : ", cparams.target) print("infile : ", cparams.infile) print("xlabel : ", cparams.xlabel) print("ylabel : ", cparams.ylabel) print("method : ", cparams.method) fit = tkFit(tol = cparams.tol, nmaxiter = cparams.nmaxiter, print_interval = cparams.print_interval, plot_interval = cparams.plot_interval) print("tol : ", fit.tol) print("nmaxiter : ", fit.nmaxiter) print("print_interval: ", fit.print_interval) print("plot_interval : ", fit.plot_interval) Aop_params = app.Aop_params AT0_params = app.AT0_params AT32_params = app.AT32_params VB_params = app.VB_params fit.Aop_varname = [ "A1", "B1", "C1", "log(N)_01"] fit.Aop_unit = [ "", "", "", ""] fit.Aop_pk = [Aop_params.A1, Aop_params.B1, Aop_params.C1, Aop_params.logN01] #Eop fit.Aop_optid = [ 1, 1, 1, 1] fit.AT0_varname = [ "A2", "B2", "C2", "log(N)_02"] fit.AT0_unit = [ "", "", "", ""] fit.AT0_pk = [AT0_params.A2, AT0_params.B2, AT0_params.C2, AT0_params.logN02] fit.AT0_optid = [ 1, 1, 1, 1] fit.AT32_varname = [ "A3", "B3", "C3", "D3", "F3", "log(N)_03"] fit.AT32_unit = [ "", "", "", "", "", ""] fit.AT32_pk = [AT32_params.A3, AT32_params.B3, AT32_params.C3, AT32_params.D3, AT32_params.F3, AT32_params.logN03] fit.AT32_optid = [ 1, 1, 1, 1, 1, 1] fit.VB_varname = [ "A6", "B6", "C6", "D6", "log(N)_06"] fit.VB_unit = [ "", "", "", "", ""] fit.VB_pk = [VB_params.A6, VB_params.B6, VB_params.C6, VB_params.D6, VB_params.logN06] fit.VB_optid = [ 1, 1, 1, 1, 1] if cparams.target == 'Aop': fit.varname = fit.Aop_varname fit.unit = fit.Aop_unit fit.pk = fit.Aop_pk fit.optid = fit.Aop_optid fit.func = Aop elif cparams.target == 'A_T0': fit.varname = fit.AT0_varname fit.unit = fit.AT0_unit fit.pk = fit.AT0_pk fit.optid = fit.AT0_optid fit.func = AT0 elif cparams.target == 'A_T32': fit.varname = fit.AT32_varname fit.unit = fit.AT32_unit fit.pk = fit.AT32_pk fit.optid = fit.AT32_optid fit.func = AT32 elif cparams.target == 'VB': fit.varname = fit.VB_varname fit.unit = fit.VB_unit fit.pk = fit.VB_pk fit.optid = fit.VB_optid fit.func = VB print("") print(f"Read data from [{cparams.infile}]") fit.read_data(cparams.infile, cparams.xlabel, cparams.ylabel, usage = usage) ndata = len(fit.x) print("x=", fit.xlabel, fit.x) print("y=", fit.ylabel, fit.y) print("ndata=", ndata) print("pk=", fit.pk) yini = fit.cal_ylist(fit.pk, fit.x) print("") print("Initial guess") print(f"{'i':4} {cparams.xlabel:10} {cparams.ylabel:12} {cparams.xlabel+'(cal)':12}") for i in range(len(fit.x)): print(f"{i:4} {fit.x[i]:10.4g} {fit.y[i]:12.4g} {yini[i]:12.4g}") # plot fig, axes = plt.subplots(1, 1, figsize = cparams.figsize) axes = [axes] # axes = axes.flatten() fit.initial_plot(axes[0], yini, cparams.fontsize) pfin, ffin, success, res = fit.minimize(cparams.method, cparams.tol, cparams.nmaxiter) if success: # print("") # print(res) print("") print(f"Converged at iteration: {res.nit}") else: print(f"Function did not converge") print(res) print("Final parameters") n = len(pfin) for i in range(n): print(f" {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}") print(f" f={ffin:12.6g}") yfin = fit.cal_ylist(pfin) print("") print("Final result") print(f"{'i':4} {fit.xlabel:10} {fit.ylabel:12} {fit.ylabel + '(ini)':12} {'y(fin)':12}") for i in range(len(fit.x)): print(f"{i:4} {fit.x[i]:10.4g} {fit.y[i]:12.4g} {yini[i]:12.4g} {yfin[i]:12.4g}") outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-{target}-fit.xlsx"], ext_dict = { "target": cparams.target } ) df = pd.DataFrame(np.array([fit.x, fit.y, yini, yfin]).T, columns = [cparams.xlabel, cparams.ylabel, f"{cparams.ylabel}(ini)", f"{cparams.ylabel}(fin)"]) print("") print(f"Save results to [{outfile}]") df.to_excel(outfile, index = False) print("") print(f"Save parameters to {cparams.parameterfile}") cparams.save_parameters (cparams.parameterfile, section = 'Preferences', sort_by_keys = False, update_commandline = False, IsPrint = False) Aop_params.save_parameters (cparams.parameterfile, section = 'Aop', sort_by_keys = False, update_commandline = False, IsPrint = False) AT0_params.save_parameters (cparams.parameterfile, section = 'A_T0', sort_by_keys = False, update_commandline = False, IsPrint = False) AT32_params.save_parameters(cparams.parameterfile, section = 'A_T3/2', sort_by_keys = False, update_commandline = False, IsPrint = False) VB_params.save_parameters (cparams.parameterfile, section = 'VB', sort_by_keys = False, update_commandline = False, IsPrint = False) fit.fit_data[0].set_data(fit.x, yfin) plt.tight_layout() plt.subplots_adjust(top = 0.90, bottom = 0.10) plt.pause(0.0001) app.terminate("\nPress ENTER to exit>>", usage = usage, pause = True) def main(): app, cparams = initialize() update_vars(app, cparams) cparams.logfile = app.replace_path(cparams.infile) cparams.outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-out.xlsx"]) print(f"Open logfile [{cparams.logfile}]") app.redict(targets = ["stdout", cparams.logfile], mode = 'w') cparams.parameterfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}.in"]) cparams.parameterbkfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-back.in"]) print("") print("infile : {}".format(cparams.infile)) print("parameter file : {}".format(cparams.parameterfile)) print("parameter backup file: {}".format(cparams.parameterbkfile)) print("mode : {}".format(cparams.mode)) print("target : {}".format(cparams.target)) print("") print("Read [{}]".format(cparams.parameterfile)) cparams.read_parameters(cparams.parameterfile) # check args again to use given parameters update_vars(app, cparams) if cparams.mode == 'init': init(app) elif cparams.mode == 'fit': if cparams.target == 'mu_T': fit_mu_T(app) elif cparams.target == 'mu_TN': fit_mu_TN(app) elif cparams.target == 'mu_T_Nfix': cparams.fix_N = True fit_mu_TN(app) elif cparams.target == 'Aop' or cparams.target == 'A_T0' or cparams.target == 'A_T32' or cparams.target == 'VB': fit_parameters(app) else: app.terminate("Error in main: Invalide target [{}]".format(cparams.target), usage = usage) elif cparams.mode == 'sim': sim(app) else: app.terminate("Error in main: Invalide mode [{}]".format(cparams.mode), usage = usage) app.terminate("", usage = usage, pause = True) if __name__ == "__main__": main()