/*********************************/
/* Hopfieldネットワーク( TSP ) */
/* Coded by Y.Suganuma) */
/*********************************/
import java.io.*;
import java.util.*;
public class Hopfield_c {
static int n = 4;
static double A = 10.0, B = 10.0, C = 10.0, D = 2.0, d[][] = new double [n][n];
public static void main(String args[])
{
// 距離
double r = Math.sqrt(2.0);
for (int i1 = 0; i1 < n; i1++) {
for (int 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;
// 初期状態
double u[][] = new double [n][n];
Random rn = new Random();
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
u[i1][i2] = 0.1 * rn.nextDouble() - 0.05;
}
Console con = System.console();
con.printf("初期状態(出力):\n");
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
con.printf("%6.3f", out(u[i1][i2]));
con.printf("\n");
}
// 更新
double time = 0.0;
double dx[] = new double[n*n];
double uu[] = new double[n*n];
double g[][] = new double [4][n*n];
int k = 0;
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
uu[k] = u[i1][i2];
k++;
}
}
for (int i1 = 0; i1 < 100; i1++)
time = rkg(time, 1, uu, dx, g, n*n);
con.printf("最終状態(出力):\n");
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
con.printf("%6.3f", out(uu[i1*n+i2]));
con.printf("\n");
}
}
/**************************/
/* クロネッカのデルタ関数 */
/**************************/
static int delta(int i, int j)
{
if (i < 0)
i = n - 1;
else if (i >= n)
i = 0;
if (j < 0)
j = n - 1;
else if (j >= n)
j = 0;
int k = 0;
if (i == j)
k = 1;
return k;
}
/******************/
/* ユニットの出力 */
/******************/
static double out(double x)
{
return 0.5 * (1.0 + Math.tanh(x));
// return 1.0 / (1.0 + Math.exp(-x));
}
/****************/
/* 微係数の計算 */
/****************/
static void snx(double time, double u[], double du[])
{
for (int x = 0; x < n; x++) {
for (int i = 0; i < n; i++) {
int k = n * x + i;
du[k] = 0.0;
for (int j = 0; j < n; j++) {
if (j != i)
du[k] -= A * out(u[n*x+j]);
}
for (int y = 0; y < n; y++) {
if (y != x)
du[k] -= B * out(u[n*y+i]);
}
double N = 0.0;
for (int xx = 0; xx < n; xx++) {
for (int ii = 0; ii < n; ii++)
N += out(u[n*xx+ii]);
}
du[k] -= C * (N - n);
for (int y = 0; y < n; y++) {
int m1 = (i + 1) % n;
int 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 : 微分方程式の次数 */
/* return : time+h */
/*************************************************/
static double rkg(double time, double h, double x[], double dx[], double g[][], int n)
{
int i1;
double h2;
h2 = 0.5 * h;
snx(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];
snx(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];
snx(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];
snx(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;
}
}