import argparse import sys from pathlib import Path from shutil import copy2 import matplotlib.pyplot as plt import numpy as np from openpyxl import Workbook, load_workbook from scipy.constants import Boltzmann, elementary_charge, h, m_e from tkplot.tkplotevent import tkPlotEvent from tkutils import read_data_table, redirect_output from tklsq import ( delta_method_variance, polynomial_design_matrix, polynomial_lsq, predict_polynomial, regression_scores, ) from tklsq.tkfitdiag import diagnose_covariance K_B_EV = Boltzmann / elementary_charge def float_or_auto(value): """'*' または空文字を None、それ以外を float に変換する。""" if value is None: return None if isinstance(value, str) and value.strip() in ('', '*'): return None return float(value) def initialize(argv=None): """コマンドライン引数を読み込み、argparse.Namespace を返す。""" parser = argparse.ArgumentParser(description="Arrhenius plot と物理モデルフィット") parser.add_argument('--infile', type=str, default="Hall-T.xlsx", help='Input file') parser.add_argument( '--model', '--mode', dest='model', type=str, default='simple', help="simple, Nc, Seto, Percolation, Hopping（--mode も同じ）") parser.add_argument('--Tlabel', type=str, default='T(K)', help='T関連データ列ラベル') parser.add_argument('--Plabel', type=str, default='P', help='P関連データ列ラベル') parser.add_argument('--header_row', type=int, default=1, help='表のヘッダー行（1始まり）') parser.add_argument('--first_column', type=int, default=1, help='表の先頭列（1始まり）') parser.add_argument('--sheet', default=None, help='Excelシート名。省略時はアクティブシート') 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=float_or_auto, default=None, help="計算に用いるTの下限（'*' または省略で自動）") parser.add_argument('--Tcalmax', type=float_or_auto, default=None, help="計算に用いるTの上限（'*' または省略で自動）") 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('--diagnostics', type=int, default=0, help='Print coefficient covariance diagnostics') 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) cfg = parser.parse_args(argv) input_path = Path(cfg.infile) cfg.logfile = input_path.with_name(f"{input_path.stem}.log") cfg.output_fitting_path = input_path.with_name(f"{input_path.stem}-fit.xlsm") cfg.output_parameter_path = input_path.with_name( f"{input_path.stem}-parameters.xlsx") return cfg, parser def _excel_value(value): """NumPy値をopenpyxlで保存可能なPython値に変換する。""" if isinstance(value, np.generic): value = value.item() if isinstance(value, float) and not np.isfinite(value): return None return value def resolve_template_path(template, input_path=None): """Excelテンプレートを、代表的な相対パス候補から探索する。""" if template is None: return None template_path = Path(template).expanduser() if template_path.is_absolute(): return template_path if template_path.exists() else None candidates = [Path.cwd() / template_path] if input_path is not None: candidates.append(Path(input_path).resolve().parent / template_path) candidates.append(Path(__file__).resolve().parent / template_path) for candidate in candidates: if candidate.exists(): return candidate return None def write_columns_to_excel(output_path, headers, columns, template=None, input_path=None): """ 列配列をExcelへ保存する。 template が見つかればコピーしてVBAを保持する。見つからない場合は 例外で解析全体を止めず、通常の .xlsx ファイルへフォールバックする。 戻り値は実際に保存したファイルの Path。 """ output_path = Path(output_path) output_path.parent.mkdir(parents=True, exist_ok=True) if len(headers) != len(columns): raise ValueError('headers and columns must have the same length') template_path = resolve_template_path(template, input_path=input_path) if template is not None and template_path is None: fallback_path = output_path.with_suffix('.xlsx') print( f"WARNING: Excel template [{template}] was not found. " f"Saving data without the template to [{fallback_path}]." ) output_path = fallback_path wb = Workbook() elif template_path is not None: copy2(template_path, output_path) wb = load_workbook( output_path, keep_vba=output_path.suffix.lower() == '.xlsm') else: wb = Workbook() ws = wb.active for col_idx, header in enumerate(headers, start=1): ws.cell(row=1, column=col_idx, value=header) values = [] if columns[col_idx - 1] is None else list(columns[col_idx - 1]) for row_idx, value in enumerate(values, start=2): ws.cell(row=row_idx, column=col_idx, value=_excel_value(value)) wb.save(output_path) return output_path def print_scores(heading, y1, y2, p=0): """tklsq.regression_scores() の結果を従来形式で表示する。""" scores = regression_scores(y1, y2, p=p) print(heading) print(f" N = {int(scores['N'])}") print(f" MAE = {scores['MAE']:.8g}") print(f" RMSE = {scores['RMSE']:.8g}") print(f" max|err| = {scores['max_abs_error']:.8g}") print(f" R2 = {scores['R2']:.8g}") def normalize_model(model): """モデル名を正規化する（大文字・小文字は区別しない）。""" table = { 'simple': 'simple', 'nc': 'Nc', 'seto': 'Seto', 'percolation': 'Percolation', 'hopping': 'Hopping', # 旧デフォルト値との互換性 'simple arrhenius #log(p)=a-(eea/kb)*(1/t)': 'simple', } key = str(model).strip().lower() if key not in table: raise ValueError( f"Unknown model [{model}]. Choose simple, Nc, Seto, " "Percolation, or Hopping.") return table[key] def transformed_ordinate(T, P, model): """ 各モデルを直線化した従属変数を返す。 simple: log10(P) Nc: log10(P / T^(3/2)) Seto: log10(P * T^(1/2)) Percolation: log10(P) Hopping: log10(P * T) """ T = np.asarray(T, dtype=float) P = np.asarray(P, dtype=float) factors = { 'simple': np.ones_like(T), 'Nc': T**(-1.5), 'Seto': np.sqrt(T), 'Percolation': np.ones_like(T), 'Hopping': T, } return np.log10(P * factors[model]), factors[model] def physical_parameters(beta, cov_beta, model): """ 回帰係数を物理パラメータへ変換し、delta method で1σ誤差を求める。 Ncモデルでは P を非縮退キャリア濃度 [cm^-3] と仮定する。 Percolationモデルでは障壁高さがGaussian分布に従う近似を用いる。 """ beta = np.asarray(beta, dtype=float) kB_eV = K_B_EV slope_to_eV = 1000.0 * kB_eV * np.log(10.0) rows = [] def add(name, value, gradient, unit): std = np.nan if cov_beta is not None: gradient = np.asarray(gradient, dtype=float) variance = float(gradient @ cov_beta @ gradient) std = np.sqrt(max(variance, 0.0)) rows.append((name, float(value), float(std), unit)) if model in ('simple', 'Nc', 'Seto', 'Hopping'): grad_E = np.zeros_like(beta) grad_E[1] = -slope_to_eV energy_name = 'barrier_height' if model == 'Seto' else 'activation_energy' add(energy_name, -slope_to_eV * beta[1], grad_E, 'eV') if model == 'simple': value = 10.0**beta[0] grad = np.zeros_like(beta) grad[0] = np.log(10.0) * value add('P0', value, grad, 'same unit as P') elif model == 'Nc': # Nc = 2(2*pi*m_dos*kB*T/h^2)^(3/2) [m^-3] # P/T^(3/2) の前置因子は cm^-3 K^-3/2。 A_cm = 10.0**beta[0] mass_kg = (h**2 / (2.0 * np.pi * kB)) * ( A_cm * 1.0e6 / 2.0)**(2.0 / 3.0) mass_ratio = mass_kg / m_e grad = np.zeros_like(beta) grad[0] = (2.0 / 3.0) * np.log(10.0) * mass_ratio add('DOS_effective_mass', mass_ratio, grad, 'm_e') grad_A = np.zeros_like(beta) grad_A[0] = np.log(10.0) * A_cm add('Nc_prefactor', A_cm, grad_A, 'cm^-3 K^-3/2') elif model == 'Seto': value = 10.0**beta[0] grad = np.zeros_like(beta) grad[0] = np.log(10.0) * value add('Seto_prefactor', value, grad, 'P K^1/2') elif model == 'Percolation': # log10(P)=c0-Ebar/(1000*kB*ln10)x # +sigma_E^2/(2*(1000*kB)^2*ln10)x^2 grad_mean = np.zeros_like(beta) grad_mean[1] = -slope_to_eV add('mean_barrier_height', -slope_to_eV * beta[1], grad_mean, 'eV') if beta[2] > 0.0: sigma_E = 1000.0 * kB_eV * np.sqrt( 2.0 * np.log(10.0) * beta[2]) grad_sigma = np.zeros_like(beta) grad_sigma[2] = sigma_E / (2.0 * beta[2]) add('barrier_height_std', sigma_E, grad_sigma, 'eV') else: rows.append(('barrier_height_std', np.nan, np.nan, 'eV (undefined: c2 <= 0)')) elif model == 'Hopping': sigma0 = 10.0**beta[0] grad = np.zeros_like(beta) grad[0] = np.log(10.0) * sigma0 add('sigma0', sigma0, grad, 'P K') return rows def polynomial_derivative(x, beta): """多項式を解析微分する。""" x = np.asarray(x, dtype=float) derivative = np.zeros_like(x) for j in range(1, len(beta)): derivative += j * beta[j] * x**(j - 1) return derivative def activation_energy_curve(x, beta, cov_beta): """局所勾配から見かけの活性化エネルギーと1σ誤差を返す。""" x = np.asarray(x, dtype=float) factor = -1000.0 * K_B_EV * np.log(10.0) Ea = factor * polynomial_derivative(x, beta) if cov_beta is None: return Ea, None J = np.zeros((len(x), len(beta))) for j in range(1, len(beta)): J[:, j] = factor * j * x**(j - 1) var = delta_method_variance(J, cov_beta) return Ea, np.sqrt(np.maximum(var, 0.0)) def execute(cfg): try: cfg.model = normalize_model(cfg.model) except ValueError as exc: raise ValueError(str(exc)) from exc # 以下元コード同様にデータ読み込みと変換 print("#======================================================") print("# Analyze activation energy etc by Arrhenius plot") print("#======================================================") print(f"infile : {cfg.infile}") print(f"model : {cfg.model}") print(f"Tlabel : {cfg.Tlabel}, Plabel: {cfg.Plabel}") print(f"Ttype : {cfg.Ttype}, Ptype: {cfg.Ptype}") print(f"Fitting x range: {cfg.xmin} - {cfg.xmax}") print(f"Fitting T range: {cfg.Tmin} - {cfg.Tmax}") # データ読み込み print(f"Read [{cfg.infile}]") table = read_data_table( cfg.infile, header_row=cfg.header_row, first_column=cfg.first_column, force_numeric=True, sheet=cfg.sheet, ) labels, datalist = table.labels, table.columns label_x, xX = table.find_column(cfg.Tlabel, ignore_case=True) label_y, yY = table.find_column(cfg.Plabel, ignore_case=True) if xX is None or yY is None: raise ValueError("指定ラベルのデータが見つかりません") # T の変換 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: raise ValueError(f"Invalid Ttype [{cfg.Ttype}]") # P の変換 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: raise ValueError(f"Invalid Ptype [{cfg.Ptype}]") T = np.asarray(T, dtype=float) P = np.asarray(P, dtype=float) if len(T) != len(P): raise ValueError("T and P have different lengths") valid = np.isfinite(T) & np.isfinite(P) & (T > 0.0) & (P > 0.0) if not np.all(valid): print(f"WARNING: {np.count_nonzero(~valid)} rows with T<=0, P<=0, " "NaN, or infinity are ignored.") xX_valid = np.asarray(xX)[valid] yY_valid = np.asarray(yY)[valid] T = T[valid] P = P[valid] if len(T) == 0: raise ValueError("No valid data (T>0 and P>0)") print("T(K)=", T) print("P =", P) # フィッティング対象抽出 T1000 = 1000.0 / T log10P = np.log10(P) y_model, model_factor = transformed_ordinate(T, P, cfg.model) x_fit = [] y_fit = [] T_fit = [] P_fit_input_raw = [] for xi_orig, temperature, xi, yi, pi in zip( xX_valid, T, T1000, y_model, P): 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) T_fit.append(temperature) P_fit_input_raw.append(pi) # 多項式次数決定 if cfg.model == 'Percolation': norder = 2 else: norder = 1 if len(x_fit) < norder + 1: raise ValueError( f"{len(x_fit)} fitting points for {norder + 1} parameters") print(f"\nLeast-squares fitting with {norder}-th order polynomial of 1000/T") print(f"Physical model: {cfg.model}") print("Data to be fitted:") print(f" {'1000/T':12} {'model ordinate':16}") for xi, yi in zip(x_fit, y_fit): print(f" {xi:12.4g} {yi:12.4g}") X = polynomial_design_matrix(x_fit, order=norder) lsq_result = polynomial_lsq(x_fit, y_fit, order=norder) beta = lsq_result.beta beta_std = lsq_result.beta_std cov_beta = lsq_result.cov_beta sigma2_resid = lsq_result.sigma2_resid error_estimation_enabled = lsq_result.error_estimation_enabled print( f"Fit diagnostics: N={lsq_result.N}, p={lsq_result.p}, " f"dof={lsq_result.dof}, rank={lsq_result.rank}, " f"RSS={lsq_result.RSS:12.4g}, " f"cond={lsq_result.condition_number:12.4g}" ) if lsq_result.warning: print(lsq_result.warning) if cfg.diagnostics and cov_beta is not None: coefficient_names = [f"c{i}" for i in range(len(beta))] diagnostics = diagnose_covariance( coefficient_names, beta, cov_beta, jacobian=X, ) print("\nCoefficient correlation matrix:") print(diagnostics.corr) if diagnostics.cond_jtj is not None: print(f"cond(J^T J) = {diagnostics.cond_jtj:.8g}") if diagnostics.warning: print(diagnostics.warning) print("Fitted polynomial coefficients:") for i, coef in enumerate(beta): if beta_std is None: print(f" c{i} = {coef:g}") else: print(f" c{i} = {coef:g} +- {beta_std[i]:g}") if sigma2_resid is None: print(" Residual variance sigma2_resid = undefined") else: print(f" Residual variance sigma2_resid = {sigma2_resid:12.4g}") parameter_rows = physical_parameters(beta, cov_beta, cfg.model) print("\nPhysical parameters (standard uncertainty = 1 sigma):") for name, value, std, unit in parameter_rows: if np.isfinite(std): print(f" {name} = {value:.8g} +- {std:.3g} {unit}") else: print(f" {name} = {value:.8g} {unit}") print(f"Save parameters to {cfg.output_parameter_path}") parameter_names = [row[0] for row in parameter_rows] parameter_values = [row[1] for row in parameter_rows] parameter_stds = [row[2] for row in parameter_rows] parameter_units = [row[3] for row in parameter_rows] write_columns_to_excel( cfg.output_parameter_path, ["parameter", "value", "std(1sigma)", "unit"], [parameter_names, parameter_values, parameter_stds, parameter_units]) # 独立な測定誤差が入力されていないため、全データの標準偏差を # 測定誤差として扱わない。予測帯には残差分散を使用する。 sigma_meas = None Tcalmin = min(T) if cfg.Tcalmin is None else cfg.Tcalmin Tcalmax = max(T) if cfg.Tcalmax is None else cfg.Tcalmax if cfg.ncal < 2: raise ValueError("ncal must be >= 2") if Tcalmin <= 0.0 or Tcalmax <= Tcalmin: raise ValueError("require 0 < Tcalmin < Tcalmax") Tstep = (Tcalmax - Tcalmin) / (cfg.ncal - 1) T_plot = [Tcalmin + i * Tstep for i in range(cfg.ncal)] x_plot = [1000.0 / v for v in T_plot] xcal = np.linspace(1000 / Tcalmax, 1000 / Tcalmin, cfg.ncal) prediction = predict_polynomial(xcal, lsq_result, order=norder) ycal = prediction.y_mean # model ordinate prediction P_plot = [10**val for val in ycal] _, factor_fit = transformed_ordinate( np.asarray(T_fit), np.asarray(P_fit_input_raw), cfg.model) P_fit = 10.0**(X @ beta) / factor_fit print_scores( heading="\nScores between P(input) and P(fit)", y1=P_fit_input_raw, y2=P_fit, p=len(beta)) print_scores( heading="\nScores in the model-transformed logarithmic ordinate", y1=y_fit, y2=lsq_result.y_fit, p=len(beta)) bands = { 'y_mean': prediction.y_mean, 'sigma_param': prediction.sigma_param, 'sigma_pred': prediction.sigma_pred, # 独立な測定誤差が無いため、従来通りparamのみと同じ。 'sigma_combined': prediction.sigma_param, } print(f"Save results to [{cfg.output_fitting_path}]") xlabel = label_x ylabel = label_y X_for_cal = polynomial_design_matrix(T1000, order=norder) ymodel_cal = X_for_cal @ beta log10P_cal = ymodel_cal - np.log10(model_factor) Ea_curve, Ea_curve_std = activation_energy_curve(xcal, beta, cov_beta) if error_estimation_enabled: write_columns_to_excel(cfg.output_fitting_path, [xlabel, ylabel, "T(K)", "P", "1000/T (K^-1)", "log10(P)", "model ordinate", "log10(P)(cal)", "", "1000/T (K^-1)", "model ordinate(cal)", "sigma(param)", 'sigma(param&resid)', "Ea/apparent barrier (eV)", "Ea std(1sigma)"], [xX_valid, yY_valid, T, P, T1000, log10P, y_model, log10P_cal, [], xcal, bands['y_mean'], bands['sigma_param'], bands['sigma_pred'], Ea_curve, Ea_curve_std], template=cfg.xlsm_template, input_path=cfg.infile) else: print("WARNING: Uncertainty columns are omitted from Excel output.") write_columns_to_excel(cfg.output_fitting_path, [xlabel, ylabel, "T(K)", "P", "1000/T (K^-1)", "log10(P)", "model ordinate", "log10(P)(cal)", "", "1000/T (K^-1)", "model ordinate(cal)", "Ea/apparent barrier (eV)"], [xX_valid, yY_valid, T, P, T1000, log10P, y_model, log10P_cal, [], xcal, bands['y_mean'], Ea_curve], template=cfg.xlsm_template, input_path=cfg.infile) """ write_columns_to_excel(cfg.output_fitting_path, [xlabel, ylabel, 'T(K)', 'P', '1000/T (K^-1)', 'log10(P)', "log10(P)(cal)", "", "T (K)", "P(cal)", "sigma(param)", 'sigma(param&resid)', 'sigma(param&noise)'], [xX, yY, T, P, T1000, log10P, ycal, [], xcal, bands['y_mean'], bands['sigma_param'], bands['sigma_pred'], bands['sigma_combined']], template = cfg.xlsm_template) """ """ print(f"Save to [{cfg.outfile}]") write_columns_to_excel(cfg.outfile, [label_x, label_y, 'T(K)', 'P', '1000/T (K^-1)', 'log10(P)', 'log10(P)(cal)', '', 'T (K)', '1000/T (K^-1)', 'P(cal)', 'log10(P)(cal)'], [xX, yY, T, P, T1000, log10P, list(X @ beta), [], T_plot, x_plot, P_plot, list(ycal)]) """ # Ea plot: フィット多項式を解析微分 Ea_plot, Ea_plot_std = activation_energy_curve( xcal, beta, cov_beta) fig, axes = plt.subplots(2, 3, figsize=cfg.figsize) plot_event = tkPlotEvent(plt) axes = axes.flatten() for ax in axes: ax.tick_params(labelsize=cfg.fontsize) # (0) X-Y axes[0].plot(xX, yY, linestyle='', marker='o', markerfacecolor='black', markersize=5.0) axes[0].set_xlabel(label_x, fontsize=cfg.fontsize) axes[0].set_ylabel(label_y, fontsize=cfg.fontsize) # (1) T-P axes[1].plot(T, P, linestyle='', marker='o', markerfacecolor='black', markersize=5.0) # axes[1].plot(T_plot, P_plot, linestyle='-', color='red', linewidth=0.5) axes[1].set_xlabel('$T$ (K)', fontsize=cfg.fontsize) axes[1].set_ylabel('$P$', fontsize=cfg.fontsize) # (2) Arrhenius plot axes[2].plot(x_fit, y_fit, 'o', label='data', color = 'black', markersize=1.5) axes[2].plot(xcal, bands['y_mean'], label='fit', linewidth = 0.5, color='red') if cfg.plot_sigma_param and bands['sigma_param'] is not None: 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 and bands['sigma_pred'] is not None: 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)') if cfg.plot_sigma_combined and bands['sigma_combined'] is not None: axes[2].fill_between(xcal, bands['y_mean'] - bands['sigma_combined'], bands['y_mean'] + bands['sigma_combined'], color='purple', alpha=0.5, label='±σ(param&noise)') axes[2].set_xlabel('$1000/T$ (K$^{-1}$)', fontsize=cfg.fontsize) ordinate_labels = { 'simple': r'$\log_{10}(P)$', 'Nc': r'$\log_{10}(P/T^{3/2})$', 'Seto': r'$\log_{10}(P\sqrt{T})$', 'Percolation': r'$\log_{10}(P)$', 'Hopping': r'$\log_{10}(PT)$', } axes[2].set_ylabel(ordinate_labels[cfg.model], fontsize=cfg.fontsize) axes[2].legend(fontsize=cfg.fontsize_legend) axes[3].plot(xcal, Ea_plot, label='Ea (eV)', linestyle='-', linewidth=0.5, color='black') if Ea_plot_std is not None: axes[3].fill_between( xcal, Ea_plot - Ea_plot_std, Ea_plot + Ea_plot_std, color='#dddddd', alpha=0.7, label='±1σ') axes[3].set_xlabel('$1000/T$ (K$^{-1}$)', fontsize=cfg.fontsize) axes[3].set_ylabel('$E_a$ (eV)', fontsize=cfg.fontsize) axes[3].legend(fontsize=cfg.fontsize_legend) minEa = min(Ea_plot) maxEa = max(Ea_plot) if abs(maxEa - minEa) < 1.0e-6: minEa -= 1.0e-3 maxEa += 1.0e-3 axes[3].set_ylim([minEa, maxEa]) if axes[4]: idxs = np.arange(len(beta)) if beta_std is None: axes[4].plot(idxs, beta, 'o', label='coeff') axes[4].text(0.05, 0.95, 'parameter error disabled', transform=axes[4].transAxes, va='top', fontsize=cfg.fontsize_legend) else: axes[4].errorbar(idxs, beta, yerr=beta_std, fmt='o', capsize=3, label='coeff ±1σ') axes[4].set_xlabel('$i$', fontsize=cfg.fontsize) axes[4].set_ylabel('$c_i$', fontsize=cfg.fontsize) axes[4].tick_params(labelsize=cfg.fontsize) axes[4].legend(fontsize=cfg.fontsize_legend) # plot_event 登録 # data=None の場合、新パッケージは axis.lines を自動取得する。 # ただし Arrhenius/Ea 軸にはフィット線もあるため、クリック対象を # 実測点または主曲線に明示的に限定する。 all_data = datalist plot_event.add_data({"label": "X-Y plot", "plot_type": "2D", "axis": axes[0], "data": [axes[0].lines[0]], "xlist": all_data, "xlabels": labels}) plot_event.add_data({"label": "T-P plot", "plot_type": "2D", "axis": axes[1], "data": [axes[1].lines[0]], "xlist": all_data, "xlabels": labels}) plot_event.add_data({"label": "Arrhenius data", "plot_type": "2D", "axis": axes[2], "data": [axes[2].lines[0]]}) plot_event.add_data({"label": "Ea curve", "plot_type": "2D", "axis": axes[3], "data": [axes[3].lines[0]]}) plot_event.register_event(fig, event="button_press_event") plt.tight_layout() plt.pause(0.001) def main(argv=None): cfg = None try: cfg, _ = initialize(argv) print(f"Open logfile [{cfg.logfile}]") with redirect_output(cfg.logfile): execute(cfg) except Exception as exc: print(f"Error: {exc}", file=sys.stderr) exit_code = 1 else: exit_code = 0 if cfg is not None and cfg.pause: try: input("Press Enter to exit...") except EOFError: pass return exit_code if __name__ == "__main__": raise SystemExit(main())