import sys import numpy as np from numpy import exp, log, log10, sin, cos, tan, arcsin, arccos, arctan, sinh, cosh, tanh, sqrt, abs import openpyxl import pandas as pd import matplotlib.pyplot as plt 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.tkoptimize import mlsq_general from tklib.tktransport.tkTransport import integrate_Simpson_list, fe, fh, NC2meff_FEA from tklib.tktransport.tkmobility import tkMobility from tklib.tksci.tkFit import tkFit from tklib.tksci.tkFit_m import tkFit_m from tklib.tkgraphic.tkplotevent import tkPlotEvent from tklib.tksci.tksci import asin, acos, atan, degcos, degsin, degtan, degacos, degasin, degatan, arcsin, arccos, arctan from tklib.tksci.tksci import factorial, gamma, combination, eVTonm, nmToeV, Bn from tklib.tksci.tksci import Ea_Arrhenius from tklib.tksci.tksci import h, h_bar, hbar, e, kB, NA, c, pi, pi2, torad, todeg, basee from tklib.tksci.tksci import me, mp, mn from tklib.tksci.tksci import u0, e0, a0, R, F, g, G from tklib.tksci.tksci import HartreeToeV, RyToeV, KToeV, eVToK, JToeV, eVToJ, Debye from tklib.tksci.tkFit import tkFit from tklib.tksci.tkFit_m import tkFit_m norder = 3 infile = 'random-poly.xlsx' xmin = -1.0e100 xmax = 1.0e100 xcalmin = '*' xcalmax = '*' xlabel = 0 ylabel = 1 fontsize = 16 app = tkApplication() argv = sys.argv narg = len(argv) if narg >= 2: infile = argv[1] if narg >= 3: norder = pint(argv[2]) if narg >= 4: xlabel = argv[3] if narg >= 5: ylabel = argv[4] if narg >= 6: xmin = pfloat(argv[5]) if narg >= 7: xmax = pfloat(argv[6]) if narg >= 8: xcalmin = pfloat(argv[7], defval = xcalmin) if narg >= 9: xcalmax = pfloat(argv[8], defval = xcalmax) def main(): global xlabel, ylabel, xcalmin, xcalmax logfile = app.replace_path(infile) print(f"Open logfile [{logfile}]") app.redict(targets = ["stdout", logfile], mode = 'w') print("Least-squares method for polynomial order {}".format(norder)) print(f"infile: {infile}") print(f"xlabel: {xlabel}") print(f"ylabel: {ylabel}") print(f"norder: {norder}") print(f"Fitting data range: {xmin} - {xmax}") if xlabel == '': app.terminate("\nError: Null xlabel is not allowed\n", pause = True) if ylabel == '': app.terminate("\nError: Null ylabel is not allowed\n", pause = True) print("") print(f"Read [{infile}]") datafile = tkVariousData(infile) labels, datalist = datafile.Read_minimum_matrix(close_fp = True, force_numeric = False) _xlabel, xin = datafile.FindDataArray(xlabel, flag = 'i') _ylabel, yin = datafile.FindDataArray(ylabel, flag = 'i') if xin is None: app.terminate("\nError: xlabel [{xlabel} is not found\n", pause = True) if yin is None: app.terminate("\nError: ylabel [{ylabel} is not found\n", pause = True) x = [] y = [] for i in range(len(xin)): if xmin <= xin[i] <= xmax: x.append(xin[i]) y.append(yin[i]) ndata = len(x) print("") print("ndata=", ndata) print(f"{labels[0]:12} {labels[1]:12}") for i in range(ndata): print(f"{x[i]:12.4g} {y[i]:12.4g}") print("") print(f"Execute linear least-squares method") ci = np.polyfit(x, y, norder) print("LSQ function") print("f(x) = {}".format(ci[norder]), end = '') for i in range(1, norder+1): print(" + {} * x^{}".format(ci[norder-i], i), end = '') fit = tkFit() print("") print("Calculating fitting result") xcalmin = pfloat(xcalmin, defval = min(x)) xcalmax = pfloat(xcalmax, defval = max(x)) ncaldata = 201 xcalstep = (xcalmax - xcalmin) / (ncaldata - 1) xcal = [xcalmin + i * xcalstep for i in range(ncaldata)] print(f"Calculation range: {xcalmin} - {xcalmax}") print(f" Number of calculation points: {ncaldata}") ycal1 = np.poly1d(ci)(x) ycal2 = np.poly1d(ci)(xcal) fit.print_scores(heading = "\nScores between y(input) and y(fit)", y1 = yin, y2 = ycal1) """ print(f"{xlabel:12} {ylabel:12} {ylabel+'(fit)':12}") for i in range(ndata): print(f"{x[i]:12.4g} {y[i]:12.4g} {ycal[i]:12.4g}") """ outfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-fit.xlsx"]) print("") print(f"Save results to [{outfile}]") fit.to_excel(outfile, [xlabel, ylabel, f"{ylabel}(fit)", "", xlabel, f"{ylabel}(cal)"], [x, y, ycal1, [], xcal, ycal2]) #================================================================ # Plot #================================================================ fig, axes = plt.subplots(1, 2, figsize = (8, 6)) plot_event = tkPlotEvent(plt) input_data = axes[0].plot(xin, yin, label = 'input', linestyle = '', marker = 'o') fit_data = axes[0].plot(xcal, ycal2, label = 'fit', linestyle = '-') # axes[0].plot(x, ycal, label = 'fit', linestyle = '-') axes[0].tick_params(labelsize = fontsize) axes[0].set_xlabel(xlabel, fontsize = fontsize) axes[0].set_ylabel(ylabel, fontsize = fontsize) axes[0].legend(fontsize = fontsize) coeff_data = axes[1].plot(range(norder+1), [ci[norder-i] for i in range(norder+1)], label = 'coeff', linestyle = '', marker = 'o') axes[1].tick_params(labelsize = fontsize) axes[1].set_xlabel('$i$', fontsize = fontsize) axes[1].set_ylabel('$c_i$', fontsize = fontsize) axes[1].legend(fontsize = fontsize) plot_event.add_data({"label": "input", "plot_type": "plot", "axis": axes[0], "data": input_data}) plot_event.add_data({"label": "fit", "plot_type": "plot", "axis": axes[0], "data": fit_data}) plot_event.add_data({"label": "coefficient", "plot_type": "plot", "axis": axes[1], "data": coeff_data}) plot_event.register_event(fig, event = "button_press_event", callback = lambda event: plot_event.onclick(event)) plt.tight_layout() plt.pause(0.1) app.terminate("", pause = True) if __name__ == "__main__": main()