stochastic.cpp

Simulated Annealing for both Stochastic models and Deterministic models

#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <cmath>
#include <fstream>
#include <stdio.h>
#include <time.h>
#define PI  3.1415926/180
#define PAI  3.1415926

using namespace std;


double ObjectFunction( double x, double y )
{

            double z = 0.0;

          //  z=sin(PI*x)*sin(PI*x)+cos(PI*y)*cos(PI*y); //F1 [-5,5]:   0.015195
          //z = sin(PI*x*y)+ x*x +y*y;//F2:   50.3556
         // z = sin(PI*4*x)*y+cos(PI*4*y)*x; //F3:   12.81206
		//	 z= (x*sin(x/PAI)+y*cos(y/PAI))*PAI*0.3/360; //F4 [0,82] range=0.758
    		z = PAI* x * 0.5 * cos(((y*y)/(PAI*PAI)+ 2*(x*x)/(PAI*PAI)))/(1+(y*y)/(PAI*PAI)+ 2*(x*x)/(PAI*PAI)) ;//F5  [0,45]
                                                               // range= 2.92846

            return z;
}	


double Stochastic ( double x, double y,  int n, double BF)
{
    double randf();
	double z;
	double sum = 0.0; 
	double stocha =0.0; //the stochastic data of Z'
	int j; 

            // z=sin(PI*x)*sin(PI*x)+cos(PI*y)*cos(PI*y); //F1 [-5,5]
           // z = sin(PI*x*y)+ x*x +y*y;//F2  50.3556
           // z = sin(PI*4*x)*y+cos(PI*4*y)*x; //F3
	      // z= (x*sin(x/PAI)+y*cos(y/PAI))*PAI*0.3/360; //F4 [0,82] range=0.758
	       z = PAI* x * 0.5 * cos(((y*y)/(PAI*PAI)+ 2*(x*x)/(PAI*PAI)))/(1+(y*y)/(PAI*PAI)+ 2*(x*x)/(PAI*PAI)) ;//F5  [0,45]
                                                               // [0,45]: range= 2.92846

	for(j=1;j<=n;j++)
	{
		stocha =( z - BF * 2.92846)+ (randf() * 2 * BF * 2.92846);
        sum = sum + stocha;
	}
		return sum / n;
}

	  
double randf() 
{ 
            return (double)(rand()/(double)RAND_MAX); 
}  



void main()

{       
	    int      InLoop = 0; 

		int      OutLoop = 0;

		double   TemperaturePause;

        int      MarkovLength;               // iteration

        double   Temperature;               // the initial temperature

        double   LastTemperature = 0.0;

	    double   Keep;


        double   DecayScale;     // decay parameter
			  
	    double   StepFactor;            //strategy for neighborhood

        double   BoundFactor = 0.0;           // stochastic, unif[Z-BoundFactor*Z,Z+BoundFactor*Z]
 
        int      sample = 0;                //  stochastic, # of samples of Z    

        const double XMAX = 45;          // define the space of X and Y

        const double YMAX = 45;  


        double   PreX,NextX;           // prior and next value of x 

        double   PreY,NextY;           // prior and next value of y 
			 
        double   PreBestX,PreBestY;    // the previous optimun value       

        double   BestX,BestY;          // the optimun value

        int      AcceptPoints = 0.0;      // #.of points during the Metropolis


        double   PreAveZ;
 
		double   NextAveZ; 
   
		double 	 BestAveZ;

		double   PreBestZ;


        double randf();
        
		double ObjectFunction( double x, double y );
	    
		double Stochastic ( double x, double y,  int n, double BF);
              

         // outside loop, change the temperature 

			 /*  cout<<"----------------------Stochastic---------------------"<<endl;
			   
			   cout<<"The Initial Temperature = "; cin>>Temperature; 
			   
			   cout<<"The Pause-Temperature = "; cin>>TemperaturePause; 
			   
			   cout<<"The Inside Loop = "; cin>>MarkovLength;  

               cout<<"The Decay Scale = "; cin>>DecayScale;  
			  
			    
                            
               cout<<"Neighborhood Step-Factor= "; cin>>StepFactor; 
			   
			   cout<<"Noise Half-Percentage= "; cin>>BoundFactor;
			   
	           cout<<"define the # of sample = "; cin>>sample;cout<<endl<<endl;"*/

			   Temperature = 2.78; 
			   TemperaturePause = 0.0424; 
               MarkovLength = 20;  
               DecayScale = 0.97;  
               StepFactor = 2.25;  //neighborhood
               BoundFactor = 0; // noise
               sample = 20;   
			   
			   ofstream SaveFile1("01.txt",ios_base::out);
			   
			   SaveFile1<<"Initial T ="<<Temperature<<";  Pause T="<<TemperaturePause<<";   Neighborhood Step-Factor ="<<StepFactor<<endl;
               
			   SaveFile1<<"Decay Scale ="<<DecayScale<<";  Noise Half-Percentage ="<<BoundFactor<<";  Sample Size ="<<sample<<";   Inside Loop ="<<MarkovLength<<endl<<endl; 

			   SaveFile1<<"      AveZ                               X                                        Y                                    Z"<<endl;  
								
			  		// pick the initial point randomly, initialization
              srand((unsigned int)time((time_t *)NULL));


              PreX = XMAX * randf();

              PreY = YMAX * randf();

			  PreAveZ = Stochastic (PreX,PreY,sample,BoundFactor);

              PreBestX = BestX = PreX;

              PreBestY = BestY = PreY;

              BestAveZ = PreAveZ;
			  
			  while ( Temperature > TemperaturePause)
			   {  // inside loop, at the current temperature

                   for (int i=0;i<MarkovLength;i++)

                   {
					   do

                       {

                            NextX = PreX - StepFactor + randf() * 2 * StepFactor; // TEST here for XMAX

                            NextY = PreY - StepFactor + randf() * 2 * StepFactor;

                       }
                        while ( !(NextX >= 0 && NextX <= XMAX && NextY >= 0 && NextY <= YMAX) ); //bound
                     
                        
                        //PreAveZ = Stochastic (PreX,PreY,sample,BoundFactor);

						NextAveZ = Stochastic (NextX,NextY,sample,BoundFactor); 

					//	BestAveZ = PreAveZ;
						
						
						if( PreAveZ- NextAveZ > 0 )

                       {
						    PreX=NextX;

                            PreY=NextY;

                            AcceptPoints++;

                            PreAveZ = NextAveZ; 

						//	SaveFile1<<PreAveZ<<"              "<<PreX<<"             "<<PreY<<"            "<<ObjectFunction (PreX,PreY)<<endl; 
					
                       }


                        else

						{  
							double change = -1 * ( NextAveZ - PreAveZ ) / Temperature ;

                            Keep = randf();

							if( exp(change) > Keep ) // exp > random value 0~1

                            {    
	
								 PreX=NextX;

                                 PreY=NextY;

                                 AcceptPoints++; //accept

	                             PreAveZ = NextAveZ;
	
							//	SaveFile1<<PreAveZ<<"              "<<PreX<<"             "<<PreY<<"            "<<ObjectFunction (PreX,PreY)<<endl; 

                           }
                      }
                  
                    	SaveFile1<<PreAveZ<<"                            "<<PreX<<"                           "<<PreY<<"                            "<<ObjectFunction (PreX,PreY)<<endl; 


						if (BestAveZ > NextAveZ )  // check the current best value, compared to the PreBest X,Y

					{   // keep the original BEST value

                            PreBestX = BestX;

                            PreBestY = BestY;

                            PreBestZ = BestAveZ;

                            // the new BEST value 

                            BestX = NextX;

                            BestY = NextY;

                            BestAveZ = NextAveZ;
                       }
				  
				  InLoop++;
				   
                   }
				  
				   LastTemperature = Temperature;
                   
				   Temperature = Temperature * DecayScale;
		
                   OutLoop++;
              
				 }
			      
			     SaveFile1<<"--------------------------"<<endl;
				 
				 SaveFile1<<"BestAveZ: "<<BestAveZ<<"     BestX: "<<BestX<<"      BestY: "<<BestY<<"       Z: "<<ObjectFunction(BestX,BestY)<<endl;
				  
				 SaveFile1<<"Accept Points: "<<AcceptPoints<<endl;  
				 
			     SaveFile1.close();
	
			      cout<<"-------------------------Final Result-----------------------"<<endl;
                  cout<<"The Best value of X is  "<<BestX<<endl;
				  cout<<"The Best value of Y is  "<<BestY<<endl;
				  cout<<"The Best average value of Z is  "<<BestAveZ<<endl;
				  cout<<"The Original Z is  "<<ObjectFunction(BestX,BestY)<<endl;
                  cout<<"The number of the accept points are  "<< AcceptPoints<<endl<<endl;
}