#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cstdlib>
#include <cmath>

#include "Eoe.h"
#include "fitPar.h"
#include <TRandom3.h>
#include <TFile.h>
#include <TNtuple.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TPad.h>
#include <TH1D.h>
#include <TMinuit.h>
#include <TGraph.h>
#include <TPostScript.h>

Eoe* eoeGlobal;  
bool debug;

int fitPar (Eoe* eoe, vector<bool> freePar){
  
// set up start values, step sizes, etc. for fit

  eoeGlobal = eoe;                  // communicate to fcn via global
  int npar = eoe->npar();
  int N = eoe->N();
  vector<double> parVec = eoe->parVec();
  vector<double> sigParVec = eoe->sigParVec();

  const int maxpar = 100;
  double par[maxpar];
  string parName[maxpar];
  double stepSize[maxpar];
  double minVal[maxpar];
  double maxVal[maxpar];

  eoe->setStartPar();               // sets m_startPar, m_stepSize
  for (int i=0; i<npar; i++){
    par[i] = eoe->startPar(i);
    stepSize[i] = eoe->stepSize(i);
  }

  for (int i=0; i<npar; i++){
    stringstream ss;
    string iString;
    ss << i;
    ss >> iString;
    parName[i] = "p" + iString;
  }

  for (int i=0; i<npar; i++){    // some of these may be overridden below
    minVal[i] = 0.;
    maxVal[i] = 0.;
  }

  //   TMinuit* minuit = eoe->minuit();
  TMinuit* minuit = new TMinuit();
  minuit->mncler();

  bool verbose = true;
  if ( verbose ) {
    minuit->SetPrintLevel(0);
  }
  else {
    minuit->SetPrintLevel(-1);
    // cout << "set print level"  << endl;
    int ierr = 0;
    minuit->mnexcm("SET NOWarnings", 0, 0, ierr);
    // cout << "set no warnings" << endl;
  }

  double up = eoe->up();
  minuit->SetErrorDef(up);
  //   if ( up != 1.0 ) {  
  //   cout << "fitPar:  up = " << up << endl;
  // }

  minuit->SetFCN(fcn);
  for (int i=0; i<npar; i++){
    minuit->DefineParameter(i, parName[i].c_str(), 
    par[i], stepSize[i], minVal[i], maxVal[i]);
  }

  for (int i=0; i<npar; i++){
    if ( !freePar[i] ){
      minuit->FixParameter(i);
    }
    else {
      minuit->Release(i);
    }
  }

  // check the data
  // for (int i=0; i<N; i++){
  //    cout << "fitPar: i, y, u, v: " << i << "  " << eoe->y(i) << "  "
  // 	 << eoe->u(i) << "  " << eoe->v(i) << endl;
  // }

// do the fit and extract results (skip fit if no free parameters!)

  int status = 0;
  int nFreePar = minuit->GetNumFreePars();
  if ( nFreePar > 0 ) {
    status = minuit->Command("MIGRAD");
    // cout << "after migrad:  status = " << status << endl;
    // cout << "after hesse:  status = " << status << endl;
    if ( status != 0 ) {
      debug = true;
      status = minuit->Command("SIMPLEX");
      status = minuit->Command("MIGRAD");
      debug = false;
    }
  }

  parVec.clear();
  sigParVec.clear();
  double fitpar[nFreePar], err[nFreePar];
  for (int i=0; i<npar; i++){
    minuit->GetParameter(i, fitpar[i], err[i]);
    parVec.push_back(fitpar[i]);
    sigParVec.push_back(err[i]);
  }
  eoe->setParVec(parVec);
  eoe->setSigParVec(sigParVec);

  // do minos analysis for mu  (internally minuit par no. 1)

  double arglist[10];
  int ierr;
  // up set above, do not reset here
  // arglist[0] = eoe->up();        // "up"
  // if ( arglist[0] != 1.0 ) {  
  //   cout << "arglist[0] = " << arglist[0] << endl;
  // }
  // minuit->mnexcm("SET ERR", arglist, 1, ierr);

  // arglist[0] = 1.e-15;
  // minuit->mnexcm("SET EPS", arglist, 1, ierr);
  // arglist[0] = 2;
  // minuit->mnexcm("SET STR", arglist, 1, ierr);
  arglist[0] = 1000;      // max calls

  arglist[1] = 1;        // 1 for mu only
  minuit->mnexcm("MINOS", arglist, 2, ierr);   // for mu only

  // for debugging:  get minos for all parameters, not only mu
  // minuit->mnexcm("MINOS", arglist, 1, ierr);   // for all params
  // cout << "i, epl, emi, gcc,.." << endl;
  // for (int i=0; i<npar; i++){
  //   double epl = gMinuit->fErp[i];
  //   double emi = gMinuit->fErn[i];
  //   double gcc = gMinuit->fGlobcc[i];
  //   cout << i << "  " << epl << "  " << emi << "  " << gcc << endl;
  // }

  double eplus, eminus, eparab, gcc;
  minuit->mnerrs(0, eplus, eminus, eparab, gcc);  // 0 for par no. 1
  eoe->setMinosErrors(eparab, eplus, eminus);

  // Get covariance matrix

  /*
  cout << "covariance and correlation coefficients..." << endl;
  double covmat[npar][npar];
  minuit->mnemat(&covmat[0][0], npar);
  for (int i=0; i<npar; i++){
    for (int j=0; j<npar; j++){
      double sigma_i = sqrt(covmat[i][i]);
      double sigma_j = sqrt(covmat[j][j]);
      double rho = covmat[i][j]/(sigma_i*sigma_j);
      cout << i << "  " << j << "  " << covmat[i][j] << "  " << rho << endl;
    }
  }
  */

  // Scan of mu (off by default)
  double eps = 1.e-2;
  //   if ( eplus < eps || abs(eminus) < eps ) {
  if ( false ) {

    cout << "eplus, eminus, eparab = " << eplus << "  " << eminus << 
      "  " << eparab << endl;

    arglist[0] = 1;      // par no.
    arglist[1] = 50;     // num points
    arglist[2] = fitpar[0] - 1;    // from
    arglist[3] = fitpar[0] + 1;     // to
    minuit->mnexcm("SCAN", arglist, 4, ierr);

    TApplication* theApp = new TApplication("App", 0, 0);
    TCanvas* canvas = new TCanvas("canvas", "canvas", 10, 10, 500, 500);
    TGraph* gr = (TGraph*)minuit->GetPlot();
    // gr->SetMinimum(0);
    gr->Draw();
    canvas->Update();
    canvas->Draw();
    canvas->Show();
    canvas->Print("scan.gif", "gif");
    theApp->Run(true);

    string dummy;
    cout << "Enter any key to continue: ";
    getline(cin, dummy);
    delete theApp;
    delete canvas;
    delete gr;

  }       // plot scan

  if ( status != 0 ) {
    cout << "status = " << status << endl;
    cout << "y, u, s..." << endl;
    for (int i=0; i<N; i++){
      cout << eoe->y(i) << "  " << eoe->u(i) << "  " << eoe->s(i) << endl;
    }
    cout << "i, par, sig  ... " << endl;
    for (int i=0; i<npar; i++){
      cout << i << "  " << fitpar[i] << " " << err[i] << endl;
    }
    cout << endl;
  }

  delete minuit;
  return status;     // 0 = normal completion, 4 = MIGRAD not converged

}

// fcn must be a non-member function, uses global Eoe object 

void fcn(int& npar, double* deriv, double& f, double par[], int flag){

  int np = eoeGlobal->npar();
  vector<double> param;
  for (int i=0; i<np; i++){
    param.push_back(par[i]);
  }
  f = -2.*eoeGlobal->lnL(param);

  // if ( debug ) {
  if ( false ) {
    for (int i=0; i<np; i++){
      cout << "i, par = " << i << "  " << par[i] << endl;
    }
    cout << "f = " << f << endl;
    cout << endl;
  }

}