#*********************************/
# 準Newton法によるの最小値の計算 */
# coded by Y.Suganuma */
#*********************************/
#********************************************/
# 与えられた点xから，dx方向へk*dxだけ進んだ */
# 点における関数値及び微係数を計算する */
#********************************************/
# 関数１
snx1 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(2)
w[0] = x[0] + k * dx[0]
w[1] = x[1] + k * dx[1]
x1 = w[0] - 1.0
y1 = w[1] - 2.0
x1 * x1 + y1 * y1
else
dx[0] = -2.0 * (x[0] - 1.0)
dx[1] = -2.0 * (x[1] - 2.0)
end
}
# 関数２
snx2 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(2)
w[0] = x[0] + k * dx[0]
w[1] = x[1] + k * dx[1]
x1 = w[1] - w[0] * w[0]
y1 = 1.0 - w[0]
100.0 * x1 * x1 + y1 * y1
else
dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
dx[1] = -200.0 * (x[1] - x[0] * x[0])
end
}
# 関数３
snx3 = Proc.new { |sw, k, x, dx|
if sw == 0
w = Array.new(2)
w[0] = x[0] + k * dx[0]
w[1] = x[1] + k * dx[1]
x1 = 1.5 - w[0] * (1.0 - w[1])
y1 = 2.25 - w[0] * (1.0 - w[1] * w[1])
z1 = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1])
x1 * x1 + y1 * y1 + z1 * z1
else
dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
end
}
#***************************************************************/
# 黄金分割法(与えられた方向での最小値) */
# a,b : 初期区間 a < b */
# eps : 許容誤差 */
# val : 関数値 */
# ind : 計算状況 */
# >= 0 : 正常終了(収束回数) */
# = -1 : 収束せず */
# max : 最大試行回数 */
# w : 位置 */
# dw : 傾きの成分 */
# snx : 関数値を計算する関数の名前 */
# return : 結果(w＋y＊dwのy) */
#***************************************************************/
def gold(a, b, eps, val, ind, max, w, dw, &snx)
# 初期設定
tau = (Math.sqrt(5.0) - 1.0) / 2.0
x = 0.0
count = 0
ind[0] = -1
x1 = b - tau * (b - a)
x2 = a + tau * (b - a)
f1 = snx.call(0, x1, w, dw)
f2 = snx.call(0, x2, w, dw)
# 計算
while count < max && ind[0] < 0
count += 1
if f2 > f1
if (b-a).abs() < eps
ind[0] = 0
x = x1
val[0] = f1
else
b = x2
x2 = x1
x1 = a + (1.0 - tau) * (b - a)
f2 = f1
f1 = snx.call(0, x1, w, dw)
end
else
if (b-a).abs() < eps
ind[0] = 0
x = x2
val[0] = f2
f1 = f2
else
a = x1
x1 = x2
x2 = b - (1.0 - tau) * (b - a)
f1 = f2
f2 = snx.call(0, x2, w, dw)
end
end
end
# 収束した場合の処理
if ind[0] == 0
ind[0] = count
fa = snx.call(0, a, w, dw)
fb = snx.call(0, b, w, dw)
if fb < fa
a = b
fa = fb
end
if fa < f1
x = a
val[0] = fa
end
end
return x
end
#*******************************************************/
# DFP法 or BFGS法 */
# method : =0 : DFP法 */
# =1 : BFGS法 */
# opt_1 : =0 : １次元最適化を行わない */
# =1 : １次元最適化を行う */
# max : 最大繰り返し回数 */
# n : 次元 */
# eps : 収束判定条件 */
# step : きざみ幅 */
# y : 最小値 */
# x1 : x(初期値と答え) */
# dx : 関数の微分値 */
# H : Hesse行列の逆行列 */
# snx : 関数値とその微分を計算する関数名 */
# （微分は，その符号を変える) */
# return : >=0 : 正常終了(収束回数) */
# =-1 : 収束せず */
# =-2 : １次元最適化の区間が求まらない */
# =-3 : 黄金分割法が失敗 */
#*******************************************************/
def DFP_BFGS(method, opt_1, max, n, eps, step, y, x1, dx, h, &snx)
count = 0
sw = 0
y2 = Array.new(1)
sw1 = Array.new(1)
x2 = Array.new(n)
g = Array.new(n)
g1 = Array.new(n)
g2 = Array.new(n)
s = Array.new(n)
w = Array.new(n)
y1 = snx.call(0, 0.0, x1, dx)
snx.call(1, 0.0, x1, g1)
h1 = Array.new(n)
h2 = Array.new(n)
i = Array.new(n)
for i1 in 0 ... n
h1[i1] = Array.new(n)
h2[i1] = Array.new(n)
i[i1] = Array.new(n)
for i2 in 0 ... n
h[i1][i2] = 0.0
i[i1][i2] = 0.0
end
h[i1][i1] = 1.0
i[i1][i1] = 1.0
end
while count < max && sw == 0
count += 1
# 方向の計算
for i1 in 0 ... n
dx[i1] = 0.0
for i2 in 0 ... n
dx[i1] -= h[i1][i2] * g1[i2]
end
end
# １次元最適化を行わない
if opt_1 == 0
# 新しい点
for i1 in 0 ... n
x2[i1] = x1[i1] + step * dx[i1]
end
# 新しい関数値
y2[0] = snx.call(0, 0.0, x2, dx)
# １次元最適化を行う
else
# 区間を決める
sw1[0] = 0
f1 = y1
sp = step
f2 = snx.call(0, sp, x1, dx)
if f2 > f1
sp = -step
end
for i1 in 0 ... max
f2 = snx.call(0, sp, x1, dx)
if f2 > f1
sw1[0] = 1
break
else
sp *= 2.0
f1 = f2
end
end
# 区間が求まらない
if sw1[0] == 0
sw = -2
# 区間が求まった
else
# 黄金分割法
k = gold(0.0, sp, eps, y2, sw1, max, x1, dx, &snx)
# 黄金分割法が失敗
if sw1[0] < 0
sw = -3
# 黄金分割法が成功
else
# 新しい点
for i1 in 0 ... n
x2[i1] = x1[i1] + k * dx[i1]
end
end
end
end
if sw == 0
# 収束（関数値の変化＜ｅｐｓ）
if (y2[0]-y1).abs() < eps
sw = count
y[0] = y2[0]
for i1 in 0 ... n
x1[i1] = x2[i1]
end
# 関数値の変化が大きい
else
# 傾きの計算
snx.call(1, 0.0, x2, g2)
sm = 0.0
for i1 in 0 ... n
sm += g2[i1] * g2[i1]
end
sm = Math.sqrt(sm)
# 収束（傾き＜ｅｐｓ）
if sm < eps
sw = count
y[0] = y2[0]
for i1 in 0 ... n
x1[i1] = x2[i1]
end
# 収束していない
else
for i1 in 0 ... n
g[i1] = g2[i1] - g1[i1]
s[i1] = x2[i1] - x1[i1]
end
sm1 = 0.0
for i1 in 0 ... n
sm1 += s[i1] * g[i1]
end
# DFP法
if method == 0
for i1 in 0 ... n
w[i1] = 0.0
for i2 in 0 ... n
w[i1] += g[i2] * h[i2][i1]
end
end
sm2 = 0.0
for i1 in 0 ... n
sm2 += w[i1] * g[i1]
end
for i1 in 0 ... n
w[i1] = 0.0
for i2 in 0 ... n
w[i1] += h[i1][i2] * g[i2]
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h1[i1][i2] = w[i1] * g[i2]
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h2[i1][i2] = 0.0
for i3 in 0 ... n
h2[i1][i2] += h1[i1][i3] * h[i3][i2]
end
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h[i1][i2] = h[i1][i2] + s[i1] * s[i2] / sm1 - h2[i1][i2] / sm2
end
end
# BFGS法
else
for i1 in 0 ... n
for i2 in 0 ... n
h1[i1][i2] = i[i1][i2] - s[i1] * g[i2] / sm1
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h2[i1][i2] = 0.0
for i3 in 0 ... n
h2[i1][i2] += h1[i1][i3] * h[i3][i2]
end
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h[i1][i2] = 0.0
for i3 in 0 ... n
h[i1][i2] += h2[i1][i3] * h1[i3][i2]
end
end
end
for i1 in 0 ... n
for i2 in 0 ... n
h[i1][i2] += s[i1] * s[i2] / sm1
end
end
end
y1 = y2[0]
for i1 in 0 ... n
g1[i1] = g2[i1]
x1[i1] = x2[i1]
end
end
end
end
end
if sw == 0
sw = -1
end
return sw
end
# データの入力
str = gets()
a = str.split(" ")
fun = Integer(a[1])
n = Integer(a[3])
max = Integer(a[5])
opt_1 = Integer(a[7])
str = gets();
a = str.split(" ");
method = Integer(a[1])
eps = Float(a[3])
step = Float(a[5])
x = Array.new(n)
dx = Array.new(n)
h = Array.new(n)
y = Array.new(1)
sw = 0
str = gets();
a = str.split(" ");
for i1 in 0 ... n
x[i1] = Float(a[i1+1])
dx[i1] = 0.0
h[i1] = Array.new(n)
end
# 実行
case fun
when 1
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx1)
when 2
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx2)
when 3
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, h, &snx3)
end
# 結果の出力
if sw < 0
printf(" 収束しませんでした！")
case sw
when -1
printf("（収束回数）\n")
when -2
printf("（１次元最適化の区間）\n")
when -3
printf("（黄金分割法）\n")
end
else
printf(" 結果＝")
for i1 in 0 ... n
printf("%f ", x[i1])
end
printf(" 最小値＝%f 回数＝%d\n", y[0], sw)
end