# -*- coding: UTF-8 -*-
from math import *
import numpy as np
############################################
# 黄金分割法(与えられた方向での最小値)
# a,b : 初期区間 a < b
# eps : 許容誤差
# val : 関数値
# ind : 計算状況
# >= 0 : 正常終了(収束回数)
# = -1 : 収束せず
# max : 最大試行回数
# w : 位置
# dw : 傾きの成分
# snx : 関数値を計算する関数の名前
# return : 結果(w＋y＊dwのy)
# coded by Y.Suganuma
############################################
def gold(a, b, eps, val, ind, max, w, dw, snx) :
# 初期設定
tau = (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(x1, w, dw)
f2 = snx(x2, w, dw)
# 計算
while count < max and ind[0] < 0 :
count += 1
if f2 > f1 :
if abs(b-a) < 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(x1, w, dw)
else :
if abs(b-a) < 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(x2, w, dw)
# 収束した場合の処理
if ind[0] == 0 :
ind[0] = count
fa = snx(a, w, dw)
fb = snx(b, w, dw)
if fb < fa :
a = b
fa = fb
if fa < f1 :
x = a
val[0] = fa
return x
############################################
# DFP法 or BFGS法
# method : =0 : DFP法
# =1 : BFGS法
# opt_1 : =0 : １次元最適化を行わない
# =1 : １次元最適化を行う
# max : 最大繰り返し回数
# n : 次元
# eps : 収束判定条件
# step : きざみ幅
# y : 最小値
# x1 : x(初期値と答え)
# dx : 関数の微分値
# H : Hesse行列の逆行列
# snx : 関数値を計算する関数名
# dsnx : 関数の微分を計算する関数名(符号を変える)
# return : >=0 : 正常終了(収束回数)
# =-1 : 収束せず
# =-2 : １次元最適化の区間が求まらない
# =-3 : 黄金分割法が失敗
# coded by Y.Suganuma
############################################
def DFP_BFGS(method, opt_1, max, n, eps, step, y, x1, dx, H, snx, dsnx) :
count = 0
sw = 0
sw1 = np.empty(1, np.int)
y2 = np.empty(1, np.float)
x2 = np.empty(n, np.float)
g = np.empty(n, np.float)
g1 = np.empty(n, np.float)
g2 = np.empty(n, np.float)
s = np.empty(n, np.float)
w = np.empty(n, np.float)
I = np.identity(n, np.float)
H1 = np.empty((n, n), np.float)
H2 = np.empty((n, n), np.float)
y1 = snx(0.0, x1, dx)
dsnx(x1, g1)
while count < max and sw == 0 :
count += 1
# 方向の計算
for i1 in range(0, n) :
dx[i1] = 0.0
for i2 in range(0, n) :
dx[i1] -= H[i1][i2] * g1[i2]
# １次元最適化を行わない
if opt_1 == 0 :
# 新しい点
for i1 in range(0, n) :
x2[i1] = x1[i1] + step * dx[i1]
# 新しい関数値
y2[0] = snx(0.0, x2, dx)
# １次元最適化を行う
else :
# 区間を決める
sw1[0] = 0
f1 = y1
sp = step
f2 = snx(sp, x1, dx)
if f2 > f1 :
sp = -step
for i1 in range(0, max) :
f2 = snx(sp, x1, dx)
if f2 > f1 :
sw1[0] = 1
break
else :
sp *= 2.0
f1 = f2
# 区間が求まらない
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 range(0, n) :
x2[i1] = x1[i1] + k * dx[i1]
if sw == 0 :
# 収束（関数値の変化＜ｅｐｓ）
if abs(y2[0]-y1) < eps :
sw = count
y[0] = y2[0]
for i1 in range(0, n) :
x1[i1] = x2[i1]
# 関数値の変化が大きい
else :
# 傾きの計算
dsnx(x2, g2)
sm = 0.0
for i1 in range(0, n) :
sm += g2[i1] * g2[i1]
sm = sqrt(sm)
# 収束（傾き＜ｅｐｓ）
if sm < eps :
sw = count
y[0] = y2[0]
for i1 in range(0, n) :
x1[i1] = x2[i1]
# 収束していない
else :
for i1 in range(0, n) :
g[i1] = g2[i1] - g1[i1]
s[i1] = x2[i1] - x1[i1]
sm1 = 0.0
for i1 in range(0, n) :
sm1 += s[i1] * g[i1]
# DFP法
if method == 0 :
for i1 in range(0, n) :
w[i1] = 0.0
for i2 in range(0, n) :
w[i1] += g[i2] * H[i2][i1]
sm2 = 0.0
for i1 in range(0, n) :
sm2 += w[i1] * g[i1]
for i1 in range(0, n) :
w[i1] = 0.0
for i2 in range(0, n) :
w[i1] += H[i1][i2] * g[i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H1[i1][i2] = w[i1] * g[i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H2[i1][i2] = 0.0
for i3 in range(0, n) :
H2[i1][i2] += H1[i1][i3] * H[i3][i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H[i1][i2] = H[i1][i2] + s[i1] * s[i2] / sm1 - H2[i1][i2] / sm2
# BFGS法
else :
for i1 in range(0, n) :
for i2 in range(0, n) :
H1[i1][i2] = I[i1][i2] - s[i1] * g[i2] / sm1
for i1 in range(0, n) :
for i2 in range(0, n) :
H2[i1][i2] = 0.0
for i3 in range(0, n) :
H2[i1][i2] += H1[i1][i3] * H[i3][i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H[i1][i2] = 0.0
for i3 in range(0, n) :
H[i1][i2] += H2[i1][i3] * H1[i3][i2]
for i1 in range(0, n) :
for i2 in range(0, n) :
H[i1][i2] += s[i1] * s[i2] / sm1
y1 = y2[0]
for i1 in range(0, n) :
g1[i1] = g2[i1]
x1[i1] = x2[i1]
if sw == 0 :
sw = -1
return sw
----------------------------------
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
from function import DFP_BFGS
############################################
# 準Newton法によるの最小値の計算
# coded by Y.Suganuma
############################################
############################################
# 与えられた点xから，dx方向へk*dxだけ進んだ
# 点における関数値を計算する
############################################
# 関数１
def snx1(k, x, dx) :
w = np.empty(2, np.float)
w[0] = x[0] + k * dx[0]
w[1] = x[1] + k * dx[1]
x1 = w[0] - 1.0
y1 = w[1] - 2.0
y = x1 * x1 + y1 * y1
return y
# 関数２
def snx2(k, x, dx) :
w = np.empty(2, np.float)
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]
y = 100.0 * x1 * x1 + y1 * y1
return y
# 関数３
def snx3(k, x, dx) :
w = np.empty(2, np.float)
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])
y = x1 * x1 + y1 * y1 + z1 * z1
return y
############################################
# 微係数を計算する
############################################
# 関数１
def dsnx1(x, dx) :
dx[0] = -2.0 * (x[0] - 1.0)
dx[1] = -2.0 * (x[1] - 2.0)
# 関数２
def dsnx2(x, dx) :
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])
# 関数３
def dsnx3(x, dx) :
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]))
# データの入力
sw = 0
line = sys.stdin.readline()
s = line.split()
fun = int(s[1])
n = int(s[3])
max = int(s[5])
opt_1 = int(s[7])
line = sys.stdin.readline()
s = line.split()
method = int(s[1])
eps = float(s[3])
step = float(s[5])
x = np.empty(n, np.float)
dx = np.empty(n, np.float)
H = np.identity(n, np.float)
y = np.empty(1, np.float)
line = sys.stdin.readline()
s = line.split()
for i1 in range(0, n) :
x[i1] = float(s[i1+1])
# 実行
if fun == 1 :
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx1, dsnx1)
elif fun == 2 :
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx2, dsnx2)
else :
sw = DFP_BFGS(method, opt_1, max, n, eps, step, y, x, dx, H, snx3, dsnx3)
# 結果の出力
if sw < 0 :
print(" 収束しませんでした！", end="")
if sw == -1 :
print("（収束回数）")
elif sw == -2 :
print("（１次元最適化の区間）")
else :
print("（黄金分割法）")
else :
print(" 結果＝", end="")
for i1 in range(0, n) :
print(str(x[i1]) + " ", end="")
print(" 最小値＝" + str(y[0]) + " 回数＝" + str(sw))