/******************************************************************************************
> [KAZOE]                                                                                 <
> Initial Sofware by [Yannick Verdie, Florent Lafarge], Copyright (C) Inria (GEOMETRICA)  <
> All Rights Reserved, 2012, [Version 1.0].                                               <
******************************************************************************************/
#include "Simulation_GPU.h"

	Simulation_GPU::Simulation_GPU(const CImg<int>  pictureInput , const ParameterFile parameterInput)
	{
		ptrPictureInput = new CImg<int>  (pictureInput);       
		ptrParamsInput = new ParameterFile(parameterInput);

		srand(time(0));

		m_time = 0;
		m_stat = NULL;
		m_ratio_birth_death = 0.5;
		m_method = METHOD_DATATERM::BAYESIAN;
	}

 
	Simulation_GPU::~Simulation_GPU()
	{
		if ( ptrParamsInput) delete ptrParamsInput;
		if (my_param) delete my_param;
		cuda_delete();
	}

 
	void Simulation_GPU::setStat(GPU_Stat* stat)
	{
		m_stat = stat;
	}

 
	const double Simulation_GPU 
		::getToTalEnergy() const
	{

		int number = 0;
		double energy = 0;

		Ellipse_Shape* result = NULL; 
		getResult_Ellipse_Shape(&result, &number, &energy);
		if (result) free(result);
		 
		return energy;

	}

	 	const double Simulation_GPU 
		::getToTalTime()
	{		
			 m_time += cuda_getTime();
			 cuda_resetTime();
			 return m_time;
	}

 
	void Simulation_GPU::geometricalDecrease (int verbose, bool saveLastState)
	{

		currentTemperature    =  this -> ptrParamsInput -> StartTemperature;
		currentDecreaseFactor =  this -> ptrParamsInput -> DecreaseFactor;
		
		double energy = 0; 
		double previousenergy = 0;
		int number = 0;
		int nb_iter = 0; 
		float mindeltaEnergy = this -> ptrParamsInput ->MinDeltaEnergy;

		while ( currentTemperature > this -> ptrParamsInput -> EndTemperature )
		{
			for ( int i = 0 ; i <  this -> ptrParamsInput -> PlateauLength; i++, nb_iter++)
			{
				tryCasualMove(currentTemperature);
			}
			currentTemperature *= currentDecreaseFactor;


			if (m_stat || verbose || mindeltaEnergy != 0)
			{
			   m_time += cuda_getTime();
			   //---------------------------------
			   {	//Ellipse
					Ellipse_Shape* a = NULL;
					getResult_Ellipse_Shape(&a, &number, &energy, false/*verbose*/);

				   	//do something with the result....
				   	if (saveLastState && currentTemperature <= this -> ptrParamsInput -> EndTemperature)
					{

		 				std::ofstream m_myfile("LastPosition.txt");
						m_myfile.setf ( ios::fixed, ios::floatfield );  

						for (int nb=0;nb<number;nb++)
							m_myfile<<a[nb].x<<" "<<a[nb].y<<" "<<a[nb].angle<<" "<<a[nb].axis0<<" "<<a[nb].axis1<<endl;

						m_myfile.close();
					}

					if (this -> ptrParamsInput -> SaveCurrentStageAsFile)
					{
						system("mkdir ResultImage");
						std::stringstream ss;
						ss<<"ResultImage/step"<<nb_iter<<".bmp";
						save(const_cast<char*>(ss.str().c_str()));  
					}

					//then delete
					if(a) free(a);
			   }

			  if (verbose)
				  printf("it:%d, ener:%f, time:%f, number:%d, T%f, dT%f\n", nb_iter,energy,m_time,number,currentTemperature,currentDecreaseFactor  );

			  if (m_stat)
				  m_stat->write(nb_iter,energy,m_time,number,currentTemperature,currentDecreaseFactor,0,0);

			  if (100*abs(previousenergy - energy) < abs(mindeltaEnergy*previousenergy) && previousenergy != 0)
			  {
				  printf("delta energy low, ... STOP (energy %f, previous%f, delta %f\n",energy,previousenergy,abs(previousenergy - energy)/abs(previousenergy));
				  break;
			  }

			  //------------------------------------
			  cuda_resetTime();

			  previousenergy  = energy;
			}
		}
	}


 
    void Simulation_GPU::startTimer()
	{	
		m_time = 0;
		cuda_startTimer();
	}

 
    float Simulation_GPU::getTime()
	{	
		return cuda_getTime();
	}

 
    void Simulation_GPU::waitUntilDone()
	{	
		cuda_blockandWait();
	}
	
	
	//
	//
	//init


    void Simulation_GPU::init(METHOD_DATATERM method)
	{	
		m_method = method;

		half_size_cell = ceil(max((float)this -> ptrParamsInput->Ellipse_MaxWidth,(float)this -> ptrParamsInput->Ellipse_MaxHeight));


		if ( this -> ptrParamsInput->ProbaTranslation != 0) half_size_cell+=this -> ptrParamsInput ->MaxTranslation;

		sizeBorder = half_size_cell;
		
		switch(m_method)
		{
		case METHOD_DATATERM::BAYESIAN:
		my_param = new bayesian_Parameter(
			this -> ptrParamsInput ->Ellipse_MinWidth, this -> ptrParamsInput ->Ellipse_MaxWidth, this -> ptrParamsInput ->Ellipse_MinHeight, this -> ptrParamsInput ->Ellipse_MaxHeight,
			 this -> ptrParamsInput ->Bayesian_Alpha, this -> ptrParamsInput ->Background_Avg, this -> ptrParamsInput ->Background_Sigma,
			this -> ptrParamsInput ->Shape_Avg, this -> ptrParamsInput ->Shape_Sigma, this -> ptrParamsInput ->Intersections_Alpha, this -> ptrParamsInput ->Hard_Core, this -> ptrParamsInput ->MaxTranslation,
			this-> ptrParamsInput->ConstructDataMap, const_cast<char*>(this-> ptrParamsInput->FileName.c_str()), this -> ptrParamsInput ->Verbose );
		break;

		case METHOD_DATATERM::BHATTACHARYA:
		my_param = new bhatt_Parameter(
			this -> ptrParamsInput ->Ellipse_MinWidth, this -> ptrParamsInput ->Ellipse_MaxWidth, this -> ptrParamsInput ->Ellipse_MinHeight, this -> ptrParamsInput ->Ellipse_MaxHeight,
			this -> ptrParamsInput ->Intersections_Alpha, this -> ptrParamsInput ->Hard_Core, this -> ptrParamsInput ->MaxTranslation,
			this -> ptrParamsInput ->Bhatt_Alpha, this -> ptrParamsInput ->Bhatt_Thickness, this -> ptrParamsInput ->Bhatt_Threshold, this-> ptrParamsInput->ConstructDataMap, const_cast<char*>(this-> ptrParamsInput->FileName.c_str()), this -> ptrParamsInput ->Verbose );
		break;
		}

		m_ratio_birth_death = cuda_RJMCMC_init(const_cast<int*>(ptrPictureInput->data()), ptrPictureInput->width(), ptrPictureInput->height(), half_size_cell, my_param, sizeBorder, method );
	}



	//
	//
	//
 
    double Simulation_GPU 
        ::tryCasualMove ( double temperature)
    {
        if ( this -> ptrParamsInput -> ProbaBirthDeath
             +  this -> ptrParamsInput -> ProbaEnlargeShrink 
             +  this -> ptrParamsInput -> ProbaTranslation
             +  this -> ptrParamsInput -> ProbaRotation
			 +  this -> ptrParamsInput -> ProbaBijection
			 +  this -> ptrParamsInput -> ProbaIncreateDecreaseH > 0 )
        {
            unsigned long randValue = rand() %
              ( this -> ptrParamsInput -> ProbaBirthDeath
                +  this -> ptrParamsInput -> ProbaEnlargeShrink
                +  this -> ptrParamsInput -> ProbaTranslation
                +  this -> ptrParamsInput -> ProbaRotation
				+  this -> ptrParamsInput -> ProbaBijection
				+  this -> ptrParamsInput -> ProbaIncreateDecreaseH
				+ 1);				 //between [0 prob]

			//randValue = 0;
            unsigned int sum = 0;

            if (randValue < (sum+= this->ptrParamsInput->ProbaBirthDeath))
            {
				return this -> tryBirthDeath ( temperature );
            }
            else if (randValue <(sum+=this->ptrParamsInput->ProbaEnlargeShrink))
            {
                return this -> tryEnlargeShrink ( temperature );
            }
            else if (randValue < (sum+= this->ptrParamsInput->ProbaTranslation))
            {
                return this -> tryTranslation ( temperature );
            }
            else if (randValue < (sum+= this->ptrParamsInput->ProbaRotation))
            {
                return this -> tryRotation ( temperature );
            }
		else if (randValue < (sum+= this->ptrParamsInput->ProbaBijection))
            {
                return this -> tryBijection ( temperature );
            }
		else if (randValue < (sum+= this -> ptrParamsInput -> ProbaIncreateDecreaseH))
			{
				return this -> tryChangeHeight ( temperature );
			}
        }
        return 0.0;
    }

 
	int Simulation_GPU::
	getNbShapes()
	{
		return getNbResult();
	}


 
	const void Simulation_GPU::
	save( const char* const ptrPathandName , char* nameImageBackground,bool verbose)
	{
		 printf("saving in %s...",   ptrPathandName);
		 CImg<int>  img(*this -> ptrPictureInput);
		 string n(nameImageBackground);
		 if	(n.empty())
		 {
			 printf("save from %s\n",nameImageBackground);
			img =CImg<int> (nameImageBackground);
		 }


		int number = 0;
		double energy = 0;

		{
			CImg<int> saver(img);


			double border = this -> ptrParamsInput ->Bhatt_Thickness;
			Ellipse_Shape* result = NULL;
			getResult_Ellipse_Shape(&result, &number, &energy);

			//loop for each
			for ( int   i = 0 ; i < number ; i++)
			{
				int searchArea = (int) ceil ( result[i].getRadius() + border + 0.5 );
				int searchAreaExt = searchArea + this->ptrParamsInput->Bhatt_Thickness;

				double axis0vec[2] = { cos(result[i].angle) , sin (result[i].angle) };
				double axis1vec[2] = { axis0vec[1] , -axis0vec[0] };

				for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
					for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
					{
						double dot0 = (axis0vec[0]*x)*(axis0vec[0]*x)+(axis0vec[1]*y)*(axis0vec[1]*y);
						double dot1 = (axis1vec[0]*x)*(axis1vec[0]*x)+(axis1vec[1]*y)*(axis1vec[1]*y);
						double ratio0big = dot0 / (result[i].axis0 + border + 1.5f);	//+1.5 instead of 0.5f  (draw a little bit larger)
						double ratio1big = dot1 / (result[i].axis1 + border + 1.5f);	//+1.5 instead of 0.5f  (draw a little bit larger)
						if ( ratio0big * ratio0big + ratio1big * ratio1big <= 1. )
						{
							double ratio0 = dot0 / (result[i].axis0 + 0.5 );
							double ratio1 = dot1 / (result[i].axis1 + 0.5 );
							if ( ratio0 * ratio0 + ratio1 * ratio1 > 1. )
							{
								double coordX = result[i].x + (double)x;
								if ( coordX < 0 )  coordX = floor ( coordX );
								else               coordX = ceil ( coordX );

								double coordY = result[i].y + (double)y;
								if ( coordY < 0 )  coordY = floor ( coordY );
								else               coordY = ceil ( coordY );

								if (coordX >=0 && coordX < this -> ptrPictureInput -> width() && coordY >=0 && coordY < this -> ptrPictureInput -> height() )
								{
									saver((int)coordX , (int)coordY ,0,0) = 255;
									saver((int)coordX , (int)coordY ,0,1) = 0;
									saver((int)coordX , (int)coordY ,0,2) = 0;
								}
							}
						}
					}


			}
			
			saver . save ( ptrPathandName );
			if (result) free(result);
			}

		 printf("done\n");
	}


 

	double Simulation_GPU::
        tryBirthDeath ( double temperature)
    {

        // ----------------------------
        // Birth
        // ----------------------------
        if ( rand()%101 < 100*m_ratio_birth_death) //rand()%101  is between [0 100]
        {

				birth_Simulation_GPU_Ellipse(temperature, const_cast<int*>(ptrPictureInput->data()), ptrPictureInput->width(), ptrPictureInput->height(),sizeBorder, my_param);
        }


        // ----------------------------
        // Death
        // ----------------------------
        else
        {
				death_Simulation_GPU_Ellipse(temperature, const_cast<int*>(ptrPictureInput->data()), ptrPictureInput->width(), ptrPictureInput->height(),sizeBorder, my_param);
        }

        return 0.0;
    }


	double Simulation_GPU::
        tryTranslation ( double temperature)
    {
		Translate_Simulation_GPU_Ellipse(temperature, const_cast<int*>(ptrPictureInput->data()), ptrPictureInput->width(), ptrPictureInput->height(),sizeBorder, my_param);
		return 0.0;
    }


	double Simulation_GPU::
        tryEnlargeShrink ( double temperature)
    {
		Enlarge_Simulation_GPU_Ellipse(temperature, const_cast<int*>(ptrPictureInput->data()), ptrPictureInput->width(), ptrPictureInput->height(),sizeBorder, my_param);
		return 0.0;
    }


	double Simulation_GPU::
        tryChangeHeight ( double temperature)
    {
        return 0.0;
    }

	double Simulation_GPU::
        tryRotation ( double temperature)
    {

		  Rotate_Simulation_GPU_Ellipse(temperature, const_cast<int*>(ptrPictureInput->data()), ptrPictureInput->width(), ptrPictureInput->height(),sizeBorder, my_param);
		return 0.0;
    }

	
	double Simulation_GPU::
        tryBijection ( double temperature)
    {
        return 0.0;
    }
