"""
lsq_func_tkFit.py
非線形最小二乗法を用いてデータに対する関数フィッティングを実行するスクリプトです。
指定された入力データと関数に基づき、最適化手法（例: BFGS）を用いてパラメータをフィッティングします。
初期値の設定、フィッティング範囲の指定、グラフ表示、結果の出力（Excel、パラメータファイル）に対応しています。
関連リンク: :doc:`lsq_func_tkFit_usage`
"""
import sys
import os
import numpy as np
from numpy import exp, log, log10, sin, cos, tan, arcsin, arccos, arctan, sinh, cosh, tanh, sqrt, abs
import pandas as pd
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import matplotlib.widgets as wg
from tklib.tkutils import replace_path, getarg, getintarg, getfloatarg, pint, pfloat
from tklib.tkapplication import tkApplication
from tklib.tkparams import tkParams
from tklib.tksci.tksci import asin, acos, atan, degcos, degsin, degtan, degacos, degasin, degatan, arcsin, arccos, arctan
from tklib.tksci.tksci import factorial, gamma, combination, eVTonm, nmToeV, Bn
from tklib.tksci.tksci import h, h_bar, hbar, e, kB, NA, c, pi, pi2, torad, todeg, basee
from tklib.tksci.tksci import me, mp, mn
from tklib.tksci.tksci import u0, e0, a0, R, F, g, G
from tklib.tksci.tksci import HartreeToeV, RyToeV, KToeV, eVToK, JToeV, eVToJ, Debye
from tklib.tksci.tkFit import tkFit
from tklib.tksci.tkFit_m import tkFit_m
#=======================================
# プログラム開始
#=======================================
[ドキュメント]
def initialize(app):
"""
アプリケーションの初期パラメータを設定します。
フィッティング方法、入出力ファイル名、フィッティング関数、初期パラメータ、
フィッティング範囲、数値微分時の微小変位、最大イテレーション回数、許容誤差、
グラフ表示に関する設定など、プログラムの動作に必要な様々なパラメータを初期化します。
:param app: tkApplication - アプリケーションインスタンス。
:returns: None
"""
cparams = tkParams()
app.cparams = cparams
#nelder-mead Downhill simplex
#powell Modified Powell
#cg conjugate gradient (Polak-Ribiere method)
#bfgs BFGS法
cparams.method = "bfgs"
#==========================================
# Source parameters to be fitted
#==========================================
cparams.infile = 'peak.xlsx'
cparams.nx = None
cparams.mode = 'fit'
cparams.func = "p[0]*exp(-((x[0]-p[1])/p[2])**2)"
cparams.p0s = "1.3,0.6,0.1"
cparams.fit_range = "-1e100:1e100, -1e100:1e100, -1e100:1e100"
cparams.olabel = ''
cparams.xlabel = ''
cparams.ylabel = ''
cparams.zlabel = ''
# 数値微分する際の変数の微小変位
cparams.h = 0.01
cparams.maxiter = 100
cparams.tol = 1.0e-5
#==========================================
# Graph parameters
#==========================================
cparams.ngdata = 51
cparams.xgmin = -4.0
cparams.xgmax = 4.0
cparams.tsleep = 0.01
cparams.figsize = (10, 5)
cparams.fontsize = 16
cparams.legend_fontsize = 10
cparams.outputinterval = 5
cparams.graphupdateinterval = 10
[ドキュメント]
def update_vars(app):
"""
コマンドライン引数に基づいてアプリケーションのパラメータを更新します。
`sys.argv` を解析し、`getarg`, `getintarg`, `getfloatarg` ヘルパー関数を使用して、
`app.cparams` に格納されているデフォルトのパラメータ値を上書きします。
これにより、スクリプトの実行時に外部から設定を変更できるようになります。
:param app: tkApplication - アプリケーションインスタンス。
:returns: None
"""
cparams = app.cparams
argv = sys.argv
# n = len(argv)
#for i in range(n):
# print("i=", i, argv[i])
cparams.mode = getarg( 1, cparams.mode)
cparams.infile = getarg( 2, cparams.infile)
cparams.method = getarg( 3, cparams.method)
cparams.func = getarg( 4, cparams.func)
cparams.p0s = getarg( 5, cparams.p0s)
cparams.fit_range = getarg( 6, cparams.fit_range)
cparams.maxiter = getintarg (7, cparams.maxiter)
cparams.tol = getfloatarg(8, cparams.tol)
cparams.h = getfloatarg(9, cparams.h)
cparams.olabel = getarg(10, cparams.olabel)
cparams.xlabel = getarg(11, cparams.xlabel)
cparams.ylabel = getarg(12, cparams.ylabel)
cparams.zlabel = getarg(13, cparams.zlabel)
# app.add_argument(opt = "-s", type = "str", var_name = 'script_list_name', opt_str = "-s=script_list_name", desc = 'Script list file',
# defval = "default", optional = True);
# app.add_argument(opt = "-i", type = "float", var_name = 'idx', opt_str = "-i=integer", desc = 'Index to speify alpha',
# defval = 5, optional = False);
# app.add_argument(opt = None, type = "str", var_name = 'infile', opt_str = "infile", desc = 'Input file',
# defval = "input.txt", optional = False);
# app.add_argument(opt = None, type = "str", var_name = 'nmax', opt_str = "nmax", desc = 'Max number of interation',
# defval = 5, optional = False);
# args_opt, args_idx, args_vars = app.read_args(vars = cparams.__dict__, check_allowed_args = True)
# 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)
# app.set_usage(usage_str)
#==========================================
# functions
#==========================================
[ドキュメント]
def cal_y(xd, pk, app):
"""
与えられた関数とパラメータに基づいてy値を計算します。
`app.cparams.func` に定義された文字列形式の関数を `eval` を用いて実行し、
独立変数 `xd` とフィッティングパラメータ `pk` を用いてy値を算出します。
この関数は、`tkFit_m` クラスの`func`属性として登録され、フィッティング計算に使用されます。
:param xd: list - 独立変数（xデータ）のリスト。
:param pk: list - フィッティングパラメータのリスト。
:param app: tkApplication - アプリケーションインスタンス。
:returns: float or array-like - 計算されたy値。
"""
y = eval(app.cparams.func, globals(), {"x": xd, "p": pk})
return y
[ドキュメント]
def fit(app):
"""
データに対する非線形最小二乗フィッティングを実行します。
入力データの読み込み、フィッティング関数の設定、初期パラメータの準備を行い、
初期プロットを表示します。その後、`cparams.method` で指定された最適化手法を用いて
フィッティングを実行し、最終結果の表示、フィッティング結果のExcelファイルへの保存、
パラメータファイルの保存を行います。最後に、最適化された関数のプロットを表示します。
:param app: tkApplication - アプリケーションインスタンス。
:returns: None
"""
cparams = app.cparams
app.log_path = app.replace_path(cparams.infile)
app.redirect(["stdout", app.log_path], "w")
cparams.p0s = [float(s) for s in cparams.p0s.split(',')]
r = []
for s in cparams.fit_range.split(','):
xmin, xmax = s.split(':')
r.append([float(xmin), float(xmax)])
cparams.fit_range = r
cparams.out_file = replace_path(cparams.infile, "{dirname}/{filebody}-fitting.xlsx")
cparams.convergence_file = replace_path(cparams.infile, "{dirname}/{filebody}-convergence.xlsx")
cparams.parameter_file = replace_path(cparams.infile, "{dirname}/{filebody}.prm")
print("")
print("Fitting program to given function")
print(f"mode : {cparams.mode}")
print(f"method : {cparams.method}")
print(f"input file : {cparams.infile}")
print(f" fit_range : {cparams.fit_range}")
print(f"output file : {cparams.out_file}")
print(f"convergence file: {cparams.convergence_file}")
print(f"parameter file : {cparams.parameter_file}")
print(f"func : {cparams.func}")
print(f"Target : {cparams.olabel}")
print(f"x : {cparams.xlabel}")
print(f"y : {cparams.ylabel}")
print(f"z : {cparams.zlabel}")
print(f"tol : {cparams.tol}")
print(f"maxiter : {cparams.maxiter}")
if '***' in cparams.method:
input("\nError: Choose method\n")
exit()
if '***' in cparams.func:
input("\nError: Choose or input function\n")
exit()
fit = tkFit_m(tol = cparams.tol, nmaxiter = cparams.maxiter,
print_interval = cparams.outputinterval, plot_interval = cparams.graphupdateinterval)
read_labels = []
if cparams.xlabel != '---':
read_labels.append(cparams.xlabel)
if cparams.ylabel != '---':
read_labels.append(cparams.ylabel)
if cparams.zlabel != '---':
read_labels.append(cparams.zlabel)
print("")
print(f"Read [{cparams.infile}]")
ranges = np.array(cparams.fit_range).T
print(" For ranges:", ranges)
fit.read_data(cparams.infile, xlabels = read_labels, ylabel = cparams.olabel,
xmins = ranges[0], xmaxs = ranges[1], usage = lambda: usage(app))
xd = fit.x_list
xd_labels = fit.xlabels
nx = len(xd)
ndata = len(xd[0])
nvars = len(cparams.p0s)
print("")
print("x ranges:")
for i in range(len(xd)):
xmin = min(xd[i])
xmax = max(xd[i])
print(f" x[{i}]: {xmin} - {xmax}")
print(f" fit in: {cparams.fit_range[i][0]} - {cparams.fit_range[i][1]}")
fit.varname = [f'p({i})' for i in range(nvars)]
fit.unit = ['' for i in range(nvars)]
fit.pk = [pfloat(cparams.p0s[i]) for i in range(nvars)]
fit.optid = [1 for i in range(nvars)]
# print("")
# print(f"initial values : {cparams.p0s}")
# print( " pk:", fit.pk)
fit.func = lambda x_list, pk: cal_y(x_list, pk, app)
#=============================
# 初期値関数
#=============================
print("")
print("Calculate initial function:")
yini = fit.cal_ylist(fit.pk)
fini = fit.minimize_func(fit.pk)
print("")
fit.print_data(heading = "Initial functions", yini = yini)
#=============================
# グラフの表示
#=============================
print("")
print("plot")
cparams.fig, axes = plt.subplots(1, 2, figsize = cparams.figsize)
# axes = [axes]
# axes = axes.flatten()
fit.initial_plot(data_axis = axes[0], error_axis = axes[1], yini = yini, fmin = fini,
fig = cparams.fig, fontsize = cparams.fontsize)
if cparams.mode == 'sim':
app.terminate("", pause = True)
#=============================
# Optimization
#=============================
print("")
print(f"Nonlinear least-squares fitting by [{cparams.method}]")
fit.print_variables()
print(f" tol={fit.tol}")
print(f" nmaxiter={fit.nmaxiter}")
pfin, ffin, success, res = fit.minimize(cparams.method)
if success:
print("")
print(f"\nConverged at iteration: {res.nit}")
else:
print(f"***Warning: Function did not converge")
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("Optimized at S2={:12.6g}:".format(ffin))
fit.print_variables()
#========================================
# Final result
#========================================
yfin = fit.cal_ylist(pfin)
print("")
# fit.print_data(heading = "Final functions", yini = yini, yfin = yfin)
fit.print_scores(heading = "\nScores between y(input) and y(fit)", y1 = fit.y, y2 = yfin)
print("")
print(f"Save results to [{cparams.out_file}]")
fit.to_excel(cparams.convergence_file,
[*xd_labels, cparams.olabel, 'initial', 'final'], [*xd, fit.y, yini, yfin])
print("")
print(f"Save convergence process to [{cparams.convergence_file}]")
fit.to_excel(cparams.convergence_file,
['iter', 'MSE'], [fit.iter_list, fit.fmin_list])
print("")
print(f"Save parameters to [{cparams.parameter_file}]")
cparams.fit_range = "{}".format(cparams.fit_range)
cresults = tkParams()
cresults.initial_values = f"{cparams.p0s}"
cresults.final_values = "{}".format(list(res.x))
cresults.MAE = res.fun
cresults.iteration = res.nit
del cparams.p0s
cresults.initial_values = cresults.initial_values.replace('[', '').replace(']', '')
cresults.final_values = cresults.final_values.replace('[', '').replace(']', '')
cparams.save_parameters(cparams.parameter_file, section = 'Parameters', sort_by_keys = False, update_commandline = False, IsPrint = False)
cresults.save_parameters(cparams.parameter_file, section = 'Fitting', sort_by_keys = False, update_commandline = False, IsPrint = False)
#=============================
# グラフの表示
#=============================
print("")
print("Plot optimized")
fit.finalize_plot(yfin, iter = res.nit, fmin = ffin)
app.terminate("", pause = True)
#==========================================
# Main routine
#==========================================
[ドキュメント]
def main():
"""
プログラムのエントリポイントです。
`tkApplication` インスタンスを作成し、以下の処理を順に実行します。
1. `initialize()`: アプリケーションの初期パラメータを設定します。
2. `update_vars()`: コマンドライン引数に基づいてパラメータを更新します。
3. `fit()`: データに対する非線形最小二乗フィッティングを実行します。
:returns: None
"""
app = tkApplication()
print(f"Initialize parameters")
initialize(app)
print(f"Update parameters by command-line arguments")
update_vars(app)
fit(app)
if __name__ == "__main__":
main()