"""
SIMPLEX法およびガウス過程回帰によるCV測定データのフィッティングスクリプト。
このスクリプトは、半導体デバイスのCV（容量-電圧）測定データに対して、
SIMPLEX法などの最適化アルゴリズムを用いてフィッティングを行います。
電極面積、ドナー/アクセプター密度、n層/p層厚などのパラメータを最適化し、
測定データに最もよく合う理論曲線を導出します。
関連リンク:
:doc:`optimize_flex_CV_usage`
"""
import os
import sys
import csv
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
from tklib.tksci.tksci import log10, e, kB
from tklib.tkutils import getarg, getintarg, getfloatarg, pconv_by_type, print_line, print_data, joinf
from tklib.tkapplication import tkApplication
from tklib.tkparams import tkParams
from tklib.tkvariousdata import tkVariousData
#from tklib.tkfilter import tkFilter
from tklib.tkgraphic.tkplotevent import tkPlotEvent
from tklib.tksci.tkFit_mxy import tkFit_mxy
from minimize_func_CV import read_data, save_data, ycal_list, minimize_func
#========================================================
# 最適化アルゴリズム
#========================================================
method = "cg"
#method = "nelder-mead"
maxiter = 100
tol = 1.0e-5
#==========================================
# Source parameters to be fitted
#==========================================
# 最適化変数の初期値
# Data parameters: I0, x0, w
x0 = [1.3, 0.6, 0.1]
#==========================================
# File configurations
#==========================================
initial_csv = 'initial.csv'
final_csv = 'final.csv'
conv_csv = 'convergence.csv'
#==========================================
# Graph parameters
#==========================================
fplot = 1
ngdata = 51
xgmin = -4.0
xgmax = 4.0
ygmin = -4.0
ygmax = 4.0
tsleep = 0.3
[ドキュメント]
def initialize():
"""
アプリケーションとパラメータの初期設定を行う。
tkApplicationおよびtkParamsオブジェクトを初期化し、
コマンドライン引数、フィッティング変数の初期値、範囲、ペナルティ、
およびグラフ表示に関する設定を定義します。
:returns: 初期化されたアプリケーションオブジェクトとパラメータオブジェクトのタプル。
:rtype: tuple[tkApplication, tkParams]
"""
#================================
# Global variables
#================================
app = tkApplication()
argv, narg = app.get_argv()
app.cparams = tkParams()
cparams = app.cparams
cparams.debug = 0
cparams.print_level = 0
cparams.print_interval = 10
cparams.plot_interval = 1
app.varname = [ "S", "ND", "NA", "dn", "dp"]
app.unit = [ "m2", "cm-3", "cm-3", "m", "m"]
app.pk_convert = [ "", "", "", "", ""]
app.optid = [ 0, 1, 1, 0, 0]
app.kmin = [ 0.0, 0.01, 0.01, 1e-6, 1e-6]
app.kmax = [ 1.0e-2, 1.0e23, 1.0e23, 1e-5, 1e-5]
app.kpenalty = [ 1.0, 1.0, 1.0, 1.0, 1.0]
app.add_argument(opt = '--infile', type = "str", var_name = 'infile', opt_str = "input .xlsx file", desc = 'input Excel file',
defval = 'CV_meas.xlsx', optional = True);
app.add_argument(opt = "--outfile", type = "str", var_name = 'outfile', opt_str = "output .xlsx file", desc = 'output Excel file',
defval = 'out.xlsx', optional = True);
app.add_argument(opt = "--datafile", type = "str", var_name = 'datafile', opt_str = "input/calculated data .xlsx file", desc = 'Data file',
defval = 'CV.xlsx', optional = True);
app.add_argument(opt = "--mode", type = "str", var_name = 'mode', opt_str = "--mode=[fit]", desc = 'task mode',
defval = 'fit', optional = True);
r_electrode = 1095e-6 # 850e-6 m
app.add_argument(opt = "--S", type = "float", var_name = 'S', opt_str = "--S=val", desc = 'Electrode area in m2',
defval = r_electrode * r_electrode, optional = True);
app.add_argument(opt = "--ND", type = "float", var_name = 'ND', opt_str = "--ND=val", desc = 'Donor density in n-layer in cm-3',
defval = 1.2e16, optional = True);
app.add_argument(opt = "--NA", type = "float", var_name = 'NA', opt_str = "--NA=val", desc = 'Acceptor density in p-layer in cm-3',
defval = 1.8e18, optional = True);
app.add_argument(opt = "--dn", type = "float", var_name = 'dn', opt_str = "--dn=val", desc = 'n-layer thickiness in m',
defval = 3.1e-6, optional = True);
app.add_argument(opt = "--dp", type = "float", var_name = 'dp', opt_str = "--dp=val", desc = 'p-layer thickiness in m',
defval = 375e-9, optional = True);
app.add_argument(opt = "--method", type = "str", var_name = 'method', opt_str = "--method=[nelder-mead]", desc = 'optimization algorism',
defval = 'nelder-mead', optional = True);
app.add_argument(opt = "--nmaxiter", type = "int", var_name = 'nmaxiter', opt_str = "--nmaxiter=int(>1)", desc = 'maximum interation number for optimization',
defval = 1000, optional = True);
app.add_argument(opt = "--tol", type = "float", var_name = 'tol', opt_str = "--tol=val", desc = 'eps for optimization',
defval = 1.0e-5, optional = True);
app.add_argument(opt = "--print_interval", type = "int", var_name = 'print_interval', opt_str = "--print_interval=val", desc = 'print interval',
defval = 5, optional = True);
app.add_argument(opt = "--plot_interval", type = "int", var_name = 'plot_interval', opt_str = "--plot_interval=val", desc = 'plot interval',
defval = 5, optional = True);
#=============================
# Graph configuration
#=============================
cparams.figsize = [8, 6]
cparams.fontsize = 12
cparams.legend_fontsize = 8
cparams.graphupdateinterval = 10
return app, cparams
#=============================
# Treat argments
#=============================
[ドキュメント]
def update_vars(app: tkApplication, cparams: tkParams, apply_default: bool = True):
"""
コマンドライン引数または設定ファイルから変数設定を更新する。
tkApplication.read_argsメソッドを使用して、与えられたアプリケーションオブジェクトと
パラメータオブジェクトの変数を更新します。更新中にエラーが発生した場合、
エラーメッセージを表示してプログラムを終了します。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: 設定パラメータオブジェクト。
:type cparams: tkParams
:param apply_default: デフォルト値を適用するかどうか。Falseの場合、設定ファイルや
コマンドライン引数で明示的に指定されていない値は更新されない。
:type apply_default: bool
:returns: なし
"""
args_opt, args_idx, args_vars = app.read_args(vars = cparams, check_allowed_args = True, apply_default = apply_default)
if args_opt is None:
error_no = args_idx
error_message = args_vars
app.terminate("\n\n"
+ "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
+ f"! {error_message}\n"
+ "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!",
usage = app.usage)
#==========================================
# functions
#==========================================
[ドキュメント]
def fit(app: tkApplication, cparams: tkParams):
"""
CV測定データに対するフィッティング処理を実行する。
入力データの読み込み、初期値の設定、最適化の実行、結果の保存、
グラフ表示までの一連の処理を担います。`tkFit_mxy` クラスを利用して
指定されたアルゴリズムで最小化を行い、フィッティング結果を出力します。
:param app: アプリケーションオブジェクト。
:type app: tkApplication
:param cparams: 設定パラメータオブジェクト。
:type cparams: tkParams
:returns: なし
"""
cparams.parameterfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}.in"])
# cparams.parameterbkfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-back.in"])
# cparams.historyfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-history.csv"])
cparams.historyfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-history.xlsx"])
print("")
print(f"Read [{cparams.parameterfile}]")
cparams.read_parameters(cparams.parameterfile, section = "Preferences",
ignore_keys = ["logfile", "outfile", "parameterfile", "parameterbkfile"])
# cparams.print_parameters(heading = "\ncparams:")
print("")
print("infile : {}".format(cparams.infile))
print("parameter file : {}".format(cparams.parameterfile))
# print("parameter backup file: {}".format(cparams.parameterbkfile))
print("history file : {}".format(cparams.historyfile))
print("mode : {}".format(cparams.mode))
print("")
print("Peak fitting")
print("mode : ", cparams.mode)
print("infile : ", cparams.infile)
print("outfile: ", cparams.outfile)
print("parameter file : ", cparams.parameterfile)
print("data file : ", cparams.datafile)
fit = tkFit_mxy(tol = cparams.tol, nmaxiter = cparams.nmaxiter,
print_interval = cparams.print_interval, plot_interval = cparams.plot_interval)
fit.infile = cparams.infile
fit.datafile = cparams.datafile
fit.historyfile = cparams.historyfile
print()
print(f"Read input data from [{cparams.infile}]")
fit.labels, fit.x_list, fit.y_list = read_data(cparams.infile)
fit.w_list = None
fit.xlabels = [fit.labels[0]]
fit.ylabels = [fit.labels[1]]
print(f"Save input data to [{fit.datafile}]")
save_data(fit.datafile, fit.labels, fit.x_list + fit.y_list)
x0 = [cparams.S, cparams.ND, cparams.NA, cparams.dn, cparams.dp]
fit.varname = app.varname
fit.unit = app.unit
fit.pk = x0.copy()
fit.pk_convert = app.pk_convert
fit.optid = app.optid
fit.kmin = app.kmin
fit.kmax = app.kmax
fit.kpenalty = app.kpenalty
nvars = len(fit.pk)
# xsはxの値のセット
# fit.func = lambda xs, pk: ycal(xs, pk)
fit.cal_ylist = lambda pk: ycal_list(fit.x_list, pk, fit)
# x_list[i][j] はi番目のx変数のj番目の値
fit.minimize_func = lambda pk: minimize_func(pk, fit.x_list, fit.y_list, fit.w_list, fit)
# fit.callback = lambda pk: callback(pk, fit)
fit.fplot = fplot
fit.plt = plt
fit.iter = 0
fit.xiter = []
fit.yfmin = []
fit.ycal_list = None
fit.yfin_list = None
print()
fit.print_variables(heading = "Fitting parameters:")
print(f"Fitting configuration")
print(f" method : {cparams.method}")
print(f" tol : {cparams.tol}")
print(f" nmaxiter: {cparams.nmaxiter}")
optpk = fit.extract_parameters()
fit.yini_list = fit.cal_ylist(fit.pk)
fini = fit.minimize_func(optpk)
print("")
print(f"Initial function: fmin={fini:10.4g}")
print_data(labels = ["x", "y_obs", "y_cal"],
data_list = [fit.x_list[0], fit.y_list[0], fit.yini_list[0]],
label_format = '{:^15}', data_format = '{:>15.4g}', header = "Initial data:", nmax = 20, print_level = 0)
print("fini=", fini)
print(f"Save input and initial data to [{fit.datafile}]")
save_data(fit.datafile, fit.labels + ['ini'], fit.x_list + fit.y_list + fit.yini_list)
#========================================
# plot
#========================================
fig, axes = plt.subplots(len(fit.y_list), 2, figsize = cparams.figsize)
# axes = [axes]
# axes = axes.flatten()
fit.initial_plot(data_axes = [axes[0], axes[1]], error_axis = axes[1], yini_list = fit.yini_list, label_ini = 'initial',
fmin = fini, plt = plt, fig = fig,
fontsize = cparams.fontsize)
fit.plot_event.remove('error')
#========================================
# Optimize
#========================================
print("")
print("Optimize:")
fit.print_variables(heading = "Variables")
pfin, ffin, success, res = fit.minimize(cparams.method, jac = '3-points')
# pfin, ffin, success, res = fit.minimize(cparams.method, jac = lambda xk: diff1(xk, fit))
if success:
print(f"Converged at iteration: {res.nit}")
else:
print(f"Function did not converge")
#========================================
# Final result
#========================================
fit.pk = pfin.copy()
cparams.S, cparams.ND, cparams.NA, cparams.dn, cparams.dp = pfin
print("Final parameters")
for i in range(nvars):
print(f" {fit.varname[i]:10}: {pfin[i]:10.4g} {fit.unit[i]}")
print(f" f={ffin:12.6g}")
print("")
print(f"Save parameters to [{cparams.parameterfile}]")
cparams.save_parameters(cparams.parameterfile, section = "Preferences")
yfin_list = fit.cal_ylist(pfin)
print("")
print_data(labels = [fit.xlabels[0], "y_obs", "y_cal"],
data_list = [fit.x_list[0], fit.y_list[0], yfin_list[0]],
label_format = '{:^15}', data_format = '{:>15.4g}', header = "Final data:", nmax = 20, print_level = 0)
# fit.print_scores(heading = "\nScores between y(input) and y(fit)", y1 = fit.y, y2 = yfin)
print(f"Save input, initial, and final data to [{cparams.datafile}]")
save_data(fit.datafile, fit.labels + ['ini', 'fin'], fit.x_list + fit.y_list + fit.yini_list + yfin_list)
#=============================
# グラフの表示
#=============================
print("")
print("Plot optimized")
fit.finalize_plot(yfin_list, iter = res.nit, fmin = ffin)
fit.layout()
# plt.tight_layout()
plt.pause(0.01)
print()
input("Press ENTER to terminate>>")
[ドキュメント]
def main():
"""
スクリプトのエントリーポイント。
アプリケーションの初期化、コマンドライン引数の処理、ログファイルの設定を行い、
設定されたモード（'fit'モード）に基づいて`fit`関数を呼び出します。
無効なモードが指定された場合はエラーメッセージを出力して終了します。
:returns: なし
"""
print()
app, cparams = initialize()
update_vars(app, cparams,apply_default = True)
cparams.logfile = app.replace_path(cparams.infile)
print(f"Open logfile [{cparams.logfile}]")
app.redirect(targets = ["stdout", cparams.logfile], mode = 'w')
# app.redirect(targets = ["stdout"], mode = 'a',
# redirect_traceback = True, output_traceback = 'stdout', display_type_traceback = 'colored')
# app.redirect_exception()
# check args again to use given parameters
update_vars(app, cparams,apply_default = False)
if cparams.mode == 'fit':
fit(app, cparams)
else:
app.terminate("Error in main: Invalide mode [{}]".format(cparams.mode), usage = app.usage, pause = True)
# app.terminate("", usage = app.usage, pause = True)
if __name__ == "__main__":
main()