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

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

Copyright (c) J.-S. Caux.

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

File:  ABACUS_Scan.h

Purpose:  Declares all classes and functions used in the
          ABACUS logic of scanning with threads.

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

#ifndef ABACUS_SCAN_H
#define ABACUS_SCAN_H

#include "ABACUS.h"

namespace ABACUS {

  const int MAX_STATE_LIST_SIZE = 10000;

  // Functions in src/UTILS/Data_File_Name.cc:
  void Data_File_Name (std::stringstream& name, char whichDSF, DP c_int, DP L, int N,
		       int iKmin, int iKmax, DP kBT, DP L2, std::string defaultname);
  void Data_File_Name (std::stringstream& name, char whichDSF, int iKmin, int iKmax, DP kBT,
		       LiebLin_Bethe_State& State, LiebLin_Bethe_State& RefScanState, std::string defaultname);
  void Data_File_Name (std::stringstream& name, char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
		       DP kBT, int N2, std::string defaultname);
  void Data_File_Name (std::stringstream& name, char whichDSF, int iKmin, int iKmax, DP kBT,
		       Heis_Bethe_State& State, Heis_Bethe_State& RefScanState, std::string defaultname);
  void ODSLF_Data_File_Name (std::stringstream& name, char whichDSF, DP Delta, int N, int M,
			     int iKmin, int iKmax, DP kBT, int N2, std::string defaultname);
  void Data_File_Name (std::stringstream& name, char whichDSF, int iKmin, int iKmax, DP kBT,
		       ODSLF_Bethe_State& State, ODSLF_Bethe_State& RefScanState, std::string defaultname);

  // Coding to convert ints to strings: for application in reduced labels
  //const int ABACUScodingsize = 64; // use a multiple of 2 to accelerate divisions in labeling.
  //const char ABACUScoding[] = {'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z', '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '!', '?'};

  // From ABACUS++T_8 onwards: forbid special characters as |, :, !, ? and all capital letters in labels.
  // This is due to the dumb capitalization-preserving but capitalization-insensitive HFS+ filesystem on Mac OS X.
  const int ABACUScodingsize = 32;
  const char ABACUScoding[] = {'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v'};
  const char LABELSEP = '_'; // was _
  const char TYPESEP = 'x';  // was |
  const char EXCSEP = 'y';   // was :
  const char INEXCSEP = 'z'; // was @

  struct State_Label_Data {

    Vect<int> type; // integer type labels of the types present
    Vect<int> M; // how many particles of each type
    Vect<int> nexc; // how many excitations as compared to the reference state used
    Vect<Vect<int> > Ix2old; // which Ix2 will be excited
    Vect<Vect<int> > Ix2exc; // which Ix2 the excitation has shifted to

    State_Label_Data (const Vect<int>& type_ref, const Vect<int>& M_ref,
		      const Vect<int>& nexc_ref, const Vect<Vect<int> >& Ix2old_ref, const Vect<Vect<int> >& Ix2exc_ref)
    {
      type = type_ref; M = M_ref; nexc = nexc_ref; Ix2old = Ix2old_ref; Ix2exc = Ix2exc_ref;
    }

  };

  std::string Extract_Base_Label (std::string label); // works for labels and complabels
  std::string Extract_nexc_Label (std::string label);

  // For compressed labels: conversions between integers and char/strings.
  std::string Convert_POSINT_to_STR (int int_to_convert);
  int Convert_CHAR_to_POSINT (char char_to_convert);
  int Convert_STR_to_POSINT (std::string str_to_convert);

  State_Label_Data Read_Base_Label (std::string label);
  State_Label_Data Read_State_Label (std::string label, const Vect<Vect<int> >& OriginIx2);
  State_Label_Data Read_State_Label (std::string label, const Vect<int>& OriginIx2); // if there is only one type
  std::string Return_State_Label (State_Label_Data data, const Vect<Vect<int> >& OriginIx2);
  std::string Return_State_Label (State_Label_Data data, const Vect<int>& OriginIx2); // if there is only one type
  std::string Return_State_Label (const Vect<Vect<int> >& ScanIx2, const Vect<Vect<int> >& OriginIx2);
  std::string Return_State_Label (const Vect<int>& ScanIx2, const Vect<int>& OriginIx2); // if there is only one type
  Vect<Vect<int> > Return_Ix2_from_Label (std::string label_ref, const Vect<Vect<int> >& OriginIx2);
  Vect<int> Return_Ix2_from_Label (std::string label_ref, const Vect<int>& OriginIx2); // specialization to Lieb-Liniger

  // Functions for descending states: in SCAN/Descendents.cc
  Vect<std::string> Descendent_States_with_iK_Stepped_Up (std::string ScanIx2_label, const LiebLin_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);
  Vect<std::string> Descendent_States_with_iK_Stepped_Down (std::string ScanIx2_label, const LiebLin_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);
  Vect<std::string> Descendent_States_with_iK_Preserved (std::string ScanIx2_label, const LiebLin_Bethe_State& OriginState, bool disperse_only_current_exc_up, bool preserve_nexc_up, bool disperse_only_current_exc_down, bool preserve_nexc_down);
  Vect<std::string> Descendent_States_with_iK_Stepped_Up (std::string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);
  Vect<std::string> Descendent_States_with_iK_Stepped_Down (std::string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);
  Vect<std::string> Descendent_States_with_iK_Preserved (std::string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc_up, bool preserve_nexc_up, bool disperse_only_current_exc_down, bool preserve_nexc_down);

  // For symmetric state scanning:
  Vect<std::string> Descendent_States_with_iK_Stepped_Up_rightIx2only
    (std::string ScanIx2_label, const LiebLin_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);
  Vect<std::string> Descendent_States_with_iK_Stepped_Down_rightIx2only
    (std::string ScanIx2_label, const LiebLin_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);
  Vect<std::string> Descendent_States_with_iK_Stepped_Up_rightIx2only
    (std::string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);
  Vect<std::string> Descendent_States_with_iK_Stepped_Down_rightIx2only
    (std::string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);

  // Functions to prepare and wrapup parallel scans:
  void Prepare_Parallel_Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
				      std::string defaultname, int paralevel, Vect<int> rank_lower_paralevels,
				      Vect<int> nr_processors_lower_paralevels, int nr_processors_at_newlevel);
  void Wrapup_Parallel_Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
				     std::string defaultname, int paralevel, Vect<int> rank_lower_paralevels,
				     Vect<int> nr_processors_lower_paralevels, int nr_processors_at_newlevel);

  void Prepare_Parallel_Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
				   int paralevel, Vect<int> rank_lower_paralevels,
				   Vect<int> nr_processors_lower_paralevels, int nr_processors_at_newlevel);
  void Wrapup_Parallel_Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
				  int paralevel, Vect<int> rank_lower_paralevels,
				  Vect<int> nr_processors_lower_paralevels, int nr_processors_at_newlevel);


  void Filter_RAW_File_for_iK (const char ff_file[], int iKneeded);
  void Sort_RAW_File (const char ffsq_file[], char optionchar);
  void Sort_RAW_File (const char ffsq_file[], char optionchar, char whichDSF);

  // Functions for data interpretation:
  DP Smoothen_RAW_into_SF (std::string prefix, int iKmin, int iKmax, int DiK, DP devmax,
			   DP ommin, DP ommax, int Nom, DP gwidth, DP normalization, DP denom_sum_K);
  DP Smoothen_RAW_into_SF (std::string prefix, Vect<std::string> rawfilename, Vect<DP> weight, int iKmin, int iKmax, int DiK, DP devmax,
			   DP ommin, DP ommax, int Nom, DP gwidth, DP normalization, DP denom_sum_K);
  void Write_K_File (DP Length, int iKmin, int iKmax);
  void Write_Omega_File (int Nout_omega, DP omegamin, DP omegamax);
  void Write_Time_File (int Nt, DP tmin, DP tmax);
  // Smoothen with gaussian width scaled with two-particle bandwidth
  DP Smoothen_RAW_into_SF_LiebLin_Scaled (std::string prefix, DP L, int N, int iKmin, int iKmax, int DiK, DP ommin, DP ommax, int Nom, DP width, DP normalization);


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

  struct Scan_Info {

    DP sumrule_obtained;
    DP Nfull; // dimensionality of (sub)Hilbert space considered
    long long int Ninadm;
    long long int Ndata;
    long long int Ndata_conv;
    long long int Ndata_conv0;
    double TT; // total computation time in seconds

  public:
    Scan_Info();  // constructor, puts everything to zero
    Scan_Info (DP sr, DP Nf, long long int Ni, long long int Nd, long long int Ndc, long long int Ndc0, double t);

    void Save (const char* outfile_Cstr);
    void Load (const char* infile_Cstr);

    inline Scan_Info& operator = (const Scan_Info& ref_info)
    {
      sumrule_obtained = ref_info.sumrule_obtained;
      Nfull = ref_info.Nfull;
      Ninadm = ref_info.Ninadm;
      Ndata = ref_info.Ndata;
      Ndata_conv = ref_info.Ndata_conv;
      Ndata_conv0 = ref_info.Ndata_conv0;
      TT = ref_info.TT;

      return(*this);
    }

    inline Scan_Info& operator+= (const Scan_Info& ref_info)
    {
      if (this != &ref_info) {
	sumrule_obtained += ref_info.sumrule_obtained;
	Nfull += ref_info.Nfull;
	Ninadm += ref_info.Ninadm;
	Ndata += ref_info.Ndata;
	Ndata_conv += ref_info.Ndata_conv;
	Ndata_conv0 += ref_info.Ndata_conv0;
	TT += ref_info.TT;
      }

      return(*this);
    }

    inline Scan_Info& operator-= (const Scan_Info& ref_info)
    {
      if (this != &ref_info) {
	sumrule_obtained -= ref_info.sumrule_obtained;
	Nfull -= ref_info.Nfull;
	Ninadm -= ref_info.Ninadm;
	Ndata -= ref_info.Ndata;
	Ndata_conv -= ref_info.Ndata_conv;
	Ndata_conv0 -= ref_info.Ndata_conv0;
	TT -= ref_info.TT;
      }

      return(*this);
    }

  };

  std::ostream& operator<< (std::ostream& s, const Scan_Info& info);

  // Functions in src/SCAN/General_Scan.cc:
  template<class Tstate>
    Scan_Info General_Scan (char whichDSF, int iKmin, int iKmax, int iKmod, DP kBT, Tstate& AveragingState, Tstate& SeedScanState,
			    std::string defaultScanStatename, int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors);

  Scan_Info Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
			  int Max_Secs, DP target_sumrule, bool refine,
			  int paralevel, Vect<int> rank, Vect<int> nr_processors);
  Scan_Info Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
			  int Max_Secs, DP target_sumrule, bool refine);
  Scan_Info Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, std::string defaultScanStatename,
			  int iKmin, int iKmax,
			  int Max_Secs, DP target_sumrule, bool refine,
			  int paralevel, Vect<int> rank, Vect<int> nr_processors);
  Scan_Info Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, std::string defaultname,
			  int iKmin, int iKmax,
			  int Max_Secs, DP target_sumrule, bool refine);
  Scan_Info Scan_LiebLin_Geometric_Quench (DP c_int, DP L_1, int type_id_1, long long int id_1, DP L_2, int N,
					   int iK_UL, int Max_Secs, DP target_sumrule, bool refine);

  Scan_Info Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
		       int Max_Secs, DP target_sumrule, bool refine,
		       int paralevel, Vect<int> rank, Vect<int> nr_processors);
  Scan_Info Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
		       int Max_Secs, DP target_sumrule, bool refine);
  Scan_Info Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, std::string defaultScanStatename,
		       int iKmin, int iKmax,
		       int Max_Secs, DP target_sumrule, bool refine,
		       int paralevel, Vect<int> rank, Vect<int> nr_processors);
  Scan_Info Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, std::string defaultScanStatename,
		       int iKmin, int iKmax,
		       int Max_Secs, DP target_sumrule, bool refine,
		       int paralevel, Vect<int> rank, Vect<int> nr_processors);
  Scan_Info Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, std::string defaultScanStatename,
		       int iKmin, int iKmax,
		       int Max_Secs, DP target_sumrule, bool refine,
		       int paralevel, Vect<int> rank, Vect<int> nr_processors);

  Scan_Info Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
			int Max_Secs, DP target_sumrule, bool refine, int rank, int nr_processors);
  Scan_Info Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax, int Max_Secs, bool refine);
  Scan_Info Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKneeded, int Max_Secs, bool refine);
  Scan_Info Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int Max_Secs, bool refine);


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

  // Functions in src/SCAN/Descendents.cc:
  Vect<std::string> Descendents (const LiebLin_Bethe_State& ScanState, const LiebLin_Bethe_State& OriginState, int type_required);
  Vect<std::string> Descendents (const Heis_Bethe_State& ScanState, const Heis_Bethe_State& OriginState, int type_required);


  struct Scan_Thread {

    std::string label;
    int type;

    Scan_Thread ();
    Scan_Thread (std::string label_ref, int type_ref) {
      label = label_ref;
      type = type_ref;
    }

    Scan_Thread& operator= (const Scan_Thread& RefThread);
  };


  struct Scan_Thread_Data {

    // By convention, a Scan_Thread_Data object handles a list of threads which are yet to be descended.
    // Improvement on Scan_Thread_Set used up to ABACUS++G_7, saving data to disk instead of holding it in memory.

    int nlists = 6400; // number of threads lists, fixed to this number by convention.
    DP logscale = (1.0/64) * log(2.0); // each separate list contains threads differing by a scale factor of 2^{1/64} \approx 1.01.

    std::string thrdir_name; // directory in which threads files are saved.

    Vect<int> nthreads_total;

    Vect<int> nthreads_on_disk;
    int lowest_il_with_nthreads_neq_0;

    // In-memory storage, for adding threads efficiently without constantly writing to disk
    // These objects are saved to disk when Next_Scan_Threads are called.
    Vect<int> dim;
    Vect<int> nthreads_in_memory;
    Vect<Vect<std::string> > label;
    Vect<Vect<int> > type; // which type of descendent is needed

    Vect<std::string> filename;

    Scan_Thread_Data ();
    Scan_Thread_Data (std::string thrdir_name_ref, bool refine);
    ~Scan_Thread_Data ();

    bool Increase_Memory_Size (int il, int nr_to_add);
    void Include_Thread (DP abs_data_value_ref, std::string label_ref, int type_ref);
    void Include_Thread (int il, std::string label_ref, int type_ref);

    Vect<Scan_Thread> Extract_Next_Scan_Threads (); // returns a vector of the threads that are next in line. By defn, all threads with index il == lowest_il_with_nthreads_neq_0. These are removed from the object.
    Vect<Scan_Thread> Extract_Next_Scan_Threads (int min_nr); // as above, but returns a minimum of min_nr threads.


    void Flush_to_Disk (int il);
    void Save ();
    void Load ();
  };



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

  // To populate a list of states for scanning:

  inline void Scan_for_Possible_Bases (const Vect<int> SeedNrap, const Vect<int> Str_L,
				       int Mdown_remaining, Vect<std::string>& possible_base_label, int& nfound, int nexc_max_used,
				       int base_level_to_scan, Vect<int>& Nrapidities)
  {
    if (Mdown_remaining < 0) { ABACUSerror("Scan_for_Possible_Bases:  shouldn't be here..."); }  // reached inconsistent point

    if (base_level_to_scan == 0) {
      if (Str_L[0] != 1) ABACUSerror("Str_L[0] != 1 in ABACUS_Scan.h Scan_for_Possible_Bases.");
      Nrapidities[0] = Mdown_remaining;

      // Set label:
      std::stringstream M0out;
      M0out << Nrapidities[0];
      possible_base_label[nfound] = M0out.str();
      for (int itype = 1; itype < Nrapidities.size(); ++itype)
	if (Nrapidities[itype] > 0) {
	  possible_base_label[nfound] += TYPESEP;
	  std::stringstream typeout;
	  typeout << itype;
	  possible_base_label[nfound] += typeout.str();
	  possible_base_label[nfound] += EXCSEP;
	  std::stringstream Mout;
	  Mout << Nrapidities[itype];
	  possible_base_label[nfound] += Mout.str();
	}
      nfound++;
    }

    else {
      // Preserve the number of strings at this level as compared to SeedState:
      Nrapidities[base_level_to_scan] = SeedNrap[base_level_to_scan];
      if (Mdown_remaining - Str_L[base_level_to_scan] * Nrapidities[base_level_to_scan] >= 0)
	Scan_for_Possible_Bases (SeedNrap, Str_L, Mdown_remaining - Str_L[base_level_to_scan] * Nrapidities[base_level_to_scan],
				 possible_base_label, nfound, nexc_max_used, base_level_to_scan - 1, Nrapidities);

      // Reduce number of strings at this level as compared to SeedState:
      for (int i = 1; i <= ABACUS::min(SeedNrap[base_level_to_scan], nexc_max_used/Str_L[base_level_to_scan]); ++i) {
	Nrapidities[base_level_to_scan] = SeedNrap[base_level_to_scan] - i;
	if (Mdown_remaining - Str_L[base_level_to_scan] * Nrapidities[base_level_to_scan] >= 0)
	  Scan_for_Possible_Bases (SeedNrap, Str_L, Mdown_remaining - Str_L[base_level_to_scan] * Nrapidities[base_level_to_scan],
				   possible_base_label, nfound, nexc_max_used - i*Str_L[base_level_to_scan], base_level_to_scan - 1, Nrapidities);
      }
      // Increase the number of strings at this level as compared to SeedState:
      for (int i = 1; i <= ABACUS::min(Mdown_remaining/Str_L[base_level_to_scan], nexc_max_used/Str_L[base_level_to_scan]); ++i) {
	Nrapidities[base_level_to_scan] = SeedNrap[base_level_to_scan] + i;
	if (Mdown_remaining - Str_L[base_level_to_scan] * Nrapidities[base_level_to_scan] >= 0)
	  Scan_for_Possible_Bases (SeedNrap, Str_L, Mdown_remaining - Str_L[base_level_to_scan] * Nrapidities[base_level_to_scan],
				   possible_base_label, nfound, nexc_max_used - i*Str_L[base_level_to_scan], base_level_to_scan - 1, Nrapidities);
      }
    }
  }


  inline Vect<std::string> Possible_Bases (const Vect<int> SeedNrap, const Vect<int> Str_L, int Mdown)//const Heis_Bethe_State& SeedState)
  {
    int nexc_max_used = NEXC_MAX_HEIS;

    Vect<std::string> possible_base_label (1000);
    int nfound = 0;
    Vect<int> Nrapidities = SeedNrap;
    int Mdown_remaining = Mdown;

    Scan_for_Possible_Bases (SeedNrap, Str_L, Mdown_remaining, possible_base_label, nfound, nexc_max_used, SeedNrap.size() - 1, Nrapidities);

    // Copy results into a clean vector:
    Vect<std::string> possible_base_label_found (nfound);
    for (int i = 0; i < nfound; ++i) possible_base_label_found[i] = possible_base_label[i];

    return(possible_base_label_found);
  }



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

  template<class Tstate>
    class Scan_State_List {

  public:
    int ndef;
    Vect<Tstate> State;
    Vect<std::string> base_label;
    Vect<Scan_Info> info;  // info for base and type of State[n]
    Vect<bool> flag_for_scan; // set to true, next round of scanning will use this base/type
    Vect<bool> scan_attempted; // whether this has already been attempted

  public:
    inline Scan_State_List (char whichDSF, const Tstate& SeedScanState);

  public:
    inline Tstate& Return_State (std::string base_label_ref);  // returns a state corresponding to same base and type
    inline void Populate_List (char whichDSF, const Tstate& SeedScanState);  // creates all types of states containing up to nexc_max excitations
    inline void Include_Info (Scan_Info& info_to_add, std::string base_label_ref);
    inline void Raise_Scanning_Flags (DP threshold); // checks whether base/type should be scanned based on simpler base/type combinations

    inline void Order_in_SRC ();
    inline void Save_Info (const char* sumfile_Cstr);
    inline void Load_Info (const char* sumfile_Cstr);
  };

  // Do the explicit class specializations:

  template<>
    inline Scan_State_List<LiebLin_Bethe_State>::Scan_State_List (char whichDSF, const LiebLin_Bethe_State& SeedScanState)
    : ndef(0), State(Vect<LiebLin_Bethe_State>(MAX_STATE_LIST_SIZE)), base_label(Vect<std::string>(MAX_STATE_LIST_SIZE)),
    info(Vect<Scan_Info>(MAX_STATE_LIST_SIZE)), flag_for_scan(Vect<bool>(false, MAX_STATE_LIST_SIZE)),
    scan_attempted(Vect<bool>(false, MAX_STATE_LIST_SIZE))
      {
	State[0] = SeedScanState;
      }

  template<>
    inline Scan_State_List<XXZ_Bethe_State>::Scan_State_List (char whichDSF, const XXZ_Bethe_State& SeedScanState)
    : ndef(0), State(Vect<XXZ_Bethe_State>(MAX_STATE_LIST_SIZE)), base_label(Vect<std::string>(MAX_STATE_LIST_SIZE)),
    info(Vect<Scan_Info>(MAX_STATE_LIST_SIZE)), flag_for_scan(Vect<bool>(false, MAX_STATE_LIST_SIZE)),
    scan_attempted(Vect<bool>(false, MAX_STATE_LIST_SIZE))
      {
	State[0] = SeedScanState;
      }

  template<>
    inline Scan_State_List<XXX_Bethe_State>::Scan_State_List (char whichDSF, const XXX_Bethe_State& SeedScanState)
    : ndef(0), State(Vect<XXX_Bethe_State>(MAX_STATE_LIST_SIZE)), base_label(Vect<std::string>(MAX_STATE_LIST_SIZE)),
    info(Vect<Scan_Info>(MAX_STATE_LIST_SIZE)), flag_for_scan(Vect<bool>(false, MAX_STATE_LIST_SIZE)),
    scan_attempted(Vect<bool>(false, MAX_STATE_LIST_SIZE))
      {
	State[0] = SeedScanState;
      }

  template<>
    inline Scan_State_List<XXZ_gpd_Bethe_State>::Scan_State_List (char whichDSF, const XXZ_gpd_Bethe_State& SeedScanState)
    : ndef(0), State(Vect<XXZ_gpd_Bethe_State>(MAX_STATE_LIST_SIZE)), base_label(Vect<std::string>(MAX_STATE_LIST_SIZE)),
    info(Vect<Scan_Info>(MAX_STATE_LIST_SIZE)), flag_for_scan(Vect<bool>(false, MAX_STATE_LIST_SIZE)),
    scan_attempted(Vect<bool>(false, MAX_STATE_LIST_SIZE))
      {
	State[0] = SeedScanState;
      }

  /* IN_DEVELOPMENT
     template<>
     inline Scan_State_List<ODSLF_XXZ_Bethe_State>::Scan_State_List (char whichDSF, const ODSLF_XXZ_Bethe_State& RefState)
     : ndef(0), State(Vect<ODSLF_XXZ_Bethe_State>(MAX_STATE_LIST_SIZE)), base_label(Vect<std::string>(MAX_STATE_LIST_SIZE)),
     info(Vect<Scan_Info>(MAX_STATE_LIST_SIZE)), flag_for_scan(Vect<bool>(false, MAX_STATE_LIST_SIZE)),
     scan_attempted(Vect<bool>(false, MAX_STATE_LIST_SIZE))
     {
     if (whichDSF == 'Z' || whichDSF == 'z') State[0] = ODSLF_XXZ_Bethe_State(RefState.chain, RefState.base.Mdown);
     else if (whichDSF == 'm') State[0] = ODSLF_XXZ_Bethe_State(RefState.chain, RefState.base.Mdown - 1);
     else if (whichDSF == 'p') State[0] = ODSLF_XXZ_Bethe_State(RefState.chain, RefState.base.Mdown + 1);
     else ABACUSerror("Unknown whichDSF in Scan_State_List<ODSLF_XXZ... constructor.");
     }
  */

  template<>
    inline LiebLin_Bethe_State& Scan_State_List<LiebLin_Bethe_State>::Return_State (std::string base_label_ref)
    {
      int n = 0;
      while (n < ndef && base_label_ref.compare(base_label[n]) != 0) n++;

      if (n == ndef) {
	State[n] = State[0];
	base_label[n] = base_label_ref;
	info[n].Nfull = 1LL; // Nfull not definable for LiebLin
	ndef++;
      }

      return(State[n]);
    }

  template<>
    inline XXZ_Bethe_State& Scan_State_List<XXZ_Bethe_State>::Return_State (std::string base_label_ref)
    {
      int n = 0;
      while (n < ndef && base_label_ref.compare(base_label[n]) != 0) n++;

      if (n == ndef) {
	Heis_Base checkbase (State[0].chain, base_label_ref);
	State[n] = XXZ_Bethe_State (State[0].chain, checkbase);
	info[n].Nfull = checkbase.dimH;
	ndef++;
      }

      return(State[n]);
    }

  template<>
    inline XXX_Bethe_State& Scan_State_List<XXX_Bethe_State>::Return_State (std::string base_label_ref)
    {
      int n = 0;
      while (n < ndef && base_label_ref.compare(base_label[n]) != 0) n++;

      if (n == ndef) {
	Heis_Base checkbase (State[0].chain, base_label_ref);
	State[n] = XXX_Bethe_State (State[0].chain, checkbase);
	info[n].Nfull = checkbase.dimH;
	ndef++;
      }

      return(State[n]);
    }

  template<>
    inline XXZ_gpd_Bethe_State& Scan_State_List<XXZ_gpd_Bethe_State>::Return_State (std::string base_label_ref)
    {
      int n = 0;
      while (n < ndef && base_label_ref.compare(base_label[n]) != 0) n++;

      if (n == ndef) {
	Heis_Base checkbase (State[0].chain, base_label_ref);
	State[n] = XXZ_gpd_Bethe_State (State[0].chain, checkbase);
	info[n].Nfull = checkbase.dimH;
	ndef++;
      }

      return(State[n]);
    }

  /* IN DEVELOPMENT
     template<>
     inline ODSLF_XXZ_Bethe_State& Scan_State_List<ODSLF_XXZ_Bethe_State>::Return_State (long long int base_id_ref, long long int type_id_ref)
     {
     int n = 0;
     while (n < ndef && !(base_id_ref == State[n].base_id && type_id_ref == State[n].type_id)) n++;

     if (n == ndef) {
     State[n] = ODSLF_XXZ_Bethe_State (State[0].chain, base_id_ref, type_id_ref);
     info[n].Nfull = State[n].maxid + 1LL;
     ndef++;
     }

     return(State[n]);
     }
  */

  template<>
    inline void Scan_State_List<LiebLin_Bethe_State>::Populate_List (char whichDSF, const LiebLin_Bethe_State& SeedScanState)
    {
      // For LiebLin_Bethe_State: only one base is used, so there is only one state here.

      if (ndef != 0) ABACUSerror("Please only populate a virgin Scan_State_List.");

      std::stringstream baselabel;
      baselabel << State[0].N;
      base_label[0] = baselabel.str();
      std::stringstream label0;
      label0 << State[0].N << LABELSEP << ABACUScoding[0] << LABELSEP;//"_0_";
      State[0].Set_to_Label (label0.str(), SeedScanState.Ix2);
      info[0].Nfull = 1LL; // Nfull not definable for LiebLin
      ndef = 1;
    }

  template<>
    inline void Scan_State_List<XXZ_Bethe_State>::Populate_List (char whichDSF, const XXZ_Bethe_State& SeedScanState)
    {
      // creates all types of states containing up to nexc_max excitations

      if (ndef != 0) ABACUSerror("Please only populate a virgin Scan_State_List.");

      // We assume that SeedScanState has quantum numbers which are set according to the relevant AveragingState.

      // This function creates a list of states with other bases in the vicinity of that of SeedScanState,
      // matching the quantum numbers as closely as possible.

      Vect<int> Str_L(SeedScanState.chain.Nstrings);
      for (int i = 0; i < SeedScanState.chain.Nstrings; ++i) Str_L[i] = SeedScanState.chain.Str_L[i];

      // First of all, we create a list of the possible bases themselves.
      Vect<std::string> bases_label = Possible_Bases (SeedScanState.base.Nrap, Str_L, SeedScanState.base.Mdown);  // returns a vector of possible bases

      for (int ib = 0; ib < bases_label.size(); ++ib) {

	Heis_Base checkbase (State[0].chain, bases_label[ib]);

	State[ndef] = XXZ_Bethe_State (State[0].chain, checkbase);
	State[ndef].Set_to_Closest_Matching_Ix2_fixed_Base (SeedScanState);
	State[ndef].Set_Label_from_Ix2 (State[ndef].Ix2); // sets to trivial label for this base
	base_label[ndef] = bases_label[ib];
	info[ndef].Nfull = State[ndef].base.dimH;
	ndef++;
	if (ndef >= MAX_STATE_LIST_SIZE) ABACUSerror("Increase number of elements in ScanStateList.");
      }
    }

  template<>
    inline void Scan_State_List<XXX_Bethe_State>::Populate_List (char whichDSF, const XXX_Bethe_State& SeedScanState)
    {
      // creates all types of states containing up to nexc_max excitations

      if (ndef != 0) ABACUSerror("Please only populate a virgin Scan_State_List.");

      // We assume that SeedScanState has quantum numbers which are set according to the relevant AveragingState.

      // This function creates a list of states with other bases in the vicinity of that of SeedScanState,
      // matching the quantum numbers as closely as possible.

      Vect<int> Str_L(SeedScanState.chain.Nstrings);
      for (int i = 0; i < SeedScanState.chain.Nstrings; ++i) Str_L[i] = SeedScanState.chain.Str_L[i];

      // To take infinite rapidities into account, we use intermediate states with up to 2 less finite rapidities (1 for Szz, 2 for Spm)
      int nrinfrapmax = 0;
      if (whichDSF == 'z') nrinfrapmax = 1;
      else if (whichDSF == 'p') nrinfrapmax = ABACUS::min(2, SeedScanState.base.Mdown);

      Vect<int> Nrapmod = SeedScanState.base.Nrap;

      for (int nrinfrap = 0; nrinfrap <= nrinfrapmax; ++nrinfrap) {

	Nrapmod[0] = SeedScanState.base.Nrap[0] - nrinfrap;
	if (Nrapmod[0] < 0) ABACUSerror("Putting too many rapidities at infinity in ABACUS_Scan.h: Possible_Bases.");

	Vect<std::string> bases_label = Possible_Bases (Nrapmod, Str_L, SeedScanState.base.Mdown-nrinfrap);  // returns a vector of possible bases

	for (int ib = 0; ib < bases_label.size(); ++ib) {

	  Heis_Base checkbase (State[0].chain, bases_label[ib]);

	  State[ndef] = XXX_Bethe_State (State[0].chain, checkbase);
	  State[ndef].Set_to_Closest_Matching_Ix2_fixed_Base (SeedScanState);
	  base_label[ndef] = bases_label[ib];
	  info[ndef].Nfull = State[ndef].base.dimH;
	  ndef++;
	  if (ndef >= MAX_STATE_LIST_SIZE) ABACUSerror("Increase number of elements in ScanStateList.");
	}
      } // for nrinfrap
    }

  template<>
    inline void Scan_State_List<XXZ_gpd_Bethe_State>::Populate_List (char whichDSF, const XXZ_gpd_Bethe_State& SeedScanState)
    {
      // creates all types of states containing up to nexc_max excitations

      if (ndef != 0) ABACUSerror("Please only populate a virgin Scan_State_List.");

      // We assume that SeedScanState has quantum numbers which are set according to the relevant AveragingState.

      // This function creates a list of states with other bases in the vicinity of that of SeedScanState,
      // matching the quantum numbers as closely as possible.

      Vect<int> Str_L(SeedScanState.chain.Nstrings);
      for (int i = 0; i < SeedScanState.chain.Nstrings; ++i) Str_L[i] = SeedScanState.chain.Str_L[i];

      // First of all, we create a list of the possible bases themselves.
      Vect<std::string> bases_label = Possible_Bases (SeedScanState.base.Nrap, Str_L, SeedScanState.base.Mdown);  // returns a vector of possible bases

      for (int ib = 0; ib < bases_label.size(); ++ib) {

	Heis_Base checkbase (State[0].chain, bases_label[ib]);

	State[ndef] = XXZ_gpd_Bethe_State (State[0].chain, checkbase);
	State[ndef].Set_to_Closest_Matching_Ix2_fixed_Base (SeedScanState);
	base_label[ndef] = bases_label[ib];
	info[ndef].Nfull = State[ndef].base.dimH;
	ndef++;
	if (ndef >= MAX_STATE_LIST_SIZE) ABACUSerror("Increase number of elements in ScanStateList.");
      }
    }


  /* IN DEVELOPMENT
     template<>
     inline void Scan_State_List<ODSLF_XXZ_Bethe_State>::Populate_List ()
     {
     // creates all types of states containing up to nexc_max excitations

     if (ndef != 0) ABACUSerror("Please only populate a virgin Scan_State_List.");

     //std::cout << "In Populate_List: " << State[0] << std::endl;

     Vect<long long int> bases_id = State[0].chain.Possible_Bases (State[0].base.Mdown);  // returns a vector of possible bases

     //std::cout << "Mdown = " << State[0].base.Mdown << "\tPossible bases size: " << bases_id.size() << "\tPossible bases: " << bases_id << std::endl;

     for (int ib = 0; ib < bases_id.size(); ++ib) {
     ODSLF_Base checkbase (State[0].chain, bases_id[ib]);

     Vect<long long int> types_id = checkbase.Possible_Types ();  // returns a vector of possible types

     //std::cout << "For base_id " << bases_id[ib] << "\t found types " << types_id << std::endl;

     for (int it = 0; it < types_id.size(); ++it) {

     if (bases_id[ib] < 1000000) { // FUDGE:  consider only one-strings
     //std::cout << "Populate list:  constructing state: " << bases_id[ib] << "\t" << types_id[it] << std::endl;
     State[ndef] = ODSLF_XXZ_Bethe_State (State[0].chain, bases_id[ib], types_id[it]);
     //std::cout << "Populate list:  before setting id: " << std::endl << State[ndef] << std::endl;
     State[ndef].Set_to_id(0LL);
     //std::cout << "Populate list:  after setting id:  "  << std::endl << State[ndef] << std::endl;
     info[ndef].Nfull = State[ndef].maxid + 1LL;
     ndef++;
     if (ndef >= MAX_STATE_LIST_SIZE) ABACUSerror("Increase number of elements in ScanStateList.");
     }
     }

     }
     }

  */



  template<class Tstate>
    inline void Scan_State_List<Tstate>::Include_Info (Scan_Info& info_to_add, std::string base_label_ref)
    {

      int n = 0;
      while (n < ndef && base_label_ref.compare(base_label[n]) != 0) n++;

      if (n == ndef) {
	std::cout << "ndef = " << ndef << std::endl;
	for (int i = 0; i < ndef; ++i) std::cout << base_label[i] << "\t";
	std::cout << std::endl;
	std::cout << "base_label_ref " << base_label_ref << std::endl;
	ABACUSerror("Did not find base_label_ref in Scan_State_List::Include_Info.");
      }

      info[n] += info_to_add;
      return;
    }


  template<class Tstate>
    inline void Scan_State_List<Tstate>::Raise_Scanning_Flags (DP threshold)
    {
      flag_for_scan = true;
    }

  template<class Tstate>
    inline void Scan_State_List<Tstate>::Order_in_SRC ()
    {
      if (ndef > 0) {

	Vect_INT index(ndef);
	for (int i = 0; i < ndef; ++i) index[i] = i;

	Vect<DP> sr (ndef);
	for (int i = 0; i < ndef; ++i) sr[i] = info[i].sumrule_obtained;

	sr.QuickSort(index, 0, ndef - 1);

	Vect<Tstate> State_ordered(ndef);
	Vect<std::string> base_label_ordered(ndef);
	Vect<Scan_Info> info_ordered(ndef);
	Vect<bool> flag_for_scan_ordered(ndef);

	// Put data in proper order
	for (int i = 0; i < ndef; ++i) {
	  State_ordered[i] = State[index[ndef - 1 - i] ];
	  base_label_ordered[i] = base_label[index[ndef - 1 - i] ];
	  info_ordered[i] = info[index[ndef - 1 - i] ];
	  flag_for_scan_ordered[i] = flag_for_scan[index[ndef - 1 - i] ];
	}

	// Put back in *this object:
	for (int i = 0; i < ndef; ++i) {
	  State[i] = State_ordered[i];
	  base_label[i] = base_label_ordered[i];
	  info[i] = info_ordered[i];
	  flag_for_scan[i] = flag_for_scan_ordered[i];
	} // The rest are all simply 0.
      }
    }

  template<class Tstate>
    inline void Scan_State_List<Tstate>::Save_Info (const char* sumfile_Cstr)
    {
      std::ofstream outfile;

      outfile.open(sumfile_Cstr);
      if (outfile.fail()) ABACUSerror("Could not open outfile... ");

      outfile.setf(std::ios::fixed);
      outfile.precision(16);

      outfile << std::setw(20) << "base" << std::setw(25) << "sumrule_obtained" << std::setw(25) << "Nfull" << std::setw(10) << "Ninadm" << std::setw(10) << "Ndata" << std::setw(10) << "conv" << std::setw(10) << "conv0" << std::setw(10) << "TT.";

      for (int i = 0; i < ndef; ++i)
	if (info[i].Nfull > 0.0) {
	  int TT_hr = int(info[i].TT/3600);
	  int TT_min = int((info[i].TT - 3600.0*TT_hr)/60);
	  outfile << std::endl << std::setw(20) << base_label[i] << std::setw(25);
	  if (info[i].sumrule_obtained < 1.0) outfile << std::fixed;
	  else outfile << std::scientific;
	  outfile << std::setprecision(16) << info[i].sumrule_obtained;
	  if (info[i].Nfull < 1.0e+10) outfile << std::setw(25) << std::fixed << std::setprecision(0) << info[i].Nfull;
	  else outfile << std::setw(25) << std::scientific << std::setprecision(16) << info[i].Nfull;
	  outfile << std::setw(10) << info[i].Ninadm << std::setw(10) << info[i].Ndata << std::setw(10) << info[i].Ndata_conv << std::setw(10) << info[i].Ndata_conv0 << std::setw(10) << TT_hr << " h " << TT_min << " m " << std::fixed << std::showpoint << std::setprecision(3) << info[i].TT - 3600.0*TT_hr - 60.0*TT_min << " s";
	}
      outfile.close();
    }

  template<class Tstate>
    inline void Scan_State_List<Tstate>::Load_Info (const char* sumfile_Cstr)
    {
      std::ifstream infile;
      infile.open(sumfile_Cstr);
      if(infile.fail()) {
	std::cout << std::endl << sumfile_Cstr << std::endl;
	ABACUSerror("Could not open input file in Scan_State_List::Load_Info.");
      }

      // Load first line, containing informative text:
      char junk[256];
      infile.getline(junk, 256);

      // Now load the previous info's:
      std::string base_label_ref;
      DP sr_ref;
      DP Nfull_ref;
      long long int Ninadm_ref, Ndata_ref, conv_ref, conv0_ref;
      DP TT_ref;
      int TT_hr, TT_min;
      DP TT_sec;
      char a;

      while (infile.peek() != EOF) {
	infile >> base_label_ref >> sr_ref >> Nfull_ref >> Ninadm_ref >> Ndata_ref >> conv_ref >> conv0_ref >> TT_hr >> a >> TT_min >> a >> TT_sec >> a;
	TT_ref = 3600.0 * TT_hr + 60.0* TT_min + TT_sec;
	Scan_Info info_ref (sr_ref, Nfull_ref, Ninadm_ref, Ndata_ref, conv_ref, conv0_ref, TT_ref);
	(*this).Include_Info (info_ref, base_label_ref);
      }

      infile.close();

      return;
    }


} // namespace ABACUS

#endif