# -*- 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
############################################
# 線形連立方程式を解く（逆行列を求める）
# w : 方程式の左辺及び右辺
# n : 方程式の数
# m : 方程式の右辺の列の数
# eps : 逆行列の存在を判定する規準
# return : =0 : 正常
# =1 : 逆行列が存在しない
# coded by Y.Suganuma
############################################
def gauss(w, n, m, eps) :
nm = n + m
ind = 0
for i1 in range(0, n) :
y1 = 0.0
m1 = i1 + 1
m2 = 0
# ピボット要素の選択
for i2 in range(i1, n) :
y2 = abs(w[i2][i1])
if y1 < y2 :
y1 = y2
m2 = i2
# 逆行列が存在しない
if y1 < eps :
ind = 1
break
# 逆行列が存在する
else : # 行の入れ替え
for i2 in range(i1, nm) :
y1 = w[i1][i2]
w[i1][i2] = w[m2][i2]
w[m2][i2] = y1
# 掃き出し操作
y1 = 1.0 / w[i1][i1]
for i2 in range(m1, nm) :
w[i1][i2] *= y1
for i2 in range(0, n) :
if i2 != i1 :
for i3 in range(m1, nm) :
w[i2][i3] -= (w[i2][i1] * w[i1][i3])
return ind
############################################
# Newton法
# opt_1 : =0 : １次元最適化を行わない
# =1 : １次元最適化を行う
# max : 最大繰り返し回数
# n : 次元
# eps : 収束判定条件
# step : きざみ幅
# y : 最小値
# x : x(初期値と答え)
# dx : 関数の微分値
# H : Hesse行列の逆行列
# snx : 関数値を計算する関数名
# dsnx : 関数の微分を計算する関数名(符号を変える)
# hesse : Hesse行列の逆行列を計算する関数名
# return : >=0 : 正常終了(収束回数)
# =-1 : 収束せず
# =-2 : １次元最適化の区間が求まらない
# =-3 : 黄金分割法が失敗
# coded by Y.Suganuma
############################################
def Newton(opt_1, max, n, eps, step, y, x, dx, H, snx, dsnx, hesse) :
wk = np.empty(n, np.float)
count = 0
sw = 0
y2 = np.empty(1, np.float)
sw1 = np.empty(1, np.int)
y1 = snx(0.0, x, dx)
while count < max and sw == 0 :
# 傾きの計算
dsnx(x, wk)
# Hesse行列の逆行列の計算
sw1[0] = hesse(x, H, eps)
# 収束していない
if sw1[0] == 0 :
# 方向の計算
count += 1
for i1 in range(0, n) :
dx[i1] = 0.0
for i2 in range(0, n) :
dx[i1] += H[i1][n+i2] * wk[i2]
# １次元最適化を行わない
if opt_1 == 0 :
# 新しい点
for i1 in range(0, n) :
x[i1] += dx[i1]
# 新しい関数値
y2[0] = snx(0.0, x, dx)
# 関数値の変化が大きい
if abs(y2[0]-y1) > eps :
y1 = y2[0]
# 収束（関数値の変化＜ｅｐｓ）
else :
sw = count
y[0] = y2[0]
# １次元最適化を行う
else :
# 区間を決める
sw1[0] = 0
f1 = y1
sp = step
f2 = snx(sp, x, dx)
if f2 > f1 :
sp = -step
for i1 in range(0, max) :
f2 = snx(sp, x, 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, x, dx, snx)
# 黄金分割法が失敗
if sw1[0] < 0 :
sw = -3
# 黄金分割法が成功
else :
# 新しい点
for i1 in range(0, n) :
x[i1] += k * dx[i1]
# 関数値の変化が大きい
if abs(y1-y2[0]) > eps :
y1 = y2[0]
# 収束（関数値の変化＜ｅｐｓ）
else :
sw = count
y[0] = y2[0]
# 収束（傾き＜ｅｐｓ）
else :
sw = count
y[0] = y1
if sw == 0 :
sw = -1
return sw
----------------------------------
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
from function import Newton, gauss
############################################
# 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]))
############################################
# Hesse行列の逆行列を計算する */
############################################
# 関数１
def hesse1(x, H, eps) :
H[0][0] = 2.0
H[0][1] = 0.0
H[1][0] = 0.0
H[1][1] = 2.0
H[0][2] = 1.0
H[0][3] = 0.0
H[1][2] = 0.0
H[1][3] = 1.0
sw = gauss(H, 2, 2, eps)
return sw
# 関数２
def hesse2(x, H, eps) :
H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
H[0][1] = -400.0 * x[0]
H[1][0] = H[0][1]
H[1][1] = 200.0
H[0][2] = 1.0
H[0][3] = 0.0
H[1][2] = 0.0
H[1][3] = 1.0
sw = gauss(H, 2, 2, eps)
return sw
# 関数３
def hesse3(x, H, eps) :
x1 = 1.0 - x[1]
x2 = 1.0 - x[1] * x[1]
x3 = 1.0 - x[1] * x[1] * x[1]
H[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3
H[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 + 4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 + 6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3
H[1][0] = H[0][1]
H[1][1] = 2.0 * x[0] * x[0] + 4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] + 12.0 * x[0] * x[1] * (2.625 - x[0] * x3) + 18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1]
H[0][2] = 1.0
H[0][3] = 0.0
H[1][2] = 0.0
H[1][3] = 1.0
sw = gauss(H, 2, 2, eps)
return sw
# データの入力
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()
eps = float(s[1])
step = float(s[3])
x = np.empty(n, np.float)
dx = np.empty(n, np.float)
H = np.empty((n, 2*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 = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx1, dsnx1, hesse1)
elif fun == 2 :
sw = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx2, dsnx2, hesse2)
else :
sw = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx3, dsnx3, hesse3)
# 結果の出力
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))