import argparse import sys from contextlib import redirect_stderr, redirect_stdout 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 tklib.tkvariousdata import tkVariousData from tkplot.tkplotevent import tkPlotEvent 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('--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('--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 class Tee: """複数ストリームへ同じ文字列を書き込む簡単な tee。""" def __init__(self, *streams): self.streams = streams def write(self, text): for stream in self.streams: stream.write(text) stream.flush() return len(text) def flush(self): for stream in self.streams: stream.flush() 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): """2系列間の基本的な適合度指標を表示する。""" y1 = np.asarray(y1, dtype=float) y2 = np.asarray(y2, dtype=float) if y1.shape != y2.shape: raise ValueError(f"score arrays have different shapes: {y1.shape} and {y2.shape}") valid = np.isfinite(y1) & np.isfinite(y2) y1 = y1[valid] y2 = y2[valid] if len(y1) == 0: print(heading) print(' No finite paired data') return resid = y1 - y2 mae = float(np.mean(np.abs(resid))) rmse = float(np.sqrt(np.mean(resid**2))) max_abs = float(np.max(np.abs(resid))) ss_res = float(np.sum(resid**2)) ss_tot = float(np.sum((y1 - np.mean(y1))**2)) r2 = np.nan if ss_tot == 0.0 else 1.0 - ss_res / ss_tot print(heading) print(f" N = {len(y1)}") print(f" MAE = {mae:.8g}") print(f" RMSE = {rmse:.8g}") print(f" max|err| = {max_abs:.8g}") print(f" R2 = {r2:.8g}") def build_design_matrix(x, order): """ Construct the design matrix for polynomial regression. Args: x (array-like): Input data points. order (int): Polynomial order. Returns: np.ndarray: (N, order+1) design matrix. """ x = np.asarray(x) N = len(x) X = np.ones((N, order + 1)) for j in range(1, order + 1): # Using scalar pow for clarity X[:, j] = [pow(val, j) for val in x] return X def mlsq_error(X, y, rcond=None): """ Perform least squares fitting with optional error estimation. If the residual degrees of freedom N - p is zero or negative, the regression itself is still performed, but residual variance, parameter standard errors, covariance matrix, and confidence/prediction bands are disabled because the unbiased residual variance cannot be estimated. Args: X (np.ndarray): Design matrix (N, p). y (array-like): Observed values (N,). rcond (float or None): Cutoff for np.linalg.lstsq. Returns: beta (np.ndarray): Estimated coefficients (p,). beta_std (np.ndarray or None): Standard error of parameters. cov_beta (np.ndarray or None): Covariance matrix of parameters. sigma2_resid (float or None): Residual variance estimate. fit_info (dict): Diagnostic information. """ X = np.asarray(X, dtype=float) y = np.asarray(y, dtype=float) N, p = X.shape dof = N - p beta, residuals_lstsq, rank, svals = np.linalg.lstsq(X, y, rcond=rcond) residuals = y - X @ beta RSS = float(residuals.T @ residuals) fit_info = { 'N': N, 'p': p, 'dof': dof, 'rank': rank, 'RSS': RSS, 'singular_values': svals, 'error_estimation_enabled': True, 'warning': '', } if rank < p: fit_info['error_estimation_enabled'] = False fit_info['warning'] = ( f"WARNING: Design matrix is rank deficient (rank={rank}, p={p}). " "Least-squares coefficients were calculated by np.linalg.lstsq, " "but parameter errors/covariance and uncertainty bands are disabled." ) return beta, None, None, None, fit_info if dof <= 0: fit_info['error_estimation_enabled'] = False fit_info['warning'] = ( f"WARNING: Residual degrees of freedom N - p = {N} - {p} = {dof}. " "The residual variance is undefined, so parameter errors, " "prediction errors, confidence intervals, and uncertainty-band plots " "are disabled. Least-squares fitting itself continues." ) return beta, None, None, None, fit_info XtX_inv = np.linalg.inv(X.T @ X) sigma2_resid = RSS / dof cov_beta = sigma2_resid * XtX_inv beta_std = np.sqrt(np.diag(cov_beta)) return beta, beta_std, cov_beta, sigma2_resid, fit_info def compute_param_uncertainty(X, cov_beta): """ Compute parameter-based prediction variance at each X row. Args: X (np.ndarray): Design matrix (N, p). cov_beta (np.ndarray): Covariance matrix (p, p). Returns: y_var (np.ndarray): Variance of predictions (N,). """ # Var[y_mean] = diag(X @ cov_beta @ X.T) return np.sum((X @ cov_beta) * X, axis=1) def compute_measurement_error(y): """ Estimate measurement error from input data variance. Args: y (array-like): Observed values. Returns: float: Estimated measurement std (unbiased). """ 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=None, sigma2_resid=None, sigma_meas=None): """ Compute mean predictions and, when available, uncertainty bands on xcal. If cov_beta or sigma2_resid is None, only y_mean is returned and all uncertainty arrays are set to None. This is used for exactly determined or underdetermined fits where residual variance cannot be estimated. """ Xcal = build_design_matrix(xcal, len(beta) - 1) y_mean = Xcal @ beta if cov_beta is None or sigma2_resid is None: return { 'y_mean': y_mean, 'sigma_param': None, 'sigma_pred': None, 'sigma_combined': None } if sigma_meas is None: sigma_meas = 0.0 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) return { 'y_mean': y_mean, 'sigma_param': sigma_param, 'sigma_pred': sigma_pred, 'sigma_combined': sigma_combined } 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 = np.sum((J @ cov_beta) * J, axis=1) 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}]") datafile = tkVariousData(cfg.infile) labels, datalist = datafile.Read_minimum_matrix(close_fp=True, usage=None) 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: 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 = build_design_matrix(x_fit, norder) beta, beta_std, cov_beta, sigma2_resid, fit_info = mlsq_error(X, y_fit) error_estimation_enabled = fit_info['error_estimation_enabled'] print(f"Fit diagnostics: N={fit_info['N']}, p={fit_info['p']}, " f"dof={fit_info['dof']}, rank={fit_info['rank']}, RSS={fit_info['RSS']:12.4g}") if fit_info['warning']: print(fit_info['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) ycal = build_design_matrix(xcal, norder) @ beta # log10P 予測 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) print_scores( heading="\nScores in the model-transformed logarithmic ordinate", y1=y_fit, y2=list(X @ beta)) bands = compute_bands(xcal, beta, cov_beta, sigma2_resid, sigma_meas) print(f"Save results to [{cfg.output_fitting_path}]") xlabel = label_x ylabel = label_y X_for_cal = build_design_matrix(T1000, 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 Path(cfg.logfile).open('w', encoding='utf-8') as log_fp: tee_out = Tee(sys.stdout, log_fp) tee_err = Tee(sys.stderr, log_fp) with redirect_stdout(tee_out), redirect_stderr(tee_err): 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())