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

#include "Prof.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>

Prof* profGlobal;  
bool debug;

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

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

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

  prof->getStartPar();               // sets m_startPar, m_stepSize
  for (int i=0; i<npar; i++){
    par[i] = prof->startPar(i);
    stepSize[i] = prof->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 = prof->minuit();

  bool verbose = true;
  if ( verbose ) {
    minuit->SetPrintLevel(1);
  }
  else {
    minuit->SetPrintLevel(-1);
    // cout << "set print level"  << endl;
    int ierr = 0;
    minuit->mnexcm("SET NOWarnings", 0, 0, ierr);
    // cout << "set no warnings" << 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: " << i << "  " << prof->y(i) << "  "
  //	 << prof->u(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;
    // status = minuit->Command("HESSE");
    // 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]);
  }
  prof->setParVec(parVec);
  prof->setSigParVec(sigParVec);

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

  double arglist[10];
  int ierr;
  arglist[0] = 1.;        // "up"
  minuit->mnexcm("SET ERR", arglist, 1, ierr);
  arglist[0] = 500;      // max calls
  arglist[1] = 0;        // 0 for all, 1 for mu only?
  minuit->mnexcm("MINOS", arglist, 2, ierr);
  double eplus, eminus, eparab, gcc;
  minuit->mnerrs(1, eplus, eminus, eparab, gcc);
  prof->setMinosErrors(eparab, eplus, eminus);

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

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

}

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

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

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

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

}