#include "DataSet.hh"
#include "Model.hh"
#include "TString.h"

string DataSet::__str__() const
{
  Table T("index", "flavor", "channel", "source", "time", "zenith", "ra", "dec", "dir-error", "Erec", "Erec-error");

  for (auto& sev : data )
  {
    T <<
      sev.index <<
      sev.flav << 
      sev.chan <<
      sev.source <<
      sev.t.AsString("s") <<
      sev.zenith << 
      sev.coordinates.ra_deg() <<
      sev.coordinates.dec_deg() <<
      sev.direction_uncertainty <<
      sev.Eproxy <<
      sev.Eproxy_uncertainty;
  }
  return T.__str__();
}


void DataSet::prnt( Model& M ) 
{
  Table T("index", "channel", "source", "time", "zenith", "ra", "dec","Erec");

  for (auto& c : M.components ) {
    T.add_column( c->name );
  }

  for (auto& sev : data )
  {
    T <<
      sev.index <<
      sev.chan <<
      sev.source <<
      sev.t.AsString("s") <<
      sev.zenith << 
      sev.coordinates.ra_deg() <<
      sev.coordinates.dec_deg() <<
      sev.Eproxy;
      
    for (auto& c : M.components ) {
      T << c->dN_dOmega_dlogE( sev );
  }
  
  } 
  cout << T.__str__();
}


template<typename T>
void purge( vector<T>& v, const T& item )
{
  v.erase(std::remove(v.begin(), v.end(), item ), v.end());
}


void DataSet::report_histogram(bool org )
{
  int bb = 0;
  for (int by = 1 ; by <= H_nevents->GetNbinsY() ; by ++ )
  for (int bx = 1 ; bx <= H_nevents->GetNbinsX() ; bx ++ )
  {
      int b  = H_nevents->GetBin( bx,by );
      int n  = H_nevents->GetBinContent( b );

      if ( n== 0 ) continue;

      double minx = H_nevents->GetXaxis()->GetBinLowEdge( bx );
      double maxx = H_nevents->GetXaxis()->GetBinWidth( bx ) + minx;

      double miny = H_nevents->GetYaxis()->GetBinLowEdge( by );
      double maxy = H_nevents->GetYaxis()->GetBinWidth( by ) + miny;

      double bbb = org ? b : bb;

      print (b, bb, "|", bx, by, "|", n, nevents[bbb], " | ", minx, maxx, mean_eproxy[bbb], " | ", miny, maxy, mean_costheta[bbb] );
      bb++;
    
  }
}


void DataSet::prepare_histogram()
{
  print ("preparing histogram");

  if ( H_nevents == nullptr )
  {
      H_nevents = new TH2D("H_nevents","H_nevents",70,1,8,40,-1,1 );
  }

  H_nevents->Reset();

  int N = (H_nevents->GetNbinsX()+2) * (H_nevents->GetNbinsY()+2);

  nevents.resize(N);
  mean_eproxy.resize(N);
  mean_costheta.resize(N);

  std::fill(nevents.begin(), nevents.end(), 0);
  std::fill(mean_eproxy.begin(), mean_eproxy.end(), 0.0);
  std::fill(mean_costheta.begin(), mean_costheta.end(), 0.0);

  //print ("histogram vectors",nevents.size(),mean_eproxy.size(), mean_costheta.size());

  for (auto& sev : data )
  {
    double bx = H_nevents -> GetXaxis() -> FindBin ( sev.Eproxy );
    double by = H_nevents -> GetYaxis() -> FindBin ( cos( sev.zenith ));

    if ( bx == 0 || bx > H_nevents -> GetXaxis() -> GetNbins()  )
    {
      print ("Warning : energy out of bounds ", sev.Eproxy );
    }
  
    if ( by == 0 || by > H_nevents -> GetYaxis() -> GetNbins()  )
    {
      print ("Warning : cos-zenith out of bounds ", cos( sev.zenith ) );
    }

    int b  = H_nevents->GetBin( bx,by );

    double ee = H_nevents -> GetXaxis() -> GetBinCenter(bx);
    double zz = H_nevents -> GetYaxis() -> GetBinCenter(by);

    // we can either use the bin-center, or the actual everage of the
    // Erec & costheta inside the bin.

    mean_eproxy  [b]+= ee; // sev.Eproxy;
    mean_costheta[b]+= zz; // cos( sev.zenith );
    H_nevents->SetBinContent( b, H_nevents->GetBinContent(b)+1 );
  }


  for (int bx = 1 ; bx <= H_nevents->GetNbinsX() ; bx ++ )
    {
      for (int by = 1 ; by <= H_nevents->GetNbinsY() ; by ++ )
      {
        int b  = H_nevents->GetBin( bx,by );
        int n  = H_nevents->GetBinContent( b );
        nevents[b] = n;

        if (n == 0 )
        {
          mean_eproxy  [b] = 0;
          mean_costheta[b] = 0;
        }
        else
        {
          mean_eproxy  [b] /= n;
          mean_costheta[b] /= n;
        }

        // print ( b, nevents[b], mean_eproxy[b], mean_costheta[b]);
      }
    }


  // print ("histogram vectors",nevents.size(),mean_eproxy.size(), mean_costheta.size());
  // report_histogram(true );


  // remove all zeros from those vectors;
  purge ( nevents,       0   );
  purge ( mean_costheta, 0.0 );
  purge ( mean_eproxy,   0.0 );

  // report_histogram();

  if ( nevents.size() != mean_eproxy.size() || nevents.size() != mean_costheta.size() )
  {
    fatal("error filling histogram vectors",nevents.size(),mean_eproxy.size(), mean_costheta.size());
  }

  // print (" data sumamrized in ", nevents.size() , " non-empty bins ");
}


SkyMap DataSet::skymap()
{
  SkyMap s;

  TMarker marker;
  marker.SetMarkerSize ( 0.5 );

  for ( auto& sev : data )
  {
    marker.SetMarkerColor( 1  + sev.source  );
    marker.SetMarkerStyle( 20 + sev.chan );
    s.add( marker, sev.coordinates );
  }

  return s;
}