"""
最適化アルゴリズムを使用して関数の最小値を探索するスクリプト。
本スクリプトは、シンプレックス法、最急降下法、共役勾配法、ニュートン法などの
様々な最適化アルゴリズムを適用して、指定された目的関数の最小値を見つけます。
ユーザーは、楕円関数、一般関数、円関数の中から目的関数を選択でき、
初期点、アルゴリズム、ラインサーチモード、グラフ表示設定などを
コマンドライン引数またはスクリプト内の設定でカスタマイズ可能です。
最適化の過程はグラフとしてリアルタイムに可視化されます。
関連リンク: :doc:`optimize_usage`
"""
import sys
import numpy as np
from numpy import sin, cos, tan, pi, exp, sqrt
from tklib.tkargs import tkArgs
from tklib.tksci.tkoptimize import tkOptimize
from tklib.tkgraphic.tkplot import tkPlot
# target function
#functype = 'ellipsoid'
#functype = 'ellipsoid2'
#functype = 'circle'
functype = 'general'
#algorism: simplex, sd, cg, newton, broyden, dfp, bfgs
algorism = 'simplex'
#algorism = 'sd'
#algorism = 'cg'
#algorism = 'newton'
#algorism = 'broyden'
#algorism = 'dfp'
#algorism = 'bfgs'
# for line search '' or 'none'/'newton', 'one', 'simple', 'exact', 'golden', 'armijo'
lsmode = 'armijo'
# optimization parameters
x0 = [0.0, 0.0]
optid = [ 1, 1]
#for simplex
#dx = np.array([3.0, 3.0]) # range to make an initial simplex
dx = np.array([3.0, 3.0]) # range to make an initial simplex
# colormap for contour plot: line, hsv, cool, iwnter, gray, gist_gray, cividis etc
colormap = 'cool'
#colormap = 'hsv'
#colormap = 'Spectral'
# number of contours
nlevels = 51
if __name__ == "__main__":
if len(sys.argv) >= 2:
x0[0] = float(sys.argv[1])
if len(sys.argv) >= 3:
x0[1] = float(sys.argv[2])
if len(sys.argv) >= 4:
algorism = sys.argv[3]
if algorism == 'newton':
lsmode = 'newton'
if len(sys.argv) >= 5:
lsmode = sys.argv[4]
if len(sys.argv) >= 6:
functype = sys.argv[5]
if len(sys.argv) >= 7:
colormap = sys.argv[6]
if len(sys.argv) >= 8:
nlevels = int(sys.argv[7])
if len(sys.argv) >= 9:
initial_simplex = sys.argv[8]
if len(sys.argv) >= 10:
simplex_scale = float(sys.argv[9])
# for global optimization
nmaxiter = 200
tolx = 1.0e-5
tolf = 1.0e-5
iprintinterval = 1
print_level = 4
# for Newton-Raphson method
dump = 0.3
# line search parameters
ls_nmaxiter = 100
# for lsmode == 'exact'
ls_h = 1.0e-3 # x step to calculate first and second derivatives
# for lsmode == 'golden'
ls_alphaeps = 1.0e-2
ls_xrange = 0.5
if 'ellipsoid' in functype or 'circle' in functype:
ls_alpha = 0.03
else:
ls_alpha = 0.01
if lsmode == 'golden':
ls_alpha = 1.0
elif lsmode == 'armijo':
ls_alpha = 1.0
# graph plot parameters
fplot = 1
ngdata = 51
xgmin = -4.0
xgmax = 4.0
ygmin = -4.0
ygmax = 4.0
tsleep = 0.2
# file configuration
parameter_path = 'params.prm'
finalparameter_path = 'final_params.prm'
[ドキュメント]
def ediff1(i, x, optdata = None):
"""
目的関数がellipsoidの場合の第1次偏導関数を計算します。
軸が揃った楕円関数 f(x0, x1) = x0^2 + 9x1^2 の x_i に関する偏導関数を返します。
:param i: int - 偏導関数を計算する変数軸のインデックス (0または1)。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 指定された軸 i における目的関数の第1次偏導関数の値。
"""
if i == 0:
return 2.0 * x[0]
if i == 1:
return 18.0 * x[1]
[ドキュメント]
def gdiff1(i, x, optdata = None):
"""
目的関数がgeneralの場合の第1次偏導関数を計算します。
複雑な一般関数の x_i に関する偏導関数を返します。
:param i: int - 偏導関数を計算する変数軸のインデックス (0または1)。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 指定された軸 i における目的関数の第1次偏導関数の値。
"""
if i == 0:
return - 10.0 - 60.0 * x[0] + 4.5 * x[0]*x[0] + 12.0 * x[0]*x[0]*x[0] \
+ 3.0 * x[1] * x[1]
if i == 1:
return 30.0 - 60.0 * x[1] + 12.0 * x[1]*x[1]*x[1] \
+ 6.0 * x[0] * x[1]
[ドキュメント]
def gdiff2(i, j, x, optdata = None):
"""
目的関数がgeneralの場合の第2次偏導関数（ヘッセ行列の要素）を計算します。
複雑な一般関数の x_i, x_j に関する2次偏導関数を返します。
:param i: int - 第1変数軸のインデックス (0または1)。
:param j: int - 第2変数軸のインデックス (0または1)。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 指定された軸 i, j における目的関数の第2次偏導関数の値。
"""
if i == 0 and j == 0:
return -60.0 + 9.0 * x[0] + 36.0 * x[0]*x[0]
if i == 0 and j == 1:
return 6.0 * x[1]
if i == 1 and j == 0:
return 6.0 * x[1]
if i == 1 and j == 1:
return -60.0 +36.0 * x[1]*x[1] + 6.0 * x[0]
[ドキュメント]
def gfunc(x, optdata = None):
"""
一般的な非線形関数を評価します。
複数の項からなる複雑な2変数関数 f(x0, x1) を計算します。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 目的関数の値。
"""
return -3.0 - 10.0 * x[0] - 30.0 * x[0]*x[0] + 1.5 * x[0]*x[0]*x[0] + 3.0 * x[0]*x[0]*x[0]*x[0] \
+ 30.0 * x[1] - 30.0 * x[1]*x[1] + 3.0 * x[1]*x[1]*x[1]*x[1] \
+ 3.0 * x[0] * x[1] * x[1]
[ドキュメント]
def ediff2(i, j, x, optdata = None):
"""
目的関数がellipsoidの場合の第2次偏導関数（ヘッセ行列の要素）を計算します。
軸が揃った楕円関数 f(x0, x1) = x0^2 + 9x1^2 の x_i, x_j に関する2次偏導関数を返します。
:param i: int - 第1変数軸のインデックス (0または1)。
:param j: int - 第2変数軸のインデックス (0または1)。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 指定された軸 i, j における目的関数の第2次偏導関数の値。
"""
if i == 0 and j == 0:
return 2.0
if i == 0 and j == 1:
return 0.0
if i == 1 and j == 0:
return 0.0
if i == 1 and j == 1:
return 18.0
[ドキュメント]
def efunc(x, optdata = None):
"""
軸が揃った楕円関数を評価します。
f(x0, x1) = x0^2 + 9x1^2 を計算します。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 目的関数の値。
"""
return x[0]*x[0] + 9.0*x[1]*x[1]
ae2xx = 1.0
ae2xy = -2.0
ae2yy = 3.0
[ドキュメント]
def ediff2b(i, j, x, optdata = None):
"""
目的関数がellipsoid2の場合の第2次偏導関数（ヘッセ行列の要素）を計算します。
交差項を持つ楕円関数 f(x0, x1) = ae2xx*x0^2 + ae2xy*x0*x1 + ae2yy*x1^2 の
x_i, x_j に関する2次偏導関数を返します。
:param i: int - 第1変数軸のインデックス (0または1)。
:param j: int - 第2変数軸のインデックス (0または1)。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 指定された軸 i, j における目的関数の第2次偏導関数の値。
"""
if i == 0 and j == 0:
return 2.0 * acxx
if i == 0 and j == 1:
return ae2xy
if i == 1 and j == 0:
return ae2xy
if i == 1 and j == 1:
return 2.0 * acyy
[ドキュメント]
def ediff1b(i, x, optdata = None):
"""
目的関数がellipsoid2の場合の第1次偏導関数を計算します。
交差項を持つ楕円関数 f(x0, x1) = ae2xx*x0^2 + ae2xy*x0*x1 + ae2yy*x1^2 の
x_i に関する偏導関数を返します。
:param i: int - 偏導関数を計算する変数軸のインデックス (0または1)。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 指定された軸 i における目的関数の第1次偏導関数の値。
"""
if i == 0:
return 2.0 * ae2xx * x[0] + ae2xy * x[1]
if i == 1:
return ae2xy * x[0] + 2.0 * ae2yy * x[1]
[ドキュメント]
def efuncb(x, optdata = None):
"""
交差項を持つ楕円関数を評価します。
f(x0, x1) = ae2xx*x0^2 + ae2xy*x0*x1 + ae2yy*x1^2 を計算します。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 目的関数の値。
"""
return ae2xx * x[0]*x[0] + ae2xy * x[0]*x[1] + ae2yy * x[1]*x[1]
acxx = 3.0
[ドキュメント]
def ediff2c(i, j, x, optdata = None):
"""
目的関数がcircleの場合の第2次偏導関数（ヘッセ行列の要素）を計算します。
円関数 f(x0, x1) = acxx*(x0^2 + x1^2) の x_i, x_j に関する2次偏導関数を返します。
:param i: int - 第1変数軸のインデックス (0または1)。
:param j: int - 第2変数軸のインデックス (0または1)。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 指定された軸 i, j における目的関数の第2次偏導関数の値。
"""
if i == 0 and j == 0:
return 2.0 * acxx
if i == 0 and j == 1:
return 0.0
if i == 1 and j == 0:
return 0.0
if i == 1 and j == 1:
return 2.0 * acxx
[ドキュメント]
def ediff1c(i, x, optdata = None):
"""
目的関数がcircleの場合の第1次偏導関数を計算します。
円関数 f(x0, x1) = acxx*(x0^2 + x1^2) の x_i に関する偏導関数を返します。
:param i: int - 偏導関数を計算する変数軸のインデックス (0または1)。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 指定された軸 i における目的関数の第1次偏導関数の値。
"""
if i == 0:
return 2.0 * acxx * x[0]
if i == 1:
return 2.0 * acxx * x[1]
[ドキュメント]
def efuncc(x, optdata = None):
"""
円関数を評価します。
f(x0, x1) = acxx*(x0^2 + x1^2) を計算します。
:param x: list or numpy.ndarray - 変数 x の現在値 [x0, x1]。
:param optdata: dict or object, optional - 最適化に関連する追加データ（この関数では未使用）。
:returns: float - 目的関数の値。
"""
return acxx * (x[0]*x[0] + x[1]*x[1])
# global function variables
if functype == 'ellipsoid':
diff1 = ediff1
diff2 = ediff2
func = efunc
elif functype == 'ellipsoid2':
diff1 = ediff1b
diff2 = ediff2b
func = efuncb
elif functype == 'circle':
diff1 = ediff1c
diff2 = ediff2c
func = efuncc
else:
diff1 = gdiff1
diff2 = gdiff2
func = gfunc
[ドキュメント]
def cfunc(optdata):
"""
最適化プロセスの各イテレーションで呼び出されるコールバック関数です。
現在の最適化ステップにおけるパラメータをリカバリし、
グラフ表示が有効な場合は最適化の軌跡をプロットして可視化します。
この関数が0より大きい値を返すと、シンプレックスのイテレーションは終了します。
:param optdata: tklib.tksci.tkoptimize.tkOptimizeの内部データオブジェクト - 最適化に関する現在の状態を保持します。
:returns: int - 常に0を返し、最適化を継続させます。
"""
# print("callback at iter={}".format(optdata.iter))
x0 = optdata.RecoverParams(optdata.x0, optdata.x0_all, optdata.optid)
if optdata.fplot == 1:
graph = optdata.graph
ax = graph.axes[0]
ax2 = graph.axes[1]
if optdata.method == 'simplex':
vtx = optdata.vtx
x = []
y = []
for i in range(optdata.nvtx):
x.append(vtx[i].x[0])
y.append(vtx[i].x[1])
x.append(vtx[0].x[0])
y.append(vtx[0].x[1])
graph.plot(0, x, y, color = 'blue', marker = '', linestyle = '-')
else:
xt = ax.xlist[0]
yt = ax.ylist[0]
xt.append(x0[0])
yt.append(x0[1])
graph.set_data(0, 0, 0, xt, yt)
graph.pause(optdata.tsleep)
# graph.sleep(optdata.tsleep)
return 0
[ドキュメント]
def main():
"""
メイン処理を実行し、最適化プロセスを初期化、実行、可視化します。
コマンドライン引数を解析して最適化パラメータを設定し、tkOptimizeインスタンスを
作成して最適化関数とアルゴリズムを設定します。
初期パラメータの表示、目的関数の等高線および3D曲面のプロット、
最適化の実行、最終結果の表示と保存を行います。
"""
print("Find minimum point by steepest-descend / conjugate gradient / simplex methods")
print("")
print("x0= ({}, {})".format(x0[0], x0[1]))
print("algorism=", algorism)
print(f" For simplex:")
print(f" initial_simplex={initial_simplex}")
print(f" scale={simplex_scale}")
if len(sys.argv) >= 10:
simplex_scale = float(sys.argv[9])
print("lsmode=", lsmode)
print("functype=", functype)
print("colormap=", colormap)
print("")
for i in range(len(dx)):
dx[i] = simplex_scale
opt = tkOptimize()
opt.add_parameter("x", x0[0], dx[0], optid[0])
opt.add_parameter("y", x0[1], dx[1], optid[1])
opt.set_method(algorism, lsmode, initial_simplex = initial_simplex)
opt.set_functions(fitfunc = func, func = func, diff1func = diff1, diff2func = diff2)
opt.initialize(callback = cfunc)
# opt.initialize(
# nmaxiter = nmaxiter, tolx = tolx, tolf = tolf,
# callback = cfunc, print_level = print_level, iprintinterval = iprintinterval,
# dump = dump,
# ls_xrange = ls_xrange, ls_h = ls_h, ls_alpha = ls_alpha, ls_dump = dump,
# ls_nmaxiter = ls_nmaxiter, ls_alphaeps = ls_alphaeps
# )
var = opt.make_optdata(tsleep = tsleep)
# inf = opt.read_parameters(parameter_path, update_params = 0)
# print("parameters in [{}]".format(parameter_path))
# for key, val in inf.items():
# print("{}: {}".format(key, val))
# print("save to [{}]".format(parameter_path))
# opt.save_parameters(parameter_path)
# exit()
print("")
print("Initial parameters:")
opt.print_parameters()
# plot surface graph
if fplot == 1:
xgstep = (xgmax - xgmin) / (ngdata - 1)
ygstep = (ygmax - ygmin) / (ngdata - 1)
xg = np.empty([ngdata, ngdata])
yg = np.empty([ngdata, ngdata])
zg = np.empty([ngdata, ngdata])
for ix in range(ngdata):
for iy in range(ngdata):
xg[ix, iy] = xgmin + ix * xgstep
yg[ix, iy] = ygmin + iy * ygstep
zg[ix, iy] = func([xg[ix][iy], yg[ix][iy]])
graph = tkPlot(figsize = (10, 5))
graph.add_subplot(1, 2, 0)
graph.add_subplot(1, 2, 1, projection = '3d')
f = func(x0)
print("x0 = ({}, {}): f = {}".format(x0[0], x0[1], f))
xt = []
yt = []
xt.append(x0[0])
yt.append(x0[1])
graph.plot(0, xt, yt, color = 'blue', linestyle = '-', linewidth = 0.5,
fillstyle = 'full', marker = 'o', markersize = 5)
if colormap == 'line':
graph.plot_contour(0, xg, yg, zg, levels = nlevels,
cmap = 'line', colors = ['black'], linewidths = 0.5, aspect = 'equal')
else:
graph.plot_contour(0, xg, yg, zg, levels = nlevels, cmap = colormap, aspect = 'equal')
# ax.set_aspect('equal')
graph.set_axtitle(0, 'contour')
graph.plot_wireframe(1, xg, yg, zg, rstride = 2, cstride = 2)
graph.pause()
# graph.show()
var.graph = graph
var.fplot = fplot
# calculate initial parameters
xmin, fmin, optdata = opt.optimize()
print("")
print("Optimized parameters:")
opt.print_parameters(x0 = xmin, f = fmin)
if opt.save_parameters(finalparameter_path):
print("")
print("Final resuls are stored to [{}]".format(finalparameter_path))
else:
print("")
print("Error: Can not save the results to [{}]".format(finalparameter_path))
if fplot == 1:
print("Press enter to terminate:", end = '')
ret = input()
if __name__ == "__main__":
main()