<?php
/**************************************/
/* シンプレックス法による最小値の計算 */
/* coded by Y.Suganuma */
/**************************************/
$sw = 0;
// データの入力
fscanf(STDIN, "%*s %d %*s %d %*s %d %*s %lf", $fun, $n, $max, $k);
fscanf(STDIN, "%*s %lf %*s %lf %*s %lf", $eps, $r1, $r2);
$x = array($n);
$str = trim(fgets(STDIN));
strtok($str, " ");
for ($i1 = 0; $i1 < $n; $i1++)
$x[$i1] = floatval(strtok(" "));
// 実行
switch ($fun) {
case 1:
$sw = simplex($max, $n, $k, $eps, $r1, $r2, $y, $x, "snx1");
break;
case 2:
$sw = simplex($max, $n, $k, $eps, $r1, $r2, $y, $x, "snx2");
break;
case 3:
$sw = simplex($max, $n, $k, $eps, $r1, $r2, $y, $x, "snx3");
break;
}
// 結果の出力
if ($sw < 0)
printf(" 収束しませんでした!");
else {
printf(" 結果=");
for ($i1 = 0; $i1 < $n; $i1++)
printf("%f ", $x[$i1]);
printf(" 最小値=%f 回数=%d\n", $y, $sw);
}
/******************************************/
/* シンプレックス法 */
/* max : 最大繰り返し回数 */
/* n : 次元 */
/* k : 初期値設定係数 */
/* r1 : 縮小比率 */
/* r2 : 拡大比率 */
/* y : 最小値 */
/* x : x(初期値と答え) */
/* snx : 関数値を計算する関数名 */
/* return : >=0 : 正常終了(収束回数) */
/* =-1 : 収束せず */
/******************************************/
function simplex($max, $n, $k, $eps, $r1, $r2, &$y, &$x, $snx)
{
// 初期値の設定
$fh = -1;
$fg = -1;
$fl = -1;
$count = 0;
$sw = -1;
$yy = array($n+1);
$xg = array($n);
$xr = array($n);
$xn = array($n);
$xx = array($n+1);
for ($i1 = 0; $i1 < $n+1; $i1++)
$xx[$i1] = array($n);
for ($i1 = 0; $i1 < $n+1; $i1++) {
for ($i2 = 0; $i2 < $n; $i2++)
$xx[$i1][$i2] = $x[$i2];
if ($i1 > 0)
$xx[$i1][$i1-1] += $k;
$yy[$i1] = $snx($xx[$i1]);
}
// 最大値,最小値の計算
for ($i1 = 0; $i1 < $n+1; $i1++) {
if ($fh < 0 || $yy[$i1] > $yy[$fh])
$fh = $i1;
if ($fl < 0 || $yy[$i1] < $yy[$fl])
$fl = $i1;
}
for ($i1 = 0; $i1 < $n+1; $i1++) {
if ($i1 != $fh && ($fg < 0 || $yy[$i1] > $yy[$fg]))
$fg = $i1;
}
// 最小値の計算
while ($count < $max && $sw < 0) {
$count++;
// 重心の計算
for ($i1 = 0; $i1 < $n; $i1++)
$xg[$i1] = 0.0;
for ($i1 = 0; $i1 < $n+1; $i1++) {
if ($i1 != $fh) {
for ($i2 = 0; $i2 < $n; $i2++)
$xg[$i2] += $xx[$i1][$i2];
}
}
for ($i1 = 0; $i1 < $n; $i1++)
$xg[$i1] /= $n;
// 最大点の置き換え
for ($i1 = 0; $i1 < $n; $i1++)
$xr[$i1] = 2.0 * $xg[$i1] - $xx[$fh][$i1];
$yr = $snx($xr);
if ($yr >= $yy[$fh]) { // 縮小
for ($i1 = 0; $i1 < $n; $i1++)
$xr[$i1] = (1.0 - $r1) * $xx[$fh][$i1] + $r1 * $xr[$i1];
$yr = $snx($xr);
}
else if ($yr < ($yy[$fl]+($r2-1.0)*$yy[$fh])/$r2) { // 拡大
for ($i1 = 0; $i1 < $n; $i1++)
$xn[$i1] = $r2 * $xr[$i1] - ($r2 - 1.0) * $xx[$fh][$i1];
$yn = $snx($xn);
if ($yn <= $yr) {
for ($i1 = 0; $i1 < $n; $i1++)
$xr[$i1] = $xn[$i1];
$yr = $yn;
}
}
for ($i1 = 0; $i1 < $n; $i1++)
$xx[$fh][$i1] = $xr[$i1];
$yy[$fh] = $yr;
// シンプレックス全体の縮小
if ($yy[$fh] >= $yy[$fg]) {
for ($i1 = 0; $i1 < $n+1; $i1++) {
if ($i1 != $fl) {
for ($i2 = 0; $i2 < $n; $i2++)
$xx[$i1][$i2] = 0.5 * ($xx[$i1][$i2] + $xx[$fl][$i2]);
$yy[$i1] = $snx($xx[$i1]);
}
}
}
// 最大値,最小値の計算
$fh = -1;
$fg = -1;
$fl = -1;
for ($i1 = 0; $i1 < $n+1; $i1++) {
if ($fh < 0 || $yy[$i1] > $yy[$fh])
$fh = $i1;
if ($fl < 0 || $yy[$i1] < $yy[$fl])
$fl = $i1;
}
for ($i1 = 0; $i1 < $n+1; $i1++) {
if ($i1 != $fh && ($fg < 0 || $yy[$i1] > $yy[$fg]))
$fg = $i1;
}
// 収束判定
$e = 0.0;
for ($i1 = 0; $i1 < $n+1; $i1++) {
if ($i1 != $fl) {
$yr = $yy[$i1] - $yy[$fl];
$e += $yr * $yr;
}
}
if ($e < $eps)
$sw = 0;
}
if ($sw == 0) {
$sw = $count;
$y = $yy[$fl];
for ($i1 = 0; $i1 < $n; $i1++)
$x[$i1] = $xx[$fl][$i1];
}
return $sw;
}
/*******************/
/* 関数値の計算 */
/* y = f(x) */
/* return : y */
/*******************/
// 関数1
function snx1($x)
{
$x1 = $x[0] - 1.0;
$y1 = $x[1] - 2.0;
$y = $x1 * $x1 + $y1 * $y1;
return $y;
}
// 関数2
function snx2($x)
{
$x1 = $x[1] - $x[0] * $x[0];
$y1 = 1.0 - $x[0];
$y = 100.0 * $x1 * $x1 + $y1 * $y1;
return $y;
}
// 関数3
function snx3($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]);
$y = $x1 * $x1 + $y1 * $y1 + $z1 * $z1;
return $y;
}
?>