#include "Model.hh"

double Model::dN_dOmega_dlogE( Component* component, SearchEvent& searchevent, bool allow_zero   ) const
{
  double p;
  {
    p = component->dN_dOmega_dlogE( searchevent );
  }

  if ( !allow_zero && p <= 0 ) {
    cout << searchevent.str() << endl;
    cout << info() << endl;
    cout << "component =" << component->IsA()->GetName();
    fatal ( "p at or below 0 for the above searchevent ", p );
  }
  if ( ::isnan(p) )
  {
    cout << searchevent.__str__() << endl;
    cout << "        p = " << p << endl;
    cout << "component = " << component->IsA()->GetName() << endl;
    cout << info() << endl;

    fatal ( "p is nan for the above searchevent " );
  }
  return p;
}


double Model::extended_log_likelihood ( DataSet& dataset )
{
  TIMEFUNC
  MEMCHECK

  double L = 0;

  for ( auto searchevent : dataset.data )
  {
    double P = 0;

    for ( auto component : components )
    {
      //Timer _("extended_log_likelihood::dN_dOmega_dlogE");

      P += dN_dOmega_dlogE( component, searchevent , true );
    }


    if ( P <= 0 ) {
      cout << "for searchevent" << searchevent.str() << ": " << endl;
      cout << ":zero probability or smaller. (P = " << P << ") ==> force L += 0 to avoid error in log calculation" << endl;
      // cout << searchevent.__str__() << endl;
      // cout << info() << endl;
      L += -20;  // NB: += 0 is adding a very high probability as if log(P)=0: the P = 1
    }

    if ( P > 0 ) {
      L += log(P);
    }

  }


  for ( auto component : components )
  {
    // Timer _("ntot");

    L -= Ntot( component );
  }


  // if ( true ) {
  //   cout << "likelihood at  ";
  //   const auto& pars = parameters();
  //   for (auto& p : pars ) {
  //     cout << p->value << "  = ";
  //   }
  //   cout << L << endl;
  // }


  return L;
}


void Model::fit( DataSet& dataset )
{
  TIMEFUNC
  MEMCHECK

  init(dataset); // initialize all components for dealing with this dataset.

  vector<Parameter*> params = parameters();
  ncalls = 0.;

// root minimizer stuff is a mess.
// ROOT::Math::Minimizer* M =
// ROOT::Math::Factory::CreateMinimizer( minimizer_name, minimizer_algo );
// auto M = new ROOT::Math::GSLMinimizer( ROOT::Math::kVectorBFGS2 );
// auto M = new ROOT::Minuit2::Minuit2Minimizer();
  //TMinuitMinimizer

  auto M = new TMinuitMinimizer();

  if (!M)
  {
    fatal("failed to create minimizer");
  }

// M->SetMaxFunctionCalls(10000);   // for Minuit/Minuit2
// M->SetMaxIterations(100);        // for GSL
// M->SetTolerance(0.01);
// M->SetPrintLevel( -1 );          // -1 seems to be the lowest setting in Gslminimizer

  if ( use_binned_likelihood )
  {
    dataset.prepare_histogram();
  }



  auto Lfunc = use_binned_likelihood ? &Model::extended_log_likelihood_binnned : &Model::extended_log_likelihood;

  auto score_function = [&] ( const double * pars )
  {
    for (unsigned i = 0; i < params.size() ; i++ )
    {
      params[i]->value = pars[i];
    }
    update_components();

    double r  =  -(this->*Lfunc)( dataset );
    ncalls++;
    return r;
  };

  ROOT::Math::Functor f( score_function, params.size() );
  M->SetFunction( f );

  for (unsigned i = 0 ; i < params.size() ; i++ )
  {
    Parameter* p = params[i];

    if ( p->fit )
    {
      M -> SetVariable( i, p->fullname() , p->fit_start, p->fit_step );
      M -> SetVariableLimits(i, p->fit_min , p->fit_max );
    }
    else
    {
      M -> SetFixedVariable( i, p->fullname() , p->fit_start );
    }
  }

  M -> Minimize(); // for great justice!

  for (unsigned i = 0 ; i < params.size() ; i++ )
  {
    //params[i]->fit_result   = M->X()[i];
    params[i]->value        = M->X()[i];
    params[i]->fit_error_up = M->Errors()[i]; // depricated
    params[i]->fit_error    = M->Errors()[i];
  }
  fitted_likelihood =  -M -> MinValue();

}