# `mu_fit.py` ドキュメント
このドキュメントは、Hall効果移動度解析・フィッティングツールである `mu_fit.py` スクリプトの機能と使用方法について説明します。
## 概要
`mu_fit.py` は、Hall効果測定から得られた温度依存移動度データに対して、複数の散乱メカニズム（光学フォノン、音響フォノン、中性不純物、イオン化不純物、粒界散乱）の物理モデルに基づいたフィッティングを行います。このツールは、パラメータの最適化、誤差評価、およびパラメータ感度診断の機能を提供します。
主な機能は以下の通りです。
* 実験データの読み込み（CSVまたはExcel形式）。
* 線形最小二乗法 (LLSQ) による初期値推定。
* 非線形最適化によるモデルパラメータの抽出。
* 誤差伝播法によるモデルの不確かさ（1σ範囲）の算出。
* 共分散行列、相関係数、固有値解析を用いたパラメータ感度診断。
:::note
フィッティングの目的関数は、移動度の対数スケールにおける残差二乗和を最小化します。
:::
## 要件とインストール
### Pythonバージョン
本スクリプトは Python 3 で動作します。
### 非標準ライブラリ
以下のライブラリが必要です。
* ``pandas``: データフレーム操作用
* ``numpy``: 数値計算用
* ``matplotlib``: グラフ描画用
* ``scipy``: 最適化、線形代数計算用
これらのライブラリは ``pip`` コマンドでインストールできます。
```bash
pip install pandas numpy matplotlib scipy
```
## 使用方法
``mu_fit.py`` はコマンドライン引数を受け取り、様々なモードで実行可能です。
### CLIオプション
スクリプトは ``argparse`` モジュールを使用してコマンドライン引数を解析します。
* ``--input FILE_PATH``
* 入力ファイル名を指定します。CSVまたはExcel形式のファイルを指定可能です。
* デフォルト値: ``Hall-T1.xlsx``
* ``--temp_col INT``
* 入力ファイル内の温度データが存在する列のインデックス（0開始）を指定します。
* デフォルト値: ``0``
* ``--mu_col INT``
* 入力ファイル内の移動度データが存在する列のインデックス（0開始）を指定します。
* デフォルト値: ``2``
* ``--mode {read,llsq,fit,weight}``
* 実行モードを指定します。
* ``read``: 実験データを読み込み、プロットのみを行います。
* ``llsq``: 線形最小二乗法で初期値を推定し、結果をプロットしてJSONファイルに保存します。
* ``fit``: 非線形最適化によるフィッティングを実行し、結果をプロットします。
* ``weight``: 保存されたパラメータに基づいて、各散乱メカニズムの寄与度を可視化します。
* デフォルト値: ``read``
* ``--method TEXT``
* 非線形最適化に使用する ``scipy.optimize.minimize`` のメソッドを指定します。
* デフォルト値: ``Nelder-Mead``
* ``--eop FLOAT``
* 光学フォノンエネルギーを ``eV`` 単位で指定します。
* デフォルト値: ``0.045``
* ``--fix PARAM_LIST``
* フィッティング時に固定するパラメータをカンマ区切りで指定します（例: ``a2,VB``）。
* デフォルト値: 空文字列（全パラメータを最適化）
* ``--Tfitmin FLOAT``
* フィッティングに使用するデータの最小温度を ``K`` 単位で指定します。
* デフォルト値: ``-1e100`` (実質的に制限なし)
* ``--Tfitmax FLOAT``
* フィッティングに使用するデータの最大温度を ``K`` 単位で指定します。
* デフォルト値: ``+1e100`` (実質的に制限なし)
* ``--band_sigma FLOAT``
* 予測誤差帯の幅をシグマ倍率で指定します。例えば ``1.0`` であれば ±1σ を意味します。
* デフォルト値: ``1.0``
* ``--jac_relstep FLOAT``
* 数値微分の相対ステップサイズを指定します。ヤコビ行列の計算に使用されます。
* デフォルト値: ``1e-6``
* ``--jac_absstep FLOAT``
* 数値微分の絶対ステップサイズを指定します。ヤコビ行列の計算に使用されます。
* デフォルト値: ``1e-12``
### 基本的な実行例
1. **データの表示のみ:**
```bash
python mu_fit.py --input data.csv --mode read
```
（``data.csv`` の代わりにデフォルトの ``Hall-T1.xlsx`` を使う場合は ``--input`` は不要）
2. **線形最小二乗法による初期値推定:**
```bash
python mu_fit.py --input data.csv --mode llsq
```
このコマンドは、推定された初期パラメータを ``llsq_params.json`` に保存します。
3. **非線形最適化によるフィッティング:**
```bash
python mu_fit.py --input data.csv --mode fit --eop 0.05 --fix VB
```
このコマンドは、光学フォノンエネルギーを ``0.05 eV`` とし、粒界散乱障壁 ``VB`` を固定してフィッティングを行います。最適化されたパラメータは ``fit_params.json`` に保存され、診断結果は ``fit_fix_suggestions.json`` に保存されます。
4. **各散乱メカニズムの寄与度表示:**
```bash
python mu_fit.py --input data.csv --mode weight
```
このコマンドは、最後に保存されたパラメータ（``fit_params.json`` または ``llsq_params.json``）を使用して、各散乱メカニズムの寄与度をプロットします。
## 関数の詳細
### `load_hall_data()`
Hall効果の測定データをファイルから読み込みます。
指定されたパスのファイル（CSVまたはExcel）を ``pandas.DataFrame`` として読み込みます。ファイルが存在しない場合はエラーメッセージを表示し、``None`` を返します。
* `param str file_path`: 読み込むファイルのパス (``.csv`` または ``.xlsx``)。
* `returns`: 読み込まれたデータフレーム。失敗した場合は ``None``。
* `rtype`: `pandas.DataFrame` or `None`
### `save_params()`
フィッティングパラメータまたは診断結果をJSON形式で保存します。
辞書形式のデータを指定されたファイル名でJSONファイルとして保存します。
* `param dict params`: 保存するパラメータの辞書。
* `param str filename`: 保存先ファイル名。デフォルトは ``fit_params.json``。
* `returns`: `None`
### `load_params()`
JSONファイルからフィッティングパラメータを読み込みます。
指定されたJSONファイルからパラメータを読み込みます。ファイルが存在しない場合は、代わりに ``llsq_params.json`` を試行します。それもなければ、デフォルトの初期値（``{'aop': 1e-3, 'a1': 1e-7, 'a2': 1e-4, 'a3': 1e1, 'VB': 0.0}``）を返します。
* `param str filename`: 読み込み先ファイル名。デフォルトは ``fit_params.json``。
* `returns`: パラメータ名と値の辞書。
* `rtype`: `dict`
### `get_inv_mu_components()`
各散乱メカニズムの逆移動度の寄与と、合計の逆移動度を計算します。
Matthiessenの則に基づき、光学フォノン、音響フォノン、中性不純物、イオン化不純物の各散乱成分の逆移動度を計算します。さらに、粒界散乱の効果を指数関数的に適用して合計逆移動度を算出します。
* `param numpy.ndarray T`: 温度 ``K`` の配列。
* `param list params`: パラメータリスト ``[aop, a1, a2, a3, VB]``。
* `param float Eop`: 光学フォノンエネルギー ``eV``。
* `returns`: 各成分の辞書と、合計逆移動度の配列。
* `rtype`: (`dict`, `numpy.ndarray`)
### `solve_llsq()`
線形最小二乗法を用いて散乱係数(`a_i`)の初期値を推定します。
``VB``（粒界散乱障壁）を0と仮定し、逆移動度に対して線形フィットを行います。推定された係数は負の値にならないようクリッピングされます。
* `param numpy.ndarray T`: 温度 ``K`` の配列。
* `param numpy.ndarray mu_exp`: 実測移動度 ``cm^2/Vs`` の配列。
* `param float Eop`: 光学フォノンエネルギー ``eV``。
* `returns`: 推定された係数配列 ``[aop, a1, a2, a3]``。
* `rtype`: `numpy.ndarray`
### `fit_mask()`
指定された温度範囲に基づいてデータのマスクを作成します。
温度配列 ``T`` の中から、``Tfitmin`` と ``Tfitmax`` の範囲内にある要素に対応する真偽値マスクを生成します。
* `param numpy.ndarray T`: 温度 ``K`` の配列。
* `param float Tfitmin`: フィットに使用する最小温度。
* `param float Tfitmax`: フィットに使用する最大温度。
* `returns`: 指定範囲内のデータを示すブール値配列。
* `rtype`: `numpy.ndarray`
### `build_full_params()`
最適化対象のパラメータと固定パラメータを統合し、全5パラメータのリストを作成します。
初期値/固定値のリストと、最適化によって得られた変動パラメータのリストを組み合わせ、完全なパラメータセットを再構築します。
* `param list init_full`: 初期/固定値を含む全パラメータリスト (5要素)。
* `param list optimize_indices`: 最適化対象のパラメータのインデックスのリスト。
* `param list p_opt`: 最適化されたパラメータ値のリスト。``optimize_indices`` の長さと一致します。
* `returns`: 統合された全パラメータリスト (5要素)。
* `rtype`: `list`
### `residuals_log10()`
移動度の対数値(`log10`)における実測値とモデル値の残差を計算します。
実測移動度と、現在のパラメータで計算されるモデル移動度の常用対数を取り、その差を計算します。数値的な安定性のため、移動度が非常に小さい場合はクリッピングを行います。
* `param numpy.ndarray T`: 温度 ``K`` の配列。
* `param numpy.ndarray mu_exp`: 実測移動度の配列。
* `param list params_full`: 全5パラメータのリスト。
* `param float Eop`: 光学フォノンエネルギー ``eV``。
* `returns`: ``log10(mu_exp) - log10(mu_model)`` の残差配列。
* `rtype`: `numpy.ndarray`
### `numerical_jacobian()`
中心差分法を用いてベクトル値関数のヤコビ行列を数値的に計算します。
各パラメータについて、微小なステップで関数値を perturbed し、その差分から偏微分を近似することでヤコビ行列を構築します。
* `param callable fun_vec`: 残差ベクトルを返す関数。
* `param numpy.ndarray p`: パラメータベクトル。
* `param float rel_step`: 数値微分の相対ステップ。
* `param float abs_step`: 数値微分の絶対ステップ。
* `returns`: ヤコビ行列 (N, M)。
* `rtype`: `numpy.ndarray`
### `param_covariance()`
線形近似に基づき、最適化パラメータの共分散行列と標準誤差を推定します。
最適化後の残差とヤコビ行列を用いて、非線形最小二乗法におけるパラメータの共分散行列を計算し、それから標準誤差を導出します。自由度が0以下の場合、``(None, None)`` を返します。
* `param numpy.ndarray residual_vec`: 最適化後の残差ベクトル。
* `param numpy.ndarray J`: ヤコビ行列。
* `param int dof`: 自由度 (データ数 - パラメータ数)。
* `returns`: 共分散行列と標準誤差の配列のタプル。
* `rtype`: (`numpy.ndarray`, `numpy.ndarray`)
### `cov_to_corr()`
共分散行列を相関係数行列に変換します。
共分散行列から標準偏差を計算し、その情報を用いて相関係数行列を生成します。対角要素は ``1.0`` とし、ゼロ除算を避ける処理を含みます。
* `param numpy.ndarray cov`: 共分散行列。
* `returns`: 相関係数行列と標準偏差の配列。
* `rtype`: (`numpy.ndarray`, `numpy.ndarray`)
### `eigen_sorted_sym()`
対称行列の固有値分解を行い、固有値の降順でソートして返します。
``numpy`` の `eigh()` 関数で対称行列の固有値と固有ベクトルを計算し、固有値の降順でソートして返します。
* `param numpy.ndarray A`: 対象とする対称行列。
* `returns`: 降順ソートされた固有値と対応する固有ベクトル。
* `rtype`: (`numpy.ndarray`, `numpy.ndarray`)
### `summarize_eigenvectors()`
固有ベクトルを人間が読みやすい形式（パラメータ成分の寄与度）で要約します。
共分散行列の固有値と固有ベクトルから、どのパラメータが不確かさの主要な方向に強く寄与しているかをランク付けして表示します。
* `param numpy.ndarray evals`: 固有値の配列。
* `param numpy.ndarray evecs`: 固有ベクトルの行列。
* `param list names`: パラメータ名のリスト。
* `param int topk`: 出力する固有値の数。デフォルトは ``3``。
* `param int compk`: 各固有ベクトルについて表示する成分の数。デフォルトは ``3``。
* `returns`: ランクごとの固有値と成分寄与度を含む辞書のリスト。
* `rtype`: `list`
### `propose_fix_candidates()`
数値的な不安定性（大きな誤差や強相関）に基づき、固定を検討すべきパラメータを提案します。
パラメータの相対誤差、相関係数、および共分散行列の固有分析結果を総合的に評価し、最適化の収束を改善するために固定することが推奨されるパラメータ候補をスコア付けして提示します。
* `param list names`: パラメータ名のリスト。
* `param numpy.ndarray values`: 最適化された値。
* `param numpy.ndarray stderr`: 標準誤差。
* `param numpy.ndarray corr`: 相関係数行列。
* `param numpy.ndarray evals_cov`: 共分散行列の固有値（オプション）。
* `param numpy.ndarray evecs_cov`: 共分散行列の固有ベクトル（オプション）。
* `param float corr_thr`: 強相関とみなす閾値。デフォルトは ``0.95``。
* `param float relerr_thr`: 相対誤差が大きいとみなす閾値。デフォルトは ``0.5``。
* `param int topn`: 提案する最大数。デフォルトは ``3``。
* `returns`: 提案パラメータと理由を含む辞書のリスト。
* `rtype`: `list`
### `prediction_band_log10_mu()`
デルタ法（誤差伝播）を用いて、モデルの対数移動度における予測誤差帯を計算します。
最適化されたパラメータの共分散行列と、モデル関数のパラメータに対する数値的な勾配を用いて、モデル予測値の不確かさ（標準偏差）を計算し、予測帯の上限と下限を導出します。
* `param numpy.ndarray T`: 温度 ``K`` の配列。
* `param list params_full`: 全パラメータのリスト。
* `param float Eop`: 光学フォノンエネルギー ``eV``。
* `param list optimize_indices`: 最適化対象のインデックス。
* `param numpy.ndarray cov_opt`: 最適化パラメータの共分散行列。
* `param float nsigma`: 誤差帯のシグマ倍率。デフォルトは ``1.0``。
* `param float rel_step`: 勾配算出用の相対ステップ。デフォルトは ``1e-6``。
* `param float abs_step`: 勾配算出用の絶対ステップ。デフォルトは ``1e-12``。
* `returns`: (平均値, 下限, 上限) の各配列。
* `rtype`: (`numpy.ndarray`, `numpy.ndarray`, `numpy.ndarray`)
### `visualize_fit()`
実験データ、モデルのフィット曲線、および誤差帯をプロットして保存します。
実測移動度、フィットされたモデル移動度、およびその誤差帯を対数スケールでグラフ化し、指定されたファイル名で画像として保存し、表示します。
* `param numpy.ndarray T`: 温度 ``K``。
* `param numpy.ndarray mu_exp`: 実測移動度。
* `param numpy.ndarray mu_fit`: モデルによるフィット値（オプション）。
* `param numpy.ndarray mu_lo`: 誤差帯の下限（オプション）。
* `param numpy.ndarray mu_hi`: 誤差帯の上限（オプション）。
* `param str title`: グラフのタイトル。デフォルトは ``Hall Mobility Fit``。
* `param str save_name`: 保存ファイル名。デフォルトは ``plot.png``。
* `returns`: `None`
### `visualize_weights()`
各散乱メカニズムが全体の散乱（逆移動度）に占める割合（寄与度）を可視化します。
温度に対して、各散乱メカニズム（光学フォノン、音響フォノン、不純物、粒界）が合計の逆移動度（散乱）に占める割合をパーセンテージでプロットします。
* `param numpy.ndarray T`: 温度 ``K``。
* `param dict components`: 各成分の逆移動度。
* `param numpy.ndarray total`: 合計逆移動度。
* `param str save_name`: 保存ファイル名。デフォルトは ``weight_plot.png``。
* `returns`: `None`
### `main()`
コマンドライン引数を解析し、Hall効果移動度の解析・フィッティング処理を実行するエントリーポイント。
`argparse` を使用してコマンドライン引数を処理し、指定されたモード（データ読み込み、LLSQ初期推定、非線形フィット、寄与度可視化）に応じて、Hall効果測定データのフィッティングと結果の保存・可視化を行います。パラメータ固定、誤差推定、フィッティング温度範囲の指定など、詳細なオプションに対応します。
* `returns`: `None`
## アルゴリズムの概要
### 線形最小二乗法 (LLSQ)
``solve_llsq()`` 関数で使用され、移動度モデルの初期パラメータを効率的に推定するために利用されます。粒界散乱障壁 ``VB`` を ``0`` と仮定し、逆移動度を各散乱成分の線形結合として扱い、 ``numpy.linalg.lstsq`` を用いて係数を求めます。
### 非線形最適化 (`scipy.optimize.minimize`)
``main()`` 関数の ``fit`` モードで、``objective()`` 関数（対数スケールでの残差二乗和を最小化）を最適化するために使用されます。これにより、`[aop, a1, a2, a3, VB]` の5つのパラメータを非線形にフィッティングします。最適化メソッドは ``--method`` オプションで指定可能です。
### 誤差伝播 (Delta Method)
``prediction_band_log10_mu()`` 関数でモデルの予測誤差帯を計算するために適用されます。最適化されたパラメータの共分散行列と、モデル関数の各パラメータに対する数値的な勾配（ヤコビ行列）を用いて、モデル予測値の不確かさを推定します。これにより、フィッティングされたモデルの信頼区間（予測帯）が提供されます。
### パラメータ診断
フィッティング結果の安定性と信頼性を評価するために、以下の診断が行われます。
* **共分散行列と標準誤差**: ``param_covariance()`` 関数で計算され、各パラメータの推定の不確かさを示します。
* **相関係数行列**: ``cov_to_corr()`` 関数で計算され、パラメータ間の相関の強さを示します。特に強い相関は、フィッティングが不安定である可能性を示唆します。
* **固有値解析**: ``eigen_sorted_sym()`` と ``summarize_eigenvectors()`` 関数で共分散行列に対して行われ、不確かさの主要な方向と、それに寄与するパラメータを特定します。
* **固定候補の提案**: ``propose_fix_candidates()`` 関数は、これらの診断結果（相対誤差が大きい、強相関がある、固有方向で支配的であるなど）に基づいて、フィッティングを改善するために固定を検討すべきパラメータを提案します。
## 入出力ファイル形式
### 入力データ
* **ファイル形式**: CSV (``.csv``) または Excel (``.xlsx``)
* **内容**: 温度 (``T``) とホール移動度 (``mu_exp``) のデータ。コマンドラインオプション ``--temp_col`` と ``--mu_col`` で列インデックスを指定します。
* **例**:
```
Temperature(K),Carrier Concentration(cm-3),Mobility(cm2/Vs)
300,1.23e18,150.5
250,1.18e18,200.1
...
```
上記の場合、``--temp_col 0 --mu_col 2`` を指定します。
### 出力データ
* **JSONファイル**:
* ``llsq_params.json``: LLSQモードで推定された初期パラメータが保存されます。
* ``fit_params.json``: フィッティングモードで最適化された最終パラメータが保存されます。
* ``fit_fix_suggestions.json``: フィッティングモードでのパラメータ診断結果（固定候補の提案）が保存されます。
* **画像ファイル**:
* ``plot.png`` (または ``visualize_fit()`` で指定されるファイル名): フィット曲線と実験データをプロットしたグラフ。
* ``weight_plot.png`` (または ``visualize_weights()`` で指定されるファイル名): 各散乱メカニズムの寄与度をプロットしたグラフ。
## 定数
* ``K_B``: ボルツマン定数。値は ``8.617333262e-5`` で、単位は ``eV/K`` です。
````