commit a6a4a1bdb450c1835827fbb0390e9660aff7e3ee Author: Jean-Sébastien Caux Date: Tue Sep 22 14:20:09 2026 +0200 Initiate diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..409dfda --- /dev/null +++ b/.gitignore @@ -0,0 +1,18 @@ + +*~ +benchmarks +bin +dep +docs +lib +obj +tests +*html + +gcm.cache +ltximg + +abacus.makeinfo + +# Apple annoyances +.DS_Store \ No newline at end of file diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..fc74804 --- /dev/null +++ b/Makefile @@ -0,0 +1,186 @@ +################################################################### +# +# This software is part of Jean-Sébastien Caux's Abacus toolsuite. +# +# Copyright © Jean-Sébastien Caux. +# +################################################################### + + + +VERSION = abacus-2.0.0 + +SRC_DIR = src +SRC_EXT = cc +DEP_DIR = dep +OBJ_DIR = obj +LIB_DIR = lib +BIN_DIR = bin + +$(shell mkdir -p $(DEP_DIR) >/dev/null) +$(shell mkdir -p $(OBJ_DIR) >/dev/null) +$(shell mkdir -p $(LIB_DIR) >/dev/null) +$(shell mkdir -p $(BIN_DIR) >/dev/null) + + + +###################### +# Compilation with g++ +###################### +# NOTE: g++-15 is required +# Uncomment the lines below for g++ compilation +CXX = g++ #g++-15 +OPTIMIZATIONS = -O3 +CXXFLAGS = -std=c++23 -fmodules -Wall -Wextra -Wconversion -pedantic-errors $(OPTIMIZATIONS) #-g -pg +COMPILE = $(CXX) $(CXXFLAGS) +COMPILE_EXECS = $(COMPILE) + +######################## +# Compilation with clang +######################## +# NOTE 2025-03-14: -O3 compilation fails on Linux with clang++ 18.1.3 +# ******* NOTE FOR FUTURE: once -O3 compilation problem is solved, +# ******* also add -O3 to clang-std make option at bottom of this Makefile +# Uncomment the lines below for clang++ compilation +# clang 18: +# CXX = clang++ +# CXXFLAGS = -std=c++23 -stdlib=libc++ --gcc-toolchain=/usr/local/gcc-15-20250302 -fprebuilt-module-path=obj/ #-O3 +# COMPILE = $(CXX) $(CXXFLAGS) -x c++-module -fmodule-output +# COMPILE_EXECS = $(CXX) $(CXXFLAGS) obj/std.pcm +# clang 21 +# CXX = clang++-21 +# CXXFLAGS = -std=c++23 -stdlib=libc++ --gcc-toolchain=/usr/local/gcc-15.2.0 -fprebuilt-module-path=obj/ -O3 +# COMPILE = $(CXX) $(CXXFLAGS) -x c++-module -fmodule-output +# COMPILE_EXECS = $(CXX) $(CXXFLAGS) + + +CXXINFO = $(shell $(CXX) --version) + + +############################# +# Dependencies auto-detection +# 2025-03: not working yet +DEPFLAGS = -MT $@ -MMD -MP -MF $(DEP_DIR)/$*.Td +# COMPILE = $(CXX) $(CXXFLAGS) $(DEPFLAGS) -c +POSTCOMPILE = mv -f $(DEP_DIR)/$*.Td $(DEP_DIR)/$*.d && touch $@ + + +# All sources +SOURCES_ALL = $(shell find $(SRC_DIR) -name '*.$(SRC_EXT)') + +# Sources for executables +SOURCES_EXECS = $(shell find $(SRC_DIR)/execs -name '*.$(SRC_EXT)') + +# Library sources +SOURCES = $(filter-out $(SOURCES_EXECS), $(SOURCES_ALL)) + +# Dependencies +DEPS = $(patsubst %.$(SRC_EXT), $(DEP_DIR)/%.d, $(notdir $(SOURCES))) + +# Objects to go into library +OBJECTS = $(patsubst %.$(SRC_EXT), $(OBJ_DIR)/%.o, $(notdir $(SOURCES))) + +# Executables +EXECS = $(patsubst %.$(SRC_EXT), $(BIN_DIR)/%, $(notdir $(SOURCES_EXECS))) +#EXECS = $(BIN_DIR)/debug + + +# Default target: executables +all: $(EXECS) + echo "Abacus version $(VERSION) copyright © Jean-Sébastien Caux\ncompiled on $(shell date -Imin) with\n$(CXXINFO)\nusing flags:\n$(CXXFLAGS)" | tee abacus.makeinfo > /dev/null + + + +# Create the library +$(LIB_DIR)/lib$(VERSION).a : $(OBJECTS) + ar -cr $(LIB_DIR)/lib$(VERSION).a $(OBJECTS) + + +# Build executables +$(EXECS): $(BIN_DIR)/%: $(SRC_DIR)/execs/%.$(SRC_EXT) $(LIB_DIR)/lib$(VERSION).a + $(COMPILE_EXECS) $< -o $@ -L$(LIB_DIR) -l$(VERSION) + + + +# Autodetection of dependencies fails using DEPFLAGS above. +# Temporary backup strategy: dependencies set by hand + +$(OBJ_DIR)/conveniences.o: $(SRC_DIR)/generic/conveniences.cc + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/timer.o: $(SRC_DIR)/util/timer.cc + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/labels.o: $(SRC_DIR)/generic/labels.cc $(OBJ_DIR)/conveniences.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/spaces.o: $(SRC_DIR)/generic/spaces.cc $(OBJ_DIR)/conveniences.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/model.o: $(SRC_DIR)/generic/model.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/spaces.o $(OBJ_DIR)/labels.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/model_LiebLiniger.o: $(SRC_DIR)/LiebLiniger/repulsive/model_LiebLiniger.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/spaces.o $(OBJ_DIR)/model.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/quantumnumbers.o: $(SRC_DIR)/generic/quantumnumbers.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/labels.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/plex.o: $(SRC_DIR)/generic/plex.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/labels.o $(OBJ_DIR)/quantumnumbers.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/matrix.o: $(SRC_DIR)/math/matrix.cc $(OBJ_DIR)/conveniences.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/state.o: $(SRC_DIR)/generic/state.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/timer.o $(OBJ_DIR)/labels.o $(OBJ_DIR)/spaces.o $(OBJ_DIR)/model.o $(OBJ_DIR)/plex.o $(OBJ_DIR)/matrix.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/state_LiebLiniger.o: $(SRC_DIR)/LiebLiniger/repulsive/state_LiebLiniger.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/labels.o $(OBJ_DIR)/spaces.o $(OBJ_DIR)/model.o $(OBJ_DIR)/model_LiebLiniger.o $(OBJ_DIR)/matrix.o $(OBJ_DIR)/state.o $(OBJ_DIR)/tba_LiebLiniger.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/matrixelements_LiebLiniger.o: $(SRC_DIR)/LiebLiniger/repulsive/matrixelements_LiebLiniger.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/matrix.o $(OBJ_DIR)/state_LiebLiniger.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/calculus.o: $(SRC_DIR)/math/calculus.cc $(OBJ_DIR)/conveniences.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/tba_LiebLiniger.o: $(SRC_DIR)/LiebLiniger/repulsive/tba_LiebLiniger.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/calculus.o + $(COMPILE) -c $< -o $@ + +$(OBJ_DIR)/abacus.o: $(SRC_DIR)/abacus.cc $(OBJ_DIR)/conveniences.o $(OBJ_DIR)/timer.o $(OBJ_DIR)/labels.o $(OBJ_DIR)/spaces.o $(OBJ_DIR)/model.o $(OBJ_DIR)/model_LiebLiniger.o $(OBJ_DIR)/quantumnumbers.o $(OBJ_DIR)/plex.o $(OBJ_DIR)/matrix.o $(OBJ_DIR)/state.o $(OBJ_DIR)/state_LiebLiniger.o $(OBJ_DIR)/matrixelements_LiebLiniger.o $(OBJ_DIR)/calculus.o $(OBJ_DIR)/tba_LiebLiniger.o + $(COMPILE) -c $< -o $@ + + +########################################### +# g++ precompiled standard library elements + +.PHONY: gcm +gcm: + rm -rf gcm.cache + rm -f obj/std.o + $(CXX) -std=c++23 -fmodules -fsearch-include-path $(OPTIMIZATIONS) -c bits/std.cc + mv std.o obj/ + + +########################################### +# clang++ std precompilation +# NOTE: +# In Linux, this requires installation of llvm-18, to ensure std.cppm is available. +# Installing on other platforms might require specifying the location of std.cppm differently + +.PHONY: clang-std +clang-std: + clang++-21 -std=c++23 -stdlib=libc++ --precompile -o obj/std.pcm /usr/lib/llvm-21/share/libc++/v1/std.cppm + + +################# +# General cleanup + +.PHONY: clean +clean: + rm -f abacus.makeinfo + rm -rf gcm.cache + rm -rf $(DEP_DIR) + rm -f $(OBJ_DIR)/* + rm -f $(LIB_DIR)/lib$(VERSION).a + rm -rf $(BIN_DIR) diff --git a/README.md b/README.md new file mode 100644 index 0000000..feac3d8 --- /dev/null +++ b/README.md @@ -0,0 +1,33 @@ +# Abacus (Lieb-Liniger minimal) + +Computational tools for Bethe Ansatz-solvable models. + +This is a stripped-down version of the full Abacus stack (in development), containing only basic features for (repulsive) Lieb-Liniger. + +Licensed under the [AGPLv3](https://www.gnu.org/licenses/agpl-3.0.en.html) license. Copyright © [Jean-Sébastien Caux](https://jscaux.org). + + +## Installation + +You will need a `c++23`-enabled compiler. Compiling works for GNU `g++` version `15` and above. + +To precompile the `std` standard library module. From the base directory, simply run + +``` shell +$ make gcm +``` + +`Abacus` itself is then built by invoking + +``` shell +$ make +``` + +This will produce all executables (located in the `bin` folder), together with a library `abacus-[vn].a` (located in the `lib` folder), where `[vn]` follows semantic versioning conventions. + +##### Executables +All executables are in the `bin/` folder. Invoking them with no arguments will print out usage instructions. + + + +An outdated description of the original (and deprecated) version of ABACUS can be found in J.-S. Caux, J. Math. Phys. 50, 095214 (2009), [doi:10.1063/1.3216474](https://doi.org/10.1063/1.3216474). diff --git a/src/LiebLiniger/repulsive/matrixelements_LiebLiniger.cc b/src/LiebLiniger/repulsive/matrixelements_LiebLiniger.cc new file mode 100644 index 0000000..eab6293 --- /dev/null +++ b/src/LiebLiniger/repulsive/matrixelements_LiebLiniger.cc @@ -0,0 +1,172 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module matrixelements_LiebLiniger; + + +import std; + + +import conveniences; +import matrix; +import state_LiebLiniger; + + +/////////////////////// +// ↓ Matrix elements // +/////////////////////// + + +// Density operator ρ (x=0) + +std::complex V_ρ +(IndexU α, const LiebLinigerBetheState& bra, const LiebLinigerBetheState& ket) +{ + std::complex result { Real(1) }; + for (IndexU β {0}; β < ket.ƛ_.size(); ++β) + result *= (bra.ƛ_[β] - ket.ƛ_[α] + 1_ir)/(ket.ƛ_[β] - ket.ƛ_[α] + 1_ir); + return(result); +} + +export std::complex matrix_element_ρ +(const LiebLinigerBetheState& bra, const LiebLinigerBetheState& ket) +{ + if (bra == ket) return bra.ρ(); + else if (bra.N() != ket.N() || bra.iK() == ket.iK()) return std::complex(0); + + Matrix one_plus_U (bra.N()); + + std::vector> Vplus (bra.ƛ_.size()); + + std::vector Fn_Prod (bra.ƛ_.size()); + + std::vector rKern (bra.ƛ_.size()); + + // "Phantom" rapidity in ME expression. Choice doesn't matter, + // see 1990_Slavnov_TMP_82 after (3.8). Choose rapidity around the middle. + IndexU p = ket.N()/2-1; + + IndexU a { 0 }; + + for (a = 0; a < bra.ƛ_.size(); ++a) { + Vplus[a] = V_ρ (a, bra, ket); + Fn_Prod[a] = Real(1); + for (IndexU m {0}; m < bra.ƛ_.size(); ++m) + if (m != a) Fn_Prod[a] *= (bra.ƛ_[m] - ket.ƛ_[a])/(ket.ƛ_[m] - ket.ƛ_[a]); + rKern[a] = -ket.model_.get().dφdƛ_(ket.ƛ_[a] - ket.ƛ_[p]); + } + + for (a = 0; a < bra.ƛ_.size(); ++a) + for (IndexU b = 0; b < bra.ƛ_.size(); ++b) + one_plus_U(a,b) = (a == b ? Real(1) : Real(0)) + + Real(0.5) * ((bra.ƛ_[a] - ket.ƛ_[a])/imag(Vplus[a])) + * Fn_Prod[a] * (-ket.model_.get().dφdƛ_(ket.ƛ_[a] - ket.ƛ_[b]) - rKern[b]); + + std::complex ddalpha_sigma { std::exp(one_plus_U.lndet_LU_destroy()) }; + + std::complex ln_prod_V { Real(0) }; + for (IndexU a {0}; a < Vplus.size(); ++a) ln_prod_V += log(2_ir * imag(Vplus[a])); + + std::complex ln_prod_2 { Real(0) }; + for (IndexU a {0}; a < bra.ƛ_.size(); ++a) + for (IndexU b {0}; b < bra.ƛ_.size(); ++b) + ln_prod_2 += log((ket.ƛ_[a] - ket.ƛ_[b] + 1_ir)/(bra.ƛ_[a] - ket.ƛ_[b])); + + return (ket.K() - bra.K()) * ddalpha_sigma * + std::exp(ln_prod_V + ln_prod_2 - Real(0.5)*(bra.lnnorm_ + ket.lnnorm_)) + /(2 * imag(Vplus[p])); +} + + +// Field annihilation operator ψ (x=0) + +std::complex V_ψ +(IndexU α, const LiebLinigerBetheState& bra, const LiebLinigerBetheState& ket) +{ + std::complex result { Real(1) }; + for (IndexU β {0}; β < bra.ƛ_.size(); ++β) + result *= (bra.ƛ_[β] - ket.ƛ_[α] + 1_ir)/(ket.ƛ_[β] - ket.ƛ_[α] + 1_ir); + result /= ket.ƛ_.back() - ket.ƛ_[α] + 1_ir; + return(result); +} + +export std::complex matrix_element_ψ +(const LiebLinigerBetheState& bra, const LiebLinigerBetheState& ket) +{ + if (bra.N() + 1 != ket.N()) return std::complex(0); + + Matrix U (bra.N()); + + std::vector> Vplus (bra.N()); + + std::vector Fn_Prod (bra.N()); + + std::vector rKern (bra.N()); + + // "Phantom" rapidity in ME expression. Choice doesn't matter, + // see 1990_Slavnov_TMP_82 after (3.8). + IndexU p = ket.N()-1; + + for (IndexU a {0}; a < bra.ƛ_.size(); ++a) + { + Vplus[a] = V_ψ (a, bra, ket); + Fn_Prod[a] = (bra.ƛ_[a] - ket.ƛ_[a])/(ket.ƛ_.back() - ket.ƛ_[a]); + + for (IndexU m {0}; m < bra.ƛ_.size(); ++m) + if (m != a) Fn_Prod[a] *= (bra.ƛ_[m] - ket.ƛ_[a])/(ket.ƛ_[m] - ket.ƛ_[a]); + + rKern[a] = -ket.model_.get().dφdƛ_(ket.ƛ_[a] - ket.ƛ_[p]); + } + + for (IndexU a {0}; a < bra.ƛ_.size(); ++a) + { + for (IndexU b = 0; b < bra.ƛ_.size(); ++b) + { + U(a,b) = (a == b ? Real(2)*imag(Vplus[a]) : Real(0)) + + Fn_Prod[a] * (-ket.model_.get().dφdƛ_(ket.ƛ_[a] - ket.ƛ_[b]) - rKern[b]); + } + } + + // std::complex det_U { U.determinant() }; + std::complex det_U { std::exp(U.lndet_LU_destroy()) }; + + Real ln_prod_ƛsq_plus_1 { Real(0) }; + + for (IndexU a {0}; a+1 < ket.ƛ_.size(); ++a) + { + for (IndexU b {a+1}; b < ket.ƛ_.size(); ++b) + ln_prod_ƛsq_plus_1 += std::logl(std::pow(ket.ƛ_[a] - ket.ƛ_[b], 2) + Real(1)); + } + + std::complex ln_prod_ƛa_min_μb { Real(0) }; + + for (IndexU a {0}; a < ket.ƛ_.size(); ++a) + { + for (IndexU b {0}; b < bra.ƛ_.size(); ++b) + ln_prod_ƛa_min_μb += log(std::complex(ket.ƛ_[a] - bra.ƛ_[b])); + } + + return (det_U * std::sqrt(bra.model_.get().c_) * + std::exp(ln_prod_ƛsq_plus_1 - ln_prod_ƛa_min_μb + - Real(0.5)*(bra.lnnorm_ + ket.lnnorm_))); +} + + +// Field creation operator ψdag (x=0) + +export std::complex matrix_element_ψdag +(const LiebLinigerBetheState& bra, const LiebLinigerBetheState& ket) +{ + return conj(matrix_element_ψ(ket, bra)); +} + + +/////////////////////// +// ↑ Matrix elements // +/////////////////////// diff --git a/src/LiebLiniger/repulsive/model_LiebLiniger.cc b/src/LiebLiniger/repulsive/model_LiebLiniger.cc new file mode 100644 index 0000000..678ada2 --- /dev/null +++ b/src/LiebLiniger/repulsive/model_LiebLiniger.cc @@ -0,0 +1,127 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module model_LiebLiniger; + + +import std; + + +import conveniences; +import spaces; +import model; + +////////////////////////////// +// ↓ Class LiebLinigerModel // +////////////////////////////// + +// forward declaration for friendship +// class LiebLinigerBetheState; + + +export class LiebLinigerModel : public Model +{ + +public: // public interface + + // constructors + LiebLinigerModel (BosonicContinuum& space, Real c) + : Model(space, Model::Type::LiebLiniger, c*space.L()) + , c_ { verify_c(c) } + , cxL_ { c*space.L() } + {} + + // utilities + // LiebLinigerModel& operator= (const LiebLinigerModel& m); + bool operator== (const LiebLinigerModel& rhs) const; + + std::string get_filename_prefix () const override; + + // friendship + // template TModel> + // friend class BetheState; + // friend class LiebLinigerBetheState; + + // friend std::complex matrix_element_ρ + // (const LiebLinigerBetheState& bra, const LiebLinigerBetheState& ket); + + // friend std::complex matrix_element_ψ + // (const LiebLinigerBetheState& bra, const LiebLinigerBetheState& ket); + + // friend std::complex matrix_element_ψdag + // (const LiebLinigerBetheState& bra, const LiebLinigerBetheState& ket); + + +public: // protected: + Real c_; + Real cxL_; + +private: + static Real verify_c (Real c) + { + if (c < 0.0L) + { + throw std::invalid_argument("c must be greater than 0 in LiebLinigerModel"); + } + return c; + } + +public: // protected: + Real θ_ (Real ƛ) const override { return ƛ; } + Real θinv_ (Real ƛ) const override { return ƛ; } + Real dθdƛ_ ([[maybe_unused]] Real ƛ) const override { return Real(1); } + // Real φ_ (Real ƛ) const override { return -2*std::atan(ƛ); } + long double φ_ (long double ƛ) const override { return -2*std::atan(ƛ); } + double φ_ (double ƛ) const override { return -2*std::atan(ƛ); } + Real dφdƛ_ (Real ƛ) const override { return -2/(ƛ*ƛ + 1); } + Real θ_ ([[maybe_unused]] int nj, + [[maybe_unused]] int pj, + [[maybe_unused]] Real ƛ) const override { return Real(0); }; + Real θinv_ ([[maybe_unused]] int nj, + [[maybe_unused]] int pj, + [[maybe_unused]] Real ƛ) const override { return Real(0); }; + Real dθdƛ_ ([[maybe_unused]] int nj, + [[maybe_unused]] int pj, + [[maybe_unused]] Real ƛ) const override { return Real(0); }; + Real φ_ ([[maybe_unused]] IndexU k, + [[maybe_unused]] Real ƛ) const override { return Real(0); }; + Real dφdƛ_ ([[maybe_unused]] IndexU k, + [[maybe_unused]] Real ƛ) const override { return Real(0); }; + Real φ_ ([[maybe_unused]] IndexU j, + [[maybe_unused]] IndexU k, + [[maybe_unused]] Real ƛ) const override { return Real(0); }; + Real dφdƛ_ ([[maybe_unused]] IndexU j, + [[maybe_unused]] IndexU k, + [[maybe_unused]] Real ƛ) const override { return Real(0); }; +}; + +// LiebLinigerModel& LiebLinigerModel::operator= (const LiebLinigerModel& m) +// { +// if (space_ != m.space_ || type_ != m.type_ || Ł_ != m.Ł_ || c_ != m.c_) +// throw "Cannot change LiebLinigerModel by assignment"; + +// return *this; +// } + +bool LiebLinigerModel::operator== (const LiebLinigerModel& rhs) const +{ + return (space_ == rhs.space_ && c_ == rhs.c_); +} + +std::string LiebLinigerModel::get_filename_prefix () const +{ + std::stringstream output; + output << get_modelname_prefix() << "_c_" << c_ << "_L_" << space_.L() ; + return output.str(); +} + + +////////////////////////////// +// ↑ Class LiebLinigerModel // +////////////////////////////// diff --git a/src/LiebLiniger/repulsive/state_LiebLiniger.cc b/src/LiebLiniger/repulsive/state_LiebLiniger.cc new file mode 100644 index 0000000..aab615e --- /dev/null +++ b/src/LiebLiniger/repulsive/state_LiebLiniger.cc @@ -0,0 +1,343 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module state_LiebLiniger; + + +import std; + + +import conveniences; +import labels; +import spaces; +import model; +import model_LiebLiniger; +import state; +import tba_LiebLiniger; + + +/////////////////////////////////// +// ↓ Class LiebLinigerBetheState // +/////////////////////////////////// + +export class LiebLinigerBetheState : public BetheState { + +public: // public interface + + // constructors + LiebLinigerBetheState (LiebLinigerModel& model, int N); + + LiebLinigerBetheState + (const LiebLinigerBetheState& refstate, int Δf, bool relative_label); + + LiebLinigerBetheState + (const LiebLinigerBetheState& refstate, + std::string baselabel, bool relative_label); + + LiebLinigerBetheState + (LiebLinigerModel& model, TBASolutionLiebLiniger& tbasol); + + // LiebLinigerBetheState& operator= (const LiebLinigerBetheState& rhs); + + // physical properties + Real c () const { return model_.get().c_; } + Real L() const { return model_.get().space_.L(); }; + int N () const { return N_; }; + Real ρ () const { return N_/model_.get().space_.L(); }; + Real λF () const; + Real ln_subspace_dim_this_filling () const; + std::vector spectrum_1ph (int iKexc) const; // returns ordered set of energies of 1ph excitations + + // utilities + std::string get_filename_prefix () const override; + + // friendship + friend std::ostream& operator<< (std::ostream& s, const LiebLinigerBetheState& state); + +protected: + //const int N_; + int N_; + +private: + static int verify_N (int N) { + if (N <= 0) { + throw std::invalid_argument("N must be greater than 0 in LiebLinigerBetheState"); + } + return N; + } + static int get_N (Real L, TBASolutionLiebLiniger& tbasol) { + return int(L * tbasol.ρ() + 0.5); + } + +public: + void set_ground_state_Ix2 () override; + void set_Ix2_from_tba (TBASolutionLiebLiniger& tbasol); + void initialize () override; + bool Gaudin_g_left_edge_is_monotonic () const override { return true; }; + bool Gaudin_g_right_edge_is_monotonic () const override { return true; }; + void control_δƛ () override {}; // nothing to control here + void compute_lnnorm () override; + void compute_Momentum () override; + void compute_Energy () override; + void populate_λ () override; +}; + +LiebLinigerBetheState::LiebLinigerBetheState +(LiebLinigerModel& model, int N) + : BetheState(model, N) + , N_ { verify_N(N) } +{ + set_ground_state_Ix2(); + initialize(); + compute_Momentum(); +} + +LiebLinigerBetheState::LiebLinigerBetheState +(const LiebLinigerBetheState& refstate, int Δf, bool relative_label) + : BetheState(refstate, Δf, relative_label) + , N_ { refstate.N_ + Δf } +{ + initialize(); + compute_Momentum(); +} + +LiebLinigerBetheState::LiebLinigerBetheState +(const LiebLinigerBetheState& refstate, std::string baselabel, bool relative_label) + : LiebLinigerBetheState + (refstate, refstate.model_.get().g_f_from_label(baselabel) - refstate.f_, relative_label) +{} + +LiebLinigerBetheState::LiebLinigerBetheState +(LiebLinigerModel& model, TBASolutionLiebLiniger& tbasol) + : BetheState(model, get_N(model.space_.L(), tbasol)) + , N_ { verify_N(get_N(model.space_.L(), tbasol)) } +{ + set_Ix2_from_tba (tbasol); + tags_.push_back(std::make_pair("T", tbasol.T_)); + initialize(); + compute_Momentum(); +} + +// LiebLinigerBetheState& LiebLinigerBetheState::operator= +// (const LiebLinigerBetheState& rhs) +// { +// if (model_ != rhs.model_ || N_ != rhs.N_) { +// throw AbacusException("Cannot assign to LiebLinigerBetheState from state with different model/filling"); +// } + +// Ix2_ = rhs.Ix2_; +// tags_ = rhs.tags_; +// λ = rhs.λ; +// label_ = rhs.label_; +// baselabel_ = rhs.baselabel_; +// patternlabel_ = rhs.patternlabel_; +// δB_ = Real(1); +// converged_ = false; + +// return *this; +// } + +std::string LiebLinigerBetheState::get_filename_prefix () const { + std::stringstream output; + output << model_.get().get_filename_prefix() << "_N_" << N_; + + if (!tags_.empty()) { + for (auto tag : tags_) output << "_" << tag.first << "_" << tag.second; + } + else output << "_" << label_; + + return output.str(); +} + +void LiebLinigerBetheState::set_ground_state_Ix2 () { + for (IndexU α { 0 }; α < ƛ_.size(); ++α) Ix2_[α] = -N_ + 1 + 2*ι(α); +} + +void LiebLinigerBetheState::set_Ix2_from_tba (TBASolutionLiebLiniger& tbasol) +{ + // The counting function is defined as c(λ) = L ∫_-∞^λ dλ' ρ(λ') + // Logic: we pop an occupation each time c(\lambda) crosses value (integer - 1/2) + IndexU nfound { 0 }; + + std::vector x_found; + Real count { 0 }; + Real count_prev; + + for (IndexU i { 0 }; i < tbasol.ρ_.interval_.λ.size(); ++i) { + count_prev = count; + count += L() * tbasol.ρ_.interval_.dλ[i] * tbasol.ρ(i); + // if (count > nfound + Real(1.5)) { // more than 1 rapidity in this dλ element + // std::cerr << "In LiebLinigerBetheState::set_Ix2_from_tba with " + // << " L = " << L() + // << ", count - nfound = " << count-nfound + // << " so there is more than one rapidity in an integration element L ρ dλ " + // << "(" << std::floor(count + Real(0.5) - nfound) << " were found).\n" + // << "-> to ensure a reliable discretization, " + // << "it is advisable to improve the accuracy of the " + // << "TBASolutionLiebLiniger& tbasol argument." << std::endl; + // //throw; + // } + //if (count > nfound + Real(0.5)) { + for (int n { 0 }; n < std::floor(count + Real(0.5) - nfound); ++n) { + // N.B.: we consider that the integral carried by count takes its value at λ[i]+0.5dλ[i]. + // The found rapidity is determined by linear interpolation, and thus solves + // nfound + 0.5 = count_prev + (λ - (λ[i]-0.5dλ[i]))*(count - count_prev)/dλ[i] + // so λ = λ[i] - 0.5dλ[i] + (nfound + 0.5 - count_prev)*dλ[i]/(count - count_prev). + // We then immediately compute the counting function at this found rapidity. + x_found.push_back(tbasol.x(tbasol.ρ_.interval_.λ[i] - Real(0.5)*tbasol.ρ_.interval_.dλ[i] + + (nfound + Real(0.5) - count_prev)*tbasol.ρ_.interval_.dλ[i]/(count - count_prev))); + nfound++; + } + } + + // The quantum numbers are then the Ix2 (with correct parity) closest to the found L*x + std::vector Ix2_found; + for (Real x : x_found) { + // Logic: 2Lx between Ix2-1 and Ix2+1 gets associated to Ix2 + // Let Ix2 = 2n + 1-N%2, so 2Lx between 2n-N%2 and 2n-N%2+2 gets associated to n, + // or Lx between n-0.5*(N%2) and n + -0.5(N%2) + 1 or Lx + 0.5*(N%2) between n and n+1. + // Thus, n = floor(Lx+ 0.5*(N%2)) and Ix2 = 1-N%2 + 2*floor(Lx+0.5*(N%2)) + Ix2_found.push_back(1 - N_%2 + 2*std::floor(L() * x + Real(0.5)*(N_%2))); + } + // Check that the Ix2_found is symmetric: + for (IndexU i { 0 }; i < Ix2_found.size()/2; ++i) + if (Ix2_found[i] != -Ix2_found[Ix2_found.size() - 1 - i]) { + std::cout << "LiebLinigerBetheState::set_Ix2_from_tba yielded " + << "an asymmetric state at L = " << L() << "\n" << Ix2_found << "\n" + << "i = " << i << "\tIx2[i] = " << Ix2_found[i] + << ", Ix2[size-1 - i] = " << Ix2_found[Ix2_found.size()-1 - i] + << "\n-> to ensure a reliable discretization, " + << "it is advisable to improve the accuracy of the " + << "TBASolutionLiebLiniger& tbasol argument." << std::endl; + throw AbacusException(""); + } + Ix2_ = Ix2_found; +} + +void LiebLinigerBetheState::initialize () { + if (model_.get().c_ > 1.0L) { + for (IndexU α { 0 }; α < ƛ_.size(); ++α) ƛ_[α] = pi_r * Ix2_[α]/model_.get().cxL_; + } + else { + // For small values of c, use better approximation using approximate + // zeroes of Hermite polynomials: see Gaudin eqn 4.71. + Real f = 1.0L/std::pow(model_.get().cxL_ * N_, 0.5); + for (IndexU α { 0 }; α < ƛ_.size(); ++α) ƛ_[α] = pi_r * Ix2_[α] * f; + } + compute_B_(); + std::fill(iter_count_.begin(), iter_count_.end(), 0); + std::fill(iter_time_.begin(), iter_time_.end(), 0.0); +} + +void LiebLinigerBetheState::populate_λ() { + for (IndexU α { 0 }; α < ƛ_.size(); ++α) λ[α] = model_.get().c_ * ƛ_[α]; +} + +Real LiebLinigerBetheState::ln_subspace_dim_this_filling () const { + return ln_dim_; +} + +std::vector LiebLinigerBetheState::spectrum_1ph (int iKexc) const { + // returns ordered set of energies of 1ph excitations + LiebLinigerBetheState estate(*this, 0, false); + estate.set_g_Ix2(this->Ix2_); + std::vector spectrum; + for (IndexU α {0}; α < Ix2_.size(); ++α) if (estate.excite_g (α, iKexc)) { + estate.polish(); + spectrum.push_back(estate.E() - E()); + estate.set_g_Ix2(this->Ix2_); + } + return spectrum; +} + +void LiebLinigerBetheState::compute_lnnorm() { + // lnnorm_ = std::logl(Gaudin_det_) + charge_ * std::logl(model_.get().Ł_); + lnnorm_ = ln_Gaudin_det_ + charge_ * std::logl(model_.get().Ł_); + for (IndexU α {0}; α+1 < ƛ_.size(); ++α) + for (IndexU β {α+1}; β < ƛ_.size(); ++β) + lnnorm_ += std::logl(Real(1) + Real(1)/std::pow(ƛ_[α] - ƛ_[β], Real(2))); +} + +void LiebLinigerBetheState::compute_Momentum() { + iK_ = 0; + for (int i : Ix2_) iK_ += i; + iK_ /= 2; // because we summed Ix2, not I + K_ = twopi_r * iK_/L(); +} + +void LiebLinigerBetheState::compute_Energy() { + E_ = Real(0); + for (Real ƛ : ƛ_) E_ += ƛ*ƛ; + E_ *= model_.get().c_ * model_.get().c_; +} + +Real LiebLinigerBetheState::λF() const { // valid only if this is the ground state + return -0.5L*(3*λ[0] - λ[1]); +} + +export std::ostream& operator<< (std::ostream& s, const LiebLinigerBetheState& state) { + s << "label: " << state.label() << "\tconverged: " << state.converged() << "\tδB: " << state.δB() << "\t" << std::numeric_limits::epsilon()*1000*state.charge_ << "\n"; + for (int n : state.Ix2_) s << n << "\t"; + s << std::endl; + for (Real l : state.ƛ_) s << l << "\t"; + s << std::endl; + s << std::accumulate(state.iter_count_.begin(), state.iter_count_.end(), 0) << "\t" << state.δB_ << std::endl; + return s; +} + +/////////////////////////////////// +// ↑ Class LiebLinigerBetheState // +/////////////////////////////////// + + +export LiebLinigerBetheState thermal_state (LiebLinigerModel& model, int N, Real T) +{ + if (T > 0) { + // try to fetch the Ix2 from file + std::stringstream filename; + filename << model.get_filename_prefix() << "_N_" << N << "_T_" << T << ".Ix2"; + + std::ifstream infile; + infile.open(filename.str()); + + if (!infile.fail()) { + + std::vector Ix2; + int Ix2_read; + + for (int count { 0 }; count < N; ++count) { + infile >> Ix2_read; + Ix2.push_back(Ix2_read); + } + infile.close(); + + LiebLinigerBetheState llbs(model, N); + llbs.Ix2_ = Ix2; + llbs.tags_.push_back(std::make_pair("T", T)); + llbs.initialize(); + llbs.compute_Momentum(); + + return llbs; + } + else { // no infile, so we build the thermal state from TBA + TBASolutionLiebLiniger tbasol { + //tba_solution_LiebLiniger_at_filling(model.c_, T, N/model.space_.L(), 1.0e-6) + tba_solution_LiebLiniger_at_filling(model.c_, T, N/model.space_.L()) + }; + LiebLinigerBetheState llbs(model, tbasol); + std::ofstream outfile; + outfile.open(filename.str()); + for (int Ix2 : llbs.Ix2_) outfile << "\t" << Ix2; + outfile.close(); + return llbs; + } + } + // else T == 0 so return the default ground state + return LiebLinigerBetheState(model, N); +} diff --git a/src/LiebLiniger/repulsive/tba_LiebLiniger.cc b/src/LiebLiniger/repulsive/tba_LiebLiniger.cc new file mode 100644 index 0000000..3c0eae9 --- /dev/null +++ b/src/LiebLiniger/repulsive/tba_LiebLiniger.cc @@ -0,0 +1,416 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module tba_LiebLiniger; + + +import std; + + +import conveniences; +import calculus; + + +export class TBASolutionLiebLiniger { + +public: + Real c_; + Real T_; + Real μ_; + Function ε_; + Function dεdμ_; + Function ρ_; + Function ρh_; + + // Computational conveniences + Function predecessor_ε_; + Function precomputed_; + Function precomputed_ddλ_; + std::vector Cauchy_Exact_0; + std::vector Cauchy_Exact_1; + +public: + TBASolutionLiebLiniger(Real c, Real T, Real μ); + + Real ϕ (Real λ) { return 2*std::atan(λ/c_); } + Real Cauchy (Real λ) { return c_/(pi_r * (λ*λ + c_*c_)); } + void refine_interval (); + void reset_μ (Real new_μ); + Real Tln1pem (Real& ε); + void iterate_ε_old (); + void iterate_ε (); + void solve_for_ε_diagonal (IndexU& i, Real& Ri); + bool solve_for_ε (Real eps = Real(0)); + void converge_ε (Real eps = Real(0)); + void iterate_dεdμ_old (); + void iterate_dεdμ (); + bool solve_for_dεdμ (Real eps = Real(0)); + void populate_ρ_and_ρh (); + void solve (Real eps = Real(0)); + Real evaluate_convergence (); + Real ρ (IndexU i) { return ρ_.val_[i]; } + Real μ () { return μ_; } + + // physical properties + Real g (); /// free energy density + Real ρ (); /// total density + Real x (Real λ); /// counting function +}; + +TBASolutionLiebLiniger::TBASolutionLiebLiniger (Real c, Real T, Real μ) + : c_(c) + , T_(T) + , μ_(μ) +{ + Interval interval_( + // Choose initial λmax such that e^{-λmax^2/T} << accuracy + pi_r + std::sqrt(-T_*std::log(std::numeric_limits::epsilon())), + // and nr of points very small to begin with + 1+2*int(pi_r+std::sqrt(-T_*std::log(std::numeric_limits::epsilon())))); + + ε_ = Function(interval_); + for (IndexU i { 0 }; i < ε_.interval_.λ.size(); ++i) { + ε_.val_[i] = ε_.interval_.λ[i] * ε_.interval_.λ[i] - μ_; + } + +} + +void TBASolutionLiebLiniger::refine_interval () +{ + Real λmax { ε_.interval_.λ.back() + 0.5* ε_.interval_.dλ.back() }; + // If the endpoint values of ε are not large enough, increase λmax further + if (ε_.val_.front() < -T_*std::log(std::numeric_limits::epsilon())) + λmax *= Real(1.1); + + IndexU npts { ε_.interval_.λ.size() }; + + predecessor_ε_ = std::move(ε_); + + // Refine grid for ε, and initiate values based on predecessor + ε_ = Function(Interval(λmax, 1 + 2*int(0.7*npts))); + + for (IndexU i { 0 }; i < ε_.interval_.λ.size(); ++i) { + ε_.val_[i] = predecessor_ε_.evaluate_at(ε_.interval_.λ[i]); + } +} + +void TBASolutionLiebLiniger::reset_μ (Real new_μ) +{ + // Reset the chemical potential, + // guess new value of ε based on linear fit, + // and mark all functions as not converged + Real Δμ { new_μ - μ_ }; + μ_ = new_μ; + + Function estimated_ε (ε_.interval_); + for (IndexU i { 0 }; i < ε_.interval_.λ.size(); ++i) { + estimated_ε.val_[i] = ε_.val_[i] + dεdμ_.val_[i] * Δμ; + } + + // Downgrade ε to base λmax and number of points as in constructor + ε_ = Function + (Interval(pi_r + std::sqrt(-T_*std::log(std::numeric_limits::epsilon())), + 1+2*int(pi_r+std::sqrt(-T_*std::log(std::numeric_limits::epsilon()))))); + + // Populate based on estimated values: + for (IndexU i { 0 }; i < ε_.interval_.λ.size(); ++i) { + ε_.val_[i] = estimated_ε.evaluate_at(ε_.interval_.λ[i]); + } +} + +/// Returns the integral of absolute value difference of +/// free energy integrand between most recent ε and predecessor ε +/// evaluated over the predecessor ε's set of points +Real TBASolutionLiebLiniger::evaluate_convergence () +{ + Real diff { 0 }; + Real latest_ε; + for (IndexU i { 0 }; i < predecessor_ε_.interval_.λ.size(); ++i) { + latest_ε = ε_.evaluate_at(predecessor_ε_.interval_.λ[i]); + diff += predecessor_ε_.interval_.dλ[i] * + std::abs(Tln1pem(predecessor_ε_.val_[i]) - Tln1pem(latest_ε)); + } + return T_ * diff/twopi_r; +} + +/// Iterates ε using a diagonal Newton-style step +void TBASolutionLiebLiniger::iterate_ε () +{ + ε_.val_prev_.swap(ε_.val_); + ε_.δval_ = Real(0); + + // Set the values of Tln(1+e^{-ε/T}) from previous iteration + for (IndexU i { 0 }; i < ε_.interval_.λ.size(); ++i) { + precomputed_.val_[i] = Tln1pem(ε_.val_prev_[i]); + } + + Real Ri; // remainder, which we aim to bring down to zero + + for (IndexU i { 0 }; i < ε_.interval_.λ.size(); ++i) { + Ri = ε_.interval_.λ[i] * ε_.interval_.λ[i] - μ_; + + for (IndexU j { 0 }; j < i; ++j) { + Ri -= Cauchy_Exact_0[i-j] * precomputed_.val_[j]; + } + for (IndexU j { i+1 }; j < ε_.interval_.λ.size(); ++j) { + Ri -= Cauchy_Exact_0[j-i] * precomputed_.val_[j]; + } + + solve_for_ε_diagonal(i, Ri); + + // Simple measure of convergence: ε not ideal, use free energy weight + ε_.δval_ += std::abs(Tln1pem(ε_.val_[i]) - precomputed_.val_[i]); + } // for i +} + +Real TBASolutionLiebLiniger::Tln1pem (Real& ε) +{ + return ε > Real(0)? + (ε > -T_* ln_sqrt_real_eps ? + T_ * std::exp(-ε/T_) : T_ * std::log(1 + std::exp(-ε/T_))) + : + -ε + (ε < T_* ln_sqrt_real_eps ? + T_ * std::exp(ε/T_) : T_ * std::log(1 + std::exp(ε/T_))); +} + +void TBASolutionLiebLiniger::solve_for_ε_diagonal (IndexU& i, Real& Ri) +{ + // Solve f(ε) = Ri for ε, where f(ε) = ε + Δ Tln(1+e^{-ε/T}) + + Real Δ { 2*std::atan(ε_.interval_.dλ[i]/(2*c_)) * oneoverpi_r }; + Real δε, dfdε; + + Real prev_diff { std::numeric_limits::max() }; + Real diff { prev_diff/2 }; + + int niter { 0 }; + + while (diff < prev_diff) { + prev_diff = diff; + dfdε = 1 - Δ/(1 + std::exp(ε_.val_[i]/T_)); + δε = (Ri - (ε_.val_[i] + Δ*Tln1pem(ε_.val_[i])))/dfdε; + diff = std::abs(δε); + ε_.val_[i] += δε; + } +} + +bool TBASolutionLiebLiniger::solve_for_ε (Real eps) +{ + precomputed_ = Function(ε_.interval_); + precomputed_ddλ_ = Function(ε_.interval_); + + Cauchy_Exact_0 = std::vector(ε_.interval_.λ.size()); + + Real dλ { ε_.interval_.dλ[0] }; + Cauchy_Exact_0[0] = oneoverpi_r * Real(2) * std::atan(dλ/(2*c_)); + for (int i { 1 }; i < ε_.interval_.λ.size(); ++i) { + Cauchy_Exact_0[i] = oneoverpi_r * std::atan(c_ * dλ/(c_*c_ + (4*i*i - 1)*dλ*dλ/4)); + } + + IndexU niter { 0 }; + Real prev_δval_ { }; + do { + prev_δval_ = ε_.δval_; + iterate_ε(); + niter++; + } while (ε_.δval_ > eps && ε_.δval_ < prev_δval_); // if improving, keep going + + return ε_.δval_ < std::sqrt(std::numeric_limits::epsilon()); +} + +void TBASolutionLiebLiniger::converge_ε (Real eps) +{ + solve_for_ε(eps); + + Real previous_Gibbs { std::numeric_limits::max() }; + Real Gibbs { 0 }; + Real previous_Δ { std::numeric_limits::max() }; + Real Δ { previous_Δ/2 }; + + while (Δ > eps && + (Δ < previous_Δ || ε_.val_.size() < 1000) + && ε_.val_.size() < 8000 + ) { + refine_interval(); + solve_for_ε(eps); + previous_Gibbs = Gibbs; + Gibbs = g(); + previous_Δ = Δ; + Δ = evaluate_convergence(); + // std::cout << "converge_ε: nr pts " << ε_.interval_.dλ.size() + // << "\t" << ε_.interval_.dλ[(ε_.val_.size() - 1)/2] + // << "\t" << ε_.interval_.λ[(ε_.val_.size() - 1)/2] + // << "\t" << ε_.val_[(ε_.val_.size() - 1)/2] << "\n" + // << "\tg " << std::setprecision(20) << Gibbs + // << "\tG-pG " << std::setprecision(8) << Gibbs-previous_Gibbs + // << "\tΔ " << Δ << "\tprevious_Δ " << previous_Δ << std::endl; + }; + + // If the last iteration made things worse, pull back + if (Δ < previous_Δ) { + ε_.val_.swap(ε_.val_prev_); + } +} + +void TBASolutionLiebLiniger::iterate_dεdμ () +{ + dεdμ_.val_prev_.swap(dεdμ_.val_); + dεdμ_.δval_ = Real(0); + + for (IndexU i { 0 }; i < dεdμ_.interval_.λ.size(); ++i) { + precomputed_.val_[i] = dεdμ_.val_prev_[i] * + (ε_.val_[i] > Real(0) ? std::exp(-ε_.val_[i]/T_)/(1 + std::exp(-ε_.val_[i]/T_)) + : Real(1)/(1 + std::exp(ε_.val_[i]/T_))); + } + + Real Ri; + + for (IndexU i { 0 }; i < dεdμ_.interval_.λ.size(); ++i) { + Ri = -Real(1); + + for (IndexU j { 0 }; j < i; ++j) { + Ri += Cauchy_Exact_0[i-j] * precomputed_.val_[j]; + } + for (IndexU j { i+1 }; j < ε_.interval_.λ.size(); ++j) { + Ri += Cauchy_Exact_0[j-i] * precomputed_.val_[j]; + } + + // directly solve since this is a linear system (for fixed ε): + dεdμ_.val_[i] = Ri/(1 - + 2*std::atan(ε_.interval_.dλ[i]/(2*c_)) * oneoverpi_r * + (ε_.val_[i] > Real(0) ? + std::exp(-ε_.val_[i]/T_)/(1 + std::exp(-ε_.val_[i]/T_)) + : Real(1)/(1 + std::exp(ε_.val_[i]/T_)))); + + dεdμ_.δval_ += std::abs(dεdμ_.val_[i] - dεdμ_.val_prev_[i]); + } +} + +bool TBASolutionLiebLiniger::solve_for_dεdμ (Real eps) +{ + dεdμ_ = Function(ε_.interval_); + + IndexU niter { 0 }; + Real prev_δval_ { }; + do { + prev_δval_ = dεdμ_.δval_; + iterate_dεdμ(); + niter++; + } while (dεdμ_.δval_ > eps && dεdμ_.δval_ < prev_δval_); + + return dεdμ_.δval_ < std::sqrt(std::numeric_limits::epsilon()); + } + +void TBASolutionLiebLiniger::populate_ρ_and_ρh () +{ + ρ_ = Function(ε_.interval_); + ρh_ = Function(ε_.interval_); + + for (IndexU i { 0 }; i < ρ_.interval_.λ.size(); ++i) { + ρ_.val_[i] = -(dεdμ_.val_[i]/twopi_r) * + (ε_.val_[i] > Real(0) ? std::exp(-ε_.val_[i]/T_)/(1 + std::exp(-ε_.val_[i]/T_)) + : Real(1)/(1 + std::exp(ε_.val_[i]/T_))); + ρh_.val_[i] = -(dεdμ_.val_[i]/twopi_r) * + (ε_.val_[i] > Real(0) ? Real(1)/(1 + std::exp(-ε_.val_[i]/T_)) + : std::exp(ε_.val_[i]/T_)/(1 + std::exp(ε_.val_[i]/T_))); + } +} + +void TBASolutionLiebLiniger::solve (Real eps) +{ + converge_ε(eps); + solve_for_dεdμ(eps); + populate_ρ_and_ρh(); +} + +Real TBASolutionLiebLiniger::g () +{ + /// Free energy density + /// https://integrability.org/e_l_YY.html#l.g + Real sum { 0 }; + + for (IndexU i { 0 }; i < ε_.val_.size(); ++i) + sum += ε_.interval_.dλ[i] * Tln1pem(ε_.val_[i]); + + return -sum/twopi_r; +} + +Real TBASolutionLiebLiniger::ρ () +{ + // Integrated density + Real sum { 0 }; + for (IndexU i { 0 }; i < ρ_.val_.size(); ++i) sum += ρ_.interval_.dλ[i] * ρ_.val_[i]; + + return sum; +} + +Real TBASolutionLiebLiniger::x (Real λ) +{ + // Counting function x(λ) = (λ + ϕ*ρ(λ))/2π + Real result { λ }; + for (IndexU i { 0 }; i < ρ_.val_.size(); ++i) + result += ρ_.interval_.dλ[i] * ϕ(λ - ρ_.interval_.λ[i]) * ρ_.val_[i]; + return result/twopi_r; +} + + +export TBASolutionLiebLiniger tba_solution_LiebLiniger_at_filling +(Real c, Real T, Real target_filling) +{ + Real μ { Real(-1) }; + Real dμ { Real(0.5) }; + + Real running_eps { 1.0e-4l }; + + TBASolutionLiebLiniger tbasol(c, T, μ); + tbasol.solve(running_eps); + + Real ρ_prev { tbasol.ρ() }; + + Real dμ_prev, dρdμ; + + μ += dμ; + tbasol.reset_μ(μ); + tbasol.solve(running_eps); + + IndexU niter { 0 }; + Real Δρ_prev { }; + + do { + Δρ_prev = std::abs(tbasol.ρ() - target_filling); + dμ_prev = dμ; + dρdμ = (tbasol.ρ() - ρ_prev)/dμ; + ρ_prev = tbasol.ρ(); + dμ = (target_filling - tbasol.ρ())/dρdμ; + // stabilize iterations by avoiding μ changing by more than a factor of 2: + if (std::abs(dμ) > 2*std::abs(dμ_prev)) dμ = 2*dμ * std::abs(dμ_prev/dμ); + μ += dμ; + tbasol.reset_μ(μ); + //running_eps = std::min(running_eps/2, dμ*dμ); // not great + running_eps = running_eps/4; // to improve: make truly adaptive + tbasol.solve(running_eps); + niter++; + // std::cout << "\n*** Filling search, niter " << niter << "\trunning_eps " << running_eps + // << "\tμ " << std::setprecision(20) << μ + // << "\tdμ " << std::setprecision(8) << dμ + // << "\tρ() " << tbasol.ρ() << "\ttarget " << target_filling + // << "\tΔρ " << tbasol.ρ() - target_filling + // << "\n" << std::endl; + } while (// if the density is on target, stop + std::abs(tbasol.ρ() - target_filling) > 100*std::numeric_limits::epsilon() + && ( + // if improving, keep going + std::abs(tbasol.ρ() - target_filling) < Δρ_prev + || + // if chemical potential still too inaccurate, keep going + std::fabs(dμ) > std::sqrt(std::numeric_limits::epsilon()) + ) + ); + return tbasol; +} diff --git a/src/abacus.cc b/src/abacus.cc new file mode 100644 index 0000000..0748956 --- /dev/null +++ b/src/abacus.cc @@ -0,0 +1,23 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module abacus; + +export import conveniences; +export import timer; +export import labels; +export import spaces; +export import model; +export import model_LiebLiniger; +export import quantumnumbers; +export import plex; +export import matrix; +export import state; +export import state_LiebLiniger; +export import matrixelements_LiebLiniger; diff --git a/src/execs/LiebLiniger_typeII_dispersion.cc b/src/execs/LiebLiniger_typeII_dispersion.cc new file mode 100644 index 0000000..265e9ee --- /dev/null +++ b/src/execs/LiebLiniger_typeII_dispersion.cc @@ -0,0 +1,79 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +import std; + + +import abacus; + + +int main(int argc, char* argv[]) +{ + if (argc != 4) { // provide some info + + std::cout << "Welcome to Abacus version " << ABACUS_VERSION + << ", copyright © Jean-Sébastien Caux.\n"; + std::cout << "\nLiebLiniger_typeII_dispersion executable purpose: compute Type II dispersion and curvature for Lieb-Liniger"; + std::cout << "\nUsage: provide the following arguments:\n"; + std::cout << "c_int \t\tValue of the interaction parameter: use positive real values only\n"; + std::cout << "L \t\tLength of the system: use positive real values only\n"; + std::cout << "int N \t\tNumber of particles: use positive integer values only\n"; + std::cout << "\nEXAMPLE:\n\n"; + std::cout << "LiebLiniger_typeII_dispersion 1.0 100.0 100\n\n"; + + return 0; + } + + // correct nr of arguments + int narg { 0 }; + long double c = std::atof(argv[++narg]); + long double L = std::atof(argv[++narg]); + int N = std::atoi(argv[++narg]); + + std::cout << "# Lieb-Liniger with c = " << c << " and L, N = " << L << ", " << N << std::endl; + std::cout << "# Dispersion, velocity and curvature for Type II modes" << std::endl; + std::cout << "# (defined as hole above right Fermi edge shifting left progressively)" << std::endl; + std::cout << "# iK\tω\t\t\tdω/dk\t\t\td^2ω/dk^2" << std::endl; + + std::cout << std::setprecision(std::numeric_limits::digits10 + 1); + + + BosonicContinuum bc1(L); + LiebLinigerModel LL1(bc1, c); + + LiebLinigerBetheState gs(LL1, N); + gs.polish(); + + LiebLinigerBetheState estate(LL1, N); + estate.polish(); + + long double gs_E { gs.E() }; // baseline energy for unexcited state + + std::vector energies; + energies.push_back(estate.E()); + std::vector iKs; + iKs.push_back(estate.iK()); + + for (int α { 0 }; α < N; ++α) { + estate.Ix2_[N-1 - α] += 2; + estate.compute_Momentum(); + estate.polish(); + energies.push_back(estate.E()); + iKs.push_back(estate.iK()); + } + + + for (int α { 1 }; α < N; ++α) { + std::cout << iKs[α] << "\t" << energies[α] - gs_E << "\t" + << (energies[α+1] - energies[α-1])*L/(4*pi_r) << "\t" + << (energies[α+1] -2*energies[α] + energies[α-1])*L*L/fourpisq_r << std::endl; + } + + return 0; +} diff --git a/src/generic/conveniences.cc b/src/generic/conveniences.cc new file mode 100644 index 0000000..0fed3bc --- /dev/null +++ b/src/generic/conveniences.cc @@ -0,0 +1,215 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +/* +Unicode codes (): unused +Ł 0141 +α 03b1 β 03b2 γ 03b3 δ 03b4 ε 03b5 η 03b7 θ 03b8 ι 03b9 κ 03ba λ 03bb μ 03bc ν 03bd ξ 03be π 03c0 ρ 03c1 +σ 03c3 τ 03c4 υ 03c5 φ 03c6 χ 03c7 ψ 03c8 ω 03c9 +Δ 0394 Ξ 039e Σ 03a3 +ƛ 019b +(ϕ 03d5 ɸ 0278 Φ 03a6) +Infinity: ∞ 221e Integral: ∫ 222b +é 00e9 © 00a9 ↑ 2191 ↓ 2193 +*/ + + +/*! + @module conveniences + + @brief Module containing basic definitions and convenient utilities + + This module contains + * mathematical constants + * +*/ + +export module conveniences; + + +import std; + + +export { + + //! Codebase version + constexpr std::string_view ABACUS_VERSION { "2.0.0" }; + + //! Error handling + class AbacusException : public std::runtime_error + { + public: + AbacusException(const std::string& error) + : std::runtime_error{error} + {} + }; + + // Utilities for array/vector indices + using IndexU = std::size_t; //!< unsigned index + using IndexS = std::ptrdiff_t; //!< signed index, arithmetic (including subtraction) is allowed + + //! shortcut conversion from signed to unsigned index (capital Xi, U+039E) + constexpr IndexU Ξ (const IndexS value) { return static_cast(value); } + + //! shortcut conversion from signed to unsigned index (capital Xi, U+039E) + constexpr IndexU Ξ (const int value) { return static_cast(value); } + + //! shortcut conversion from unsigned to signed index (capital Xi, U+039E) + constexpr IndexS Ξ (const IndexU value) { return static_cast(value); } + + //! static cast to int from integral types (lowercase iota, U+03B9) + constexpr int ι (const std::integral auto value) { return static_cast(value); } + // // static cast to int from floating-point types (lowercase iota, U+03B9) + // constexpr int ι (const std::floating_point auto value) { return static_cast(value); } + + //////////////////////////// + // Mathematical conveniences + //////////////////////////// + + // // For array/vector indices: natural numbers (here including zero) + // using Natural = std::uint_fast16_t; // should be unsigned long long int on 64 bit machines + // // For general integral arithmetic: integer numbers + // using Integer = std::int_fast16_t; // should be long long int on 64 bit machines + + // conversion between signed and unsigned integral types (Ξ == capital Xi, U+039E) + // constexpr Natural Ξ (const Integer value) { return static_cast(value); } + // constexpr Integer Ξ (const Natural value) { return static_cast(value); } + + //! Typedef for floating point numbers + // using Real = float; // does not compile + // using Real = double; // can be faster than long double, at the cost of some precision + using Real = long double; + + // constants + const Real pi_r { std::numbers::pi_v }; + const Real twopi_r { Real(2) * pi_r }; + const Real fourpisq_r { twopi_r * twopi_r }; + const Real piover2_r { pi_r/2 }; + const Real oneoverpi_r { Real(1)/pi_r }; + const Real oneovertwopi_r { oneoverpi_r/2 }; + + // Precision-related constants + const Real real_eps { std::numeric_limits::epsilon() }; + const Real ln_real_eps { std::log(std::numeric_limits::epsilon()) }; + const Real sqrt_real_eps { std::sqrt(std::numeric_limits::epsilon()) }; + const Real ln_sqrt_real_eps { std::log(std::sqrt(std::numeric_limits::epsilon())) }; + const Real ten_real_eps { 10*std::numeric_limits::epsilon() }; + constexpr auto full_precision_width { std::numeric_limits::digits10 + 8 }; + const Real double_eps { std::numeric_limits::epsilon() }; + const Real sqrt_double_eps { std::sqrt(std::numeric_limits::epsilon()) }; + const Real ten_double_eps { 10*std::numeric_limits::epsilon() }; + + + // binomial coefficients + Real ln_choose (std::size_t n, std::size_t m) + { + if (n < m) throw std::invalid_argument("n < m " + std::to_string(n) + ", " + + std::to_string(m) + " in ln_choose"); + return std::lgamma(Real(n + 1)) - std::lgamma(Real(m + 1)) - + std::lgamma(Real(n - m + 1)); + } + + unsigned long long int choose_lli (std::size_t n, std::size_t m) + { + Real ln_c { ln_choose(n, m) }; + if (ln_c >= std::log(std::numeric_limits::max())) + throw std::invalid_argument("n, m = " + std::to_string(n) + ", " + std::to_string(m) + + " too high for binomial coefficient"); + return static_cast(std::exp(ln_c) + 0.5); + } + + // complex + using namespace std::complex_literals; + + constexpr std::complex operator""_ir(unsigned long long d) + { + return std::complex { Real(0), static_cast(d) }; + } + + constexpr std::complex operator""_ir(long double d) + { + return std::complex { Real(0), static_cast(d) }; + } + + + // vector + template + std::ostream& operator<< (std::ostream& s, const std::vector& v) { + for (T e : v) s << e << "\t"; + s << std::endl; + return s; + } + + template + IndexU index_in_ordered (T Ix2, const std::vector& Ix2_) + { + // checks whether Ix2 is in Ix2_, which is assumed ordered + // If found, return its index, otherwise sentinel value Ix2_.size() + + if (Ix2_.size() < 1) return Ix2_.size(); + + int index { int(Ix2_.size() - 1)/2 }; // midway in array (left middle if size is even) + int lower { 0 }; + int upper { int(Ix2_.size()) - 1 }; + + do + { + if (Ix2 == Ix2_[Ξ(index)]) return Ξ(index); + + if (upper == lower) return Ix2_.size(); + + if (Ix2 < Ix2_[Ξ(index)]) // between lower and index-1 + upper = index-1; + + else // between index+1 and upper + lower = index+1; + + index = lower + (upper - lower)/2; // midway in what's left (left middle if size is even) + + } while (lower <= upper); + + return Ix2_.size(); + } + + + template + bool is_in_ordered (T Ix2, const std::vector& Ix2_) + { + // checks whether Ix2 is in Ix2_, which is assumed ordered + return (index_in_ordered(Ix2, Ix2_) < Ix2_.size()); + } + + template + IndexU interval_index_in_ordered (T v, const std::vector& vec) + { + // returns i such that vec[i-1] < v <= vec[i], or vec.size() if v > vec.back() + + if (v <= vec[0]) return 0; + if (v > vec.back()) return vec.size(); + + int index = int(vec.size() - 1)/2; // midway in array (left middle if size is even) + int lower { 1 }; + int upper { int(vec.size()) - 1 }; + + do + { + if (v > vec[Ξ(index-1)] && v <= vec[Ξ(index)]) return Ξ(index); + + if (v <= vec[Ξ(index-1)]) // result must be between lower and index-1 + upper = index - 1; + else // here v > vec[index], so result must be between index+1 and upper + lower = index+1; + + index = lower + (upper - lower)/2; // midway in what's left (left middle if size is even) + } while (true); + + return vec.size(); + } + +} diff --git a/src/generic/labels.cc b/src/generic/labels.cc new file mode 100644 index 0000000..d9ee992 --- /dev/null +++ b/src/generic/labels.cc @@ -0,0 +1,542 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +/* +Labels + +Notations: +- we use ~ as the exc separator, _ as the level separator +- l==level, f==filling, p==plexlabel, and using g to denote the ground level + +Definitions: +- plexlabel: label for a given plex, format: [p] == (excs only, can be empty, see userguide) +- levellabel: label for a given level, format: [ll] == [l]~[f]~[p] +- strlabel: label for higher strings, format: [sl] == [ll0]_[ll1]..._[lln] + with an entry for all higher levels with filling > 0 + (higher strings then occupy levels 0,..., str_size-1) +- glabel: label for the ground level, + format: [lg] == [fg]~[pg] if [pg] not empty, + [fg] otherwise +- label: label for a BetheState, format: + label == compress([lg]_[sl]) if [sl] is not empty, + compress([lg]) otherwise + +Remarks: +- the number of excitations in the ground level == count([lg],'~')/2 = (count([pg],'~')+1)/2 +- the number of excitations in a higher level == (count([ll],'~')-1)/2 = (count([p],'~')/2-1 +*/ + + +////////////// +// ↓ Labels // +////////////// + + +export module labels; + + +import std; + + +import conveniences; + + + +export const char LABEL_LEVEL_SEPARATOR { '~' }; +export const char LABEL_EXC_SEPARATOR { '_' }; + +const std::map compress_map { + { '0', 'g' }, { '1', 'h' }, { '2', 'i' }, { '3', 'j' }, + { '4', 'k' }, { '5', 'l' }, { '6', 'm' }, { '7', 'n' }, + { '8', 'o' }, { '9', 'p' }, { 'a', 'q' }, { 'b', 'r' }, + { 'c', 's' }, { 'd', 't' }, { 'e', 'u' }, { 'f', 'v' } +}; + +const std::map inflate_map { + { '~', "~"}, { '_', "_" }, { '0', "0" }, { '1', "1" }, + { '2', "2" }, { '3', "3" }, { '4', "4" }, { '5', "5" }, + { '6', "6" }, { '7', "7" }, { '8', "8" }, { '9', "9" }, + { 'a', "a" }, { 'b', "b" }, { 'c', "c" }, { 'd', "d" }, + { 'e', "e" }, { 'f', "f" }, { 'g', "_0" }, { 'h', "_1" }, + { 'i', "_2" }, { 'j', "_3" }, { 'k', "_4" }, { 'l', "_5" }, + { 'm', "_6" }, { 'n', "_7" }, { 'o', "_8" }, { 'p', "_9" }, + { 'q', "_a" }, { 'r', "_b" }, { 's', "_c" }, { 't', "_d" }, + { 'u', "_e" }, { 'v', "_f" } +}; + +export std::string compress (std::string_view label) +{ + std::string compressed_label; + + IndexU i { 0 }; + while (i < label.size()) { + compressed_label += + (label[i] == LABEL_EXC_SEPARATOR) ? compress_map.at(label[++i]) : label[i]; + i++; + } + return compressed_label; +} + +std::string inflate (std::string_view label) +{ + std::string inflated_label; + for (char c : label) inflated_label += inflate_map.at(c); + return inflated_label; +} + + +/// Given quantum number vectors Lx2 and Ix2, +/// return a plexlabel specifying vectors giving the indices of Lx2 elements not in Ix2, +/// and Ix2 values not in Lx2. +/// Assumed precondition: Lx2 and Ix2 are ordered +/// +export std::string plexlabel (const std::vector& Ix2, const std::vector& Lx2) +{ + if (Lx2.size() != Ix2.size()) + throw AbacusException("Unequal sized vectors in plexlabel(Ix2, Lx2)"); + + IndexU Ii { 0 }; + IndexU Oi { 0 }; + + std::vector hi; ///< hole indices + std::vector pIx2; ///< particle excitation quantum numbers + + while (Ii < Ix2.size() && Oi < Lx2.size()) + { + while (Ii < Ix2.size() && Oi < Lx2.size() && Ix2[Ii] == Lx2[Oi]) { Ii++; Oi++; } + + if (Ii < Ix2.size() && Oi < Lx2.size()) + { + if (Ix2[Ii] < Lx2[Oi]) pIx2.push_back(Ix2[Ii++]); + else hi.push_back(Oi++); + } + if (Ii == Ix2.size()) while (Oi < Lx2.size()) hi.push_back(Oi++); + if (Oi == Lx2.size()) while (Ii < Ix2.size()) pIx2.push_back(Ix2[Ii++]); + } + + IndexU nex { hi.size() }; + std::stringstream plexlabel_; + + if (nex > 0) + { + plexlabel_ << std::hex << hi[0] << LABEL_EXC_SEPARATOR; + + for (IndexU i { 1 }; i < hi.size(); ++i) + plexlabel_ << hi[i] - hi[Ξ(Ξ(i)-1)] - 1 << LABEL_EXC_SEPARATOR; + + // can be positive or negative; even: positive, odd: negative + plexlabel_ << (pIx2[nex-1] >= 0 ? 2*pIx2[nex-1] : -2*pIx2[nex-1]-1); + + for (IndexS i { Ξ(nex-2) }; i >= 0; --i) + plexlabel_ << LABEL_EXC_SEPARATOR << (pIx2[Ξ(i+1)] - pIx2[Ξ(i)])/2 - 1; // always positive + } + + return plexlabel_.str(); +} + +export IndexU count_nex_in_plexlabel (std::string_view plexlabel) +{ + // Let nex represent the number of excitations between Ix2 and Lx2. + // The plexlabel then contains 2nex - 1 separators. + return Ξ(std::count(plexlabel.begin(), plexlabel.end(), + LABEL_EXC_SEPARATOR) + 1)/2; +} + + +/*! + @brief Utility class to manage the conversion from plexlabels + to quantum numbers at a given level. + + Closely related to ParsedLabel, which does this + for all levels, starting from a (full) label. +*/ +class ParsedPlexlabel +{ +public: + // constructors + ParsedPlexlabel () : nex_ { 0 } {} + ParsedPlexlabel (std::string_view plexlabel); + + // data access + IndexU nex () const { return nex_; } + + // utilities + bool is_compatible (const std::vector& Lx2) const; + void set_Ix2 (std::vector& Ix2, const std::vector& Lx2) const; + +public: + IndexU nex_; + std::vector hi_; + std::vector pIx2_; +}; + +ParsedPlexlabel::ParsedPlexlabel (std::string_view plexlabel) + : nex_ { count_nex_in_plexlabel(plexlabel) } + , hi_ { std::vector(nex_) } + , pIx2_ { std::vector(nex_) } +{ + if (nex_ > 0) + { + std::stringstream plexlabel_stream; + plexlabel_stream << plexlabel; + + std::string token; + std::vector tokens; + + while (std::getline(plexlabel_stream, token, LABEL_EXC_SEPARATOR)) + tokens.push_back(std::move(token)); + + hi_[0] = std::stoi(tokens[0], nullptr, 16); + pIx2_[nex_-1] = std::stoi(tokens[nex_], nullptr, 16); + + // p[nex_-1-a] <-> tokens[nex_+a] so index(token) = 2nex_-1 - index(p) + pIx2_[nex_-1] = (pIx2_[nex_-1] % 2 ? -(pIx2_[nex_-1]/2)-1 : pIx2_[nex_-1]/2); + + for (IndexU i { 0 }; i+1 < nex_; ++i) + hi_[i+1] = std::stoi(tokens[i+1], nullptr, 16) + hi_[i] + 1; + + for (int i { int(nex_)-2 }; i >= 0; --i) { + pIx2_[Ξ(i)] = pIx2_[Ξ(i)+1] + - 2*std::stoi(tokens[Ξ(Ξ(2*nex_)-1-i)], nullptr, 16) - 2; + } + } +} + +bool ParsedPlexlabel::is_compatible (const std::vector& Lx2) const +{ + // check validity of excitations, with conditions: + // - hi in [0, Lx2.size()[ + // - hi strictly increasing + // - pIx2 strictly increasing + // - pIx2 not in Lx2 + + // std::cout << "\nIn ParsedPlexlabel: nex_ " << nex_ << "\tLx2.size " << Lx2.size() << "\n"; + if (nex_ == 0) return true; + if (nex_ > Lx2.size()) return false; + + // std::cout << "\nIn ParsedPlexlabel::is_compatible\nhi = "; + // for (IndexU i : hi_) std::cout << i << "\t"; + // std::cout << "\npIx2 = "; + // for (int i : pIx2_) std::cout << i << "\t"; + // std::cout << std::endl; + // std::cout << "\nLx2 "; + // for (int i : Lx2) std::cout << i << "\t"; + // std::cout << std::endl; + + for (IndexU i { 0 }; i < nex_;++i) + if (hi_[i] >= Lx2.size() || + (i+1 < nex_ && + (hi_[i+1] <= hi_[i] || + pIx2_[i+1] <= pIx2_[i])) || + is_in_ordered(pIx2_[i], Lx2)) + return false; + + return true; +} + +void ParsedPlexlabel::set_Ix2 (std::vector& Ix2, const std::vector& Lx2) const +{ + // Given a parsed plexlabel, sets the Ix2 relative to Lx2. + // The incoming value of Ix2 is disregarded. + // parsed is assumed to be compatible. + // Lx2 is assumed to be consistent. + + Ix2 = Lx2; + + for (IndexU i { 0 }; i < nex_ ;++i) + { + Ix2[hi_[i]] = pIx2_[i]; + } + + std::sort(Ix2.begin(), Ix2.end()); +} + + +/*! + @brief Utility class to manage the conversion from labels to quantum numbers. + + Closely related to ParsedPlexlabel, which does this for an individual level. + Note: this will also work to parse a baselabel, which is simply a label with only trivial plexlabels at all levels +*/ +export class ParsedLabel +{ + +public: + // constructors + ParsedLabel (std::string_view label); + + // friendship + friend std::ostream& operator<< (std::ostream& s, const ParsedLabel& parsed); + + +public: + int filling_g_; + std::string plexlabel_g_; + ParsedPlexlabel parsed_plexlabel_g_; + + std::vector level_; // level; only contains entry if filling is > 0 + std::vector filling_; // filling; only contains entry if filling is > 0 + + std::vector plexlabel_; // plexlabel; only contains entry if filling is > 0 + std::vector parsed_plexlabel_; // only contains entry if filling is > 0 + + std::string baselabel_; // label, but remove the plexlabels + std::string patternlabel_; // label, but replace the plexlabels by level's nex +}; + +ParsedLabel::ParsedLabel (std::string_view label) +{ + std::stringstream label_stream; + // label_stream << label; + label_stream << inflate(label); + + std::stringstream baselabel_stream; + baselabel_stream << std::hex; + + std::stringstream patternlabel_stream; + patternlabel_stream << std::hex; + + std::string tmp; + std::vector levellabels; + + while (std::getline(label_stream, tmp, LABEL_LEVEL_SEPARATOR)) + levellabels.push_back(std::move(tmp)); + + // Here, let LABEL_EXC_SEPARATOR == _ + // ground level: levellabel[0] is either [fg]_[pg] if pg nontrivial, or [fg] otherwise + std::stringstream ll_ss { levellabels[0] }; // label at level + std::string f_sh; // filling (as string hex) at level + + std::getline(ll_ss, f_sh, LABEL_EXC_SEPARATOR); // f_sh: ground filling (string hex) + filling_g_ = std::stoi(f_sh, nullptr, 16); + baselabel_stream << filling_g_; // stream maps back to hex + patternlabel_stream << filling_g_; // stream maps back to hex + + std::getline(ll_ss, plexlabel_g_); + parsed_plexlabel_g_ = ParsedPlexlabel(plexlabel_g_); + + if (parsed_plexlabel_g_.nex() > 0) + patternlabel_stream << LABEL_EXC_SEPARATOR << parsed_plexlabel_g_.nex(); + + // higher string levels: + // levellabel is [l]_[f]_[p] if [p] nontrivial, [l]_[f] otherwise + // baselabel is [l]_[f] if [f] is nontrivial, empty otherwise + // patternlabel is [l]_[f]_[nex] if [f] is nontrivial, empty otherwise + + std::string level_sh; // level (as string hex) + int li; // level converted to decimal + int fi; // f converted to decimal + std::string plexlabel_sh; + + for (IndexU llj { 1 }; llj < levellabels.size(); ++llj) + { + // start at 1, 0 was for ground + std::stringstream ll { levellabels[llj] }; + level_sh.clear(); + std::getline(ll, level_sh, LABEL_EXC_SEPARATOR); // level_sh: level (as string hex) + li = std::stoi(level_sh, nullptr, 16); + level_.push_back(li); + + f_sh.clear(); + std::getline(ll, f_sh, LABEL_EXC_SEPARATOR); // f_sh: filling (as string hex) + fi = std::stoi(f_sh, nullptr, 16); // always > 0 + filling_.push_back(fi); + + baselabel_stream << LABEL_LEVEL_SEPARATOR << li + << LABEL_EXC_SEPARATOR << fi; + patternlabel_stream << LABEL_LEVEL_SEPARATOR << li + << LABEL_EXC_SEPARATOR << fi; + + plexlabel_sh.clear(); + if (std::getline(ll, plexlabel_sh)) + { + patternlabel_stream << LABEL_EXC_SEPARATOR << count_nex_in_plexlabel(plexlabel_sh); + plexlabel_.push_back(std::move(plexlabel_sh)); + parsed_plexlabel_.emplace_back(plexlabel_.back()); + } + else + { + plexlabel_.push_back(""); + parsed_plexlabel_.emplace_back(""); + } + } + + baselabel_ = baselabel_stream.str(); + patternlabel_ = patternlabel_stream.str(); +} + + +export std::ostream& operator<< (std::ostream& s, const ParsedLabel& parsed) +{ + s << "Parsed label:\n"; + s << "g: " << parsed.filling_g_ << "\t" << parsed.plexlabel_g_ << "\n"; + + for (IndexU i { 0 }; i < parsed.level_.size(); ++i) + { + s << i << ": " << parsed.filling_[i] << "\t" << parsed.plexlabel_[i] << "\t"; + s << "\n"; + } + + return s; +} + +////////////// +// ↑ Labels // +////////////// + + +///////////////////// +/// Further utilities +///////////////////// + +// baselabel from fillings + +std::string baselabel_from_fillings +(int g_f, const std::vector& str_f) +{ + std::stringstream baselabel_stream; + baselabel_stream << std::hex << g_f; + for (IndexU j { 0 }; j < str_f.size(); ++j) + if (str_f[j] > 0) + baselabel_stream << LABEL_LEVEL_SEPARATOR << j + << LABEL_EXC_SEPARATOR << str_f[j]; + return baselabel_stream.str(); +} + +// partition rapidities into strings + +// forward declaration +void append_possible_baselabels +( + std::vector& possible_baselabels, + const std::vector str_l_, + int g_fa, + const std::vector fa, + int g_f, + std::vector str_f, + int target_Δ_string_weight + ); + +void process_candidate_baselabel +( + std::vector& possible_baselabels, + const std::vector& str_l_, + int g_fa, + const std::vector& fa, + int g_f, + const std::vector& f, + int target_Δ_string_weight +) +{ + std::string candidate_baselabel { baselabel_from_fillings(g_f, f) }; + if (!candidate_baselabel.empty()) { + if (g_fa - g_f == target_Δ_string_weight) { + possible_baselabels.emplace_back(candidate_baselabel); + } + append_possible_baselabels + (possible_baselabels, str_l_, g_fa, fa, g_f, f, target_Δ_string_weight); + } +} + +/// Append possible baselabels +/// +/// principles: +/// - filling modifications can only be made at or above the highest level with Δf \neq 0 +/// - filling modifications can be positive or negative, but all f must be \geq 0 +/// - any nonzero Δf can only change further in the same direction +/// +void append_possible_baselabels +( + std::vector& possible_baselabels, + const std::vector str_l_, ///< str_l_: string lengths (fixed) + int g_fa, ///< g_fa: anchor ground filling, assumed one-strings + const std::vector fa, ///< fa: anchor fillings (fixed) + int g_f, ///< g_f: current ground filling, assumed one-strings + std::vector f, ///< str_f: current fillings + int target_Δ_string_weight ///< target change of weight in higher strings (sum of nr * str_l) + ///< note: the string weight is just g_fa - g_f + ) +{ + // find the highest level with Δf \neq 0: + int idx { int(str_l_.size()) - 1 }; + while (idx >= 0 && f[idx] - fa[idx] == 0) idx--; + + // std::cout << "idx = " << idx << "\n"; + + std::vector f_mod = f; + int g_f_mod; + std::string candidate_baselabel; + + // do a further change (increasing or decreasing filling) at all levels above idx + for (IndexU j { Ξ(idx)+1 }; j < str_l_.size(); ++j) { + + if (g_f >= str_l_[j]) { // deplete ground by adding a string at level j + f_mod = f; + f_mod[j] += 1; + g_f_mod = g_f - str_l_[j]; + process_candidate_baselabel + (possible_baselabels, str_l_, g_fa, fa, g_f_mod, f_mod, target_Δ_string_weight); + } + + if (f[j] >= 1) { // demote a higher string back to ground + f_mod = f; + f_mod[j] -= 1; + g_f_mod = g_f + str_l_[j]; + process_candidate_baselabel + (possible_baselabels, str_l_, g_fa, fa, g_f_mod, f_mod, target_Δ_string_weight); + } + } // for j + + // if idx >= 0, do a further change at this level and call self recursively + if (idx >= 0) { + f_mod = f; + f_mod[idx] += (f[idx]-fa[idx]) > 0 ? 1 : -1; + g_f_mod = g_f - str_l_[idx] * ((f[idx]-fa[idx]) > 0 ? 1 : -1); + if (g_f_mod >= 0 && f_mod[idx] >= 0) { + process_candidate_baselabel + (possible_baselabels, str_l_, g_fa, fa, g_f_mod, f_mod, target_Δ_string_weight); + } + } +} + +export std::vector list_possible_baselabels +( + std::vector possible_baselabels, + const std::vector& str_l_, // str_l_: string lengths (fixed) + int g_fa, // g_fa: anchor ground filling, assumed one-strings + const std::vector& str_fa // fa: anchor fillings (fixed) + ) +{ + // overloaded function to initiate the base descendents listing process + + int string_weight { 0 }; + for (IndexU j { 0 }; j < str_l_.size(); ++j) string_weight += str_l_[j] * str_fa[j]; + + // add other equal-string-weight baselabels + append_possible_baselabels + (possible_baselabels, str_l_, g_fa, str_fa, g_fa, str_fa, 0); + + // and modified string charge ones, moving one step up/down at a time + // (so that the earlier the base descendent is listed, the more similar the ground filling is + for (int Δ_string_weight { 1 }; + // upper limit in next line: deplete all strings, or all ground rapidities + Δ_string_weight <= std::max(string_weight, g_fa); + ++Δ_string_weight + ) { + if (Δ_string_weight <= g_fa) { // up to ground depletion + append_possible_baselabels + (possible_baselabels, str_l_, g_fa, str_fa, g_fa, str_fa, Δ_string_weight); + } + if (Δ_string_weight <= string_weight) { // down to strings depletion + append_possible_baselabels + (possible_baselabels, str_l_, g_fa, str_fa, g_fa, str_fa, -Δ_string_weight); + } + } + return possible_baselabels; +} diff --git a/src/generic/model.cc b/src/generic/model.cc new file mode 100644 index 0000000..0cd9cf7 --- /dev/null +++ b/src/generic/model.cc @@ -0,0 +1,149 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module model; + + +import std; + + +import conveniences; +import labels; +import spaces; + + +/////////////////// +// ↓ Class Model // +/////////////////// + +export template +class Model +{ + +public: + enum Type + { + LiebLiniger, + LiebLiniger_Attractive, + SpinHalf_XXX_AntiFerro, + SpinHalf_XXZ_Axial_AntiFerro, + SpinHalf_XXZ_Planar_AntiFerro, + SpinHalf_XX_AntiFerro, + SpinHalf_XXX_Ferro, + SpinHalf_XXZ_Axial_Ferro, + SpinHalf_XXZ_Planar_Ferro, + SpinHalf_XX_Ferro, + }; + +public: // public interface + + // constructors + Model (TSpace& space, Model::Type type, Real Ł) + : space_ { space } + , type_ { type } + , Ł_ { Ł } + { } + + // Model& operator= (const Model& m); + + // utilities + std::string get_modelname_prefix () const; + std::string get_modelname_prefix (int n) const; + virtual std::string get_filename_prefix () const = 0; + int g_f_from_label (std::string label) const; + +public: // protected: + //const TSpace& space_; + //TSpace& space_; + TSpace space_; + //const Type type_; + Type type_; + //const Real Ł_; // scaled length + Real Ł_; // scaled length + + // kinetic and scattering phase functions + // kinetic phase function for the ground string: + virtual Real θ_ (Real ƛ) const = 0; + virtual Real θinv_ (Real ƛ) const = 0; + virtual Real dθdƛ_ (Real ƛ) const = 0; + + // ground string - ground string scattering phase function: + //virtual Real φ_ (Real ƛ) const = 0; + virtual long double φ_ (long double ƛ) const = 0; + virtual double φ_ (double ƛ) const = 0; + virtual Real dφdƛ_ (Real ƛ) const = 0; + + // kinetic phase functions for higher strings: + // For these 3, use integer length and parity arguments + virtual Real θ_ (int nj, int pj, Real ƛ) const = 0; + virtual Real θinv_ (int nj, int pj, Real ƛ) const = 0; + virtual Real dθdƛ_ (int nj, int pj, Real ƛ) const = 0; + + // ground string - higher string scattering phase functions: + virtual Real φ_ (IndexU k, Real ƛ) const = 0; + virtual Real dφdƛ_ (IndexU k, Real ƛ) const = 0; + // string-string scattering phase functions: + virtual Real φ_ (IndexU j, IndexU k, Real ƛ) const = 0; + virtual Real dφdƛ_ (IndexU j, IndexU k, Real ƛ) const = 0; + + // virtual ~Model () {} +}; + +// template +// Model& Model::operator= (const Model& m) +// { +// if (space_ != m.space_ || type_ != m.type_ || Ł_ != m.Ł_) +// throw "Cannot change Model by assignment"; + +// return *this; +// } + +template +std::string Model::get_modelname_prefix () const +{ + switch (type_) + { + case LiebLiniger: return "LiebLiniger"; + case LiebLiniger_Attractive: return "LiebLiniger-a"; + case SpinHalf_XXX_AntiFerro: return "XXX"; + case SpinHalf_XXZ_Axial_AntiFerro: return "XXZ-a"; + case SpinHalf_XXZ_Planar_AntiFerro: return "XXZ-p"; + case SpinHalf_XX_AntiFerro: return "XX"; + case SpinHalf_XXX_Ferro: return "XXX-f"; + case SpinHalf_XXZ_Axial_Ferro: return "XXZ-a-f"; + case SpinHalf_XXZ_Planar_Ferro: return "XXZ-p-f"; + case SpinHalf_XX_Ferro: return "XX-f"; + default: return "undefined"; + } +} + +template +std::string Model::get_modelname_prefix (int n) const +{ + return get_modelname_prefix(static_cast::Type>(n)); +} + +template +int Model::g_f_from_label (std::string label) const +{ + ParsedLabel parsed(label); + return parsed.filling_g_; +} + + + + +// Template specialization concept +export template +concept DerivedModel = std::is_base_of, TModel>::value; + + +/////////////////// +// ↑ Class Model // +/////////////////// diff --git a/src/generic/plex.cc b/src/generic/plex.cc new file mode 100644 index 0000000..0c38a37 --- /dev/null +++ b/src/generic/plex.cc @@ -0,0 +1,198 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module plex; + + +import std; + + +import conveniences; +import labels; +import quantumnumbers; + + +////////////////// +// ↓ Class Plex // +////////////////// + + +export class Plex : public QuantumNumbers { // quantum numbers and rapidities at a given base level + +public: // public interface + + // constructors + Plex () = default; + Plex (int f); + Plex (int l, int p, int f); + Plex (int l, int p, const std::vector& Ox2, int Δf); + // virtual ~Plex () = default; + + // // operators + // Plex& operator= (const Plex& rhs); + + // data access + void print() const; + + // manipulation + bool excite (IndexU α, int δI); + void boost (int δI); + + // // friendship + // template TModel> + // friend class BetheState; + // friend class XXXBetheState; + // friend class XXZAxialBetheState; + // friend class XXZPlanarBetheState; + + // friend std::complex V_ρ + // ( + // IndexU α, + // const LiebLinigerBetheState& bra, + // const LiebLinigerBetheState& ket + // ); + + // friend std::complex matrix_element_ρ + // ( + // const LiebLinigerBetheState& bra, + // const LiebLinigerBetheState& ket + // ); + + // friend std::complex V_ψ + // ( + // IndexU α, + // const LiebLinigerBetheState& bra, + // const LiebLinigerBetheState& ket + // ); + + // friend std::complex matrix_element_ψ + // ( + // const LiebLinigerBetheState& bra, + // const LiebLinigerBetheState& ket + // ); + + // friend std::complex matrix_element_ψdag + // ( + // const LiebLinigerBetheState& bra, + // const LiebLinigerBetheState& ket + // ); + + +public: // protected: + int l_; // string length + int p_; // parity + std::vector ƛ_; // (rescaled) rapidities + std::vector δƛ_; // latest iterative change + std::vector B_; // Bethe function + + void shift_ƛ_with_δƛ_(); + +private: + int check_Δf(const std::vector& Ox2, int Δf) + { + if (int(Ox2.size()) + Δf < 0) throw "Ox2 size + Δf < 0 in QuantumNumbers constructor"; + return Δf; + } +}; + +// constructors +Plex::Plex (int f) : Plex(1, 1, f) {} + +Plex::Plex (int l, int p, int f) + : QuantumNumbers(f) + , l_ { l }, p_ { p } + , ƛ_ { std::vector(f) } + , δƛ_ { std::vector(f) } + , B_ { std::vector(f) } +{} + +Plex::Plex (int l, int p, const std::vector& Ox2, int Δf) + : QuantumNumbers(Ox2, check_Δf(Ox2, Δf)) + , l_ { l }, p_ { p } + , ƛ_ { std::vector(int(Ox2.size()) + Δf) } + , δƛ_ { std::vector(int(Ox2.size()) + Δf) } + , B_ { std::vector(int(Ox2.size()) + Δf) } +{} + +// // operators +// Plex& Plex::operator= (const Plex& rhs) +// { +// if (this == &rhs) return *this; + +// l_ = rhs.l_; +// p_ = rhs.p_; +// ƛ_ = rhs. ƛ_ ; +// δƛ_ = rhs.δƛ_; +// B_ = rhs.B_; // Bethe function + +// return *this; +// } + +// data access +void Plex::print () const +{ + std::cout << l_ << "\t" << p_ << "\t" << f_ << "\n"; + for (int Ix2 : Ix2_) std::cout << Ix2 << "\t"; + std::cout << "\n"; + for (Real ƛ : ƛ_) std::cout << ƛ << "\t"; + std::cout << "\n"; + for (Real δƛ : δƛ_) std::cout << δƛ << "\t"; + std::cout << "\n"; + for (Real B : B_) std::cout << B << "\t"; + std::cout << "\n"; +} + +// manipulation +bool Plex::excite (IndexU α, int δI) { + // shifts the ground Ix2_[α] by δIx2 if this is allowed + if ( + α < Ix2_.size() && + Ix2_[α] + δI*2 >= Ix2_min_ && Ix2_[α] + δI*2 <= Ix2_max_ + && !is_in_ordered(Ix2_[α] + δI*2, Ix2_) + ) { + Ix2_[α] += δI*2; + std::sort(Ix2_.begin(), Ix2_.end()); + set_plexlabel_from_Ix2_(); + return true; + } + return false; +} + +void Plex::boost (int δI) { + // To accelerate fixed-momentum scans, this boosts both Ix2_ and Ox2_ + // to the momentum value which is required. + if (δI > 0 && + Ix2_.back() + δI*2 >= Ix2_min_ + && Ix2_.back() + δI*2 <= Ix2_max_ + ) { + // std::cout << " boosting " << Ix2_.back(); + Ix2_.back() += δI*2; + Ox2_.back() += δI*2; + // std::cout << " to " << Ix2_.back() << "\n"; + set_plexlabel_from_Ix2_(); + } + if (δI < 0 && + Ix2_.front() + δI*2 >= Ix2_min_ + && Ix2_.front() + δI*2 <= Ix2_max_ + ) { + // std::cout << " boosting " << Ix2_.back(); + Ix2_.front() += δI*2; + Ox2_.front() += δI*2; + // std::cout << " to " << Ix2_.back() << "\n"; + set_plexlabel_from_Ix2_(); + } +} + +void Plex::shift_ƛ_with_δƛ_() { + for (IndexU α { 0 }; α < ƛ_.size(); ++α) ƛ_[α] += δƛ_[α]; +} + +////////////////// +// ↑ Class Plex // +////////////////// diff --git a/src/generic/quantumnumbers.cc b/src/generic/quantumnumbers.cc new file mode 100644 index 0000000..cffa4a8 --- /dev/null +++ b/src/generic/quantumnumbers.cc @@ -0,0 +1,514 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module quantumnumbers; + + +import std; + + +import conveniences; +import labels; +// import descendents; + + + +//////////////////////////// +// ↓ Class QuantumNumbers // +//////////////////////////// + +export class QuantumNumbers // quantum numbers for an individual level +{ + +public: + // constructors + QuantumNumbers (); + QuantumNumbers (int f); + QuantumNumbers (const std::vector& Ox2); + QuantumNumbers (const std::vector& Ox2, int Δf); // filling adapted + // virtual ~QuantumNumbers () = default; + + // manipulation + bool set_Ix2 (const std::vector& Ix2_req); + + // checks + bool is_symmetrical () const; + bool zero_occupied () const; + + // void append_descendent_plexlabels + // (DescendentType type, std::vector& list); + +public: // protected: + std::string plexlabel_; + int f_; // filling, int (and not IndexU) due to many uses in arithmetic + int Ix2_min_; + int Ix2_max_; + Real ln_dim_; // ln of dimensionality of this subspace + std::vector Ix2_; // (doubled) quantum numbers + std::vector Ox2_; // origin (doubled) quantum numbers (may be boosted) + std::vector Lx2_; // origin (doubled) quantum numbers (unboosted, for labelling) + + void set_Ix2_limits_(int Ix2_min, int Ix2_max); + void set_plexlabel_from_Ix2_(); + // IndexU lowest_excitation_ (Chirality chirality); + // std::variant> descendents + // (DescendentType type); + // std::string descendent_a_ (Chirality chirality); + // std::string descendent_b_leading_ (Chirality chirality); + // std::vector descendents_b_subleading_ (Chirality chirality); + // std::string descendent_d_leading_ (Chirality chirality); + // std::vector descendents_d_subleading_ (Chirality chirality); + +private: + int check_filling (const std::vector& Ox2, int Δf) + { + if (int(Ox2.size()) + Δf < 0) throw "Negative filling in QuantumNumbers"; + return IndexU(int(Ox2.size()) + Δf); + } +}; + +QuantumNumbers::QuantumNumbers () + : QuantumNumbers(0) {} + +QuantumNumbers::QuantumNumbers (int f) + : f_ { f } + , Ix2_min_ { std::numeric_limits::min() } + , Ix2_max_ { std::numeric_limits::max() } + , ln_dim_ { f_ > 0 ? std::numeric_limits::infinity() : Real(0) } + , Ix2_ { std::vector(f) } + , Ox2_ { std::vector(f) } + , Lx2_ { std::vector(f) } +{ + for (int a { 0 }; a < f; ++a) { + Ix2_[Ξ(a)] = -(f - 1) + 2*a; + Ox2_[Ξ(a)] = -(f - 1) + 2*a; + } + Lx2_ = Ox2_; + set_plexlabel_from_Ix2_(); +} + +QuantumNumbers::QuantumNumbers (const std::vector& Ox2) + : QuantumNumbers(Ox2, 0) {} + +QuantumNumbers::QuantumNumbers (const std::vector& Ox2, int Δf) // filling adapted + : f_ { check_filling(Ox2, Δf) } + , Ix2_min_ { std::numeric_limits::min() } + , Ix2_max_ { std::numeric_limits::max() } + , ln_dim_ { f_ > 0 ? std::numeric_limits::infinity() : Real(0) } + , Ix2_ { std::vector(f_) } + , Ox2_ { std::vector(f_) } + , Lx2_ { std::vector(f_) } +{ + // Set Ix2_ and Ox2_ to filling-adapted Ox2, adding/removing from the middle + + if (Δf == 0) + { + // Ix2_ = Ox2; // 2025-01-31 + Ox2_ = Ox2; + } + else if (Δf > 0) + { // we add Δf particles in the middle + // shifting by Δf gives the correct evenness of the quantum numbers + + if (Ox2.size() == 0) + { + for (IndexU a { 0 }; a < Ξ(Δf); ++a) Ox2_[a] = -Δf + 1 + 2*ι(a); + } + else + { + for (IndexU a { 0 }; a <= Ox2.size()/2; ++a) Ox2_[a] = Ox2[a] - Δf; + + for (IndexU a { 1 }; a <= Ξ(Δf); ++a) + Ox2_[Ox2.size()/2 + a] = ι(Ox2_[Ox2.size()/2] + 2*a); + + for (IndexU a { Ox2.size()/2 + 1 }; a < Ox2.size(); ++a) + Ox2_[a + Ξ(Δf)] = Ox2[a] + Δf; + } + } + else + { // Δf < 0, we remove |Δf| particles from the middle, + // and shift towards the middle (remembering that Δf < 0) + if (Ox2_.size() > 0) + { + for (IndexU a { 0 }; a < Ox2_.size(); ++a) + { + Ox2_[a] = Ox2[a + (a > Ox2.size()/2 ? Ξ(-Δf) : 0)] + + (a > Ox2.size()/2 ? Δf : -Δf); + } + } + } + // finally, set the Ix2_ to the obtained Ox2_ + // for (IndexU a { 0 }; a < Ox2_.size(); ++a) Ix2_[a] = Ox2_[a]; + Ix2_ = Ox2_; + Lx2_ = Ox2_; + set_plexlabel_from_Ix2_(); +} + + +bool QuantumNumbers::set_Ix2 (const std::vector& Ix2_req) +{ + // check vector length + if (Ix2_req.size() != Ξ(f_)) return false; + + // check parity of each required quantum nr + // Ix2 + f_ must be odd + for (int Ix2_given : Ix2_req) if (!((Ix2_given + f_) % 2)) return false; + + Ix2_ = Ix2_req; + set_plexlabel_from_Ix2_(); // 2025-01-31 + return true; +} + +bool QuantumNumbers::is_symmetrical () const { + for (IndexU α { 0 }; α < (Ix2_.size() + 1)/2; ++α) + if (Ix2_[α] != -Ix2_[Ξ(Ix2_.size() - 1)]) return false; + return true; +} + +bool QuantumNumbers::zero_occupied () const { + return is_in_ordered(0, Ix2_); +} + +void QuantumNumbers::set_Ix2_limits_(int Ix2_min, int Ix2_max) +{ + if (Ix2_min > Ix2_max) + { + std::cout << Ix2_min << "\t" << Ix2_max << std::endl; + throw std::invalid_argument("Improper Ix2_min,max in QuantumNumbers::set_Ix2_limits_"); + } + + if (Ix2_max - Ix2_min + 2 < f_) + throw std::invalid_argument("Insufficient limits in QuantumNumbers::set_Ix2_limits_"); + + Ix2_min_ = Ix2_min; + Ix2_max_ = Ix2_max; + + // Ensure proper parity (by narrowing): if filling is even(odd), Ix2 must odd(even) + if (!((Ix2_min_ + f_) % 2)) Ix2_min_++; + if (!((Ix2_max_ + f_) % 2)) Ix2_max_--; + + if (f_ > 0) + { + ln_dim_ = std::lgamma(Real((Ix2_max_ - Ix2_min_)/2 + 2)) - + std::lgamma(Real((Ix2_max_ - Ix2_min_)/2 + 2 - f_)) - std::lgamma(Real(f_+1)); + } + else + { + ln_dim_ = Real(0); + } + if (std::isnan(ln_dim_)) ln_dim_ = std::numeric_limits::infinity(); +} + +void QuantumNumbers::set_plexlabel_from_Ix2_() +{ + // plexlabel_ = plexlabel(Ix2_, Ox2_); + plexlabel_ = plexlabel(Ix2_, Lx2_); +} + +// IndexU QuantumNumbers::lowest_excitation_ (Chirality chirality) +// { +// if (chirality == Chirality::right) // look for leftmost R moving +// { +// for (IndexU α { 0 }; α < Ix2_.size(); ++α) +// if (Ix2_[α] > Ox2_[α]) return α; +// } +// else // look for rightmost L moving +// { +// for (int a { int(Ix2_.size()) - 1}; a >= 0; --a) +// if (Ix2_[Ξ(a)] < Ox2_[Ξ(a)]) return Ξ(a); +// } +// return Ix2_.size(); // not found, sentinel value +// } + +// void QuantumNumbers::append_descendent_plexlabels +// (DescendentType type, std::vector& list) +// { +// Move move { std::get(type) }; +// Chirality chirality { std::get(type) }; + +// std::string tmp; + +// switch (move) +// { +// case Move::a_lead: +// tmp = descendent_a_(chirality); +// if (!tmp.empty()) list.push_back(tmp); +// break; + +// case Move::b_lead: +// tmp = descendent_b_leading_(chirality); +// if (!tmp.empty()) list.push_back(tmp); +// break; + +// case Move::b_sub: +// for (std::string label : descendents_b_subleading_(chirality)) +// if (!label.empty()) list.push_back(label); +// break; + +// case Move::d_lead: +// tmp = descendent_d_leading_(chirality); +// if (!tmp.empty()) list.push_back(tmp); +// break; + +// case Move::d_sub: +// for (std::string label : descendents_d_subleading_(chirality)) +// if (!label.empty()) list.push_back(label); +// break; +// } +// } + +// std::string QuantumNumbers::descendent_a_ (Chirality chirality) +// { +// // returns the plexlabel of the a-type descendent, if it exists, +// // otherwise empty string +// std::string descendent_plexlabel; + +// if (Ix2_.size() == 0) return descendent_plexlabel; // no leading possible + +// IndexU α { lowest_excitation_(chirality) }; + +// if (α == Ix2_.size()) return descendent_plexlabel; // no leading found + +// if (chirality == Chirality::right) // a^+ move on leading +// { +// if (!is_in_ordered(Ix2_[α], Ox2_) +// && !is_in_ordered(Ix2_[α] + 2, Ox2_) +// && ((α + 1 == Ix2_.size()) || (Ix2_[α] + 2 != Ix2_[α+1])) +// && Ix2_[α] + 2 <= Ix2_max_ +// ) +// { +// Ix2_[α] += 2; +// // descendent_plexlabel = plexlabel(Ix2_, Ox2_); +// descendent_plexlabel = plexlabel(Ix2_, Lx2_); +// Ix2_[α] -= 2; +// } +// } +// else // a^- move on leading +// { +// if (!is_in_ordered(Ix2_[α] - 2, Ox2_) +// && !is_in_ordered(Ix2_[α], Ox2_) +// && (α == 0 || (Ix2_[α] - 2 != Ix2_[Ξ(Ξ(α) - 1)])) +// && Ix2_[α] - 2 >= Ix2_min_ +// ) +// { +// Ix2_[α] -= 2; +// // descendent_plexlabel = plexlabel(Ix2_, Ox2_); +// descendent_plexlabel = plexlabel(Ix2_, Lx2_); +// Ix2_[α] += 2; +// } +// } +// return descendent_plexlabel; +// } + +// std::string QuantumNumbers::descendent_b_leading_ (Chirality chirality) +// { +// // returns the plexlabel of the leading b-type descendent, if it exists, +// // otherwise empty string +// std::string descendent_plexlabel; + +// if (Ix2_.size() == 0) return descendent_plexlabel; // no leading possible + +// IndexU α { lowest_excitation_(chirality) }; + +// if (α == Ix2_.size()) return descendent_plexlabel; // no leading found + +// if (chirality == Chirality::right) // b^+ move on leading +// { +// if (is_in_ordered(Ix2_[α], Ox2_) +// && !is_in_ordered(Ix2_[α] + 2, Ox2_) +// && ((α+1 == Ix2_.size()) || (Ix2_[α] + 2 != Ix2_[α+1])) +// && (Ix2_[α] + 2 <= Ix2_max_) +// ) +// { +// Ix2_[α] += 2; +// // descendent_plexlabel = plexlabel(Ix2_, Ox2_); +// descendent_plexlabel = plexlabel(Ix2_, Lx2_); +// Ix2_[α] -= 2; +// } +// } +// else // b^- move on leading +// { +// if (!is_in_ordered(Ix2_[α] - 2, Ox2_) +// && is_in_ordered(Ix2_[α], Ox2_) +// && (α == 0 || (Ix2_[α] - 2 != Ix2_[Ξ(Ξ(α)-1)])) +// && (Ix2_[α] - 2 >= Ix2_min_) +// ) +// { +// Ix2_[α] -= 2; +// // descendent_plexlabel = plexlabel(Ix2_, Ox2_); +// descendent_plexlabel = plexlabel(Ix2_, Lx2_); +// Ix2_[α] += 2; +// } +// } +// return descendent_plexlabel; +// } + +// std::vector QuantumNumbers::descendents_b_subleading_ +// (Chirality chirality) +// { +// // return a vector of plexlabels for subleading b-type descendents. +// std::vector descendent_plexlabels; + +// if (Ix2_.size() == 0) // no subleading possible, return empty +// return descendent_plexlabels; + +// IndexU α { lowest_excitation_(chirality) }; + +// if (chirality == Chirality::right) +// { // b^+ move sought on subleading, left of leading +// for (IndexU β { 0 }; β < α; ++β) +// { +// if (Ix2_[β] == Ox2_[β] // unmoved up to now +// && !is_in_ordered(Ix2_[β] + 2, Ox2_) +// && ((β+1 == Ix2_.size()) || (Ix2_[β] + 2 != Ix2_[β+1])) +// && Ix2_[β] + 2 <= Ix2_max_ +// ) +// { +// Ix2_[β] += 2; +// // descendent_plexlabels.push_back(plexlabel(Ix2_, Ox2_)); +// descendent_plexlabels.push_back(plexlabel(Ix2_, Lx2_)); +// Ix2_[β] -= 2; +// } +// } +// } +// else // b^- move sought on subleading, right of leading +// { +// for (IndexU β { α < Ix2_.size() ? α+1 : 0 }; β < Ix2_.size(); ++β) +// { +// if (Ix2_[β] == Ox2_[β] // unmoved up to now +// && !is_in_ordered(Ix2_[β] - 2, Ox2_) +// && (β == 0 || (Ix2_[β] - 2 != Ix2_[Ξ(Ξ(β)-1)])) +// && Ix2_[β] - 2 >= Ix2_min_ +// ) +// { +// Ix2_[β] -= 2; +// // descendent_plexlabels.push_back(plexlabel(Ix2_, Ox2_)); +// descendent_plexlabels.push_back(plexlabel(Ix2_, Lx2_)); +// Ix2_[β] += 2; +// } +// } +// } +// return descendent_plexlabels; +// } + +// std::string QuantumNumbers::descendent_d_leading_ (Chirality chirality) +// { +// // returns the plexlabel of the leading d-type descendent, if it exists, +// // otherwise empty string +// IndexU α { lowest_excitation_(chirality) }; +// std::string descendent_plexlabel; + +// if (α == Ix2_.size()) // no leading found, return empty +// return descendent_plexlabel; + +// if (chirality == Chirality::right) // d^+ move on leading +// { +// // shift right while both Ox2 and Ix2 are empty to the right: +// // Oα is the index of the leading quantum number in Ox2_ +// IndexU Oα { index_in_ordered(Ix2_[α], Ox2_) }; +// if (Oα + 1 < Ox2_.size()) // leading Ix2 found in Ox2_, +// { +// // and there exists at least one other Ox2 to the right +// // Size of empty block in d^{+,n} is n = (Ox2_[Oα + 1] - Ox2_[Oα])/2 - 1 +// if (α+1 == Ix2_.size() || // rightmost or +// // empty block in Ix2_ and empty target position +// (Ix2_[α+1] > Ix2_[α] + Ox2_[Oα + 1] - Ox2_[Oα]) +// ) +// { +// // no need to check for Ix2_max since it's occupied in Ox2_ and thus within limits +// Ix2_[α] += Ox2_[Oα + 1] - Ox2_[Oα]; +// // descendent_plexlabel = plexlabel(Ix2_, Ox2_); +// descendent_plexlabel = plexlabel(Ix2_, Lx2_); +// Ix2_[α] -= Ox2_[Oα + 1] - Ox2_[Oα]; +// } +// } +// } +// else // d^- move on leading +// { +// // shift left while both Ox2 and Ix2 are empty to the left: +// // Oα is the index of the leading quantum number in Ox2_ +// IndexU Oα { index_in_ordered(Ix2_[α], Ox2_) }; +// if (Oα > 0 && Oα < Ox2_.size()) // leading Ix2 found in Ox2_, +// { +// // and there exists at least one other Ox2 to the left +// // Size of empty block in d^{+,n} is n = (Ox2_[Oα] - Ox2_[Oα-1])/2 - 1 +// if (α == 0 || +// // empty block in Ix2_ and empty target position +// (Ix2_[Ξ(Ξ(α)-1)] < Ix2_[α] - Ox2_[Oα] + Ox2_[Ξ(Ξ(Oα)-1)]) +// ) +// { +// Ix2_[α] -= Ox2_[Oα] - Ox2_[Ξ(Ξ(Oα)-1)]; +// // descendent_plexlabel = plexlabel(Ix2_, Ox2_); +// descendent_plexlabel = plexlabel(Ix2_, Lx2_); +// Ix2_[α] += Ox2_[Oα] - Ox2_[Ξ(Ξ(Oα)-1)]; +// } +// } +// } +// return descendent_plexlabel; +// } + +// std::vector QuantumNumbers::descendents_d_subleading_ +// (Chirality chirality) +// { +// // returns a vector of plexlabels for subleading d-type descendent. +// std::vector descendent_plexlabels; + +// if (Ix2_.size() == 0) // no subleading possible, return empty +// return descendent_plexlabels; + +// IndexU α { lowest_excitation_(chirality) }; + +// if (chirality == Chirality::right) +// { // d^+ move sought on subleading, left of leading +// for (IndexU β { 0 }; β < α; ++β) +// { +// if (Ix2_[β] == Ox2_[β] // undisplaced +// && (β + 1 < Ox2_.size())) +// { // there exists at least one other Ox2 to the right +// // Size of empty block in d^{+,n} is n = (Ox2_[Oβ + 1] - Ox2_[Oβ])/2 - 1 +// if (Ix2_[β+1] > Ox2_[β + 1] // empty block in Ix2_ and empty target position +// ) +// { +// // no need to check for Ix2_max since it's occupied in Ox2_ +// // and thus within limits +// Ix2_[β] += Ox2_[β + 1] - Ox2_[β]; +// // descendent_plexlabels.push_back(plexlabel(Ix2_, Ox2_)); +// descendent_plexlabels.push_back(plexlabel(Ix2_, Lx2_)); +// Ix2_[β] -= Ox2_[β + 1] - Ox2_[β]; +// } +// } +// } +// } +// else +// { // d^- move sought on subleading, subleading sought to the right of leading +// for (IndexU β { α < Ix2_.size() ? α+1 : 0 }; β < Ix2_.size(); ++β) +// { +// if (Ix2_[β] == Ox2_[β] // undisplaced +// && β > 0) // there exists at least one other Ox2 to the left +// { +// // Size of empty block in d^{+,n} is n = (Ox2_[Oβ] - Ox2_[Oβ-1])/2 - 1 +// // empty block in Ix2_ and empty target position +// if (Ix2_[Ξ(Ξ(β)-1)] < Ox2_[Ξ(Ξ(β)-1)] +// ) +// { +// Ix2_[β] -= Ox2_[β] - Ox2_[Ξ(Ξ(β)-1)]; +// // descendent_plexlabels.push_back(plexlabel(Ix2_, Ox2_)); +// descendent_plexlabels.push_back(plexlabel(Ix2_, Lx2_)); +// Ix2_[β] += Ox2_[β] - Ox2_[Ξ(Ξ(β)-1)]; +// } +// } +// } +// } +// return descendent_plexlabels; +// } + + +//////////////////////////// +// ↑ Class QuantumNumbers // +//////////////////////////// diff --git a/src/generic/spaces.cc b/src/generic/spaces.cc new file mode 100644 index 0000000..3118ae9 --- /dev/null +++ b/src/generic/spaces.cc @@ -0,0 +1,232 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module spaces; + + +import std; + + +import conveniences; + + + +/////////////////////// +// ↓ Class FockSpace // +/////////////////////// + +/// Abstract base class for derived physical Fock spaces +/// + +export class FockSpace +{ + +public: + /** + Explicit list of supported Fock spaces + */ + enum Type { + BosonicContinuum, ///< For use with the Lieb-Liniger model + SpinHalfChain, ///< For use with Heisenberg chains + }; + +public: // public interface + + // constructors + FockSpace (FockSpace::Type type) : type_ { type } {} + + // physical properties + virtual Real ln_dim () const = 0; + virtual Real ln_subspace_dim_at_filling (int filling) const = 0; + +protected: + Type type_; +}; + +/////////////////////// +// ↑ Class FockSpace // +/////////////////////// + + + +export template +concept DerivedFockSpace = std::is_base_of::value; + + + +////////////////////////////// +// ↓ Class BosonicContinuum // +////////////////////////////// + +export class BosonicContinuum : public FockSpace +{ + +public: + enum Operator { psi, rho, psidag, }; + +public: // public interface + + // constructors + BosonicContinuum (Real L) + : FockSpace(FockSpace::Type::BosonicContinuum) + , L_ { verify_L(L) } + {} + + // physical properties + Real L () const { return L_; } + Real ln_dim () const override + { + return std::numeric_limits::infinity(); + } + Real ln_subspace_dim_at_filling ([[maybe_unused]]int filling) const override + { + return std::numeric_limits::infinity(); + } + + // utilities + bool operator== (const BosonicContinuum& space) const; + +protected: + Real L_; + +private: + static Real verify_L (Real L) + { + if (L < 0.0L) + { + throw std::invalid_argument("L must be greater than 0 in BosonicContinuum"); + } + return L; + } +}; + +bool BosonicContinuum::operator== (const BosonicContinuum& space) const +{ + return (type_ == space.type_ && L_ == space.L_); +} + + +constexpr std::string_view get_operator_name (BosonicContinuum::Operator op) +{ + switch (op) + { + case BosonicContinuum::Operator::psi: return "psi"; + case BosonicContinuum::Operator::rho: return "rho"; + case BosonicContinuum::Operator::psidag: return "psidag"; + default: + throw std::invalid_argument + ("Unknown BosonicContinuum::Operator op in get_operator_name"); + } +} + +export constexpr int Δf (BosonicContinuum::Operator op) +{ + switch (op) + { + case BosonicContinuum::Operator::psi: return -1; + case BosonicContinuum::Operator::rho: return 0; + case BosonicContinuum::Operator::psidag: return 1; + default: + throw std::invalid_argument + ("Unknown BosonicContinuum::Operator op in Δf"); + } +} + +export std::ostream& operator<< (std::ostream& s, BosonicContinuum::Operator op) +{ + return s << get_operator_name(op); +} + +////////////////////////////// +// ↑ Class BosonicContinuum // +////////////////////////////// + + + +/////////////////////////// +// ↓ Class SpinHalfChain // +/////////////////////////// + +export class SpinHalfChain : public FockSpace +{ + +public: + enum Operator { Sm, Sz, Sp, }; + +public: // public interface + + // constructors + SpinHalfChain (int N) + : FockSpace(FockSpace::Type::SpinHalfChain) + , N_ { verify_N(N) } + {} + + // physical properties + int N () const { return N_; } + Real ln_dim () const override { return N_ * std::logl(Real(2)); }; + Real ln_subspace_dim_at_filling (int filling) const override + { + return std::lgamma(Real(N_) + 1) + - std::lgamma(Real(N_) + filling + 1) - std::lgamma(filling + 1); + } + + // utilities + bool operator== (const SpinHalfChain& space) const + { + return (type_ == space.type_ && N_ == space.N_); + } + +protected: + int N_; + +private: + static int verify_N (const int N) + { + if (N < 4 || N%2) + { + throw std::invalid_argument("N must be even and >= 4 in HeisenbergModel"); + } + return N; + } +}; + +constexpr std::string_view get_operator_name (SpinHalfChain::Operator op) +{ + switch (op) + { + case SpinHalfChain::Operator::Sm: return "Sm"; + case SpinHalfChain::Operator::Sz: return "Sz"; + case SpinHalfChain::Operator::Sp: return "Sp"; + default: + throw std::invalid_argument + ("Unknown SpinHalfChain::Operator op in get_operator_name"); + } +} + +export constexpr int Δf (SpinHalfChain::Operator op) +{ + switch (op) + { + case SpinHalfChain::Operator::Sm: return 1; + case SpinHalfChain::Operator::Sz: return 0; + case SpinHalfChain::Operator::Sp: return -1; + default: + throw std::invalid_argument + ("Unknown SpinHalfChain::Operator op in get_operator_name"); + } +} + +export std::ostream& operator<< (std::ostream& s, SpinHalfChain::Operator op) +{ + return s << get_operator_name(op); +} + +/////////////////////////// +// ↑ Class SpinHalfChain // +/////////////////////////// diff --git a/src/generic/state.cc b/src/generic/state.cc new file mode 100644 index 0000000..88ec3ce --- /dev/null +++ b/src/generic/state.cc @@ -0,0 +1,877 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module state; + + +import std; + + +import conveniences; +import timer; +import labels; +import spaces; +import model; +// import descendents; +import plex; +import matrix; + + +export enum IterationMethod +{ + ouroboros, + diagonal, + tridiagonal, + pentadiagonal, + Newton, + maxIterationMethod, +}; + +constexpr std::string_view getIterationMethodName (IterationMethod method) +{ + switch (method) + { + case ouroboros: return "ouroboros"; + case diagonal: return "diagonal"; + case tridiagonal: return "tridiagonal"; + case pentadiagonal: return "pentadiagonal"; + case Newton: return "Newton"; + default: return "undefined"; + } +} + +constexpr std::string_view getIterationMethodName (int nr) +{ + return getIterationMethodName(static_cast(nr)); +} + + + +//////////////////////// +// ↓ Class BetheState // +//////////////////////// + +export template TModel> +class BetheState : public Plex { + + // Notes + // - ground strings are kept separate from higher strings + // (to ease Lieb-Liniger, and accelerate others) + +public: // public data interface + + // constructors, from scratch + BetheState (TModel& model, int g_f); + + // constructors based on a preexisting state, adapted in filling + BetheState (const BetheState& refstate, int Δf, bool relative_label); + + // Identification helpers, for use in filename + std::vector> tags_; // e.g. T_0 + + // accessors + std::vector λ; // ground rapidities + // virtual Real ƛ (IndexU α) const; // accessor for rescaled rapidity + + // subspace dimensionality + Real ln_subspace_dim_this_base () const; + Real ln_subspace_dim_this_pattern () const; + Real ln_subspace_dim_for_pattern (std::string patternlabel) const; + + // manipulation + bool set_g_Ix2 (const std::vector& Ix2_req); + bool excite_g (IndexU α, int δI); + void boost (int δI); + // std::vector descendents (DescendentType type); + // // base descendents + // IndexU nr_base_descendents () { return 0; } + // std::vector base_descendent_baselabels () { return {}; } + + // checks + bool is_symmetrical () const { return Plex::is_symmetrical(); } + + // solve Bethe equations and compute properties + bool draft (std::string_view label); + void polish (); + bool converged() const { return converged_; }; + + // physical properties + virtual std::string output() const; + void print_dim_info() const; + Real E() const { return E_; }; + int iK() const { return iK_; }; + Real K() const { return K_; }; + + // utilities: + virtual std::string get_filename_prefix () const = 0; + bool operator== (const BetheState& rhs) const { + // return (model_ == rhs.model_ && label_ == rhs.label_); + return (model_.get() == rhs.model_.get() && has_identical_Ix2(rhs)); + } + bool has_identical_Ix2 (const BetheState& rhs) const; + + // identification: + std::string_view label() const { return label_; }; + std::string_view baselabel() const { return baselabel_; }; + std::string_view patternlabel() const { return patternlabel_; }; + + // construction, convergence, information: + Real δB() const { return δB_; }; + int charge() const { return charge_; }; + int iter_count (IterationMethod method) const { return iter_count_[method]; }; + std::string iter_details () const; + + +public: + + std::string label_; + std::string baselabel_; + std::string patternlabel_; + + std::reference_wrapper model_; + int charge_; // number of Bethe rapidities (all-inclusive, an n-string counts as n) + int string_charge_; // number of higher string centers (including ground level and all up) + + // one-strings are handled via the inherited Plex + + // precomputed phase shifts to accelerate iterating Bethe equations, + // by exploiting symmetry + Matrix φ_; + Matrix dφdƛ_; + + Real δB_ { Real(0) }; // sum of abs(B_) + + Matrix Gaudin_g_; // ground level only + Matrix Gaudin_; // all-inclusive + Real ln_Gaudin_det_; + Real lnnorm_; + + std::array iter_count_; + std::array iter_time_; + bool converged_ { false }; + + Real E_; + int iK_; + Real K_; + +protected: // member functions + // setting quantum numbers + void set_label_from_plexlabels_(); + std::string label_for_modified_g_plexlabel_ + (std::string modified_g_plexlabel) const; + std::string label_for_modified_str_plexlabel_ + (IndexU jmod, std::string modified_str_plexlabel) const; + // std::vector descendents_g_ (DescendentType type); + bool is_compatible_ (const ParsedLabel& parsed); + + // Scattering sums for ground level rapidities: + void compute_φ_ (); + Real Σφ_ (IndexU α) const; + void compute_dφdƛ_ (); + Real Σdφdƛ_ (IndexU α) const; + + // Bethe equations and Gaudin matrix + void compute_B_(); + void build_Gaudin_g_(); + void build_Gaudin_(); + void compute_Gaudin_det_(); + void build_Gaudin_g_diagonal(); + void build_Gaudin_g_tridiagonal (); + void build_Gaudin_g_pentadiagonal (); + void apply_g_Thomas_algorithm (std::vector& scratch); + void apply_g_pentadiagonal_algorithm (std::vector>& scratch); + + void reset_iter_info() { + std::fill(iter_count_.begin(), iter_count_.end(), 0); + std::fill(iter_time_.begin(), iter_time_.end(), 0.0); + } + +public: // member functions + virtual void set_ground_state_Ix2() = 0; + virtual void initialize() = 0; + void iterate_Bethe_equations_ouroboros (); + void iterate_Bethe_equations_g_diagonal (); + void iterate_Bethe_equations_diagonal (); + void iterate_Bethe_equations_tridiagonal (std::vector& scratch); + void iterate_Bethe_equations_pentadiagonal (std::vector>& scratch); + void iterate_Bethe_equations_g_Newton (); + void iterate_Bethe_equations_Newton (); + void iterate_Bethe_equations (IterationMethod method, std::vector>& scratch); + virtual bool Gaudin_g_left_edge_is_monotonic () const; + virtual bool Gaudin_g_right_edge_is_monotonic () const; + virtual void control_δƛ (); + void shift_all_ƛ_with_δƛ (); + void approach_solution_to_Bethe_equations + ( + IterationMethod method=IterationMethod::diagonal + ); + virtual void compute_lnnorm () = 0; + virtual void compute_Momentum () = 0; + virtual void compute_Energy () = 0; + virtual void populate_λ () = 0; +}; + +// constructors +template TModel> +BetheState::BetheState (TModel& model, int g_f) + : Plex(g_f) + , λ { std::vector(g_f) } + , model_ { model } + , charge_ { g_f } + , string_charge_ { g_f } + , φ_ { Matrix(g_f) } + , dφdƛ_ { Matrix(g_f) } + , Gaudin_g_ { Matrix(g_f) } + , Gaudin_ { Matrix(string_charge_) } + // , iter_count_ { std::vector(IterationMethod::maxIterationMethod) } + // , iter_time_ { std::vector(IterationMethod::maxIterationMethod) } +{ + set_label_from_plexlabels_(); +} + +template TModel> +BetheState::BetheState +(const BetheState& refstate, int Δf, bool relative_label) + : Plex(refstate.l_, refstate.p_, relative_label ? refstate.Ix2_ : refstate.Ox2_, Δf) + , λ { std::vector(int(refstate.Ox2_.size()) + Δf) } + , model_ { refstate.model_ } + , charge_ { refstate.charge_ + Δf } + , string_charge_ { refstate.string_charge_ + Δf } + , φ_ { Matrix(λ.size()) } + , dφdƛ_ { Matrix(λ.size()) } + , Gaudin_g_ { Matrix(λ.size()) } + , Gaudin_ { Matrix(string_charge_) } + // , iter_count_ { std::vector(IterationMethod::maxIterationMethod) } + // , iter_time_ { std::vector(IterationMethod::maxIterationMethod) } +{ + set_label_from_plexlabels_(); +} + +template TModel> +Real BetheState::ln_subspace_dim_this_base () const { + return ln_dim_; +} + +template TModel> +Real BetheState::ln_subspace_dim_this_pattern () const { + Real ln_dim { 0 }; + + // ground level dimension: + // putting nex particles in (Ix2_max-Ix2_min)/2 + 1 - f_ slots, + // putting nex holes in f_ slots + IndexU nex { count_nex_in_plexlabel(plexlabel_) }; + if (nex > 0) { + ln_dim += std::lgamma(Real((Ix2_max_ - Ix2_min_)/2 + 2 - f_)) + - std::lgamma(Real((Ix2_max_ - Ix2_min_)/2 + 2 - f_ - nex)) + + std::lgamma(Real(f_ + 1)) - std::lgamma(Real(f_ - nex + 1)) - + 2*std::lgamma(Real(nex + 1)); + } + return ln_dim; +} + +template TModel> +Real BetheState::ln_subspace_dim_for_pattern (std::string patternlabel) const { + // read the nex from the patternlabel + int nex_g; + + std::stringstream label_stream; + label_stream << patternlabel; + std::string tmp; + + std::getline(label_stream, tmp, LABEL_EXC_SEPARATOR); + tmp.clear(); + std::getline(label_stream, tmp); + nex_g = tmp.empty() ? 0 : stoi(tmp, nullptr, 16); + + Real ln_dim { 0 }; + + // ground level dimension: + // putting nex particles in (Ix2_max-Ix2_min)/2 + 1 - f_ slots, + // putting nex holes in f_ slots + if (nex_g > 0) { + ln_dim += std::lgamma(Real((Ix2_max_ - Ix2_min_)/2 + 2 - f_)) + - std::lgamma(Real((Ix2_max_ - Ix2_min_)/2 + 2 - f_ - nex_g)) + + std::lgamma(Real(f_ + 1)) - std::lgamma(Real(f_ - nex_g + 1)) + - 2*std::lgamma(Real(nex_g + 1)); + } + if (std::isnan(ln_dim)) ln_dim = std::numeric_limits::infinity(); + return ln_dim; +} + + +// public manipulation +template TModel> +bool BetheState::draft (std::string_view label) { + // checks compatibility of label, and if compatible, + // sets + // - label (and plexlabels) + // - momentum + + // std::string label_given { label }; + // ParsedLabel parsed(label_given); + ParsedLabel parsed(label); + + if (is_compatible_(parsed)) { + parsed.parsed_plexlabel_g_.set_Ix2(Ix2_, Lx2_); + set_plexlabel_from_Ix2_(); + set_label_from_plexlabels_(); + compute_Momentum(); + return true; + } + return false; +} + +template TModel> +bool BetheState::set_g_Ix2 (const std::vector& Ix2_req) { + if (Plex::set_Ix2(Ix2_req)) { + set_label_from_plexlabels_(); + return true; + } + return false; +} + +template TModel> +bool BetheState::excite_g (IndexU α, int δI) { + // shifts the ground Ix2_[α] by δI * 2 if this is allowed + if (Plex::excite(α, δI)) { + set_label_from_plexlabels_(); + return true; + } + return false; +} + +template TModel> +void BetheState::boost (int δI) { + Plex::boost(δI); + set_label_from_plexlabels_(); + compute_Momentum(); +} + +// template TModel> +// std::vector BetheState::descendents (DescendentType type) +// { +// Chirality chirality { std::get(type) }; + +// IndexU g_exc_α { lowest_excitation_(chirality) }; + +// // all descendents imply modifications at the ground level only +// return descendents_g_(type); +// } + +template TModel> +void BetheState::polish () { + approach_solution_to_Bethe_equations (); +} + +template TModel> +bool BetheState::has_identical_Ix2 (const BetheState& rhs) const { + if (Ix2_ != rhs.Ix2_) return false; + return true; +} + +template TModel> +std::string BetheState::iter_details () const { + std::stringstream output { }; + for (int method { 0 }; method < IterationMethod::maxIterationMethod; ++method) { + if (iter_count_[method] > 0) + output << getIterationMethodName(method) << ": " << iter_count_[method] + << " iterations in " + << iter_time_[method] << " seconds\t"; + } + return output.str(); +} + +// protected member functions + +template TModel> +void BetheState::set_label_from_plexlabels_ () { + std::stringstream label_stream; + label_stream << std::hex << f_; + std::stringstream baselabel_stream; + baselabel_stream << std::hex << f_; + std::stringstream patternlabel_stream; + patternlabel_stream << std::hex << f_; + IndexU nex { 0 }; + + if (!Plex::plexlabel_.empty()) { + label_stream << LABEL_EXC_SEPARATOR << Plex::plexlabel_; + nex = count_nex_in_plexlabel(Plex::plexlabel_); + if (nex > 0) patternlabel_stream << LABEL_EXC_SEPARATOR << nex; + } + + label_ = compress(label_stream.str()); + baselabel_ = baselabel_stream.str(); + patternlabel_ = patternlabel_stream.str(); +} + +template TModel> +std::string BetheState::label_for_modified_g_plexlabel_ +(std::string modified_g_plexlabel) const { + // Give the state's label if the ground plexlabel was modified + // Used when computing descendent labels + + std::stringstream label_stream; + label_stream << std::hex << f_ << LABEL_EXC_SEPARATOR << modified_g_plexlabel; + + return compress(label_stream.str()); +} + +// template TModel> +// std::vector BetheState::descendents_g_ (DescendentType type) +// { +// std::vector desc_plexlabels; + +// append_descendent_plexlabels(type, desc_plexlabels); + +// // translate the changed plexlabels into full labels +// std::vector desc_g_labels; +// for (std::string plexlabel: desc_plexlabels) +// if (!plexlabel.empty()) +// desc_g_labels.push_back(label_for_modified_g_plexlabel_(plexlabel)); + +// return (desc_g_labels); +// } + +template TModel> +bool BetheState::is_compatible_ (const ParsedLabel& parsed) +{ + // check if parsed label is compatible with state + if (f_ != parsed.filling_g_ || // wrong ground filling + !parsed.parsed_plexlabel_g_.is_compatible(Lx2_)) + return false; + return true; +} + +template TModel> +void BetheState::compute_φ_ () { + if (δB_ > ten_double_eps) { // use faster phase + for (IndexU α { 0 }; α + 1 < f_; ++α) { + for (IndexU β { α+1 }; β < f_; ++β) { + φ_(α, β) = model_.get().φ_(static_cast(ƛ_[α] - ƛ_[β])); + } + } + } + else { + for (IndexU α { 0 }; α + 1 < f_; ++α) { + for (IndexU β { α+1 }; β < f_; ++β) { + φ_(α, β) = model_.get().φ_(ƛ_[α] - ƛ_[β]); + } + } + } +} + +template TModel> +Real BetheState::Σφ_ (IndexU α) const { + // Scattering sum for a ground level rapidity + Real sum { Real(0) }; + for (IndexU β { 0 }; β < α; ++β) sum -= φ_ (β, α); + for (IndexU β { α+1 }; β < f_; ++β) sum += φ_ (α, β); + return sum; +} + +template TModel> +void BetheState::compute_dφdƛ_ () { + for (IndexU α { 0 }; α + 1 < f_; ++α) { + for (IndexU β { α+1 }; β < f_; ++β) { + dφdƛ_(α, β) = model_.get().dφdƛ_(ƛ_[α] - ƛ_[β]); + } + } +} + +template TModel> +Real BetheState::Σdφdƛ_ (IndexU α) const { + // Derivative of scattering sum for a ground level rapidity + Real sum = Real(0.0); + for (IndexU β { 0 }; β < α; ++β) sum += dφdƛ_ (β, α); + for (IndexU β { α+1 }; β < f_; ++β) sum += dφdƛ_ (α, β); + + return sum; +} + +template TModel> +void BetheState::compute_B_ () { + // Expectation: δB ~ N^2 epsilon + // (each B has factor N and precision epsilon; there are N of them) + δB_ = Real(0.0); + for (IndexU α { 0 }; α < f_; ++α) { + B_[α] = model_.get().Ł_ * model_.get().θ_(ƛ_[α]) - Σφ_(α) - pi_r * Ix2_[α]; + δB_ += std::fabs(B_[α]); + } + δB_ /= model_.get().Ł_ * charge(); // scale so that δB_ ~ \epsilon +} + +template TModel> +void BetheState::build_Gaudin_g_ () { + compute_dφdƛ_(); + for (IndexU α { 0 }; α < f_; ++α) { + Gaudin_g_(α, α) = model_.get().Ł_ * model_.get().dθdƛ_(ƛ_[α]) - Σdφdƛ_(α); + for (IndexU β { α + 1 }; β < f_; ++β) { + Gaudin_g_(α, β) = dφdƛ_(α, β); + Gaudin_g_(β, α) = Gaudin_g_(α, β); + } + } +} + +template TModel> +void BetheState::build_Gaudin_ () +{ + compute_dφdƛ_(); + for (IndexU α { 0 }; α < f_; ++α) { + Gaudin_(α, α) = model_.get().Ł_ * model_.get().dθdƛ_(ƛ_[α]) - Σdφdƛ_(α); + for (IndexU β { α + 1 }; β < f_; ++β) { + Gaudin_(α, β) = dφdƛ_(α, β); + Gaudin_(β, α) = Gaudin_(α, β); + } + } +} + +template TModel> +void BetheState::compute_Gaudin_det_ () { + build_Gaudin_(); + Gaudin_ /= model_.get().Ł_; + + ln_Gaudin_det_ = real(Gaudin_.lndet_LU_destroy()); +} + +template TModel> +void BetheState::build_Gaudin_g_diagonal () { + Gaudin_g_.setZero(); + IndexU α { 0 }; + + for (α = 0; α < f_; ++α) { + Gaudin_g_(α, α) = model_.get().Ł_ * model_.get().dθdƛ_(ƛ_[α]) - Σdφdƛ_(α); + } +} + +template TModel> +void BetheState::build_Gaudin_g_tridiagonal () { + build_Gaudin_g_diagonal(); + for (IndexU α { 1 }; α < f_; ++α) + Gaudin_g_(α, Ξ(α)-1) = model_.get().dφdƛ_(ƛ_[α] - ƛ_[Ξ(Ξ(α)-1)]); + for (IndexU α { 0 }; α < f_ - 1; ++α) + Gaudin_g_(α, α+1) = model_.get().dφdƛ_(ƛ_[α] - ƛ_[α+1]); +} + +template TModel> +void BetheState::build_Gaudin_g_pentadiagonal () { + build_Gaudin_g_tridiagonal(); + for (IndexU α { 2 }; α < f_; ++α) + Gaudin_g_(α, Ξ(α)-2) = model_.get().dφdƛ_(ƛ_[α] - ƛ_[Ξ(Ξ(α)-2)]); + for (IndexU α { 0 }; α < f_ - 2; ++α) + Gaudin_g_(α, α+2) = model_.get().dφdƛ_(ƛ_[α] - ƛ_[α+2]); +} + +template TModel> +void BetheState::apply_g_Thomas_algorithm (std::vector& scratch) { + // See https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm + // a[i] is Gaudin[i][i-1], b[i] is Gaudin[i][i], c[i] is Gaudin[i][i+1] + for (IndexU α { 0 }; α < f_; ++α) δƛ_[α] = -B_[α]; + scratch[0] = Gaudin_g_(0,1)/Gaudin_g_(0,0); + δƛ_[0] = δƛ_[0]/Gaudin_g_(0,0); + for (IndexU ix { 1 }; ix < f_; ix++) { + if (ix+1 < f_) + scratch[ix] = Gaudin_g_(ix,ix+1) / + (Gaudin_g_(ix,ix) - Gaudin_g_(ix,Ξ(Ξ(ix)-1)) * scratch[Ξ(Ξ(ix)-1)]); + δƛ_[ix] = (δƛ_[ix] - Gaudin_g_(ix,Ξ(Ξ(ix)-1)) * δƛ_[Ξ(Ξ(ix)-1)]) / + (Gaudin_g_(ix,ix) - Gaudin_g_(ix,Ξ(Ξ(ix)-1)) * scratch[Ξ(Ξ(ix)-1)]); + } + for (IndexS ix { f_ - 2 }; ix >= 0; ix--) + δƛ_[Ξ(ix)] -= scratch[Ξ(ix)] * δƛ_[Ξ(ix) + 1]; +} + +template TModel> +void BetheState::apply_g_pentadiagonal_algorithm +(std::vector>& scratch) { + // see https://www.hindawi.com/journals/mpe/2015/232456/ + // but be careful for the error it contains (see below) + // For i = 0, ..., n-1 (shift by 1) in the paper, identification with Gaudin matrix elements: + // a[i] = G[i][i+1] for i = 0, ..., n-2 + // b[i] = G[i][i+2] for i = 0, ..., n-3 + // d[i] = G[i][i] + // c[i] = G[i][i-1] for i = 1, ..., n-1 + // e[i] = G[i][i-2] for i = 2, ..., n-1 + // The input y[i] = -B_[i] + + int n { f_ }; + Real ga, mu; // forward substitutions for gamma, mu + + // alpha vector: scratch[0] + // beta vector: scratch[1] + // z vector: scratch[2] + + // From the algorithm (with shifted indices) + // mu[0] = d[0]; + // al[0] = a[0]/mu[0]; + // be[0] = b[0]/mu[0]; + // z[0] = y[0]/mu[0]; + // ga[1] = c[1]; + // mu[1] = d[1] - al[0] * ga[1]; + // al[1] = (a[1] - be[0] * ga[1])/mu[1]; + // be[1] = b[1]/mu[1]; + // z[1] = (y[1] - z[0] * ga[1])/mu[1]; + + // With direct substitutions + mu = Gaudin_g_(0,0); + scratch[0][0] = Gaudin_g_(0,1)/mu; + scratch[1][0] = Gaudin_g_(0,2)/mu; + scratch[2][0] = -B_[0]/mu; + ga = Gaudin_g_(1,0); + mu = Gaudin_g_(1,1) - scratch[0][0] * ga; + scratch[0][1] = (Gaudin_g_(1,2) - scratch[1][0] * ga)/mu; + scratch[1][1] = Gaudin_g_(1,3)/mu; + scratch[2][1] = (-B_[1] - scratch[2][0] * ga)/mu; + + // for (i = 2; i < n-4; ++i) { + // ga[i] = c[i] - al[i-2] * e[i]; + // mu[i] = d[i] - be[i-2] * e[i] - al[i-1] * ga[i]; + // al[i] = (a[i] - be[i-1] * ga[i])/mu[i]; + // be[i] = b[i]/mu[i]; + // z[i] = (y[i] - z[i-2] * e[i] - z[i-1] * ga[i])/mu[i]; + // } + for (IndexU i {0}; i < IndexU(n-4); ++i) { + ga = Gaudin_g_(i+2,i+1) - scratch[0][i] * Gaudin_g_(i+2,i); + mu = Gaudin_g_(i+2,i+2) - scratch[1][i] * Gaudin_g_(i+2,i) - scratch[0][i+1] * ga; + scratch[0][i+2] = (Gaudin_g_(i+2,i+3) - scratch[1][i+1] * ga)/mu; + scratch[1][i+2] = Gaudin_g_(i+2,i+4)/mu; + scratch[2][i+2] = (-B_[i+2] - scratch[2][i] * Gaudin_g_(i+2,i) - scratch[2][i+1] * ga)/mu; + } + + // ga[n-2] = c[n-2] - al[n-4] * e[n-2]; + // mu[n-2] = d[n-2] - be[n-4] * e[n-2] - al[n-3] * ga[n-2]; + // al[n-2] = (a[n-2] - be[n-3] * ga[n-2])/mu[n-2]; + // ga[n-1] = c[n-1] - al[n-3] * e[n-1]; + // mu[n-1] = d[n-1] - be[n-3] * e[n-1] - al[n-2] * ga[n-1]; + // z[n-2] = (y[n-2] - !!z[n-4] * e[n-2] - z[n-3] * ga[n-2])/mu[n-2]; // error in paper: at !!, z[n-4] instead of z[n-3] + // z[n-1] = (y[n-1] - !!z[n-3] * e[n-1] - z[n-2] * ga[n-1])/mu[n-1]; // same error as above: n-3 instead of n-2 + ga = Gaudin_g_(n-2,n-3) - scratch[0][n-4] * Gaudin_g_(n-2,n-4); + mu = Gaudin_g_(n-2,n-2) - scratch[1][n-4] * Gaudin_g_(n-2,n-4) - scratch[0][n-3] * ga; + scratch[0][n-2] = (Gaudin_g_(n-2,n-1) - scratch[1][n-3] * ga)/mu; + scratch[2][n-2] = (-B_[n-2] - scratch[2][n-4] * Gaudin_g_(n-2,n-4) - scratch[2][n-3] * ga)/mu; + ga = Gaudin_g_(n-1,n-2) - scratch[0][n-3] * Gaudin_g_(n-1,n-3); + mu = Gaudin_g_(n-1,n-1) - scratch[1][n-3] * Gaudin_g_(n-1,n-3) - scratch[0][n-2] * ga; + scratch[2][n-1] = (-B_[n-1] - scratch[2][n-3] * Gaudin_g_(n-1,n-3) - scratch[2][n-2] * ga)/mu; + + δƛ_[n-1] = scratch[2][n-1]; + δƛ_[n-2] = scratch[2][n-2] - scratch[0][n-2] * δƛ_[n-1]; + + for (IndexS i {n-3}; i >= 0; --i) { + δƛ_[Ξ(i)] = scratch[2][Ξ(i)] - scratch[0][Ξ(i)] * δƛ_[Ξ(i)+1] - scratch[1][Ξ(i)] * δƛ_[Ξ(i)+2]; + } +} + +template TModel> +void BetheState::iterate_Bethe_equations_ouroboros () +{ + Real f = 1.0L/model_.get().Ł_; + + for (IndexU α { 0 }; α < f_; ++α) + δƛ_[α] = model_.get().θinv_(model_.get().θ_(ƛ_[α]) -B_[α] * f) - ƛ_[α]; +} + +template TModel> +void BetheState::iterate_Bethe_equations_g_diagonal () { + build_Gaudin_g_diagonal(); + for (IndexU α { 0 }; α < f_; ++α) δƛ_[α] = -B_[α]/Gaudin_g_(α, α); +} + +template TModel> +void BetheState::iterate_Bethe_equations_diagonal () +{ + iterate_Bethe_equations_g_diagonal(); +} + +template TModel> +void BetheState::iterate_Bethe_equations_tridiagonal (std::vector& scratch) { + build_Gaudin_g_tridiagonal(); + apply_g_Thomas_algorithm(scratch); +} + +template TModel> +void BetheState::iterate_Bethe_equations_pentadiagonal +( + std::vector>& scratch + ) { + build_Gaudin_g_pentadiagonal(); + apply_g_pentadiagonal_algorithm(scratch); +} + +template TModel> +void BetheState::iterate_Bethe_equations_g_Newton () +{ + Timer t; + + build_Gaudin_g_(); + Gaudin_g_ /= model_.get().Ł_; + + std::vector RHS(f_); + for (IndexU α { 0 }; α < f_; ++α) RHS[α] = -B_[α]/model_.get().Ł_; + std::vector indx(f_); + Real d; + Gaudin_g_.ludcmp(indx, d); + Gaudin_g_.lubksb(indx, RHS); + δƛ_ = RHS; +} + +template TModel> +void BetheState::iterate_Bethe_equations_Newton () +{ + iterate_Bethe_equations_g_Newton(); +} + +template TModel> +void BetheState::iterate_Bethe_equations +( + IterationMethod method, + std::vector>& scratch + ) +{ + compute_dφdƛ_(); + + switch (method) + { + case ouroboros: // Method: straight iter + iterate_Bethe_equations_ouroboros(); + break; + + case diagonal: // Method: diagonal Newton + iterate_Bethe_equations_diagonal(); + break; + + case tridiagonal: // Method: Newton but using only tridiagonal form, + iterate_Bethe_equations_tridiagonal(scratch[0]); + break; + + case pentadiagonal: // Method: Newton but using only pentadiagonal form, + iterate_Bethe_equations_pentadiagonal(scratch); + break; + + case Newton: // Method: Newton + iterate_Bethe_equations_g_Newton(); + break; + + default: + throw std::invalid_argument("Not a known iteration method"); + + } // switch (method) + + control_δƛ(); + shift_all_ƛ_with_δƛ(); + compute_φ_(); + compute_B_(); + ++iter_count_[method]; +} + +template TModel> +bool BetheState::Gaudin_g_left_edge_is_monotonic () const { + // Check the sign of the leftmost ground rapidity's Gaudin diagonal after the update + Real sum = Real(0.0); + for (IndexU α { 1 }; α < f_; ++α) sum += model_.get().dφdƛ_(ƛ_[0] + δƛ_[0] - (ƛ_[α] + δƛ_[α])); + return model_.get().Ł_ * model_.get().dθdƛ_(ƛ_[0] + δƛ_[0]) - sum > 0; +} + +template TModel> +bool BetheState::Gaudin_g_right_edge_is_monotonic () const { + // Check the sign of the rightmost ground rapidity's Gaudin diagonal after the update + Real sum = Real(0.0); + IndexU α_r = Ξ(std::ssize(ƛ_)-1); + for (IndexU α { 0 }; α < α_r; ++α) sum += model_.get().dφdƛ_(ƛ_[α_r] + δƛ_[α_r] - (ƛ_[α] + δƛ_[α])); + return model_.get().Ł_ * model_.get().dθdƛ_(ƛ_[α_r] + δƛ_[α_r]) - sum > 0; +} + +template TModel> +void BetheState::control_δƛ() { + // Check if the rightmost rapidity's Gaudin diagonal remains positive after the update. + // If not, pull back. + int ctr { 0 }; + int max_ctr { 6 }; + if (f_ > 0) { + while (!Gaudin_g_left_edge_is_monotonic() && ctr++ < max_ctr) { // pull back the left edge rapidity + δƛ_.front() *= 0.5; + } + if (ctr == max_ctr) δƛ_.front() = 0; // give up and reset + + ctr = 0; + while (!Gaudin_g_right_edge_is_monotonic() && ctr++ < max_ctr) { // pull back the right edge rapidity + δƛ_.back() *= 0.5; + } + if (ctr == max_ctr) δƛ_.back() = 0; // give up and reset + } +} + +template TModel> +void BetheState::shift_all_ƛ_with_δƛ() { + shift_ƛ_with_δƛ_(); +} + +template TModel> +void BetheState::approach_solution_to_Bethe_equations +( + IterationMethod method + ) { + iter_count_[method] = 0; + iter_time_[method] = 0.0; + Timer timer; + std::vector> scratch { + std::vector(f_), + std::vector(f_), + std::vector(f_), + }; + Real previous_δB_; + iterate_Bethe_equations(method, scratch); // do at least one iteration + do { + previous_δB_ = δB_; + iterate_Bethe_equations(method, scratch); + } while (((δB_ > sqrt_real_eps && iter_count_[method] < 1000) + || δB_ < previous_δB_)); + + converged_ = δB_ < sqrt_real_eps; + + if (converged_) { + compute_Gaudin_det_(); + compute_lnnorm(); + compute_Energy(); + populate_λ(); + } + iter_time_[method] += timer.elapsed(); +} + +template TModel> +void BetheState::print_dim_info() const { + std::cout << "Ix2 max: " << Ix2_max_ << "\t" + << "Dimensionality: " << ln_dim_ << "\t"; + if (ln_dim_ < std::logl(std::numeric_limits::max())) + std::cout << static_cast(std::exp(ln_dim_)+0.5l); + else std::cout << std::numeric_limits::infinity(); + std::cout << "\n"; +} + +template TModel> +std::string BetheState::output() const { + std::stringstream s; + s << "State with label " << label_ << " (baselabel " << baselabel_ << ")\n"; + s << "Converged: " << converged() << "\tIterations: " + << iter_details() << "\n\t\tTotal iterations: " + << std::accumulate(iter_count_.begin(), iter_count_.end(), 0) << " in " + << std::accumulate(iter_time_.begin(), iter_time_.end(), float(0)) << " seconds\n" + << "\t\tResulting δB: " << δB_ << "\n"; + s << "\nQuantum numbers:\nGround level:\n"; + for (int n : Ix2_) s << n << "\t"; + s << "\n"; + return s.str(); +} + + +// Template specialization concept +export template +concept BetheStateOf = std::is_base_of, TState>::value; + + +//////////////////////// +// ↑ Class BetheState // +//////////////////////// diff --git a/src/math/calculus.cc b/src/math/calculus.cc new file mode 100644 index 0000000..043a620 --- /dev/null +++ b/src/math/calculus.cc @@ -0,0 +1,126 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +export module calculus; + + +import std; + + +import conveniences; + + +export class Interval { +public: + std::vector λ; + std::vector dλ; + +public: + Interval() {}; + Interval(Real λmax, int npts); + Interval(int npts); /// for improper integrals from -∞ to ∞ +}; + +Interval::Interval(Real λmax, int npts) + : λ(std::vector(npts)) + , dλ(std::vector(npts)) +{ + for (int i { 0 }; i < npts; ++i) { + λ[Ξ(i)] = -λmax + (2*i + 1)* λmax/npts; + dλ[Ξ(i)] = 2*λmax/npts; + } +} + +Interval::Interval(int npts) + : λ(std::vector(npts)) + , dλ(std::vector(npts)) +{ + /// Use coordinate mapping λ = tan θ, with θ uniformly spaced in ]-π/2, π/2 [. + // for (int i { 0 }; i < npts; ++i) + // λ[Ξ(i)] = std::tan((2*i + 2 - npts)* pi_r/(2*npts)); + + /// Use λ = ξ e^(ξ^2) with ξ in ]-10, 10[ + // Real ξ; + // for (int i { 0 }; i < npts; ++i) { + // ξ = (2*i + 2 - npts)* Real(10)/npts; + // λ[Ξ(i)] = ξ * std::exp(ξ*ξ); + // } + + /// Use λ = ξ/(1 - ξ^2) with ξ in ]-1, 1[ + Real ξ; + for (int i { 0 }; i < npts; ++i) { + ξ = (2*i + 1 - npts)* Real(1)/npts; + λ[Ξ(i)] = ξ/(Real(1) - ξ*ξ); + } + + // Fill in the differential elements + for (IndexU i { 1 }; i < npts - 1; ++i) + dλ[i] = (λ[i+1] - λ[i-1])/2; + dλ[0] = dλ[1]; + dλ[Ξ(npts-1)] = dλ[Ξ(npts-2)]; +} + + +export class Function { + +public: + Interval interval_; + std::vector val_; + std::vector val_prev_; + Real δval_; // Σ abs(val - val_prev) + +public: + Function() {}; + Function(const Interval& interval); + + Real evaluate_at (Real coordinate); /// returns linearly interpolated value +}; + +Function::Function(const Interval& interval) + : interval_(interval) + , val_(std::vector(interval.λ.size(), Real(0))) + , val_prev_(std::vector(interval.λ.size(), Real(0))) + , δval_(std::numeric_limits::max()) +{} + +Real Function::evaluate_at (Real coordinate) +{ + if (coordinate <= interval_.λ.front() || + coordinate > interval_.λ.back()) return Real(0); + + // find index i such that λ[i-1] < coordinate <= λ[i] + IndexU i { interval_index_in_ordered(coordinate, interval_.λ) }; + + // if (interval_.λ[i-1] >= coordinate || interval_.λ[i] < coordinate) { + // std::cout << i << "\t" << interval_.λ[i-1] << " +class Matrix { + +private: + std::size_t dim_; + std::vector element_; + +public: + Matrix (std::size_t dim); + + inline std::size_t size() const { return dim_; }; + + inline T operator() (IndexU i, IndexU j) const { return element_[i*dim_ + j]; }; + inline T& operator() (IndexU i, IndexU j) { return element_[i*dim_ + j]; }; + + void setZero(); + + Matrix& operator*= (const T& a); + Matrix& operator/= (const T& a); + + void ludcmp (std::vector& indx, T& d); + void lubksb (std::vector& indx, std::vector& b); + std::complex lndet_LU_destroy (); +}; + +template +Matrix::Matrix (std::size_t dim) + : dim_(dim) + , element_(std::vector(dim_*dim_)) +{} + +template +void Matrix::setZero () +{ + std::fill(element_.begin(), element_.end(), T(0)); +} + + +template +Matrix& Matrix::operator*= (const T& a) +{ + std::transform(element_.begin(), element_.end(), element_.begin(), [a](T el) { return el * a; }); + return *this; +} + +template +Matrix& Matrix::operator/= (const T& a) +{ + T oneovera = T(1)/a; + std::transform(element_.begin(), element_.end(), element_.begin(), [oneovera](T el) { return el * oneovera; }); + return *this; +} + + +template +void Matrix::ludcmp (std::vector& indx, T& d) +{ + IndexU i, j, k; + IndexU imax { 0 }; + IndexU idim_, jdim_, imaxdim_; + T big, dum, sum, temp; + + IndexU n { size() }; + + std::vector vv(n); + d = T(1); + for (i = 0; i < n; i++) { + big = T(0); + idim_ = i*dim_; + for (j = 0; j < n; j++) { + if ((std::fabs(temp = element_[idim_ + j])) > std::fabs(big)) big = temp; + } + if (big == T(0)) throw; + vv[i] = T(1)/big; + } + + for (j = 0; j < n; j++) { + for (i = 0; i < j; i++) { + idim_ = i*dim_; + sum = element_[idim_ + j]; + for (k = 0; k < i; k++) sum -= element_[idim_ + k] * element_[k*dim_ + j]; + element_[idim_ + j] = sum; + } + big = T(0); + for (i = j; i < n; i++) { + idim_ = i*dim_; + sum = element_[idim_ + j]; + for (k = 0; k < j; k++) sum -= element_[idim_ + k] * element_[k*dim_ + j]; + element_[idim_ + j] = sum; + if ((std::fabs(dum = vv[i]*sum)) >= std::fabs(big)) { + big = dum; + imax = i; + } + } + jdim_ = j*dim_; + if (j != imax) { + imaxdim_ = imax*dim_; + for (k = 0; k < n; k++) { + dum = element_[imaxdim_ + k]; + element_[imaxdim_ + k] = element_[jdim_ + k]; + element_[jdim_ + k] = dum; + } + d = -d; + vv[imax] = vv[j]; + } + indx[j] = imax; + if (j !=n-1) { + dum = T(1)/(element_[jdim_ + j]); + for (i = j + 1; i < n; i++) element_[i*dim_ + j] *= dum; + } + } +} + +template +void Matrix::lubksb (std::vector& indx, std::vector& b) +{ + // int i, ip, j; + // int ii { 0 }; + // int idim_; + // T sum; + + // int n { int(size()) }; + // for (i = 0; i < n; i++) { + // ip = indx[i]; + // sum = b[ip]; + // b[ip] = b[i]; + // idim_ = i*dim_; + // if (ii != 0) + // for (j = ii-1; j < i; j++) sum -= element_[idim_ + j] * b[j]; + // else if (sum != T(0)) + // ii = i + 1; + // b[i] = sum; + // } + // for (i = n - 1; i >= 0; i--) { + // sum = b[i]; + // idim_ = i*dim_; + // for (j = i + 1; j < n; j++) sum -= element_[idim_ + j] * b[j]; + // b[i] = sum/element_[idim_ + i]; + // } + + std::size_t i, ip, j; + std::size_t ii { 0 }; + std::size_t idim_; + T sum; + std::size_t n { size() }; + for (i = 0; i < n; i++) { + ip = indx[i]; + sum = b[ip]; + b[ip] = b[i]; + idim_ = i*dim_; + if (ii != 0) + for (j = ii-1; j < i; j++) sum -= element_[idim_ + j] * b[j]; + else if (sum != T(0)) + ii = i + 1; + b[i] = sum; + } + for (i = n; i--;) { + sum = b[i]; + idim_ = i*dim_; + for (j = i + 1; j < n; j++) sum -= element_[idim_ + j] * b[j]; + b[i] = sum/element_[idim_ + i]; + } +} + +template +std::complex Matrix::lndet_LU_destroy () +{ + std::vector indx(size()); + T d; + std::complex lndet { 0 }; + + (*this).ludcmp(indx, d); + + lndet = log(std::complex(d)); + + for (IndexU j { 0 }; j < size(); j++) { + lndet += log(std::complex(element_[j*dim_ + j])); + } + + return lndet; +} + + +export template +std::ostream& operator<< (std::ostream& s, const Matrix& M) +{ + for (IndexU i { 0 }; i < M.size(); ++i) { + for (IndexU j { 0 }; j < M.size(); ++j) s << M(i,j) << "\t"; + s << std::endl; + } + return s; +} + +//////////////////// +// ↑ Class Matrix // +//////////////////// + + + +// Additional utilities + + +inline Real SIGN (const Real& a, const Real& b) +{ + return b >= Real(0) ? (a >= Real(0) ? a : -a) : (a >= Real(0) ? -a : a); +} + +Real pythag(const Real a, const Real b) +{ + Real absa, absb; + absa = std::abs(a); + absb = std::abs(b); + if (absa > absb) return absa * std::sqrt(Real(1) + (absb/absa)*(absb/absa)); + else return absb * std::sqrt(Real(1) + (absa/absb)*(absa/absb)); +} + + +// Singular value decomposition +// Numerical Recipes section 2.6 +// Restricted to square matrix due to use of Matrix class +export void svdcmp(Matrix& a, std::vector& w, Matrix& v) +{ + bool flag; + int i, its, j, jj, k, l, nm; + Real anorm, c, f, g, h, s, scale, x, y, z; + + int m { int(a.size()) }; + int n { int(a.size()) }; + + std::vector rv1(n); + g = Real(0); scale = Real(0); anorm = Real(0); + + for (i = 0; i < n; i++) { + l = i+2; + rv1[i] = scale*g; + g = Real(0); s = Real(0); scale = Real(0); + if (i < m) { + for (k = i; k < m; k++) scale += std::abs(a(k,i)); + if (scale != Real(0)) { + for (k = i; k < m; k++) { + a(k,i) /= scale; + s += a(k,i) * a(k,i); + } + f = a(i,i); + g = -SIGN(std::sqrt(s), f); + h = f*g-s; + a(i,i) = f-g; + for (j = l-1; j < n; j++) { + for (s=Real(0), k=i; k < m; k++) s += a(k,i)*a(k,j); + f=s/h; + for (k=i; k < m; k++) a(k,j) += f*a(k,i); + } + for (k=i; k < m; k++) a(k,i) *= scale; + } + } + + w[i] = scale*g; + g = Real(0); s = Real(0); scale = Real(0); + + if (i+1 <= m && i != n) { + for (k = l-1; k < n; k++) scale += std::abs(a(i,k)); + if (scale != Real(0)) { + for (k=l-1; k < n; k++) { + a(i,k) /= scale; + s += a(i,k)*a(i,k); + } + f = a(i, l-1); + g = -SIGN(std::sqrt(s), f); + h = f*g - s; + a(i, l-1) = f-g; + for (k=l-1; k < n; k++) rv1[k] = a(i,k)/h; + for (j = l-1; j < m; j++) { + for (s=Real(0), k=l-1; k < n; k++) s += a(j,k)*a(i,k); + for (k=l-1; k < n; k++) a(j,k) += s*rv1[k]; + } + for (k=l-1; k < n; k++) a(i,k) *= scale; + } + } + anorm = std::max(anorm, std::abs(w[i]) + std::abs(rv1[i])); + } + + for (i = n-1; i >= 0; i--) { // accumulation of right-hand transforms + if (i < n-1) { + if (g != Real(0)) { + for (j = l; j < n; j++) v(j,i) = (a(i,j)/a(i,l))/g; + for (j=l; j < n; j++) { + for (s=Real(0), k=l; k < n; k++) s += a(i,k)*v(k,j); + for (k=l; k < n; k++) v(k,j) += s*v(k,i); + } + } + for (j=l; j < n; j++) { + v(i,j) = Real(0); + v(j,i) = Real(0); + } + } + v(i,i) = Real(1); + g = rv1[i]; + l = i; + } + + for (i = std::min(m,n) -1; i >= 0; i--) { // accumulation of left-hand transforms + l = i+1; + g = w[i]; + for (j=l; j < n; j++) a(i,j) = Real(0); + if (g != Real(0)) { + g = Real(1)/g; + for (j=l; j < n; j++) { + for (s=Real(0), k=l; k < m; k++) s += a(k,i) * a(k,j); + f = (s/a(i,i))*g; + for (k=i; k < m; k++) a(k,j) += f*a(k,i); + } + for (j=i; j < m; j++) a(j,i) *= g; + } + else for (j=i; j < m; j++) a(j,i) = Real(0); + ++a(i,i); + } + + for (k=n-1; k >= 0; k--) { + for (its=0; its < 30; its++) { + flag = true; + for (l=k; l >= 0; l--) { + nm = l-1; + if (std::abs(rv1[l]) + anorm == anorm) { + flag = false; + break; + } + if (std::abs(w[nm]) + anorm == anorm) break; + } + if (flag) { + c = Real(0); + s = Real(1); + for (i=l-1; i < k+1; i++) { + f = s*rv1[i]; + rv1[i] = c*rv1[i]; + if (std::abs(f) + anorm == anorm) break; + g = w[i]; + h = pythag(f,g); + w[i] = h; + h = Real(1)/h; + c = g*h; + s = -f*h; + for (j=0; j < m; j++) { + y = a(j, nm); + z = a(j,i); + a(j, nm) = y*c + z*s; + a(j,i) = z*c - y*s; + } + } + } + + z = w[k]; + if (l == k) { + if (z < Real(0)) { + w[k] = -z; + for (j=0; j < n; j++) v(j,k) = -v(j,k); + } + break; + } + + if (its == 29) throw AbacusException("No convergence in 30 svdmp iterations"); + x = w[l]; + nm = k-1; + y = w[nm]; + g = rv1[nm]; + h = rv1[k]; + f = ((y-z) * (y+z) + (g-h) * (g+h))/(2*h*y); + g = pythag(f, Real(1)); + f = ((x-z)*(x+z) + h*((y/(f+SIGN(g,f)))-h))/x; + c = Real(1); s = Real(1); // Next QR transformation + for (j=l; j <= nm; j++) { + i = j+1; + g = rv1[i]; + y = w[i]; + h = s*g; + g = c*g; + z = pythag(f,h); + rv1[j] = z; + c = f/z; + s = h/z; + f = x*c + g*s; + g = g*c - x*s; + h = y*s; + y *= c; + for (jj = 0; jj < n; jj++) { + x = v(jj, j); + z = v(jj, i); + v(jj, j) = x*c + z*s; + v(jj, i) = z*c - x*s; + } + z = pythag(f,h); + w[j] = z; + if (z) { + z = Real(1)/z; + c = f*z; + s = h*z; + } + f = c*g + s*y; + x = c*y - s*g; + for (jj=0; jj < m; jj++) { + y = a(jj, j); + z = a(jj, i); + a(jj, j) = y*c + z*s; + a(jj, i) = z*c - y*s; + } + } + rv1[l] = Real(0); + rv1[k] = f; + w[k] = x; + } + } +} + + +// Householder reduction of a real symmetric matrix +// NR section 11.2 +export void tred2 (Matrix& a, std::vector& d, std::vector& e) +{ + int l, k, j, i; + Real scale, hh, h, g, f; + + int n { int(a.size()) }; + for (i = n-1; i > 0; i--) { + l = i - 1; + h = scale = Real(0); + if (l > 0) { + for (k = 0; k < l + 1; k++) scale += std::abs(a(i,k)); + if (scale == Real(0)) e[i] = a(i,l); + else { + for (k = 0; k < l + 1; k++) { + a(i,k) /= scale; + h += a(i,k) * a(i,k); + } + f = a(i,l); + g = (f >= Real(0) ? -std::sqrt(h) : std::sqrt(h)); + e[i] = scale * g; + h -= f * g; + a(i,l) = f - g; + f = Real(0); + for (j = 0; j < l + 1; j++) { + a(j,i) = a(i,j)/h; + g = Real(0); + for (k = 0; k < j + 1; k++) g += a(j,k) * a(i,k); + for (k = j + 1; k < l + 1; k++) g += a(k,j) * a(i,k); + e[j] = g/h; + f += e[j] * a(i,j); + } + hh = f/(h +h); + for (j = 0; j < l + 1; j++) { + f = a(i,j); + e[j] = g = e[j] - hh * f; + for (k = 0; k < j + 1; k++) a(j,k) -= (f * e[k] + g * a(i,k)); + } + } + } + else e[i] = a(i,l); + d[i] = h; + } + d[0] = Real(0); + e[0] = Real(0); + + for (i = 0; i < n; i++) { + l = i; + if (d[i] != Real(0)) { + for (j = 0; j < l; j++) { + g = Real(0); + for (k = 0; k < l; k++) g += a(i,k) * a(k,j); + for (k = 0; k < l; k++) a(k,j) -= g * a(k,i); + } + } + d[i] = a(i,i); + a(i,i) = Real(1); + for (j = 0; j < l; j++) a(j,i) = a(i,j) = Real(0); + } +} // tred2 + + + +// QL algorithm with implicit shifts +// NR section 11.3 +export void tqli (std::vector& d, std::vector& e, Matrix& z) +{ + int m, l, iter, i, k; + Real s, r, p, g, f, dd, c, b; + + int n { int(d.size()) }; + for (i = 1; i < n; i++) e[i-1] = e[i]; + e[n-1] = Real(0); + for (l = 0; l < n; l++) { + iter = 0; + do { + for (m = l; m < n-1; m++) { + dd = std::abs(d[m]) + std::abs(d[m+1]); + if (std::abs(e[m]) + dd == dd) break; + } + if (m != l) { + if (iter++ == 30) { + std::cout << "Too many iterations in tqli" << std::endl; + std::exit(1); + } + g = (d[l + 1] - d[l])/(2 * e[l]); + r = pythag(g, Real(1)); + g = d[m] - d[l] + e[l]/(g + SIGN(r, g)); + s = c = Real(1); + p = Real(0); + for (i = m - 1; i >= l; i--) { + f = s * e[i]; + b = c * e[i]; + e[i + 1] = (r = pythag(f,g)); + if (r == Real(0)) { + d[i + 1] -= p; + e[m] = Real(0); + break; + } + s = f/r; + c = g/r; + g = d[i + 1] - p; + r = (d[i] - g) * s + 2 * c * b; + d[i + 1] = g + (p = s * r); + g = c * r - b; + for (k = 0; k < n; k++) { + f = z(k, i + 1); + z(k, i + 1) = s * z(k, i) + c * f; + z(k, i) = c * z(k, i) - s * f; + } + } + if (r == Real(0) && i >= l) continue; + d[l] -= p; + e[l] = g; + e[m] = Real(0); + } + } while (m != l); + } +} // tqli + + + + +// Failed version using T** +// +// //////////////////// +// // ↓ Class Matrix // +// //////////////////// + +// export template +// class Matrix { + +// private: +// std::size_t dim_; +// T** element_; + +// public: +// Matrix (std::size_t dim); +// // ~Matrix () = default; +// ~Matrix (); + +// inline std::size_t size() { return dim_; }; + +// // inline T* operator[] (const IndexU i); +// // inline const T* operator[] (const IndexU i) const; + +// inline T operator() (IndexU i, IndexU j) const { return element_[i][j]; }; +// inline T& operator() (IndexU i, IndexU j) { return element_[i][j]; }; + +// void setZero(); + +// // Matrix& operator= (const Matrix& rhs); + +// Matrix& operator*= (const T& a); +// Matrix& operator/= (const T& a); + +// void ludcmp (std::vector& indx, T& d); +// void lubksb (std::vector& indx, std::vector& b); +// std::complex lndet_LU_destroy (); +// }; + +// template +// Matrix::Matrix (std::size_t dim) +// : dim_(dim) +// , element_(new T*[dim]) +// { +// // element_[0] = new T[dim_*dim_]; +// // for (IndexU i { 0 }; i + 1 < dim_; i++) element_[i+1] = element_[i] + dim_; +// for (IndexU i { 0 }; i < dim_; ++i) element_[i] = new T[dim]; +// } + +// template +// Matrix::~Matrix() +// { +// // if (element_ != 0) { +// // delete[] (element_[0]); +// // delete[] (element_); +// // } +// for (IndexU i { 0 }; i < dim_; ++i) delete[] element_[i]; +// delete[] element_; +// } + +// // template +// // inline T* Matrix::operator[] (const IndexU i) +// // { +// // return element_[i]; +// // } + +// // template +// // inline const T* Matrix::operator[] (const IndexU i) const +// // { +// // return element_[i]; +// // } + +// template +// void Matrix::setZero () +// { +// for (IndexU i { 0 }; i < dim_; ++i) +// for (IndexU j { 0 }; j < dim_; ++j) element_[i][j] = T(0); +// } + +// // template +// // Matrix& Matrix::operator= (const Matrix& rhs) +// // { +// // if (this != &rhs) { +// // if (dim_ != rhs.dim_) { +// // throw; +// // } + +// // for (int i = 0; i < dim_; ++i) +// // for (int j = 0; j < dim_; ++j) element_[i][j] = rhs.element_[i][j]; +// // } +// // return *this; +// // } + + + +// template +// Matrix& Matrix::operator*= (const T& a) +// { +// for (IndexU i { 0 }; i < dim_; ++i) +// for (IndexU j { 0 }; j < dim_; ++j) element_[i][j] *= a; +// return *this; +// } + +// template +// Matrix& Matrix::operator/= (const T& a) +// { +// T oneovera = T(1)/a; +// for (IndexU i { 0 }; i < dim_; ++i) +// for (IndexU j { 0 }; j < dim_; ++j) element_[i][j] *= oneovera; +// return *this; +// } + + +// template +// void Matrix::ludcmp (std::vector& indx, T& d) +// { +// IndexU i, j, k; +// IndexU imax { 0 }; +// T big, dum, sum, temp; + +// IndexU n { size() }; + +// std::vector vv(n); +// d = T(1); +// for (i = 0; i < n; i++) { +// big = T(0); +// for (j = 0; j < n; j++) { +// if ((std::fabs(temp = element_[i][j])) > std::fabs(big)) big = temp; +// } +// if (big == T(0)) throw; +// vv[i] = T(1)/big; +// } + +// for (j = 0; j < n; j++) { +// for (i = 0; i < j; i++) { +// sum = element_[i][j]; +// for (k = 0; k < i; k++) sum -= element_[i][k] * element_[k][j]; +// element_[i][j] = sum; +// } +// big = T(0); +// for (i = j; i < n; i++) { +// sum = element_[i][j]; +// for (k = 0; k < j; k++) sum -= element_[i][k] * element_[k][j]; +// element_[i][j] = sum; +// if ((std::fabs(dum = vv[i]*sum)) >= std::fabs(big)) { +// big = dum; +// imax = i; +// } +// } +// if (j != imax) { +// for (k = 0; k < n; k++) { +// dum = element_[imax][k]; +// element_[imax][k] = element_[j][k]; +// element_[j][k] = dum; +// } +// d = -d; +// vv[imax] = vv[j]; +// } +// indx[j] = imax; +// if (j !=n-1) { +// dum = T(1)/(element_[j][j]); +// for (i = j + 1; i < n; i++) element_[i][j] *= dum; +// } +// } +// } + +// template +// void Matrix::lubksb (std::vector& indx, std::vector& b) +// { +// int i, ip, j; +// int ii { 0 }; +// T sum; + +// int n { int(size()) }; +// for (i = 0; i < n; i++) { +// ip = indx[i]; +// sum = b[ip]; +// b[ip] = b[i]; +// if (ii != 0) +// for (j = ii-1; j < i; j++) sum -= element_[i][j] * b[j]; +// else if (sum != T(0)) +// ii = i + 1; +// b[i] = sum; +// } +// for (i = n - 1; i >= 0; i--) { +// sum = b[i]; +// for (j = i + 1; j < n; j++) sum -= element_[i][j] * b[j]; +// b[i] = sum/element_[i][i]; +// } +// } + +// template +// std::complex Matrix::lndet_LU_destroy () +// { +// std::vector indx(size()); +// T d; +// std::complex lndet { 0 }; + +// (*this).ludcmp(indx, d); + +// lndet = log(std::complex(d)); + +// for (IndexU j { 0 }; j < size(); j++) { +// lndet += log(std::complex(element_[j][j])); +// } + +// return lndet; +// } + + +// template +// std::ostream& operator<< (std::ostream& s, const Matrix& M) +// { +// for (IndexU i { 0 }; i < M.dim(); ++i) { +// for (IndexU j { 0 }; j < M.dim(); ++j) s << M[i][j] << "\t"; +// s << std::endl; +// } +// return s; +// } + +// //////////////////// +// // ↑ Class Matrix // +// //////////////////// + + + +// // Version using vector of vector (not working) + +// //////////////////// +// // ↓ Class Matrix // +// //////////////////// + +// export template +// class Matrix { + +// private: +// std::size_t dim_; +// std::vector> element_; + +// public: +// Matrix (std::size_t dim); +// ~Matrix () = default; + +// inline std::size_t size() { return dim_; }; +// inline std::vector> rows () const { return element_; }; +// inline std::vector& operator[] (const std::size_t i); +// inline T operator() (IndexU i, IndexU j) const { return element_[i][j]; }; +// inline T& operator() (IndexU i, IndexU j) { return element_[i][j]; }; + +// void setZero(); + +// Matrix& operator*= (const T& a); +// Matrix& operator/= (const T& a); + +// void ludcmp (std::vector& indx, T& d); +// void lubksb (std::vector& indx, std::vector& b); +// std::complex lndet_LU_destroy (); +// }; + +// template +// inline std::vector& Matrix::operator[] (const std::size_t i) +// { +// return element_[i]; +// } + +// template +// Matrix::Matrix (std::size_t dim) +// : dim_(dim) +// , element_(std::vector>(dim_, std::vector(dim_))) +// { +// } + +// template +// Matrix& Matrix::operator*= (const T& a) +// { +// for (auto &row : element_) +// for (auto &col : row) col *= a; +// return *this; +// } + +// template +// Matrix& Matrix::operator/= (const T& a) +// { +// for (auto &row : element_) +// for (auto &col : row) col /= a; +// return *this; +// } + +// template +// void Matrix::setZero () +// { +// for (auto &row : element_) +// for (auto &col : row) col = T(0); +// } + +// template +// void Matrix::ludcmp (std::vector& indx, T& d) +// { +// IndexU i, j, k; +// IndexU imax { 0 }; +// T big, dum, sum, temp; + +// IndexU n { size() }; + +// std::vector vv(n); +// d = T(1); +// for (i = 0; i < n; i++) { +// big = T(0); +// for (j = 0; j < n; j++) { +// if ((std::fabs(temp = element_[i][j])) > std::fabs(big)) big = temp; +// } +// if (big == T(0)) throw; +// vv[i] = T(1)/big; +// } + +// for (j = 0; j < n; j++) { +// for (i = 0; i < j; i++) { +// sum = element_[i][j]; +// for (k = 0; k < i; k++) sum -= element_[i][k] * element_[k][j]; +// element_[i][j] = sum; +// } +// big = T(0); +// for (i = j; i < n; i++) { +// sum = element_[i][j]; +// for (k = 0; k < j; k++) sum -= element_[i][k] * element_[k][j]; +// element_[i][j] = sum; +// if ((std::fabs(dum = vv[i]*sum)) >= std::fabs(big)) { +// big = dum; +// imax = i; +// } +// } +// if (j != imax) { +// for (k = 0; k < n; k++) { +// dum = element_[imax][k]; +// element_[imax][k] = element_[j][k]; +// element_[j][k] = dum; +// } +// d = -d; +// vv[imax] = vv[j]; +// } +// indx[j] = imax; +// if (j !=n-1) { +// dum = T(1)/(element_[j][j]); +// for (i = j + 1; i < n; i++) element_[i][j] *= dum; +// } +// } +// } + +// template +// void Matrix::lubksb (std::vector& indx, std::vector& b) +// { +// int i, ip, j; +// int ii { 0 }; +// T sum; + +// int n { int(size()) }; +// for (i = 0; i < n; i++) { +// ip = indx[i]; +// sum = b[ip]; +// b[ip] = b[i]; +// if (ii != 0) +// for (j = ii-1; j < i; j++) sum -= element_[i][j] * b[j]; +// else if (sum != T(0)) +// ii = i + 1; +// b[i] = sum; +// } +// for (i = n - 1; i >= 0; i--) { +// sum = b[i]; +// for (j = i + 1; j < n; j++) sum -= element_[i][j] * b[j]; +// b[i] = sum/element_[i][i]; +// } +// } + +// template +// std::complex Matrix::lndet_LU_destroy () +// { +// std::vector indx(size()); +// T d; +// std::complex lndet { 0 }; + +// (*this).ludcmp(indx, d); + +// lndet = log(std::complex(d)); + +// for (IndexU j { 0 }; j < size(); j++) { +// lndet += log(std::complex(element_[j][j])); +// } + +// return lndet; +// } + + +// template +// std::ostream& operator<< (std::ostream& s, const Matrix& M) +// { +// for (auto row : M.rows()) { +// for (auto col : row) s << col << "\t"; +// s << std::endl; +// } +// return s; +// } + +// //////////////////// +// // ↑ Class Matrix // +// //////////////////// diff --git a/src/util/timer.cc b/src/util/timer.cc new file mode 100644 index 0000000..4a82a9d --- /dev/null +++ b/src/util/timer.cc @@ -0,0 +1,55 @@ +/***************************************************************** + +This software is part of Jean-Sébastien Caux's Abacus toolsuite. + +Copyright © Jean-Sébastien Caux. + +*****************************************************************/ + + +///////////// +// ↓ Timer // +///////////// + +// Timer utility from https://www.learncpp.com/cpp-tutorial/timing-your-code/ + + +export module timer; + + +import std; + + +export class Timer +{ +private: + // Type aliases to make accessing nested type easier + using Clock = std::chrono::steady_clock; + using Second = std::chrono::duration>; + + std::chrono::time_point m_beg { Clock::now() }; + +public: + void reset () + { + m_beg = Clock::now(); + } + + double elapsed () const + { + return std::chrono::duration_cast(Clock::now() - m_beg).count(); + } +}; + + +export std::string timestamp () +{ + std::time_t current_time { std::time(nullptr) }; + char timestr[100]; + std::strftime(timestr, sizeof(timestr), "%Y-%m-%d %H:%M:%S", std::gmtime(¤t_time)); + return timestr; +} + +///////////// +// ↑ Timer // +/////////////