#include <math.h>
#include <stdlib.h>

/* 連続値ホップフィールドモデルの諸定数 */
#define TAU		1
#define DELTA_T		0.02
#define MU0		0.1
#define f(u)		(1/(1+exp(-2*u/MU0)))

/* 問題の定義 */
#define MAX_WEIGHT	14
int a[]={	/* 品物の重量 */
	3, 5, 4, 7, 2, 4
};
int c[]={	/* 品物の価値 */
	10, 13, 10, 16, 2, 3
};
/* 解=1 0 1 1 0 0 */
		/* 品物の数 */
#define N	(sizeof(a)/sizeof(a[0]))

/* エネルギー関数のための係数 */
#define A		0.3

/* ニューロンの定義 */
typedef struct {
	double w[N];		/* 結合係数 */
	double theta;		/* 閾値 */
	double u;		/* 内部状態 */
	double x;		/* 出力 */
} NeuronRec;

NeuronRec neuron[N];

#define rnd()	((double)rand()/RAND_MAX)
		/* [0, 1] の乱数 */

/* エネルギー関数 */
double Energy()
{
	register int i, j;
	double E;

	E=0.0;
	for(i=0; i<N; i++) {
		for(j=0; j<i; j++)
			E+=-neuron[i].w[j]*neuron[i].x*neuron[j].x;
		E+=neuron[i].theta*neuron[i].x;
	}
	/* 簡単のため，積分項を無視 */
	return(E);
}

/* Hopfield ネットワークの定義 */
Initialize()
{
	register int i, j;

	for(i=0; i<N; i++) {
		neuron[i].w[i]=0.0;
		for(j=0; j<i; j++)
			neuron[i].w[j]=neuron[j].w[i]=-2*A*a[i]*a[j];
		neuron[i].theta=-A*MAX_WEIGHT*a[i]-c[i];
	}
}

main()
{
	register int i, j, t;
	int changed;
	double x, du;

	/* 結合荷重と閾値の設定 */
	Initialize();

	/* ネットワークの初期状態の設定 */
	for(i=0; i<N; i++)
		neuron[i].u=rnd();
	for(t=0; ; t++) {
		changed=0;
		/* 出力の計算 */
		for(i=0; i<N; i++) {
			x=f(neuron[i].u);
			if(x!=neuron[i].x) {
				neuron[i].x=x;
				changed=1;
			}
		}

		/* 出力の表示 */
		printf("%3d: ", t);
		for(i=0; i<N; i++)
			printf("%lf ", neuron[i].x);
		printf(" => %lf\n", Energy());

		if(!changed) break;

		/* 状態の計算 */
		for(i=0; i<N; i++) {
			du=-neuron[i].u/TAU;
			for(j=0; j<N; j++)
				du+=neuron[i].w[j]*neuron[j].x;
			du-=neuron[i].theta;
			neuron[i].u+=du*DELTA_T;
		}
	}
}