using System;
using System.Windows.Forms;

namespace aiNnFingers
{
	public class BPN
	{
		public BPN()
		{

		}

		public RichTextBox richTextBox;

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

		//typedef int           BOOL;
		//typedef int           INT;
		//typedef double        REAL;

		//#define FALSE         0
		//#define TRUE          1
		//#define NOT           !
		//#define AND           &&
		//#define OR            ||

		//#define MIN_REAL      -HUGE_VAL
		const double MIN_REAL = double.MinValue;
		//#define MAX_REAL      +HUGE_VAL
		const double MAX_REAL = double.MaxValue;
		//#define MIN(x,y)      ((x)<(y) ? (x) : (y))
		//#define MAX(x,y)      ((x)>(y) ? (x) : (y))

		//#define LO            0.1
		const double LO = 0.1;
		//#define HI            0.9
		const double HI = 0.9;
		//#define BIAS          1
		const double BIAS = 1.0;

		//#define sqr(x)        ((x)*(x))


		public struct LAYER 
		{                     /* A LAYER OF A NET:                     */
			public int           Units;         /* - number of units in this layer       */
			public double[]         Output;        /* - output of ith unit                  */
			public double[]         Error;         /* - error term of ith unit              */
			public double[][]        Weight;        /* - connection weights to ith unit      */
			public double[][]        WeightSave;    /* - saved weights for stopped training  */
			public double[][]        dWeight;       /* - last weight deltas for momentum     */
		}

		public struct NET 
		{                     /* A NET:                                */
			public LAYER[]       Layer;         /* - layers of this net                  */
			public LAYER        InputLayer;    /* - input layer                         */
			public LAYER        OutputLayer;   /* - output layer                        */
			public double          Alpha;         /* - momentum factor                     */
			public double          Eta;           /* - learning rate                       */
			public double          Gain;          /* - gain of sigmoid function            */
			public double          Error;         /* - total net error                     */
		}


		/******************************************************************************
				R A N D O M S   D R A W N   F R O M   D I S T R I B U T I O N S
		******************************************************************************/


		Random r;
		const int RAND_MAX = Int32.MaxValue;
		void InitializeRandoms()
		{
			//srand(4711);
			r = new Random(4711);
		}


		int RandomEqualINT(int Low, int High)
		{
			//return rand() % (High-Low+1) + Low;
			int retVal = r.Next() % (High-Low+1) + Low;
			return retVal;
		}      


		double RandomEqualREAL(double Low, double High)
		{
			double retVal = ((double) r.Next() / (double) RAND_MAX) * (High-Low) + Low; //268 times
			return retVal;
		}              

		/******************************************************************************
					A P P L I C A T I O N - S P E C I F I C   C O D E
		******************************************************************************/


		//#define NUM_LAYERS    3
		const int NUM_LAYERS = 3;
		//#define N             30
		const int N = 30;
		//#define M             1
		const int M = 1;
		//int                   Units[NUM_LAYERS] = {N, 10, M};
		int [] Units = new int [3] {N, 10, M};

		//#define FIRST_YEAR    1700
		const int FIRST_YEAR = 1700;
		//#define NUM_YEARS     280
		const int NUM_YEARS = 280;

		//#define TRAIN_LWB     (N)
		const int TRAIN_LWB = N;
		//#define TRAIN_UPB     (179)
		const int TRAIN_UPB = 179;
		//#define TRAIN_YEARS   (TRAIN_UPB - TRAIN_LWB + 1)
		const int TRAIN_YEARS = TRAIN_UPB - TRAIN_LWB + 1;
		//#define TEST_LWB      (180)
		const int TEST_LWB = 180;
		//#define TEST_UPB      (259)
		const int TEST_UPB = 259;
		//#define TEST_YEARS    (TEST_UPB - TEST_LWB + 1)
		const int TEST_YEARS = TEST_UPB - TEST_LWB + 1;
		//#define EVAL_LWB      (260)
		const int EVAL_LWB = 260;
		//#define EVAL_UPB      (NUM_YEARS - 1)
		const int EVAL_UPB = NUM_YEARS - 1;
		//#define EVAL_YEARS    (EVAL_UPB - EVAL_LWB + 1)
		const int EVAL_YEARS = EVAL_UPB - EVAL_LWB + 1;

		double []                  Sunspots_ = new double[NUM_YEARS];
		double []                 Sunspots = new double[NUM_YEARS] {

								0.0262,  0.0575,  0.0837,  0.1203,  0.1883,  0.3033,  
								0.1517,  0.1046,  0.0523,  0.0418,  0.0157,  0.0000,  
								0.0000,  0.0105,  0.0575,  0.1412,  0.2458,  0.3295,  
								0.3138,  0.2040,  0.1464,  0.1360,  0.1151,  0.0575,  
								0.1098,  0.2092,  0.4079,  0.6381,  0.5387,  0.3818,  
								0.2458,  0.1831,  0.0575,  0.0262,  0.0837,  0.1778,  
								0.3661,  0.4236,  0.5805,  0.5282,  0.3818,  0.2092,  
								0.1046,  0.0837,  0.0262,  0.0575,  0.1151,  0.2092,  
								0.3138,  0.4231,  0.4362,  0.2495,  0.2500,  0.1606,  
								0.0638,  0.0502,  0.0534,  0.1700,  0.2489,  0.2824,  
								0.3290,  0.4493,  0.3201,  0.2359,  0.1904,  0.1093,  
								0.0596,  0.1977,  0.3651,  0.5549,  0.5272,  0.4268,  
								0.3478,  0.1820,  0.1600,  0.0366,  0.1036,  0.4838,  
								0.8075,  0.6585,  0.4435,  0.3562,  0.2014,  0.1192,  
								0.0534,  0.1260,  0.4336,  0.6904,  0.6846,  0.6177,  
								0.4702,  0.3483,  0.3138,  0.2453,  0.2144,  0.1114,  
								0.0837,  0.0335,  0.0214,  0.0356,  0.0758,  0.1778,  
								0.2354,  0.2254,  0.2484,  0.2207,  0.1470,  0.0528,  
								0.0424,  0.0131,  0.0000,  0.0073,  0.0262,  0.0638,  
								0.0727,  0.1851,  0.2395,  0.2150,  0.1574,  0.1250,  
								0.0816,  0.0345,  0.0209,  0.0094,  0.0445,  0.0868,  
								0.1898,  0.2594,  0.3358,  0.3504,  0.3708,  0.2500,  
								0.1438,  0.0445,  0.0690,  0.2976,  0.6354,  0.7233,  
								0.5397,  0.4482,  0.3379,  0.1919,  0.1266,  0.0560,  
								0.0785,  0.2097,  0.3216,  0.5152,  0.6522,  0.5036,  
								0.3483,  0.3373,  0.2829,  0.2040,  0.1077,  0.0350,  
								0.0225,  0.1187,  0.2866,  0.4906,  0.5010,  0.4038,  
								0.3091,  0.2301,  0.2458,  0.1595,  0.0853,  0.0382,  
								0.1966,  0.3870,  0.7270,  0.5816,  0.5314,  0.3462,  
								0.2338,  0.0889,  0.0591,  0.0649,  0.0178,  0.0314,  
								0.1689,  0.2840,  0.3122,  0.3332,  0.3321,  0.2730,  
								0.1328,  0.0685,  0.0356,  0.0330,  0.0371,  0.1862,  
								0.3818,  0.4451,  0.4079,  0.3347,  0.2186,  0.1370,  
								0.1396,  0.0633,  0.0497,  0.0141,  0.0262,  0.1276,  
								0.2197,  0.3321,  0.2814,  0.3243,  0.2537,  0.2296,  
								0.0973,  0.0298,  0.0188,  0.0073,  0.0502,  0.2479,  
								0.2986,  0.5434,  0.4215,  0.3326,  0.1966,  0.1365,  
								0.0743,  0.0303,  0.0873,  0.2317,  0.3342,  0.3609,  
								0.4069,  0.3394,  0.1867,  0.1109,  0.0581,  0.0298,  
								0.0455,  0.1888,  0.4168,  0.5983,  0.5732,  0.4644,  
								0.3546,  0.2484,  0.1600,  0.0853,  0.0502,  0.1736,  
								0.4843,  0.7929,  0.7128,  0.7045,  0.4388,  0.3630,  
								0.1647,  0.0727,  0.0230,  0.1987,  0.7411,  0.9947,  
								0.9665,  0.8316,  0.5873,  0.2819,  0.1961,  0.1459,  
								0.0534,  0.0790,  0.2458,  0.4906,  0.5539,  0.5518,  
								0.5465,  0.3483,  0.3603,  0.1987,  0.1804,  0.0811,  
								0.0659,  0.1428,  0.4838,  0.8127 

							};

		double                  Mean;
		double                  TrainError;
		double                  TrainErrorPredictingMean;
		double                  TestError;
		double                  TestErrorPredictingMean;

		//FILE*                 f;


		void NormalizeSunspots()
		{
			int  Year;
			double Min, Max;
				
			Min = MAX_REAL;
			Max = MIN_REAL;
			for (Year=0; Year<NUM_YEARS; Year++) {
				Min = Math.Min(Min, Sunspots[Year]);
				Max = Math.Max(Max, Sunspots[Year]);
			}
			Mean = 0;
			for (Year=0; Year<NUM_YEARS; Year++) {
				Sunspots_[Year] = 
				Sunspots [Year] = ((Sunspots[Year]-Min) / (Max-Min)) * (HI-LO) + LO;
				Mean += Sunspots[Year] / NUM_YEARS;
			}
		}


		void InitializeApplication(ref NET Net)
		{
			int  Year, i;
			double Out, Err;

			Net.Alpha = 0.5;
			Net.Eta   = 0.05;
			Net.Gain  = 1;

			NormalizeSunspots();
			TrainErrorPredictingMean = 0;
			for (Year=TRAIN_LWB; Year<=TRAIN_UPB; Year++) 
			{
				for (i=0; i<M; i++) 
				{
					Out = Sunspots[Year+i];
					Err = Mean - Out;
					TrainErrorPredictingMean += 0.5 * (Err * Err);
				}
			}
			TestErrorPredictingMean = 0;
			for (Year=TEST_LWB; Year<=TEST_UPB; Year++) 
			{
				for (i=0; i<M; i++) 
				{
					Out = Sunspots[Year+i];
					Err = Mean - Out;
					TestErrorPredictingMean += 0.5 * (Err * Err);
				}
			}
			//f = fopen("BPN.txt", "w");
		}


		void FinalizeApplication(ref NET Net)
		{
			//fclose(f);
		}


		/******************************************************************************
								I N I T I A L I Z A T I O N
		******************************************************************************/


		void GenerateNetwork(ref NET Net)
		{
			int l,i;

			Net.Layer = new LAYER[NUM_LAYERS]; //(LAYER**) calloc(NUM_LAYERS, sizeof(LAYER*));
			   
			for (l=0; l<NUM_LAYERS; l++) 
			{
				Net.Layer[l] = new LAYER(); //(LAYER*) malloc(sizeof(LAYER));
			      
				Net.Layer[l].Units      = Units[l];
				Net.Layer[l].Output     = new double[Units[l] + 1]; //(double*)  calloc(Units[l]+1, sizeof(double));
				Net.Layer[l].Error      = new double[Units[l] + 1]; //(double*)  calloc(Units[l]+1, sizeof(double));
				Net.Layer[l].Weight     = new double[Units[l] + 1][]; //(double**) calloc(Units[l]+1, sizeof(double*));
				Net.Layer[l].WeightSave = new double[Units[l] + 1][]; //(double**) calloc(Units[l]+1, sizeof(double*));
				Net.Layer[l].dWeight    = new double[Units[l] + 1][]; //(double**) calloc(Units[l]+1, sizeof(double*));
				Net.Layer[l].Output[0]  = BIAS;
			      
				if (l != 0) 
				{
					for (i=1; i<=Units[l]; i++) 
					{
						Net.Layer[l].Weight[i]     = new double[Units[l-1] + 1]; //(double*) calloc(Units[l-1]+1, sizeof(double));
						Net.Layer[l].WeightSave[i] = new double[Units[l-1] + 1]; //(double*) calloc(Units[l-1]+1, sizeof(double));
						Net.Layer[l].dWeight[i]    = new double[Units[l-1] + 1]; //(double*) calloc(Units[l-1]+1, sizeof(double));
					}
				}
			}
			Net.InputLayer  = Net.Layer[0];
			Net.OutputLayer = Net.Layer[NUM_LAYERS - 1];
			Net.Alpha       = 0.9;
			Net.Eta         = 0.25;
			Net.Gain        = 1;
		}


		void RandomWeights(ref NET Net)
		{
			int l,i,j;
			   
			for (l=1; l<NUM_LAYERS; l++) 
			{
				for (i=1; i<=Net.Layer[l].Units; i++) 
				{
					for (j=0; j<=Net.Layer[l-1].Units; j++) 
					{
						Net.Layer[l].Weight[i][j] = RandomEqualREAL(-0.5, 0.5);
					}
				}
			}
		}


		//void SetInput(ref NET Net, double[] Input)
		void SetInput(ref NET Net, double[] Input, int Year)
		{
			int i;
			   
			for (i=1; i<=Net.InputLayer.Units; i++) 
			{
				//Net.InputLayer.Output[i] = Input[i-1];
				Net.InputLayer.Output[i] = Input[Year+i-1];
			}
		}


		void GetOutput(ref NET Net, double[] Output)
		{
			int i;
			   
			for (i=1; i<=Net.OutputLayer.Units; i++) 
			{
				Output[i-1] = Net.OutputLayer.Output[i];
			}
		}


		/******************************************************************************
					S U P P O R T   F O R   S T O P P E D   T R A I N I N G
		******************************************************************************/


		void SaveWeights(ref NET Net)
		{
			int l,i,j;

			for (l=1; l<NUM_LAYERS; l++) 
			{
				for (i=1; i<=Net.Layer[l].Units; i++) 
				{
					for (j=0; j<=Net.Layer[l-1].Units; j++) 
					{
						Net.Layer[l].WeightSave[i][j] = Net.Layer[l].Weight[i][j];
					}
				}
			}
		}


		void RestoreWeights(ref NET Net)
		{
			int l,i,j;

			for (l=1; l<NUM_LAYERS; l++) 
			{
				for (i=1; i<=Net.Layer[l].Units; i++) 
				{
					for (j=0; j<=Net.Layer[l-1].Units; j++) 
					{
						Net.Layer[l].Weight[i][j] = Net.Layer[l].WeightSave[i][j];
					}
				}
			}
		}


		/******************************************************************************
							P R O P A G A T I N G   S I G N A L S
		******************************************************************************/


		void PropagateLayer(ref NET Net, LAYER Lower, LAYER Upper)
		{
			int  i,j;
			double Sum;

			for (i=1; i<=Upper.Units; i++) 
			{
				Sum = 0;
				for (j=0; j<=Lower.Units; j++) 
				{
					Sum += Upper.Weight[i][j] * Lower.Output[j];
				}
				Upper.Output[i] = 1 / (1 + Math.Exp(-Net.Gain * Sum));
			}
		}


		void PropagateNet(ref NET Net)
		{
			int l;
			   
			for (l=0; l<NUM_LAYERS-1; l++) 
			{
				PropagateLayer(ref Net, Net.Layer[l], Net.Layer[l+1]);
			}
		}


		/******************************************************************************
						B A C K P R O P A G A T I N G   E R R O R S
		******************************************************************************/


		void ComputeOutputError(ref NET Net, double[] Target, int targetYear)
		{
			int  i;
			double Out, Err;
			   
			Net.Error = 0;
			for (i=1; i<=Net.OutputLayer.Units; i++) 
			{
				Out = Net.OutputLayer.Output[i];
				//Err = Target[i-1]-Out;
				Err = Target[targetYear+i-1]-Out;
				Net.OutputLayer.Error[i] = Net.Gain * Out * (1-Out) * Err;
				Net.Error += 0.5 * (Err * Err);
			}
		}


		void BackpropagateLayer(ref NET Net, LAYER Upper, LAYER Lower)
		{
			int  i,j;
			double Out, Err;
			   
			for (i=1; i<=Lower.Units; i++) 
			{
				Out = Lower.Output[i];
				Err = 0;
				for (j=1; j<=Upper.Units; j++) 
				{
					Err += Upper.Weight[j][i] * Upper.Error[j];
				}
				Lower.Error[i] = Net.Gain * Out * (1-Out) * Err;
			}
		}


		void BackpropagateNet(ref NET Net)
		{
			int l;
			   
			for (l=NUM_LAYERS-1; l>1; l--) 
			{
				BackpropagateLayer(ref Net, Net.Layer[l], Net.Layer[l-1]);
			}
		}


		void AdjustWeights(ref NET Net)
		{
			int  l,i,j;
			double Out, Err, dWeight;
			   
			for (l=1; l<NUM_LAYERS; l++) 
			{
				for (i=1; i<=Net.Layer[l].Units; i++) 
				{
					for (j=0; j<=Net.Layer[l-1].Units; j++) 
					{
						Out = Net.Layer[l-1].Output[j];
						Err = Net.Layer[l].Error[i];
						dWeight = Net.Layer[l].dWeight[i][j];
						Net.Layer[l].Weight[i][j] += Net.Eta * Err * Out + Net.Alpha * dWeight;
						Net.Layer[l].dWeight[i][j] = Net.Eta * Err * Out;
					}
				}
			}
		}


		/******************************************************************************
							S I M U L A T I N G   T H E   N E T
		******************************************************************************/


		//void SimulateNet(ref NET Net, double Input, double Output, double[] Target, bool Training)
		//void SimulateNet(ref NET Net, double[] Input, double[] Output, double[] Target, bool Training)
		void SimulateNet(ref NET Net, double[] Input, int inputYear, double[] Output, double[] Target, int targetYear, bool Training)
		{
			//SetInput(ref Net, Input);
			SetInput(ref Net, Input, inputYear);
			PropagateNet(ref Net);
			GetOutput(ref Net, Output);
			   
			//ComputeOutputError(ref Net, Target);
			ComputeOutputError(ref Net, Target, targetYear);
			if (Training) 
			{
				BackpropagateNet(ref Net);
				AdjustWeights(ref Net);
			}
		}


		void TrainNet(ref NET Net, int Epochs)
		{
			int  Year, n;
			double [] Output = new double[M];

			for (n=0; n<Epochs*TRAIN_YEARS; n++) 
			{
				Year = RandomEqualINT(TRAIN_LWB, TRAIN_UPB);
				//SimulateNet(ref Net, &(Sunspots[Year-N]), Output, &(Sunspots[Year]), true);
				//SimulateNet(ref Net, (Sunspots[Year-N]), Output, (Sunspots[Year]), true);
				SimulateNet(ref Net, Sunspots, Year-N, Output, Sunspots, Year, true);
			}
		}


		int iter = 0;
		void TestNet(ref NET Net)
		{
			int  Year;
			double [] Output = new double[M];

			TrainError = 0;
			for (Year=TRAIN_LWB; Year<=TRAIN_UPB; Year++) 
			{
				//SimulateNet(ref Net, &(Sunspots[Year-N]), Output, &(Sunspots[Year]), false);
				//SimulateNet(ref Net, (Sunspots[Year-N]), Output, (Sunspots[Year]), false);
				SimulateNet(ref Net, Sunspots, Year-N, Output, Sunspots, Year, false);
				TrainError += Net.Error;
			}
			TestError = 0;
			for (Year=TEST_LWB; Year<=TEST_UPB; Year++) 
			{
				//SimulateNet(ref Net, &(Sunspots[Year-N]), Output, &(Sunspots[Year]), false);
				//SimulateNet(ref Net, (Sunspots[Year-N]), Output, (Sunspots[Year]), false);
				SimulateNet(ref Net, Sunspots, Year-N, Output, Sunspots, Year, false);
				TestError += Net.Error;
			}
			iter = iter + 1;
			double out1 = TrainError / TrainErrorPredictingMean;
			double out2 = TestError / TestErrorPredictingMean;
			richTextBox.AppendText("\n" +iter.ToString() + " NMSE is " + out1.ToString("F3") + " on Training Set and " + out2.ToString("F3") + " on Test Set");
			//fprintf(f, "\nNMSE is %0.3f on Training Set and %0.3f on Test Set",
			//			TrainError / TrainErrorPredictingMean,
			//			TestError / TestErrorPredictingMean);
		}


		void EvaluateNet(ref NET Net)
		{
			int  Year;
			double [] Output = new double[M];
			double [] Output_ = new double[M];

			//fprintf(f, "\n\n\n");
			richTextBox.AppendText("\n\n\n");
			richTextBox.AppendText("Year    Sunspots    Open-Loop Prediction    Closed-Loop Prediction\n");
			richTextBox.AppendText("\n");
			//fprintf(f, "Year    Sunspots    Open-Loop Prediction    Closed-Loop Prediction\n");
			//fprintf(f, "\n");
			for (Year=EVAL_LWB; Year<=EVAL_UPB; Year++) 
			{
				//SimulateNet(ref Net, &(Sunspots [Year-N]), Output,  &(Sunspots [Year]), false);
				//SimulateNet(ref Net, (Sunspots [Year-N]), Output,  (Sunspots [Year]), false);
				SimulateNet(ref Net, Sunspots, Year-N, Output,  Sunspots, Year, false);
				//SimulateNet(ref Net, &(Sunspots_[Year-N]), Output_, &(Sunspots_[Year]), false);
				//SimulateNet(ref Net, (Sunspots_[Year-N]), Output_, (Sunspots_[Year]), false);
				SimulateNet(ref Net, Sunspots_, Year-N, Output_, Sunspots_, Year, false);
				Sunspots_[Year] = Output_[0];
				int out1 = FIRST_YEAR + Year;
				double out2 = Sunspots[Year];
				double out3 = Output[0];
				double out4 = Output_[0];
				richTextBox.AppendText(out1.ToString() + "\t\t" + out2.ToString("F3") + "\t\t" + out3.ToString("F3") + "\t\t" + out4.ToString("F3") + "\n");
				//fprintf(f, "%d       %0.3f                   %0.3f                     %0.3f\n",
				//		FIRST_YEAR + Year,
				//		Sunspots[Year],
				//		Output [0],
				//		Output_[0]);
			}
		}


		/******************************************************************************
											M A I N
		******************************************************************************/


		public void main()
		{
			NET  Net = new NET();
			bool Stop;
			double MinTestError;

			InitializeRandoms();
			GenerateNetwork(ref Net);
			RandomWeights(ref Net);
			InitializeApplication(ref Net);

			Stop = false;
			MinTestError = MAX_REAL;
			do 
			{
				TrainNet(ref Net, 10);
				TestNet(ref Net);
				if (TestError < MinTestError)
				{
					richTextBox.AppendText(" - saving Weights ...");
					//fprintf(f, " - saving Weights ...");
					MinTestError = TestError;
					SaveWeights(ref Net);
				}
				else if (TestError > 1.2 * MinTestError) 
				{
					richTextBox.AppendText(" - stopping Training and restoring Weights ...");
					//fprintf(f, " - stopping Training and restoring Weights ...");
					Stop = true;
					RestoreWeights(ref Net);
				}
			} while (!Stop);

			TestNet(ref Net);
			EvaluateNet(ref Net);
			   
			FinalizeApplication(ref Net);
		}
	}
}