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

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

Copyright (c) J.-S. Caux.

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

File:  ABACUS_ODSLF.h

Purpose:  Declares lattice spinless fermion classes and functions.

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

#ifndef ABACUS_ODSLF_H
#define ABACUS_ODSLF_H

#include "ABACUS.h"

namespace ABACUS {

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

  // Objects in class ODSLF_Base are a checked vector containing the number of rapidities of allowable types for a given state

  class ODSLF_Base {

  public:
    int Mdown;         // total number of down spins
    Vect<int> Nrap;    // Nrap[i] contains the number of rapidities of type i, i = 0, Nstrings - 1.
    int Nraptot;       // total number of strings in this state
    Vect<DP> Ix2_infty;  // Ix2_infty[i] contains the max of BAE function for the (half-)integer I[i], i = 0, Nstrings - 1.
    Vect<int> Ix2_max;    // Ix2_max[i] contains the integer part of 2*I_infty, with correct parity for base.
    long long int id;  // identification number

  public:
    ODSLF_Base ();
    ODSLF_Base (const ODSLF_Base& RefBase);  // copy constructor
    ODSLF_Base (const Heis_Chain& RefChain, int M);  // constructs configuration with all Mdown in one-string of +1 parity
    ODSLF_Base (const Heis_Chain& RefChain, const Vect<int>& Nrapidities);  // sets to Nrapidities vector, and checks consistency
    ODSLF_Base (const Heis_Chain& RefChain, long long int id_ref);
    inline int& operator[] (const int i);
    inline const int& operator[] (const int i) const;
    ODSLF_Base& operator= (const ODSLF_Base& RefBase);
    bool operator== (const ODSLF_Base& RefBase);
    bool operator!= (const ODSLF_Base& RefBase);

    void Compute_Ix2_limits(const Heis_Chain& RefChain);  // computes the Ix2_infty and Ix2_max

    void Scan_for_Possible_Types (Vect<long long int>& possible_type_id, int& nfound, int base_level, Vect<int>& Nexcitations);
    Vect<long long int> Possible_Types ();  // returns a vector of possible types

  };

  inline int& ODSLF_Base::operator[] (const int i)
  {
    return Nrap[i];
  }

  inline const int& ODSLF_Base::operator[] (const int i) const
  {
    return Nrap[i];
  }


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

  // Objects in class ODSLF_Ix2_Config carry all the I's of a given state

  class ODSLF_Ix2_Config {

  public:
    int Nstrings;
    Vect<int> Nrap;
    int Nraptot;

    int** Ix2;

  public:
    ODSLF_Ix2_Config ();
    ODSLF_Ix2_Config (const Heis_Chain& RefChain, int M);   // constructor, puts I's to ground state
    ODSLF_Ix2_Config (const Heis_Chain& RefChain, const ODSLF_Base& base);   // constructor, putting I's to lowest-energy config
    // consistent with Heis_Base configuration for chain RefChain
    ODSLF_Ix2_Config& operator= (const ODSLF_Ix2_Config& RefConfig);
    inline int* operator[] (const int i);
    inline const int* operator[] (const int i) const;
    ~ODSLF_Ix2_Config();
  };

  inline int* ODSLF_Ix2_Config::operator[] (const int i)
  {
    return Ix2[i];
  }

  inline const int* ODSLF_Ix2_Config::operator[] (const int i) const
  {
    return Ix2[i];
  }

  std::ostream& operator<< (std::ostream& s, const ODSLF_Ix2_Config& RefConfig);


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

  // Objects in class ODSLF_Lambda carry all rapidities of a state

  class ODSLF_Lambda {

  private:
    int Nstrings;
    Vect<int> Nrap;
    int Nraptot;
    DP** lambda;

  public:
    ODSLF_Lambda ();
    ODSLF_Lambda (const Heis_Chain& RefChain, int M);   // constructor, puts all lambda's to zero
    ODSLF_Lambda (const Heis_Chain& RefChain, const ODSLF_Base& base);   // constructor, putting I's to lowest-energy config
    // consistent with Heis_Base configuration for chain RefChain
    ODSLF_Lambda& operator= (const ODSLF_Lambda& RefConfig);
    inline DP* operator[] (const int i);
    inline const DP* operator[] (const int i) const;
    ~ODSLF_Lambda();

  };

  inline DP* ODSLF_Lambda::operator[] (const int i)
  {
    return lambda[i];
  }

  inline const DP* ODSLF_Lambda::operator[] (const int i) const
  {
    return lambda[i];
  }


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

  // Objects in class ODSLF_Ix2_Offsets carry Young tableau representations of the Ix2 configurations

  class ODSLF_Ix2_Offsets {

  public:
    ODSLF_Base base;
    Vect<Young_Tableau> Tableau;  // vector of pointers to tableaux at each level
    long long int type_id;
    long long int id;        // id number of offset
    long long int maxid;     // max id number allowable

  public:
    ODSLF_Ix2_Offsets ();
    ODSLF_Ix2_Offsets (const ODSLF_Ix2_Offsets& RefOffset);  // copy constructor
    ODSLF_Ix2_Offsets (const ODSLF_Base& RefBase, long long int req_type_id);
    ODSLF_Ix2_Offsets (const ODSLF_Base& RefBase, Vect<int> nparticles);  // sets all tableaux to empty ones, with nparticles[] at each level
    ODSLF_Ix2_Offsets& operator= (const ODSLF_Ix2_Offsets& RefOffset);
    bool operator<= (const ODSLF_Ix2_Offsets& RefOffsets);
    bool operator>= (const ODSLF_Ix2_Offsets& RefOffsets);

  public:
    void Set_to_id (long long int idnr);
    void Compute_id ();
    void Compute_type_id ();

  public:
    bool Add_Boxes_From_Lowest (int Nboxes, bool odd_sectors); // adds Nboxes in minimal energy config, all boxes in either even or odd sectors

  };

  inline long long int ODSLF_Ix2_Offsets_type_id (Vect<int>& nparticles)
  {
    long long int type_id_here = 0ULL;

    for (int i = 0; i < nparticles.size(); ++i)
      type_id_here += nparticles[i] * pow_ulli(10ULL, i);

    return(type_id_here);
  }


  //****************************************************************************
  // Objects in class ODSLF_Ix2_Offsets_List carry a vector of used Ix2_Offsets

  class ODSLF_Ix2_Offsets_List {

  public:
    int ndef;
    Vect<ODSLF_Ix2_Offsets> Offsets;

  public:
    ODSLF_Ix2_Offsets_List ();
    ODSLF_Ix2_Offsets& Return_Offsets (ODSLF_Base& RefBase, Vect<int> nparticles);  // returns the Ix2_Offsets corresponding to nparticles[]/base
    ODSLF_Ix2_Offsets& Return_Offsets (ODSLF_Base& RefBase, long long int req_type_id);
  };


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

  // Objects in class ODSLF_Bethe_State carry all information about an eigenstate

  // Derived classes include XXZ_Bethe_State, XXX_Bethe_State, XXZ_gpd_Bethe_State
  // These contain subclass-specific functions and data.

  class ODSLF_Bethe_State {

  public:
    Heis_Chain chain;
    ODSLF_Base base;
    ODSLF_Ix2_Offsets offsets;
    ODSLF_Ix2_Config Ix2;
    ODSLF_Lambda lambda;
    ODSLF_Lambda BE;            // Bethe equation for relevant rapidity, in the form BE = theta - (1/N)\sum ... - \pi I/N = 0
    DP diffsq;             // sum of squares of rapidity differences in last iteration
    int conv;              // convergence status
    int iter;              // number of iterations necessary for convergence
    int iter_Newton;              // number of iterations necessary for convergence (Newton method)
    DP E;                  // total energy
    int iK;                // K = 2.0*PI * iK/Nsites
    DP K;                  // total momentum
    DP lnnorm;             // ln of norm of reduced Gaudin matrix
    long long int base_id;
    long long int type_id;
    long long int id;
    long long int maxid;
    int nparticles;

  public:
    ODSLF_Bethe_State ();
    ODSLF_Bethe_State (const ODSLF_Bethe_State& RefState);  // copy constructor
    ODSLF_Bethe_State (const ODSLF_Bethe_State& RefState, long long int type_id_ref);  // new state with requested type_id
    ODSLF_Bethe_State (const Heis_Chain& RefChain, int M);  // constructor to ground-state configuration
    ODSLF_Bethe_State (const Heis_Chain& RefChain, const ODSLF_Base& base);  // constructor to lowest-energy config with base
    ODSLF_Bethe_State (const Heis_Chain& RefChain, long long int base_id_ref, long long int type_id_ref);
    virtual ~ODSLF_Bethe_State () {};

  public:
    int Charge () { return(base.Mdown); };
    void Set_Ix2_Offsets (const ODSLF_Ix2_Offsets& RefOffset);  // sets the Ix2 to given offsets
    void Set_to_id (long long int id_ref);
    void Set_to_id (long long int id_ref, ODSLF_Bethe_State& RefState);
    int Nparticles (); // counts the number of particles in state once Ix2 offsets set (so type_id is correctly set)
    bool Check_Symmetry ();  // checks whether the I's are symmetrically distributed
    void Compute_diffsq ();  // \sum BE[j][alpha]^2
    void Iterate_BAE ();     // Finds new set of lambda[j][alpha] from previous one by simple iteration
    void Iterate_BAE_Newton ();     // Finds new set of lambda[j][alpha] from previous one by a Newton step
    void Find_Rapidities (bool reset_rapidities);  // Finds the rapidities
    void Find_Rapidities_Twisted (bool reset_rapidities, DP twist);  // Finds the rapidities with twist added to RHS of logBE
    void BAE_smackdown (DP max_allowed);
    void Solve_BAE_smackdown (DP max_allowed, int maxruns);
    void Solve_BAE (int j, int alpha, DP req_prec, int itermax);
    void Solve_BAE_interp (DP interp_prec, int max_iter_interp);
    void Solve_BAE_straight_iter (DP interp_prec, int max_iter_interp);
    void Solve_BAE_Newton (DP Newton_prec, int max_iter_Newton);
    void Compute_lnnorm ();
    void Compute_Momentum ();
    void Compute_All (bool reset_rapidities);     // solves BAE, computes E, K and lnnorm
    bool Boost_Momentum (int iKboost);

    // Virtual functions, all defined in the derived classes

  public:
    virtual void Set_Free_lambdas() { ABACUSerror("ODSLF_Bethe_State::..."); }  // sets the rapidities to solutions of BAEs without scattering terms
    virtual bool Check_Admissibility(char option) { ABACUSerror("ODSLF_Bethe_State::..."); return(false); }
    // verifies that we don't have a symmetrical Ix2 config with a Ix2 == 0 for a string of even length >= 2.
    virtual void Compute_BE (int j, int alpha) { ABACUSerror("ODSLF_Bethe_State::..."); }
    virtual void Compute_BE () { ABACUSerror("ODSLF_Bethe_State::..."); }
    virtual DP Iterate_BAE(int i, int alpha) { ABACUSerror("ODSLF_Bethe_State::..."); return(0.0);}
    virtual bool Check_Rapidities() { ABACUSerror("ODSLF_Bethe_State::..."); return(false); }
    virtual void Compute_Energy () { ABACUSerror("ODSLF_Bethe_State::..."); }
    virtual void Build_Reduced_Gaudin_Matrix (SQMat<std::complex<DP> >& Gaudin_Red) { ABACUSerror("ODSLF_Bethe_State::..."); }
  };

  inline bool Force_Descent (char whichDSF, ODSLF_Bethe_State& ScanState, ODSLF_Bethe_State& RefState, int desc_type_required, int iKmod, DP Chem_Pot)
  {
    ABACUSerror("Need to implement Force_Descent properly for ODSLF.");

    bool force_descent = false;

    // Force descent if quantum nr distribution is symmetric:
    if (RefState.Check_Symmetry()) force_descent = true;

    return(force_descent);
  }

  std::ostream& operator<< (std::ostream& s, const ODSLF_Bethe_State& state);


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

  // Objects in class XXZ_Bethe_State carry all extra information pertaining to XXZ gapless

  class ODSLF_XXZ_Bethe_State : public ODSLF_Bethe_State {

  public:
    ODSLF_Lambda sinhlambda;
    ODSLF_Lambda coshlambda;
    ODSLF_Lambda tanhlambda;

  public:
    ODSLF_XXZ_Bethe_State ();
    ODSLF_XXZ_Bethe_State (const ODSLF_XXZ_Bethe_State& RefState);  // copy constructor
    ODSLF_XXZ_Bethe_State (const Heis_Chain& RefChain, int M);  // constructor to ground-state configuration
    ODSLF_XXZ_Bethe_State (const Heis_Chain& RefChain, const ODSLF_Base& base);  // constructor to lowest-energy config with base
    ODSLF_XXZ_Bethe_State (const Heis_Chain& RefChain, long long int base_id_ref, long long int type_id_ref);  // constructor to lowest-energy config with bas

  public:
    ODSLF_XXZ_Bethe_State& operator= (const ODSLF_XXZ_Bethe_State& RefState);

  public:
    void Set_Free_lambdas(); // sets the rapidities to solutions of BAEs without scattering terms
    void Compute_sinhlambda();
    void Compute_coshlambda();
    void Compute_tanhlambda();
    bool Check_Admissibility(char option);           // verifies that we don't have a symmetrical Ix2 config with a Ix2 == 0 for a string of even length >= 2.
    void Compute_BE (int j, int alpha);
    void Compute_BE ();
    DP Iterate_BAE(int i, int j);
    bool Check_Rapidities();  // checks that all rapidities are not nan
    void Compute_Energy ();
    void Build_Reduced_Gaudin_Matrix (SQMat<std::complex<DP> >& Gaudin_Red);

    // XXZ specific functions:
  public:

  };


  //****************************************************************************
  /*
  // Objects in class ODSLF_XXX_Bethe_State carry all extra information pertaining to XXX antiferromagnet

  class ODSLF_XXX_Bethe_State : public ODSLF_Bethe_State {

  public:
  ODSLF_XXX_Bethe_State ();
  ODSLF_XXX_Bethe_State (const ODSLF_XXX_Bethe_State& RefState);  // copy constructor
  ODSLF_XXX_Bethe_State (const Heis_Chain& RefChain, int M);  // constructor to ground-state configuration
  ODSLF_XXX_Bethe_State (const Heis_Chain& RefChain, const ODSLF__Base& base);  // constructor to lowest-energy config with base
  ODSLF_XXX_Bethe_State (const Heis_Chain& RefChain, long long int base_id_ref, long long int type_id_ref);  // constructor to lowest-energy config with base

  public:
  ODSLF_XXX_Bethe_State& operator= (const ODSLF_XXX_Bethe_State& RefState);

  public:
  void Set_Free_lambdas(); // sets the rapidities to solutions of BAEs without scattering terms
  bool Check_Admissibility(char option);           // verifies that we don't have a symmetrical Ix2 config with a Ix2 == 0 for a string of even length >= 2.
  void Compute_BE (int j, int alpha);
  void Compute_BE ();
  DP Iterate_BAE(int i, int j);
  bool Check_Rapidities();  // checks that all rapidities are not nan
  void Compute_Energy ();
  //void Compute_Momentum ();
  void Build_Reduced_Gaudin_Matrix (SQMat<std::complex<DP> >& Gaudin_Red);

  // XXX specific functions
  public:
  bool Check_Finite_rap ();

  };
  */
  //****************************************************************************
  /*
  // Objects in class ODSLF_XXZ_gpd_Bethe_State carry all extra information pertaining to XXZ gapped antiferromagnets

  class ODSLF_XXZ_gpd_Bethe_State : public ODSLF__Bethe_State {

  public:
  Lambda sinlambda;
  Lambda coslambda;
  Lambda tanlambda;

  public:
  ODSLF_XXZ_gpd_Bethe_State ();
  ODSLF_XXZ_gpd_Bethe_State (const ODSLF_XXZ_gpd_Bethe_State& RefState);  // copy constructor
  ODSLF_XXZ_gpd_Bethe_State (const Heis_Chain& RefChain, int M);  // constructor to ground-state configuration
  ODSLF_XXZ_gpd_Bethe_State (const Heis_Chain& RefChain, const ODSLF_Base& base);  // constructor to lowest-energy config with base
  ODSLF_XXZ_gpd_Bethe_State (const Heis_Chain& RefChain, long long int base_id_ref, long long int type_id_ref);  // constructor to lowest-energy config with base

  public:
  ODSLF_XXZ_gpd_Bethe_State& operator= (const ODSLF_XXZ_gpd_Bethe_State& RefState);

  public:
  void Set_Free_lambdas(); // sets the rapidities to solutions of BAEs without scattering terms
  void Compute_sinlambda();
  void Compute_coslambda();
  void Compute_tanlambda();
  int Weight();   // weight function for contributions cutoff
  bool Check_Admissibility(char option);           // verifies that we don't have a symmetrical Ix2 config with a Ix2 == 0 for a string of even length >= 2.
  void Compute_BE (int j, int alpha);
  void Compute_BE ();
  DP Iterate_BAE(int i, int j);
  void Iterate_BAE_Newton();
  bool Check_Rapidities();  // checks that all rapidities are not nan and are in interval ]-PI/2, PI/2]
  void Compute_Energy ();
  //void Compute_Momentum ();
  void Build_Reduced_Gaudin_Matrix (SQMat<std::complex<DP> >& Gaudin_Red);

  // XXZ_gpd specific functions
  public:

  };
  */
  //***********************************************

  // Function declarations
  /*
  // in M_vs_H.cc
  DP Ezero (DP Delta, int N, int M);
  DP H_vs_M (DP Delta, int N, int M);
  DP HZmin (DP Delta, int N, int M, Vect_DP& Ezero_ref);
  int M_vs_H (DP Delta, int N, DP HZ);

  DP X_avg (char xyorz, DP Delta, int N, int M);
  */
  DP Chemical_Potential (const ODSLF_Bethe_State& RefState);
  DP Sumrule_Factor (char whichDSF, ODSLF_Bethe_State& RefState, DP Chem_Pot, int iKmin, int iKmax);
  void Evaluate_F_Sumrule (std::string prefix, char whichDSF, const ODSLF_Bethe_State& RefState, DP Chem_Pot, int iKmin, int iKmax);

  std::complex<DP> ln_Sz_ME (ODSLF_XXZ_Bethe_State& A, ODSLF_XXZ_Bethe_State& B);
  std::complex<DP> ln_Smin_ME (ODSLF_XXZ_Bethe_State& A, ODSLF_XXZ_Bethe_State& B);

  DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, ODSLF_XXZ_Bethe_State& LeftState,
				     ODSLF_XXZ_Bethe_State& RefState, DP Chem_Pot, std::stringstream& DAT_outfile);


} // namespace ABACUS

#endif