/*
# mrcImageNeuralNetParticleLearning : $Revision$  
# $Date$ 
# Created by $Author$
# Usage : mrcImageNeuralNetParticleLearning
# Attention
#   $Loccker$
#  	$State$ 
#
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#define GLOBAL_DECLARATION
#include "../inc/config.h"

#undef DEBUG
#include "genUtil.h"
#include "mrcImage.h"
#include "Memory.h"
#include "String.h"

/*
#include "pickup.h"
*/

#include "NeuralNet.h"
#include "lmrcNeuralNet.h"

int
main(int argc, char* argv[]) 
{
	mrcImageNeuralNetParticleLearningInfo info;
	NeuralNetInfo linfo;
    mrcImage*   weightImage;
	mrcImage*   positive;
	mrcImage*   negative;
	double***   weight;
	int mode;
	int i,j;
	int x,y;
	int k;
	double E;
	int counter;
	double error;

	init0(&info);
    argCheck(&info, argc, argv);
    init1(&info);

    DEBUGPRINT("program start!\n");
  	// Information Set 		 	
	linfo.numOfLayer = info.flagOut + 1;		
	linfo.numOfCellOfLayer = (int*)memoryAllocate(sizeof(int)*linfo.numOfLayer, "in main");
	linfo.bias             = (int*)memoryAllocate(sizeof(int)*linfo.numOfLayer, "in main");
	fseek(info.fptOutList, 0L, SEEK_SET);
	for(i=1; i<linfo.numOfLayer; i++) {
		char s[1024];			
		fgets(s, 1023, info.fptOutList);
		linfo.numOfCellOfLayer[i] = stringGetNthIntegerData(s, 2, " ,\t");
		linfo.bias[i]             = stringGetNthIntegerData(s, 3, " ,\t");
	}
	linfo.numOfCellOfLayer[0] = 0;
	linfo.bias[0]             = info.Bias;

	positive = (mrcImage*)memoryAllocate(sizeof(mrcImage)*info.flagInPositive, "in main for positive");
	for(i=0; i<info.flagInPositive; i++) {
		mrcFileRead(&positive[i], info.InPositive[i], "in main for positive", 0);
		if(0==linfo.numOfCellOfLayer[0]) {
			linfo.numOfCellOfLayer[0] = positive[i].PixelPerImage;
		} else if(linfo.numOfCellOfLayer[0] != positive[i].PixelPerImage){
			fprintf(stderr, "The size of image %d is different from image 0\n", i);  			
			exit(EXIT_FAILURE);
		}
	}
	
	negative = (mrcImage*)memoryAllocate(sizeof(mrcImage)*info.flagInNegative, "in main for negative");
	for(i=0; i<info.flagInNegative; i++) {
		mrcFileRead(&negative[i], info.InNegative[i], "in main for negative", 0);
		if(linfo.numOfCellOfLayer[0] != negative[i].PixelPerImage){
			fprintf(stderr, "The size of image %d is different from image 0\n", i);  			
			exit(EXIT_FAILURE);
		}
	}
	linfo.numOfPositive = info.flagInPositive;
	linfo.numOfNegative = info.flagInNegative;
   	linfo.numOfCellOfLayer[0] = positive[0].HeaderN.x*positive[0].HeaderN.y;

    linfo.numOfSeparation     = 3;
	
    // this is the memoty of weight  
	weight = (double***)memoryAllocate(sizeof(double**)*(linfo.numOfLayer-1), "in  lmrcImageNeuralNetParticleLearning");
	for(i=0; i<linfo.numOfLayer-1; i++) {
		weight[i] = (double**)memoryAllocate(sizeof(double*)*(linfo.numOfCellOfLayer[i]+2), "in  lmrcImageNeuralNetParticleLearning");	
		for(j=0; j<linfo.numOfCellOfLayer[i]+2; j++) {
			weight[i][j]=(double*)memoryAllocate(sizeof(double)*(linfo.numOfCellOfLayer[i+1]+2), "lmrcImageNeuralNetParticleLearning");
		}
	}
	/*
		weight[l][m][n] : the weight between the m'th cell in the l'th layer and the n'th cell in the l+1'th layer   

			l+1      ..... n ..... 
                          /
                         / weight[l][m][n]
                        / 
			l        . m .........
 


			0        .............     : input-layer
 	*/


	// main's caluclation

	RandomToWeight(&linfo,weight);

	DEBUGPRINT("RandomToWeight OK!!\n");

	error = info.EPS*1e6;
	counter = 0;
    while(counter < info.MaxCounter &&  info.EPS < error) {
		error = 0;	
		linfo.Lambda              = info.Lambda;
		for(k=0;k<linfo.numOfPositive;k++) {
			E = lmrcImageNeuralNetParticleLearning(1, weight, &positive[k], &linfo, mode);
			error += E;
		}

		linfo.Lambda              = -1.0*info.Lambda;
		for(k=0;k<linfo.numOfNegative;k++) {
			E = lmrcImageNeuralNetParticleLearning(0, weight, &negative[k], &linfo, mode);
			error += E;
		}
		error = sqrt(error/(linfo.numOfPositive+linfo.numOfNegative));
		fprintf(stderr, "count: %4d error: %15.6f\n", counter, error);
		counter++;
	}

	weightImage = (mrcImage*)memoryAllocate(sizeof(mrcImage)*info.flagOut, "in main for weight");
	for(i=0; i<linfo.numOfLayer-1; i++) {		
		weightImage[i].HeaderN.x = linfo.numOfCellOfLayer[i]; 
		weightImage[i].HeaderN.y = linfo.numOfCellOfLayer[i+1]; 
		weightImage[i].HeaderN.z = 1;
		weightImage[i].HeaderMode = mrcFloatImage;
		mrcInit(&weightImage[i], NULL);

		for(x=0; x<weightImage[i].HeaderN.x; x++) {
		for(y=0; y<weightImage[i].HeaderN.y; y++) {
			mrcPixelDataSet(&weightImage[i], x, y, 0, weight[i][x][y], mrcPixelRePart);
		}
		}
	}

	for(i=0; i<info.flagOut; i++) {
		mrcFileWrite(&weightImage[i], info.Out[i], "in main for weight", 0);
	}

	return 0;
}

void
additionalUsage()
{
	fprintf(stderr, "----- Additional Usage -----\n");
}

void 
RandomToWeight(lmrcImageNeuralNetParticleLearningInfo* linfo, double*** weight)
{
	int l,m,n;
	srand((unsigned)time(NULL));

	for( l=0; l<linfo->numOfLayer-1 ; l++)
	{
		for(  m=0; m<linfo->numOfCellOfLayer[l]+2; m++)
		{
			for(  n=0; n<linfo->numOfCellOfLayer[l+1]+2; n++)
			{
				weight[l][m][n] = (double)rand() / (double)RAND_MAX  - 0.5 ;
				DEBUGPRINT4("firstweight[%d][%d][%d]=%f\n",l,m,n,weight[l][m][n]);
			}		
		}
	}
}

double
lmrcImageNeuralNetParticleLearning(
	int pattern,
	double*** weight,
	mrcImage* in,
	lmrcImageNeuralNetParticleLearningInfo* linfo,
	int mode)
{

	int       i,j;
	double**  cell;
	int       counter;
	int       index;
	float     x,y,z;
	double    error;

	cell = (double**)memoryAllocate(sizeof(double*)*(linfo->numOfLayer), "in  lmrcImageNeuralNetParticleLearning");
	for(i=0; i<linfo->numOfLayer; i++) {	
		cell[i] = (double*)memoryAllocate(sizeof(double)*(linfo->numOfCellOfLayer[i]+2), "in lmrcImageNeuralNetParticleLearning");
	}

	/*
		cell[l][n] : the value of the n'th cell in the l'th layer 	
			
			l		..... n .......
						  |
						cell[l][n]
	*/

	// Get A Sample Image
	index = 0;
	for(x=0; x<in->HeaderN.x; x++) {
	for(y=0; y<in->HeaderN.y; y++) {
	for(z=0; z<in->HeaderN.z; z++) {
		mrcPixelDataGet(in, x, y, z, &cell[0][index], mrcPixelRePart, mrcPixelHowNearest);
		cell[0][index] = sigmoidalFunction((cell[0][index] - in->HeaderAMean)/(in->HeaderAMax - in->HeaderAMin)*4);
		DEBUGPRINT4("x=%lf , y=%lf , cell[0][%d]=%f\n", x, y, index, cell[0][index]);
		index++;
	}
	}
	}
	cell[0][linfo->numOfCellOfLayer[0]]   = linfo->bias[0];
	cell[0][linfo->numOfCellOfLayer[0]+1] = in->HeaderAMean;
    linfo->numOfCellOfLayer[0] = index+1;

	counter = 0;
	if(pattern==1) {
		DEBUGPRINT("positive start\n");			
	} else if(pattern==0) {
		DEBUGPRINT("negative start\n");			
	}	

	OutputLayerCalculation(cell, weight, linfo, mode); 	
    Back_propagation(cell, weight, pattern, linfo, mode);
	OutputLayerCalculation(cell, weight, linfo, mode); 	
	
	// Calc. Error 
	error = 0;
    for(j=0 ; j<linfo->numOfCellOfLayer[linfo->numOfLayer-1] ; j++) {
		error += SQR(pattern - cell[linfo->numOfLayer-1][j]);
	}
	fprintf(stderr, "error=%f\n",error);

	for(i=0; i<linfo->numOfLayer; i++) {	
		free(cell[i]);
	}
	free(cell);
	return error;
}
				
void
OutputLayerCalculation(double** cell, double*** weight, lmrcImageNeuralNetParticleLearningInfo* linfo, int mode)
{
	int l,m,n;
	
    // Calclation
	DEBUGPRINT("Calculation start\n");
    for( l=0 ; l<linfo->numOfLayer-1 ; l++) {
		for( n=0 ; n<linfo->numOfCellOfLayer[l+1] ; n++ ) {
			cell[l+1][n] = 0;
			for( m=0 ; m<linfo->numOfCellOfLayer[l] ; m++ ) {
		   		cell[l+1][n]+=cell[l][m]*weight[l][m][n];
			}
			cell[l+1][n]=sigmoidalFunction(cell[l+1][n]);
			DEBUGPRINT6("l=%d , m=%d , n=%d ,cell[%d][%d]=%f\n",l,m,n,l+1,n,cell[l+1][n]);
		}
	}
}

double
sigmoidalFunction(double cell)
{
	return (1.0/(1.0+exp(cell)));
}

/*
	Estimate weights between layers using back-propagation

	Input: cell 
		   linfo

	Output: weight
*/

void Back_propagation
(double** cell, double*** weight, double Teacher, lmrcImageNeuralNetParticleLearningInfo* linfo, int mode)
{

	int i,j,l,m,n;
	double    sum;
	double**  squarederror;
    double*** weightadjust;
	
	DEBUGPRINT("Start Back Propagation\n");
	squarederror = (double**)memoryAllocate(sizeof(double**)*(linfo->numOfLayer), "in ImrcImageNeuralNetParticleLearning");
	for(i=0 ; i<linfo->numOfLayer ; i++) {
		squarederror[i] = (double*)memoryAllocate(sizeof(double)*(linfo->numOfCellOfLayer[i]+2), "in ImrcImageNeuralNetParticl			eLearning");
	}
	/*
	    squarederror[l][n] : the sqaurederror of n'th cell in the l'th layer

		l      . . . . n . . . .
		               |
				 squarederror[l][n]

	*/			
	
    
	weightadjust = (double***)memoryAllocate(sizeof(double**)*(linfo->numOfLayer-1), "in  lmrcImageNeuralNetParticleLearning");
	for(i=0; i<linfo->numOfLayer-1; i++) {
		weightadjust[i] = (double**)memoryAllocate(sizeof(double*)*(linfo->numOfCellOfLayer[i]+2), "in  lmrcImageNeuralNetParticleLearning");	
		for(j=0; j<linfo->numOfCellOfLayer[i]; j++) {
			weightadjust[i][j]=(double*)memoryAllocate(sizeof(double)*(linfo->numOfCellOfLayer[i+1]+2), "lmrcImageNeuralNetParticleLearning");
		}
	}
    /*
	    weightadjust[l][m][n] : the weightadjust adjust weight of l'th layer between m'th cell in the l'th layer and n'th cell in the l-1'th layer

		l+1    . . . . n . . . .
                      / 
					 /  weightadjust[l][m][n]
					/ 
		l      . . m . . . . . .
    */


    /*  OutputLayer' caluclation  */

	for(i=0 ; i<linfo->numOfCellOfLayer[linfo->numOfLayer-1] ; i++) {
		squarederror[linfo->numOfLayer-1][i] = pow(cell[linfo->numOfLayer-1][i] - Teacher, 2.0 );
		DEBUGPRINT1("Top: %lf\n", squarederror[linfo->numOfLayer-1][i]);
    }
		
	for( m=0 ; m<linfo->numOfCellOfLayer[linfo->numOfLayer-2] ; m++) {
		for( n=0 ; n<linfo->numOfCellOfLayer[linfo->numOfLayer-1] ; n++) {
			weightadjust[linfo->numOfLayer-2][m][n] = -1.0 * linfo->Lambda * squarederror[linfo->numOfLayer-1][n] * cell[linfo->numOfLayer-2][m] ;
			DEBUGPRINT3("Adjust from Top: wa %lf cell %lf sqr %lf\n", 
					weightadjust[linfo->numOfLayer-2][m][n],  
					cell[linfo->numOfLayer-2][m],
					 squarederror[linfo->numOfLayer-1][n]);
		}	
	}	

		
	for( l=linfo->numOfLayer-2 ; l>0 ; l--) {
		for( m=0 ; m<linfo->numOfCellOfLayer[l] ; m++) {   
		    sum=0;
			for( n=0 ; n<linfo->numOfCellOfLayer[l+1] ; n++) {
	       		sum += weight[l][m][n] * squarederror[l+1][n] ;
			}	
			squarederror[l][m] = sum * (cell[l][m] * (1.0 - cell[l][m]));
		}
	
		for( m=0 ; m<linfo->numOfCellOfLayer[l-1] ; m++) {
			for( n=0 ; n<linfo->numOfCellOfLayer[l] ; n++) {
				weightadjust[l-1][m][n] = -1.0 * linfo->Lambda * squarederror[l][n] * cell[l-1][m] ;
			    DEBUGPRINT4("weightadjust[%d][%d][%d]=%f\n",l-1,m,n,weightadjust[l-1][m][n]);
			}
		}	
	}		

	DEBUGPRINT("weightadjust finish\n");
    for( l=0 ; l<linfo->numOfLayer-1 ; l++) {
    	for( m=0 ; m<linfo->numOfCellOfLayer[l] ; m++) {
			for ( n=0 ; n<linfo->numOfCellOfLayer[l+1] ; n++) {
				weight[l][m][n] += weightadjust[l][m][n];
				DEBUGPRINT4("weight[%d][%d][%d]=%f\n",l,m,n,weight[l][m][n]);
			}
		}
	}
	DEBUGPRINT("weight finish\n");

	/* release heap-areas */
	for(i=0 ; i<linfo->numOfLayer ; i++) {
		free(squarederror[i]);
	}
	free(squarederror);

	for(i=0; i<linfo->numOfLayer-1; i++) {
		for(j=0; j<linfo->numOfCellOfLayer[i]; j++) {
			free(weightadjust[i][j]);
		}
		free(weightadjust[i]);
	}
	free(weightadjust);

	DEBUGPRINT("End Back Propagation\n");
}


