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

/* 層の定義 */
#define R	3		/* 層の数 */
#define N1	2		/* 入力層のユニット数 */
#define NR	1		/* 出力層のユニット数 */
#define MAX_N	5		/* 層中の最大ユニット数 */

int n[R]={N1, 2, NR};		/* 各層のユニット数 */
typedef struct {
	double w[MAX_N][MAX_N];
	double theta[MAX_N];
	double output[MAX_N];
	double delta[MAX_N];
/* モーメント法 */
	double w_t1[MAX_N][MAX_N];
	double theta_t1[MAX_N];
} LayerRec;
LayerRec L[R];

/* 教師データの定義 */
#define P	4		/* 教師データの数 */
typedef struct {
	double x[N1];
	double t[NR];
} TeachData;

TeachData T[P]={		/* XOR */
	0, 0, 0,
	0, 1, 1,
	1, 0, 1,
	1, 1, 0,
};

/* 学習時の定数の定義 */
#define EPSILON		1.0e-3	/* 学習を停止する許容誤差 */
#define INIT_W		0.5	/* 結合荷重の初期乱数の範囲 */
#define INIT_THETA	1.0	/* 閾値の初期乱数の範囲 */
#define ETA_W		0.5	/* 結合荷重の学習係数 */
#define ETA_THETA	1.0	/* 閾値の学習係数 */
/* モーメント法 */
#define ALPHA_MIN	0.2	/* モーメント法の例 0.2 */
#define ALPHA_MAX	0.9
#define ALPHA_DELTA	0.001	/* 修正モーメント法の例 0.001 */
double Alpha;

#define f(u)		(1/(1+exp(-(u))))
#define f_dash(o)	(o*(1.0-o))
#define rnd()		(2.0*(double)rand()/RAND_MAX-1.0)
			/* [-1, 1] の乱数 */
#define square(x)	((x)*(x))

#define BeginTable()	printf("     t     E w1-11 w1-12 θ1-1 w1-21 w1-22 θ1-2 w2-11 w2-12 θ2-1\n")
#define BeginResultTable()	printf("    p    x1    x2    tp     o\n")
#define EndTable()	printf("\n")
#define BeginRow()
#define EndRow()	printf("\n")
#define PrintValue(v)	printf("%5g ", v)

/* ニューラルネットワークの初期化 */
Initialize()
{
	int i, j, k;

	for(k=1; k<R; k++)
	    for(i=0; i<n[k]; i++) {
		for(j=0; j<n[k-1]; j++) {
			L[k].w[i][j]=INIT_W*rnd();
			L[k].w_t1[i][j]=0;
		}
		L[k].theta[i]=INIT_THETA*rnd();
		L[k].theta_t1[i]=0;
	    }
}

/* 入力→出力 */
Forward(int p)
{
	int i, j, k;
	double u;

	k=0;
	for(i=0; i<n[k]; i++)
	    L[k].output[i]=T[p].x[i];
	for(k=1; k<R; k++)
	    for(i=0; i<n[k]; i++) {
		u=0.0;
		for(j=0; j<n[k-1]; j++)
			u+=L[k].w[i][j]*L[k-1].output[j];
		u-=L[k].theta[i];
		L[k].output[i]=f(u);
	    }
}

/* 出力→入力 */
Backward(int p)
{
	int i, j, k, l;
	double delta;

	k=R-1;
	for(i=0; i<n[k]; i++)
	    L[k].delta[i]= -(T[p].t[i]-L[k].output[i])*f_dash(L[k].output[i]);
	/* 中間層の誤差 */
	for(k=R-2; k>=1; k--)
	    for(i=0; i<n[k]; i++) {
		delta=0.0;
		for(l=0; l<n[k+1]; l++)
			delta+=L[k+1].w[l][i]*L[k+1].delta[l];
		L[k].delta[i]=delta*f_dash(L[k].output[i]);
	    }
}

/* 学習 */
Learn()
{
	int i, j, k;
	double delta;

	for(k=R-1; k>=1; k--)
	    for(i=0; i<n[k]; i++) {
		for(j=0; j<n[k-1]; j++) {
			delta=-ETA_W*L[k].delta[i]*L[k-1].output[j]+Alpha*L[k].w_t1[i][j];
			L[k].w[i][j]+=delta;
			L[k].w_t1[i][j]=delta;
		}
		delta=ETA_THETA*L[k].delta[i]+Alpha*L[k].theta_t1[i];
		L[k].theta[i]+=delta;
		L[k].theta_t1[i]=delta;
	    }
}

/* 誤差関数 */
double MSE()
{
	int i, k, p;
	double E;

	k=R-1;	/* 出力層 */
	E=0.0;
	for(p=0; p<P; p++) {
		Forward(p);
		for(i=0; i<n[k]; i++)
			E+=square(T[p].t[i]-L[k].output[i]);
	}
	return(E/P);
}

/* ニューラルネットワークの表示 */
PrintNetwork()
{
	int i, j, k;

	for(k=1; k<R; k++) {
	    for(i=0; i<n[k]; i++) {
		for(j=0; j<n[k-1]; j++)
			PrintValue(L[k].w[i][j]);
		PrintValue(L[k].theta[i]);
	    }
	}
}

Test(int p)
{
	int i, k;

	BeginRow();
	PrintValue((double)p);
	Forward(p);
	k=0;
	for(i=0; i<n[k]; i++)
		PrintValue(L[k].output[i]);
	k=R-1;
	for(i=0; i<n[k]; i++)
		PrintValue(T[p].t[i]);
	for(i=0; i<n[k]; i++)
		PrintValue(L[k].output[i]);
	EndRow();
}

main()
{
	int p, loop;
	double E;

	BeginTable();
	Initialize();
	E=MSE();
	Alpha=ALPHA_MIN;
	for(loop=0; E>EPSILON; loop++) {
		if(loop%100==0) {
			BeginRow();
			PrintValue((double)loop);
			PrintValue(E);
			PrintNetwork();
			EndRow();
		}
		for(p=0; p<P; p++) {
			Forward(p);	/* 出力の計算 */
			Backward(p);	/* 誤差逆伝播 */
			Learn();	/* 各層の学習 */
		}
		E=MSE();
		if(Alpha<ALPHA_MAX)
			Alpha+=ALPHA_DELTA;
	}
	BeginRow();
	PrintValue((double)loop);
	PrintValue(E);
	PrintNetwork();
	EndRow();
	EndTable();

	/* テスト結果の出力 */
	BeginResultTable();
	for(p=0; p<P; p++)
		Test(p);
	EndTable();
}