#*******************************/
# 共役勾配法による最小値の計算 */
# 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
#*******************************************************/
# 共役勾配法 */
# opt_1 : =0 : １次元最適化を行わない */
# =1 : １次元最適化を行う */
# max : 最大繰り返し回数 */
# n : 次元 */
# eps : 収束判定条件 */
# step : きざみ幅 */
# y : 最小値 */
# x : x(初期値と答え) */
# dx : 関数の微分値 */
# snx : 関数値とその微分を計算する関数名 */
# （微分は，その符号を変える) */
# return : >=0 : 正常終了(収束回数) */
# =-1 : 収束せず */
# =-2 : １次元最適化の区間が求まらない */
# =-3 : 黄金分割法が失敗 */
#*******************************************************/
def conjugate(opt_1, max, n, eps, step, y, x, dx, &snx)
b = 0.0
s1 = 1.0
count = 0
sw = 0
y2 = Array.new(1)
sw1 = Array.new(1)
g = Array.new(n)
y1 = snx.call(0, 0.0, x, dx)
while count < max && sw == 0
# 傾きの計算
snx.call(1, 0.0, x, g)
# 傾きのチェック
s2 = 0.0
for i1 in 0 ... n
s2 += g[i1] * g[i1]
end
# 収束していない
if Math.sqrt(s2) > eps
# 方向の計算
count += 1
if count%n == 1
b = 0.0
else
b = s2 / s1
end
for i1 in 0 ... n
dx[i1] = g[i1] + b * dx[i1]
end
# １次元最適化を行わない
if opt_1 == 0
# 新しい点
for i1 in 0 ... n
x[i1] += step * dx[i1]
end
# 新しい関数値
y2[0] = snx.call(0, 0.0, x, dx)
# 関数値の変化が大きい
if (y2[0]-y1).abs() > eps
y1 = y2[0]
s1 = s2
# 収束（関数値の変化＜ｅｐｓ）
else
sw = count
y[0] = y2[0]
end
# １次元最適化を行う
else
# 区間を決める
sw1[0] = 0
f1 = y1
sp = step
f2 = snx.call(0, sp, x, dx)
if f2 > f1
sp = -step
end
for i1 in 0 ... max
f2 = snx.call(0, sp, x, 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, x, dx, &snx)
# 黄金分割法が失敗
if sw1[0] < 0
sw = -3
# 黄金分割法が成功
else
# 新しい点
for i1 in 0 ... n
x[i1] += k * dx[i1]
end
# 関数値の変化が大きい
if (y1-y2[0]).abs() > eps
y1 = y2[0]
s1 = s2
# 収束（関数値の変化＜ｅｐｓ）
else
sw = count
y[0] = y2[0]
end
end
end
end
# 収束（傾き＜ｅｐｓ）
else
sw = count
y[0] = y1
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(" ")
eps = Float(a[1])
step = Float(a[3])
x = Array.new(n)
dx = 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
end
# 実行
case fun
when 1
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, &snx1)
when 2
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, &snx2)
when 3
sw = conjugate(opt_1, max, n, eps, step, y, x, dx, &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