#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Hall効果移動度解析・フィッティングツール (mu_fit.py)
====================================================
このモジュールは、Hall効果測定から得られた温度依存移動度データを、
各種散乱メカニズム（光学フォノン、音響フォノン、中性不純物、イオン化不純物、粒界散乱）
の物理モデルに基づいてフィッティングし、パラメータの最適化および誤差評価を行います。
主な機能:
* 実験データの読み込み（CSV/Excel）
* 線形最小二乗法（LLSQ）による初期値推定
* 非線形最適化によるパラメータ抽出
* 誤差伝播法によるモデルの不確かさ（1σ範囲）の算出
* 共分散行列・相関係数・固有値解析によるパラメータ感度診断
関連ドキュメント:
* :doc:`usage`
* :doc:`physics_models`
.. note::
フィッティングの目的関数は、移動度の対数スケールでの残差二乗和を最小化します。
"""
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import argparse
import os
import json
# 定数: ボルツマン定数 (eV/K)
K_B = 8.617333262e-5
# ---------------------------------------------------------
# 1. データ入出力関連の関数
# ---------------------------------------------------------
[ドキュメント]
def load_hall_data(file_path):
"""
Hall効果の測定データをファイルから読み込む。
指定されたパスのファイル（CSVまたはExcel）をPandas DataFrameとして読み込む。
ファイルが存在しない場合はエラーメッセージを表示し、Noneを返す。
:param str file_path: 読み込むファイルのパス (.csv または .xlsx)。
:returns: 読み込まれたデータフレーム。失敗した場合は None。
:rtype: pandas.DataFrame or None
"""
if not os.path.exists(file_path):
print(f"エラー: ファイル '{file_path}' が見つかりません。")
return None
if file_path.endswith('.csv'):
return pd.read_csv(file_path)
return pd.read_excel(file_path)
[ドキュメント]
def save_params(params, filename='fit_params.json'):
"""
フィッティングパラメータまたは診断結果をJSON形式で保存する。
辞書形式のデータを指定されたファイル名でJSONファイルとして保存する。
:param dict params: 保存するパラメータの辞書。
:param str filename: 保存先ファイル名。
:returns: None
"""
with open(filename, 'w', encoding='utf-8') as f:
json.dump(params, f, indent=4, ensure_ascii=False)
print(f"パラメータを保存しました: {filename}")
[ドキュメント]
def load_params(filename='fit_params.json'):
"""
JSONファイルからフィッティングパラメータを読み込む。
指定されたJSONファイルからパラメータを読み込む。ファイルが存在しない場合は、
`llsq_params.json`を試行し、それもなければデフォルトの初期値を返す。
:param str filename: 読み込み先ファイル名。
:returns: パラメータ名と値の辞書。
:rtype: dict
"""
if not os.path.exists(filename):
alt_filename = 'llsq_params.json'
if os.path.exists(alt_filename):
filename = alt_filename
else:
print(f"警告: パラメータファイルが見つかりません。デフォルト値を使用します。")
return {'aop': 1e-3, 'a1': 1e-7, 'a2': 1e-4, 'a3': 1e1, 'VB': 0.0}
with open(filename, 'r', encoding='utf-8') as f:
return json.load(f)
# ---------------------------------------------------------
# 2. 物理モデルと解析関連の関数
# ---------------------------------------------------------
[ドキュメント]
def get_inv_mu_components(T, params, Eop):
"""
各散乱メカニズムの逆移動度の寄与と、合計の逆移動度を計算する。
Matthiessenの則に基づき、光学フォノン、音響フォノン、中性不純物、イオン化不純物の
各散乱成分の逆移動度を計算し、さらに粒界散乱の効果を指数関数的に適用して合計逆移動度を算出する。
:param numpy.ndarray T: 温度 [K] の配列。
:param list params: パラメータリスト [aop, a1, a2, a3, VB]。
:param float Eop: 光学フォノンエネルギー [eV]。
:returns: 各成分の辞書と、合計逆移動度の配列。
:rtype: (dict, numpy.ndarray)
"""
aop, a1, a2, a3, VB = params
# Optical phonon factor (Bose-Einstein)
f_op = 1.0 / (np.exp(Eop / (K_B * T)) - 1.0)
f_ac = T**1.5
f_ni = np.ones_like(T)
f_ii = T**-1.5
components = {
'Optical Phonon': aop * f_op,
'Acoustic Phonon': a1 * f_ac,
'Neutral Impurity': a2 * f_ni,
'Ionized Impurity': a3 * f_ii
}
inv_mu_bulk = np.maximum(sum(components.values()), 1e-30)
# exp(VB/kT) can overflow for small T; clip exponent safely.
expo = VB / (K_B * T)
expo = np.clip(expo, -700.0, 700.0)
exp_factor = np.exp(expo)
inv_mu_total = inv_mu_bulk * exp_factor
components['Grain Boundary'] = inv_mu_total - inv_mu_bulk
return components, inv_mu_total
[ドキュメント]
def solve_llsq(T, mu_exp, Eop):
"""
線形最小二乗法を用いて散乱係数(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
"""
f_op = 1.0 / (np.exp(Eop / (K_B * T)) - 1.0)
f_ac = T**1.5
f_ni = np.ones_like(T)
f_ii = T**-1.5
X = np.column_stack([f_op, f_ac, f_ni, f_ii])
y = 1.0 / mu_exp
coeffs, _, _, _ = np.linalg.lstsq(X, y, rcond=None)
# 負の値にならないようクリップ
return np.maximum(coeffs, 0)
[ドキュメント]
def fit_mask(T, Tfitmin=-1e100, Tfitmax=+1e100):
"""
指定された温度範囲に基づいてデータのマスクを作成する。
温度配列`T`の中から、`Tfitmin`と`Tfitmax`の範囲内にある要素に対応する真偽値マスクを生成する。
:param numpy.ndarray T: 温度 [K] の配列。
:param float Tfitmin: フィットに使用する最小温度。
:param float Tfitmax: フィットに使用する最大温度。
:returns: 指定範囲内のデータを示すブール値配列。
:rtype: numpy.ndarray
"""
return (T >= Tfitmin) & (T <= Tfitmax)
[ドキュメント]
def build_full_params(init_full, optimize_indices, p_opt):
"""
最適化対象のパラメータと固定パラメータを統合し、全5パラメータのリストを作成する。
初期値/固定値のリストと、最適化によって得られた変動パラメータのリストを組み合わせ、
完全なパラメータセットを再構築する。
:param list init_full: 初期/固定値を含む全パラメータリスト。
:param list optimize_indices: 最適化対象のインデックス。
:param list p_opt: 最適化されたパラメータ値のリスト。
:returns: 統合された全パラメータリスト。
:rtype: list
"""
p = list(init_full)
for idx, val in zip(optimize_indices, p_opt):
p[idx] = val
return p
[ドキュメント]
def residuals_log10(T, mu_exp, params_full, Eop):
"""
移動度の対数値(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
"""
_, inv_total = get_inv_mu_components(T, params_full, Eop)
mu_model = 1.0 / np.maximum(inv_total, 1e-300)
mu_model = np.maximum(mu_model, 1e-300)
mu_exp_safe = np.maximum(mu_exp, 1e-300)
return np.log10(mu_exp_safe) - np.log10(mu_model)
[ドキュメント]
def numerical_jacobian(fun_vec, p, rel_step=1e-6, abs_step=1e-12):
"""
中心差分法を用いてベクトル値関数のヤコビ行列を数値的に計算する。
各パラメータについて、微小なステップで関数値を perturbed し、その差分から偏微分を近似することで
ヤコビ行列を構築する。
:param callable fun_vec: 残差ベクトルを返す関数。
:param numpy.ndarray p: パラメータベクトル。
:param float rel_step: 数値微分の相対ステップ。
:param float abs_step: 数値微分の絶対ステップ。
:returns: ヤコビ行列 (N, M)。
:rtype: numpy.ndarray
"""
p = np.asarray(p, dtype=float)
f0 = np.asarray(fun_vec(p), dtype=float)
n = f0.size
m = p.size
J = np.zeros((n, m), dtype=float)
for j in range(m):
dp = rel_step * (abs(p[j]) + 1.0) + abs_step
p1 = p.copy(); p1[j] += dp
p2 = p.copy(); p2[j] -= dp
f1 = np.asarray(fun_vec(p1), dtype=float)
f2 = np.asarray(fun_vec(p2), dtype=float)
J[:, j] = (f1 - f2) / (2.0 * dp)
return J
[ドキュメント]
def param_covariance(residual_vec, J, dof):
"""
線形近似に基づき、最適化パラメータの共分散行列と標準誤差を推定する。
最適化後の残差とヤコビ行列を用いて、非線形最小二乗法におけるパラメータの共分散行列を計算し、
それから標準誤差を導出する。
:param numpy.ndarray residual_vec: 最適化後の残差ベクトル。
:param numpy.ndarray J: ヤコビ行列。
:param int dof: 自由度 (データ数 - パラメータ数)。
:returns: 共分散行列と標準誤差の配列のタプル。
:rtype: (numpy.ndarray, numpy.ndarray)
"""
if dof <= 0:
return None, None
s2 = float(np.sum(residual_vec**2) / dof)
JTJ = J.T @ J
try:
cov = s2 * np.linalg.inv(JTJ)
except np.linalg.LinAlgError:
cov = s2 * np.linalg.pinv(JTJ)
stderr = np.sqrt(np.maximum(np.diag(cov), 0.0))
return cov, stderr
# ============================================================
# covariance / correlation / eigen diagnostics
# ============================================================
[ドキュメント]
def cov_to_corr(cov):
"""
共分散行列を相関係数行列に変換する。
共分散行列から標準偏差を計算し、その情報を用いて相関係数行列を生成する。
対角要素は1.0とし、ゼロ除算を避ける処理を含む。
:param numpy.ndarray cov: 共分散行列。
:returns: 相関係数行列と標準偏差の配列。
:rtype: (numpy.ndarray, numpy.ndarray)
"""
cov = np.asarray(cov, dtype=float)
std = np.sqrt(np.maximum(np.diag(cov), 0.0))
denom = np.outer(std, std)
with np.errstate(invalid='ignore', divide='ignore'):
corr = cov / denom
np.fill_diagonal(corr, 1.0)
z = denom == 0
corr[z] = np.nan
np.fill_diagonal(corr, 1.0)
return corr, std
[ドキュメント]
def eigen_sorted_sym(A):
"""
対称行列の固有値分解を行い、固有値の降順でソートして返す。
numpyの`eigh`関数で対称行列の固有値と固有ベクトルを計算し、
固有値の降順でソートして返す。
:param numpy.ndarray A: 対象とする対称行列。
:returns: 降順ソートされた固有値と対応する固有ベクトル。
:rtype: (numpy.ndarray, numpy.ndarray)
"""
A = np.asarray(A, dtype=float)
evals, evecs = np.linalg.eigh(A)
idx = np.argsort(evals)[::-1] # descending
return evals[idx], evecs[:, idx]
[ドキュメント]
def summarize_eigenvectors(evals, evecs, names, topk=3, compk=3):
"""
固有ベクトルを人間が読みやすい形式（パラメータ成分の寄与度）で要約する。
共分散行列の固有値と固有ベクトルから、どのパラメータが不確かさの主要な方向に
強く寄与しているかをランク付けして表示する。
:param numpy.ndarray evals: 固有値の配列。
:param numpy.ndarray evecs: 固有ベクトルの行列。
:param list names: パラメータ名のリスト。
:param int topk: 出力する固有値の数。
:param int compk: 各固有ベクトルについて表示する成分の数。
:returns: ランクごとの固有値と成分寄与度を含む辞書のリスト。
:rtype: list
"""
names = list(names)
out = []
k = min(topk, len(evals))
for r in range(k):
v = evecs[:, r]
order = np.argsort(np.abs(v))[::-1]
comps = []
for j in order[:min(compk, len(order))]:
comps.append((names[j], float(v[j])))
out.append({
'rank': r + 1,
'eigenvalue': float(evals[r]),
'components': comps
})
return out
[ドキュメント]
def propose_fix_candidates(names, values, stderr, corr, evals_cov=None, evecs_cov=None,
corr_thr=0.95, relerr_thr=0.5, topn=3):
"""
数値的な不安定性（大きな誤差や強相関）に基づき、固定を検討すべきパラメータを提案する。
パラメータの相対誤差、相関係数、および共分散行列の固有分析結果を総合的に評価し、
最適化の収束を改善するために固定することが推奨されるパラメータ候補をスコア付けして提示する。
: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: 強相関とみなす閾値。
:param float relerr_thr: 相対誤差が大きいとみなす閾値。
:param int topn: 提案する最大数。
:returns: 提案パラメータと理由を含む辞書のリスト。
:rtype: list
"""
names = list(names)
values = np.asarray(values, dtype=float)
stderr = np.asarray(stderr, dtype=float)
tiny = 1e-30
relerr = stderr / (np.abs(values) + tiny)
score = relerr.copy()
reasons = {n: [] for n in names}
for i, n in enumerate(names):
if not np.isfinite(relerr[i]):
continue
if relerr[i] >= relerr_thr:
reasons[n].append(f"relative stderr = {relerr[i]:.3g} (>= {relerr_thr})")
else:
reasons[n].append(f"relative stderr = {relerr[i]:.3g}")
corr = np.asarray(corr, dtype=float)
for i in range(len(names)):
for j in range(i+1, len(names)):
cij = corr[i, j]
if not np.isfinite(cij):
continue
if abs(cij) >= corr_thr:
bonus = (abs(cij) - corr_thr) * 5.0 + 0.5
score[i] += bonus
score[j] += bonus
reasons[names[i]].append(f"strong corr with {names[j]}: {cij:+.3f} (>= {corr_thr})")
reasons[names[j]].append(f"strong corr with {names[i]}: {cij:+.3f} (>= {corr_thr})")
if evals_cov is not None and evecs_cov is not None and len(evals_cov) > 0:
v = evecs_cov[:, 0]
w = np.abs(v)
w = w / (w.max() + tiny)
for i, n in enumerate(names):
if w[i] >= 0.5:
add = 0.7 * w[i]
score[i] += add
reasons[n].append(f"dominant in most-uncertain eigen-direction: |v|/max={w[i]:.2f}")
items = []
for i, n in enumerate(names):
items.append({
'param': n,
'score': float(score[i]),
'value': float(values[i]),
'stderr': float(stderr[i]),
'relerr': float(relerr[i]) if np.isfinite(relerr[i]) else None,
'reasons': reasons[n],
})
items.sort(key=lambda d: d['score'], reverse=True)
filtered = []
for it in items:
if (it['relerr'] is not None and it['relerr'] >= relerr_thr) or any('strong corr' in r for r in it['reasons']):
filtered.append(it)
if not filtered:
filtered = items[:min(topn, len(items))]
return filtered[:min(topn, len(filtered))]
[ドキュメント]
def prediction_band_log10_mu(T, params_full, Eop, optimize_indices, cov_opt, nsigma=1.0,
rel_step=1e-6, abs_step=1e-12):
"""
デルタ法（誤差伝播）を用いて、モデルの対数移動度における予測誤差帯を計算する。
最適化されたパラメータの共分散行列と、モデル関数のパラメータに対する数値的な勾配を用いて、
モデル予測値の不確かさ（標準偏差）を計算し、予測帯の上限と下限を導出する。
: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: 誤差帯のシグマ倍率。
:param float rel_step: 勾配算出用の相対ステップ。
:param float abs_step: 勾配算出用の絶対ステップ。
:returns: (平均値, 下限, 上限) の各配列。
:rtype: (numpy.ndarray, numpy.ndarray, numpy.ndarray)
"""
_, inv0 = get_inv_mu_components(T, params_full, Eop)
mu0 = 1.0 / np.maximum(inv0, 1e-300)
y0 = np.log10(np.maximum(mu0, 1e-300))
m = len(optimize_indices)
if cov_opt is None or m == 0:
return y0, None, None
grads = np.zeros((len(T), m), dtype=float)
p0_opt = np.array([params_full[i] for i in optimize_indices], dtype=float)
for j in range(m):
dp = rel_step * (abs(p0_opt[j]) + 1.0) + abs_step
p_plus = p0_opt.copy(); p_plus[j] += dp
params_plus = list(params_full)
for idx, val in zip(optimize_indices, p_plus):
params_plus[idx] = val
_, invp = get_inv_mu_components(T, params_plus, Eop)
mup = 1.0 / np.maximum(invp, 1e-300)
yp = np.log10(np.maximum(mup, 1e-300))
p_minus = p0_opt.copy(); p_minus[j] -= dp
params_minus = list(params_full)
for idx, val in zip(optimize_indices, p_minus):
params_minus[idx] = val
_, invm = get_inv_mu_components(T, params_minus, Eop)
mum = 1.0 / np.maximum(invm, 1e-300)
ym = np.log10(np.maximum(mum, 1e-300))
grads[:, j] = (yp - ym) / (2.0 * dp)
var_y = np.einsum('ni,ij,nj->n', grads, cov_opt, grads)
sig_y = np.sqrt(np.maximum(var_y, 0.0))
ylow = y0 - nsigma * sig_y
yhigh = y0 + nsigma * sig_y
return y0, ylow, yhigh
# ---------------------------------------------------------
# 3. 可視化関連の関数
# ---------------------------------------------------------
[ドキュメント]
def visualize_fit(T, mu_exp, mu_fit=None, mu_lo=None, mu_hi=None,
title='Hall Mobility Fit', save_name='plot.png'):
"""
実験データ、モデルのフィット曲線、および誤差帯をプロットして保存する。
実測移動度、フィットされたモデル移動度、およびその誤差帯を対数スケールでグラフ化し、
指定されたファイル名で画像として保存し、表示する。
: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: グラフのタイトル。
:param str save_name: 保存ファイル名。
:returns: None
"""
plt.figure(figsize=(8, 6))
plt.scatter(T, mu_exp, color='red', label='Experimental', alpha=0.6)
idx = np.argsort(T)
T_sorted = T[idx]
if mu_fit is not None:
mu_fit_sorted = np.asarray(mu_fit)[idx]
plt.plot(T_sorted, mu_fit_sorted, color='blue', label='Model Fit', linewidth=2)
if mu_lo is not None and mu_hi is not None:
mu_lo_sorted = np.asarray(mu_lo)[idx]
mu_hi_sorted = np.asarray(mu_hi)[idx]
plt.fill_between(T_sorted, mu_lo_sorted, mu_hi_sorted, alpha=0.25, label='Model ±1σ', color='skyblue')
plt.xlabel('Temperature (K)')
plt.ylabel('Mobility (cm²/Vs)')
plt.title(title)
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig(save_name)
plt.show()
[ドキュメント]
def visualize_weights(T, components, total, save_name='weight_plot.png'):
"""
各散乱メカニズムが全体の散乱（逆移動度）に占める割合（寄与度）を可視化する。
温度に対して、各散乱メカニズム（光学フォノン、音響フォノン、不純物、粒界）が
合計の逆移動度（散乱）に占める割合をパーセンテージでプロットする。
:param numpy.ndarray T: 温度 [K]。
:param dict components: 各成分の逆移動度。
:param numpy.ndarray total: 合計逆移動度。
:param str save_name: 保存ファイル名。
:returns: None
"""
plt.figure(figsize=(10, 6))
idx = np.argsort(T)
T_sorted = T[idx]
for name, inv_mu in components.items():
weight = (inv_mu / total) * 100
plt.plot(T_sorted, weight[idx], marker='o', markersize=4, label=name, linestyle='-', alpha=0.8)
plt.xlabel('Temperature (K)')
plt.ylabel('Contribution to Scattering (%)')
plt.title('Scattering Mechanism Contributions')
plt.ylim(-5, 105)
plt.grid(True, linestyle='--', alpha=0.5)
plt.legend(loc='best')
plt.tight_layout()
plt.savefig(save_name)
plt.show()
# ---------------------------------------------------------
# 4. メイン処理
# ---------------------------------------------------------
[ドキュメント]
def main():
"""
コマンドライン引数を解析し、Hall効果移動度の解析・フィッティング処理を実行するエントリーポイント。
`argparse`を使用してコマンドライン引数を処理し、指定されたモード
（データ読み込み、LLSQ初期推定、非線形フィット、寄与度可視化）に応じて、
Hall効果測定データのフィッティングと結果の保存・可視化を行う。
パラメータ固定、誤差推定、フィッティング温度範囲の指定など、詳細なオプションに対応する。
:returns: None
"""
parser = argparse.ArgumentParser(description='Hall効果移動度フィッティング（パラメータ固定＆誤差推定付き）')
parser.add_argument('--input', type=str, default='Hall-T1.xlsx', help='入力ファイル名')
parser.add_argument('--temp_col', type=int, default=0, help='温度列(0開始)')
parser.add_argument('--mu_col', type=int, default=2, help='移動度列(0開始)')
parser.add_argument('--mode', type=str, choices=['read', 'llsq', 'fit', 'weight'], default='read')
parser.add_argument('--method', type=str, default='Nelder-Mead')
parser.add_argument('--eop', type=float, default=0.045, help='光学フォノンエネルギー(eV)')
parser.add_argument('--fix', type=str, default='', help='固定するパラメータをカンマ区切りで指定 (例: a2,VB)')
parser.add_argument('--Tfitmin', type=float, default=-1e100, help='フィットに使う最小温度[K]')
parser.add_argument('--Tfitmax', type=float, default=+1e100, help='フィットに使う最大温度[K]')
parser.add_argument('--band_sigma', type=float, default=1.0, help='誤差帯の幅 (nsigma). default 1.0 (±1σ)')
parser.add_argument('--jac_relstep', type=float, default=1e-6, help='数値微分の相対ステップ')
parser.add_argument('--jac_absstep', type=float, default=1e-12, help='数値微分の絶対ステップ')
args = parser.parse_args()
print()
print(f"input file={args.input}")
print(f"mode={args.mode}")
print(f"method={args.method}")
print(f"eop={args.eop} eV")
print(f"fix={args.fix}")
print(f"Tfitmin={args.Tfitmin}, Tfitmax={args.Tfitmax}")
df = load_hall_data(args.input)
if df is None:
return
T = df.iloc[:, args.temp_col].values.astype(float)
mu_exp = df.iloc[:, args.mu_col].values.astype(float)
labels = ['aop', 'a1', 'a2', 'a3', 'VB']
fixed_list = [s.strip() for s in args.fix.split(',') if s.strip()] if args.fix else []
if args.mode == 'read':
visualize_fit(T, mu_exp, title='Experimental Data')
elif args.mode == 'llsq':
c = solve_llsq(T, mu_exp, args.eop)
params_dict = {'aop': float(c[0]), 'a1': float(c[1]), 'a2': float(c[2]), 'a3': float(c[3]), 'VB': 0.0}
save_params(params_dict, filename='llsq_params.json')
_, inv_fit = get_inv_mu_components(T, list(c) + [0.0], args.eop)
visualize_fit(T, mu_exp, 1 / inv_fit, title='LLSQ Initial Fit (VB=0)')
elif args.mode == 'fit':
mfit = fit_mask(T, args.Tfitmin, args.Tfitmax)
if np.sum(mfit) < 3:
print("エラー: フィットに使える点が少なすぎます。")
return
T_fit = T[mfit]
mu_fit_exp = mu_exp[mfit]
p_base = load_params()
init_full = [float(p_base.get(l, 0.0)) for l in labels]
optimize_indices = [i for i, l in enumerate(labels) if l not in fixed_list]
def objective(p_opt):
p_current = build_full_params(init_full, optimize_indices, p_opt)
if any(np.array(p_current)[:-1] < 0):
return 1e18
r = residuals_log10(T_fit, mu_fit_exp, p_current, args.eop)
return float(np.sum(r**2))
init_opt = [init_full[i] for i in optimize_indices]
res = minimize(objective, init_opt, method=args.method, options={'maxiter': 5000})
final_values = build_full_params(init_full, optimize_indices, res.x)
final_params = {l: float(v) for l, v in zip(labels, final_values)}
m = len(optimize_indices)
dof = int(np.sum(mfit) - m)
cov_opt, stderr_opt = None, None
if m > 0 and dof > 0:
def r_of_popt(p_opt):
p_full = build_full_params(init_full, optimize_indices, p_opt)
return residuals_log10(T_fit, mu_fit_exp, p_full, args.eop)
rvec = r_of_popt(res.x)
J = numerical_jacobian(r_of_popt, np.asarray(res.x), rel_step=args.jac_relstep, abs_step=args.jac_absstep)
cov_opt, stderr_opt = param_covariance(rvec, J, dof)
# 診断出力・可視化
if cov_opt is not None:
opt_names = [labels[i] for i in optimize_indices]
opt_vals = [final_values[i] for i in optimize_indices]
corr_opt, _ = cov_to_corr(cov_opt)
evals_cov, evecs_cov = eigen_sorted_sym(cov_opt)
suggestions = propose_fix_candidates(opt_names, opt_vals, stderr_opt, corr_opt, evals_cov, evecs_cov)
save_params({'suggestions': suggestions}, filename='fit_fix_suggestions.json')
y0, ylo, yhi = prediction_band_log10_mu(T, final_values, args.eop, optimize_indices, cov_opt, nsigma=args.band_sigma)
mu_lo = 10**ylo if ylo is not None else None
mu_hi = 10**yhi if yhi is not None else None
visualize_fit(T, mu_exp, 10**y0, mu_lo, mu_hi, title=f'Final Fit (±{args.band_sigma}σ)')
elif args.mode == 'weight':
p_dict = load_params()
p_list = [float(p_dict.get(l, 0.0)) for l in labels]
comp, total = get_inv_mu_components(T, p_list, args.eop)
visualize_weights(T, comp, total)
if __name__ == '__main__':
main()