#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
fit_vbo_xps.py
概要: XPS価電子帯スペクトルをフィッティングするためのスクリプト。
詳細説明:
フィルム/基板混合XPS価電子帯スペクトルを、参照フィルムスペクトルと参照基板スペクトルの
重み付き合計、および多項式背景の和としてフィッティングします。
各参照スペクトルには独立したエネルギーシフトが適用されます。
フィッティングはNelder-Mead最適化手法を使用し、堅牢な損失関数とソフトバウンド制約をサポートします。
モデル:
I_mix(E) = A_f * I_film(E - dE_f) + A_s * I_sub(E - dE_s)
+ polynomial_background(E)
フィッティングされたオフセットは主に以下として報告されます。
dE_shift = dE_f - dE_s
参照VBM位置が提供される場合、VBOは以下として計算されます。
VBO_BE = (VBM_film_ref + dE_f) - (VBM_sub_ref + dE_s)
結合エネルギー軸については、符号の慣例に注意が必要です。
スクリプトは結合エネルギー差とその負の値の両方を出力します。
例:
python fit_vbo_xps.py ^
--film Bi2OS2100nm_20260508-1_Bi2OS2_pSi_HAXPES_fit_YBPS100_VBwideY2-deconvoluted_point13.xlsm ^
--substrate n-Si_20260508-1_Bi2OS2_pSi_HAXPES_fit_YNS_VBnarrowY-Si-deconvoluted_point10.xlsm ^
--mix TOA90_20260508-1_Bi2OS2_pSi_HAXPES_fit_YBPS25-90_VBwideY2-deconvoluted_point13.xlsm ^
--fit-xmin -2.5 --fit-xmax 8 ^
--show 1 --amp-sub0 1.0
関連リンク: :doc:`fit_vbo_xps_usage`
"""
from __future__ import annotations
import argparse
import os
import sys
import traceback
from dataclasses import dataclass
from typing import Tuple, List
try:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.optimize import minimize
except Exception:
print("ERROR: required library import failed.")
traceback.print_exc()
print("\nInstall example:")
print(" pip install numpy pandas matplotlib scipy openpyxl")
input("\nPress ENTER to terminate>>\n")
sys.exit(1)
[ドキュメント]
@dataclass
class Spectrum:
"""
概要: XPSスペクトルデータを保持するデータクラス。
詳細説明:
スペクトルの名前、エネルギー（または結合エネルギー）軸のデータ、
および強度（またはカウント）軸のデータを格納します。
Attributes:
name (str): スペクトルの識別名（通常はファイル名）。
x (np.ndarray): エネルギー（または結合エネルギー）軸のデータ。
y (np.ndarray): 強度（またはカウント）軸のデータ。
"""
name: str
x: np.ndarray
y: np.ndarray
[ドキュメント]
def read_spectrum(path: str, sheet: str | int = 0, x_col: int = 0, y_col: int = 2,
skip_first_blank: int = 1) -> Spectrum:
"""
概要: XLSX/XLSM/CSVファイルからXPSスペクトルデータを読み込む。
詳細説明:
指定されたファイルパス、シート、X列、Y列からデータを読み込みます。
数値に変換できないデータは除外され、NaN値も処理されます。
X軸データは昇順にソートされ、重複するX値のYデータは平均化されます。
:param path: 読み込むファイルへのパス。
:param sheet: Excelファイルの場合のシート名またはゼロベースのシートインデックス。
CSV/TXTファイルの場合は無視されます。デフォルトは0。
:param x_col: X軸データのゼロベースの列インデックス。デフォルトは0。
:param y_col: Y軸データ（強度）のゼロベースの列インデックス。
デフォルトは2（供給されるファイルではデコンボリューションされたデータ）。
:param skip_first_blank: 最初の行が空の場合にスキップするかどうか（1: スキップ、0: スキップしない）。
供給されるファイルでは通常空の最初の行があるため、デフォルトは1。
:returns: 読み込まれたデータを含む :class:`Spectrum` オブジェクト。
:raises ValueError: 有効なデータポイントが少なすぎる場合。
"""
ext = os.path.splitext(path)[1].lower()
if ext in [".csv", ".txt", ".dat"]:
df = pd.read_csv(path, header=None)
else:
df = pd.read_excel(path, sheet_name=sheet, header=None, engine="openpyxl")
if skip_first_blank:
# The supplied files have an empty first row. Numeric conversion below also removes it.
pass
x = pd.to_numeric(df.iloc[:, x_col], errors="coerce").to_numpy(dtype=float)
y = pd.to_numeric(df.iloc[:, y_col], errors="coerce").to_numpy(dtype=float)
m = np.isfinite(x) & np.isfinite(y)
x = x[m]
y = y[m]
if len(x) < 5:
raise ValueError(f"too few valid data points in {path}")
# sort by increasing x and merge duplicate x by averaging y
order = np.argsort(x)
x = x[order]
y = y[order]
xu, inv = np.unique(x, return_inverse=True)
if len(xu) != len(x):
yu = np.zeros_like(xu, dtype=float)
cnt = np.zeros_like(xu, dtype=float)
for i, k in enumerate(inv):
yu[k] += y[i]
cnt[k] += 1
x, y = xu, yu / cnt
return Spectrum(os.path.basename(path), x, y)
[ドキュメント]
def crop(sp: Spectrum, xmin: float | None, xmax: float | None) -> Spectrum:
"""
概要: Spectrumオブジェクトのデータを指定されたX軸範囲でトリミングする。
詳細説明:
`xmin` と `xmax` の値に基づいて、`sp` オブジェクトのX軸とY軸のデータをフィルタリングし、
新しい :class:`Spectrum` オブジェクトを返します。
:param sp: トリミング対象の :class:`Spectrum` オブジェクト。
:param xmin: X軸の最小値。`None` の場合は下限なし。
:param xmax: X軸の最大値。`None` の場合は上限なし。
:returns: トリミングされたデータを含む新しい :class:`Spectrum` オブジェクト。
"""
m = np.ones_like(sp.x, dtype=bool)
if xmin is not None:
m &= sp.x >= xmin
if xmax is not None:
m &= sp.x <= xmax
return Spectrum(sp.name, sp.x[m], sp.y[m])
[ドキュメント]
def make_interp(sp: Spectrum, fill_value: float = 0.0):
"""
概要: Spectrumオブジェクトから線形補間関数を作成する。
詳細説明:
`scipy.interpolate.interp1d` を使用して、指定された :class:`Spectrum` オブジェクトの
X軸とY軸データに基づいて線形補間関数を構築します。
補間範囲外の値は `fill_value` で埋められます。
:param sp: 補間関数を作成する基になる :class:`Spectrum` オブジェクト。
:param fill_value: 補間範囲外のX値に対して返される値。デフォルトは0.0。
:returns: 作成された線形補間関数 (:class:`scipy.interpolate.interp1d` オブジェクト)。
"""
return interp1d(sp.x, sp.y, kind="linear", bounds_error=False, fill_value=fill_value,
assume_sorted=True)
[ドキュメント]
def robust_scale(y: np.ndarray) -> float:
"""
概要: データのロバストなスケール（強度範囲）を計算する。
詳細説明:
入力配列 `y` の95パーセンタイルと5パーセンタイルの差を計算することで、
外れ値に強い強度スケールを推定します。
計算されたスケールが非数値または非正の場合、代わりに絶対値の最大値を使用します。
:param y: スケールを計算する対象の強度データ配列。
:returns: 計算されたロバストなスケール値。
"""
s = np.nanpercentile(y, 95) - np.nanpercentile(y, 5)
if not np.isfinite(s) or s <= 0:
s = np.nanmax(np.abs(y))
return float(s if s > 0 else 1.0)
[ドキュメント]
def poly_background(x: np.ndarray, coeff: np.ndarray, x0: float, xscale: float) -> np.ndarray:
"""
概要: 指定された係数を持つ多項式背景を生成する。
詳細説明:
X軸データ `x` に基づいて多項式背景を計算します。
多項式の評価には、安定性のために `x0` と `xscale` を用いてX軸が標準化されます。
例えば、`coeff` が `[c0, c1, c2]` の場合、背景は `c0 + c1*(x-x0)/xscale + c2*((x-x0)/xscale)^2` となります。
:param x: 多項式背景を計算するX軸データ。
:param coeff: 多項式の係数配列（例: `[c0, c1, c2]` は `c0 + c1*xr + c2*xr^2` に対応）。
:param x0: X軸のオフセット（中心点）。
:param xscale: X軸のスケール因子。
:returns: 計算された多項式背景のY値配列。
"""
xr = (x - x0) / xscale
bg = np.zeros_like(xr)
for i, c in enumerate(coeff):
bg += c * xr**i
return bg
[ドキュメント]
def normalize_intensity(y: np.ndarray, mode: str = "max") -> np.ndarray:
"""
概要: スペクトル強度データを特定のモードで正規化する。
詳細説明:
入力 `y` 配列の強度を、視覚的な比較のために指定されたモードで正規化します。
- `"max"`: 強度の絶対値の最大値で割って正規化します（デフォルト）。
- `"minmax"`: (y - y_min) / (y_max - y_min) で正規化し、値を0から1の範囲にスケーリングします。
- `"area"`: 強度の絶対値の積分（台形法）で割って正規化します。
:param y: 正規化する強度データ配列。
:param mode: 正規化モード。"max", "minmax", または "area" のいずれか。デフォルトは"max"。
:returns: 正規化された強度データ配列。
"""
yy = np.asarray(y, dtype=float)
if len(yy) == 0:
return yy.copy()
if mode == "minmax":
ymin = np.nanmin(yy)
ymax = np.nanmax(yy)
scale = ymax - ymin
if not np.isfinite(scale) or scale <= 0:
scale = 1.0
return (yy - ymin) / scale
if mode == "area":
# For area normalization, keep the sign of the spectrum and normalize
# by the integral of abs(y) over the available point index axis.
scale = np.trapz(np.abs(yy))
if not np.isfinite(scale) or scale <= 0:
scale = 1.0
return yy / scale
# default: maximum-height normalization without baseline subtraction.
scale = np.nanmax(np.abs(yy))
if not np.isfinite(scale) or scale <= 0:
scale = 1.0
return yy / scale
[ドキュメント]
def fit_mixture(film: Spectrum, sub: Spectrum, mix: Spectrum, args):
"""
概要: フィルムと基板の混合XPSスペクトルを最適化手法でフィッティングする。
詳細説明:
`scipy.optimize.minimize` のNelder-Mead法を使用して、混合スペクトルをフィルム、基板、
および多項式背景の線形結合でフィッティングします。
エネルギーシフト、振幅、背景係数が最適化されます。
ソフトバウンド制約と堅牢な損失関数をサポートします。
フィッティングモデルは以下の通りです。
`I_mix(E) = A_f * I_film(E - dE_f) + A_s * I_sub(E - dE_s) + poly_background(E)`
:param film: 参照フィルムスペクトルデータ。
:param sub: 参照基板スペクトルデータ。
:param mix: フィッティング対象の混合スペクトルデータ。
:param args: コマンドライン引数を格納する `argparse.Namespace` オブジェクト。
フィッティングパラメータ（初期値、制約、損失関数設定など）が含まれます。
:returns: フィッティング結果を含む辞書。
- `result` (:class:`scipy.optimize.OptimizeResult`): `scipy.optimize.minimize` の結果オブジェクト。
- `x` (:class:`np.ndarray`): フィッティング範囲内のX軸データ。
- `y` (:class:`np.ndarray`): フィッティング範囲内の観測Y軸データ。
- `yfit` (:class:`np.ndarray`): フィッティングモデルによって計算されたY軸データ。
- `film_part` (:class:`np.ndarray`): フィッティングされたフィルム成分。
- `sub_part` (:class:`np.ndarray`): フィッティングされた基板成分。
- `bg_part` (:class:`np.ndarray`): フィッティングされた背景成分。
- `params` (`dict`): 最適化されたパラメータと統計情報。
- `dE_film` (`float`): フィルムスペクトルのエネルギーシフト。
- `dE_sub` (`float`): 基板スペクトルのエネルギーシフト。
- `A_film` (`float`): フィルムスペクトルの振幅。
- `A_sub` (`float`): 基板スペクトルの振幅。
- `bg_coeff` (:class:`np.ndarray`): 背景多項式の係数。
- `rmse` (`float`): 二乗平均平方根誤差。
- `nrmse` (`float`): 正規化された二乗平均平方根誤差。
- `dE_shift_film_minus_sub` (`float`): フィルムと基板のエネルギーシフト差。
- `objective` (`float`): 最終的な目的関数の値。
- `success` (`bool`): 最適化が成功したかどうか。
- `message` (`str`): 最適化メッセージ。
- `interp` (`Tuple[interp1d, interp1d]`): フィルムと基板の補間関数。
- `x0` (`float`): 背景多項式用のX軸中心点。
- `xscale` (`float`): 背景多項式用のX軸スケール。
:raises ValueError: フィッティング範囲に有効なデータポイントが少なすぎる場合。
"""
x = mix.x.copy()
y = mix.y.copy()
if args.fit_xmin is not None or args.fit_xmax is not None:
m = np.ones_like(x, dtype=bool)
if args.fit_xmin is not None:
m &= x >= args.fit_xmin
if args.fit_xmax is not None:
m &= x <= args.fit_xmax
x, y = x[m], y[m]
if len(x) < 10:
raise ValueError("fit range contains too few points")
ff = make_interp(film, fill_value=args.outside_fill)
fs = make_interp(sub, fill_value=args.outside_fill)
yscale = robust_scale(y)
x0 = 0.5 * (x.min() + x.max())
xscale = 0.5 * (x.max() - x.min())
if xscale <= 0:
xscale = 1.0
# Initial parameters are given explicitly.
# No linear pre-fit is used because this problem has only a few parameters
# and Nelder-Mead is often more stable for this kind of spectral fitting.
bg0 = np.zeros(args.bg_order + 1, dtype=float)
bg0[0] = args.bg0
if args.bg_order >= 1:
bg0[1] = args.bg1
if args.bg_order >= 2:
bg0[2] = args.bg2
p0 = np.r_[
args.shift_film0,
args.shift_sub0,
np.log(max(args.amp_film0, args.amp_min)),
np.log(max(args.amp_sub0, args.amp_min)),
bg0,
]
def unpack(p):
d_f, d_s = p[0], p[1]
a_f, a_s = np.exp(p[2]), np.exp(p[3])
bgc = p[4:]
return d_f, d_s, a_f, a_s, bgc
def model(p, xx=x):
d_f, d_s, a_f, a_s, bgc = unpack(p)
return a_f * ff(xx - d_f) + a_s * fs(xx - d_s) + poly_background(xx, bgc, x0, xscale)
def resid(p):
yyfit = model(p)
r = (yyfit - y) / yscale
if args.loss_on_log:
yy = np.maximum(y - np.nanmin(y) + 1.0, 1.0)
mm = np.maximum(yyfit - np.nanmin(y) + 1.0, 1.0)
r = np.log(mm / yy)
return r
eval_counter = {"n": 0, "best": np.inf, "best_p": None}
def objective(p):
eval_counter["n"] += 1
d_f, d_s, a_f, a_s, bgc = unpack(p)
penalty = 0.0
# Soft bounds for unconstrained Nelder-Mead.
if abs(d_f) > args.shift_limit:
penalty += ((abs(d_f) - args.shift_limit) / max(args.shift_limit, 1.0e-12))**2 * args.penalty
if abs(d_s) > args.shift_limit:
penalty += ((abs(d_s) - args.shift_limit) / max(args.shift_limit, 1.0e-12))**2 * args.penalty
if a_f < args.amp_min:
penalty += (np.log(args.amp_min / max(a_f, 1.0e-300)))**2 * args.penalty
if a_s < args.amp_min:
penalty += (np.log(args.amp_min / max(a_s, 1.0e-300)))**2 * args.penalty
if a_f > args.amp_max:
penalty += (np.log(a_f / args.amp_max))**2 * args.penalty
if a_s > args.amp_max:
penalty += (np.log(a_s / args.amp_max))**2 * args.penalty
r = resid(p)
if args.loss == "linear":
val = float(np.mean(r * r))
elif args.loss == "logcosh":
val = float(np.mean(np.log(np.cosh(r / args.f_scale))) * args.f_scale**2)
else:
# Pseudo-Huber-like robust loss.
z = r / args.f_scale
val = float(np.mean(2.0 * args.f_scale**2 * (np.sqrt(1.0 + z*z) - 1.0)))
total = val + penalty
if total < eval_counter["best"]:
eval_counter["best"] = total
eval_counter["best_p"] = np.array(p, dtype=float).copy()
return total
callback_counter = {"iter": 0}
def callback(xk):
callback_counter["iter"] += 1
if args.print_every <= 0:
return
it = callback_counter["iter"]
if it == 1 or it % args.print_every == 0:
val = objective(xk)
d_f, d_s, a_f, a_s, bgc = unpack(xk)
print(
f"Iter {it:5d} Eval {eval_counter['n']:6d} "
f"S={val:.8g} Best={eval_counter['best']:.8g} "
f"dEf={d_f:.6g} dEs={d_s:.6g} "
f"Af={a_f:.6g} As={a_s:.6g}"
)
print("\nNelder-Mead convergence")
print(" Iter/Eval output is from the callback; S is the objective value including penalties.")
print(" Parameters: dEf, dEs [eV], Af, As")
print(
f"Initial S={objective(p0):.8g} "
f"dEf={args.shift_film0:.6g} dEs={args.shift_sub0:.6g} "
f"Af={args.amp_film0:.6g} As={args.amp_sub0:.6g}"
)
res = minimize(
objective,
p0,
method="Nelder-Mead",
callback=callback,
options={
"maxiter": args.max_iter,
"maxfev": args.max_nfev,
"xatol": args.xatol,
"fatol": args.fatol,
"disp": bool(args.verbose),
},
)
p = res.x
d_f, d_s, a_f, a_s, bgc = unpack(p)
print(
f"Final S={float(res.fun):.8g} "
f"nit={getattr(res, 'nit', -1)} nfev={getattr(res, 'nfev', -1)} "
f"dEf={d_f:.6g} dEs={d_s:.6g} Af={a_f:.6g} As={a_s:.6g}"
)
yfit = model(p)
film_part = a_f * ff(x - d_f)
sub_part = a_s * fs(x - d_s)
bg_part = poly_background(x, bgc, x0, xscale)
rmse = float(np.sqrt(np.mean((yfit - y)**2)))
nrmse = rmse / yscale
return {
"result": res,
"x": x, "y": y, "yfit": yfit,
"film_part": film_part, "sub_part": sub_part, "bg_part": bg_part,
"params": {"dE_film": d_f, "dE_sub": d_s, "A_film": a_f, "A_sub": a_s,
"bg_coeff": bgc, "rmse": rmse, "nrmse": nrmse,
"dE_shift_film_minus_sub": d_f - d_s,
"objective": float(res.fun),
"success": bool(res.success),
"message": str(res.message)},
"interp": (ff, fs),
"x0": x0, "xscale": xscale,
}
[ドキュメント]
def save_outputs(film: Spectrum, sub: Spectrum, mix: Spectrum, fit, args) -> Tuple[str, str, str, str]:
"""
概要: フィッティング結果と関連プロットをファイルに保存する。
詳細説明:
フィッティングされた混合スペクトル、個々の成分、残差を含むCSVファイル、
フィッティングパラメータとVBM情報をまとめたテキストサマリーファイル、
およびフィッティング結果と入力スペクトルの比較プロット画像をPNG形式で出力します。
:param film: 参照フィルムスペクトルデータ。
:param sub: 参照基板スペクトルデータ。
:param mix: フィッティング対象の混合スペクトルデータ。
:param fit: `fit_mixture` 関数から返されたフィッティング結果辞書。
:param args: コマンドライン引数を格納する `argparse.Namespace` オブジェクト。
出力パス、保存/表示設定、正規化モードなどが含まれます。
:returns: 生成されたCSV、サマリー、フィッティングプロットPNG、正規化プロットPNGのファイルパスのタプル。
"""
stem = args.output_stem
os.makedirs(os.path.dirname(stem) or ".", exist_ok=True)
csv_path = stem + "_fit.csv"
png_path = stem + "_fit.png"
norm_png_path = stem + "_input_normalized.png"
summary_path = stem + "_summary.txt"
x = fit["x"]
df = pd.DataFrame({
"E": x,
"I_mix_obs": fit["y"],
"I_mix_fit": fit["yfit"],
"I_film_part": fit["film_part"],
"I_substrate_part": fit["sub_part"],
"I_background": fit["bg_part"],
"residual": fit["yfit"] - fit["y"],
})
df.to_csv(csv_path, index=False)
p = fit["params"]
vbm_info = []
if args.vbm_film_ref is not None and args.vbm_sub_ref is not None:
vbm_f = args.vbm_film_ref + p["dE_film"]
vbm_s = args.vbm_sub_ref + p["dE_sub"]
vbo_be = vbm_f - vbm_s
vbm_info.append(f"VBM_film_fit_BE = {vbm_f:.8g} eV")
vbm_info.append(f"VBM_sub_fit_BE = {vbm_s:.8g} eV")
vbm_info.append(f"VBO_BE = VBM_film_BE - VBM_sub_BE = {vbo_be:.8g} eV")
vbm_info.append(f"VBO_energy_axis_opposite_sign = {-vbo_be:.8g} eV")
with open(summary_path, "w", encoding="utf-8") as f:
f.write("XPS film/substrate mixture fitting summary\n")
f.write("========================================\n")
f.write(f"film reference : {film.name}\n")
f.write(f"substrate reference: {sub.name}\n")
f.write(f"mixed spectrum : {mix.name}\n")
f.write(f"fit range : {x.min():.8g} to {x.max():.8g} eV\n")
f.write(f"dE_film : {p['dE_film']:.8g} eV\n")
f.write(f"dE_substrate : {p['dE_sub']:.8g} eV\n")
f.write(f"dE_film - dE_sub : {p['dE_shift_film_minus_sub']:.8g} eV\n")
f.write(f"A_film : {p['A_film']:.8g}\n")
f.write(f"A_substrate : {p['A_sub']:.8g}\n")
f.write(f"RMSE : {p['rmse']:.8g}\n")
f.write(f"NRMSE : {p['nrmse']:.8g}\n")
f.write(f"objective : {p['objective']:.8g}\n")
f.write(f"success : {p['success']}\n")
f.write(f"message : {p['message']}\n")
for i, c in enumerate(p["bg_coeff"]):
f.write(f"BG_coeff[{i}] : {c:.8g}\n")
if vbm_info:
f.write("\n")
f.write("\n".join(vbm_info) + "\n")
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(x, fit["y"], label="mixed obs", lw=1.5)
ax.plot(x, fit["yfit"], label="fit", lw=1.5)
ax.plot(x, fit["film_part"], label="film part", lw=1.0)
ax.plot(x, fit["sub_part"], label="substrate part", lw=1.0)
ax.plot(x, fit["bg_part"], label="background", lw=1.0)
ax.set_xlabel("Energy / Binding energy (eV)")
ax.set_ylabel("Intensity")
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
if args.save:
fig.savefig(png_path, dpi=200)
if args.show:
plt.show()
else:
plt.close(fig)
fig2, ax2 = plt.subplots(figsize=(8, 5))
ax2.plot(film.x, normalize_intensity(film.y, args.norm_mode), label="film input normalized", lw=1.2)
ax2.plot(sub.x, normalize_intensity(sub.y, args.norm_mode), label="substrate input normalized", lw=1.2)
ax2.plot(mix.x, normalize_intensity(mix.y, args.norm_mode), label="mixed input normalized", lw=1.2)
if args.fit_xmin is not None:
ax2.axvline(args.fit_xmin, ls="--", lw=0.8, alpha=0.5)
if args.fit_xmax is not None:
ax2.axvline(args.fit_xmax, ls="--", lw=0.8, alpha=0.5)
ax2.set_xlabel("Energy / Binding energy (eV)")
ax2.set_ylabel(f"Normalized intensity ({args.norm_mode})")
ax2.set_title("Input spectra comparison")
ax2.legend()
ax2.grid(True, alpha=0.3)
fig2.tight_layout()
if args.save:
fig2.savefig(norm_png_path, dpi=200)
if args.show:
plt.show()
else:
plt.close(fig2)
return csv_path, summary_path, png_path, norm_png_path
[ドキュメント]
def build_parser():
"""
概要: コマンドライン引数パーサーを構築する。
詳細説明:
XPS価電子帯スペクトルフィッティングスクリプトの実行に必要なすべてのコマンドライン引数
（入力ファイルパス、フィッティング範囲、初期パラメータ、最適化設定、出力設定など）を定義し、
設定された `ArgumentParser` オブジェクトを返します。
:returns: 設定された引数パーサーオブジェクト (:class:`argparse.ArgumentParser`)。
"""
p = argparse.ArgumentParser(description="Fit mixed XPS valence-band spectrum with film/substrate references.")
p.add_argument("--film", required=True, help="フィルムのみの参照スペクトル (xlsx/xlsm/csv)。")
p.add_argument("--substrate", required=True, help="基板のみの参照スペクトル (xlsx/xlsm/csv)。")
p.add_argument("--mix", required=True, help="フィルム/基板混合スペクトル (xlsx/xlsm/csv)。")
p.add_argument("--sheet", default=0, help="シート名またはゼロベースのシートインデックス [デフォルト: 0]。")
p.add_argument("--x-col", type=int, default=0, help="ゼロベースのX列インデックス [デフォルト: 0]。")
p.add_argument("--y-col", type=int, default=2, help="ゼロベースの強度列インデックス [デフォルト: 2; 提供ファイル: デコンボリューションされたデータ]。")
p.add_argument("--fit-xmin", type=float, default=None, help="フィッティングに使用するX軸の最小値。")
p.add_argument("--fit-xmax", type=float, default=None, help="フィッティングに使用するX軸の最大値。")
p.add_argument("--ref-xmin", type=float, default=None, help="参照スペクトルで保持するX軸の最小値。")
p.add_argument("--ref-xmax", type=float, default=None, help="参照スペクトルで保持するX軸の最大値。")
p.add_argument("--shift-film0", type=float, default=0.0, help="フィルムスペクトルの初期エネルギーシフト。")
p.add_argument("--shift-sub0", type=float, default=0.0, help="基板スペクトルの初期エネルギーシフト。")
p.add_argument("--shift-limit", type=float, default=5.0, help="各エネルギーシフトの絶対限界 [eV]。")
p.add_argument("--bg-order", type=int, default=1, choices=[0, 1, 2], help="多項式背景の次数。")
p.add_argument("--outside-fill", type=float, default=0.0, help="参照範囲外の補間値。")
p.add_argument("--amp-film0", type=float, default=1.0, help="フィルム振幅の初期値。")
p.add_argument("--amp-sub0", type=float, default=0.1, help="基板振幅の初期値。")
p.add_argument("--amp-min", type=float, default=1e-12, help="最小正振幅; ペナルティによって強制されます。")
p.add_argument("--amp-max", type=float, default=1e12, help="最大正振幅; ペナルティによって強制されます。")
p.add_argument("--bg0", type=float, default=0.0, help="初期の定数背景。")
p.add_argument("--bg1", type=float, default=0.0, help="初期の線形背景係数。")
p.add_argument("--bg2", type=float, default=0.0, help="初期の二次背景係数。")
p.add_argument("--loss", default="soft_l1", choices=["linear", "soft_l1", "logcosh"], help="Nelder-Meadの目的関数の損失タイプ。")
p.add_argument("--f-scale", type=float, default=0.05, help="堅牢な損失のスケーリング係数。")
p.add_argument("--loss-on-log", type=int, default=0, choices=[0, 1], help="線形残差ではなく、ログ強度比をフィッティングするかどうか。")
p.add_argument("--max-iter", type=int, default=5000, help="Nelder-Meadの最大反復回数。")
p.add_argument("--max-nfev", type=int, default=10000, help="最大関数評価回数。")
p.add_argument("--xatol", type=float, default=1e-8, help="Nelder-MeadのX許容誤差。")
p.add_argument("--fatol", type=float, default=1e-10, help="Nelder-Meadの目的関数許容誤差。")
p.add_argument("--penalty", type=float, default=1.0e6, help="ソフトバウンド制約のペナルティ強度。")
p.add_argument("--print-every", type=int, default=10, help="Nelder-Meadの収束情報をN回のコールバックごとに表示; 0で無効。")
p.add_argument("--norm-mode", default="max", choices=["max", "minmax", "area"], help="入力比較プロットの正規化モード。")
p.add_argument("--vbm-film-ref", type=float, default=None, help="同じX軸上のフィルム参照のVBM。")
p.add_argument("--vbm-sub-ref", type=float, default=None, help="同じX軸上の基板参照のVBM。")
p.add_argument("--output-stem", default="xps_vbo", help="出力ファイルのステム名。")
p.add_argument("--save", type=int, default=1, choices=[0, 1], help="PNGプロットを保存するかどうか。")
p.add_argument("--show", type=int, default=1, choices=[0, 1], help="プロットウィンドウを表示するかどうか。")
p.add_argument("--verbose", type=int, default=0, choices=[0, 1], help="scipyオプティマイザの詳細出力。")
return p
[ドキュメント]
def main():
"""
概要: スクリプトのメイン実行関数。
詳細説明:
コマンドライン引数を解析し、必要なスペクトルデータを読み込み、
`fit_mixture` 関数を呼び出してフィッティングを実行します。
その後、最適化されたパラメータと計算されたVBM/VBO情報をコンソールに出力し、
`save_outputs` 関数を呼び出して結果をファイルに保存します。
:param: なし
:returns: なし
"""
parser = build_parser()
args = parser.parse_args()
# argparse gives sheet as string; convert simple integer strings to int.
try:
args.sheet = int(args.sheet)
except Exception:
pass
film = crop(read_spectrum(args.film, args.sheet, args.x_col, args.y_col), args.ref_xmin, args.ref_xmax)
sub = crop(read_spectrum(args.substrate, args.sheet, args.x_col, args.y_col), args.ref_xmin, args.ref_xmax)
mix = read_spectrum(args.mix, args.sheet, args.x_col, args.y_col)
print(f"film : {film.name}, n={len(film.x)}, x=[{film.x.min():g}, {film.x.max():g}]")
print(f"substrate: {sub.name}, n={len(sub.x)}, x=[{sub.x.min():g}, {sub.x.max():g}]")
print(f"mix : {mix.name}, n={len(mix.x)}, x=[{mix.x.min():g}, {mix.x.max():g}]")
print("Optimizer: scipy.optimize.minimize(method='Nelder-Mead')")
fit = fit_mixture(film, sub, mix, args)
p = fit["params"]
print("\nOptimized parameters")
print(f" dE_film = {p['dE_film']:.8g} eV")
print(f" dE_substrate = {p['dE_sub']:.8g} eV")
print(f" dE_film-dE_sub = {p['dE_shift_film_minus_sub']:.8g} eV")
print(f" A_film = {p['A_film']:.8g}")
print(f" A_substrate = {p['A_sub']:.8g}")
print(f" RMSE = {p['rmse']:.8g}")
print(f" NRMSE = {p['nrmse']:.8g}")
print(f" objective = {p['objective']:.8g}")
print(f" success = {p['success']}")
print(f" message = {p['message']}")
if args.vbm_film_ref is not None and args.vbm_sub_ref is not None:
vbm_f = args.vbm_film_ref + p["dE_film"]
vbm_s = args.vbm_sub_ref + p["dE_sub"]
vbo_be = vbm_f - vbm_s
print("\nVBM / offset")
print(f" VBM_film_fit_BE = {vbm_f:.8g} eV")
print(f" VBM_sub_fit_BE = {vbm_s:.8g} eV")
print(f" VBO_BE = {vbo_be:.8g} eV (= film BE - substrate BE)")
print(f" opposite sign = {-vbo_be:.8g} eV")
else:
print("\nNote: supply --vbm-film-ref and --vbm-sub-ref to convert shift difference to VBO.")
csv_path, summary_path, png_path, norm_png_path = save_outputs(film, sub, mix, fit, args)
print("\nOutput")
print(f" fit csv : {csv_path}")
print(f" summary : {summary_path}")
if args.save:
print(f" fit plot : {png_path}")
print(f" normalized plot: {norm_png_path}")
if __name__ == "__main__":
main()