/* 
   gcc -g simple-cnn_8x8.cpp -o 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"である
};

+2.0

double learning_input[][_image_Size][_image_Size] = {
  {
	// "A"の画像イメージ(1つ目)
	{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"の画像イメージ(1つ目)
	{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"の画像イメージ(1つ目)
	{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"の画像イメージ(1つ目)
	{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"の画像イメージ(2つ目)
	{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"の画像イメージ(2つ目)
	{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"の画像イメージ(2つ目)
	{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"の画像イメージ(2つ目)
	{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}
  }
};

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

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

// フィルタの初期化
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() ;
	  }
	}
  }
} 

// フィルタの適用 
double convolution_calculation(double filter[][_filter_Size]
							   ,double e[][_image_Size],int i,int j)
{
  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 ;
}


// 畳み込みの計算    
void convolution(double filter[][_filter_Size],
				 double e[][_image_Size],
				 double convolution_output[][_image_Size])
{
  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) ;
	}
  }
}


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

// プーリングの計算
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) ;
	}
  }
}


// 中間層の重みの初期化
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() ;
	}
  }
}

// 出力層の重みの初期化
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() ;
	}
  }
}

// 順方向の計算      
double forward(double weight_hidden[_hidden_Node_Number][_input_Node_Number + 1],
			   double weight_output[_hidden_Node_Number + 1],double hidden_input[],
			   double e[])
{
  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][_input_Node_Number] ;
	
	hidden_input[i] = sigmoid_function(u) ;
  }

  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) ;
}


// 出力層の重み学習
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[_hidden_Node_Number] += _alpha * (-1.0) * d ; //しきい値の学習
}

// 中間層の重み学習 
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++){  

	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++){
	  weight_hidden[j][i] += _alpha * e[i] * dj ;
	}
	weight_hidden[j][_input_Node_Number] += _alpha*(-1.0)*dj ;
  }
}

// 結果の出力      
void print_weight(double weight_hidden[_hidden_Node_Number][_input_Node_Number + 1]
				 ,double weight_output[_output_Node_Number][_hidden_Node_Number + 1])
{
  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") ;
}


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++){
		for(int n = 0; n < _pool_Outsize; n++){
		  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) ;

	for(int j = 0; j < _input_Node_Number; j++){
	  printf("%lf ",e[i][j]) ;
	}
	printf("\n") ;

	for(int j = _input_Node_Number; j < _input_Node_Number+_output_Node_Number; j++){
	  printf("%lf ",e[i][j]) ;
	}
	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 ;
}