import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.optimize import minimize from numpy import sqrt, log, exp import argparse import os import traceback from typing import Dict, Any, List, Tuple plt.rcParams["font.family"] = "MS Gothic" # --- デフォルト値 --- EF0_def = 0.0 BG0_def = 0.0 w0_def = 0.20 I00_def = 200000.0 DOSa0_def = 0.0 # フィッティングパラメータの範囲、初期値、最適化ID、変換タイプ、EPSを定義 fit_params = { "EF": { "optid": 1, "convert": "", "eps": 1.0e-5, "x0": EF0_def, "min": -0.5, "max": 0.5, "kpenalty": 1e8 }, "BG": { "optid": 1, "convert": "", "eps": 1.0e-3, "x0": BG0_def, "min": None, "max": None, "kpenalty": 1e8 }, "w": { "optid": 1, "convert": "log", "eps": 1.0e-5, "x0": w0_def, "min": 0.01, "max": 0.5, "kpenalty": 1e8 }, "I0": { "optid": 1, "convert": "log", "eps": 1.0e+0, "x0": I00_def, "min": 0.0, "max": None, "kpenalty": 1e8 }, "DOSa": { "optid": 0, "convert": "", "eps": 1.0e-5, "x0": DOSa0_def, "min": None, "max": None, "kpenalty": 1e8 }, } PRINT_INTERVAL_def = 5 NMAXITER_def = 1000 # パラメータ名のリストをグローバルに定義 ALL_PARAM_NAMES = list(fit_params.keys()) OPT_PARAM_NAMES = [name for name, info in fit_params.items() if info["optid"] == 1] FIXED_PARAM_NAMES = [name for name, info in fit_params.items() if info["optid"] == 0] def terminate(): input("\nPress ENTER to terminate>>\n") exit() # ------------------------------- # 変換関数 # ------------------------------- def convert_params(params_dict): """ 全パラメータ辞書から、最適化対象のみを抽出し、logスケールに変換して配列として返す。 """ opt_x0 = [] for name in OPT_PARAM_NAMES: value = params_dict[name] info = fit_params[name] if info["convert"] == "log": if value <= 0: raise ValueError(f"Parameter '{name}' must be positive for log transformation.") opt_x0.append(log(value)) else: opt_x0.append(value) return np.array(opt_x0) def recover_params(opt_params_array, all_initial_params_dict): """ 最適化パラメータ配列と、固定パラメータの初期値から、全パラメータの物理量リストを再構成。 戻り値は [EF, BG, w, I0, DOSa] の順序のリスト。 """ all_params_dict = all_initial_params_dict.copy() opt_params_idx = 0 for name in OPT_PARAM_NAMES: value_raw = opt_params_array[opt_params_idx] info = fit_params[name] if info["convert"] == "log": value_recovered = exp(value_raw) else: value_recovered = value_raw all_params_dict[name] = value_recovered opt_params_idx += 1 return [all_params_dict[name] for name in ALL_PARAM_NAMES] def kG2fwhm(kG): """kGからFWHM wに変換 (配列対応)""" return np.where(kG <= 0, 1e10, 2.0 * sqrt(log(2.0) / kG)) def fwhm2kG(w): """FWHM wからkGに変換 (配列対応)""" return np.where(w <= 0, 1e20, log(2.0) / (0.5 * w)**2) # ------------------------------- # 引数初期化 # ------------------------------- def initialize(): parser = argparse.ArgumentParser(description="光電子分光スペクトル解析") parser.add_argument("filename", help="入力ファイル（csv, xlsx, txt）") parser.add_argument("--mode", choices=["fit", "sim", "init"], default="fit", help="モード: fit, sim, or init (初期値推定)") parser.add_argument("--method", default="Nelder-Mead", help="最適化アルゴリズム (デフォルト Nelder-Mead") parser.add_argument("--T", type=float, default=300.0, help="温度 [K]") parser.add_argument("--EF", type=float, default=None, help=f"EF [eV] (デフォルト: {fit_params['EF']['x0']})") parser.add_argument("--BG", type=float, default=None, help=f"BG (デフォルト: {fit_params['BG']['x0']})") parser.add_argument("--w", type=float, default=None, help=f"装置関数fwhm (デフォルト: {fit_params['w']['x0']})") parser.add_argument("--I0", type=float, default=None, help=f"シグナル強度I0 (デフォルト: {fit_params['I0']['x0']})") parser.add_argument("--DOSa", type=float, default=None, help=f"DOSの線形傾き (デフォルト: {fit_params['DOSa']['x0']})") parser.add_argument("--Efit_min", type=float, default=None, help="フィット範囲下限 [eV]") parser.add_argument("--Efit_max", type=float, default=None, help="フィット範囲上限 [eV]") parser.add_argument("--Escale", type=float, default=-1.0, help="強度スケーリング係数") parser.add_argument("--nmaxiter", type=int, default=NMAXITER_def, help=f"最適化の最大反復回数 (デフォルト: {NMAXITER_def})") parser.add_argument("--print_interval", type=int, default=PRINT_INTERVAL_def, help=f"ステップ出力間隔 (デフォルト: {PRINT_INTERVAL_def})") return parser.parse_args() # ------------------------------- # CSVパラメータファイルの読み書き # ------------------------------- def get_prm_csv_path(filename): """パラメータCSVファイルのパスを生成""" base, ext = os.path.splitext(filename) return base + "_parameters.csv" def load_params(prm_csv_path): """ CSVファイルからパラメータ設定を読み込む。 """ if os.path.exists(prm_csv_path): try: df = pd.read_csv(prm_csv_path, index_col='Parameter') for name, info in fit_params.items(): if name in df.index: for col in ['optid', 'convert', 'eps', 'min', 'max', 'kpenalty']: df.loc[name, col] = info[col] if info[col] is not None else np.nan print(f"Status: Parameters loaded from '{os.path.basename(prm_csv_path)}'.") return df except Exception as e: print(f"Warning: Failed to read CSV: {e}. Using code defaults.") pass # ファイルが存在しない、または読み込み失敗した場合、デフォルト値でDataFrameを作成 data = [] for name, info in fit_params.items(): row = { 'Parameter': name, 'optid': info['optid'], 'convert': info['convert'], 'eps': info['eps'], 'min': info['min'] if info['min'] is not None else np.nan, 'max': info['max'] if info['max'] is not None else np.nan, 'kpenalty': info['kpenalty'], 'Value': info['x0'], } data.append(row) df = pd.DataFrame(data).set_index('Parameter') print(f"Status: Initialized default parameters.") return df def save_params(prm_csv_path, df_params): """ DataFrame形式のパラメータ設定をCSVファイルに保存。 """ try: df_params.to_csv(prm_csv_path, index=True, float_format='%.6e') print(f"Status: Parameters saved to '{os.path.basename(prm_csv_path)}'.") except Exception as e: print(f"Error: Failed to save parameters to CSV: {e}") traceback.print_exc() def update_params(args, df_params): """ コマンドライン引数、CSVファイルから全パラメータの初期値を決定し、辞書で返す。 """ initial_params_dict = {} for name, info in fit_params.items(): # 1. CSVから 'Value' を取得 val = df_params.loc[name, 'Value'] # 2. コマンドライン引数で上書き arg_val = getattr(args, name, None) if arg_val is not None: val = arg_val initial_params_dict[name] = float(val) if initial_params_dict["w"] <= 0 or initial_params_dict["I0"] <= 0: raise ValueError("w (FWHM) and I0 must be positive.") w0 = initial_params_dict["w"] initial_params_dict["kG"] = fwhm2kG(w0) return initial_params_dict def print_params(message, params_dict): """全パラメータを出力""" print(message) w = params_dict.get('w', None) kG = fwhm2kG(w) if w is not None and w > 0 else 0.0 for name in ALL_PARAM_NAMES: val = params_dict.get(name) print(f" {name}={val:.6f}") if w is not None: print(f" w ={w:.6f} (kG={kG:.6f})") # ------------------------------- # スペクトルデータ読み込み # ------------------------------- def load_spectrum(filename, Escale = -1.0): """ データファイルを読み込み、エネルギーEと強度IをNumPy配列として返す。 """ ext = os.path.splitext(filename)[1].lower() if ext == ".csv": df = pd.read_csv(filename) elif ext == ".xlsx": df = pd.read_excel(filename) elif ext in [".txt", ".dat"]: df = pd.read_csv(filename, skiprows=1) else: raise ValueError(f"未対応のファイル形式: {ext}") df.columns = ["Energy (eV)", "Measured Intensity"] E = df["Energy (eV)"].values.astype(float) * Escale I = df["Measured Intensity"].values.astype(float) return E, I # ------------------------------- # モデルスペクトル（カーネル定義） # ------------------------------- def fermi(E, EF, T): kB = 8.617333e-5 return 1.0 / (np.exp((E - EF) / (kB * T)) + 1.0) def DOS(E, DOSa): """状態密度 (ここでは線形傾きを仮定: (1.0 - DOSa * E))""" return (1.0 - DOSa * E) def app_func(E_kernel, kG): """ガウス関数カーネル (プロット用)""" return np.exp(-kG * E_kernel**2) def gaussian_kernel_func(dE, w): """FWHM wに基づくガウス関数カーネル値""" kG = fwhm2kG(w) return np.exp(-kG * dE**2) def convolute_by_func(E, I_phys, w): """ 要素ごとのループを使って、I_physにガウス関数カーネルを畳み込む。 両端でのカーネルカットオフを、その位置での重み合計(wtot)で除算して補正する。 """ if w <= 1e-10: return I_phys ndata = len(E) dx = E[1] - E[0] di = int( (w * 5.0) / dx + 1.1 ) I_conv = np.zeros(ndata, dtype=float) for j in range(ndata): wtot = 0.0 for k in range(-di, di + 1): idx_phys = j - k if idx_phys < 0 or idx_phys >= ndata: continue dE = k * dx G_kernel = gaussian_kernel_func(dE, w) weight = G_kernel * dx wtot += weight I_conv[j] += I_phys[idx_phys] * weight if wtot > 1e-10: I_conv[j] /= wtot else: I_conv[j] = 0.0 return I_conv def model_spectrum(E, EF, BG, kG, I0, DOSa, T): w = kG2fwhm(kG) D = DOS(E, DOSa) f = fermi(E, EF, T) I_phys = D * f conv = convolute_by_func(E, I_phys, w) return I0 * conv + BG def calS(opt_params, E, I, T, initial_params_dict, termination_state): """ 目的関数 (残差二乗和 S + ペナルティ項)。 termination_state が True の場合は、キャッシュされた値を返す。 """ # 早期停止が確定している場合、キャッシュされた最終エラー値を返す if termination_state[0]: return termination_state[1] # 1. 全パラメータを元のスケールに復元 EF, BG, w, I0, DOSa = recover_params(opt_params, initial_params_dict) kG = fwhm2kG(w) if w > 0 else 1e-10 # 2. モデル計算 try: I_model = model_spectrum(E, EF, BG, kG, I0, DOSa, T) S = np.sum((I - I_model) ** 2) except Exception as e: S = 1e30 # 3. ペナルティ項の計算 (変更なし) penalty = 0.0 p_names = ["EF", "BG", "w", "I0", "DOSa"] p_values = [EF, BG, w, I0, DOSa] for name, value in zip(p_names, p_values): p_info = fit_params.get(name) if not p_info: continue kpenalty = p_info["kpenalty"] p_min = p_info["min"] p_max = p_info["max"] if (name == 'w' and value <= 0) or (name == 'I0' and value <= 0): penalty += kpenalty * 1e6 continue if p_min is not None and value < p_min: penalty += kpenalty * (p_min - value)**2 if p_max is not None and value > p_max: penalty += kpenalty * (value - p_max)**2 return S + penalty # ------------------------------- # 初期値推定モード # ------------------------------- def init_params(args, df_params): print("\n--- Guess initial parameters ---") print(f"Efit_min={args.Efit_min}, Efit_max={args.Efit_max}, T={args.T}K") E, I = load_spectrum(args.filename, args.Escale) if args.Efit_min is None or args.Efit_max is None: raise ValueError("mode='init' requires both --Efit_min and --Efit_max to be set.") mask = (E >= args.Efit_min) & (E <= args.Efit_max) E_fit, I_fit = E[mask], I[mask] if len(E_fit) < 2: raise ValueError("Fit range is too small or invalid.") idx_max = np.argmin(np.abs(E - args.Efit_max)) BG_init = I[idx_max] idx_min = np.argmin(np.abs(E - args.Efit_min)) I_min = I[idx_min] I0_init = I_min - BG_init if I0_init <= 0: I0_init = fit_params['I0']['x0'] target_intensity = BG_init + (I0_init / 2.0) intensity_diff = np.abs(I_fit - target_intensity) idx_ef = np.argmin(intensity_diff) EF_init = E_fit[idx_ef] kB = 8.617333e-5 w_thermal = 3.5 * kB * args.T w_default = fit_params['w']['x0'] w_init = max(w_default, w_thermal * 1.5) if fit_params['w']['max'] is not None and w_init > fit_params['w']['max']: w_init = fit_params['w']['max'] DOSa_init = fit_params['DOSa']['x0'] initial_params_dict = { "EF": EF_init, "BG": BG_init, "w": w_init, "I0": I0_init, "DOSa": DOSa_init } print_params("\nGuessed parameters:", initial_params_dict) # 推定値をdf_paramsのValue列に書き込み for name, value in initial_params_dict.items(): if name in df_params.index: df_params.loc[name, 'Value'] = value prm_csv_path = get_prm_csv_path(args.filename) save_params(prm_csv_path, df_params) print(f"\nGuessed parameters saved to [{os.path.basename(prm_csv_path)}]") # プロット kG_init = fwhm2kG(w_init) DOSa_init = initial_params_dict['DOSa'] I_model = model_spectrum(E, EF_init, BG_init, kG_init, I0_init, DOSa_init, args.T) plt.figure() plt.plot(E, I, 'o', label="measured") plt.plot(E, I_model, '-', label="guessed") if args.Efit_min: plt.axvline(args.Efit_min, color='gray', linestyle='--', linewidth = 0.5) if args.Efit_max: plt.axvline(args.Efit_max, color='gray', linestyle='--', linewidth = 0.5) plt.axvline(EF_init, color='red', linestyle=':', label="EF (guessed)", linewidth = 0.5) plt.axhline(BG_init + I0_init, color='green', linestyle=':', label="Max Intensity (guessed)", linewidth = 0.5) plt.axhline(BG_init, color='blue', linestyle=':', label="BG (guessed)", linewidth = 0.5) plt.xlabel("Energy (eV)") plt.ylabel("Intensity") plt.legend() plt.title("Guess result") plt.tight_layout() plt.pause(0.1) # ------------------------------- # フィッティング処理 # ------------------------------- def fit(args, prm_csv_path, df_params): E, I = load_spectrum(args.filename, args.Escale) # フィット範囲 mask = np.ones_like(E, dtype=bool) if args.Efit_min is not None: mask &= E >= args.Efit_min if args.Efit_max is not None: mask &= E <= args.Efit_max E_fit, I_fit = E[mask], I[mask] # 全パラメータの初期値設定 initial_params_dict = update_params(args, df_params) # 最適化対象のみを抽出し、logスケールに変換 log_params0 = convert_params(initial_params_dict) print_params("\nInitial parameters:", initial_params_dict) history = [] step_count = 0 # --- 外部停止フラグ [is_terminated (bool), last_error (float)] --- # minimizeの引数として渡すためリストを使用 termination_state = [False, 0.0] # ------------------------------------------------------------- # --------------------------------------------- # アニメーション初期化 # --------------------------------------------- fig, ax = plt.subplots() ax.plot(E, I, 'o', label="measured", markersize=2) kG0 = initial_params_dict['kG'] DOSa0 = initial_params_dict['DOSa'] EF0 = initial_params_dict['EF'] BG0 = initial_params_dict['BG'] I00 = initial_params_dict['I0'] I_model_init = model_spectrum(E, EF0, BG0, kG0, I00, DOSa0, args.T) # 初期スペクトルをプロット ax.plot(E, I_model_init, '--', color='gray', label="initial") # フィッティングカーブ用のラインを初期化 line_fit, = ax.plot(E, I_model_init, '-', color='red', label="optimizing...") if args.Efit_min: ax.axvline(args.Efit_min, color='gray', linestyle='--', linewidth = 0.5) if args.Efit_max: ax.axvline(args.Efit_max, color='gray', linestyle='--', linewidth = 0.5) ax.set_xlabel("Energy (eV)") ax.set_ylabel("Intensity") ax.legend() plt.ion() plt.pause(0.1) # --------------------------------------------- prev_opt_params = np.copy(log_params0) def callback(opt_params): """最適化ステップごとのコールバック関数 (eps停止条件をチェック)""" nonlocal step_count, prev_opt_params, termination_state # 停止が確定していたら、計算負荷の高い残りの処理をスキップしTrueを返す if termination_state[0]: return True step_count += 1 # 1. 全パラメータを元のスケールに復元 EF, BG, w, I0, DOSa = recover_params(opt_params, initial_params_dict) # 履歴に保存するデータを作成 (物理量のリスト) # calSの計算はここで行う (calS内でtermination_stateのチェックは行わない) err = calS(opt_params, E_fit, I_fit, args.T, initial_params_dict, [False, 0.0]) # フラグは無視 all_param_values = [EF, BG, w, I0, DOSa] history.append([step_count, *all_param_values, err]) # --------------------------------------------- # 2. 停止条件のチェック (eps判定) # --------------------------------------------- if step_count > 1: diff_satisfied = True current_param_values = np.array(all_param_values) prev_all_params_list = recover_params(prev_opt_params, initial_params_dict) prev_param_values = np.array(prev_all_params_list) param_names = ALL_PARAM_NAMES for i, name in enumerate(param_names): if fit_params[name]['optid'] == 0: continue eps = fit_params[name]['eps'] delta = np.abs(current_param_values[i] - prev_param_values[i]) if delta >= eps: diff_satisfied = False break # 早期停止条件を満たした場合、フラグを立てて最後のエラー値をキャッシュ if diff_satisfied: termination_state[0] = True termination_state[1] = err print(f"[Step {step_count:4d}] EF={EF:.5f}, BG={BG:.5f}, w={w:.5f}, I0={I0:.5f}, DOSa={DOSa:.5f}, Error={err:.6e}") print(f"\nOptimization terminated successfully: All fitted parameters changed by less than their respective EPS value (Step {step_count}).") return True # 即座に停止シグナルを返す # 次ステップのために現在のパラメータを保存 prev_opt_params[:] = opt_params[:] # 1. コンソール出力: print_intervalごと if step_count == 1 or (step_count % args.print_interval == 0): print(f"[Step {step_count:4d}] EF={EF:.5f}, BG={BG:.5f}, w={w:.5f}, I0={I0:.5f}, DOSa={DOSa:.5f}, Error={err:.6e}") # 2. アニメーション更新 if step_count % (args.print_interval // 10 + 1) == 0: kG_c = fwhm2kG(w) I_model_current = model_spectrum(E, EF, BG, kG_c, I0, DOSa, args.T) line_fit.set_ydata(I_model_current) ax.set_title(f"Step {step_count} / Error: {err:.2e}") fig.canvas.draw() fig.canvas.flush_events() plt.pause(0.001) return False # 継続 # calSの引数に termination_state を追加 options = {'maxiter': args.nmaxiter} result = minimize(calS, log_params0, args=(E_fit, I_fit, args.T, initial_params_dict, termination_state), method=args.method, callback=callback, options=options) # --------------------------------------------- # 終了後の処理 # --------------------------------------------- plt.ioff() # 最適解を元のスケールに復元 EF_fit, BG_fit, w_fit, I0_fit, DOSa_fit = recover_params(result.x, initial_params_dict) kG_fit = fwhm2kG(w_fit) final_params_dict = { "EF": EF_fit, "BG": BG_fit, "w": w_fit, "I0": I0_fit, "DOSa": DOSa_fit } print_params("Optimized", final_params_dict) # パラメータCSVの更新: Valueカラムに最適化結果を上書き for name, value in final_params_dict.items(): if name in df_params.index: df_params.loc[name, 'Value'] = value print() save_params(prm_csv_path, df_params) # 最終モデル計算 I_model_final = model_spectrum(E, EF_fit, BG_fit, kG_fit, I0_fit, DOSa_fit, args.T) # 最終結果の表示 line_fit.set_ydata(I_model_final) line_fit.set_color('blue') line_fit.set_label('optimized') ax.legend() ax.set_title("Fitting completed") fig.canvas.draw() plt.pause(0.1) # Excel出力 (NumPy配列から直接DataFrameを作成) df_output = pd.DataFrame({ "Energy (eV)": E, "Measured": I, "Initial": I_model_final, "Final": I_model_init, }) excel_path = os.path.splitext(args.filename)[0] + "_fit.xlsx" df_output.to_excel(excel_path, index=False) print(f"Fitting results saved to [{excel_path}]") # 履歴保存 (変更なし) history_columns = ["Step"] + ALL_PARAM_NAMES + ["Error"] hist_array = np.array(history) hist_df = pd.DataFrame(hist_array, columns=history_columns) hist_df['kG'] = fwhm2kG(hist_df['w']) output_cols = ["Step"] + ALL_PARAM_NAMES + ["Error"] output_hist_df = hist_df[output_cols] hist_path = os.path.splitext(args.filename)[0] + "_history.xlsx" output_hist_df.to_excel(hist_path, index=False) print(f"History saved to [{hist_path}]") # ------------------------------- # シミュレーションモード # ------------------------------- def simulate(args, df_params): E, I = load_spectrum(args.filename, args.Escale) initial_params_dict = update_params(args, df_params) EF0 = initial_params_dict['EF'] BG0 = initial_params_dict['BG'] w0 = initial_params_dict['w'] I00 = initial_params_dict['I0'] DOSa0 = initial_params_dict['DOSa'] kG0 = initial_params_dict['kG'] if w0 <= 0 or I00 <= 0: raise ValueError("w (FWHM) and I0 must be positive for simulation.") print() print("Simulation parameters:") print_params("使用パラメータ", initial_params_dict) I_model = model_spectrum(E, EF0, BG0, kG0, I00, DOSa0, args.T) D = DOS(E, DOSa0) f = fermi(E, EF0, args.T) E_kernel = np.linspace(-len(E) // 2 * (E[1]-E[0]), (len(E) // 2 - 1 + (len(E) % 2)) * (E[1]-E[0]), len(E)) kG_for_plot = fwhm2kG(w0) G = app_func(E_kernel, kG_for_plot) G /= np.max(G) plt.figure() plt.plot(E, I, label="obs", linestyle = '', marker = 'o', markersize = 2.0) plt.plot(E, I_model, label="sim", linestyle = '-') plt.plot(E, 0.5*I00 * D, label="DOS(E)", linestyle = '-', linewidth = 0.5) plt.plot(E, 0.5*I00 * f, label="fFD(E)", linestyle = '-', linewidth = 0.5) plt.plot(E, 0.5*I00 * D * f, label="DOS(E)*fFD(E)", linestyle = '-', linewidth = 0.5) plt.plot(E, 0.5*I00 * G, label="g_app(E)", linestyle = '-', linewidth = 0.5) if args.Efit_min: plt.axvline(args.Efit_min, color='gray', linestyle='--', linewidth = 0.5) if args.Efit_max: plt.axvline(args.Efit_max, color='gray', linestyle='--', linewidth = 0.5) plt.xlabel("Energy (eV)") plt.ylabel("Intensity") plt.legend() plt.tight_layout() plt.pause(0.1) # ------------------------------- # メイン関数 # ------------------------------- def main(): args = initialize() if (args.mode != "init" and (args.w is not None and args.w <= 0 or args.I0 is not None and args.I0 <= 0)): raise ValueError("w (FWHM) and I0 must be positive in command line arguments.") prm_csv_path = get_prm_csv_path(args.filename) df_params = load_params(prm_csv_path) if args.mode == "init": init_params(args, df_params) elif args.mode == "fit": fit(args, prm_csv_path, df_params) elif args.mode == "sim": simulate(args, df_params) else: raise ValueError("modeは 'fit', 'sim', 'init' のいずれかを指定してください") if __name__ == "__main__": main() terminate()