"""
Simplex最適化アルゴリズムの実装。
このモジュールは、Nelder-Mead (Simplex) 法を用いた多変数関数の最小化機能を提供します。
初期シンプレックスの構築、反射、拡張、収縮、縮小の各ステップを通じて、
関数の局所最小値を見つけます。
:doc:`optimize-simplex_usage`
"""
import csv
import numpy as np
from numpy import sin, cos, tan, pi, exp, sqrt
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm
import time
import sys
#from goto import goto, label
[ドキュメント]
def func(x):
"""
最適化対象の目的関数を定義します。
この関数は2変数の二次関数であり、f(x0, x1) = (x0 - 1.0)^2 + (x1 - 3.0)^2 を計算します。
最小値はx = [1.0, 3.0]で0.0です。
:param x: 評価する変数のnumpy.ndarray。x[0]とx[1]を含む必要があります。
:type x: numpy.ndarray
:returns: 関数の評価値。
:rtype: float
"""
return pow(x[0] - 1.0, 2) + pow(x[1] - 3.0, 2)
[ドキュメント]
def main():
"""
Simplex最適化アルゴリズムを実行するメイン関数です。
最適化の初期条件（初期点、ステップサイズ、許容誤差、最大反復回数など）を設定し、
`simplex`関数を呼び出して関数の最小化を行います。
"""
x = np.array([0.0, 0.0])
dx = np.array([1.0, 1.0])
tolx = 1.0e-5
tolf = 1.0e-5
nvar = 2
fmin = 1.0e99
nmaxiter = 100
lprint = 3
iprintinterval = 10
simplex(func, x, dx, tolx, tolf, nmaxiter, lprint, iprintinterval)
"""
(1) FUNC (EXTERNAL) THE NAME OF THE FUNCTION TO BE MINIMIZED *
(INPUT) *
(2) X (D) THE INITIAL VALUES FOR THE VARIABLES (INPUT) *
THE VALUES WHICH GIVE THE MINIMUM (OUTPUT) *
(3) DX (D) THE EXTENTION OF THE FIRST SIMPLEX TO EACH DIRECTION*
(INPUT) *
(4) TOLX (D) THE TOLERANCE FOR THE VARIABLES (INPUT) *
(5) TOLF (D) THE TOLERANCE FOR THE CHANGE IN FUNCTION VALUES *
(INPUT) *
(6) N (I) THE NUMBER OF VARIABLES (INPUT) *
(7) FMIN (D) THE MINIMUM VALUE OF THE FUNCTION (OUTPUT) *
(8) ITMAX (I) THE MAXIMUM NUMBER OF ITERATIONS (INPUT) *
(9) LPRINT (I) THE VERBOSITY OF THE PRINT OUT (INPUT) *
COPYRIGHT: Y. OYANAGI, JUNE 30, 1989 V.1 *
"""
#double CAlgorism::Simplex(double (*func)(double x[]), double x[], double dx[], double tolx, double tolf,
# int n, double *fmin, int itmax, int lprint, int iprintinterval)
[ドキュメント]
def simplex(func, x, dx, tolx, tolf, itmax, lprint, iprintinterval):
"""
Nelder-Mead (Simplex) 法を用いて多変数関数を最小化します。
この関数は、反射、拡張、収縮、縮小といった幾何学的な操作を繰り返すことで、
与えられた目的関数の局所最小値を探索します。
初期シンプレックスを構築し、各反復でシンプレックスの頂点を更新して、
関数値が最も低い点へと移動させます。
:param func: 最小化される目的関数。1つのnumpy.ndarrayを引数にとり、floatを返す関数である必要があります。
:type func: callable
:param x: 最適化の開始点となる変数の初期値。最適化完了時には最小値を与える変数が格納されます。
:type x: numpy.ndarray
:param dx: 最初のシンプレックスを各次元に拡張するためのステップサイズ。
:type dx: numpy.ndarray
:param tolx: 変数の値に関する収束判定許容誤差。シンプレックスのサイズがこの値より小さくなると収束と判断されます。
:type tolx: float
:param tolf: 関数値の変化に関する収束判定許容誤差。シンプレックスの関数値の差がこの値より小さくなると収束と判断されます。
:type tolf: float
:param itmax: 最大反復回数。この回数を超えると探索を終了します。
:type itmax: int
:param lprint: 出力の詳細度を設定します。
0: 収束/反復上限に関する情報のみ。
1: 主要な反復情報（ITER, FL, FS, FH）を出力します。
2以上: 詳細な操作（REFLECTION, EXPANSION, CONTRACTION, REDUCTION）を出力します。
:type lprint: int
:param iprintinterval: lprintが1以上のとき、主要な反復情報を出力する間隔。
:type iprintinterval: int
:returns:
- fmin (float): 最小化された関数値。
- x (numpy.ndarray): 最小値を与える変数。
:rtype: tuple[float, numpy.ndarray]
.. note::
元のコードにはC/Fortran由来と思われる1ベースインデックスの記述や、
Pythonでは使用されない`goto`文のコメントアウトが存在します。
これらの記述は既存のロジックを変更しないというルールに基づき、そのまま残されています。
特に、配列のインデックス処理(`f[il - 1]`, `xx[ih][j+1]`)は、
一部で直感的ではないように見えるかもしれませんが、元のコードの意図を尊重しています。
"""
n = len(x)
alp = 1.0
bet = 0.5
gam = 2.0
eps = 1.0e-50
f = np.empty((n+1))
xr = np.empty((n))
xx = np.empty((n+1, n))
alp1 = 2.;
bet1 = .5;
gam1 = -1.;
# SET INITIAL SIMPLEX
for i in range(0, n + 1):
for j in range(0, n):
xx[i][j] = x[j]
for i in range(0, n):
xx[i][i] = x[i] + dx[i]
for i in range(0, n + 1):
for j in range(0, n):
x[j] = xx[i][j]
f[i] = func(x)
for iter in range(0, itmax):
ih = 1
fh = f[0]
il = 1
fl = f[0]
# FIND THE HIGHEST (IH, FH), THE SECOND HIGHEST (FS), */
# AND THE LOWEST (IL, FL) POINTS */
for i in range(0, n+1):
if f[i] > fh:
ih = i
fh = f[i]
if f[i] < fl:
il = i
fl = f[i]
# START OF THE MAIN LOOP
fs = fl
for i in range(0, n+1):
if i != ih and f[i] > fs:
fs = f[i]
if lprint >= 1:
if(iter-1) % iprintinterval == 0:
print("ITER={:3d} FL, FS, FH={},{},{}".format(iter, fl, fs, fh))
# CONVERGENCE CHECK
if tolx > eps:
xint = 0.0
i__1 = n;
for i1 in range(0, n):
i__3 = n + 1;
for i2 in range(i1 + 1, n+1):
for j in range(0, n):
xdif = abs(xx[i1][j] - xx[i2][j])
if xdif > xint:
xint = xdif
if tolf > 1e-50:
if fh - fl < tolf and xint < tolx:
# goto .L500 #convergence
if lprint >= 0:
print("(SUBR.SMPLX) CONVERGENCE AT ITER={}".format(iter))
for j in range(n):
x[j] = xx[il][j]
fmin = f[il - 1]
return fmin, x
elif xint < tolx:
# goto .L500
if lprint >= 0:
print("(SUBR.SMPLX) CONVERGENCE AT ITER={}".format(iter))
for j in range(n):
x[j] = xx[il][j]
fmin = f[il - 1]
return fmin, x
else:
if tolf > eps and fh - fl < tolf:
# goto .L500
if lprint >= 0:
print("(SUBR.SMPLX) CONVERGENCE AT ITER={}".format(iter))
for j in range(n):
x[j] = xx[il][j]
fmin = f[il - 1]
return fmin, x
# CONVERGENCE
# label. L500
# CUtils::printf("lprint=%ld\n", *lprint);
# if lprint >= 0:
# print("(SUBR.SMPLX) CONVERGENCE AT ITER={}".format(iter))
# for j in range(n):
# x[j] = xx[il][j]
# fmin = f[il - 1]
# return fmin, x
# REFLECTION
# X(*) IS THE CENTROID
for j in range(n):
for j in range(n):
s = 0.0
i__1 = n + 1;
for i in range(0, n+1):
s += xx[i][j]
x[j] = (s - xx[ih][j]) / n
# XR(*) IS THE NEW POINT BY REFELCTION */
for j in range(n-1): # range(n)
xr[j] = alp1 * x[j+1] - xx[ih][j+1]
fr = func(xr);
if lprint >= 2:
print("REFLECTION APPLIED. FR ={} XR(*)={}".format(fr, xr[0]))
# SEE WHERE FR LIES
if fr <= fs:
if fr < fl:
# EXPANSION X(*) WILL BE XE IN THE NEXT LOOP
for j in range(n):
x[j+1] = xr[j - 1] * 2.0 + gam1 * x[j+1]
fe = func(x)
if lprint >= 2:
# CUtils::printf("EXPANSION APPLIED. FE ={} XE(*)={}".format(fe, x[0]))
print("EXPANSION APPLIED. FE ={} XE(*)={}".format(fe, x[0]))
if fe < fr:
for j in range(n):
xx[ih][j] = x[j];
f[ih - 1] = fe;
if lprint >= 2:
print("XH(*) IS REPLACED BY XE(*)")
else:
for j in range(n):
xx[ih][j+1] = xr[j - 1]
f[ih - 1] = fr;
if lprint >= 2:
print("XH(*) IS REPLACED BY XR(*)")
else:
for j in range(n):
xx[ih][j+1] = xr[j];
f[ih - 1] = fr
if lprint >= 2:
print("XH(*) IS REPLACED BY XR(*)")
# CONTRACTION X(*) IS STILL THE CENTROID, BUT WILL SOON BE XC
else:
if fr < fh:
for j in range(n):
xx[ih][j+1] = xr[j];
f[ih - 1] = fr
if lprint >= 2:
print("XH(*) IS REPLACED BY XR(*)")
for j in range(n):
x[j] = xx[ih][j] * 0.5 + bet1 * x[j]
fc = func(x)
if lprint >= 2:
print("CONTRACTION APPLIED. FC ={} XC(*)={}".format(fc, x[0]))
if fc < fh:
for j in range(n):
xx[ih][j] = x[j];
f[ih - 1] = fc;
if lprint >= 2:
print("XH(*) IS REPLACED BY XC(*)")
else:
# REDUCTION
if iter == itmax:
# goto .L401
if lprint >= 0:
print("(SUBR.SMPLX) ITERATION TERMINATED DUE TO ITMAX")
for j in range(n):
x[j] = xx[il][j]
fmin = f[il - 1];
return fmin, x
for i in range(n+1):
for j in range(n):
x[j] = (xx[i][j] + xx[il][j]) * 0.5
xx[i][j] = x[j];
f[i] = func(x);
if lprint >= 2:
# CUtils::printf("REDUCTION APPLIED AROUND THE POINT: {}".format(xx[il][j]))
print("REDUCTION APPLIED AROUND THE POINT: {}".format(xx[il][j]))
# END OF THE MAIN LOOP
# label .L401
#CUtils::printf("lprint=%ld\n", *lprint);
if lprint >= 0:
print("(SUBR.SMPLX) ITERATION TERMINATED DUE TO ITMAX")
for j in range(n):
x[j] = xx[il][j]
fmin = f[il - 1];
return fmin, x
"""
static double Simplex(double (*func)(double x[]), CSimplexCtrl& ctrl)
{
double *x = new double[ctrl.nParam];
double *dx = new double[ctrl.nParam];
int count = 0;
int i;
for(i = 0 ; i < ctrl.nParam ; i++)
{
if(ctrl.id[i] == 0) continue;
x[count] = ctrl.x[i];
dx[count] = ctrl.dx[i];
count++;
}
double ret = Simplex(func, x, dx, ctrl.tolx, ctrl.tolf, count, &(ctrl.fmin),
ctrl.itmax, ctrl.lprint, ctrl.iPrintInterval);
count = 0;
for(i = 0 ; i < ctrl.nParam ; i++)
{
if(ctrl.id[i] == 0) continue;
ctrl.x[i] = x[count];
ctrl.dx[i] = dx[count];
count++;
}
delete []x;
delete []dx;
return ret;
};
"""
if __name__ == "__main__":
main()