Hopfieldネットワーク（ TSP ）

<?php

/*********************************/
/* Hopfieldネットワーク（ TSP ） */
/*      Coded by Y.Suganuma)     */
/*********************************/

$n = 4;
$A = 10.0;
$B = 10.0;
$C = 10.0;
$D = 2.0;
$d = array(4);
for ($i1 = 0; $i1 < 4; $i1++)
	$d[$i1] = array(4);
					// 距離
$r = sqrt(2.0);
for ($i1 = 0; $i1 < $n; $i1++) {
	for ($i2 = 0; $i2 < $n; $i2++) {
		if ($i1 == $i2)
			$d[$i1][$i2] = 0.0;
		else
			$d[$i1][$i2] = 1.0;
	}
}
$d[0][2] = $r;
$d[1][3] = $r;
$d[2][0] = $r;
$d[3][1] = $r;
				// 初期状態
mt_srand();

$u = array(16);
for ($i1 = 0; $i1 < $n; $i1++) {
	for ($i2 = 0; $i2 < $n; $i2++)
		$u[$i1*$n+$i2] = 0.1 * (mt_rand() / mt_getrandmax()) - 0.05;
}

printf("初期状態（出力）：\n");
for ($i1 = 0; $i1 < $n; $i1++) {
	for ($i2 = 0; $i2 < $n; $i2++)
		printf("%6.3f", out($u[$i1*$n+$i2]));
	printf("\n");
}
				// 更新
$time = 0.0;
$dx   = array(16);
$g    = array(4);
for ($i1 = 0; $i1 < 4; $i1++)
	$g[$i1] = array(16);
for ($i1 = 0; $i1 < 100; $i1++)
	$time = rkg($time, 1, $u, $dx, $g, $n*$n, "snx");

printf("最終状態（出力）：\n");
for ($i1 = 0; $i1 < $n; $i1++) {
	for ($i2 = 0; $i2 < $n; $i2++)
		printf("%6.3f", out($u[$i1*$n+$i2]));
	printf("\n");
}

/******************/
/* ユニットの出力 */
/******************/
function out($x)
{
	return 0.5 * (1.0 + tanh($x));
//	return 1.0 / (1.0 + exp(-$x));
}

/****************/
/* 微係数の計算 */
/****************/
function snx($time, $u, &$du)
{
	global $n, $A, $B, $C, $D, $d;

	for ($x = 0; $x < $n; $x++) {
		for ($i = 0; $i < $n; $i++) {
			$k      = $n * $x + $i;
			$du[$k] = 0.0;
			for ($j = 0; $j < $n; $j++) {
				if ($j != $i)
					$du[$k] -= $A * out($u[$n*$x+$j]);
			}
			for ($y = 0; $y < $n; $y++) {
				if ($y != $x)
					$du[$k] -= $B * out($u[$n*$y+$i]);
			}
			$N = 0.0;
			for ($xx = 0; $xx < $n; $xx++) {
				for ($ii = 0; $ii < $n; $ii++)
					$N += out($u[$n*$xx+$ii]);
			}
			$du[$k] -= $C * ($N - $n);
			for ($y = 0; $y < $n; $y++) {
				$m1 = ($i + 1) % $n;
				$m2 = $i - 1;
				if ($m2 < 0)
					$m2 = $n - 1;
				$du[$k] -= $D * $d[$x][$y] * (out($u[$n*$y+$m1]) + out($u[$n*$y+$m2]));
			}
		}
	}
}

/*******************************************/
/* ルンゲ・クッタ法  dx/dt=f(t,x)          */
/*      time : 現在の時間                  */
/*      h : 時間刻み幅                     */
/*      x : 現在の状態                     */
/*      dx : 微係数（f(t,x):snxで計算する）*/
/*      g : 作業域（g[4][n]）              */
/*      n : 微分方程式の次数               */
/*      sub : 微係数を計算する関数の名前   */
/*      return : time+h                    */
/*******************************************/
function rkg($time, $h, &$x, $dx, $g, $n, $sub)
{
	$h2 = 0.5 * $h;

	$sub($time, $x, $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[0][$i1] = $h * $dx[$i1];

	$time += $h2;
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[1][$i1] = $x[$i1] + 0.5 * $g[0][$i1];
	$sub($time, $g[1], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[1][$i1] = $h * $dx[$i1];

	for ($i1 = 0; $i1 < $n; $i1++)
		$g[2][$i1] = $x[$i1] + 0.5 * $g[1][$i1];
	$sub($time, $g[2], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[2][$i1] = $h * $dx[$i1];

	$time += $h2;
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[3][$i1] = $x[$i1] + $g[2][$i1];
	$sub($time, $g[3], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[3][$i1] = $h * $dx[$i1];

	for ($i1 = 0; $i1 < $n; $i1++)
		$x[$i1] = $x[$i1] + ($g[0][$i1] + 2.0 * $g[1][$i1] + 2.0 * $g[2][$i1] + $g[3][$i1]) / 6.0;

	return $time;
}

?>