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

/* 問題の定義 */
#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

/* ボルツマンマシンの動作の定義 */
#define T_0		10
#define T_MIN		2
#define T_next(t)	t*0.5

/* ニューロンの定義 */
typedef struct {
	double w[N];		/* 結合係数 */
	double theta;		/* 閾値 */
	int 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, loop;
	double T, u, p;

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

	/* ネットワークの初期状態の設定 */
	for(i=0; i<N; i++)
		neuron[i].x=rnd()>0.0;
	for(T=T_0; T>T_MIN; T=T_next(T))
		for(i=0; i<N; i++) {	/* 順番にユニットを選択 */
			printf("%.2lfT(%d) ", T, i);
			/* 状態の出力 */
			for(j=0; j<N; j++)
				printf("%d ", neuron[j].x);
			printf(" => %lf\n", Energy());
			/* 出力の計算 */
			u=0.0;
			for(j=0; j<N; j++)
				u+=neuron[i].w[j]*neuron[j].x;
			u-=neuron[i].theta;
			p=1/(1+exp(-u/T));
			if(p>rnd())
				neuron[i].x=1;
      			else
				neuron[i].x=0;
	    	}
}