#include "math.h"
#include "stdlib.h"
#include "stdio.h"
#include "conio.h"
#include "iostream.h"
#include <time.h>

#define I 9     // 輸入層神經元個數 (decision variables)
#define H 4     // 隱含層神經元個數
#define O 1     // 輸出層神經元個數
#define D 8  // 訓練樣本個數

double x[D][I+1],y[D][O+1],h[D][H+1],out_y[D][O+1],wh[H+1][I+1],wo[O+1][H+1];
double ETA=1.1;
void   Input_Output(double x[D][I+1], double y[D][O+1]);
double Randw(double a, double b);
void   Init_W(void);//初始化權向量
double Sigmoid(double a);
void   Forward();
void   Backward();
double Err();
void   Out_W();
//void   BP();

void   Init_xy(void)
{
	double xx[D][I+1]={{1,1,0,0,1,0,0,1,1,1},{1,0,0,1,0,0,1,1,1,1},{1,1,1,1,0,0,1,0,0,1},
	{1,1,1,1,1,0,0,1,0,0},{1,1,1,1,0,1,0,0,1,0},{1,1,0,0,1,1,1,1,0,0},{1,0,1,0,0,1,0,1,1,1},
	{1,0,0,1,1,1,1,0,0,1}};
	for(int i=0;i<D;i++)
		for(int j=0;j<=I;j++)
			x[i][j]=xx[i][j];
	y[0][1]=0;	y[1][1]=0;	y[2][1]=0;	y[3][1]=0;
	y[4][1]=1;	y[5][1]=1;	y[6][1]=1;	y[7][1]=1;
}

double Randw(double a, double b) // Uniform Distribution
{
  double y;
  if(a>b) {
	printf("\nThe first parameter should be less than the second!");
	exit(1);
  }
  y = (double)rand()/(RAND_MAX);
  return (a+(b-a)*y); 
}

void   Init_W(void)
{
	int i,j;
	for(i=1; i<=H; i++)
		for(j=0; j<=I; j++)
			wh[i][j]= Randw(-0.3,0.3);
		for(i=1; i<=O; i++)
			for(j=0; j<=H; j++)
				wo[i][j]= Randw(-0.3,0.3);
}

double Sigmoid(double a)//sigmoid函數***************************************OK:)
{
	double b;
	a=(-1)*a;
	b=1./(1+exp(a));
	return b;
}

void Forward()
{
	int i,j,k;
	for( k=0; k<D; k++)
	{	/*計算隱含層輸出*/
		h[k][0]=1.;//
		for( i=1; i<=H; i++ )
		{
			h[k][i] = 0.;
			for( j=0; j<=I; j++ )
				h[k][i] += wh[i][j] * x[k][j];
			h[k][i] = Sigmoid(h[k][i]);//Sigmoid函數
		}
		/*計算輸出層輸出*/
		for( i=1; i<=O; i++ )
		{
			out_y[k][i] = 0.;
			for( j=0; j<=H; j++)
				out_y[k][i] += wo[i][j] * h[k][j];
			out_y[k][i] = Sigmoid(out_y[k][i]);
		}
	}
}

void Backward()
{
	int i,j,k;
	double temp;
	double DETA_WH[H+1][I+1],DETA_WO[O+1][H+1];
	double E_WO[O+1][H+1],E_WH[H+1][I+1],eo[D][O+1],eh[D][H+1];
	//差分置零
	for( i=1; i<=H; i++ )
		for( j=0; j<=I; j++ )
			DETA_WH[i][j] = 0.;
	for( i=1; i<=O; i++ )
		for( j=0; j<=H; j++ )
			DETA_WO[i][j] = 0.;
	//調整隱含層權向量
	for( i=1; i<=O; i++ )
	{	
		for( k=0; k<D; k++)
			eo[k][i] = (1-out_y[k][i])*out_y[k][i]*(out_y[k][i] - y[k][i])/D;
		for( j=0; j<=H; j++ )
		{
			E_WO[i][j]=0.;
			for(k=0;k<D;k++)
				E_WO[i][j] += eo[k][i] * h[k][j];
			DETA_WO[i][j] = ETA*E_WO[i][j];
			wo[i][j] -= DETA_WO[i][j];
		}
	}
	//調整隱含層權向量
	for( j=1; j<=H; j++)
	{
		for( k=0; k<D; k++)
		{
			for( i=1,temp=0.; i<=O; i++)
				temp += eo[k][i] * wo[i][j];
			eh[k][j] = (1-h[k][j])*h[k][j]*temp;
		}
		for(i=0; i<=I; i++)
		{
			E_WH[j][i]=0.;
			for( k=0; k<D; k++)
				E_WH[j][i] += eh[k][j]*x[k][i];
			DETA_WH[j][i] = ETA * E_WH[j][i];
			wh[j][i] -= DETA_WH[j][i];
		}
	}
}

double   Err()
{
	int i,k;
	double e[D][O+1];
	double E=0.;
	for(k=0;k<D;k++)
	{
		for( i=1; i<=O; i++ )
		{	
			e[k][i] = y[k][i] - out_y[k][i];
//			EE += fabs(e[k][i])/D;//平均誤差??????????????????
			E += 0.5*e[k][i]*e[k][i]/D;//平方和的1/2
		}
	}
	return E;
}

void Out_W()
{
	int i,j;
	for(i=1; i<=H; i++)
	{
		for(j=0; j<=I; j++)
			cout<<wh[i][j]<<'\t';
		cout<<'\n';
	}
	for(i=1; i<=O; i++)
	{
		for(j=0; j<=H; j++)
			cout<<wo[i][j]<<'\t';
		cout<<'\n';
	}
}

/////////////////////主函數////////////////////////
void main()
{
	int i,step=0;
	double E,Emax=0.000001;
	double t1,t2;
	Init_xy();
	Init_W();
	t1=(double)clock();
	do{
		step++;
		Forward();
		E=Err();
		if(step%1000==0) cout<<"\nBP_"<<step<<"    "<<E<<endl;
		if (E<=Emax)
			break;
		else
			Backward();
	}while(step<50000000);
	t2=(double)clock();
	cout<<"E="<<E<<'\n';
	cout<<"step="<<step<<endl;
	for(i=0;i<D;i++)
		cout<<out_y[i][1]<<'\n';
	cout<<"time:"<<(t2-t1)/CLK_TCK<<"秒"<<endl;
}