#include "JppGrad.hh"
#include "stringutil.hh"

JppGrad::JppGrad( vector< string > shower_pdfs,
		  double sigma_blur          ,// = 0, 
		  double R_bg                ,// background rate in [Hz]
		  double t_start             ,// start of hit time residual window 
		  double t_end,
		  uint verbose,
		  uint debug ): _background_rate( R_bg ),
				_window_start( t_start ),
				_window_end( t_end ),
				_verbose( verbose ),
				_debug( debug ),
				_dx_dtes(),
				_dXs( 20, vector< vector < double > >( 3, vector< double >( 6, 0.0 ) ) ),
				_sum_dnpes_a( 7, vector<double>( 3, 0 ) ),
				_sum_dXs_a  ( 7, vector<double>( 3, 0 ) ),
  _sum_dTs_a  ( 7, vector<double>( 6, 0 ) ),
  _sum_dTes_a ( 7, vector<double>( 7, 0 ) ),
  _sum_dnpes_n( 7, vector<double>( 3, 0 ) ),
  _sum_dXs_n  ( 7, vector<double>( 3, 0 ) ),
  _sum_dTs_n  ( 7, vector<double>( 6, 0 ) ),
  _sum_dTes_n ( 7, vector<double>( 7, 0 ) ),
  _dL_a( 7, 0.0 ),
  _dL_n( 7, 0.0 ){
  
  for( uint ii = 0; ii != 7; ii++ ){
    uint iii = ii >= 4 ? ii + 3 : ii; //Skip over the cart direction dims returned by jacobian objects
    kk_kki[ii] = iii;
    jj_jji[ii] = iii;
  };


  bool build = !file_exists("totalpdf.dat");


  if (build)
  { 
    cout << "Building JppGrad pdf from direct and scattered constituents..." << endl;

    cout << shower_pdfs[0] << endl;
    cout << shower_pdfs[1] << endl;
    cout << sigma_blur << endl;
    cout << R_bg << endl;
    cout << t_start << endl;
    cout << t_end << endl;

    Timer tim("constructPDF", true);

    _pdf_gradfuncs.resize(2);
    _npe_gradfuncs.resize(1);
      
    print ("loading scattered photon table.");
    _pdf_gradfuncs[0].load( shower_pdfs[0].c_str() );
    
    print ("loading direct photon table.");
    _pdf_gradfuncs[1].load( shower_pdfs[1].c_str() );
    
    print ("computing pdfs");

    for( vector< JppShowerGradPdf_t >::iterator ipdf = _pdf_gradfuncs.begin(); ipdf != _pdf_gradfuncs.end(); ++ipdf ) 
    {
      std::cout << "x" << std::flush;

      ipdf->transform(JppShowerGradPdf_t::transformer_type::getDefaultTransformer()); // get rid of R-dependent weight function
      ipdf->setExceptionHandler( new JppShowerGradPdf_t::JDefaultResult( JMATH::zero ) );

      auto x = new PDG_Function1D_t::JDefaultResult( JMATH::zero ) ;
      for (JppShowerGradPdf_t::super_iterator p = ipdf->super_begin(); p != ipdf->super_end(); ++p) 
      { 
        p.getValue().setExceptionHandler(x);
      }

    };

    print ("adding pdfs");
    
    add_jpp_shower_pdfs( _pdf_gradfuncs[0], _pdf_gradfuncs[1] );
    
    print("Making npe table.");
    _npe_gradfuncs[0] = JppShowerGradNpe_t( _pdf_gradfuncs[0] ); //Make NPE table for total only
    
    print("smearing combined table with sigma=", sigma_blur );
    if ( sigma_blur > 0 ){
      _pdf_gradfuncs[0].blur( sigma_blur ); // only blur total : save some time
    }

    print("background rate:", R_bg, "Hz");
    print("time window:", t_start, t_end );
    print("JppGrad constructed in", tim.value(), "sec");

    _pdf_gradfunc = &_pdf_gradfuncs[ 0 ]; //Could remove pdf_funcs vector if we really don't need components
    _npe_gradfunc = &_npe_gradfuncs[ 0 ];

    print ("saving (?)");
    _pdf_gradfunc->store("totalpdf.dat");
  } 
  else  // don't build, but load
  {
    print("loading totalpdf.dat");
    _pdf_gradfunc = new JppShowerGradPdf_t;
    _pdf_gradfunc->load("totalpdf.dat");
    _pdf_gradfunc->transform(JppShowerGradPdf_t::transformer_type::getDefaultTransformer()); // get rid of R-dependent weight function
    _pdf_gradfunc->setExceptionHandler( new JppShowerGradPdf_t::JDefaultResult( JMATH::zero ) );

    print("building npe table");
    _npe_gradfunc = new JppShowerGradNpe_t( *_pdf_gradfunc );
  }
  print ("done with JppGrad ctor.");

};




void JppGrad::eval_lik_gradient(double &r, vector<double> &ds, const Trk &trk, vector<Hit> hits, vector<Pmt> unhits)
  {
    TIMEFUNC

    // Timer time0("time0");

    // out debug : 1.66632e-07
    // set elong : 0.000437399
    // fill ds   : 3.44014e-07
    // hit loop  : 0.00382506
    // unhit loop: 1.06737e-07
    //----------inner----------
    // thit norm      : 1.72007e-07
    // thit paramst   : 3.31191e-06
    // thit jacob     : 0.000140122
    // thit A         : 0.00121552
    // thit dnpedtdts : 0.00123817
    // thit dnpedts   : 0.00122572
    // thit dte       : 1.07735e-06
    // thit dhitL     : 2.48028e-07
    // thit sanity    : 2.61082e-07
    // thit debug     : 1.55114e-07
    // tunhit norm    : 0
    // tunhit params  : 0
    // tunhit jacob   : 0
    // tunhit dnpedts : 0
    // tunhit dunhitL : 0

    // double t0 = time0.time_since_last_start();

    if (_debug)
    {
      bad_dXs = 0;
      bad_int_dXs = 0;
      bad_dtdXs = 0;
      call_dXs = 0;
      call_int_dXs = 0;
      call_dtdXs = 0;
      pdf_call_count = 0;
      npe_call_count = 0;
      bad_pdf_call_count = 0;
      bad_npe_call_count = 0;
      if (_debug > 1)
      {
        _sum_dnpes_a.assign(7, vector<double>(3, 0));
        _sum_dXs_a.assign(7, vector<double>(3, 0));
        _sum_dTs_a.assign(7, vector<double>(6, 0));
        _sum_dTes_a.assign(7, vector<double>(7, 0));

        _sum_dnpes_n.assign(7, vector<double>(3, 0));
        _sum_dXs_n.assign(7, vector<double>(3, 0));
        _sum_dTs_n.assign(7, vector<double>(6, 0));
        _sum_dTes_n.assign(7, vector<double>(7, 0));

        _dL_a.assign(7, 0.0);
        _dL_n.assign(7, 0.0);
      };
    };

    // double t1 = time0.time_since_last_start();

    { 
      Timer _ ("set_elon");
      set_elongation(trk);
    }

   
    ds.assign(_N, 0.0);
    r = 0.0;


    foreach (hit, hits)
    {
      TIMER(hitloop)

      // Here we access the pdf for the first and only time
      set_paramst(hit); // 0.96%, 1.74e-7

      if (hit.pattern_flags == 2)
      {
        continue; // Hit flagged as bad
      }

      set_jacobians(hit); // 39.1%, 7.1e-6

      // double t7 = time0.time_since_last_start();
      // thit_jacob += t7 - t6;

      const vector<double> dnpe_dtdTs = this->dnpe_dtdTs(hit);

      // double t8 = time0.time_since_last_start();
      // thit_dnpedtdts += t8 - t7;
      vector<double> dnpe_dTs;

      if (npe(hit) < _fudge2)
        dnpe_dTs = this->dnpe_dTs(hit); // 26.0%, 4.7e-6
      else
        dnpe_dTs.assign(7, 0.0);

      // double t9 = time0.time_since_last_start();
      // thit_dnpedts += t9 - t8;

      _A = dnpe_dt(hit);

      // double t10 = time0.time_since_last_start();
      // thit_A += t10 - t9;
      r -= get_hit_L(_A, npe(hit), true);

      // double t10b = time0.time_since_last_start();
      // thit_hitL += t10b - t10;

      for (uint ii = 0; ii != _N; ii++) // gradient
      {
        // double t11 = time0.time_since_last_start();

        double dt_e = 0;
        for (vector<vector<vector<double>>>::iterator idTs = _dX_dTes.begin(); idTs != _dX_dTes.end(); idTs++)
        {
          dt_e += (*idTs)[5][ii] / (float)_elong_trks.size();
        };

        // double t12 = time0.time_since_last_start();
        //	thit_dte += t12 - t11;

        ds[ii] -= dhit_L(dnpe_dtdTs[ii], dnpe_dTs[ii], _A, dt_e);

        // double t13 = time0.time_since_last_start();
        // thit_dhitL += t13 - t12;

        if (ds[ii] < -1e6 or ds[ii] != ds[ii])
        {
          cout << "ds[ii]: " << ds[ii] << endl;
          cout << "bad LL grad ii: " << ii << endl;
          cout << "dnpe_dtdTs:  " << dnpe_dtdTs[ii] << endl;
          cout << "dnpe_dTs: " << dnpe_dTs[ii] << endl;
          cout << "dt_e: " << dt_e << endl;
          cout << "_A: " << _A << endl;
          cout << "trk: " << trk << endl;
          cout << "hit: " << hit << endl;
          if (_debug)
            cin.ignore();
        };

        // double t14 = time0.time_since_last_start();
        // thit_sanity += t14 - t13;
        if (_debug > 1)
        {
          print("debug");
          fill_components(trk, hit); // Must come first as the trk is set

          double dnpe_dtdT = this->num_dnpe_dtdTs(trk, hit, 1e-5)[ii];
          double dnpe_dT;
          if (npe(hit) < _fudge2)
            dnpe_dT = this->num_dnpe_dTs(trk, hit, 1e-5)[ii];
          else
            dnpe_dT = 0.0;

          double num_dt_e = 0.0;
          for (vector<Trk>::iterator itrk = _elong_trks.begin(); itrk != _elong_trks.end(); ++itrk)
          {
            num_dt_e += num_dX_dTs(*itrk, hit, 1e-5)[5][ii] / (float)_elong_trks.size();
          };

          _dL_n[ii] -= dhit_L(dnpe_dtdT, dnpe_dT, _A, num_dt_e);
        }
        // double t15 = time0.time_since_last_start();
        // thit_debug += t15 - t14;
      }
    }

    // double t16 = time0.time_since_last_start();
    /*
    double tunhit_norm = 0.0;
    double tunhit_params = 0.0;
    double tunhit_jacob= 0.0;
    double tunhit_dnpedts = 0.0;
    double tunhit_dunhitL = 0.0;
    double tunhit_unhitL = 0.0;
    */
    foreach (unhit, unhits)
    {
      TIMER(unhit)
      
      set_params(unhit); // has to be first

      if (unhit.flag == 2)
        continue; // Hit flagged as bad

      // double t19 = time0.time_since_last_start();
      // tunhit_params += t19 - t18;
      set_jacobians(unhit);
      // double t20 = time0.time_since_last_start();
      // tunhit_jacob += t20 - t19;
      const vector<double> dnpe_int_dTs = this->dnpe_int_dTs(unhit);
      // double t21 = time0.time_since_last_start();
      // tunhit_dnpedts += t21 - t20;

      r -= -npe(unhit);

      // double t21c = time0.time_since_last_start();
      // tunhit_unhitL += t21c - t21;

      for (uint ii = 0; ii != _N; ii++)
      {

        // double t21a = time0.time_since_last_start();

        double dL_unhit = -dnpe_int_dTs[ii];
        ds[ii] -= dL_unhit;

        // double t21b = time0.time_since_last_start();
        // tunhit_dunhitL += t21b - t21a;

        if (ds[ii] != ds[ii])
        {
          cout << "nan encountered in ll grad ii: " << ii << endl;
          cout << "trk: " << trk << endl;
          cout << "unhit: " << unhit << endl;
          if (_debug)
            cin.ignore();
        };

        if (_debug > 1)
          _dL_n[ii] -= -this->num_dnpe_dTs(trk, unhit, 1e-5)[ii]; // Yes, 2 negative signs
      };

      if (_debug > 1)
        fill_components(trk, unhit);
    }

    if (_debug)
    {
      cout << "Good pdf calls: " << pdf_call_count - bad_npe_call_count << ", bad: " << bad_npe_call_count << endl;
      cout << "Good npe calls: " << npe_call_count - bad_pdf_call_count << ", bad: " << bad_pdf_call_count << endl;
      // cout << "Good calls to dnpe_int_dXs: " << call_int_dXs - bad_int_dXs << ", bad: " << bad_int_dXs << endl;
    }

    // double t23 = time0.time_since_last_start();

    /*
      cout << "out debug : " << t1 - t0 << endl;
      cout << "set elong : " << t2 - t1 << endl;
      cout << "fill ds   : " << t3 - t2 << endl;
      cout << "hit loop  : " << t16 - t3  << endl;
      cout << "unhit loop: " << t23 - t16 << endl;
      cout << "grand tot : " << t23 - t0 << endl;

      cout << "----------inner----------" << endl;

      cout <<  "thit norm      : " << thit_norm << endl;
      cout <<  "thit paramst   : " << thit_paramst << endl;
      cout <<  "thit jacob     : " << thit_jacob << endl;
      cout <<  "thit A         : " << thit_A << endl;
      cout <<  "thit dnpedtdts : " << thit_dnpedtdts << endl;
      cout <<  "thit dnpedts   : " << thit_dnpedts << endl;
      cout <<  "thit dte       : " << thit_dte << endl;
      cout <<  "thit dhitL     : " << thit_dhitL << endl;
      cout <<  "thit dhitL     : " << thit_hitL << endl;
      cout <<  "thit sanity    : " << thit_sanity << endl;
      cout <<  "thit debug     : " << thit_debug << endl;

      cout <<  "tunhit norm    : " << tunhit_norm << endl;
      cout <<  "tunhit params  : " << tunhit_params << endl;
      cout <<  "tunhit jacob   : " << tunhit_jacob << endl;
      cout <<  "tunhit dnpedts : " << tunhit_dnpedts << endl;
      cout <<  "tunhit dunhitL : " << tunhit_dunhitL << endl;
      cout <<  "tunhit unhitL : " << tunhit_unhitL << endl;

      cin.ignore();
    */
  };