# Ne_T_fit.py ドキュメント
このドキュメントは、Pythonプログラム ``Ne_T_fit.py`` のソースコードを解析し、Sphinx（MyST）でビルド可能なMarkdown形式で記述されたものです。
## 概要
``Ne_T_fit.py`` は、半導体内の電子・正孔濃度、イオン化不純物濃度を温度の関数として計算するシミュレーターであり、実験的なキャリア密度データにモデルをフィッティングする機能も提供します。
フィッティングはデフォルトでネルダー・ミード法を使用し、ドナー/アクセプターのイオン化エネルギー、濃度、バンドギャップをフィッティングパラメータとして利用します。
ログ空間での最適化と、パラメータの不確かさ（共分散、相関、標準誤差）推定もサポートしています。
## 非標準ライブラリ
``Ne_T_fit.py`` が使用する非標準ライブラリは以下の通りです。
- ``numpy``: 数値計算
- ``pandas``: データフレームの操作、データ読み込み (CSV, Excel)
- ``argparse``: コマンドライン引数の解析
- ``matplotlib.pyplot``: データのプロット
- ``scipy``: 数値積分 (``integrate``), 最適化 (``optimize``), 特殊関数 (``special``), 物理定数 (``constants``)
## 主要な機能
### グローバル定数とフラグ
- ``KbT_const``: ボルツマン定数 ``k`` と素電荷 ``e`` の比 ``k / e`` を計算します。単位は ``[eV/K]`` です。
- ``FD_CACHE_DISABLE``: フェルミ・ディラック積分 ``F_half()`` のキャッシュ機構を制御するブールフラグです。``True`` に設定するとキャッシュが無効になり、より高精度な計算が可能になりますが、計算速度は低下します。数値ヤコビアンの計算時に一時的に ``True`` に設定されます。
### コア物理計算
#### ``_F_half_cached(eta_rounded_milli: int) -> float``
フェルミ・ディラック積分 $F_{1/2}(\eta)$ の結果をキャッシュする内部関数です。
高速化のため、$\eta$ を1000倍して整数に丸めた値をキャッシュキーとして使用します。
- **パラメータ**:
- ``eta_rounded_milli``: ``eta`` を1000倍して整数に丸めた値。
- **戻り値**:
- フェルミ・ディラック積分 $F_{1/2}(\eta)$ の結果。
#### ``F_half(eta: float) -> float``
係数 $2/\sqrt{\pi}$ を含まない1/2次のフェルミ・ディラック積分を計算します。
デフォルトでは高速化のために $\eta$ を $10^{-3}$ の精度で量子化し、結果をキャッシュします。
``FD_CACHE_DISABLE`` が ``True`` の場合は量子化を行わず、より正確な計算（ただし低速）を行います。これは特に数値ヤコビアンの計算時に重要となります。
- **パラメータ**:
- ``eta``: 無次元化されたフェルミ準位 $\eta = (E_f - E_c) / k_B T$ または $(E_v - E_f) / k_B T$。
- **戻り値**:
- フェルミ・ディラック積分 $F_{1/2}(\eta)$ の結果。
#### ``m2Nc(T: float, me_eff: float) -> float``
伝導帯の有効状態密度 $N_c$ ``[cm^-3]`` を計算します。
電子の有効質量 ``me_eff`` と温度 ``T`` を用いて有効状態密度を計算し、単位を ``m^-3`` から ``cm^-3`` に変換します。
- **パラメータ**:
- ``T``: 温度 ``[K]``。
- ``me_eff``: 電子の有効質量比 ($m_e^* / m_0$)。
- **戻り値**:
- 伝導帯の有効状態密度 $N_c$ ``[cm^-3]``。
#### ``m2Nv(T: float, mh_eff: float) -> float``
価電子帯の有効状態密度 $N_v$ ``[cm^-3]`` を計算します。
正孔の有効質量 ``mh_eff`` と温度 ``T`` を用いて有効状態密度を計算し、単位を ``m^-3`` から ``cm^-3`` に変換します。
- **パラメータ**:
- ``T``: 温度 ``[K]``。
- ``mh_eff``: 正孔の有効質量比 ($m_h^* / m_0$)。
- **戻り値**:
- 価電子帯の有効状態密度 $N_v$ ``[cm^-3]``。
#### ``Ne(Ef: float, Ec: float, T: float, Nc: float) -> float``
伝導帯における電子濃度 $n$ ``[cm^-3]`` を計算します。
フェルミ準位 ``Ef`` と伝導帯端 ``Ec`` の相対位置から無次元化されたフェルミ準位 $\eta$ を計算し、フェルミ・ディラック積分を用いて電子濃度を算出します。
- **パラメータ**:
- ``Ef``: フェルミ準位 ``[eV]``。
- ``Ec``: 伝導帯端のエネルギー ``[eV]``。
- ``T``: 温度 ``[K]``。
- ``Nc``: 伝導帯の有効状態密度 ``[cm^-3]``。
- **戻り値**:
- 電子濃度 $n$ ``[cm^-3]``。
#### ``Nh(Ef: float, Ev: float, T: float, Nv: float) -> float``
価電子帯における正孔濃度 $p$ ``[cm^-3]`` を計算します。
フェルミ準位 ``Ef`` と価電子帯端 ``Ev`` の相対位置から無次元化されたフェルミ準位 $\eta$ を計算し、フェルミ・ディラック積分を用いて正孔濃度を算出します。
- **パラメータ**:
- ``Ef``: フェルミ準位 ``[eV]``。
- ``Ev``: 価電子帯端のエネルギー ``[eV]``。
- ``T``: 温度 ``[K]``。
- ``Nv``: 価電子帯の有効状態密度 ``[cm^-3]``。
- **戻り値**:
- 正孔濃度 $p$ ``[cm^-3]``。
#### ``NDp(Ef: float, Ec: float, ED: float, ND: float, T: float, g: float = 1.0) -> float``
イオン化されたドナー密度 $N_D^+$ ``[cm^-3]`` を計算します。
ドナー準位 ``ED`` は伝導帯端 ``Ec`` から ``ED`` ``[eV]`` 下にあります (``ED > 0`` でより深い)。
ドナーの中性占有率 $f_D$ を用い、$N_D^+ = N_D (1 - f_D)$ でイオン化ドナー密度を算出します。
``special.expit`` を使用してオーバーフローを回避します。
- **パラメータ**:
- ``Ef``: フェルミ準位 ``[eV]``。
- ``Ec``: 伝導帯端のエネルギー ``[eV]``。
- ``ED``: ドナーのイオン化エネルギー（伝導帯端からの深さ） ``[eV]``。
- ``ND``: ドナーの全密度 ``[cm^-3]``。
- ``T``: 温度 ``[K]``。
- ``g``: ドナーの縮退因子 (デフォルト: ``1.0``)。
- **戻り値**:
- イオン化されたドナー密度 $N_D^+$ ``[cm^-3]``。
#### ``NAm(Ef: float, Ev: float, EA: float, NA: float, T: float, g: float = 1.0) -> float``
イオン化されたアクセプター密度 $N_A^-$ ``[cm^-3]`` を計算します。
アクセプター準位 ``EA`` は価電子帯端 ``Ev`` から ``EA`` ``[eV]`` 上にあります (``EA > 0`` でより深い)。
アクセプターのイオン化された状態の占有率 $f_A$ を用い、$N_A^- = N_A f_A$ でイオン化アクセプター密度を算出します。
``special.expit`` を使用してオーバーフローを回避します。
- **パラメータ**:
- ``Ef``: フェルミ準位 ``[eV]``。
- ``Ev``: 価電子帯端のエネルギー ``[eV]``。
- ``EA``: アクセプターのイオン化エネルギー（価電子帯端からの深さ） ``[eV]``。
- ``NA``: アクセプターの全密度 ``[cm^-3]``。
- ``T``: 温度 ``[K]``。
- ``g``: アクセプターの縮退因子 (デフォルト: ``1.0``)。
- **戻り値**:
- イオン化されたアクセプター密度 $N_A^-$ ``[cm^-3]``。
#### ``deltaQ(Ef: float, Ec: float, Ev: float, ED: float, EA: float, ND: float, NA: float, T: float, Nc: float, Nv: float, gD: float = 1.0, gA: float = 1.0) -> float``
電荷中性条件の残差 $p - n + N_D^+ - N_A^-$ を計算します。
フェルミ準位 ``Ef`` を解くための根を探索する関数です。
使用する縮退因子 ``gD`` / ``gA`` は、``Ef`` の計算時と $N_D^+$ / $N_A^-$ の事後評価時で一貫している必要があります。
- **パラメータ**:
- ``Ef``: フェルミ準位 ``[eV]``。
- ``Ec``: 伝導帯端のエネルギー ``[eV]``。
- ``Ev``: 価電子帯端のエネルギー ``[eV]``。
- ``ED``: ドナーのイオン化エネルギー ``[eV]``。
- ``EA``: アクセプターのイオン化エネルギー ``[eV]``。
- ``ND``: ドナーの全密度 ``[cm^-3]``。
- ``NA``: アクセプターの全密度 ``[cm^-3]``。
- ``T``: 温度 ``[K]``。
- ``Nc``: 伝導帯の有効状態密度 ``[cm^-3]``。
- ``Nv``: 価電子帯の有効状態密度 ``[cm^-3]``。
- ``gD``: ドナーの縮退因子 (デフォルト: ``1.0``)。
- ``gA``: アクセプターの縮退因子 (デフォルト: ``1.0``)。
- **戻り値**:
- 電荷中性条件の残差。
#### ``solve_Ef(Ec: float, Ev: float, ED: float, EA: float, ND: float, NA: float, T: float, Nc: float, Nv: float, gD: float = 1.0, gA: float = 1.0, low0: float = -5.0, high0: float = 5.0, max_expand: int = 50) -> float``
Brentq法と自動ブラケット拡張を用いてフェルミ準位 ``Ef`` を解きます。
``deltaQ()`` 関数のゼロ点を見つけることで、``Ef`` を計算します。
初期の探索範囲 (``low0``, ``high0``) で根が見つからない場合、自動的に探索範囲を拡張します。
- **パラメータ**:
- ``Ec``: 伝導帯端のエネルギー ``[eV]``。
- ``Ev``: 価電子帯端のエネルギー ``[eV]``。
- ``ED``: ドナーのイオン化エネルギー ``[eV]``。
- ``EA``: アクセプターのイオン化エネルギー ``[eV]``。
- ``ND``: ドナーの全密度 ``[cm^-3]``。
- ``NA``: アクセプターの全密度 ``[cm^-3]``。
- ``T``: 温度 ``[K]``。
- ``Nc``: 伝導帯の有効状態密度 ``[cm^-3]``。
- ``Nv``: 価電子帯の有効状態密度 ``[cm^-3]``。
- ``gD``: ドナーの縮退因子 (デフォルト: ``1.0``)。
- ``gA``: アクセプターの縮退因子 (デフォルト: ``1.0``)。
- ``low0``: ``Ef`` 探索の初期下限値 ``[eV]``。
- ``high0``: ``Ef`` 探索の初期上限値 ``[eV]``。
- ``max_expand``: 探索範囲を拡張する最大回数。
- **戻り値**:
- 計算されたフェルミ準位 ``Ef`` ``[eV]``。
#### `class Params`
半導体の物理パラメータを格納するデータクラスです。
バンドギャップ、ドナー/アクセプターのイオン化エネルギー、濃度、有効質量といった主要な物理パラメータをカプセル化します。
- **属性**:
- ``Eg``: バンドギャップ ``[eV]``。
- ``ED``: ドナーのイオン化エネルギー ``[eV]``。
- ``EA``: アクセプターのイオン化エネルギー ``[eV]``。
- ``ND``: ドナー濃度 ``[cm^-3]``。
- ``NA``: アクセプター濃度 ``[cm^-3]``。
- ``me``: 電子の有効質量比 ($m_e^* / m_0$)。
- ``mh``: 正孔の有効質量比 ($m_h^* / m_0$)。
#### ``calc_at_T(T: float, p: Params, args=None) -> Dict[str, float]``
指定された温度 ``T`` と物理パラメータ ``p`` におけるキャリア統計関連の値を計算します。
有効状態密度 ($N_c$, $N_v$) を計算し、電荷中性条件を満たすフェルミ準位 (``Ef``) を求め、そこから電子・正孔濃度、イオン化ドナー・アクセプター濃度を算出します。
結果は辞書形式で返されます。
- **パラメータ**:
- ``T``: 温度 ``[K]``。
- ``p``: 半導体の物理パラメータ (``Params`` オブジェクト)。
- ``args``: コマンドライン引数を格納するオブジェクト（ドナー/アクセプターの縮退因子 ``gD`` / ``gA`` を取得するために使用）。
- **戻り値**:
- 計算結果を格納した辞書。キーは ``"T"``, ``"Ef"``, ``"n"``, ``"p"``, ``"NDp"``, ``"NAm"``, ``"Nc"``, ``"Nv"`` です。
### フィッティングユーティリティ
#### ``FIT_NAMES``
フィッティングの対象となるパラメータ名がリストされた定数です。
``["ED", "ND", "EA", "NA", "Eg"]``。
#### ``parse_fix_list(s: str | None) -> set[str]``
カンマ区切りの文字列から固定パラメータ名のセットをパースします。
``--fix`` オプションで与えられた文字列を処理し、空白を除去して有効なパラメータ名のセットを生成します。
- **パラメータ**:
- ``s``: 固定パラメータ名のカンマ区切り文字列、または ``None``。
- **戻り値**:
- 固定パラメータ名を含む文字列のセット。
#### ``sanitize_args_for_fit(args, fix: set[str]) -> None``
フィッティングのためにCLI引数のパラメータ値をサニタイズします。
- ``ND`` と ``NA`` が自由変数として最適化される場合、``log10`` 変換を避けるために0より大きい値に強制します。
- ``Eg``, ``ED``, ``EA`` については物理的に妥当な範囲に調整します。
- ``gD``, ``gA`` は ``log()`` で使用されるため0より大きい値に強制します。
- **パラメータ**:
- ``args``: CLI引数を格納するオブジェクト。
- ``fix``: 固定パラメータ名のセット。
#### ``read_data(path: str, temp_col: int = 0, ne_col: int = 3, sheet_name: str | int | None = None) -> Tuple[np.ndarray, np.ndarray, pd.DataFrame]``
CSVまたはExcelファイルから実験データを読み込みます。
指定されたパスからファイル（CSV, TSV, Excel）を読み込み、温度 (``temp_col``) とキャリア濃度 (``ne_col``) の列を抽出します。
NaN値は除去されます。Excelファイルの場合、シート名を指定できます。
- **パラメータ**:
- ``path``: 入力ファイルのパス。
- ``temp_col``: 温度データの列インデックス (0始まり) (デフォルト: ``0``)。
- ``ne_col``: キャリア濃度データの列インデックス (0始まり) (デフォルト: ``3``)。
- ``sheet_name``: Excelファイルの場合のシート名またはシートインデックス (デフォルト: ``None``, 最初のシート)。
- **戻り値**:
- 温度のNumpy配列、キャリア濃度のNumpy配列、読み込んだDataFrame。
#### ``save_params(params: dict, filename: str = "fit_params.json")``
パラメータ辞書をJSONファイルとして保存します。
フィッティング結果などのパラメータを、指定されたファイル名でJSON形式で保存します。
- **パラメータ**:
- ``params``: 保存するパラメータを格納した辞書。
- ``filename``: 保存するJSONファイル名 (デフォルト: ``"fit_params.json"``)。
#### ``load_params(filename: str = "fit_params.json") -> dict``
JSONファイルからパラメータ辞書を読み込みます。
指定されたJSONファイルからパラメータを読み込み、辞書として返します。
ファイルが存在しない場合は空の辞書を返します。
- **パラメータ**:
- ``filename``: 読み込むJSONファイル名 (デフォルト: ``"fit_params.json"``)。
- **戻り値**:
- 読み込んだパラメータを格納した辞書。
#### ``_cli_overrides(argv: List[str]) -> set[str]``
コマンドラインで明示的に指定された引数の宛先名（``argparse`` の ``dest`` 名）のセットを返します。
``sys.argv`` をパースし、``--opt value`` または ``--opt=value`` 形式で提供されたオプションに対応する ``argparse`` の宛先名を特定します。これにより、設定ファイルからのロード時にCLIで指定された値が上書きされないようにします。
- **パラメータ**:
- ``argv``: コマンドライン引数のリスト（通常は ``sys.argv[1:]``）。
- **戻り値**:
- コマンドラインで指定された引数の宛先名を含むセット。
#### ``_apply_loaded_params_to_args(args, params: dict, cli_overrides: set[str]) -> None``
読み込まれたパラメータを ``args`` オブジェクトに適用します。
パラメータ辞書 ``params`` の値を ``args`` オブジェクトの対応する属性に設定します。
ただし、``cli_overrides`` に含まれる、つまりCLIで明示的に指定されたパラメータは上書きしません。
既知の数値パラメータのみを対象とします。
- **パラメータ**:
- ``args``: ``argparse.Namespace`` オブジェクト。
- ``params``: 読み込まれたパラメータを格納した辞書。
- ``cli_overrides``: CLIで明示的に指定されたパラメータ名のセット。
#### ``pack_free_params(args, fix: set[str]) -> Tuple[np.ndarray, List[str]]``
フィッティングのために、自由パラメータを内部表現のベクトル ``theta0`` にパックします。
固定されていないパラメータのみを抽出し、``ND``, ``NA`` は ``log10`` 空間に変換してベクトルに格納します。
エネルギー系のパラメータは線形空間のままです。
- **パラメータ**:
- ``args``: CLI引数を格納するオブジェクト。
- ``fix``: 固定パラメータ名のセット。
- **戻り値**:
- 自由パラメータの初期値ベクトル (``numpy.ndarray``) とその名前のリスト。
#### ``unpack_params(theta: np.ndarray, args, fix: set[str]) -> Params``
内部表現のパラメータベクトル ``theta`` と固定値を物理的な ``Params`` オブジェクトにアンパックします。
最適化された自由パラメータ ``theta`` の値と、CLI引数 ``args`` から取得した固定パラメータの値を結合し、物理的な意味を持つ ``Params`` オブジェクトを構築します。
``log10`` 変換された ``ND``, ``NA`` は元のスケールに戻されます。
- **パラメータ**:
- ``theta``: 自由パラメータの内部表現ベクトル (``numpy.ndarray``)。
- ``args``: CLI引数を格納するオブジェクト（固定パラメータの値、および ``me``, ``mh`` を取得）。
- ``fix``: 固定パラメータ名のセット。
- **戻り値**:
- 物理パラメータを格納した ``Params`` オブジェクト。
#### ``residuals_log10(T: np.ndarray, Ne_obs: np.ndarray, p: Params, args=None) -> np.ndarray``
モデルと観測値の間の残差ベクトルを ``log10`` 空間で計算します。
モデル計算で得られた電子濃度 ``n_model`` と観測された電子濃度 ``Ne_obs`` の比の ``log10`` を残差として返します。
これにより相対誤差が強調され、広い濃度範囲でのフィッティングに適します。
``Ne_obs`` は0より大きい必要があります。
- **パラメータ**:
- ``T``: 温度のNumpy配列 ``[K]``。
- ``Ne_obs``: 観測された電子濃度データのNumpy配列 ``[cm^-3]``。
- ``p``: 物理パラメータ (``Params`` オブジェクト)。
- ``args``: コマンドライン引数を格納するオブジェクト（縮退因子を取得するため）。
- **戻り値**:
- ``log10`` 空間での残差ベクトル (``numpy.ndarray``)。
#### ``objective(theta: np.ndarray, T: np.ndarray, Ne_obs: np.ndarray, args, fix: set[str]) -> float``
最適化アルゴリズムのための目的関数です。``log10`` 空間での残差の二乗和を返します。
自由パラメータ ``theta`` から物理パラメータを再構築し、``residuals_log10()`` を計算し、その二乗和を目的値とします。物理的に不可能なパラメータ値（例: バンドギャップが負）や不自然な関係（例: ドナー準位がバンドギャップより深い）に対してはペナルティ項を加えて最適化をガイドします。
- **パラメータ**:
- ``theta``: 自由パラメータの内部表現ベクトル (``numpy.ndarray``)。
- ``T``: 温度のNumpy配列 ``[K]``。
- ``Ne_obs``: 観測された電子濃度データのNumpy配列 ``[cm^-3]``。
- ``args``: コマンドライン引数を格納するオブジェクト。
- ``fix``: 固定パラメータ名のセット。
- **戻り値**:
- ``log10`` 空間での残差の二乗和（目的関数の値）。
#### ``numerical_jacobian(theta: np.ndarray, T: np.ndarray, Ne_obs: np.ndarray, args, fix: set[str]) -> np.ndarray``
残差ベクトルに関するパラメータ ``theta`` の数値ヤコビアンを計算します。中央差分法を使用します。
**重要**:
- このコードベースでは $F_{1/2}(\eta)$ の計算を高速化するために $\eta$ を $10^{-3}$ 単位で量子化してキャッシュします。
- 有限差分によるヤコビアン計算の場合、この量子化は小さなパラメータ摂動に対して残差を区分的に定数とし、``J≈0``、したがって ``(J^T J)≈0`` を引き起こす可能性があります。
- したがって、ここでは一時的に量子化キャッシュ (``FD_CACHE_DISABLE=True``) を無効にします。
- **パラメータ**:
- ``theta``: 自由パラメータの内部表現ベクトル (``numpy.ndarray``)。
- ``T``: 温度のNumpy配列 ``[K]``。
- ``Ne_obs``: 観測された電子濃度データのNumpy配列 ``[cm^-3]``。
- ``args``: コマンドライン引数を格納するオブジェクト。
- ``fix``: 固定パラメータ名のセット。
- **戻り値**:
- ヤコビアン行列 (``numpy.ndarray``)。
#### ``estimate_param_errors(theta_hat: np.ndarray, T: np.ndarray, Ne_obs: np.ndarray, args, fix: set[str], free_names: List[str]) -> Dict[str, Dict[str, float]]``
最適化されたパラメータの1-シグマ誤差を推定します。
フィッティング後のパラメータ ``theta_hat`` における残差と数値ヤコビアンを用いて、自由パラメータの共分散行列を推定します。共分散行列から標準誤差（1-シグマ誤差）を導出し、``log10`` 変換されたパラメータ (``ND``, ``NA``) については物理空間での誤差を伝播させて算出します。
診断目的で相関行列なども含めて返します。
- **パラメータ**:
- ``theta_hat``: 最適化された自由パラメータの内部表現ベクトル (``numpy.ndarray``)。
- ``T``: フィッティングに使用された温度のNumpy配列 ``[K]``。
- ``Ne_obs``: フィッティングに使用された観測電子濃度データのNumpy配列 ``[cm^-3]``。
- ``args``: コマンドライン引数を格納するオブジェクト。
- ``fix``: 固定パラメータ名のセット。
- ``free_names``: 自由パラメータ名のリスト。
- **戻り値**:
- 各パラメータの推定値、標準誤差、および共分散・相関行列などの診断情報を含む辞書。
### 共分散/相関固有値診断と提案
#### ``cov_to_corr(Cov: np.ndarray) -> Tuple[np.ndarray, np.ndarray]``
共分散行列から相関行列と標準偏差を計算します。
共分散行列 ``Cov`` の対角要素から標準偏差を抽出し、それらを用いて相関行列 ``Corr`` を計算します。
対角要素は ``1.0``、標準偏差が ``0`` の場合は対応する相関値を ``NaN`` とします。
- **パラメータ**:
- ``Cov``: 共分散行列 (``numpy.ndarray``)。
- **戻り値**:
- 相関行列 (``numpy.ndarray``) と標準偏差の配列 (``numpy.ndarray``)。
#### ``eigen_sorted_sym(A: np.ndarray, descending: bool = True) -> Tuple[np.ndarray, np.ndarray]``
対称行列の固有値分解を行い、結果をソートして返します。
対称行列 ``A`` の固有値と固有ベクトルを計算し、固有値を降順（``descending=True`` の場合）または昇順にソートして返します。
固有ベクトルは対応する固有値の順序に合わせて並べ替えられます。
- **パラメータ**:
- ``A``: 対称行列 (``numpy.ndarray``)。
- ``descending``: 固有値を降順にソートするかどうか (デフォルト: ``True``)。
- **戻り値**:
- ソートされた固有値の配列 (``numpy.ndarray``) と対応する固有ベクトルの行列 (``numpy.ndarray``)。
#### ``summarize_eigenvectors(evals: np.ndarray, evecs: np.ndarray, names: List[str], topk: int = 3, compk: int = 4)``
主要な固有ベクトルを人間が読める形式で要約します。
固有値と固有ベクトル、パラメータ名のリストを受け取り、上位 ``topk`` 個の固有ベクトルについて、その主要な ``compk`` 個の成分と対応するパラメータ名を報告する辞書のリストを生成します。
これにより、どのパラメータが特定の不確かさの方向に強く寄与しているかを理解しやすくなります。
- **パラメータ**:
- ``evals``: 固有値の配列。
- ``evecs``: 固有ベクトルの行列。
- ``names``: パラメータ名のリスト。
- ``topk``: 要約する上位の固有ベクトルの数 (デフォルト: ``3``)。
- ``compk``: 各固有ベクトルで報告する主要成分の数 (デフォルト: ``4``)。
- **戻り値**:
- 要約された固有ベクトル情報を含む辞書のリスト。
#### ``propose_fix_candidates(names: List[str], theta_hat: np.ndarray, stderr: np.ndarray, corr: np.ndarray, evals_cov: np.ndarray, evecs_cov: np.ndarray, corr_thr: float = 0.95, relerr_thr: float = 0.5, topn: int = 3)``
フィッティングパラメータとして固定すべき/制約すべき候補をヒューリスティックに提案します。
- 自由パラメータの相対標準誤差が閾値 ``relerr_thr`` を超えるもの。
- 他のパラメータとの相関が閾値 ``corr_thr`` を超えるもの。
- 共分散行列の最も不確かな固有方向（最大固有値に対応する固有ベクトル）に強く寄与するものを特定し、提案リストを生成します。
- 各候補について、スコア、推定値、標準誤差、相対誤差、および提案理由を報告します。
- **パラメータ**:
- ``names``: 自由パラメータ名のリスト。
- ``theta_hat``: 最適化された自由パラメータの内部表現ベクトル (``numpy.ndarray``)。
- ``stderr``: 自由パラメータの標準誤差の配列 (``numpy.ndarray``)。
- ``corr``: 自由パラメータの相関行列 (``numpy.ndarray``)。
- ``evals_cov``: 共分散行列の固有値の配列。
- ``evecs_cov``: 共分散行列の固有ベクトルの行列。
- ``corr_thr``: 強相関と見なす相関係数の閾値 (デフォルト: ``0.95``)。
- ``relerr_thr``: 相対誤差が高いと見なす閾値 (デフォルト: ``0.5``)。
- ``topn``: 上位何件の候補を提案するか (デフォルト: ``3``)。
- **戻り値**:
- 固定/制約候補の提案を含む辞書のリスト。
#### ``prediction_band_log10_Ne(T: np.ndarray, theta_hat: np.ndarray, args, fix: set[str], Cov: np.ndarray, nsigma: float = 1.0) -> Tuple[np.ndarray, np.ndarray, np.ndarray]``
``log10(Ne)`` における予測誤差帯（デルタ法）を計算します。
最適化されたパラメータ ``theta_hat`` とその共分散行列 ``Cov`` を用いて、デルタ法により予測値 ``y0`` (``log10(Ne)``) の標準誤差を計算し、指定された ``nsigma`` 幅の誤差帯 (``ylo``, ``yhi``) を提供します。
数値ヤコビアンの計算と同様に ``F_half()`` 関数のキャッシュを一時的に無効にします。
- **パラメータ**:
- ``T``: 温度のNumpy配列 ``[K]``。
- ``theta_hat``: 最適化された自由パラメータの内部表現ベクトル (``numpy.ndarray``)。
- ``args``: コマンドライン引数を格納するオブジェクト。
- ``fix``: 固定パラメータ名のセット。
- ``Cov``: 自由パラメータの共分散行列 (``numpy.ndarray``)。
- ``nsigma``: 誤差帯の幅（標準偏差の倍数） (デフォルト: ``1.0``)。
- **戻り値**:
- 予測中央値 (``y0``)、下限 (``ylo``)、上限 (``yhi``) のNumpy配列 (``log10`` 空間)。共分散行列が利用できない場合は ``None`` を返します。
### CLI / main 関数
#### ``initialize()``
コマンドライン引数をパースし、``argparse.Namespace`` オブジェクトとして返します。
シミュレーションおよびフィッティングに必要な各種パラメータ（I/O、物理定数、温度範囲、フィッティングオプション、診断出力など）を定義し、ユーザーからの入力を処理します。
- **戻り値**:
- パースされたコマンドライン引数を格納した ``argparse.Namespace`` オブジェクト。
#### ``run_read(args)``
入力ファイルからデータを読み込み、プロットして表示します。
``read_data()`` 関数を使用して指定された入力ファイルから温度と電子濃度データを読み込み、``matplotlib`` を使用して ``Ne(T)`` データをプロットし、表示します。
- **パラメータ**:
- ``args``: コマンドライン引数を格納する ``argparse.Namespace`` オブジェクト。
#### ``run_sim(args)``
指定されたパラメータでキャリア統計をシミュレーションし、結果をプロットします。
- オプションでJSONファイルから初期パラメータをロードします（CLI指定が優先）。
- 入力データファイルが指定されている場合は、その温度範囲を使用してシミュレーションを行います。そうでない場合は、``--Tmin``, ``--Tmax``, ``--nT`` で定義された温度範囲を使用します。
- ``calc_at_T()`` 関数で電子濃度を計算し、結果をデータ（あれば）と共に ``1000/T`` vs ``Ne`` の半対数プロットとして表示します。
- **パラメータ**:
- ``args``: コマンドライン引数を格納する ``argparse.Namespace`` オブジェクト。
#### ``run_fit(args)``
実験データに対して半導体パラメータをフィッティングします。
- 入力ファイルから実験データを読み込み、指定された温度範囲でフィッティングを行います。
- JSONファイルから初期パラメータをロードします（CLI指定が優先）。
- 指定されたパラメータを固定し、残りの自由パラメータを ``nelder-mead`` などの最適化アルゴリズムでフィッティングします。``ND``, ``NA`` は ``log10`` 空間で最適化されます。
- フィット結果、パラメータの不確かさ（標準誤差、共分散、相関）、および固定すべきパラメータの提案をコンソールに出力し、オプションでJSONファイルに保存します。
- モデル計算値と実験データ、および予測誤差帯をプロットします。
- **パラメータ**:
- ``args``: コマンドライン引数を格納する ``argparse.Namespace`` オブジェクト。
#### ``replace_ext(path: str, target_ext: str) -> str``
ファイルパスの拡張子を新しいものに置き換えます。
与えられたファイルパスの拡張子を ``target_ext`` で指定されたものに置き換えた新しいファイルパス文字列を返します。
ファイル名自体は変更しません。
- **パラメータ**:
- ``path``: 元のファイルパス文字列。
- ``target_ext``: 新しい拡張子（例: ``".json"``）。
- **戻り値**:
- 拡張子が置き換えられた新しいファイルパス文字列。
#### ``main()``
プログラムのエントリポイントです。
``initialize()`` 関数を呼び出してコマンドライン引数をパースします。
``args.save``, ``args.diag_save``, ``args.suggest_save`` が指定されていない場合、``args.input`` に基づいてデフォルトのファイル名を生成します。
``--mode`` 引数の値に応じて、``run_read()``, ``run_fit()``, ``run_sim()`` のいずれかの関数を実行します。
## 入出力仕様
### コマンドライン引数
``Ne_T_fit.py`` は ``argparse`` モジュールを使用してコマンドライン引数を処理します。
- **入出力ファイル**:
- ``--input <FILE>``: 入力ファイル名。CSVまたはExcel形式をサポートします (デフォルト: ``Hall-T1.xlsx``)。
- ``--sheet <NAME_OR_INDEX>``: Excelファイルの場合のシート名またはシートインデックス (デフォルト: 最初のシート)。
- ``--temp_col <INT>``: 温度データが含まれる列のインデックス（0から始まる） (デフォルト: ``0``)。
- ``--ne_col <INT>``: キャリア濃度データが含まれる列のインデックス（0から始まる） (デフォルト: ``3``)。
- ``--save <FILE>``: フィッティング結果のパラメータを保存するJSONファイル名 (デフォルト: ``<INPUT_BASENAME>_params.json``)。
- ``--load <FILE>``: シミュレーションまたはフィッティングの初期パラメータを読み込むJSONファイル名 (デフォルト: 空文字列 ``""``)。
- **プログラムモード**:
- ``--mode {read,sim,fit}``: 実行モードを指定します。
- ``read``: 入力データを読み込み、プロットして表示します。
- ``sim``: 指定されたパラメータでキャリア濃度をシミュレーションし、プロットします。
- ``fit``: 実験データに対してパラメータをフィッティングし、結果をプロットします。
(デフォルト: ``read``)
- **物理パラメータ (シミュレーション初期値/フィット対象)**:
- ``--Eg <FLOAT>``: バンドギャップ ``[eV]`` (デフォルト: ``1.12``)。
- ``--me <FLOAT>``: 電子の有効質量比 ($m_e^* / m_0$) (デフォルト: ``1.08``)。
- ``--mh <FLOAT>``: 正孔の有効質量比 ($m_h^* / m_0$) (デフォルト: ``0.55``)。
- ``--ND <FLOAT>``: ドナー濃度 ``[cm^-3]`` (デフォルト: ``1e17``)。
- ``--NA <FLOAT>``: アクセプター濃度 ``[cm^-3]`` (デフォルト: ``1e15``)。
- ``--ED <FLOAT>``: ドナーのイオン化エネルギー（伝導帯端からの深さ） ``[eV]`` (デフォルト: ``0.045``)。
- ``--EA <FLOAT>``: アクセプターのイオン化エネルギー（価電子帯端からの深さ） ``[eV]`` (デフォルト: ``0.045``)。
- ``--gD <FLOAT>``: ドナーの縮退因子 (デフォルト: ``1.0``)。
- ``--gA <FLOAT>``: アクセプターの縮退因子 (デフォルト: ``1.0``)。
- **温度/濃度範囲 (シミュレーション/プロット)**:
- ``--Tmin <FLOAT>``: シミュレーションまたはプロットの最小温度 ``[K]`` (デフォルト: ``50``)。
- ``--Tmax <FLOAT>``: シミュレーションまたはプロットの最大温度 ``[K]`` (デフォルト: ``600``)。
- ``--nT <INT>``: シミュレーション時の温度点の数 (入力データがない場合) (デフォルト: ``50``)。
- ``--Nmin <FLOAT>``: プロットの最小キャリア濃度 ``[cm^-3]`` (デフォルト: ``1e10``)。
- ``--Nmax <FLOAT>``: プロットの最大キャリア濃度 ``[cm^-3]`` (デフォルト: ``1e22``)。
- **フィッティングオプション**:
- ``--method <STR>``: フィットモードで使用する最適化手法 (デフォルト: ``nelder-mead``)。
- ``--fix <STR>``: 固定するパラメータ名をカンマ区切りで指定します。対象は ``ED``, ``ND``, ``EA``, ``NA``, ``Eg`` です (デフォルト: ``Eg,NA,EA``)。
- ``--Tfitmin <FLOAT>``: フィットに使用する最小温度 ``[K]`` (デフォルト: ``-1e100``)。
- ``--Tfitmax <FLOAT>``: フィットに使用する最大温度 ``[K]`` (デフォルト: ``+1e100``)。
- **不確かさ計算/ヤコビアン制御**:
- ``--band_sigma <FLOAT>``: 誤差帯の幅（標準偏差の倍数） (デフォルト: ``1.0``)。
- ``--jac_relstep <FLOAT>``: 数値微分の相対ステップサイズ (デフォルト: ``1e-6``)。
- ``--jac_absstep <FLOAT>``: 数値微分の絶対ステップサイズ (デフォルト: ``1e-12``)。
- **診断出力**:
- ``--diag_save <FILE>``: 診断出力（共分散行列、固有値など）を保存するJSONファイル名 (デフォルト: ``<INPUT_BASENAME>_diagnostics_opt.json``)。
- ``--suggest_save <FILE>``: 固定候補の提案を保存するJSONファイル名 (デフォルト: ``<INPUT_BASENAME>_fix_suggestions.json``)。
### 入力データ形式
入力ファイルはCSV、TSV、またはExcel形式をサポートします。
データは少なくとも2つの数値列を含んでいる必要があり、``--temp_col`` と ``--ne_col`` で指定された列がそれぞれ温度とキャリア濃度を表します。
- **CSV/TSVの例**:
```text
Temperature, Electron_Density, ...
10.0, 1.2e10, ...
20.0, 5.5e11, ...
...
```
- **Excelの例**:
- シート名は ``--sheet`` オプションで指定可能です。指定しない場合は最初のシートが使用されます。
- ``temp_col`` および ``ne_col`` は0から始まる列インデックスです。
### 出力データ形式
- **標準出力**: シミュレーション結果、フィッティングの進行状況、フィットされたパラメータ、不確かさの推定、診断情報、固定候補の提案などがコンソールに出力されます。
- **プロット**: ``matplotlib`` を用いて、``Ne(T)`` (または ``Ne(1000/T)``) の半対数プロットがリアルタイムで表示されます。フィットモードでは、実験データ、モデル計算値、および予測誤差帯がプロットされます。
- **JSONファイル**:
- ``--save`` オプションで指定されたファイル (デフォルト: ``<INPUT_BASENAME>_params.json``) に、フィットされた物理パラメータ (``Eg``, ``ED``, ``EA``, ``ND``, ``NA``, ``me``, ``mh`` など) がJSON形式で保存されます。
- ``--diag_save`` オプションで指定されたファイル (デフォルト: ``<INPUT_BASENAME>_diagnostics_opt.json``) に、フィッティングの診断情報（共分散行列、相関行列、固有値分解結果など）がJSON形式で保存されます。
- ``--suggest_save`` オプションで指定されたファイル (デフォルト: ``<INPUT_BASENAME>_fix_suggestions.json``) に、固定すべきパラメータの提案がJSON形式で保存されます。
## 使用例
### 1. データの読み込みと表示 (``read`` モード)
入力ファイルを読み込み、``Ne(T)`` データをプロットして表示します。
```bash
python Ne_T_fit.py --mode read --input Hall_data.csv --temp_col 0 --ne_col 1
2. パラメータを指定してシミュレーション (sim モード)
指定された物理パラメータでキャリア濃度を計算し、Ne(1000/T) をプロットします。
--input が指定されている場合は、その温度範囲でプロットし、データがあればオーバーレイされます。
python Ne_T_fit.py --mode sim --Eg 1.12 --ND 1e16 --ED 0.045 --NA 1e14 --EA 0.1 --Tmin 100 --Tmax 500 --nT 100
JSONファイルから初期パラメータをロードしてシミュレーションを行うことも可能です。
python Ne_T_fit.py --mode sim --load my_initial_params.json --input Hall_data.csv
3. 実験データにフィッティング (fit モード)
実験データ Hall_data.csv にモデルをフィッティングし、フィット結果を fit_result.json に保存します。
Eg, NA, EA を固定し、ND, ED を自由パラメータとしてフィッティングする例です。
python Ne_T_fit.py --mode fit --input Hall_data.csv --temp_col 0 --ne_col 1 --fix Eg,NA,EA --save fit_result.json
特定の温度範囲のみをフィッティングに使用することも可能です。
python Ne_T_fit.py --mode fit --input Hall_data.csv --Tfitmin 150 --Tfitmax 400 --save fit_partial_T.json
フィッティングの診断情報と固定候補の提案をファイルに保存します。
python Ne_T_fit.py --mode fit --input Hall_data.csv --diag_save diagnostics.json --suggest_save suggestions.json