最適化(シンプレックス法)

<?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;
}

?>