"""
半導体の散乱モデルの最小二乗法
概要:
半導体のホール移動度に対する散乱モデルの最小二乗フィッティングを行うスクリプトです。
詳細説明:
本スクリプトは、温度およびキャリア濃度依存性を持つホール移動度データを、
複数の散乱メカニズム（光学フォノン散乱、音響フォノン散乱、イオン化不純物散乱、障壁散乱など）を考慮したモデルでフィッティングします。
これにより、各散乱パラメータを決定し、その結果をグラフ表示するとともに、
Excelファイルおよびパラメータファイルに保存します。
動作モードは 'init', 'fit', 'sim' があり、'fit' モードではさらに 'mu_T', 'mu_TN', 'mu_T_Nfix',
'Aop', 'A_T0', 'A_T32', 'VB' のターゲットを指定して特定のフィッティングを実行できます。
関連リンク:
:doc:`mu-TN-Fit_usage`
"""
import sys
import numpy as np
from numpy import sin, cos, tan, pi, exp, log, sqrt
from scipy.optimize import minimize
import pandas as pd
from matplotlib import pyplot as plt
import matplotlib.widgets as wg
from tklib.tkutils import getarg, getintarg, getfloatarg, pconv_by_type
from tklib.tkapplication import tkApplication
from tklib.tkparams import tkParams
from tklib.tkvariousdata import tkVariousData
from tklib.tksci.tksci import log10, e, kB
from tklib.tksci.tkFit import tkFit
from tklib.tksci.tkFit_m import tkFit_m
usage_str = '''
""
f"(i) usage: python {sys.argv[0]} fit mode target infile xlabel ylabel"
'''[1:-1]
[ドキュメント]
def initialize():
"""
アプリケーションとパラメータを初期化します。
詳細説明:
tkApplicationとtkParamsのインスタンスを作成し、アプリケーション全体のデフォルト設定や
フィッティングパラメータの初期値を設定します。
主に、モード、ターゲット、入力ファイル、データカラム、温度範囲、
最適化アルゴリズム、フィッティングパラメータの初期値などを設定します。
:returns:
tuple[tkApplication, tkParams]: 初期化されたアプリケーションオブジェクトと共通パラメータオブジェクトのタプル。
"""
#================================
# Global variables
#================================
app = tkApplication(usage_str = usage_str)
argv, narg = app.get_argv()
app.cparams = tkParams()
cparams = app.cparams
cparams.debug = 0
cparams.print_level = 0
# mode: 'init', 'fit', 'sim'
# cparams.mode = 'init'
cparams.mode = 'fit'
# cparams.target = 'Aop'
# cparams.target = 'A_T0'
cparams.target = 'A_T32'
# cparams.target = 'VB'
cparams.fix_N = False
# data file to fit (.xlsx)
cparams.infile = 'IO.xlsx'
cparams.parameterfile = None
cparams.parameterbkfile = None
cparams.xlabel = 2
cparams.ylabel = 6
cparams.sample_label = 0
cparams.formula_label = 1
cparams.T_label = 6
cparams.mu_label = 7
cparams.Nlabel = 8
cparams.Tbase = 300.0
# Temperature range
cparams.Tmin = 50 # K
cparams.Tmax = 1000
cparams.Tstep = 10
cparams.nT = int((cparams.Tmax - cparams.Tmin) / cparams.Tstep + 1)
# 最適化パラメータの初期値
#Aop
app.Aop_params = tkParams()
app.Aop_params.varname = ["Eop_0", "A1", "B1", "C1", "logN01"]
app.Aop_params.varunit = ["eV", "", "", "", ""]
app.Aop_params.value = [0.046, 0.0143, 20.06, 0.0, 19.58]
app.Aop_params.optid = [0, 1, 0, 0, 0]
#AT0
app.AT0_params = tkParams()
app.AT0_params.varname = ["A2", "B2", "C2", "logN02"]
app.AT0_params.varunit = ["", "", "", ""]
app.AT0_params.value = [0.0135, 16.15, 0.00647, 18.994]
app.AT0_params.optid = [1, 1, 1, 1]
#A_T3/2
app.AT32_params = tkParams()
app.AT32_params.varname = ["A3", "B3", "C3", "D3", "F3", "logN03"]
app.AT32_params.varunit = ["", "", "", "", "", ""]
app.AT32_params.value = [8.00e-36, 15.0, 0.0, 0.0, 0.0, 18.5]
app.AT32_params.optid = [1, 1, 1, 1, 1, 1]
#VB
app.VB_params = tkParams()
app.VB_params.varname = ["A6", "B6", "C6", "D6", "logN06"]
app.VB_params.varunit = ["", "", "", "", ""]
app.VB_params.value = [4.64e-2, 10.80, 0.0, 0.031, 18.37]
app.VB_params.optid = [1, 1, 1, 1, 1]
#scipy.optimize.minimizeで使うアルゴリズム
#nelder-mead Downhill simplex
#cg conjugate gradient (Polak-Ribiere method)
#bfgs BFGS法
#newton-cg Newton-CG
cparams.method = "nelder-mead"
cparams.nmaxiter = 3000
cparams.tol = 1.0e-3
cparams.h_diff = 1.0e-3
cparams.print_interval = 5
cparams.plot_interval = 1
#=============================
# Graph configuration
#=============================
cparams.fig = None
cparams.figsize = [6, 6]
cparams.fontsize = 14
cparams.legend_fontsize = 12
cparams.graphupdateinterval = 10
return app, cparams
#=============================
# Treat argments
#=============================
[ドキュメント]
def usage(app):
"""
スクリプトの使用方法を標準出力に表示します。
詳細説明:
`app` オブジェクトに格納されている使用方法文字列を整形して表示します。
:param app: (tkApplication) アプリケーションオブジェクト。
:returns: None
"""
cparams = app.cparams
for s in app.usage_str.split('\n'):
cmd = 'print({})'.format(s.rstrip())
eval(cmd)
[ドキュメント]
def update_vars(app, cparams):
"""
コマンドライン引数に基づいてパラメータを更新します。
詳細説明:
スクリプト起動時に指定されたコマンドライン引数を解析し、
共通パラメータ `cparams` の各値を上書きします。
引数が指定されていない場合は、現在の `cparams` の値が維持されます。
:param app: (tkApplication) アプリケーションオブジェクト。
:param cparams: (tkParams) 更新対象の共通パラメータオブジェクト。
:returns: None
"""
argv = app.argv
# if len(argv) <= 1:
# app.terminate(usage = usage)
cparams.mode = getarg(1, cparams.mode)
cparams.target = getarg(2, cparams.target)
cparams.infile = getarg(3, cparams.infile)
cparams.xlabel = getarg(4, cparams.xlabel)
cparams.ylabel = getarg(5, cparams.ylabel)
cparams.Nlabel = getarg(6, cparams.Nlabel)
[ドキュメント]
def Aop(N, Eop_0, A1, B1, C1, logN01):
"""
光学フォノン散乱に関連するパラメータ (Eop_0, Aop) を計算します。
詳細説明:
キャリア濃度 `N` に基づいて、光学フォノン散乱のパラメータである `Eop_0` と `Aop` を計算します。
`Aop` はフェルミ分布関数に似た形でキャリア濃度に依存します。
:param N: (float) キャリア濃度 (cm^-3)。
:param Eop_0: (float) 光学フォノン散乱のエネルギーパラメータ。
:param A1: (float) 光学フォノン散乱の振幅パラメータ。
:param B1: (float) 光学フォノン散乱の濃度依存性パラメータ。
:param C1: (float) 光学フォノン散乱のオフセットパラメータ。
:param logN01: (float) 光学フォノン散乱のキャリア濃度基準値の常用対数。
:returns: (tuple[float, float]) `(Eop_0, Aop)` のタプル。`Eop_0` はそのまま、`Aop` は計算された値。
"""
logN = log10(N)
return Eop_0, A1 / (1 + exp(-B1 * (logN - logN01))) + C1
[ドキュメント]
def Aop1(N, Eop_0, A1, B1, C1, logN01):
"""
光学フォノン散乱のパラメータ Aop を計算します（Eop_0 を含まない）。
詳細説明:
キャリア濃度 `N` に基づいて、光学フォノン散乱のパラメータ `Aop` を計算します。
`Aop` はフェルミ分布関数に似た形でキャリア濃度に依存します。
`Aop` 関数と異なり、`Eop_0` は返しません。
:param N: (float) キャリア濃度 (cm^-3)。
:param Eop_0: (float) (未使用) 光学フォノン散乱のエネルギーパラメータ。
:param A1: (float) 光学フォノン散乱の振幅パラメータ。
:param B1: (float) 光学フォノン散乱の濃度依存性パラメータ。
:param C1: (float) 光学フォノン散乱のオフセットパラメータ。
:param logN01: (float) 光学フォノン散乱のキャリア濃度基準値の常用対数。
:returns: (float) 計算された光学フォノン散乱のパラメータ `Aop`。
"""
logN = log10(N)
return A1 / (1 + exp(-B1 * (logN - logN01))) + C1
[ドキュメント]
def AT0(N, A2, B2, C2, logN02):
"""
音響フォノン散乱に関連するパラメータ AT0 を計算します。
詳細説明:
キャリア濃度 `N` に基づいて、音響フォノン散乱のパラメータ `AT0` を計算します。
これは、キャリア濃度に依存する飽和関数として表現されます。
:param N: (float) キャリア濃度 (cm^-3)。
:param A2: (float) 音響フォノン散乱の振幅パラメータ。
:param B2: (float) 音響フォノン散乱の濃度依存性パラメータ。
:param C2: (float) 音響フォノン散乱のオフセットパラメータ。
:param logN02: (float) 音響フォノン散乱のキャリア濃度基準値の常用対数。
:returns: (float) 計算された音響フォノン散乱のパラメータ `AT0`。
"""
logN = log10(N)
return A2 /(1.0 + exp(B2 * (logN - logN02))) + C2
[ドキュメント]
def AT32(N, A3, B3, C3, D3, F3, logN03):
"""
イオン化不純物散乱に関連するパラメータ AT32 を計算します。
詳細説明:
キャリア濃度 `N` に基づいて、イオン化不純物散乱のパラメータ `AT32` を計算します。
複数の項が組み合わされており、キャリア濃度に対する複雑な依存性を示します。
:param N: (float) キャリア濃度 (cm^-3)。
:param A3: (float) イオン化不純物散乱の第1項パラメータ。
:param B3: (float) イオン化不純物散乱の濃度依存性パラメータ。
:param C3: (float) イオン化不純物散乱の第2項振幅パラメータ。
:param D3: (float) イオン化不純物散乱の第2項濃度依存性パラメータ。
:param F3: (float) イオン化不純物散乱のオフセットパラメータ。
:param logN03: (float) イオン化不純物散乱のキャリア濃度基準値の常用対数。
:returns: (float) 計算されたイオン化不純物散乱のパラメータ `AT32`。
"""
logN = log10(N)
K1 = 1.0 + exp(B3 * (logN - logN03))
a1 = A3 * N * N / K1
K2 = 1.0 + exp(-B3 * (logN - logN03))
a2 = C3 / K2 / (1.0 + D3 * N) + F3
return a1 + a2
[ドキュメント]
def VB(N, A6, B6, C6, D6, logN06):
"""
障壁散乱（またはキャリア濃度依存の有効質量）に関連するパラメータ VB を計算します。
詳細説明:
キャリア濃度 `N` に基づいて、障壁散乱のパラメータ `VB` を計算します。
これはキャリア濃度の対数に線形に依存する項を含む関数として表現されます。
:param N: (float) キャリア濃度 (cm^-3)。
:param A6: (float) 障壁散乱の振幅パラメータ。
:param B6: (float) 障壁散乱の濃度依存性パラメータ。
:param C6: (float) 障壁散乱のオフセットパラメータ。
:param D6: (float) 障壁散乱の対数濃度依存性パラメータ。
:param logN06: (float) 障壁散乱のキャリア濃度基準値の常用対数。
:returns: (float) 計算された障壁散乱のパラメータ `VB`。
"""
logN = log10(N)
return A6 * (1.0 - D6 * logN) / (1.0 + exp(B6 * (logN - logN06))) + C6
[ドキュメント]
def cal_mu_from_params(T, _Eop, _Aop, _AT0, _AT32, _VB):
"""
各散乱パラメータからホール移動度を計算します。
詳細説明:
温度 `T` と各散乱メカニズムのパラメータ（`_Eop`, `_Aop`, `_AT0`, `_AT32`, `_VB`）を
用いて、ホール移動度 `mu` を計算します。
各散乱メカニズムによる移動度の逆数を合計し、それに温度依存の係数を乗じます。
:param T: (float) 温度 (K)。
:param _Eop: (float) 光学フォノン散乱のエネルギーパラメータ。
:param _Aop: (float) 光学フォノン散乱の振幅パラメータ。
:param _AT0: (float) 音響フォノン散乱のパラメータ。
:param _AT32: (float) イオン化不純物散乱のパラメータ。
:param _VB: (float) 障壁散乱のパラメータ。
:returns: (float) 計算されたホール移動度 `mu`。
"""
ekBT = e / kB / T
muinv = _Aop / (exp(_Eop * ekBT) - 1.0)
muinv += _AT0
muinv += _AT32 * pow(T, -1.5)
mu = 1.0 / muinv * exp(-_VB * ekBT)
return mu
[ドキュメント]
def cal_mulist(pk, T = None):
"""
温度のリストに対してホール移動度のリストを計算します。
詳細説明:
与えられたパラメータ `pk` と温度リスト `T` を使用して、
各温度におけるホール移動度を `cal_mu_from_params` 関数を呼び出して計算し、
そのリストを返します。
:param pk: (list[float]) 全散乱パラメータのリスト。
:param T: (list[float] | None) 温度のリスト (K)。デフォルトは `slef.T`。
:returns: (list[float]) 各温度に対応するホール移動度のリスト。
"""
if T is None:
T = slef.T
return [cal_mu_from_params(T[i], *pk) for i in range(len(T))]
[ドキュメント]
def pk2params(*pk):
"""
全パラメータリストを個別の散乱モデルパラメータ群に分割します。
詳細説明:
結合されたパラメータリスト `pk` を、光学フォノン散乱 (`Aop_params`)、
音響フォノン散乱 (`AT0_params`)、イオン化不純物散乱 (`AT32_params`)、
障壁散乱 (`VB_params`) の各モデルに対応する部分リストに分割します。
:param *pk: (float) 可変長引数として渡される全フィッティングパラメータ。
:returns: (tuple[list[float], list[float], list[float], list[float]]) 各散乱モデルのパラメータリストを含むタプル。
"""
return pk[0:5], pk[5:9], pk[9:15], pk[15:20]
[ドキュメント]
def cal_mu_from_TN(x_list, pk):
"""
温度とキャリア濃度からホール移動度を計算します。
詳細説明:
温度 `T` とキャリア濃度 `N`、および全フィッティングパラメータ `pk` を使用して、
各散乱メカニズム（光学フォノン、音響フォノン、イオン化不純物、障壁）によるパラメータを計算し、
それらを用いて最終的なホール移動度を算出します。
:param x_list: (list[float, float]) 温度 `T` とキャリア濃度 `N` を含むリスト。
:param pk: (list[float]) 全フィッティングパラメータのリスト。
:returns: (float) 計算されたホール移動度 `mu`。
"""
T, N = x_list
ekBT = e / kB / T
Aop_params, AT0_params, AT32_params, VB_params = pk2params(*pk)
_Aop, _Eop = Aop(N, *Aop_params)
_AT0 = AT0(N, *AT0_params)
_AT32 = AT32(N, *AT32_params)
_VB = VB(N, *VB_params)
#def VB(N, A6, B6, C6, D6, logN06):
# print("params=", _Aop, _Eop, _AT0, _AT32, _VB)
return cal_mu_from_params(T, _Eop, _Aop, _AT0, _AT32, _VB)
[ドキュメント]
def fit_mu_TN(app, fit):
"""
温度とキャリア濃度依存のホール移動度データに対する最小二乗フィッティングを実行します。
詳細説明:
Excelファイルから温度 `T`、キャリア濃度 `N`、ホール移動度 `mu` のデータを読み込み、
複数の散乱モデルを含む複合関数 `cal_mu_from_TN` を用いて、非線形最小二乗法でフィッティングを行います。
フィッティング結果はグラフ表示され、Excelファイルとパラメータファイルに保存されます。
また、指定された温度でキャリア濃度を固定するオプションもサポートします。
:param app: (tkApplication) アプリケーションオブジェクト。
:param fit: (tkFit_m) フィッティング処理を行う `tkFit_m` オブジェクト。
:returns: None
"""
cparams = app.cparams
fit = tkFit_m(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
# fit = tkFit_TN(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
# print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
print("")
print("Fitting to mu-T Hall data by nonlinear least-squares method")
print("mode : ", cparams.mode)
print("target : ", cparams.target)
print("infile : ", cparams.infile)
print("parameter file : ", cparams.parameterfile)
print("xlabel : ", cparams.xlabel)
print("ylabel : ", cparams.ylabel)
print("Nlabel : ", cparams.Nlabel)
print("method : ", cparams.method)
print("tol : ", cparams.tol)
print("nmaxiter: ", cparams.nmaxiter)
print("")
print(f"Read data from [{cparams.infile}]")
print(f"Extract data for {cparams.Tmin} <= T <= {cparams.Tmax}")
# cparams.T_label = cparams.xlabel
# cparams.mu_label = cparams.ylabel
fit.read_data(cparams.infile, xlabels = [cparams.xlabel, cparams.Nlabel], ylabel = cparams.mu_label,
xmins = [cparams.Tmin, None], xmaxs = [cparams.Tmax, None], usage = usage)
# fit.read_data(cparams.infile, cparams.T_label, cparams.mu_label, usage = usage)
# N_label, N_in = fit.datafile.FindDataArray(cparams.Nlabel, flag = 'i')
sample_label, sample_in = fit.datafile.FindDataArray(cparams.sample_label, flag = 'i')
# formula_label, formula_in = fit.datafile.FindDataArray(cparams.formula_label, flag = 'i')
fit.sample = sample_in
fit.sample_label = sample_label
# x_list[0]: T
# x_list[1]: N
print("ndata_all:", fit.ndata_all)
print("ndata :", fit.ndata)
print("xlabel[0](T): ", fit.xlabels[0])
print("xlabel[1](N): ", fit.xlabels[1])
print("ylabel : ", fit.ylabel)
print(f"{'sample':10}: {fit.xlabels[0]:10} {fit.xlabels[1]:10} {fit.ylabel:10}")
for i in range(fit.ndata):
print(f"{sample_in[i]:10}: {fit.x_list[0][i]:10.4g} {fit.x_list[1][i]:10.4g} {fit.y[i]:10.4g}")
# Calculate Nbase at Tbase for different samples
data_set = {}
for i in range(fit.ndata):
index = fit.included_index[i]
s = fit.sample[index]
if data_set.get(s, None) is None:
data_set[s] = [[fit.x_list[1][i]], [fit.x_list[1][i]]]
else:
data_set[s][0].append(fit.x_list[0][i])
data_set[s][1].append(fit.x_list[1][i])
fit.Nbases = []
Nbases_dict = {}
for i in range(fit.ndata):
index = fit.included_index[i]
s = fit.sample[index]
Nbase = np.interp([cparams.Tbase], data_set[s][0], data_set[s][1])[0]
Nbases_dict[s] = Nbase
fit.Nbases.append(Nbase)
# print("s=", sample_in)
# print("Nbase=", fit.Nbases)
print("")
if cparams.fix_N:
fit.fix_N = True
fit.Tbase = cparams.Tbase
fit.x_list[1] = fit.Nbases
for s in Nbases_dict.keys():
print(f"sample [{s}]: N is fixed to N={Nbases_dict[s]:10.4g} cm-3 at T={cparams.Tbase} K")
else:
print(f"Nbase={cparams.Nbase:10.4g} cm-3 from T={cparams.Tbase} K")
print(f"N values are taken from T-dependent Hall data")
#=============================================
# Build pk
#=============================================
fit.func = cal_mu_from_TN
fit.varname = app.Aop_params.varname.copy()
fit.unit = app.Aop_params.varunit.copy()
fit.pk = app.Aop_params.value.copy()
fit.optid = app.Aop_params.optid.copy()
fit.varname += app.AT0_params.varname
fit.unit += app.AT0_params.varunit
fit.pk += app.AT0_params.value
fit.optid += app.AT0_params.optid
fit.varname += app.AT32_params.varname
fit.unit += app.AT32_params.varunit
fit.pk += app.AT32_params.value
fit.optid += app.AT32_params.optid
fit.varname += app.VB_params.varname
fit.unit += app.VB_params.varunit
fit.pk += app.VB_params.value
fit.optid += app.VB_params.optid
# print("pk=", len(fit.pk), fit.pk)
# print("optid=", len(fit.optid), fit.optid)
print("")
print("# of variables:", len(fit.pk))
fit.print_variables()
print("")
yini = fit.cal_ylist(fit.pk, fit.x_list)
fit.print_data(heading = "Initial functions", yini = yini)
#========================================
# plot
#========================================
fig, axes = plt.subplots(1, 1, figsize = cparams.figsize)
axes = [axes]
# axes = axes.flatten()
fit.initial_plot(axes[0], yini = yini, fontsize = cparams.fontsize)
#========================================
# Optimize
#========================================
pfin, ffin, success, res = fit.minimize(cparams.method)
if success:
print("")
print(f"Converged at iteration: {res.nit}")
else:
print(f"Function did not converge")
print(res)
pk = pfin
n = len(pk)
print("Final parameters")
for i in range(n):
print(f" {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}")
print(f" f={ffin:12.6g}")
app.Aop_params.value, app.AT0_params.value, app.AT32_params.value, app.VB_params.value = pk2params(*pfin)
#========================================
# Final result
#========================================
yfin = fit.cal_ylist(pfin)
fit.print_data(heading = "Final functions", yini = yini, yfin = yfin)
outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-mu_TN-fit.xlsx"])
print("")
print(f"Save results to [{outfile}]")
fit.to_excel(outfile, ['N(cm-3)', 'T(K)', 'mu', 'mu(ini)', 'mu(fin)'], [fit.x_list[1], fit.x_list[0], fit.y, yini, yfin])
# Delete parameters in cparams so as not to save to Preferences section
for key in ["Eop_0", "A1", "B1", "C1", "logN01", "A2", "B2", "C2", "logN02", "A3", "B3", "C3", "D3", "F3",
"logN03", "A6", "B6", "C6", "D6", "logN06", "Aop", "AT0", "AT32", "VB"]:
if cparams.get(key, None) is not None:
cparams.__dict__[key] = None
print("")
print(f"Save parameters to {cparams.parameterfile}")
# cparams.print_parameters()
cparams.save_parameters(cparams.parameterfile, section = 'Preferences', sort_by_keys = False, update_commandline = False, IsPrint = False)
fit.save_parameters (cparams.parameterfile, section = 'Aop', heading = "Save Aop_parameters",
keys = app.Aop_params.varname, values = app.Aop_params.value, optid = app.Aop_params.optid)
fit.save_parameters (cparams.parameterfile, section = 'AT0', heading = "Save AT0_parameters",
keys = app.AT0_params.varname, values = app.AT0_params.value, optid = app.AT0_params.optid)
fit.save_parameters (cparams.parameterfile, section = 'AT32', heading = "Save AT32_parameters",
keys = app.AT32_params.varname, values = app.AT32_params.value, optid = app.AT32_params.optid)
fit.save_parameters (cparams.parameterfile, section = 'VB', heading = "Save VB_parameters",
keys = app.VB_params.varname, values = app.VB_params.value, optid = app.VB_params.optid)
print("")
print(f"Scattering parameters at Tbase={cparams.Tbase:10.4g} K")
print("a=", app.Aop_params.value)
for s in Nbases_dict.keys():
Nbase = Nbases_dict[s]
mu_params = tkParams()
mu_params.Eop, mu_params.Aop = Aop(Nbase, *app.Aop_params.value)
mu_params.AT0 = AT0 (Nbase, *app.AT0_params.value)
mu_params.AT32 = AT32(Nbase, *app.AT32_params.value)
mu_params.VB = VB (Nbase, *app.VB_params.value)
print(f" sample [{s}]: N is fixed to N={Nbases_dict[s]:10.4g} cm-3 at T={cparams.Tbase} K")
print(f" Aop ={mu_params.Aop:10.4g}")
print(f" Eop ={mu_params.Eop:10.4g}")
print(f" AT0 ={mu_params.AT0:10.4g}")
print(f" AT32={mu_params.AT32:10.4g}")
print(f" VB ={mu_params.VB:10.4g}")
mu_params.save_parameters(cparams.parameterfile, section = s, sort_by_keys = False, update_commandline = False, IsPrint = False)
fit.finalize_plot(yfin, iter = res.nit, fmin = ffin)
app.terminate("", pause = True)
[ドキュメント]
def fit_mu_T(app, fit = None):
"""
特定のサンプルの温度依存ホール移動度データに対する最小二乗フィッティングを実行します。
詳細説明:
Excelファイルからデータを読み込み、特定のサンプルについて温度 `T` とホール移動度 `mu` のデータを抽出します。
その際、基準温度 `Tbase` におけるキャリア濃度 `Nbase` を補間し、この `Nbase` を用いて
各散乱メカニズムのパラメータを固定値またはフィッティング対象として、温度依存性のホール移動度をフィッティングします。
結果はグラフ表示され、Excelファイルとパラメータファイルに保存されます。
:param app: (tkApplication) アプリケーションオブジェクト。
:param fit: (tkFit | None) フィッティング処理を行う `tkFit` オブジェクト。デフォルトは `None`。
:returns: None
"""
cparams = app.cparams
cparams.target_sample = '02A-01'
mu_params = tkParams()
print("")
print("Fitting to mu-T Hall data by nonlinear least-squares method")
print("mode : ", cparams.mode)
print("target : ", cparams.target)
print("infile : ", cparams.infile)
print("parameter file : ", cparams.parameterfile)
print("xlabel : ", cparams.xlabel)
print("ylabel : ", cparams.ylabel)
print("method : ", cparams.method)
print("tol : ", cparams.tol)
print("nmaxiter: ", cparams.nmaxiter)
fit = tkFit(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
print("")
print(f"Read data from [{cparams.infile}]")
fit.read_data(cparams.infile, cparams.T_label, cparams.mu_label, usage = usage)
sample_label, sample_in = fit.datafile.FindDataArray(cparams.sample_label, flag = 'i')
formula_label, formula_in = fit.datafile.FindDataArray(cparams.formula_label, flag = 'i')
N_label, N_in = fit.datafile.FindDataArray(cparams.Nlabel, flag = 'i')
print("")
print(f"Extract data for sample [{cparams.target_sample}] and {cparams.Tmin} <= T <= {cparams.Tmax}")
fit.T = []
fit.mu = []
fit.N = []
for i in range(len(fit.x)):
if sample_in[i] == cparams.target_sample and cparams.Tmin <= fit.x[i] <= cparams.Tmax:
fit.T.append(fit.x[i])
fit.mu.append(fit.y[i])
fit.N.append(N_in[i])
fit.Tlabel = fit.xlabel
fit.mulabel = fit.ylabel
fit.Nlabel = N_label
fit.x = fit.T
fit.y = fit.mu
ndata = len(fit.T)
print("ndata=", ndata)
print(f"{fit.Tlabel:10} {fit.mulabel:10} {fit.Nlabel:10}")
for i in range(ndata):
print(f"{fit.T[i]:10.4g} {fit.mu[i]:10.4g} {fit.N[i]:10.4g}")
cparams.Nbase = np.interp([cparams.Tbase], fit.T, fit.N)[0]
mu_params.Eop, mu_params.Aop = Aop(cparams.Nbase, *app.Aop_params.value)
mu_params.AT0 = AT0 (cparams.Nbase, *app.AT0_params.value)
mu_params.AT32 = AT32(cparams.Nbase, *app.AT32_params.value)
mu_params.VB = VB (cparams.Nbase, *app.VB_params.value)
fit.varname = [ "Eop", "Aop", "AT0", "AT32", "VB"]
fit.unit = [ "", "", "", "", "eV"]
fit.pk = [mu_params.Eop, mu_params.Aop, mu_params.AT0, mu_params.AT32, mu_params.VB]
fit.optid = [ 0, 1, 1, 1, 1]
fit.func = lambda T, *pk: cal_mu_from_params(T, *pk)
print("")
print(f"Scattering parameters at Nbase={cparams.Nbase:10.4g} cm-3")
fit.print_variables()
optpk = fit.extract_parameters()
yini = cal_mulist(fit.pk, fit.T)
print("")
print("Initial guess")
print(f"{'i':4} {'T':10} {'mu':12} {'mu(cal)':12}")
for i in range(len(fit.T)):
print(f"{i:4} {fit.T[i]:10.4g} {fit.mu[i]:12.4g} {yini[i]:12.4g}")
#========================================
# plot
#========================================
fig, axes = plt.subplots(1, 1, figsize = cparams.figsize)
axes = [axes]
# axes = axes.flatten()
fit.initial_plot(axes[0], yini = yini, fig = fig, fontsize = cparams.fontsize)
# fit.initial_plot(data_axis = axes[0], error_axis = axes[1], yini = yini, label_ini = 'initial',
# fmin = fini, fig = fig,
# fontsize = cparams.fontsize)
#========================================
# Optimize
#========================================
pfin, ffin, success, res = fit.minimize(cparams.method)
if success:
# print("")
# print(res)
print("")
print(f"Converged at iteration: {res.nit}")
else:
print(f"Function did not converge")
print(res)
pk = fit.recover_parameters(optpk, set_member = True)
print("Final parameters")
n = len(pfin)
for i in range(n):
print(f" {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}")
print(f" f={ffin:12.6g}")
#========================================
# Final result
#========================================
yfin = fit.cal_ylist(pfin)
print("")
print("Final result")
print(f"{'i':4} {'T':10} {'mu':12} {'mu(ini)':12} {'mu(fin)':12}")
for i in range(len(fit.x)):
print(f"{i:4} {fit.T[i]:10.4g} {fit.mu[i]:12.4g} {yini[i]:12.4g} {yfin[i]:12.4g}")
outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-mu-fit.xlsx"])
df = pd.DataFrame(np.array([fit.T, fit.mu, yini, yfin]).T, columns = ['T(K)', 'mu', 'mu(ini)', 'mu(fin)'])
print("")
print(f"Save results to [{outfile}]")
df.to_excel(outfile, index = False)
print("")
print(f"Save parameters to {cparams.parameterfile}")
cparams.save_parameters (cparams.parameterfile, section = 'Preferences', sort_by_keys = False, update_commandline = False, IsPrint = False)
mu_params.save_parameters (cparams.parameterfile, section = 'mu', sort_by_keys = False, update_commandline = False, IsPrint = False)
fit.fit_data[0].set_data(fit.x, yfin)
plt.tight_layout()
plt.subplots_adjust(top = 0.90, bottom = 0.10)
plt.pause(0.0001)
app.terminate("", usage = usage, pause = True)
[ドキュメント]
def fit_parameters(app, fit = None):
"""
単一の散乱モデルパラメータに対する最小二乗フィッティングを実行します。
詳細説明:
ユーザーが指定したターゲット（例: `Aop`, `A_T0`, `A_T32`, `VB`）に基づいて、
ExcelファイルからXYデータを読み込み、対応する単一の散乱モデル関数を用いて
非線形最小二乗法でフィッティングを行います。
フィッティング結果はグラフ表示され、Excelファイルとパラメータファイルに保存されます。
:param app: (tkApplication) アプリケーションオブジェクト。
:param fit: (tkFit | None) フィッティング処理を行う `tkFit` オブジェクト。デフォルトは `None`。
:returns: None
"""
cparams = app.cparams
print("")
print("Fitting to mu-T Hall data by nonlinear least-squares method")
print("mode : ", cparams.mode)
print("target : ", cparams.target)
print("infile : ", cparams.infile)
print("xlabel : ", cparams.xlabel)
print("ylabel : ", cparams.ylabel)
print("method : ", cparams.method)
fit = tkFit(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
print("tol : ", fit.tol)
print("nmaxiter : ", fit.nmaxiter)
print("print_interval: ", fit.print_interval)
print("plot_interval : ", fit.plot_interval)
if cparams.target == 'Aop':
fit.varname = app.Aop_params.varname
fit.unit = app.Aop_params.varunit
fit.pk = app.Aop_params.value
fit.optid = app.Aop_params.optid
fit.func = Aop1
elif cparams.target == 'A_T0':
fit.varname = app.AT0_params.varname
fit.unit = app.AT0_params.varunit
fit.pk = app.AT0_params.value
fit.optid = app.AT0_params.optid
fit.func = AT0
elif cparams.target == 'A_T32':
fit.varname = app.AT32_params.varname
fit.unit = app.AT32_params.varunit
fit.pk = app.AT32_params.value
fit.optid = app.AT32_params.optid
fit.func = AT32
elif cparams.target == 'VB':
fit.varname = app.VB_params.varname
fit.unit = app.VB_params.varunit
fit.pk = app.VB_params.value
fit.optid = app.VB_params.optid
fit.func = VB
print("")
print(f"Read data from [{cparams.infile}]")
fit.read_data(cparams.infile, cparams.xlabel, cparams.ylabel, usage = usage)
ndata = len(fit.x)
print("x=", fit.xlabel, fit.x)
print("y=", fit.ylabel, fit.y)
print("ndata=", ndata)
print("pk=", fit.pk)
yini = fit.cal_ylist(fit.pk, fit.x)
print("")
print("Initial guess")
print(f"{'i':4} {cparams.xlabel:10} {cparams.ylabel:12} {cparams.xlabel+'(cal)':12}")
for i in range(len(fit.x)):
print(f"{i:4} {fit.x[i]:10.4g} {fit.y[i]:12.4g} {yini[i]:12.4g}")
# plot
fig, axes = plt.subplots(1, 1, figsize = cparams.figsize)
axes = [axes]
# axes = axes.flatten()
fit.initial_plot(data_axis = axes[0], yini = yini, fig = fig, fontsize = cparams.fontsize)
# fit.initial_plot(data_axis = axes[0], error_axis = axes[1], yini = yini, label_ini = 'initial',
# fmin = fini, fig = fig,
# fontsize = cparams.fontsize)
pfin, ffin, success, res = fit.minimize(cparams.method, cparams.tol, cparams.nmaxiter)
if success:
# print("")
# print(res)
print("")
print(f"Converged at iteration: {res.nit}")
else:
print(f"Function did not converge")
print(res)
print("Final parameters")
n = len(pfin)
for i in range(n):
print(f" {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}")
print(f" f={ffin:12.6g}")
yfin = fit.cal_ylist(pfin)
print("")
print("Final result")
print(f"{'i':4} {fit.xlabel:10} {fit.ylabel:12} {fit.ylabel + '(ini)':12} {'y(fin)':12}")
for i in range(len(fit.x)):
print(f"{i:4} {fit.x[i]:10.4g} {fit.y[i]:12.4g} {yini[i]:12.4g} {yfin[i]:12.4g}")
outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-{target}-fit.xlsx"], ext_dict = { "target": cparams.target } )
df = pd.DataFrame(np.array([fit.x, fit.y, yini, yfin]).T,
columns = [cparams.xlabel, cparams.ylabel, f"{cparams.ylabel}(ini)", f"{cparams.ylabel}(fin)"])
print("")
print(f"Save results to [{outfile}]")
df.to_excel(outfile, index = False)
print("")
print(f"Save parameters to {cparams.parameterfile}")
# cparams.print_parameters()
cparams.save_parameters(cparams.parameterfile, section = 'Preferences', sort_by_keys = False, update_commandline = False, IsPrint = False)
fit.save_parameters (cparams.parameterfile, section = 'Aop', heading = "Save Aop_parameters",
keys = app.Aop_params.varname, values = app.Aop_params.value, optid = app.Aop_params.optid)
fit.save_parameters (cparams.parameterfile, section = 'AT0', heading = "Save AT0_parameters",
keys = app.AT0_params.varname, values = app.AT0_params.value, optid = app.AT0_params.optid)
fit.save_parameters (cparams.parameterfile, section = 'AT32', heading = "Save AT32_parameters",
keys = app.AT32_params.varname, values = app.AT32_params.value, optid = app.AT32_params.optid)
fit.save_parameters (cparams.parameterfile, section = 'VB', heading = "Save VB_parameters",
keys = app.VB_params.varname, values = app.VB_params.value, optid = app.VB_params.optid)
fit.fit_data[0].set_data(fit.x, yfin)
plt.tight_layout()
plt.subplots_adjust(top = 0.90, bottom = 0.10)
plt.pause(0.0001)
app.terminate("\nPress ENTER to exit>>", usage = usage, pause = True)
[ドキュメント]
def main():
"""
スクリプトのメインエントリポイントです。
詳細説明:
アプリケーションと共通パラメータを初期化し、コマンドライン引数でパラメータを更新します。
ログファイルの準備と既存パラメータファイルの読み込みを行った後、
`cparams.mode` と `cparams.target` に応じて適切なフィッティング関数を呼び出します。
サポートされるモードは 'init', 'fit', 'sim' で、'fit' モードではさらに
'mu_T', 'mu_TN', 'mu_T_Nfix', 'Aop', 'A_T0', 'A_T32', 'VB' のターゲットが選択可能です。
:returns: None
"""
app, cparams = initialize()
update_vars(app, cparams)
cparams.logfile = app.replace_path(cparams.infile)
cparams.outfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-out.xlsx"])
print(f"Open logfile [{cparams.logfile}]")
app.redict(targets = ["stdout", cparams.logfile], mode = 'w')
cparams.parameterfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}.in"])
cparams.parameterbkfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-back.in"])
print("")
print("infile : {}".format(cparams.infile))
print("parameter file : {}".format(cparams.parameterfile))
print("parameter backup file: {}".format(cparams.parameterbkfile))
print("mode : {}".format(cparams.mode))
print("target : {}".format(cparams.target))
print("")
print("Read [{}]".format(cparams.parameterfile))
cparams.read_parameters(cparams.parameterfile, section = "Preferences")
fit = tkFit(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
# fit = tkFit_TN(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
# print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
fit.read_parameters(cparams.parameterfile, section = "Aop",
keys = app.Aop_params.varname, values = app.Aop_params.value, optid = app.Aop_params.optid)
app.Aop_params.print_parameters(heading = "Aop")
fit.read_parameters(cparams.parameterfile, section = "AT0",
keys = app.Aop_params.varname, values = app.Aop_params.value, optid = app.Aop_params.optid)
app.AT0_params.print_parameters(heading = "AT0")
fit.read_parameters(cparams.parameterfile, section = "AT32",
keys = app.Aop_params.varname, values = app.Aop_params.value, optid = app.Aop_params.optid)
app.AT32_params.print_parameters(heading = "AT32")
fit.read_parameters(cparams.parameterfile, section = "VB",
keys = app.Aop_params.varname, values = app.Aop_params.value, optid = app.Aop_params.optid)
app.VB_params.print_parameters(heading = "VB")
# check args again to use given parameters
update_vars(app, cparams)
cparams.print_parameters(heading = "\ncparams:")
app.Aop_params.print_parameters(heading = "\nAop_params:")
app.AT0_params.print_parameters(heading = "\nAT0_params:")
app.AT32_params.print_parameters(heading = "\nAT32_params:")
app.VB_params.print_parameters(heading = "\nVB_params:")
if cparams.mode == 'init':
# init関数は存在しないが、既存のコードを変更しないため、このまま。
# 本来はinit関数を定義するか、ここを削除すべき。
init(app)
elif cparams.mode == 'fit':
if cparams.target == 'mu_T':
fit_mu_T(app, fit)
elif cparams.target == 'mu_TN':
fit_mu_TN(app, fit)
elif cparams.target == 'mu_T_Nfix':
cparams.fix_N = True
fit_mu_TN(app, fit)
elif cparams.target == 'Aop' or cparams.target == 'A_T0' or cparams.target == 'A_T32' or cparams.target == 'VB':
fit_parameters(app, fit)
else:
app.terminate("Error in main: Invalide target [{}]".format(cparams.target), usage = usage)
elif cparams.mode == 'sim':
# sim関数は存在しないが、既存のコードを変更しないため、このまま。
# 本来はsim関数を定義するか、ここを削除すべき。
sim(app)
else:
app.terminate("Error in main: Invalide mode [{}]".format(cparams.mode), usage = usage)
app.terminate("", usage = usage, pause = True)
if __name__ == "__main__":
main()