MLP_cannonball8

//+------------------------------------------------------------------+
//|                                               MLP_cannonball.mq4 |
//|                      Copyright © 2008, MetaQuotes Software Corp. |
//|                                        http://www.metaquotes.net |
//+------------------------------------------------------------------+
#property copyright "Copyright © 2008, MetaQuotes Software Corp."
#property link      "http://www.metaquotes.net"
/* This program implements a simple Multi Layer Perceptron to solve the canonball problem (regression)
   Adapted by Sebastien Marcel (2003-2004) from D. Collobert
   The goal is to estimate the position X of a bullet at time t (fixed) given the initial speed (v) and angle (a)
	  
   y
   
  / \
   |
   |
   |
   |               . X
   |            .
   |          .
   |        .
   |      .
   |     .
   |    .
   |   .
   |  .
   | .
   |.   (v,a)
   --------------------------> x
   
*/

#define PI 3.14159265358979323846
#define RAND_MAX 32767.0

//
// Datasets
#define N_PATTERNS_TRAIN 500
#define N_PATTERNS_TEST 500
#define N_PATTERNS_TRAIN_TEST 1000 //N_PATTERNS_TRAIN + N_PATTERNS_TEST

//  w1[N_HU1][IN1]=>w1[4][3]
//     H  In
//     |  |
//  w1[0][0] w1[0][1] w1[0][2]  
//  w1[1][0] w1[1][1] w1[1][2]          
//  w1[2][0] w1[2][1] w1[2][2]          
//  w1[3][0] w1[3][1] w1[3][2]
//  forward train:
//  aHidden[0] = xBias;
//  aHidden[1] = xBias*w1[1][0] + (In0)*w1[1][1] + (In1)*w1[1][2];
//  aHidden[2] = xBias*w1[2][0] + (In0)*w1[2][1] + (In1)*w1[2][2];
//  aHidden[3] = xBias*w1[3][0] + (In0)*w1[3][1] + (In1)*w1[3][2];
//  yHidden[0] = xBias;
//  yHidden[1] = f(aHidden[1]);
//  yHidden[2] = f(aHidden[2]);
//  yHidden[3] = f(aHidden[3]);
//
//  w2[OUT][N_HU1]=>w2[2][4]
//
//  w2[0][0] w2[0][1] w2[0][2] w2[0][3]
//  w2[1][0] w2[1][1] w2[1][2] w2[1][3]
//     |  |
//    Out H
// forward train:
//	aOutput[0] = 0.0;
//	aOutput[1] = 0.0;
//	aOutput[0] = aOutput[0] + yHidden[0]*w2[0][0];
//	aOutput[1] = aOutput[1] + yHidden[0]*w2[1][0];
//	aOutput[0] = aOutput[0] + yHidden[1]*w2[0][1];
//	aOutput[1] = aOutput[1] + yHidden[1]*w2[1][1];
//	aOutput[0] = aOutput[0] + yHidden[2]*w2[0][2];
//	aOutput[1] = aOutput[1] + yHidden[2]*w2[1][2];
//	aOutput[0] = aOutput[0] + yHidden[3]*w2[0][3];
//	aOutput[1] = aOutput[1] + yHidden[3]*w2[1][3];
//

//	out0 = f(aOutput[0]);
//	out1 = f(aOutput[1]);

#define IN 18 //2*9
#define IN1 19 //IN+1
#define OUT 2
#define N_HU 27 //3*9
#define N_HU1 28 //N_HU + 1
//*****************************
// Inputs
//
	double DataIn0[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn1[N_PATTERNS_TRAIN_TEST]; // angle of bullet
	double DataIn2[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn3[N_PATTERNS_TRAIN_TEST]; // angle of bullet	
	double DataIn4[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn5[N_PATTERNS_TRAIN_TEST]; // angle of bullet	
	double DataIn6[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn7[N_PATTERNS_TRAIN_TEST]; // angle of bullet	

	double DataIn8[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn9[N_PATTERNS_TRAIN_TEST]; // angle of bullet
	double DataIn10[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn11[N_PATTERNS_TRAIN_TEST]; // angle of bullet	
	double DataIn12[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn13[N_PATTERNS_TRAIN_TEST]; // angle of bullet	
	double DataIn14[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn15[N_PATTERNS_TRAIN_TEST]; // angle of bullet	
	double DataIn16[N_PATTERNS_TRAIN_TEST]; // speed of bullet
	double DataIn17[N_PATTERNS_TRAIN_TEST]; // angle of bullet	
			
// X[N_PATTERNS_TRAIN + N_PATTERNS_TEST];
//*****************************
// Targets
//
	double DataOut0[N_PATTERNS_TRAIN_TEST]; // x position at time t
	double DataOut1[N_PATTERNS_TRAIN_TEST]; // y position at time t
// Y[N_PATTERNS_TRAIN + N_PATTERNS_TEST];

//*****************************
// to create the data
//
double vmax;
double gravity;

//*****************************
// to train the MLP
//
double lambda =  0.01;     // learning rate
double mu = 0.6;           // inertia momentum rate
double mse_min = 0.001;    // minimum Mean Squared Error
int max_iterations = 500; // maximum number of iterations

//*****************************
// MLP data


// weights between hidden neurons and inputs
double w1[N_HU1][IN1];
double w1old[N_HU1][IN1];

// weights between outputs and hidden neurons
double w2[OUT][N_HU1];
double w2old[OUT][N_HU1];

// values of integration function
double aHidden[N_HU1];
double aOutput[OUT];

// values of transfert function
double yHidden[N_HU1];

// constant value of bias
double xBias = 1.0;


int init()
{
  vmax = 1.0/10;
  gravity = 1.0/200;
  
  return(0);
}
//+------------------------------------------------------------------+
//| script program start function                                    |
//+------------------------------------------------------------------+
int start()
  {
//----
   initrand();
   main();
//----
   return(0);
  }
//+------------------------------------------------------------------+

//*****************************
// Sigmoid transfer function
//
double f(double x)
{
	return (1.0 /(1.0 + MathExp(-x)));
}

//*****************************
// Derivative of Sigmoid
//
double f_prime(double x)
{
   double z = MathExp(-x);
	double one_plus_z = 1.0 + z;
	
	return (z / (one_plus_z * one_plus_z));
}



//*****************************
// Compute the x position of bullet at time t
//
double xpos(double t, double v, double teta)
{
	return( v * t * MathCos(teta*PI/180.0));
}

//*****************************
// Compute the y position of bullet at time t
//
double ypos(double t, double v, double teta)
{
	return( -0.5 * gravity * t * t + v * t * MathSin(teta*PI/180.0));
}

//*****************************
// Compute the MSE
//
double MSE(
            double a0, 
            double a1, 
            double b0,
            double b1
          )
{
	return (0.5*((a0 - b0)*(a0 - b0) + (a1 - b1)*(a1 - b1)));
}

//*****************************
// Create the MLP  
//
void createMLP()
{
	int i, j;
	// weights between hidden neurons and inputs
	for (i = 0; i<IN1; i++)
	{
		for (j = 0; j< N_HU1; j++)
		{	 
			w1[j][i] = Random(-1.0, 1.0);
			w1old[j][i] = w1[j][i];
		}
	}
   // weights between outputs and hidden neurons
	for (i = 0; i<OUT; i++)
	{
		for (j = 0; j< N_HU1; j++)
		{	 
			w2[i][j] = Random(-1.0, 1.0);
			w2old[i][j] = w2[i][j];
		}
	}
}

//*****************************
// Create the datasets
//
void createDatasets()
{
	int i;

	i = 0;
	while (i < (N_PATTERNS_TRAIN + N_PATTERNS_TEST))
	{  
		DataIn0[i]  = (Random(0.0, vmax)) / vmax;
		DataIn1[i]  = (Random(0.0, 90.0)) / 90.0;
		DataIn2[i]  = (Random(0.0, vmax)) / vmax;
		DataIn3[i]  = (Random(0.0, 90.0)) / 90.0;
		DataIn4[i]  = (Random(0.0, vmax)) / vmax;
		DataIn5[i]  = (Random(0.0, 90.0)) / 90.0;		
		DataIn6[i]  = (Random(0.0, vmax)) / vmax;
		DataIn7[i]  = (Random(0.0, 90.0)) / 90.0;	

		DataIn8[i]  = (Random(0.0, vmax)) / vmax;
		DataIn9[i]  = (Random(0.0, 90.0)) / 90.0;
		DataIn10[i] = (Random(0.0, vmax)) / vmax;
		DataIn11[i] = (Random(0.0, 90.0)) / 90.0;
		DataIn12[i] = (Random(0.0, vmax)) / vmax;
		DataIn13[i] = (Random(0.0, 90.0)) / 90.0;		
		DataIn14[i] = (Random(0.0, vmax)) / vmax;
		DataIn15[i] = (Random(0.0, 90.0)) / 90.0;	
		DataIn16[i] = (Random(0.0, vmax)) / vmax;
		DataIn17[i] = (Random(0.0, 90.0)) / 90.0;	
				
		DataOut0[i] = xpos(10, DataIn0[i] * vmax, (DataIn1[i])*90.0);
		DataOut1[i] = ypos(10, DataIn0[i] * vmax, (DataIn1[i])*90.0);

		if (DataOut1[i] > 0.0) i = i + 1;
   }
}

//*****************************
// Forward the input in the MLP
//
void forward(
             double In0,
             double In1,
             double In2,
             double In3,
             double In4,
             double In5, 
             double In6,
             double In7,                                        

             double In8,
             double In9,
             double In10,
             double In11,
             double In12,
             double In13, 
             double In14,
             double In15,                                        
             double In16,
             double In17,                                        

             double &out0,
             double &out1
             )
{
	int i;

	aHidden[0] = xBias;
	yHidden[0] = xBias;
	
	for (i = 1; i <= N_HU; i++)
	{
		aHidden[i] = xBias*w1[i][0] + 
		                              (In0)*w1[i][1] + 
		                              (In1)*w1[i][2] + 
		                              (In2)*w1[i][3] + 
		                              (In3)*w1[i][4] +
		                              (In4)*w1[i][5] + 
		                              (In5)*w1[i][6] +
		                              (In6)*w1[i][7] + 
		                              (In7)*w1[i][8] +

		                              (In8)*w1[i][9] + 
		                              (In9)*w1[i][10] + 
		                              (In10)*w1[i][11] + 
		                              (In11)*w1[i][12] +
		                              (In12)*w1[i][13] + 
		                              (In13)*w1[i][14] +
		                              (In14)*w1[i][15] + 
		                              (In15)*w1[i][16] +
		                              (In16)*w1[i][17] + 
		                              (In17)*w1[i][18] ;
		                              
		yHidden[i] = f(aHidden[i]);
	}

	aOutput[0] = 0.0;
	aOutput[1] = 0.0;

	for (i = 0; i <= N_HU; i++)
	{
		aOutput[0] = aOutput[0] + yHidden[i]*w2[0][i];
		aOutput[1] = aOutput[1] + yHidden[i]*w2[1][i];
	}

	out0 = f(aOutput[0]);
	out1 = f(aOutput[1]);
} 

//*****************************
// Backward (back-propagate the gradient of error)
//
void backward(
               double Yestimate0, 
               double Yestimate1, 
               double Outwant0,
               double Outwant1,
               double In0, 
               double In1, 
               double In2, 
               double In3,  
               double In4, 
               double In5,   
               double In6, 
               double In7,  

               double In8, 
               double In9, 
               double In10, 
               double In11,  
               double In12, 
               double In13,   
               double In14, 
               double In15,     
               double In16, 
               double In17,     
                                                                           
               double lambda, 
               double mu
               )
{
	double wnew, Goutput[OUT], Ghidden[N_HU1];
	int i,j;

	// 
	Goutput[0] = (Yestimate0 - Outwant0) * f_prime(aOutput[0]);
	Goutput[1] = (Yestimate1 - Outwant1) * f_prime(aOutput[1]); 

	//
	for (i = 0; i<=N_HU; i++)
		Ghidden[i] = (Goutput[0]*w2[0][i] + Goutput[1]*w2[1][i])*f_prime(aHidden[i]);

	//
	for (j = 0; j<OUT;j++)
	{
		for (i = 0; i<=N_HU; i++)
		{
			wnew = w2[j][i] - lambda * yHidden[i] * Goutput[j] + mu * (w2[j][i] - w2old[j][i]);
			w2old[j][i] = w2[j][i];
			w2[j][i] = wnew;
		}
	}

	//
	for (i = 0; i<=N_HU; i++)
	{
		wnew = w1[i][0] - lambda * xBias * Ghidden[i] + mu * (w1[i][0] - w1old[i][0]);
		w1old[i][0] = w1[i][0];
		w1[i][0] = wnew;
	   
		wnew = w1[i][1] - lambda * In0 * Ghidden[i] + mu * (w1[i][1] - w1old[i][1]);
		w1old[i][1] = w1[i][1];
		w1[i][1] = wnew;
	   
		wnew = w1[i][2] - lambda * In1 * Ghidden[i] + mu * (w1[i][2] - w1old[i][2]);
		w1old[i][2] = w1[i][2];
		w1[i][2] = wnew;

		wnew = w1[i][3] - lambda * In2 * Ghidden[i] + mu * (w1[i][3] - w1old[i][3]);
		w1old[i][3] = w1[i][3];
		w1[i][3] = wnew;
	   
		wnew = w1[i][4] - lambda * In3 * Ghidden[i] + mu * (w1[i][4] - w1old[i][4]);
		w1old[i][4] = w1[i][4];
		w1[i][4] = wnew;		

		wnew = w1[i][5] - lambda * In4 * Ghidden[i] + mu * (w1[i][5] - w1old[i][5]);
		w1old[i][5] = w1[i][5];
		w1[i][5] = wnew;
	   
		wnew = w1[i][6] - lambda * In5 * Ghidden[i] + mu * (w1[i][6] - w1old[i][6]);
		w1old[i][6] = w1[i][6];
		w1[i][6] = wnew;		

		wnew = w1[i][7] - lambda * In6 * Ghidden[i] + mu * (w1[i][7] - w1old[i][7]);
		w1old[i][7] = w1[i][7];
		w1[i][7] = wnew;
	   
		wnew = w1[i][8] - lambda * In7 * Ghidden[i] + mu * (w1[i][8] - w1old[i][8]);
		w1old[i][8] = w1[i][8];
		w1[i][8] = wnew;	

		wnew = w1[i][9] - lambda * In8 * Ghidden[i] + mu * (w1[i][9] - w1old[i][9]);
		w1old[i][9] = w1[i][9];
		w1[i][9] = wnew;				

		wnew = w1[i][10] - lambda * In9 * Ghidden[i] + mu * (w1[i][10] - w1old[i][10]);
		w1old[i][10] = w1[i][10];
		w1[i][10] = wnew;		

		wnew = w1[i][11] - lambda * In10 * Ghidden[i] + mu * (w1[i][11] - w1old[i][11]);
		w1old[i][11] = w1[i][11];
		w1[i][11] = wnew;		
		
		wnew = w1[i][12] - lambda * In11 * Ghidden[i] + mu * (w1[i][12] - w1old[i][12]);
		w1old[i][12] = w1[i][12];
		w1[i][12] = wnew;		

		wnew = w1[i][13] - lambda * In12 * Ghidden[i] + mu * (w1[i][13] - w1old[i][13]);
		w1old[i][13] = w1[i][13];
		w1[i][13] = wnew;		

		wnew = w1[i][14] - lambda * In13 * Ghidden[i] + mu * (w1[i][14] - w1old[i][14]);
		w1old[i][14] = w1[i][14];
		w1[i][14] = wnew;		

		wnew = w1[i][15] - lambda * In14 * Ghidden[i] + mu * (w1[i][15] - w1old[i][15]);
		w1old[i][15] = w1[i][15];
		w1[i][15] = wnew;		

		wnew = w1[i][16] - lambda * In15 * Ghidden[i] + mu * (w1[i][16] - w1old[i][16]);
		w1old[i][16] = w1[i][16];
		w1[i][16] = wnew;		

		wnew = w1[i][17] - lambda * In16 * Ghidden[i] + mu * (w1[i][17] - w1old[i][17]);
		w1old[i][17] = w1[i][17];
		w1[i][17] = wnew;		

		wnew = w1[i][18] - lambda * In17 * Ghidden[i] + mu * (w1[i][18] - w1old[i][18]);
		w1old[i][18] = w1[i][18];
		w1[i][18] = wnew;		
	}
}

//*****************************
// Stochastic gradient training
//
double train(int P0,int P)
{
   double Yestimate0;
   double Yestimate1; 
	double mse_;
	double mse_total;

	mse_total = 0.0;
   
	// For each train patterns
	for(int p = P0 ; p < P ; p++)
	{
		// Forward current train pattern into the MLP
		forward(
		          //X[p], 
               DataIn0[p],
               DataIn1[p],
               DataIn2[p],
               DataIn3[p],  
               DataIn4[p],
               DataIn5[p],        
               DataIn6[p],
               DataIn7[p],  
               DataIn8[p],
               DataIn9[p],
               DataIn10[p],
               DataIn11[p],  
               DataIn12[p],
               DataIn13[p],        
               DataIn14[p],
               DataIn15[p],  
               DataIn16[p],
               DataIn17[p],  
                                                                     
		          //&Yestimate
               Yestimate0,
               Yestimate1		          
		        );
      
		// Computes the MSE
		mse_ = MSE(
		             //Y[p], 
            DataOut0[p], 
            DataOut1[p], 
		             //Yestimate
            Yestimate0,
            Yestimate1
		          );
      
		// Accumulate the MSE
		mse_total = mse_total + mse_;

		// Backward MSE gradient
		backward(
		          //Yestimate, 
                Yestimate0, 
                Yestimate1, 
		          //Y[p], 
                DataOut0[p], //Outwant
                DataOut1[p], //Outwant
		          //X[p], 
                DataIn0[p], 
                DataIn1[p], 
                DataIn2[p], 
                DataIn3[p],   
                DataIn4[p], 
                DataIn5[p],    
                DataIn6[p], 
                DataIn7[p], 
                   
                DataIn8[p], 
                DataIn9[p], 
                DataIn10[p], 
                DataIn11[p],   
                DataIn12[p], 
                DataIn13[p],    
                DataIn14[p], 
                DataIn15[p],                                                             	          
                DataIn16[p], 
                DataIn17[p],                                                             	          
                
		          lambda, 
		          mu
		          ); 
	}

	// Return normalized train MSE
	double P1=P-P0;
	return (mse_total /P1);
}
/*
		//*****************************
		// test
		//*****************************		
		mse_test = 0.0;
		// For each test pattern
		// (N_PATTERNS_TRAIN + N_PATTERNS_TEST)                       N_PATTERNS_TRAIN                            0
		// +-----------------------------------------------------------+------------------------------------------+
		//
		for(int i = N_PATTERNS_TRAIN ; i < (N_PATTERNS_TRAIN + N_PATTERNS_TEST) ; i++)
		{
			Yestimate0 = 0.0;
			Yestimate1 = 0.0;
			// forward current pattern into the MLP
			//forward(X[i], &Yestimate);
		    forward(
		          //X[p], 
               DataIn0[i],
               DataIn1[i],
               DataIn2[i],
               DataIn3[i],               
               DataIn4[i],
               DataIn5[i], 
               DataIn6[i],
               DataIn7[i],  

               DataIn8[i],
               DataIn9[i],
               DataIn10[i],
               DataIn11[i],               
               DataIn12[i],
               DataIn13[i], 
               DataIn14[i],
               DataIn15[i],  
               DataIn16[i],
               DataIn17[i],  
                                            
		          //&Yestimate
               Yestimate0,
               Yestimate1		          
		        );			
			// computes MSE
			//mse_ = MSE(Y[i], Yestimate);
			mse_ = MSE(
		             //Y[p], 
            DataOut0[i], 
            DataOut1[i], 
		             //Yestimate
            Yestimate0,
            Yestimate1
		          );
			// accumulate MSE
			mse_test += mse_;
  		}       
		// Normalize the MSE
		mse_test /= N_PATTERNS_TEST;
*/
double test(int P0,int P)
{
   double Yestimate0;
   double Yestimate1; 
	double mse_;
	double mse_total;

	mse_total = 0.0;
   
	// For each train patterns
	for(int p = P0 ; p < P ; p++)
	{
			Yestimate0 = 0.0;
			Yestimate1 = 0.0;	
		// Forward current train pattern into the MLP
		forward(
		          //X[p], 
               DataIn0[p],
               DataIn1[p],
               DataIn2[p],
               DataIn3[p],  
               DataIn4[p],
               DataIn5[p],        
               DataIn6[p],
               DataIn7[p],  
               DataIn8[p],
               DataIn9[p],
               DataIn10[p],
               DataIn11[p],  
               DataIn12[p],
               DataIn13[p],        
               DataIn14[p],
               DataIn15[p],  
               DataIn16[p],
               DataIn17[p],  
                                                                     
		          //&Yestimate
               Yestimate0,
               Yestimate1		          
		        );
      
		// Computes the MSE
		mse_ = MSE(
		             //Y[p], 
            DataOut0[p], 
            DataOut1[p], 
		             //Yestimate
            Yestimate0,
            Yestimate1
		          );
      
		// Accumulate the MSE
		mse_total = mse_total + mse_;

		// Backward MSE gradient
		backward(
		          //Yestimate, 
                Yestimate0, 
                Yestimate1, 
		          //Y[p], 
                DataOut0[p], //Outwant
                DataOut1[p], //Outwant
		          //X[p], 
                DataIn0[p], 
                DataIn1[p], 
                DataIn2[p], 
                DataIn3[p],   
                DataIn4[p], 
                DataIn5[p],    
                DataIn6[p], 
                DataIn7[p], 
                   
                DataIn8[p], 
                DataIn9[p], 
                DataIn10[p], 
                DataIn11[p],   
                DataIn12[p], 
                DataIn13[p],    
                DataIn14[p], 
                DataIn15[p],                                                             	          
                DataIn16[p], 
                DataIn17[p],                                                             	          
                
		          lambda, 
		          mu
		          ); 
	}

	// Return normalized train MSE
	double P1=P-P0;
	return (mse_total /P1);
}


//
int main()
{
	double mse_train;
	double mse_test;
   double mse_;
   double Yestimate0; 
   double Yestimate1; 
	//*****************************
	createMLP();
	//*****************************
	createDatasets();
	//*****************************	
	//for (int i = 0; i< N_PATTERNS_TRAIN; i++) printf(" TRN: x=[%f %f] y=[%f %f]\n", X[i].v, X[i].a, Y[i].x, Y[i].y);
	//for (int i = N_PATTERNS_TRAIN; i<(N_PATTERNS_TRAIN + N_PATTERNS_TEST); i++) printf(" TST: x=[%f %f] y=[%f %f]\n", X[i].v, X[i].a, Y[i].x, Y[i].y);
	Print("Stochastic gradient training:");
	//
	int pf_train = FileOpen("mse_train_mt4_i18.txt",FILE_CSV|FILE_WRITE,' ');//FileOpen("mse_train.txt", "w");
	int pf_test = FileOpen("mse_test_mt4_i18.txt",FILE_CSV|FILE_WRITE,' ');//FileOpen("mse_test.txt", "w");

	int iter;
	
	for(iter = 1 ; iter <= max_iterations ; iter++)
	{
		//*****************************
		// train
		//*****************************		
		mse_train = train(0,N_PATTERNS_TRAIN);
		FileWrite(pf_train,mse_train);
		//*****************************
		// test
		//*****************************		
      mse_test = test(N_PATTERNS_TRAIN,(N_PATTERNS_TRAIN + N_PATTERNS_TEST));
		FileWrite(pf_test,mse_test);
		Print("->",iter); 
		//fflush(stdout);
		if (mse_train < mse_min) break;
   }
	//Print("\n");
	//
	FileClose(pf_test);
	FileClose(pf_train);
	//*****************************
	// Print info about training
	//
	Print("Number of iterations = ", iter);
	Print("Final MSE train = ", mse_train);
	Print("Final MSE test = ", mse_test);
	Print("End of program.");

	return (0);
}

//use this first function to seed the random number generator,
//call this before any of the other functions
void initrand()
{
    //srand((unsigned)(time(0)));
    MathSrand(TimeLocal());
}
//generates a psuedo-random double between 0.0 and 0.999...
double randdouble0()
{
    return (MathRand()/((RAND_MAX)+1));
}
//generates a psuedo-random double between min and max
double randdouble2(double min, double max)
{
    if (min>max)
    {
        return (randdouble0()*(min-max)+max);    
    }
    else
    {
        return (randdouble0()*(max-min)+min);
    }
}


// Random generator
double Random(double inf, double sup)
{
        return(randdouble2(inf,sup));
}