#include "Math/GSLMinimizer.h"
#include "Math/Functor.h"

// Pointers to arrays holding data
Float_t *CURR_Y,*YHAT1,*LOSS1,*LOSS2,*minY, *maxY, *ORIGSEQ;
double *bestY;
Int_t ndata = 1000;


//Convert yhat to probability
double getProbability(Float_t yhat){
	return (1 - 1/(1+exp(-yhat)));
}

// The Minimizer function
double MyFunc(const double *par) {
  double totalProfit = 0;
  double totalProb=0;
  double totalEstPrem=0;
  double totalCurrEstPrem=0;
  double totalCurrProb=0;

  for (int i = 0; i < ndata; i++) {
  	
  	Float_t yhat= YHAT1[i] + 1.4858 * (par[i] - CURR_Y[i])/CURR_Y[i];
  	totalProfit += getProbability(yhat) * (par[i] - LOSS1[i] - LOSS2[i]);
    totalProb += getProbability(yhat);
    totalEstPrem += getProbability(yhat) * par[i];
    
    totalCurrEstPrem += getProbability(YHAT1[i]) * CURR_Y[i];
    totalCurrProb += getProbability(YHAT1[i]);
  }



  //If Constraints are not met, return a large number
  if (totalProb / totalCurrProb > 1.1 || 
  	  totalProb / totalCurrProb < 0.9 || 
  	  totalEstPrem/totalCurrEstPrem > 1.1 || 
  	  totalEstPrem/totalCurrEstPrem < 0.9	 ) return 0;

  //Return -totalProfit if all constraints are met
  return -totalProfit;
}

void GetData() {
  TFile *fin = new TFile("OptimizeThis.root");
	if (fin->IsZombie()) {
	   cout << "unable to read tree " << endl; 
	   break;
    }  
  Float_t effdate, polnum, curr_man_prem , yhat1, lossest,expenses,minPrem,maxPrem,origSeq;
  TTree *T = (TTree *) fin->Get("tree");

  T->SetBranchAddress("CURR_MAN_PREM",&curr_man_prem);
  T->SetBranchAddress("YHAT1",&yhat1);
  T->SetBranchAddress("LOSSEST",&lossest);
  T->SetBranchAddress("EXPENSES",&expenses);
  T->SetBranchAddress("minPrem",&minPrem);
  T->SetBranchAddress("maxPrem",&maxPrem);
  T->SetBranchAddress("_N_",&origSeq);  

  ndata = T->GetEntries();
  std::cout << "Found " << ndata << " entries in the data" << std::endl;

  for (Int_t i = 0; i < ndata; i++) {
    T->GetEntry(i);
    CURR_Y[i] = curr_man_prem;
    YHAT1[i] = yhat1;
    LOSS1[i] = lossest;
    LOSS2[i] = expenses;
    minY[i] = minPrem;
    maxY[i] = maxPrem;
    ORIGSEQ[i]=origSeq;
  }
  

}



int optimize() 
{

  // First get the data
  CURR_Y = new double[ndata];
  YHAT1 = new double[ndata];
  LOSS1 = new double[ndata];
  LOSS2 = new double[ndata];
  minY = new double[ndata];
  maxY = new double[ndata];
  ORIGSEQ = new double[ndata];

  GetData();


  // Setup the minimizer
  ROOT::Math::GSLMinimizer min( ROOT::Math::kVectorBFGS );
   
  min.SetMaxFunctionCalls(1000000);
  min.SetMaxIterations(20);
  min.SetTolerance(0.1);
  

  ROOT::Math::Functor f(&MyFunc,ndata);
  min.SetFunction(f);
  for (Int_t i = 0; i < ndata; i++) {
  	char tempString[20];sprintf(tempString,"p%i",i);
    min.SetLimitedVariable(i,tempString,floor((maxY[i]+minY[i])/20 +0.5)*10, 10,minY[i],maxY[i]);
  }

  TStopwatch StopWatch; 
  StopWatch.Reset();
  StopWatch.Start();

  
  min.Minimize(); 
  
  bestY = min.X();
  
  FILE * csvFile = fopen("optimizeResults.csv","w");
  	if (csvFile==0) {
	   cout << "unable to open output CSV file" << endl; 
	   return -1;
    } 
  
  
	fprintf(csvFile,"origSeqNo,effdate,polnum,curr_man_prem,yhat1,lossEst,expenses,minPrem,maxPrem,optimalPremium\n");

  for (Int_t i = 0; i < ndata; i++) {
  	fprintf(csvFile,"%9.0f,%10.2f,%12.9f,%10.2f,%10.2f,%10.2f,%10.2f,%10.2f\n",ORIGSEQ[i],CURR_Y[i],YHAT1[i],LOSS1[i],LOSS2[i],minY[i],maxY[i],bestY[i]);
  }
 
  fclose(csvFile);
  
  StopWatch.Stop();
  cout << "Total CPU Time = " << StopWatch.CpuTime() << " s\n"; 
  cout << "Total Real Time = " << StopWatch.RealTime() << " s\n";   
   
 
  return 1;
}