シンプレックス法

#*************************************/
# シンプレックス法による最小値の計算 */
#      coded by Y.Suganuma           */
#*************************************/

#******************/
# 関数値の計算    */
#      y = f(x)   */
#      return : y */
#******************/
			# 関数１
snx1 = Proc.new { |x|
	x1 = x[0] - 1.0
	y1 = x[1] - 2.0
	x1 * x1 + y1 * y1
}
			# 関数２
snx2 = Proc.new { |x|
	x1 = x[1] - x[0] * x[0]
	y1 = 1.0 - x[0]
	100.0 * x1 * x1 + y1 * y1
}
			# 関数３
snx3 = Proc.new { |x|
	x1 = 1.5 - x[0] * (1.0 - x[1])
	y1 = 2.25 - x[0] * (1.0 - x[1] * x[1])
	z1 = 2.625 - x[0] * (1.0 - x[1] * x[1] * x[1])
	x1 * x1 + y1 * y1 + z1 * z1
}

#*****************************************/
# シンプレックス法                       */
#      max : 最大繰り返し回数            */
#      n : 次元                          */
#      k : 初期値設定係数                */
#      r1 : 縮小比率                     */
#      r2 : 拡大比率                     */
#      y : 最小値                        */
#      x : x(初期値と答え)               */
#      snx : 関数値を計算する関数名      */
#      return : >=0 : 正常終了(収束回数) */
#               =-1 : 収束せず           */
#*****************************************/
def simplex(max, n, k, eps, r1, r2, y, x, &snx)
					# 初期値の設定
	fh    = -1
	fg    = -1
	fl    = -1
	count = 0
	sw    = -1
	yy    = Array.new(n+1)
	xg    = Array.new(n)
	xr    = Array.new(n)
	xn    = Array.new(n)
	xx    = Array.new(n+1)
	for i1 in 0 ... n+1
		xx[i1] = Array.new(n)
	end
	for i1 in 0 ... n+1
		for i2 in 0 ... n
			xx[i1][i2] = x[i2]
		end
		if i1 > 0
			xx[i1][i1-1] += k
		end
		yy[i1] = snx.call(xx[i1])
	end
					# 最大値，最小値の計算
	for i1 in 0 ... n+1
		if fh < 0 || yy[i1] > yy[fh]
			fh = i1
		end
		if fl < 0 || yy[i1] < yy[fl]
			fl = i1
		end
	end
	for i1 in 0 ... n+1
		if i1 != fh && (fg < 0 || yy[i1] > yy[fg])
			fg = i1
		end
	end
					# 最小値の計算
	while count < max && sw < 0
		count += 1
							# 重心の計算
		for i1 in 0 ... n
			xg[i1] = 0.0
		end
		for i1 in 0 ... n+1
			if i1 != fh
				for i2 in 0 ... n
					xg[i2] += xx[i1][i2]
				end
			end
		end
		for i1 in 0 ... n
			xg[i1] /= n
		end
							# 最大点の置き換え
		for i1 in 0 ... n
			xr[i1] = 2.0 * xg[i1] - xx[fh][i1]
		end
		yr = snx.call(xr)
		if yr >= yy[fh]   # 縮小
			for i1 in 0 ... n
				xr[i1] = (1.0 - r1) * xx[fh][i1] + r1 * xr[i1]
			end
			yr = snx.call(xr)
		elsif yr < (yy[fl]+(r2-1.0)*yy[fh])/r2   # 拡大
			for i1 in 0 ... n
				xn[i1] = r2 * xr[i1] - (r2 - 1.0) * xx[fh][i1]
			end
			yn = snx.call(xn)
			if yn <= yr
				for i1 in 0 ... n
					xr[i1] = xn[i1]
				end
				yr = yn
			end
		end
		for i1 in 0 ... n
			xx[fh][i1] = xr[i1]
		end
		yy[fh] = yr
							# シンプレックス全体の縮小
		if yy[fh] >= yy[fg]
			for i1 in 0 ... n+1
				if i1 != fl
					for i2 in 0 ... n
						xx[i1][i2] = 0.5 * (xx[i1][i2] + xx[fl][i2])
					end
					yy[i1] = snx.call(xx[i1])
				end
			end
		end
							# 最大値，最小値の計算
		fh = -1
		fg = -1
		fl = -1
		for i1 in 0 ... n+1
			if fh < 0 || yy[i1] > yy[fh]
				fh = i1
			end
			if fl < 0 || yy[i1] < yy[fl]
				fl = i1
			end
		end
		for i1 in 0 ... n+1
			if i1 != fh && (fg < 0 || yy[i1] > yy[fg])
				fg = i1
			end
		end
							# 収束判定
		e = 0.0
		for i1 in 0 ... n+1
			if i1 != fl
				yr  = yy[i1] - yy[fl]
				e  += yr * yr
			end
		end
		if e < eps
			sw = 0
		end
	end

	if sw == 0
		sw   = count
		y[0] = yy[fl]
		for i1 in 0 ... n
			x[i1] = xx[fl][i1]
		end
	end

	return sw
end
				# データの入力
sw  = 0
str = gets()
a   = str.split(" ")
fun = Integer(a[1])
n   = Integer(a[3])
max = Integer(a[5])
k   = Float(a[7])
str = gets();
a   = str.split(" ");
eps = Float(a[1])
r1  = Float(a[3])
r2  = Float(a[5])
x   = 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])
end
				# 実行
case fun
	when 1
		sw = simplex(max, n, k, eps, r1, r2, y, x, &snx1)
	when 2
		sw = simplex(max, n, k, eps, r1, r2, y, x, &snx2)
	when 3
		sw = simplex(max, n, k, eps, r1, r2, y, x, &snx3)
end
				# 結果の出力
if sw < 0
	printf("   収束しませんでした！")
else
	printf("   結果＝")
	for i1 in 0 ... n
		printf("%f ", x[i1])
	end
	printf(" 最小値＝%f  回数＝%d\n", y[0], sw)
end