/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  ABACUS_Integ.h

Purpose:  Declares combinatorics-related classes and functions.

***********************************************************/

#ifndef ABACUS_INTEG_H
#define ABACUS_INTEG_H

#include "ABACUS.h"

namespace ABACUS {


  //**********************  Class Domain ************************

  template<class T>
    class Domain {

  private:
    Vect<T> bdry;

  public:
    Domain () : bdry(Vect<T>(2))
      {
	bdry[0] = T(0);
	bdry[1] = T(0);
      }

  public:
    Domain (T xmin_ref, T xmax_ref) : bdry(Vect<T>(2))
      {
	if (xmax_ref < xmin_ref) ABACUSerror("Use xmax > xmin in Domain.");

	bdry[0] = xmin_ref;
	bdry[1] = xmax_ref;
      }

  public:
    inline T xmin (int i)
    {
      if (i > bdry.size()/2) ABACUSerror("i index too high in Domain::xmin.");

      return(bdry[2*i]);
    }

  public:
    inline T xmax (int i)
    {
      if (i > bdry.size()/2) ABACUSerror("i index too high in Domain::xmax.");

      return(bdry[2*i + 1]);
    }

  public:
    inline int Ndomains ()
    {
      return(bdry.size()/2);
    }

  public:
    void Include (T xmin_ref, T xmax_ref) {

      // Determine the indices of xmin_ref & xmax_ref
      int xmin_reg = -1;
      int xmax_reg = -1;
      for (int i = 0; i < bdry.size(); ++i) {
	if ((i+1) % 2 && bdry[i] <= xmin_ref) xmin_reg++;
	if (i % 2 && bdry[i] < xmin_ref) xmin_reg++;
	if ((i+1) % 2 && bdry[i] <= xmax_ref) xmax_reg++;
	if (i % 2 && bdry[i] < xmax_ref) xmax_reg++;
      }

      Vect<T> new_bdry(bdry.size() + 2 * ((xmin_reg % 2 && xmax_reg % 2) - (xmax_reg - xmin_reg)/2));

      int ishift = 0;
      for (int i = 0; i <= xmin_reg; ++i) new_bdry[i] = bdry[i];
      if (xmin_reg % 2) {
	new_bdry[xmin_reg + 1] = xmin_ref;
	ishift++;
	if (xmax_reg % 2) {
	  new_bdry[xmin_reg + 2] = xmax_ref;
	  ishift++;
	}
      }
      else if ((xmin_reg + 1) % 2 && xmax_reg % 2) {
	new_bdry[xmin_reg + 1] = xmax_ref;
	ishift++;
      }
      for (int i = xmin_reg + ishift + 1; i < new_bdry.size(); ++i)
	new_bdry[i] = bdry[xmax_reg - xmin_reg - ishift + i];

      bdry = new_bdry;

      return;
    }

  public:
    void Exclude (T xmin_ref, T xmax_ref) {

      // Determine the indices of xmin_ref & xmax_ref
      int xmin_reg = -1;
      int xmax_reg = -1;
      for (int i = 0; i < bdry.size(); ++i) {
	if ((i+1) % 2 && bdry[i] <= xmin_ref) xmin_reg++;
	if (i % 2 && bdry[i] < xmin_ref) xmin_reg++;
	if ((i+1) % 2 && bdry[i] <= xmax_ref) xmax_reg++;
	if (i % 2 && bdry[i] < xmax_ref) xmax_reg++;
      }

      Vect<T> new_bdry(bdry.size()
		       + 2 * (((xmin_reg + 1) % 2 && (xmax_reg + 1) % 2) - (xmax_reg - xmin_reg)/2));

      int ishift = 0;
      for (int i = 0; i <= xmin_reg; ++i) new_bdry[i] = bdry[i];
      if ((xmin_reg + 1) % 2) {
	new_bdry[xmin_reg + 1] = xmin_ref;
	ishift++;
	if ((xmax_reg + 1) % 2) {
	  new_bdry[xmin_reg + 2] = xmax_ref;
	  ishift++;
	}
      }
      else if (xmin_reg % 2 && (xmax_reg + 1) % 2) {
	new_bdry[xmin_reg + 1] = xmax_ref;
	ishift++;
      }
      for (int i = xmin_reg + ishift + 1; i < new_bdry.size(); ++i)
	new_bdry[i] = bdry[xmax_reg - xmin_reg - ishift + i];

      bdry = new_bdry;

      return;
    }

  };

  template<class T>
    std::ostream& operator<< (std::ostream& s, Domain<T> dom)
    {
      for (int i = 0; i < dom.Ndomains(); ++i) {
	if (i > 0) s << std::endl;
	s << dom.xmin(i) << "\t" << dom.xmax(i);
      }
      return(s);
    }


  // ********************************* struct I_table ************************

  struct I_table {

    DP (*function) (DP, DP);
    int Nvals;
    DP rhomin;
    DP rhomax;
    DP alpha;
    DP logalpha;
    DP prec;
    DP* rho_tbl;
    DP* I_tbl;

    I_table (DP (*function) (DP, DP), DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec);
    DP Return_val (DP req_rho);
    void Save ();
    bool Load (DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec);

  };


  // ********************************* struct Integral_table ************************

  struct Integral_table {

    DP (*function) (DP, DP, int);
    int Nvals;
    DP rhomin;
    DP rhomax;
    DP alpha;
    DP logalpha;
    DP prec;
    int maxnrpts;
    DP* rho_tbl;
    DP* I_tbl;

    Integral_table (DP (*function) (DP, DP, int), const char* filenameprefix_ref, DP rhomin_ref,
		    DP rhomax_ref, int Nvals_ref, DP req_prec, int max_nr_pts);
    DP Return_val (DP req_rho);
    void Save (const char* filenameprefix);
    bool Load (const char* filenameprefix, DP rhomin_ref, DP rhomax_ref, int Nvals_ref, DP req_prec, int max_nr_pts);

  };


  // ******************************** Recursive integration functions ******************************

  DP Integrate_Riemann (DP (*function) (Vect_DP),
			Vect_DP& args, int arg_to_integ,
			DP xmin, DP xmax,
			int Npts);

  DP Integrate_Riemann_using_table (DP (*function) (Vect_DP, I_table),
				    Vect_DP& args, int arg_to_integ,
				    I_table Itable,
				    DP xmin, DP xmax, int Npts);

  DP Integrate_rec (DP (*function) (Vect_DP),
		    Vect_DP& args, int arg_to_integ,
		    DP xmin, DP xmax,
		    DP req_prec, int max_rec_level);

  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table),
				Vect_DP& args, int arg_to_integ,
				I_table Itable,
				DP xmin, DP xmax,
				DP req_prec, int max_rec_level);

  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table),
				Vect_DP& args, int arg_to_integ,
				I_table Itable,
				DP xmin, DP xmax,
				DP req_prec, int max_rec_level,
				std::ofstream& outfile);

  DP Integrate_rec_using_table_and_file (DP (*function) (Vect_DP, I_table, std::ofstream&),
					 Vect_DP& args, int arg_to_integ,
					 I_table Itable,
					 DP xmin, DP xmax,
					 DP req_prec, int max_rec_level,
					 std::ofstream& outfile);



  // ******************************** Recursive version:  optimal ******************************

  struct data_pt {
    DP x;
    DP f;
    DP dx;
  };

  struct Integral_result {
    DP integ_est;
    DP abs_prec;
    DP rel_prec;
    int n_vals;
  };

  std::ostream& operator<< (std::ostream& s, const Integral_result& res);

  class Integral_data {

  private:
    data_pt* data;
    DP* abs_d2f_dx;  // second derivative * dx
    DP max_abs_d2f_dx;  //

  public:
    Integral_result integ_res;

  public:
    DP xmin;
    DP xmax;

  public:
    Integral_data (DP (*function_ref) (Vect_DP),
		   Vect_DP& args, int arg_to_integ_ref,
		   DP xmin_ref, DP xmax_ref);

    Integral_data (DP (*function_ref) (Vect_DP, I_table),
		   Vect_DP& args, int arg_to_integ_ref,
		   I_table Itable,
		   DP xmin_ref, DP xmax_ref);

    Integral_data (DP (*function_ref) (Vect_DP, Integral_table),
		   Vect_DP& args, int arg_to_integ_ref,
		   Integral_table Itable,
		   DP xmin_ref, DP xmax_ref);

    void Save (std::ofstream& outfile);

    void Improve_estimate (DP (*function) (Vect_DP),
			   Vect_DP& args, int arg_to_integ,
			   int Npts_max);

    void Improve_estimate (DP (*function) (Vect_DP, I_table),
			   Vect_DP& args, int arg_to_integ,
			   I_table Itable,
			   int Npts_max);

    void Improve_estimate (DP (*function) (Vect_DP, Integral_table),
			   Vect_DP& args, int arg_to_integ,
			   Integral_table Itable,
			   int Npts_max);

    ~Integral_data ();

  };

  Integral_result Integrate_optimal (DP (*function) (Vect_DP),
				     Vect_DP& args, int arg_to_integ,
				     DP xmin, DP xmax,
				     DP req_rel_prec, DP req_abs_prec,
				     int max_nr_pts);

  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable),
						 Vect_DP& args, int arg_to_integ,
						 I_table Itable,
						 DP xmin, DP xmax,
						 DP req_rel_prec, DP req_abs_prec,
						 int max_nr_pts);

  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, Integral_table Itable),
						 Vect_DP& args, int arg_to_integ,
						 Integral_table Itable,
						 DP xmin, DP xmax,
						 DP req_rel_prec, DP req_abs_prec,
						 int max_nr_pts);

  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable),
						 Vect_DP& args, int arg_to_integ,
						 I_table Itable,
						 DP xmin, DP xmax,
						 DP req_rel_prec, DP req_abs_prec,
						 int max_nr_pts,
						 std::ofstream& outfile);


  // ********************** Recursive version:  optimal, complex implementation ********************

  // NB:  function returns complex values but takes real arguments

  struct data_pt_CX {
    DP x;
    std::complex<DP> f;
    DP dx;
  };

  struct Integral_result_CX {
    std::complex<DP> integ_est;
    DP abs_prec;
    DP rel_prec;
    int n_vals;
  };

  class Integral_data_CX {

  private:
    data_pt_CX* data;
    DP* abs_d2f_dx;  // second derivative * dx
    DP max_abs_d2f_dx;  //

  public:
    Integral_result_CX integ_res;

  public:
    DP xmin;
    DP xmax;

  public:
    Integral_data_CX (std::complex<DP> (*function_ref) (Vect_DP),
		      Vect_DP& args, int arg_to_integ_ref,
		      DP xmin_ref, DP xmax_ref);

    void Save (std::ofstream& outfile);

    void Improve_estimate (std::complex<DP> (*function) (Vect_DP),
			   Vect_DP& args, int arg_to_integ,
			   int Npts_max);

    ~Integral_data_CX ();

  };

  Integral_result_CX Integrate_optimal (std::complex<DP> (*function) (Vect_DP),
					Vect_DP& args, int arg_to_integ,
					DP xmin, DP xmax,
					DP req_rel_prec, DP req_abs_prec,
					int max_nr_pts);

} // namespace ABACUS

#endif