/**
She works!!!

Change POINT_NUM, DEG,CENTER_POINT


Press ctrl-c to shrink the window.

***/



#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <obstack.h>
#include <gmp.h>
#include <signal.h>
//#include <iomanip.h>



#define TRYS 1000
#define MAX_DEPTH 100


#define DEG 5 /*Actualy degree-1*/
#define POINT_NUM 2
#define SAMPLE_SIZE 50
#define ITERATIONS 100000

#define SPECTRAL 1
#define CENTER_POINT 0


/*The Precision*/
#define PREC 512



//The co-efficents for the polynomials
mpf_t RAND_COEFS[SAMPLE_SIZE][DEG];
//The evaluation of their integrals
mpf_t RAND_AVES[SAMPLE_SIZE];
//fittness
mpf_t FIT;
mpf_t BEST_FIT;
//New element
mpf_t NEW_ELT[POINT_NUM];
mpf_t BEST[POINT_NUM];

mpf_t VAL,TEMP1,TEMP2,TEMP3;

mpf_t WINDOW_SIZE;


void cnt(int sig)
{  
  mpf_div_ui(WINDOW_SIZE,WINDOW_SIZE,2);
  printf("\n");
  mpf_out_str(NULL,10,10,WINDOW_SIZE);
  printf("\n");
}


//Get the fittest member of SWARM_POP, swap him to SWARM_POP[0]
int findFittest()
{
  int i;

  if(mpf_cmp(FIT,BEST_FIT)<0)
    {
      mpf_swap(FIT,BEST_FIT);
      for(i=0;i<POINT_NUM;i++)
	mpf_swap(NEW_ELT[i],BEST[i]);
      return 1;
    }
  return 0;
}


//integrate the randomly generated polynomials from 0 to 1.
//store their value in RAND_AVES
//F(1)-F(0)
//
// ax^2+bx+c
// (a/3) * x^3 +(b/2)x^2 + c*x 
//
void integrate()
{
  int i,j;
  for(i=0;i<SAMPLE_SIZE;i++)
    {
      mpf_set_d(RAND_AVES[i],0.0);
      //RAND_AVES[i]=0.0;
      for(j=0;j<DEG;j++)
	{
	  mpf_set_d(TEMP2,(double)(j+1.0));
	  mpf_div(TEMP1,RAND_COEFS[i][j],TEMP2);
	  mpf_add(RAND_AVES[i],RAND_AVES[i],TEMP1);
	  //RAND_AVES[i]+=(RAND_COEFS[i][j]/((double)j+1.0));
	}
    }
}

void polyGen()
{
  int i,j;
  for(i=0;i<SAMPLE_SIZE;i++)
    {
      for(j=0;j<DEG;j++)
	mpf_set_d(RAND_COEFS[i][j], (((double)rand())/ ((double)RAND_MAX)*2.0) -1.0);
    }
}
void evalPointGen()
{
  int i,j;


  mpf_set_d(TEMP2,0.5);
  mpf_set_d(TEMP3,1.0);

  for(j=0;j<POINT_NUM;j++)
    {
      //mpf_set_d(TEMP1,(double)rand());
      //mpf_set_d(TEMP2,(double)RAND_MAX);
      //mpf_div(SWARM_POP[i][j],TEMP1,TEMP2);  
      mpf_set_d(NEW_ELT[j],(double) ((double)rand()/ (double)RAND_MAX)); 
      
      mpf_sub(TEMP1, NEW_ELT[j],TEMP2);
      mpf_mul(TEMP1,TEMP1,WINDOW_SIZE);
      mpf_add(NEW_ELT[j],BEST[j],TEMP1);
      if(mpf_sgn(NEW_ELT[j])<0)
	mpf_set_d(NEW_ELT[j],0.000000001);
      if(mpf_cmp(NEW_ELT[j],TEMP3)>0)
	mpf_set_d(NEW_ELT[j],.9999999);
    }
}

//Evaluate each member of the population to see how close they were to the
//actual answers


void evalSingle()
{
  int i,j,l;
  int realPointNum;
  mpf_set_d(FIT,0.0);
  for(i=0;i<SAMPLE_SIZE;i++)
    {
      mpf_set_d(VAL,0.0);
      for(j=0;j<DEG;j++)
	{
	  for(l=0;l<POINT_NUM;l++)
	    {
	      mpf_pow_ui(TEMP1,NEW_ELT[l],(long)j);
	      mpf_mul(TEMP1, TEMP1,RAND_COEFS[i][j]);
	      mpf_add(VAL,VAL,TEMP1);
	      if(SPECTRAL)
		{
		  mpf_set_d(TEMP2,1.0);
		  mpf_sub(TEMP2,TEMP2,NEW_ELT[l]);
		  mpf_pow_ui(TEMP1,TEMP2,(long)j);
		  mpf_mul(TEMP1, TEMP1,RAND_COEFS[i][j]);
		  mpf_add(VAL,VAL,TEMP1);
		}
	    }
	  if(CENTER_POINT)
	    {
	      mpf_set_ui(TEMP2,.5);
	      mpf_pow_ui(TEMP1,TEMP2,(long)j);
	      mpf_mul(TEMP1, TEMP1,RAND_COEFS[i][j]);
	      mpf_add(VAL,VAL,TEMP1);
	    }
	}
      realPointNum=CENTER_POINT+POINT_NUM+(POINT_NUM*SPECTRAL);
      mpf_set_d(TEMP1,(double)realPointNum);
      mpf_div(TEMP1,VAL,TEMP1);
      mpf_sub(TEMP1,TEMP1,RAND_AVES[i]);
      //VAL=(VAL/(double)POINT_NUM)-RAND_AVES[i];
      mpf_abs(TEMP1,TEMP1);
      mpf_add(FIT,FIT,TEMP1);
    }
}





int main()
{
  int i,j,k;
  int bestCount;


  signal(SIGINT,cnt);

  //Initalize everything to zero and precision PREC
  mpf_set_default_prec(PREC);
  mpf_init(WINDOW_SIZE);
  
  
  for(i=0;i<SAMPLE_SIZE;i++)
    {
      mpf_init(RAND_AVES[i]);
      for(j=0;j<DEG;j++)
	mpf_init(RAND_COEFS[i][j]);
    }
  
  for(i=0;i<POINT_NUM;i++)
    {mpf_init(NEW_ELT[i]);
    mpf_init(BEST[i]);
    }
  mpf_init(FIT);
  mpf_init(BEST_FIT);
  
  mpf_init(VAL);mpf_init(TEMP1);mpf_init(TEMP2);mpf_init(TEMP3);
  

  mpf_set_d(WINDOW_SIZE,1.0);
  
  mpf_set_d(BEST_FIT,99999999);
  for(i=0;i<POINT_NUM;i++)
    mpf_set_d(BEST[i],.5);
  
  srand(time(NULL));
  //init population
  polyGen();
  integrate();
  for(i=0;i<MAX_DEPTH;i++)
    {
      
      for(j=0;j<TRYS;j++)
	{
	  evalPointGen();
	  evalSingle();
	  if(findFittest())
	    {
	      j=0;
	      mpf_out_str(NULL,10,10,BEST_FIT);
	      printf("::\n");
	      for(i=0;i<POINT_NUM;i++)
		{
		  mpf_out_str(NULL,10,30,BEST[i]);
		  if(SPECTRAL)
		    {
		      gmp_printf("\n");
		      mpf_set_d(TEMP1,1.0);
		      mpf_sub(TEMP2,TEMP1,BEST[i]);
		      mpf_out_str(NULL,10,30,TEMP2);
		    }  
		  gmp_printf("\n");
		}
	      if(CENTER_POINT)
		gmp_printf("%f",.5);
	      printf("\n\n");
	    }
	}
      cnt(1);      
    }
  printf("Best Points:\n");

  if(CENTER_POINT)
    {
      mpf_set_d(TEMP1,.5);
      mpf_out_str(NULL,10,30,TEMP1);
    }
  for(i=0;i<POINT_NUM;i++)
    {
      mpf_out_str(NULL,10,30,BEST[i]);
      if(SPECTRAL)
	{
	  mpf_set_d(TEMP1,1);
	  mpf_sub(TEMP2,TEMP1,BEST[i]);
	  mpf_out_str(NULL,10,30,TEMP2);
	}
      printf("\n\n");
    }
  
  return 0;
}