/* 
   gcc -g destructive_simple-cnn_8x8.cpp -o destructive_simple-cnn_8x8
*/

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


// 畳み込み
static const int _Image_Size = 8; // 入力画像の1辺のピクセル数
static const int _Filter_Size = 3; //畳み込みフィルタのサイズ
static const int _Filter_Number = 2; // フィルタの個数
static const int _Pool_Outsize = 3; //プーリング層の出力のサイズ

// 全結合
static const int _Input_Node_Number = 18;   /*入力層のセル数*/
static const int _hidden_node_number = 6;  /*中間層のセル数*/
static const int _output_node_number =4;   /*出力層のセル数*/
static const int _alpha =  1;    /*学習係数*/
static const int _max_data_number = 100;   /*学習データの最大個数*/
static const int _big_number = 100;  /*誤差の初期値*/
static const double _limit =0.001; /*誤差の上限値*/
static const int _seed = 199;  /*乱数のシード*/

//  教師データ

double learning_output[][_output_node_number] ={
  {1,0,0,0}, 
  {0,1,0,0}, // "B"である
  {0,0,1,0}, // "C"である
  {0,0,0,1}, // "D"である
  {1,0,0,0}, // "A"である
  {0,1,0,0}, // "B"である
  {0,0,1,0}, // "C"である
  {0,0,0,1}  // "D"である
};

double learning_input[][_Image_Size][_Image_Size] = {
  {
	// "A"の画像イメージ
	{0,0,0,1,1,0,0,0},
	{0,0,0,1,1,0,0,0},
	{0,0,1,1,1,1,0,0},
	{0,0,1,0,1,0,0,0},
	{0,0,1,0,1,0,0,0},
	{0,1,1,1,1,1,1,0},
	{0,1,0,0,0,0,1,0},
	{0,1,0,0,0,0,1,0}
  },

  {
	// "B"の画像イメージ
	{0,1,1,1,1,0,0,0},
	{0,1,0,0,0,1,0,0},
	{0,1,0,0,0,1,0,0},
	{0,1,1,1,1,1,0,0},
	{0,1,0,0,0,1,0,0},
	{0,1,0,0,0,0,1,0},
	{0,1,0,0,0,1,0,0},
	{0,1,1,1,1,0,0,0}
  },

  {
	// "C"の画像イメージ
	{0,0,0,1,1,1,0,0},
	{0,0,1,0,0,0,1,0},
	{0,1,0,0,0,0,0,0},
	{0,1,0,0,0,0,0,0},
	{0,1,0,0,0,0,0,0},
	{0,1,0,0,0,0,0,0},
	{0,0,1,0,0,0,1,0},
	{0,0,1,1,1,1,0,0}
  },

  {
	// "D"の画像イメージ
	{0,1,1,1,1,1,0,0},
	{0,1,0,0,0,1,1,0},
	{0,1,0,0,0,0,1,0},
	{0,1,0,0,0,0,1,0},
	{0,1,0,0,0,0,1,0},
	{0,1,0,0,0,0,1,0},
	{0,1,0,0,0,1,1,1},
	{0,1,1,1,1,1,0,0}
  },

  {
	// "A"の画像イメージ
	{0,0,0,0,0,1,0,0},
	{0,0,0,0,1,1,0,0},
	{0,0,0,1,1,0,1,0},
	{0,0,0,1,0,0,1,0},
	{0,0,1,0,0,0,1,0},
	{0,0,1,1,1,1,1,0},
	{0,1,0,0,0,0,1,0},
	{0,0,0,0,0,0,1,0}
  },

  {
	// "B"の画像イメージ
	{0,0,0,1,1,1,1,0},
	{0,0,0,1,0,0,0,1},
	{0,0,0,1,0,0,0,1},
	{0,0,0,1,1,1,1,1},
	{0,0,0,1,0,0,0,1},
	{0,0,0,1,0,0,0,1},
	{0,0,0,1,1,0,0,1},
	{0,0,0,1,1,1,1,0}
  },

  {
	// "C"の画像イメージ
	{0,1,1,1,1,0,0,0},
	{1,0,0,0,0,0,0,0},
	{1,0,0,0,0,0,0,0},
	{1,0,0,0,0,0,0,0},
	{1,0,0,0,0,0,0,0},
	{0,1,0,0,1,0,0,0},
	{0,0,1,1,0,0,0,0},
	{0,0,0,0,0,0,0,0}
  },

  {
	// "D"の画像イメージ
	{0,0,0,1,1,0,0,0},
	{0,0,1,0,0,1,0,0},
	{0,0,1,0,0,0,1,0},
	{0,0,1,0,0,0,1,0},
	{0,1,0,0,0,0,1,0},
	{0,1,0,0,0,1,1,0},
	{0,1,1,1,1,1,0,0},
	{0,0,0,0,0,0,0,0}
  }
};

/************************
 drand()関数           
-1から1の間の乱数を生成
************************/
double drand(void)
{
  double rndno ;/*生成した乱数*/
  
  while((rndno = (double)rand()/RAND_MAX) == 1.0) ;
  rndno = rndno * 2 -1 ;/*-1から1の間の乱数を生成*/
  return rndno;
}

/*********************
init_filter()関数
ィルタの初期化
*********************/
void init_filter(double filter[_Filter_Number][_Filter_Size][_Filter_Size])
{
  for(int i = 0; i < _Filter_Number; i++){
	for(int j = 0; j < _Filter_Size; j++){
	  for(int k = 0; k < _Filter_Size;k++){
		filter[i][j][k]=drand() ;
	  }
	}
  }
} 


void convolution(double filter[_Filter_Size][_Filter_Size]
		  ,double e[][_Image_Size]
		  ,double convolution_output[][_Image_Size]) ; /*畳み込みの計算*/

double convolution_calculation(double filter[][_Filter_Size]
				,double e[][_Image_Size],int i,int j) ;/*  フィルタの適用  */

void pooling(double convolution_output[][_Image_Size],double pooling_output[][_Pool_Outsize]) ; 
/*プーリングの計算*/ 

double pooling_calculation(double convolution_output[][_Image_Size]
				   ,int x,int y) ;/* 平均値プーリング */

//全結合層関連
void init_weight_hidden_layer(double weight_hidden[_hidden_node_number][_Input_Node_Number + 1]) ;

//中間層の重みの初期化
void init_weight_output_layer(double weight_output[_output_node_number][_hidden_node_number + 1]) ;

/*出力層の重みの初期化*/
double forward(double weight_hidden[_hidden_node_number][_Input_Node_Number + 1]
			   ,double [_hidden_node_number + 1],double hidden_input[]
			   ,double e[_Input_Node_Number+_output_node_number]) ; /*順方向の計算*/

void output_layer_learning(double weight_output[_hidden_node_number + 1],double hidden_input[]
			,double e[_Input_Node_Number + _output_node_number],double o,int k) ;

/*出力層の重みの学習*/
void hidden_layer_learning(double weight_hidden[_hidden_node_number][_Input_Node_Number + 1]
			,double weight_output[_hidden_node_number + 1],double hidden_input[]
			,double e[_Input_Node_Number+_output_node_number],double o,int k) ;

/*中間層の重みの学習*/
void print_weight(double weight_hidden[_hidden_node_number][_Input_Node_Number + 1]
				 ,double weight_output[_output_node_number][_hidden_node_number + 1]) ; /*結果の出力*/

double sigmoid_function(double u) ; /*シグモイド関数*/





 
/*******************/
/*    main()関数   */
/*******************/
int main()
{
  /*畳み込み演算関連*/
  double filter[_Filter_Number][_Filter_Size][_Filter_Size] ;/*畳み込みフィルタ*/

  double convolution_output[_Image_Size][_Image_Size] ;/*畳み込み出力*/
  double pooling_output[_Pool_Outsize][_Pool_Outsize] ;/*出力データ*/
  
  /*全結合層関連*/
  double weight_hidden[_hidden_node_number][_Input_Node_Number + 1] ;/*中間層の重み*/

  double weight_output[_output_node_number][_hidden_node_number + 1] ;/*出力層の重み*/

  double e[_max_data_number][_Input_Node_Number + _output_node_number] ;/*学習データセット*/

  double hidden_input[_hidden_node_number + 1] ;/*中間層の出力*/

  double o[_output_node_number]  ;/*出力*/

  double err = _big_number ;/*誤差の評価*/

  //int i,j,m,n ;/*繰り返しの制御*/

  int n_of_e ;/*学習データの個数*/
  int count = 0 ;/*繰り返し回数のカウンタ*/
  
  /*乱数の初期化*/
  srand(_seed) ;
  
  /*畳み込みフィルタの初期化*/
  init_filter(filter) ;
  
  /*重みの初期化*/
  init_weight_hidden_layer(weight_hidden) ;

  init_weight_output_layer(weight_output) ;

  print_weight(weight_hidden,weight_output) ;
  
  /*学習データの読み込み*/
  //  n_of_e=getdata(image,t) ;
  n_of_e = 8;

  printf("学習データの個数:%d\n",n_of_e) ;
  
  /*畳み込み処理*/
  for(int i = 0; i < n_of_e; i++){/*学習データ毎の繰り返し*/
	for(int j = 0; j < _Filter_Number; j++){/*フィルタ毎の繰り返し*/

	  /*畳み込みの計算*/
	  convolution(filter[j], learning_input[i], convolution_output) ;

	  /*プーリングの計算*/
	  pooling(convolution_output,pooling_output) ;

	  /*プーリング出力を全結合層の入力へコピー*/
	  for(int m = 0; m<_Pool_Outsize; m++){   // _Pool_Outsize は "3"
		for(int n = 0; n < _Pool_Outsize; n++){ // _Pool_Outsize は "3"
		  e[i][j * _Pool_Outsize * _Pool_Outsize + _Pool_Outsize * m + n] 
			= pooling_output[m][n] ;
		}
	  }

	  for(int m = 0; m < _output_node_number; m++){
		e[i][_Pool_Outsize*_Pool_Outsize*_Filter_Number+m] 
		  = learning_output[i][m] ;/*教師データ*/
	  }
	}
  }
  
  /*学習*/
  while( err > _limit){
	/*i個の出力層に対応*/
	for(int i = 0; i < _output_node_number; i++){
	  err=0.0 ;
	  for(int j = 0; j < n_of_e; j++){

		/*順方向の計算*/
		o[i] = forward(weight_hidden, weight_output[i], hidden_input, e[j]) ;

		/*出力層の重みの調整*/
		output_layer_learning(weight_output[i], hidden_input, e[j], o[i], i) ;

		/*中間層の重みの調整*/
		hidden_layer_learning(weight_hidden, weight_output[i], hidden_input, e[j], o[i], i) ;

		/*誤差の積算*/
		err += (o[i] - e[j][_Input_Node_Number+i]) * (o[i] - e[j][_Input_Node_Number+i]);
	  }

	  count++ ;
	  /*誤差の出力*/
	  
	  printf("%d\t%lf\n",count,err) ;
	}
  }/*学習終了*/
  
  /*重みの出力*/
  print_weight(weight_hidden,weight_output) ;
  
  /*学習データに対する出力*/
  for(int i = 0; i < n_of_e; i++){
	printf("%d\n",i) ;


	printf("学習データ\n");
	for(int j = 0; j < _Input_Node_Number; j++){
	  printf("%lf ",e[i][j]) ;/*学習データ*/
	}

	printf("\n") ;
	printf("教師データ\n");
	for(int j = _Input_Node_Number; j < _Input_Node_Number+_output_node_number; j++){/*教師データ*/
	  printf("%lf ",e[i][j]) ;
	}
	printf("\n") ;

	printf("ネットワーク出力\n");
	for(int j = 0; j < _output_node_number; j++){/*ネットワーク出力*/
	  printf("%lf ",forward(weight_hidden, weight_output[j], hidden_input, e[i])) ;
	}
	printf("\n") ;
  }
  
  return 0 ;
}


/*********************
init_weight_hidden_layer()関数
中間層の重みの初期化
*********************/
void init_weight_hidden_layer(double weight_hidden[_hidden_node_number][_Input_Node_Number + 1])
{
  /*　乱数による重みの初期化*/
  for(int i = 0; i < _hidden_node_number; i++){
	for(int j = 0; j < _Input_Node_Number + 1; j++){
	  weight_hidden[i][j] = drand() ;
	}
  }
}

/**********************
init_weight_output_layer()関数
出力層の重みの初期化
**********************/
void init_weight_output_layer(double weight_output[_output_node_number][_hidden_node_number + 1])
{
  /*　乱数による重みの初期化*/
  for(int i = 0; i < _output_node_number; i++){
	for(int j = 0; j < _hidden_node_number + 1; j++){
	  weight_output[i][j]=drand() ;
	}
  }
}

/*********************
convolution()関数        
畳み込みの計算    
*********************/
void convolution(double filter[][_Filter_Size],
				 double e[][_Image_Size],
				 double convolution_output[][_Image_Size])
{
  //  int i = 0,j = 0 ;/*繰り返しの制御用*/
  int startpoint = _Filter_Size / 2 ;/*畳み込み範囲の下限*/
  
  for(int i = startpoint; i < _Image_Size-startpoint; i++){
	for(int j = startpoint; j < _Image_Size-startpoint; j++){
	  convolution_output[i][j] = convolution_calculation(filter,e,i,j) ;
	}
  }
}

/**********************/
/*  convolution_calculation()関数    */
/*  フィルタの適用    */
/**********************/
double convolution_calculation(double filter[][_Filter_Size]
				,double e[][_Image_Size],int i,int j)
{
  // int m,n ;/*繰り返しの制御用*/
  double sum = 0 ;/*和の値*/
  
  for(int m = 0; m < _Filter_Size; m++){
	for(int n = 0; n < _Filter_Size; n++){
	  sum += e[i - _Filter_Size / 2 + m][j - _Filter_Size / 2 + n] * filter[m][n];
	}
  }
  
  return sum ;
}

/**********************/
/*  pooling()関数        */
/* プーリングの計算   */
/**********************/
void pooling(double convolution_output[][_Image_Size],
			 double pooling_output[][_Pool_Outsize]) 
{

  for(int i = 0; i < _Pool_Outsize; i++){
	for(int j = 0; j < _Pool_Outsize; j++){
	  pooling_output[i][j] = pooling_calculation(convolution_output, i * 2 + 1,j * 2 + 1) ;
	}
  }
}

/*********************
pooling_calculation()関数
平均値プーリング   
*********************/
double pooling_calculation(double convolution_output[][_Image_Size], int x, int y)
{
  double ave=0.0 ;/*平均値*/
  
  for(int m = x; m <= x + 1 ; m++){
	for(int n = y; n <= y + 1; n++){
	  ave += convolution_output[m][n] ;
	}
  }
  
  return ave/4.0 ;
}

/*********************
forward()関数    
順方向の計算      
*********************/
double forward(double weight_hidden[_hidden_node_number][_Input_Node_Number + 1],
			   double weight_output[_hidden_node_number + 1],double hidden_input[],
			   double e[])
{
  /*hidden_input計算*/
  for(int i = 0; i < _hidden_node_number; i++){
	/*重み付き和の計算*/
	double u = 0 ;/*重み付き和を求める*/ 

	for(int j = 0; j < _Input_Node_Number; j++){
	  u += e[j] * weight_hidden[i][j] ;
	}

	//u-=weight_hidden[i][j] ;/*しきい値の処理*/
	u -= weight_hidden[i][_Input_Node_Number] ;/*しきい値の処理*/
	
	hidden_input[i] = sigmoid_function(u) ;
  }
  /*出力outputの計算*/
  double output = 0.0 ;

  for(int i = 0; i < _hidden_node_number; i++){
	output += hidden_input[i]*weight_output[i] ;
  }

  output -= weight_output[_hidden_node_number] ;/*しきい値の処理*/
  
  return sigmoid_function(output) ;
}

/*********************
output_layer_learning()関数
出力層の重み学習
*********************/
void output_layer_learning(double weight_output[_hidden_node_number + 1]
			,double hidden_input[],double e[],double o,int k)
{
  /*重み計算に利用*/
  double d = ( e[_Input_Node_Number+k] - o) * o * (1-o) ;/*誤差の計算*/
  for(int i = 0; i < _hidden_node_number; i++){
	weight_output[i]+=_alpha*hidden_input[i]*d ;/*重みの学習*/
  }

  //weight_output[i]+=_alpha*(-1.0)*d ;/*しきい値の学習*/
  weight_output[_hidden_node_number] += _alpha * (-1.0) * d ;/*しきい値の学習*/
}

/**********************/
/*  hidden_layer_learning()関数      */
/*  中間層の重み学習  */
/**********************/
void hidden_layer_learning(double weight_hidden[_hidden_node_number][_Input_Node_Number + 1],double weight_output[_hidden_node_number + 1]
			,double hidden_input[],double e[],double output,int k)
{
  for(int j = 0; j < _hidden_node_number; j++){/*中間層の各セルjを対象*/
	/*中間層の重み計算に利用*/
	double dj = hidden_input[j] * (1 - hidden_input[j]) * weight_output[j] * (e[_Input_Node_Number+k] - output) * output * (1-output) ;
	
	for(int i = 0; i < _Input_Node_Number; i++){/*i番目の重みを処理*/
	  weight_hidden[j][i] += _alpha * e[i] * dj ;
	}

	//	weight_hidden[j][i] += _alpha*(-1.0)*dj ;/*しきい値の学習*/
	weight_hidden[j][_Input_Node_Number] += _alpha*(-1.0)*dj ;/*しきい値の学習*/
  }
}

/**********************/
/*  print_weight()関数 */
/*   結果の出力       */
/**********************/
void print_weight(double weight_hidden[_hidden_node_number][_Input_Node_Number + 1]
				 ,double weight_output[_output_node_number][_hidden_node_number + 1])
{
  //int i,j ;/*繰り返しの制御*/
  
  for(int i = 0; i < _hidden_node_number; i++){
	for(int j = 0; j < _Input_Node_Number + 1; j++){
	  printf("%lf ",weight_hidden[i][j]) ;
	}
  }
  printf("\n") ;

  for(int i = 0; i < _output_node_number; i++){
	for(int j = 0; j < _hidden_node_number + 1; j++){
	  printf("%lf ",weight_output[i][j]) ;
	}
  }
  printf("\n") ;
}

/*******************/
/* s()関数         */
/* シグモイド関数  */
/*******************/
double sigmoid_function(double u)
{
  return 1.0/(1.0+exp(-u)) ;
}