#ifndef MODEL_HH_INCLUDED
#define MODEL_HH_INCLUDED

#include "Component.hh"
//#include "Minuit2/Minuit2Minimizer.h"
#include "TMinuitMinimizer.h"

/*!
  A model to describe a set of neutrino telescope data. A model has a list
  of Components, which together make up the full description of model.
  A model owns its Components.  
*/

struct Model
{
  string name;
  vector<Component*> components;

  double fitted_likelihood = -1;
  int ncalls = -1;
  bool use_binned_likelihood;


  Model(string name = "unnamed") :
    name(name),
    use_binned_likelihood(false) {}


  Model(const Model& other) : Model( "copy_of_" + other.name )
  {
    components.clear();

    for ( auto c : other.components )
    {
      components.push_back( c->clone() );
    }
  }

  Model clone( string newname = "" )
  {
    Model r(*this);
    if (newname != "" ) r.name = newname;
    else                r.name = "clone_of_" + name;
    return r;
  }

  virtual ~Model()
  {
    // this should be done, but some components are still owned by python
    // for( auto c : components ) delete c;
    components.clear();
  }


  Table table_info() const
  {
    vector<string> names;
    const auto& pars = parameters();

    for (auto p : get_all_datamembers( *pars[0] ) )
    {
      names.push_back( p[1] );
    }

    Table T( split("component param description | truth value error | fit fit_start fit_min fit_max fit_step") );

    for (auto par : pars )
    {
       T << par->component_name << par->name << par->description << "|" 
         << par->true_value     << par->value << par->fit_error << "|"
         << par->fit << par->fit_start << par->fit_min << par->fit_max << par -> fit_step;
    }

    return T;
  }

  string info() const
  {
    Table T = table_info();
    string r = "  Model:" + name + "  ";
    r += "lik=" + str(fitted_likelihood ) + "  ";
    r += "iterations=" + str(ncalls ) + " ";
    T.title = r;
    return T.str();
  }




  /*! return the (first) component with a certain name */

  Component* get_component(string name)
  {
    for ( auto c : components ) if ( c->name == name ) return c;
    print ("no component found with name ", name );
    return nullptr;
  }

  /* return the (first) component with a certain class name */

  Component* get_component_type(string typname)
  {
    for ( auto c : components ) if ( c->IsA()->GetName()  == typname ) return c;
    print ("no component found with type ", typname );
    return nullptr;
  }


  void remove_component(string name)
  {
    vector<Component*> newlist;

    for (unsigned i = 0; i < components.size(); i++)
    {
      auto& c = *components[i];
      if ( c.name == name ) delete &c;
      else newlist.push_back( &c );
    }
    components = newlist;
  }

  void add_component( Component& c )
  {
    components.push_back( &c );
  }


  /*! call update on all of our components; do this after each parameter change */
  void update_components(){
    for( auto c : components ) c -> update();
  }


  /*! set_true values of each paramters accoring to src */

  void set_truth( Model& src ) 
  {
      for ( auto couple : hzip ( parameters(), src.parameters() )) 
        {
          if (couple[0]->name != couple[1]->name ) fatal( "non-identical paramters in set_truth" );
          couple[0]->true_value = couple[1]->value;
        }
  }



  /* initialize each component so that they may do pre-computations. This will be
     called by the fit method. */

  void init(DataSet& dset )
  {
    for ( auto c : components ) c->init( dset );
  }


  DataSet generate_dataset( bool fluctuate_number_of_events = true, int seed = 0 )
  {
    DataSet r;
    fill_dataset( r, fluctuate_number_of_events, seed );
    return r;
  }

  Component* get_background_component() {

    Component* r = get_component_type("BackgroundComponent");
    if (!r) {
      fatal ("no background component defined while generating events");
    }
    return r;
  }


  void fill_dataset_with_component( DataSet& dset, 
                                    Component* c,
                                    bool fluctuate_number_of_events = false )
  {
      c->deinit() ; // reset any caches.

      unsigned ngen_exact = fluctuate_number_of_events ? gRandom->Poisson( c->norm ) : c->norm;

      print ("generating ", ngen_exact,"for", c->name );

      for (unsigned i = 0 ; i < ngen_exact ; i++)
      {
        auto ev = c->generate_event();

        if ( ev.chan == undefined_channel ) {
          fatal ("component", c->name , "generated event without channel");
        }

        // each event must have a non-zero background probability
        // (check e.g. that no events violate the zenith-angle cut / no weird energies are generated)

        const double p = get_background_component() -> dN_dOmega_dlogE( ev );

        if ( ! (std::isfinite(p) && p>0) ) {
          print ("Event has non-finite probablity for background!", p );
          print( ev.__str__() );
          fatal("exiting");
        }

        // another sanity check.

        if ( !ev.check_yourself_before_you_wreck_yourself() )
        {
          print(ev.__str__() );
          fatal("Component", name, "generated an invalid searchevent ")
        }

        // if all good , put in the dataset.
        ev.index = dset.data.size();
        ev.source = c->component_id;
        dset.data.push_back( ev );
      }


  }


  void fill_dataset( DataSet& dset, 
                     bool fluctuate_number_of_events, 
                     int seed = 0 )
  {
    dset.clear();
    dset.seed = seed;
    if ( seed ) gRandom->SetSeed( seed );
    
    for ( auto c : components ) 
      {
        fill_dataset_with_component( dset, c, fluctuate_number_of_events );                                                        
      }
  }




  /*! add some error checking to calling component->dN_dOmega_dlogE */

  double dN_dOmega_dlogE( Component* component, SearchEvent& searchevent, bool allow_zero = false  ) const;


  /*! call component->Ntot with error checking */

  double Ntot( Component* component )
  {
    double n = component->Ntot();
    if (::isnan(n)) fatal (" nan for ntot of ", component );
    return n;
  }

  /*! Experimental binned likelihhood for diffuse analyses as function
      for E_rec,costheta_rec. The idea is that all events in a given (Erec,ctrec)
      bin will heave the ~same dN_dOmega_dlogE for all components. It is up
      to the user to ensure/verify this is a good assumption.*/

  double extended_log_likelihood_binnned( DataSet& dataset )
  {

    static int z = 0;
    if (!z)
    {
      print ("Hi there, this is your friendly neighbourhood binned likelihood function talking!");
      z = 1;
    }

    double L = 0;
    SearchEvent searchevent;

    for (unsigned ii = 0 ; ii < dataset.nevents.size() ; ii ++  )
    {
      const int nevents = dataset.nevents[ii];

      searchevent.Eproxy = dataset.mean_eproxy[ii];
      searchevent.zenith = acos( dataset.mean_costheta[ii] );

      double P = 0;
      for ( auto component : components )
      {
        P += dN_dOmega_dlogE( component, searchevent );
      }

      L += nevents * log ( P );
    }

    for ( auto component : components )
    {
      L -= Ntot( component );
    }

    return L;
  }

  double extended_log_likelihood ( DataSet& dataset );


  /*! all parameters in the model */
  vector<Parameter*> parameters() const
  {
    vector<Parameter*> r;

    for ( auto component : components )
    {
      const auto& v = component->parameters();
      r.insert( r.end(), v.begin(), v.end() );
    }
    return r;
  }


  /* define a tree with a branch for every parameter */
  void add_to_tree( TTree* T )
  {
    T->Branch( ( name + ".fitted_likelihood" ).c_str() , &fitted_likelihood, "fitted_likelihood/D" ) ;
    T->Branch( ( name + ".ncalls"            ).c_str() , &ncalls,            "ncalls/I" ) ;

    for ( auto p : parameters() )
    {
      T->Branch( ( name + "." + p->fullname()          )  .c_str() , & ( p->value ),      (p->fullname() + "/D"     ) .c_str() );
      T->Branch( ( name + "." + p->fullname() + "_error") .c_str() , & ( p->fit_error  ), (p->fullname() + "_error/D").c_str() );
      T->Branch( ( name + "." + p->fullname() + "_true")  .c_str() , & ( p->true_value ), (p->fullname() + "_true/D") .c_str() );
    }
  }

  void fit( DataSet& dataset );

  ClassDef( Model, 1 );
};

namespace cling
{
inline std::string printValue(const Model* m)
{
  return m->info();
}
}

// inline void test_comp()
// {
//   Model m("dus");

//   GaussianPointSource g;
//   FlatBackground f;

//   m.components.push_back( &f );
//   m.components.push_back( &g );

//   cout << m.info() << endl;
//   Model mm = m.clone();
//   cout << mm.get_component("gauspoint")->IsA()->GetName() << endl;
// }

#endif