import chardet 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": None, "max": None, "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=str, default=None, help="フィット範囲下限 [eV]") parser.add_argument("--Efit_max", type=str, 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 read_haxpes_spring8(file_path): """ エネルギー(E)と平均信号強度(I)のリストを返します。 """ try: with open(file_path, 'rb') as f: raw_data = f.read() except FileNotFoundError: print(f"Error: File not found: {file_path}") return [], [] except Exception as e: print(f"Error reading/decoding file: {e}") return [], [] detected = chardet.detect(raw_data) encoding = detected['encoding'] confidence = detected['confidence'] # print(f"Detected encoding: {encoding} (Confidence: {confidence:.2f})") if encoding is None: print("Warning: Failed to detect character encoding. Use python default") # encoding = 'utf-8' lines = [] try: content = raw_data.decode(encoding) lines = content.splitlines() except Exception as e: print(f"Error decoding data: {e}") return [], [] # --------------------------------------------------------- # 2. Typeチェック # --------------------------------------------------------- if not lines or '[Info]' not in lines[0]: print("Error: Invalid file format. Header '[Info]' missing in the first line.") return [], [] # --------------------------------------------------------- # 3. ヘッダー情報の解析 # --------------------------------------------------------- dim1_size = 0 # データ行数 dim2_size = 0 # データ列数 data_start_idx = -1 for i, line in enumerate(lines): line = line.strip() if line.startswith('Dimension 1 size='): dim1_size = int(line.split('=')[1]) elif line.startswith('Dimension 2 size='): dim2_size = int(line.split('=')[1]) elif line == '[Data 1]': data_start_idx = i + 1 break if dim1_size == 0: print("Error: Failed to parse dimension1 size.") return [], [] if dim2_size == 0: dim2_size = 1 if data_start_idx == -1: print("Error: Failed to find [Data 1] tag.") return [], [] # --------------------------------------------------------- # 4. データ領域の処理 # --------------------------------------------------------- E_list = [] I_list = [] # 必要な行数だけ切り出す data_block = lines[data_start_idx : data_start_idx + dim1_size] for row_idx, line in enumerate(data_block): parts = line.split() # 列数チェック expected_cols = 1 + dim2_size if len(parts) < expected_cols: print(f"Warning: Line {data_start_idx + row_idx + 1} insufficient columns.") continue try: # 1列目: エネルギー energy = float(parts[0]) # 2列目以降: 信号強度 intensities = [float(x) for x in parts[1 : expected_cols]] # 平均化 avg_intensity = sum(intensities) / len(intensities) E_list.append(energy) I_list.append(avg_intensity) except ValueError: print(f"Warning: Conversion error at line {data_start_idx + row_idx + 1}.") continue return np.array(E_list, dtype='float64'), np.array(I_list, dtype='float64') def read_hapes_xy(file_path, read_ndata = True): try: with open(file_path, 'rb') as f: raw_data = f.read() except FileNotFoundError: print(f"Error: File not found: {file_path}") return [], [] except Exception as e: print(f"Error reading/decoding file: {e}") return [], [] detected = chardet.detect(raw_data) encoding = detected['encoding'] confidence = detected['confidence'] # print(f"Detected encoding: {encoding} (Confidence: {confidence:.2f})") if encoding is None: print("Warning: Failed to detect character encoding. Use python default") # encoding = 'utf-8' lines = [] try: content = raw_data.decode(encoding) lines = content.splitlines() except Exception as e: print(f"Error decoding data: {e}") return [], [] if read_ndata: idx = 2 else: idx = 1 E_list = [] I_list = [] while True: if lines[idx] is None or lines[idx] == "": break _aa = lines[idx].split() if len(_aa) < 2: break print("_aa=", _aa) E_list.append(float(_aa[0])) I_list.append(float(_aa[1])) idx += 1 return np.array(E_list, dtype='float64'), np.array(I_list, dtype='float64') def read_haxpes_txt(file_path): try: with open(file_path, 'r') as f: line1 = f.readline() line2 = f.readline() line3 = f.readline() except FileNotFoundError: print(f"Error: File not found: {file_path}") return [], [] except Exception as e: print(f"Error reading/decoding file: {e}") return [], [] if '[Info]' in line1: return read_haxpes_spring8(file_path) if re.match(r'\s*\d+\s*$', line2): return read_hapes_xy(file_path, read_ndata = True) return read_hapes_xy(file_path, read_ndata = False) def load_spectrum(filename, Escale = -1.0): """ データファイルを読み込み、エネルギーEと強度IをNumPy配列として返す。 """ ext = os.path.splitext(filename)[1].lower() if ext == ".txt": E, I = read_haxpes_txt(filename) else: 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 # eV/K ke = (E - EF) / (kB * T) # 数値的に安全な範囲に制限 ke = np.clip(ke, -30, 30) return 1.0 / (np.exp(ke) + 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 kernel_func(dE, plist): w = plist[0] kG = fwhm2kG(w) ret = np.exp(-kG * dE**2) # print("w=", w, ret[0:5]) return ret def convolute_by_func_original(E, I_phys, w): if w <= 1e-10: return I_phys ndata = len(E) dx = np.abs(E[1] - E[0]) di = int( (w * 5.0) / abs(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 * abs(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 # ========================================== # 3. NumPy Optimized (修正版) # ========================================== def convolute_by_func_np(E, I_phys, w, kernel_range = 5.0): if w <= 1e-10: return I_phys ndata = len(E) dx = np.abs(E[1] - E[0]) di = int( (w * kernel_range) / dx + 1.1 ) if di < 1: di = 1 k_indices = np.arange(-di, di + 1) dE_array = k_indices * dx weights = kernel_func(dE_array, [w]) * dx I_conv_unnorm = np.convolve(I_phys, weights, mode='same') W_norm = np.convolve(np.ones(ndata), weights, mode='same') # 配列サイズ補正 if I_conv_unnorm.shape[0] > ndata: start = (I_conv_unnorm.shape[0] - ndata) // 2 I_conv_unnorm = I_conv_unnorm[start : start + ndata] W_norm = W_norm[start : start + ndata] I_conv = np.divide(I_conv_unnorm, W_norm, out=np.zeros_like(I_conv_unnorm), where=W_norm > 1e-10) 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 + ペナルティ項)。 """ # 早期停止が確定している場合 if termination_state[0]: return termination_state[1] # 1. 全パラメータを元のスケールに復元 EF, BG, w, I0, DOSa = recover_params(opt_params, initial_params_dict) # ------------------------------------------------------------- # ガード処理 # ------------------------------------------------------------- # [Check 1] パラメータ自体が NaN や Inf になっていないか確認 if not np.all(np.isfinite([EF, BG, w, I0, DOSa])): return 1e30 # [Check 2] w (FWHM) の物理的妥当性と kG のオーバーフロー防止 # w が 0 に近すぎると、kG = log(2)/(w/2)^2 が無限大になり、exp(-kG*x^2) の計算で不具合が出る # また、I0 が負である場合も物理的にありえないため除外 min_w_limit = 1e-5 # これより細い幅は計測限界以下として弾く if w < min_w_limit or I0 <= 1e-10: return 1e30 # [Check 3] kG の算出 kG = fwhm2kG(w) # 万が一 kG が巨大になりすぎた場合のガード（double精度の限界考慮） if kG > 1e20: return 1e30 # ------------------------------------------------------------- # 2. モデル計算 (エラーチェック済みなので直接実行) # ------------------------------------------------------------- I_model = model_spectrum(E, EF, BG, kG, I0, DOSa, T) # [Check 4] 計算結果に NaN/Inf が含まれていないか最終確認 # (convolveの結果などで万が一発生した場合用) if not np.all(np.isfinite(I_model)): return 1e30 # 残差二乗和の計算 S = np.sum((I - I_model) ** 2) # ------------------------------------------------------------- # 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"] # w, I0 の正値制約ペナルティ (念のため残すが、上流で弾いているので基本到達しない) 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 type(args.Efit_min) is str: try: args.Efit_min = float(args.Efit_min) except: args.Efit_min = min(E) if type(args.Efit_max) is str: try: args.Efit_max = float(args.Efit_max) except: args.Efit_max = max(E) print() print("Fitting:") print(f"Escale : {args.Escale}") print(f"Efit_min: {args.Efit_min}") print(f"Efit_max: {args.Efit_max}") 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) if type(args.Efit_min) is str: try: args.Efit_min = float(args.Efit_min) except: args.Efit_min = min(E) if type(args.Efit_max) is str: try: args.Efit_max = float(args.Efit_max) except: args.Efit_max = max(E) print() print("Fitting:") print(f"Escale : {args.Escale}") print(f"Efit_min: {args.Efit_min}") print(f"Efit_max: {args.Efit_max}") print(f"method : {args.method}") print(f"nmaxiter: {args.nmaxiter}") # フィット範囲 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 not early_stop: diff_satisfied = False # 早期停止条件を満たした場合、フラグを立てて最後のエラー値をキャッシュ 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ごと is_print = (step_count % args.print_interval == 0) if step_count == 1 or is_print: 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: if is_print: 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) if type(args.Efit_min) is str: try: args.Efit_min = float(args.Efit_min) except: args.Efit_min = min(E) if type(args.Efit_max) is str: try: args.Efit_max = float(args.Efit_max) except: args.Efit_max = max(E) print() print("Fitting:") print(f"Escale : {args.Escale}") print(f"Efit_min: {args.Efit_min}") print(f"Efit_max: {args.Efit_max}") 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(): global convolute_by_func, early_stop convolute_by_func = convolute_by_func_np # convolute_by_func = convolute_by_func_original # early_stop = False early_stop = True args = initialize() if args.method.lower() == 'simplex': args.method = 'nelder-mead' if args.method.lower() == 'nelder-mead': early_stop = False 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()