import sys import argparse import numpy as np from numpy import exp, log import matplotlib.pyplot as plt from tklib.tkutils import replace_path, pint, pfloat from tklib.tkvariousdata import tkVariousData from tklib.tkapplication import tkApplication from tklib.tkparams import tkParams from tklib.tksci.tkFit import tkFit from tklib.tkgraphic.tkplotevent import tkPlotEvent from tklib.tksci.tksci import kB, e def initialize(): """ argparse の結果 args を受け取り、app.cfg に設定。 """ app = tkApplication() cfg = tkParams() parser = argparse.ArgumentParser(description="Arrhenius plot と多項式フィット") parser.add_argument('--infile', type=str, default="Hall-T.xlsx", help='Input file') parser.add_argument('--model', type=str, default='simple Arrhenius #log(P)=A-(eEa/kB)*(1/T)', help="モデル選択") parser.add_argument('--Tlabel', type=str, default='T(K)', help='T関連データ列ラベル') parser.add_argument('--Plabel', type=str, default='P', help='P関連データ列ラベル') parser.add_argument('--Ttype', type=str, choices=['T(K)', 'T(C)', '1/T', '1000/T'], default='T(K)', help='Tデータ変換方法') parser.add_argument('--Ptype', type=str, choices=['P', 'log10(P)', 'log_e(P)'], default='P', help='Pデータ変換方法') parser.add_argument('--xmin', type=float, default=-1.0e100, help='フィットする x の下限') parser.add_argument('--xmax', type=float, default=1.0e100, help='フィットする x の上限') parser.add_argument('--Tmin', type=float, default=-1.0e100, help='フィットする T の下限') parser.add_argument('--Tmax', type=float, default=1.0e100, help='フィットする T の上限') parser.add_argument('--Tcalmin', type=str, default='*', help="計算用する の下限（'*' で自動）") parser.add_argument('--Tcalmax', type=str, default='*', help="計算用する の上限（'*' で自動）") parser.add_argument('--ncal', type=int, default=201, help='Number of points in xcal grid') parser.add_argument('--xlsm_template', type=str, default="StandardGraph.xlsm", help='Excel template path') parser.add_argument('--plot_ci', type=int, default=1, help='Plot confidence interval') parser.add_argument('--plot_sigma_param', type=int, default=1, help='Plot parameter uncertainty band') parser.add_argument('--plot_sigma_pred', type=int, default=1, help='Plot prediction uncertainty band') parser.add_argument('--plot_sigma_combined', type=int, default=0, help='Plot combined uncertainty band') parser.add_argument('--figsize', type=float, nargs=2, default=[12, 8], help='プロット figsize, 例: --figsize 8 8') parser.add_argument('--fontsize', type=int, default=16, help='Font size for plots') parser.add_argument('--fontsize_legend', type=int, default=12, help='Legend font size') group = parser.add_mutually_exclusive_group() group.add_argument('--pause', dest='pause', action='store_true', help='Pause on terminate') group.add_argument('--no-pause', dest='pause', action='store_false', help='Do not pause on terminate') parser.set_defaults(pause=True) args = parser.parse_args() app.cfg = cfg for key in vars(args): setattr(cfg, key, getattr(args, key)) cfg.logfile = app.replace_path(cfg.infile) cfg.output_fitting_path = app.replace_path(cfg.infile, template = "{dirname}/{filebody}-fit.xlsm") cfg.output_parameter_path = app.replace_path(cfg.infile, template = ["{dirname}", "{filebody}-parameters.xlsx"]) return app, cfg, parser def build_design_matrix(x, order): x = np.asarray(x) N = len(x) X = np.ones((N, order + 1)) for j in range(1, order + 1): X[:, j] = [pow(val, j) for val in x] return X def mlsq_error(X, y): y = np.asarray(y) XtX = X.T @ X XtX_inv = np.linalg.inv(XtX) beta = XtX_inv @ (X.T @ y) # Residuals and variance residuals = y - X @ beta N, p = X.shape RSS = float(residuals.T @ residuals) # 修正箇所: N <= p の場合の例外処理[cite: 1] if N > p: sigma2_resid = RSS / (N - p) cov_beta = sigma2_resid * XtX_inv beta_std = np.sqrt(np.diag(cov_beta)) else: print(f"\n[Warning] Data points (N={N}) <= Parameter count (p={p}).") print("Error estimation is skipped to avoid divergence.") sigma2_resid = 0.0 cov_beta = None beta_std = np.zeros(p) return beta, beta_std, cov_beta, sigma2_resid def compute_param_uncertainty(X, cov_beta): return np.sum((X @ cov_beta) * X, axis=1) def compute_measurement_error(y): y = np.asarray(y) var_unbiased = np.var(y, ddof=1) if len(y) > 1 else 0.0 return np.sqrt(var_unbiased) def compute_bands(xcal, beta, cov_beta, sigma2_resid, sigma_meas): Xcal = build_design_matrix(xcal, len(beta) - 1) y_mean = Xcal @ beta # 修正箇所: cov_beta が None の場合は誤差を 0 にする[cite: 1] if cov_beta is not None: var_param = compute_param_uncertainty(Xcal, cov_beta) sigma_param = np.sqrt(var_param) sigma_pred = np.sqrt(var_param + sigma2_resid) sigma_combined = np.sqrt(var_param + sigma_meas**2) else: n = len(xcal) sigma_param = np.zeros(n) sigma_pred = np.zeros(n) sigma_combined = np.zeros(n) return { 'y_mean': y_mean, 'sigma_param': sigma_param, 'sigma_pred': sigma_pred, 'sigma_combined': sigma_combined } def execute(app): cfg = app.cfg app.redict(targets=["stdout", cfg.logfile], mode='w') cfg.xmin = pfloat(cfg.xmin) cfg.xmax = pfloat(cfg.xmax) cfg.Tmin = pfloat(cfg.Tmin) cfg.Tmax = pfloat(cfg.Tmax) datafile = tkVariousData(cfg.infile) labels, datalist = datafile.Read_minimum_matrix(close_fp=True, usage=app.usage) label_x, xX = datafile.FindDataArray(cfg.Tlabel, flag='i') label_y, yY = datafile.FindDataArray(cfg.Plabel, flag='i') if xX is None or yY is None: app.terminate("Error: 指定ラベルのデータが見つかりません", pause=cfg.pause) # T, P 変換処理 (省略せず維持) if cfg.Ttype == 'T(K)': T = xX elif cfg.Ttype == 'T(C)': T = [v + 273.15 for v in xX] elif cfg.Ttype == '1/T': T = [1.0 / v for v in xX] elif cfg.Ttype == '1000/T': T = [1000.0 / v for v in xX] else: app.terminate(f"Invalid Ttype [{cfg.Ttype}]", pause=cfg.pause) if cfg.Ptype == 'P': P = yY elif cfg.Ptype == 'log10(P)': P = [10**v for v in yY] elif cfg.Ptype == 'log_e(P)': P = [np.exp(v) for v in yY] else: app.terminate(f"Invalid Ptype [{cfg.Ptype}]", pause=cfg.pause) T1000 = [1000.0 / v for v in T] log10P = [log(v)/log(10.0) for v in P] x_fit, y_fit = [], [] for xi_orig, xi, yi in zip(xX, T1000, log10P): if cfg.xmin <= xi_orig <= cfg.xmax and cfg.Tmin <= (1000.0/xi if xi!=0 else float('inf')) <= cfg.Tmax: x_fit.append(xi) y_fit.append(yi) norder = {'percolation': 2, '3rd order': 3, '4th order': 4}.get(cfg.model, 1) X = build_design_matrix(x_fit, norder) beta, beta_std, cov_beta, sigma2_resid = mlsq_error(X, y_fit) # 修正箇所: 標準偏差の表示を条件分岐[cite: 1] print("Fitted polynomial coefficients:") for i, coef in enumerate(beta): if cov_beta is not None: print(f" c{i} = {coef:g} +- {beta_std[i]:g}") else: print(f" c{i} = {coef:g} (Error estimation skipped)") fit = tkFit() fit.to_excel(cfg.output_parameter_path, ["coeff", "std"], [beta, beta_std]) sigma_meas = compute_measurement_error(y_fit) Tcalmin = pfloat(cfg.Tcalmin, defval=min(T)) Tcalmax = pfloat(cfg.Tcalmax, defval=max(T)) xcal = np.linspace(1000 / Tcalmax, 1000 / Tcalmin, cfg.ncal) bands = compute_bands(xcal, beta, cov_beta, sigma2_resid, sigma_meas) # Excel 保存 (データ構造維持) xlabel, ylabel = label_x, label_y X_for_cal = build_design_matrix(T1000, norder) ycal_all = X_for_cal @ beta fit.to_excel(cfg.output_fitting_path, [xlabel, ylabel, "", "1000/T (K^-1)", "log10(P)", "log10(P)(cal)", "", "1000/T (K^-1)", "log10(P)(cal)", "sigma(param)", 'sigma(param&resid)', 'sigma(param&noise)'], [xX, yY, [], T1000, log10P, ycal_all, [], xcal, bands['y_mean'], bands['sigma_param'], bands['sigma_pred'], bands['sigma_combined']], template = cfg.xlsm_template) # Plot diff_plot = np.gradient(bands['y_mean'], xcal) Ea_plot = [-d * kB / e * 1000.0 * log(10.0) for d in diff_plot] fig, axes = plt.subplots(2, 3, figsize=cfg.figsize) axes = axes.flatten() axes[0].plot(xX, yY, 'ko', markersize=5) axes[1].plot(T, P, 'ko', markersize=5) # Arrhenius Plot axes[2].plot(x_fit, y_fit, 'ko', label='data', markersize=1.5) axes[2].plot(xcal, bands['y_mean'], 'r-', label='fit', linewidth=0.5) # 修正箇所: cov_beta がある場合のみ信頼区間を描画[cite: 1] if cov_beta is not None: if cfg.plot_sigma_param: axes[2].fill_between(xcal, bands['y_mean'] - bands['sigma_param'], bands['y_mean'] + bands['sigma_param'], color='#0000cc', alpha=0.5, label='±σ(param)') if cfg.plot_sigma_pred: axes[2].fill_between(xcal, bands['y_mean'] - bands['sigma_pred'], bands['y_mean'] + bands['sigma_pred'], color='#ddddFF', alpha=0.5, label='±σ(param&resid)') axes[2].legend(fontsize=cfg.fontsize_legend) axes[3].plot(xcal, Ea_plot, 'k-', linewidth=0.5, label='Ea (eV)') # パラメータプロット if axes[4]: idxs = np.arange(len(beta)) # 修正箇所: cov_beta がある場合のみエラーバーを描画[cite: 1] if cov_beta is not None: axes[4].errorbar(idxs, beta, yerr=beta_std, fmt='o', capsize=3, label='coeff ±1σ') else: axes[4].plot(idxs, beta, 'o', label='coeff') axes[4].legend(fontsize=cfg.fontsize_legend) plt.tight_layout() plt.show() app.terminate("", usage=None, pause=cfg.pause) def main(): app, cfg, parser = initialize() execute(app) if __name__ == "__main__": main()