From 8fcb948ecca7985b0b9deeda212f516b6accb018 Mon Sep 17 00:00:00 2001 From: "matthias@quintern.xyz" Date: Sun, 26 Jan 2025 19:58:40 +0100 Subject: [PATCH] add dmrg --- .gitignore | 17 + .gitmodules | 6 + lambda-lanczos | 1 + spectra | 1 + src/.clangd | 2 +- src/.doxygen_config | 4 +- src/Makefile | 70 ++-- src/algorithm.cpp | 3 + src/algorithm.hpp | 20 + src/algorithms/dmrg.cpp | 288 +++++++++++++ src/algorithms/dmrg.hpp | 52 +++ src/algorithms/tebd.cpp | 201 ++++----- src/algorithms/tebd.hpp | 65 ++- src/benchmark.cpp | 5 +- src/main.cpp | 122 ++++-- src/model.cpp | 2 +- src/model.hpp | 5 +- src/models/bose_hubbard.cpp | 44 +- src/models/bose_hubbard.hpp | 38 +- src/models/transverse_field_ising.cpp | 22 +- src/models/transverse_field_ising.hpp | 22 +- src/mpl_colors.py | 154 ------- src/mps.cpp | 32 +- src/mps.hpp | 85 +++- src/output | 562 -------------------------- src/run.cpp | 390 ------------------ src/type.hpp | 2 + src/wasm/emscripten.hpp | 27 ++ src/wasm/range.hpp | 88 ++++ src/wasm/result.hpp | 81 ++++ src/wasm/run.cpp | 107 +++++ src/wasm/run.hpp | 166 ++++++++ src/wasm/util.hpp | 0 33 files changed, 1322 insertions(+), 1362 deletions(-) create mode 100644 .gitignore create mode 100644 .gitmodules create mode 160000 lambda-lanczos create mode 160000 spectra create mode 100644 src/algorithm.cpp create mode 100644 src/algorithm.hpp create mode 100644 src/algorithms/dmrg.cpp create mode 100644 src/algorithms/dmrg.hpp delete mode 100644 src/mpl_colors.py delete mode 100644 src/output delete mode 100644 src/run.cpp create mode 100644 src/wasm/emscripten.hpp create mode 100644 src/wasm/range.hpp create mode 100644 src/wasm/result.hpp create mode 100644 src/wasm/run.cpp create mode 100644 src/wasm/run.hpp create mode 100644 src/wasm/util.hpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..17afe23 --- /dev/null +++ b/.gitignore @@ -0,0 +1,17 @@ +*.pdf +*.html +*.js +*benchmarks +*build* +docs +include +*.so +*.py +mps +src-alloc-debug/ +src-eigen-opt/ +src/.old +src/pch.hpp.gch +test-python/ +web/ + diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..a1b4d1f --- /dev/null +++ b/.gitmodules @@ -0,0 +1,6 @@ +[submodule "lambda-lanczos"] + path = lambda-lanczos + url = https://github.com/mrcdr/lambda-lanczos +[submodule "src/spectra"] + path = spectra + url = https://github.com/yixuan/spectra diff --git a/lambda-lanczos b/lambda-lanczos new file mode 160000 index 0000000..d45badd --- /dev/null +++ b/lambda-lanczos @@ -0,0 +1 @@ +Subproject commit d45baddc4239dd1de3f0091ce495920b57a3fd5e diff --git a/spectra b/spectra new file mode 160000 index 0000000..a29f37f --- /dev/null +++ b/spectra @@ -0,0 +1 @@ +Subproject commit a29f37fe3ed1c2895b07b449fcb8f09bc0940e40 diff --git a/src/.clangd b/src/.clangd index 15a9781..59f5e49 100644 --- a/src/.clangd +++ b/src/.clangd @@ -1,3 +1,3 @@ CompileFlags: # Tweak the parse settings - Add: [ -std=c++23, -Wall, -Wpedantic, -IFastor, -I., -I.., -I../include, -Imodels, -I../models, -Ialgorithms, -DEMSCRIPTEN ] + Add: [ -std=c++23, -Wall, -Wpedantic, -IFastor, -I., -I.., -Imodels, -I../models, -Ialgorithms, -iprefix ~/Projekte/science/cpp-mps/ -iwithprefix include -iwithprefix lambda-lanczos/include -iwithprefix spectra/include -DEMSCRIPTEN ] # https://clangd.llvm.org/config diff --git a/src/.doxygen_config b/src/.doxygen_config index c86914e..235674e 100644 --- a/src/.doxygen_config +++ b/src/.doxygen_config @@ -965,7 +965,7 @@ WARN_LOGFILE = # spaces. See also FILE_PATTERNS and EXTENSION_MAPPING # Note: If this tag is empty the current directory is searched. -INPUT = . +INPUT = . models algorithms wasm # This tag can be used to specify the character encoding of the source files # that Doxygen parses. Internally Doxygen uses the UTF-8 encoding. Doxygen uses @@ -2087,7 +2087,7 @@ PAPER_TYPE = a4 # If left blank no extra packages will be included. # This tag requires that the tag GENERATE_LATEX is set to YES. -EXTRA_PACKAGES = {amsmath} +EXTRA_PACKAGES = {amsmath} {bbm} # The LATEX_HEADER tag can be used to specify a user-defined LaTeX header for # the generated LaTeX document. The header should contain everything until the diff --git a/src/Makefile b/src/Makefile index bc0af1a..a158d64 100755 --- a/src/Makefile +++ b/src/Makefile @@ -9,19 +9,34 @@ WEBEXEC = ~/Projekte/web/html-quintern/src/script/mps/mps-wasm.js # LOG_LEVEL = LOG_LEVEL_0 CXX = /usr/bin/g++ CXXFLAGS = -std=c++23 -MMD -MP -Wpedantic -Wall #-O3: optimierungsstufe 3, -MMD .d files, -MP leeres Target, Wextra alle Warnungen -CXXDEBUGFLAGS = -g -fsanitize=address -fno-omit-frame-pointer -O1 +CXXDEBUGFLAGS = -g -fno-omit-frame-pointer -O1 +# CXXDEBUGFLAGS += -fsanitize=address WEBCXX = em++ -WEBCXXFLAGS = --cache "$(WEBOBJECT_DIR)/emscripten-cache" -Wno-c++11-narrowing -Wno-unused-main -DEMSCRIPTEN -WEBCXXDEBUGFLAGS = -gsource-map -g2 -O0 -WEBLDFLAGS = -s EXPORTED_FUNCTIONS="['_free']" -s ASYNCIFY -s EXPORTED_RUNTIME_METHODS="['cwrap', 'getValue']" -lembind -s INITIAL_MEMORY=2GB -s ALLOW_MEMORY_GROWTH=1 +WEBCXXFLAGS = --cache "$(WEBOBJECT_DIR)/emscripten-cache" -Wno-c++11-narrowing -Wno-unused-main -DEMSCRIPTEN +# the TEST___FLAGS will be made available with the getBuildOptions() function in webassembly for benchmarking purposes +# TESTCXXFLAGS = -g -ffast-math +# TESTCXXFLAGS = -g -fno-math-errno +TESTCXXFLAGS = -msimd128 # concurrency, TODO: check if does something +# TESTCXXFLAGS += -g +# TESTCXXFLAGS += -s INVOKE_RUN=0 +# TESTLDFLAGS = -ffast-math +# TESTLDFLAGS += --closure 1 +TESTLDFLAGS = -s ASSERTIONS=1 -s NO_DISABLE_EXCEPTION_CATCHING -g2 +WEBCXXDEBUGFLAGS = -gsource-map -g3 -O0 -s ASSERTIONS=2 +WEBLDFLAGS = -lembind -s ASYNCIFY -s EXPORTED_RUNTIME_METHODS="['cwrap', 'getValue']" -s INITIAL_MEMORY=2GB -s ALLOW_MEMORY_GROWTH=1 +# WEBLDFLAGS += -s STACK_SIZE=131072 -sINVOKE_RUN=0 #-sSTACK_OVERFLOW_CHECK=1 #-s ASSERTIONS=2 # TODO remove +WEBLDFLAGS += -s STACK_SIZE=331072 -sINVOKE_RUN=0 #-sSTACK_OVERFLOW_CHECK=1 #-s ASSERTIONS=2 # TODO remove + +WEBLDRELEASEFLAGS = -sDISABLE_EXCEPTION_CATCHING # treated as error by em++ -LDFLAGS = #-lgzutil -IFLAGS = -I../include/ +LDFLAGS = +# IFLAGS = -I../include/ -I../../lambda-lanczos/include -I../spectra/include +IFLAGS = -iprefix $$HOME/Projekte/science/cpp-mps/ -iwithprefix include -iwithprefix lambda-lanczos/include -iwithprefix spectra/include SRCDIRS = $(wildcard */) IFLAGS += $(foreach dir,$(SRCDIRS), -I$(dir)) -IFLAGS += $(foreach dir,$(SRCDIRS), -I../$(dir)) +# IFLAGS += $(foreach dir,$(SRCDIRS), -I../$(dir)) IFLAGS += -I../ -I./ CXXFLAGS += $(IFLAGS) @@ -33,39 +48,43 @@ CXXFLAGS += $(IFLAGS) SRC = $(wildcard *.cpp) $(wildcard */*.cpp) # OBJECTS = $(SRC:%.cpp=$(OBJECT_DIR)/%.o) OBJECTS = $($(notdir SRC):%.cpp=$(OBJECT_DIR)/%.o) -WEBOBJECTS = $($(notdir SRC):%.cpp=$(WEBOBJECT_DIR)/%.o) OBJECT_DIRS = $(foreach dir,$(SRCDIRS), $(OBJECT_DIR)/$(dir)) -WEBOBJECT_DIRS = $(foreach dir,$(SRCDIRS), $(WEBOBJECT_DIR)/$(dir)) DEPENDS = ${OBJECTS:.o=.d} +WEBOBJECTS = $($(notdir SRC):%.cpp=$(WEBOBJECT_DIR)/%.o) +WEBOBJECT_DIRS = $(foreach dir,$(SRCDIRS), $(WEBOBJECT_DIR)/$(dir)) +WEBDEPENDS = ${WEBOBJECTS:.o=.d} + FMT_MESSAGE="\e[1;34m%s\e[0m %s %s %s\n" +# default: CXXFLAGS += -include pch.hpp + +# include the makefiles generated by the -M flag +-include $(DEPENDS) +-include $(WEBDEPENDS) + default: @printf $(FMT_MESSAGE) "CXXFLAGS = " "$(CXXFLAGS)" @printf $(FMT_MESSAGE) "LDFLAGS = " "$(LDFLAGS)" -default: CXXFLAGS += -include pch.hpp default: $(EXEC) -web: CXXFLAGS += $(WEBCXXFLAGS) +web: CXXFLAGS += $(WEBCXXFLAGS) $(TESTCXXFLAGS) $(TESTLDFLAGS) -D BUILD_OPTIONS_STRING="\"$(TESTCXXFLAGS) $(TESTLDFLAGS)\"" +# web: LDFLAGS += $(WEBLDFLAGS) $(TESTLDFLAGS) web: $(WEBEXEC) webDebug: WEBLDFLAGS += -sASSERTIONS webDebug: CXXFLAGS += $(WEBCXXDEBUGFLAGS) webDebug: web -webO0: CXXFLAGS += $(WEBCXXFLAGS) -O0 +webO0: TESTCXXFLAGS += -O0 webO0: web -webO1: CXXFLAGS += $(WEBCXXFLAGS) -O1 +webO1: TESTCXXFLAGS += -O1 webO1: web -webO2: CXXFLAGS += $(WEBCXXFLAGS) -O2 +webO2: TESTCXXFLAGS += -O2 webO2: web -webO3: CXXFLAGS += $(WEBCXXFLAGS) -O3 +webO3: TESTCXXFLAGS += -O3 webO3: web -webrun: CXXFLAGS += $(WEBCXXFLAGS) -webrun: $(WEBHTML) - cd $(WEBOBJECT_DIR) && emrun $(WEBHTML) - .PHONY: release release: CXXFLAGS += -O3 release: default @@ -84,9 +103,7 @@ $(EXEC): $(OBJECT_DIRS) $(OBJECTS) $(CXX) $(OBJECTS) -o $@ $(CXXFLAGS) $(LDFLAGS) $(WEBHTML) $(WEBEXEC): $(WEBOBJECT_DIRS) $(WEBOBJECTS) - $(WEBCXX) $(WEBOBJECTS) -o $@ $(CXXFLAGS) $(LDFLAGS) $(WEBLDFLAGS) -# include the makefiles generated by the -M flag --include $(DEPENDS) + $(WEBCXX) $(WEBOBJECTS) -o $@ $(CXXFLAGS) $(LDFLAGS) $(WEBLDFLAGS) # rule for all ../build/*.o files $(OBJECT_DIR)/%.o: $(shell echo $<) %.cpp @@ -106,7 +123,13 @@ $(OBJECT_DIRS) $(WEBOBJECT_DIRS): # with debug flags .PHONY += debug run gdb pch clean docs -debug: CXXFLAGS += $(CXXDEBUGFLAGS) +MASSIF_OUT=../massif-mps.out +profile: $(EXEC) + valgrind --tool=massif --stacks=yes --massif-out-file=$(MASSIF_OUT) $(EXEC) + ms_print $(MASSIF_OUT) | less + +debug: CXXFLAGS += $(CXXDEBUGFLAGS) +# debug: LDFLAGS += -static-libasan debug: default # make with debug flags and run afterwards @@ -126,6 +149,7 @@ pch: # remove all object and dependecy files clean: -rm -r $(OBJECT_DIR) +webclean: -rm -r $(WEBOBJECT_DIR) docs: doxygen .doxygen_config diff --git a/src/algorithm.cpp b/src/algorithm.cpp new file mode 100644 index 0000000..550b002 --- /dev/null +++ b/src/algorithm.cpp @@ -0,0 +1,3 @@ +#include "algorithm.hpp" +template bool postUpdateNoop(const MPS&, unsigned); + diff --git a/src/algorithm.hpp b/src/algorithm.hpp new file mode 100644 index 0000000..0e40ed5 --- /dev/null +++ b/src/algorithm.hpp @@ -0,0 +1,20 @@ +#pragma once +#include "mps.hpp" +#include "type.hpp" + +#include + +/** + * @brief Type for a function to after one algorithm iteration. + * @details + * The return value states wether the results have converged, thus + * if the function returns `true`, the simulation will stop. + */ +template +using postUpdateFunction_t = std::function&, unsigned)>; + +/** + * @brief Default post update function that does nothing + */ +template +bool postUpdateNoop(const MPS&, unsigned) { return false; }; diff --git a/src/algorithms/dmrg.cpp b/src/algorithms/dmrg.cpp new file mode 100644 index 0000000..a1a7597 --- /dev/null +++ b/src/algorithms/dmrg.cpp @@ -0,0 +1,288 @@ +#include "dmrg.hpp" + +#include "lanczos.hpp" +#include "type.hpp" +#include "util.hpp" + + +#include +using lambda_lanczos::LambdaLanczos; + +#include +#include +#include +#include + +namespace dmrg { + void HeffTheta( + const e::Tensor& LP, // vL wL* vL* + const e::Tensor& RP, // vR* wR* vR + const e::Tensor& W1, // wL wC i i* + const e::Tensor& W2, // wC wR j j* + const std::vector& thetaIn, + std::vector& thetaOut + ) { + assert(thetaIn.size() == LP.dimension(0) * W1.dimension(2) * W2.dimension(2) * RP.dimension(2)); + assert(thetaOut.size() == LP.dimension(0) * W1.dimension(2) * W2.dimension(2) * RP.dimension(2)); + e::array thetaShape = {LP.dimension(0), W1.dimension(2), W2.dimension(2), RP.dimension(2)}; + // Reshape theta into a 4D tensor + e::TensorMap> thetaTensor(thetaIn.data(), thetaShape); // vL i j vR + // e::Tensor x = theta_tensor; + auto LPT = LP.contract(thetaTensor, e::array({iP{2, 0}})); // vL wL* [vL*], [vL] i j vR + auto LPTW1 = LPT.contract(W1, e::array({iP{1, 0}, iP{2, 3}})); // vL [wL*] [i] j vR, [wL] wC i [i*] + auto LPTW1W2 = LPTW1.contract(W2, e::array({iP{3, 0}, iP{1, 3}})); // vL [j] vR [wC] i, [wC] wR j [j*] + auto LPTW1W2RP = LPTW1W2.contract(RP, e::array({iP{1, 0}, iP{3, 1}})); // vL [vR] i [wR] j, [vR*] [wR*] vR + + // assigning to map will assign to underlying std::vector + e::TensorMap> thetaOutMap(thetaOut.data(), thetaShape); // vL i j vR + thetaOutMap = LPTW1W2RP; + // std::println("thetaO: {}", thetaOut); + } + + + template + DMRGEngine::DMRGEngine(MPS& psi, const std::vector>>& H_mpo, + unsigned chiMax, T eps) + : H_mpo(H_mpo), chiMax(chiMax), eps(eps) { + int L = psi.getL(); + LPs.resize(L); + RPs.resize(L); + initialize(psi); + } + + template + void DMRGEngine::initialize(MPS& psi) { + // Initialize LPs and RPs + int D = H_mpo[0]->dimension(0); + int chi = psi.B(0).dimension(0); + e::Tensor LP(chi, D, chi), RP(chi, D, chi); + LP.setZero(); + RP.setZero(); + LP.chip(0, 1) = identityTensor(chi, chi); + RP.chip(D - 1, 1) = identityTensor(chi, chi); + + LPs[0] = LP; + RPs.back() = RP; + + for (int i = psi.getL() - 1; i > 1; --i) { + update_RP(i, psi); + } + } + + template + void DMRGEngine::sweep(MPS& psi) { + // Sweep left to right + for (unsigned i = 0; i < psi.getL() - 2; ++i) { + update_bond(i, psi); + } + // Sweep right to left + for (unsigned i = psi.getL() - 2; i > 0; --i) { + update_bond(i, psi); + } + } + + template + void DMRGEngine::update_bond(int i, MPS& psi) { + int j = i + 1; + // std::println("update({})", i); + + e::Index d1 = H_mpo[i]->dimension(2); + e::Index d2 = H_mpo[j]->dimension(2); + e::Index chi1 = LPs[i].dimension(0); + e::Index chi2 = RPs[j].dimension(2); + e::Index dim = chi1*d1*d2*chi2; + + // 1) Find GS (smallest eigenvalue / vector) + // Approaches 1 and 2 are wrong, because we need a more complicated contraction series + // instead of a reshape + simple matrix product for Heff @ theta + + // For i) and ii): create effective hamiltonian + // auto LW = LPs[i].contract(*H_mpo[i], e::array{iP{1, 0}}); // vL [wL] vL*, [wL] wR i i* + // auto LWW = LW.contract(*H_mpo[j], e::array{iP{2, 1}}); // vL vL* [wR] i i*, [wL] wR j j* + // auto LWWR = LWW.contract(RPs[j], e::array{iP{4, 1}}); // vL vL* i i* [wR] j j*, vR* [wR*] vR + // e::Tensor Heff = LWWR;// vL vL* i i* j j*, vR* vR + // printTensorDims(Heff, "Heff"); + + // printTensor(LPs[i], "LP[i]"); + // printTensor(RPs[j], "RP[j]"); + // printTensor(*H_mpo[i], "H_mpo[i]"); + // printTensor(*H_mpo[j], "H_mpo[j]"); + // e::Map> HeffMat(Heff.data(), dim, dim); // vL vL* i i*, j j* vR* vR + // std::cout << "HeffMat" << std::endl; + // std::cout << HeffMat << std::endl; + + // i) Eigen + // Realllyyyy inefficient, it calculates all eigenvalues and we just need the eigenvector to the + // smallest absolute eigenvalue + // e::SelfAdjointEigenSolver> es(HeffMat); + // const e::VectorX& S = es.eigenvalues(); + // const e::MatrixX& U = es.eigenvectors(); + // std::println("Eigenvalues with eigen:"); + // std::cout << S << std::endl; + // std::println("Eigenvector with eigen:"); + // std::cout << U.rows() << "x" << U.cols() << std::endl; + // e::Map> U1(U.col(0).data(), chi1*d1, d2*chi2); + // std::cout << "theta2 (eigen)" << std::endl; + // std::cout << U1 << std::endl; + // find the smallest eigenvalue + // e::TensorMap> theta(U.data(), chi1, d1, d2, chi2); + // auto [Ai, Sj, Bj] = split_truncate_theta(theta, chiMax, eps); + // auto [Ai, Sj, Bj] = split_truncate_theta(U1, chi1, d1, d2, chi2, chiMax, eps); + + // ii) Lanczos + // std::println("Eigenvector with Lanczos:"); + // auto [eval, theta2Vec] = lanczos(HeffMat, psi.getTheta2(i)); + // e::Map> theta2(theta2Vec.data(), chi1*d1, d2*chi2); + // std::cout << theta2.rows() << "x" << theta2.cols() << std::endl; + // std::cout << eval << std::endl; + // std::cout << theta2 << std::endl; + // auto [Ai, Sj, Bj] = split_truncate_theta(theta2, chi1, d1, d2, chi2, chiMax, eps); + // auto theta0 = psi.getTheta2(i); + // printTensor(theta0, "theta0"); + // e::Map> theta0v(theta0.data(), theta0.size()); + // std::println("Heff @ theta0"); + // e::VectorX v = HeffMat * theta0v; + // std::cout << v << std::endl; + + // iii) Lanczos with custom contraction function + e::Tensor theta0 = psi.getTheta2(i); + // Use Theta2 as the initial guess for the GS + auto initF = [&theta0](std::vector& v) { + v.resize(theta0.size()); + for (unsigned i = 0; i < theta0.size(); i++) { + v.at(i) = theta0.data()[i]; // TODO find nicer solution + } + // v = std::vector(theta0.data(), theta0.size()); + }; + auto matmul = [this, i, j](const std::vector& thetaIn, std::vector& thetaOut) { + HeffTheta(this->LPs.at(i), this->RPs.at(j), *(this->H_mpo.at(i)), *(this->H_mpo.at(j)), thetaIn, thetaOut); + }; + + // TODO set eps + LambdaLanczos engine(matmul, dim, false, 1); // Find 1 maximum eigenvalue + engine.init_vector = initF; + std::vector eigenvalues; + std::vector> eigenvectors; + engine.run(eigenvalues, eigenvectors); + e::Map> theta2(eigenvectors.front().data(), chi1*d1, d2*chi2); + // std::println("Lanczos ev: {}", eigenvalues); + // e::TensorMap> theta2T(eigenvectors.front().data(), chi1, d1, d2, chi2); + // printTensor(theta2T, "theta2"); + + auto [Ai, Sj, Bj] = split_truncate_theta(theta2, chi1, d1, d2, chi2, chiMax, eps); + + + + // 2) Truncate + // printTensor(Ai, "Ai"); + // printTensor(Sj, "Sj"); + // printTensor(Bj, "Bj"); + // std::cout << "a: " << a.dimensions() << "\nA: " << Ai.dimensions() << "\nS: " << Sj.dimensions() << "\nB: " << Bj.dimensions() << std::endl; + // renormalize to Bi + e::Tensor SiInv = 1. / psi.S(i); + e::Tensor SiDiv = diagonal(SiInv); + e::Tensor Gi = SiDiv.contract(Ai, e::array{iP{1, 0}}); // vL [vL], [vL] i vC + // printTensor(SiDiv, "SiDiv"); + // printTensor(Gi, "Gi"); + + // e::Tensor SS = Sj.contract(Sj.conjugate(), e::array{iP{0,0}}); + // std::println("Bond {:03}, Sj contracted with its conjugate = {}", i, SS(0)); + + // psi.B(i) = Gi.contract(diagonal(Sj), e::array{iP{2, 0}}); // vL i [vC], [vC] vC + // psi.B(i+1) = std::move(Bj); + // psi.S(i+1) = std::move(Sj); + typename MPS::BTensor_t Bi = Gi.contract(diagonal(Sj), e::array{iP{2, 0}}); // vL i [vC], [vC] vC + // e::Tensor equal = (psi.B(i) == Bi).all(); + // printTensor(psi.B(i), "old Bi"); + // printTensor(Bi, "new Bi"); + psi.B(i) = std::move(Bi); + psi.B(j) = std::move(Bj); + psi.S(j) = std::move(Sj); + + // 3) Update + update_LP(i, psi); + update_RP(j, psi); + } + + template + void DMRGEngine::update_RP(int i, MPS& psi) { + // std::println("update_RP({})", i); + int j = i - 1; + const e::Tensor& RP = RPs[i]; // vR* wR* vR + const e::Tensor& B = psi.B(i); // vL i vR + e::Tensor Bc = B.conjugate(); // vL* i* vR* + const e::Tensor& W = *H_mpo[i]; // wL wR i i* +#ifdef DMRGDEBUG + printTensorDims(RP, "RP"); + printTensorDims(B, "B"); + printTensorDims(Bc, "Bc"); + printTensorDims(W, "W"); + + auto RP2 = B.contract(RP, e::array{iP{2, 0}}); // vL i [vR], [vR*] wR* vR + e::Tensor RP2T = RP2; + printTensorDims(RP2T, "RP2T"); + auto RP3 = RP2T.contract(W, e::array{iP{1, 3}, iP{2, 1}}); // vL [i] [wR*] vR, wL [wR] i [i*] + e::Tensor RP3T = RP3; + printTensorDims(RP3T, "RP3T"); + auto RP4 = RP3.contract(Bc, e::array{iP{1, 2}, iP{3, 1}}); // vL [vR] wL [i], vL* [i*] [vR*] + e::Tensor RP4T = RP4; + printTensor(RP4T, "RP4T"); +#else + auto RP2 = B.contract(RP, e::array{iP{2, 0}}); // vL i [vR], [vR*] wR* vR + auto RP3 = RP2.contract(W, e::array{iP{1, 3}, iP{2, 1}}); // vL [i] [wR*] vR, wL [wR] i [i*] + auto RP4 = RP3.contract(Bc, e::array{iP{1, 2}, iP{3, 1}}); // vL [vR] wL [i], vL* [i*] [vR*] +#endif + RPs[j] = RP4; + } + + template + void DMRGEngine::update_LP(int i, MPS& psi) { + // std::println("update_LP({})", i); + int j = i + 1; + const e::Tensor& LP = LPs[i]; // vL wL vL* + const e::Tensor& B = psi.B(i); // vL i vR + const e::Tensor& W = *H_mpo[i]; // wL wR i i* + + e::Tensor SjInvV = 1. / psi.S(j); + e::Tensor SjInv = diagonal(SjInvV); + auto G = B.contract(SjInv, e::array{iP{2, 0}}); // vL i [vR], [vR*] vR + e::Tensor Si = diagonal(psi.S(i)); + auto A = Si.contract(G, e::array{iP{1, 0}}); // vL [vL*], [vL] i vR + auto Ac = A.conjugate(); // vL* i* vR* + +#ifdef DMRGDEBUG + auto LP2 = LP.contract(A, e::array{iP{2, 0}}); // vL wL* [vL*], [vL] i vR + e::Tensor LP2T = LP2; + printTensorDims(LP2T, "LP2T"); + printTensorDims(W, "W"); + auto LP3 = W.contract(LP2T, e::array{iP{0, 1}, iP{3, 2}}); // [wL] wR i [i*], vL [wL*] [i] vR + e::Tensor LP3T = LP3; + printTensorDims(LP3T, "LP3T"); + e::Tensor AcT = Ac; + printTensorDims(AcT, "AcT"); + auto LP4 = AcT.contract(LP3, e::array{iP{0, 2}, iP{1, 1}}); // [vL*] [i*] vR*, wR [i] [vL] vR + e::Tensor LP4T = LP4; + printTensor(LP4T, "LP4T"); + LPs[j] = LP4T; +#else + auto LP2 = LP.contract(A, e::array{iP{2, 0}}); // vL wL* [vL*], [vL] i vR + auto LP3 = W.contract(LP2, e::array{iP{0, 1}, iP{3, 2}}); // [wL] wR i [i*], vL [wL*] [i] vR + auto LP4 = Ac.contract(LP3, e::array{iP{0, 2}, iP{1, 1}}); // [vL*] [i*] vR*, wR [i] [vL] vR + LPs[j] = LP4; +#endif + } + + template + void run(MPS& psi, const std::vector>>& H_mpos, unsigned chiMax, T eps, unsigned nSweeps, postUpdateFunction_t postUpdateFunc) { + auto engine = DMRGEngine(psi, H_mpos, chiMax, eps); + for (unsigned i = 0; i < nSweeps; i++) { + engine.sweep(psi); + if (postUpdateFunc(psi, i)) break; + } + } + + template class dmrg::DMRGEngine; + template void run(MPS& psi, const std::vector>>& H_mpos, unsigned chiMax, dtype eps, unsigned nSweeps, postUpdateFunction_t postUpdateFunc); + +} // namespace dmrg + diff --git a/src/algorithms/dmrg.hpp b/src/algorithms/dmrg.hpp new file mode 100644 index 0000000..59d675b --- /dev/null +++ b/src/algorithms/dmrg.hpp @@ -0,0 +1,52 @@ +#pragma once +#include "algorithm.hpp" +#include "mps.hpp" +#include + + +/** + * @brief Returns a \f$L\times L \f$ identity matrix in tensor form, where \f$L = \text{min}(N, M)\f$ + */ +template +e::Tensor identityTensor(e::Index N, e::Index M) { + e::Tensor t(N, M); + t.setZero(); + e::Index min = std::min(N, M); + for (e::Index i = 0; i < min; i++) { + t(i, i) = 1; + } + return t; +} + + +/** + * @brief Implementation of the Density Matrix Renormalization Group (DMRG) algorithm + */ +namespace dmrg { + /** + * @brief Helper class that keeps the left and right blocks of the MPS + * @note Always call with the same mps (psi) that it was initialized with! + * It would have been better if the engine would take ownership of the MPS during sweeps, + * but this would make calling an update function in between more difficult. + */ + template + class DMRGEngine { + public: + std::vector>> H_mpo; + std::vector> LPs, RPs; + int chiMax; + T eps; + + DMRGEngine(MPS& psi, const std::vector>>& H_mpo, + unsigned chiMax = 100, T eps = 1e-12); + + void initialize(MPS& psi); + void sweep(MPS& psi); + void update_bond(int i, MPS& psi); + void update_RP(int i, MPS& psi); + void update_LP(int i, MPS& psi); + }; + + template + void run(MPS& psi, const std::vector>>& H_mpos, unsigned chiMax, T eps, unsigned nSweeps, postUpdateFunction_t postUpdateFunc); +} diff --git a/src/algorithms/tebd.cpp b/src/algorithms/tebd.cpp index b11c61e..e1a2971 100644 --- a/src/algorithms/tebd.cpp +++ b/src/algorithms/tebd.cpp @@ -1,10 +1,11 @@ #include "tebd.hpp" #include "type.hpp" #include "util.hpp" +#include "model.hpp" #include -#include #include #include +#include template std::vector> forEachUnique(std::function f, const std::vector>& inRange) { @@ -25,107 +26,109 @@ std::vector> forEachUnique(std::function -e::Tensor getTimeEvolutionMatrix(const e::MatrixX& bond, T dt, unsigned localDim) { - e::MatrixX U = (-bond * dt).exp(); - // std::cout << "U = \n" << U << "\n\n"; - return bondMatrix2Tensor(U, localDim); -} - -template -Bonds getTimeEvolutionMatrices(const std::vector>>& bonds, T dt, unsigned int localDim) { - std::function(const e::MatrixX&)> f = std::bind(getTimeEvolutionMatrix, std::placeholders::_1, dt, localDim); - std::vector>> U_Bonds = forEachUnique(f, bonds); - return U_Bonds; -} - - -template -void updateBond(MPS& psi, size_t i, const e::Tensor& U, unsigned chiMax, T eps) { - // std::println("BOND {}", i); - auto theta = psi.getTheta2(i); // vL i j vR - // apply U: i j i* j* - // contraction gives vL vR i j -> shuffle to vL i j vR - // std::cout << theta.dimension(2) << " - " << U.dimension(2) << " / " << theta.dimensions() << std::endl; - e::Tensor a = theta.contract(U, e::array({iP{1, 2}, iP{2, 3}})); - e::Tensor b = a.shuffle(e::array{0, 2, 3, 1}); - // e::Tensor b = a.shuffle(e::array{2, 0, 1, 3}); - // printTensor(b, "Utheta"); - // SVD and truncation - auto [Ai, Sj, Bj] = split_truncate_theta(b, chiMax, eps); - // printTensor(Ai, "Ai"); - // printTensor(Sj, "Sj"); - // printTensor(Bj, "Bj"); - // std::cout << "a: " << a.dimensions() << "\nA: " << Ai.dimensions() << "\nS: " << Sj.dimensions() << "\nB: " << Bj.dimensions() << std::endl; - // renormalize to Bi - e::Tensor SiInv = 1. / psi.S(i); - e::Tensor SiDiv = diagonal(SiInv); - e::Tensor Gi = SiDiv.contract(Ai, e::array{iP{1, 0}}); // vL [vL], [vL] i vC - // printTensor(SiDiv, "SiDiv"); - // printTensor(Gi, "Gi"); - - // e::Tensor SS = Sj.contract(Sj.conjugate(), e::array{iP{0,0}}); - // std::println("Bond {:03}, Sj contracted with its conjugate = {}", i, SS(0)); - - // psi.B(i) = Gi.contract(diagonal(Sj), e::array{iP{2, 0}}); // vL i [vC], [vC] vC - // psi.B(i+1) = std::move(Bj); - // psi.S(i+1) = std::move(Sj); - MPS::BTensor_t Bi = Gi.contract(diagonal(Sj), e::array{iP{2, 0}}); // vL i [vC], [vC] vC - // e::Tensor equal = (psi.B(i) == Bi).all(); - // printTensor(psi.B(i), "old Bi"); - // printTensor(Bi, "new Bi"); - psi.B(i) = Bi; - psi.B(i+1) = Bj; - psi.S(i+1) = Sj; -} - - -template -void runTEBD_brickwall(MPS& psi, const Bonds& bonds, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc) { - for (unsigned n = 0; n < nSteps; n++) { - // std::print("i = {:03}\r", n); - // std::flush(std::cout); - for (unsigned row = 0; row < 2; row++) { - for (unsigned i = row; i < bonds.size(); i+=2) { - // std::cout << "i = " << i << std::endl; - updateBond(psi, i, *(bonds.at(i)), chiMax, eps); - } - } - postUpdateFunc(psi, n); +namespace tebd { + template + e::Tensor getTimeEvolutionMatrix(const e::MatrixX& bond, T dt, unsigned localDim) { + e::MatrixX U = (-bond * dt).exp(); + // std::cout << "U = \n" << U << "\n\n"; + return bondMatrix2Tensor(U, localDim); } -} -// TODO -/** - * @brief Run TEBD using a second-order decomposition - * @todo I think this is currently only correct for even system sizes => check decomposition - */ -template -void runTEBD_secondOrder(MPS& psi, const Bonds& bonds, const Bonds& bondsHalfDT, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc) { - for (unsigned n = 0; n < nSteps; n++) { - // std::print("i = {:03}\r", n); - // std::flush(std::cout); - for (unsigned row = 0; row < 2; row++) { - const Bonds* useBonds = &bonds; - if (row == 0 && (n != 0 || n != nSteps-1)) { - useBonds = &bondsHalfDT; - } - for (unsigned i = row; i < bonds.size(); i+=2) { - // std::cout << "i = " << i << std::endl; - updateBond(psi, i, *(useBonds->at(i)), chiMax, eps); - } - } - postUpdateFunc(psi, n); + template + Bonds getTimeEvolutionMatrices(const std::vector>>& bonds, T dt, unsigned int localDim) { + std::function(const e::MatrixX&)> f = std::bind(getTimeEvolutionMatrix, std::placeholders::_1, dt, localDim); + std::vector>> U_Bonds = forEachUnique(f, bonds); + return U_Bonds; } -} -// Instantiations + template + void updateBond(MPS& psi, size_t i, const e::Tensor& U, unsigned chiMax, T eps) { + // std::println("BOND {}", i); + auto theta = psi.getTheta2(i); // vL i j vR + // apply U: i j i* j* + // contraction gives vL vR i j -> shuffle to vL i j vR + // std::cout << theta.dimension(2) << " - " << U.dimension(2) << " / " << theta.dimensions() << std::endl; + e::Tensor a = theta.contract(U, e::array({iP{1, 2}, iP{2, 3}})); + e::Tensor b = a.shuffle(e::array{0, 2, 3, 1}); + // e::Tensor b = a.shuffle(e::array{2, 0, 1, 3}); + // printTensor(b, "Utheta"); + // SVD and truncation + auto [Ai, Sj, Bj] = split_truncate_theta(b, chiMax, eps); + // printTensor(Ai, "Ai"); + // printTensor(Sj, "Sj"); + // printTensor(Bj, "Bj"); + // std::cout << "a: " << a.dimensions() << "\nA: " << Ai.dimensions() << "\nS: " << Sj.dimensions() << "\nB: " << Bj.dimensions() << std::endl; + // renormalize to Bi + e::Tensor SiInv = 1. / psi.S(i); + e::Tensor SiDiv = diagonal(SiInv); + e::Tensor Gi = SiDiv.contract(Ai, e::array{iP{1, 0}}); // vL [vL], [vL] i vC + // printTensor(SiDiv, "SiDiv"); + // printTensor(Gi, "Gi"); + + // e::Tensor SS = Sj.contract(Sj.conjugate(), e::array{iP{0,0}}); + // std::println("Bond {:03}, Sj contracted with its conjugate = {}", i, SS(0)); + + // psi.B(i) = Gi.contract(diagonal(Sj), e::array{iP{2, 0}}); // vL i [vC], [vC] vC + // psi.B(i+1) = std::move(Bj); + // psi.S(i+1) = std::move(Sj); + MPS::BTensor_t Bi = Gi.contract(diagonal(Sj), e::array{iP{2, 0}}); // vL i [vC], [vC] vC + // e::Tensor equal = (psi.B(i) == Bi).all(); + // printTensor(psi.B(i), "old Bi"); + // printTensor(Bi, "new Bi"); + psi.B(i) = Bi; + psi.B(i+1) = Bj; + psi.S(i+1) = Sj; + } + + + template + void runBrickwall(MPS& psi, const Bonds& bonds, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc) { + for (unsigned n = 0; n < nSteps; n++) { + // std::print("i = {:03}\r", n); + // std::flush(std::cout); + for (unsigned row = 0; row < 2; row++) { + for (unsigned i = row; i < bonds.size(); i+=2) { + // std::cout << "i = " << i << std::endl; + // std::println("Updating bond {:02}", i); + updateBond(psi, i, *(bonds.at(i)), chiMax, eps); + } + } + if (postUpdateFunc(psi, n)) break; + } + } + + // TODO + /** + * @brief Run TEBD using a second-order decomposition + * @todo I think this is currently only correct for even system sizes => check decomposition + */ + template + void runSecondOrder(MPS& psi, const Bonds& bonds, const Bonds& bondsHalfDT, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc) { + for (unsigned n = 0; n < nSteps; n++) { + // std::print("i = {:03}\r", n); + // std::flush(std::cout); + for (unsigned row = 0; row < 2; row++) { + const Bonds* useBonds = &bonds; + if (row == 0 && (n != 0 || n != nSteps-1)) { + useBonds = &bondsHalfDT; + } + for (unsigned i = row; i < bonds.size(); i+=2) { + // std::cout << "i = " << i << std::endl; + updateBond(psi, i, *(useBonds->at(i)), chiMax, eps); + } + } + if (postUpdateFunc(psi, n)) break; + } + } + + + // Instantiations + template e::Tensor getTimeEvolutionMatrix(const e::MatrixX& bond, dtype dt, unsigned localDim); + template Bonds getTimeEvolutionMatrices(const std::vector>>& bonds, dtype dt, unsigned int localDim); + template void updateBond(MPS& psi, size_t i, const e::Tensor& U, unsigned chiMax, dtype eps); + template void runBrickwall(MPS& psi, const Bonds& bonds, unsigned chiMax, dtype eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc); + template void runSecondOrder(MPS& psi, const Bonds& bonds, const Bonds&, unsigned chiMax, dtype eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc); + +}; // namespace tebd template std::vector>> forEachUnique(std::function (const e::MatrixX &)> f, const std::vector>>&); -template e::Tensor getTimeEvolutionMatrix(const e::MatrixX& bond, dtype dt, unsigned localDim); -template Bonds getTimeEvolutionMatrices(const std::vector>>& bonds, dtype dt, unsigned int localDim); -template void updateBond(MPS& psi, size_t i, const e::Tensor& U, unsigned chiMax, dtype eps); -template void postUpdateNoop(const MPS&, unsigned); -template void runTEBD_brickwall(MPS& psi, const Bonds& bonds, unsigned chiMax, dtype eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc); -template void runTEBD_secondOrder(MPS& psi, const Bonds& bonds, const Bonds&, unsigned chiMax, dtype eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc); - diff --git a/src/algorithms/tebd.hpp b/src/algorithms/tebd.hpp index 622671a..500918d 100644 --- a/src/algorithms/tebd.hpp +++ b/src/algorithms/tebd.hpp @@ -1,9 +1,8 @@ #pragma once -#include "model.hpp" +#include "algorithm.hpp" #include "mps.hpp" #include #include -#include /** * @brief Apply the function f on every element if inRange, but reuse a cached result when two pointers point to the same object @@ -15,40 +14,40 @@ template using Bonds = std::vector>>; /** - * @brief Apply time evolution on a bond - * @details - * \f[U = \exp{ B \, \del t} \f] + * @brief Implementation of the Time Evolved block Decimation (TEBD) algorithm */ -template -e::Tensor getTimeEvolutionMatrix(const e::MatrixX& bond, T dt, unsigned localDim); +namespace tebd { + /** + * @brief Apply time evolution on a bond + * @details + * \f[U = \exp{ B \, \partial t} \f] + */ + template + e::Tensor getTimeEvolutionMatrix(const e::MatrixX& bond, T dt, unsigned localDim); -template -Bonds getTimeEvolutionMatrices(const std::vector>>& bonds, T dt, unsigned int localDim); + template + Bonds getTimeEvolutionMatrices(const std::vector>>& bonds, T dt, unsigned int localDim); -/** - * @brief Apply U on bond i, perform SVD and truncate, update MPS - */ -template -void updateBond(MPS& psi, size_t i, const e::Tensor& U, unsigned chiMax, T eps); + /** + * @brief Apply U on bond i, perform SVD and truncate, update MPS + */ + template + void updateBond(MPS& psi, size_t i, const e::Tensor& U, unsigned chiMax, T eps); -template -using postUpdateFunction_t = std::function&, unsigned)>; -template -void postUpdateNoop(const MPS&, unsigned) {}; + /** + * @brief Brickwall TEBS scheme + * @details + * Apply U on every second bond starting from the first, + * then apply U on every second bond starting from the second + * + * First order scheme + * + * @param postUpdateFunc function that can be used to track expectation values after each update + */ + template + void runBrickwall(MPS& psi, const Bonds& bonds, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc); -/** - * @brief Brickwall TEBS scheme - * @details - * Apply U on every second bond starting from the first, - * then apply U on every second bond starting from the second - * - * First order scheme - * - * @param postUpdateFunc function that can be used to track expectation values after each update - */ -template -void runTEBD_brickwall(MPS& psi, const Bonds& bonds, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc); - -template -void runTEBD_secondOrder(MPS& psi, const Bonds& bonds, const Bonds& bondsHalfDT, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc); + template + void runSecondOrder(MPS& psi, const Bonds& bonds, const Bonds& bondsHalfDT, unsigned chiMax, T eps, unsigned nSteps, postUpdateFunction_t postUpdateFunc); +}; // namespace tebd diff --git a/src/benchmark.cpp b/src/benchmark.cpp index 325ca87..a08ac9b 100644 --- a/src/benchmark.cpp +++ b/src/benchmark.cpp @@ -41,12 +41,13 @@ Results _runTEBD(unsigned L, double J, double g, unsigned nSteps, dtype dt, dtyp result.M.push_back(psi.getSiteExpectationValueSum(sigma_z).real()); result.E.push_back(model.energy(psi).real()); result.S.push_back(psi.getEntanglementEntropy(middleBond).real()); + return false; }); // initial values now F(psi, 0); - const Bonds bonds = getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim); - runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, F); + const Bonds bonds = tebd::getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim); + tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, F); return result; }; diff --git a/src/main.cpp b/src/main.cpp index 146952c..a5556dc 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -7,11 +7,29 @@ #include "tebd.hpp" #include "util.hpp" #include "type.hpp" +#include "dmrg.hpp" +#include "wasm/run.hpp" #include "benchmark.hpp" using namespace std::complex_literals; +void testWasmCode() { + unsigned L = 6; + unsigned nMax = 3; + auto model = std::make_shared>(L, dtype(1.), dtype(0.05), dtype(0.5), nMax); + auto psi = initMPS_nParticles(L, nMax, 1); + + double eps = 1e-10; + unsigned chiMax = 20; + double dtReal = 0.00; + double dtImag = 0.05; + unsigned nSteps = 300; + ems::val callback; + BoseHubbardResults result(callback, 0); + runTEBD(psi, model, nSteps, dtReal, dtImag, eps, chiMax, TEBD_TYPE::BRICKWALL, result); +} + void testSVD() { e::Tensor t4(2, 3, 4, 5); auto data = t4.data(); @@ -30,19 +48,20 @@ void testSVD() { void testBoseHubbardTEBD() { - unsigned L = 8; - unsigned nMax = 5; - auto model = BoseHubbardModel1D(L, dtype(1.), dtype(2.50), dtype(0.5), nMax); + unsigned L = 6; + unsigned nMax = 3; + auto model = BoseHubbardModel1D(L, dtype(1.), dtype(0.05), dtype(0.5), nMax); dtype eps = 1e-10; unsigned chiMax = 20; - dtype dt = 0.05; + dtype dt = 0.05i; + unsigned nSteps = 300; std::println("dt={}", dt); // std::cout << "leftT:\n" << *(model.getH_Bonds().front()) << std::endl; // std::cout << "leftT2:\n" << bondMatrix2Tensor(*(model.getH_BondsMatrices().front()), model.localDim) << std::endl; std::cout << "middleT:\n" << *(model.getH_Bonds().at(1)) << std::endl; // - auto bonds = getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim); + auto bonds = tebd::getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim); std::cout << "N:\n" << model.getLadderOperators().getN() << std::endl; std::cout << "a:\n" << model.getLadderOperators().getA() << std::endl; std::cout << "aDagger:\n" << model.getLadderOperators().getADagger() << std::endl; @@ -50,36 +69,36 @@ void testBoseHubbardTEBD() { // std::cout << "H_Bond:\n" << *(bonds.front()) << std::endl; // MPS psi = initMPS_spinUp(L); - MPS psi = initMPS_nParticles(L, nMax, 2); + MPS psi = initMPS_nParticles(L, nMax, 1); - unsigned nSteps = 300; // postUpdateFunction_t fUpdate = postUpdateNoop; e::TensorMap> particleNumberOperator(model.getLadderOperators().getN().data(), model.localDim, model.localDim); postUpdateFunction_t fUpdate = [&model, &particleNumberOperator](const MPS& psi, unsigned i) { std::println("{:03}: N={:.2f} {}", i, psi.getSiteExpectationValueSum(particleNumberOperator).real() / model.L, psi.getChi()); + return false; }; // auto fUpdate = postUpdateFunction_t([](const MPS& psi, unsigned i){ // std::println("i={:03}: ={}", i, psi.collapse(0)); // }); - std::println("E = {}", model.energy(psi).real()); + // std::println("E = {}", model.energy(psi).real()); // std::println("m = {}", psi.getSiteExpectationValues(sigma_z)); - std::println("chiS: {:b}", psi.getChi()); - psi.print(); - std::println("={}", psi.collapse(0)); - - runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, fUpdate); - std::println("E = {}", model.energy(psi).real()); - std::println(" = {}", psi.getSiteExpectationValueSum(particleNumberOperator).real() / model.L); - // std::println("m = {}", psi.getSiteExpectationValues(sigma_z)); - std::println("chiS: {:b}", psi.getChi()); + // std::println("chiS: {:b}", psi.getChi()); // psi.print(); - std::println("={}", psi.collapse(0)); + // std::println("={}", psi.collapse(0)); + + tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, fUpdate); + // std::println("E = {}", model.energy(psi).real()); + // std::println(" = {}", psi.getSiteExpectationValueSum(particleNumberOperator).real() / model.L); + // std::println("m = {}", psi.getSiteExpectationValues(sigma_z)); + // std::println("chiS: {:b}", psi.getChi()); + // psi.print(); + // std::println("={}", psi.collapse(0)); } void testTFI_TEBD() { - unsigned L = 8; + unsigned L = 40; std::println("Make model"); auto model = TransverseFieldIsingModel1D(L, dtype(1.), dtype(1.5)); // auto& mpo = *(model.getH_MPO().front()); @@ -109,21 +128,21 @@ void testTFI_TEBD() { e::Tensor sigma_z(2, 2); sigma_z.setValues({{1, 0}, {0, -1}}); - dtype eps = 1e-10; - unsigned chiMax = 20; - dtype dt = 0.05i; + dtype eps = 1e-20; + unsigned chiMax = 30; + dtype dt = 0.05; std::println("dt={}", dt); std::cout << "leftT:\n" << *(model.getH_Bonds().front()) << std::endl; std::cout << "leftT2:\n" << bondMatrix2Tensor(*(model.getH_BondsMatrices().front()), model.localDim) << std::endl; std::cout << "middleT:\n" << *(model.getH_Bonds().at(1)) << std::endl; - auto bonds = getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim); + auto bonds = tebd::getTimeEvolutionMatrices(model.getH_BondsMatrices(), dt, model.localDim); std::cout << "U_Bond:\n" << *(bonds.front()) << std::endl; // std::cout << "H_Bond:\n" << *(bonds.front()) << std::endl; // MPS psi = initMPS_spinUp(L); MPS psi = initMPS_spinRight(L); - unsigned nSteps = 100; + unsigned nSteps = 200; postUpdateFunction_t fUpdate = postUpdateNoop; // auto fUpdate = postUpdateFunction_t([](const MPS& psi, unsigned i){ // std::println("i={:03}: ={}", i, psi.collapse(0)); @@ -134,7 +153,7 @@ void testTFI_TEBD() { psi.print(); std::println("={}", psi.collapse(0)); - runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, fUpdate); + tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, fUpdate); // updateBond(psi, 0, *(bonds.at(0)), chiMax, eps); // updateBond(psi, 2, *(bonds.at(2)), chiMax, eps); // updateBond(psi, 4, *(bonds.at(4)), chiMax, eps); @@ -148,7 +167,7 @@ void testTFI_TEBD() { std::println("E = {}", model.energy(psi).real()); std::println("m = {}", psi.getSiteExpectationValues(sigma_z)); std::println("chiS: {:b}", psi.getChi()); - psi.print(); + // psi.print(); std::println("={}", psi.collapse(0)); // for (unsigned i = 0; i < L; i++) { @@ -157,10 +176,59 @@ void testTFI_TEBD() { } +void testTFI_DMRG() { + unsigned L = 8; + std::println("Make model"); + auto model = TransverseFieldIsingModel1D(L, dtype(1.), dtype(1)); + e::Tensor sigma_z(2, 2); + sigma_z.setValues({{1, 0}, {0, -1}}); + + std::cout << *(model.getH_BondsMatrices().at(0)) << std::endl; + printTensor(*model.getH_Bonds().at(0), "H0"); + std::cout << *(model.getH_BondsMatrices().at(1)) << std::endl; + printTensor(*model.getH_Bonds().at(1), "H1"); + std::cout << *(model.getH_BondsMatrices().at(6)) << std::endl; + printTensor(*model.getH_Bonds().at(6), "H6"); + + + dtype eps = 1e-20; + unsigned chiMax = 10; + MPS psi = initMPS_spinUp(L); + for (unsigned i = 0; i < 7; i++) { + std::println("{}", psi.getBondExpectationValue(*(model.getH_Bonds()[i]), i)); + } + // MPS psi = initMPS_spinRight(L); + // MPS psi = initMPS_spinDown(L); + dtype E = model.energy(psi); + std::println("E{} = {}", -1, E.real()); + std::println("M{} = {}", -1, psi.getSiteExpectationValueSum(sigma_z)); + + printTensor(*(model.getH_MPO().at(0)), "H_mpo[0]"); + printTensor(*(model.getH_MPO().at(1)), "H_mpo[1]"); + + auto engine = dmrg::DMRGEngine(psi, model.getH_MPO(), chiMax, eps); + std::println("Initialized"); + E = model.energy(psi); + std::println("E{} = {}", 0, E); + // engine.update_bond(0); + for (unsigned i = 0; i < 4; i++) { + engine.sweep(psi); + std::println("\nAfter sweep #{}", i+1); + E = model.energy(psi); + std::println("E = {:.2f}", E); + std::println("M = {:.2f}", psi.getSiteExpectationValueSum(sigma_z)); + } + + // engine.psi.print(); + +} + int main() { + testTFI_DMRG(); // testSVD(); return 0; - testBoseHubbardTEBD(); + // testBoseHubbardTEBD(); // testTFI_TEBD(); + // testWasmCode(); // return runBenchmark_oneAtATime(); // return runBenchmark_all(); diff --git a/src/model.cpp b/src/model.cpp index 9c71c19..d0cd674 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -12,6 +12,7 @@ T Model1D::energy(const MPS& psi) const { return E; } + template e::Tensor constructMPO(const std::initializer_list>>& matrices) { size_t rows = matrices.size(); @@ -42,7 +43,6 @@ e::Tensor constructMPO(const std::initializer_list e::MatrixX MPO_Tensor2Matrix(const e::Tensor& mpo) { size_t dim = mpo.dimension(0); diff --git a/src/model.hpp b/src/model.hpp index c4c78f4..db1b207 100644 --- a/src/model.hpp +++ b/src/model.hpp @@ -17,12 +17,13 @@ using namespace std::complex_literals; template e::Tensor constructMPO(const std::initializer_list>>& matrices); + template e::MatrixX MPO_Tensor2Matrix(const e::Tensor& mpo); /** - * @brief Reshape a N*localDim x N*localDim matrix into a N x localDim x N x localDim tensor + * @brief Reshape a \f$N \, d \times N \, d\f$ matrix into a \f$N \times d \times N \times d\f] tensor */ template e::Tensor bondMatrix2Tensor(const e::MatrixX& bondMatrix, unsigned localDim); @@ -39,7 +40,7 @@ class Model1D { virtual ~Model1D() {}; inline const std::vector>>& getH_BondsMatrices() const { return H_BondsMatrices; } inline const std::vector>>& getH_Bonds() const { return H_Bonds; } - inline const std::vector>>& getH_MPO() const { return H_MPO; } + inline std::vector>>& getH_MPO() { return H_MPO; } const size_t L; const unsigned localDim; T energy(const MPS& psi) const; diff --git a/src/models/bose_hubbard.cpp b/src/models/bose_hubbard.cpp index 6c803b8..4c7caae 100644 --- a/src/models/bose_hubbard.cpp +++ b/src/models/bose_hubbard.cpp @@ -1,6 +1,7 @@ #include "bose_hubbard.hpp" #include + template LadderOperators::LadderOperators(unsigned dim) : identity(dim, dim), aDagger(dim, dim), a(dim, dim), n(dim, dim) { identity.setZero(); @@ -22,34 +23,37 @@ LadderOperators::LadderOperators(unsigned dim) : identity(dim, dim), aDagger( template void BoseHubbardModel1D::initH_Bonds() { - T tHalf = this->t * static_cast(0.5); - T tFull = this->t; + T UHalf = this->U * static_cast(0.5); + T UQuarter = this->U * static_cast(0.25); + T muHalf = this->mu * static_cast(0.5); /// @todo TODO figure out which is right. Doesnt t need to be halfed at the edges and not the middle? e::MatrixX N_NminusI = this->lo.getN() * (this->lo.getN() - this->lo.getIdentity()); - e::MatrixX singleParticlePart( - ( - e::kroneckerProduct(this->lo.getIdentity(), N_NminusI) - - e::kroneckerProduct(N_NminusI, this->lo.getIdentity()) - ) * (this->U/static_cast(2)) - - e::kroneckerProduct(this->lo.getN(), this->lo.getIdentity()) * (this->mu/static_cast(2)) - - e::kroneckerProduct(this->lo.getIdentity(), this->lo.getN()) * (this->mu/static_cast(2)) + + e::MatrixX twoParticlePart = + - this->t * e::kroneckerProduct(this->lo.getADagger(), this->lo.getA()) + - this->t * e::kroneckerProduct(this->lo.getA(), this->lo.getADagger() ); - // TODO this is wrong, needs A and Adagger; check which t or t/half or identities are required auto left = std::make_shared>( - singleParticlePart - - tFull * e::kroneckerProduct(this->lo.getADagger(), this->lo.getA()) - - tFull * e::kroneckerProduct(this->lo.getA(), this->lo.getADagger()) + twoParticlePart + + e::kroneckerProduct(N_NminusI, this->lo.getIdentity()) * UHalf + + e::kroneckerProduct(this->lo.getIdentity(), N_NminusI) * UQuarter + - e::kroneckerProduct(this->lo.getN(), this->lo.getIdentity()) * (this->mu) + - e::kroneckerProduct(this->lo.getIdentity(), this->lo.getN()) * (muHalf) ); auto right = std::make_shared>( - singleParticlePart - - tFull * e::kroneckerProduct(this->lo.getADagger(), this->lo.getA()) - - tFull * e::kroneckerProduct(this->lo.getA(), this->lo.getADagger()) + twoParticlePart + + e::kroneckerProduct(N_NminusI, this->lo.getIdentity()) * UQuarter + + e::kroneckerProduct(this->lo.getIdentity(), N_NminusI) * UHalf + - e::kroneckerProduct(this->lo.getN(), this->lo.getIdentity()) * (muHalf) + - e::kroneckerProduct(this->lo.getIdentity(), this->lo.getN()) * (this->mu) ); auto middle = std::make_shared>( - singleParticlePart - - tFull * e::kroneckerProduct(this->lo.getADagger(), this->lo.getA()) - - tFull * e::kroneckerProduct(this->lo.getA(), this->lo.getADagger()) + twoParticlePart + + e::kroneckerProduct(N_NminusI, this->lo.getIdentity()) * UQuarter + + e::kroneckerProduct(this->lo.getIdentity(), N_NminusI) * UQuarter + - e::kroneckerProduct(this->lo.getN(), this->lo.getIdentity()) * (muHalf) + - e::kroneckerProduct(this->lo.getIdentity(), this->lo.getN()) * (muHalf) ); // The matrices are N*N x N*N -> reshape to N x N x N x N tensors unsigned N = this->localDim; @@ -90,7 +94,7 @@ void BoseHubbardModel1D::initH_MPO() { zero.setZero(); auto W = std::make_shared>(std::move(constructMPO({ - {this->lo.getIdentity(), this->lo.getADagger(), this->lo.getA(), (this->U / 2. * (this->lo.getN() * this->lo.getN() - this->lo.getN()) - this->mu * this->lo.getN())}, + {this->lo.getIdentity(), this->lo.getADagger(), this->lo.getA(), this->U / 2. * (this->lo.getN() * this->lo.getN() - this->lo.getN()) - this->mu * this->lo.getN()}, {zero, zero, zero, (-this->t * this->lo.getA())}, {zero, zero, zero, (-this->t * this->lo.getADagger())}, {zero, zero, zero, this->lo.getIdentity()} diff --git a/src/models/bose_hubbard.hpp b/src/models/bose_hubbard.hpp index 28a4b3d..de6df2c 100644 --- a/src/models/bose_hubbard.hpp +++ b/src/models/bose_hubbard.hpp @@ -1,6 +1,7 @@ #pragma once #include "model.hpp" + /** * @brief Helper class constructing identity, particle number and both ladder operators in one loop. */ @@ -19,14 +20,44 @@ class LadderOperators { e::MatrixX n; }; + /** * @brief Bose-Hubbard Model in 1 dimension + * Latex test \f$ 1 +2 \otimes 3\f$ * @details + * The Hamiltonian: + * \f[ + * H = + * -t \sum_{i}\left( + * \hat{a}^\dagger_i \hat{a}_{i+1} + * + \hat{a}_i \hat{a}^\dagger_{i+1} + * \right) + * + \frac{U}{2} \sum_i \hat{n}_i \left(\hat{n}_i-\mathbbm{1}\right) + * - \mu \sum_i \hat{n}_i + * \f] + * * The middle bonds are: - * TODO \f[ B_i = J \, \sigma^x_i \otimes \sigma^x_{i+1} + \frac{g}{2} \,\sigma^z_i \otimes I + \frac{g}{2} \, I \otimes \sigma^z_{i+1} \f] + * \f[ + * -t \left( + * \hat{a}^\dagger_i \otimes \hat{a}_{i+1} + * + \hat{a}_i \otimes \hat{a}^\dagger_{i+1} + * \right) + * + \frac{U}{4} \Big( + * \hat{n}_{i} \left(\hat{n}_{i}-1_{i}\right) \otimes \mathbbm{1}_{i+1} + * + \mathbbm{1}_{i} \otimes \hat{n}_{i+1} \left(\hat{n}_{i+1}-1_{i+1} \right) + * \Big) + * - \frac{\mu}{2} \left( + * \hat{n}_i \otimes \mathbbm{1}_{i+1} + * \mathbbm{1}_{i} \otimes \hat{n}_{i+1} + * \right) + * \f] * * The MPO is: - * \f[ W = \begin{pmatrix} I & \hat{a}^\dagger & \hat{a} & \frac{U}{2}\hat{n}(\hat{n} - I) - \mu\hat{n} \\ 0 & 0 & 0 & -t \hat{a} \\ 0 & 0 & 0 & -t \hat{a}^\dagger\\ 0 & 0 & 0 & I \end{pmatrix} \f] + * \f[ W = \begin{pmatrix} + * \mathbbm{1} & \hat{a}^\dagger & \hat{a} & \frac{U}{2}\hat{n}(\hat{n} - \mathbbm{1}) - \mu\hat{n} \\ + * 0 & 0 & 0 & -t \hat{a} \\ 0 & 0 & 0 & -t \hat{a}^\dagger\\ + * 0 & 0 & 0 & \mathbbm{1} \end{pmatrix} + * \f] */ template class BoseHubbardModel1D : public Model1D { @@ -43,7 +74,6 @@ class BoseHubbardModel1D : public Model1D { virtual std::string format() const override { return std::format("BoseHubbardModel1D(L={},U={},t={},mu={})", this->L, this->U, this->t, this->mu); } - protected: void initH_Bonds() override; void initH_MPO() override; @@ -53,5 +83,3 @@ class BoseHubbardModel1D : public Model1D { T mu; LadderOperators lo; }; - - diff --git a/src/models/transverse_field_ising.cpp b/src/models/transverse_field_ising.cpp index e677cd4..e5b2887 100644 --- a/src/models/transverse_field_ising.cpp +++ b/src/models/transverse_field_ising.cpp @@ -8,22 +8,22 @@ PauliMatrices TransverseFieldIsingModel1D::pm; template void TransverseFieldIsingModel1D::initH_Bonds() { - T gHalf = this->g * static_cast(0.5); + T BHalf = this->B * static_cast(0.5); // TODO: improve code or directly initialize the tensors auto left = std::make_shared>( - e::kroneckerProduct(pm.getSigmaX(), pm.getSigmaX()) * (-J) - - this->g * e::kroneckerProduct(pm.getSigmaZ(), pm.getIdentity()) - - gHalf * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaZ()) + e::kroneckerProduct(pm.getSigmaZ(), pm.getSigmaZ()) * (-J) + - this->B * e::kroneckerProduct(pm.getSigmaX(), pm.getIdentity()) + - BHalf * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaX()) ); auto right = std::make_shared>( - e::kroneckerProduct(pm.getSigmaX(), pm.getSigmaX()) * (-J) - - gHalf * e::kroneckerProduct(pm.getSigmaZ(), pm.getIdentity()) - - this->g * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaZ()) + e::kroneckerProduct(pm.getSigmaZ(), pm.getSigmaZ()) * (-J) + - BHalf * e::kroneckerProduct(pm.getSigmaX(), pm.getIdentity()) + - this->B * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaX()) ); auto middle = std::make_shared>( - e::kroneckerProduct(pm.getSigmaX(), pm.getSigmaX()) * (-J) - - gHalf * e::kroneckerProduct(pm.getSigmaZ(), pm.getIdentity()) - - gHalf * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaZ()) + e::kroneckerProduct(pm.getSigmaZ(), pm.getSigmaZ()) * (-J) + - BHalf * e::kroneckerProduct(pm.getSigmaX(), pm.getIdentity()) + - BHalf * e::kroneckerProduct(pm.getIdentity(), pm.getSigmaX()) ); auto leftT = std::make_shared>(2, 2, 2, 2); auto middleT = std::make_shared>(2, 2, 2, 2); @@ -63,7 +63,7 @@ void TransverseFieldIsingModel1D::initH_MPO() { e::Matrix2 zero(2, 2); zero.setZero(); auto W = std::make_shared>(std::move(constructMPO({ - {this->pm.getIdentity(), this->pm.getSigmaZ(), (this->g * this->pm.getSigmaX())}, + {this->pm.getIdentity(), this->pm.getSigmaZ(), (-this->B) * this->pm.getSigmaX()}, {zero, zero, (-this->J * this->pm.getSigmaZ())}, {zero, zero, this->pm.getIdentity()} }))); diff --git a/src/models/transverse_field_ising.hpp b/src/models/transverse_field_ising.hpp index b7c414b..64d9b39 100644 --- a/src/models/transverse_field_ising.hpp +++ b/src/models/transverse_field_ising.hpp @@ -4,25 +4,33 @@ /** * @brief Transverse Field Ising Model in 1 dimension * @details + * The nearest neighbor interaction is in z direction, the magnetic field in x direction. + * The Hamiltonian is: + * \f[ H = -J \sum_{\langle i,j\rangle} \sigma_{z,i} \sigma_{z,j} - B \sum_i \sigma_{x,j} \f] + * * The middle bonds are: - * \f[ B_i = J \, \sigma^x_i \otimes \sigma^x_{i+1} + \frac{g}{2} \,\sigma^z_i \otimes I + \frac{g}{2} \, I \otimes \sigma^z_{i+1} \f] + * \f[ B_i = J \, \sigma^z_i \otimes \sigma^z_{i+1} + \frac{B}{2} \,\sigma^x_i \otimes I + \frac{B}{2} \, I \otimes \sigma^x_{i+1} \f] * * The MPO is: - * \f[ W = \begin{pmatrix} I & \sigma_z & g \, \sigma_x \\ 0 & 0 & J\,\sigma_z \\ 0 & 0 & I \end{pmatrix} \f] + * \f[ W = \begin{pmatrix} + * \mathbbm{1} & \sigma_z & -B \, \sigma_x \\ + * 0 & 0 & -J\,\sigma_z \\ + * 0 & 0 & \mathbbm{1} + * \end{pmatrix} + * \f] */ template class TransverseFieldIsingModel1D : public Model1D { public: - TransverseFieldIsingModel1D(size_t L, T J, T g) : Model1D(L, 2), J(J), g(g) { - // std::println("TFIModel: L={}, J={}, g={}", L, J, g); + TransverseFieldIsingModel1D(size_t L, T J, T B) : Model1D(L, 2), J(J), B(B) { this->initH_Bonds(); this->initH_MPO(); } inline const T& getJ() const { return J; } - inline const T& getg() const { return g; } + inline const T& getB() const { return B; } virtual std::string format() const override { - return std::format("TransverseFieldIsingModel1D(L={},J={},g={})", this->L, this->J, this->g); + return std::format("TransverseFieldIsingModel1D(L={},J={},B={})", this->L, this->J, this->B); } protected: void initH_Bonds() override; @@ -30,7 +38,7 @@ class TransverseFieldIsingModel1D : public Model1D { void initH_MPO() override; private: T J; - T g; + T B; static PauliMatrices pm; }; diff --git a/src/mpl_colors.py b/src/mpl_colors.py deleted file mode 100644 index 3e2ff2e..0000000 --- a/src/mpl_colors.py +++ /dev/null @@ -1,154 +0,0 @@ -import matplotlib as mpl -import matplotlib.pyplot as plt -from cycler import cycler - -# Palette -P = { - "bg0": "#282828", - "bg0-hard": "#1d2021", - "bg0-soft": "#32302f", - "bg1": "#3c3836", - "bg2": "#504945", - "bg3": "#665c54", - "bg4": "#7c6f64", - "fg0": "#fbf1c7", - "fg0-hard": "#f9f5d7", - "fg0-soft": "#f2e5bc", - "fg1": "#ebdbb2", - "fg2": "#d5c4a1", - "fg3": "#bdae93", - "fg4": "#a89984", - "dark-red": "#cc241d", - "dark-green": "#98971a", - "dark-yellow": "#d79921", - "dark-blue": "#458588", - "dark-purple": "#b16286", - "dark-aqua": "#689d6a", - "dark-orange": "#d65d0e", - "dark-gray": "#928374", - "light-red": "#fb4934", - "light-green": "#b8bb26", - "light-yellow": "#fabd2f", - "light-blue": "#83a598", - "light-purple": "#d3869b", - "light-aqua": "#8ec07c", - "light-orange": "#f38019", - "light-gray": "#a89984", - "alt-red": "#9d0006", - "alt-green": "#79740e", - "alt-yellow": "#b57614", - "alt-blue": "#076678", - "alt-purple": "#8f3f71", - "alt-aqua": "#427b58", - "alt-orange": "#af3a03", - "alt-gray": "#7c6f64", -} - -COLORS_LIGHT = [ - P["light-blue"], - P["light-orange"], - P["light-green"], - P["light-red"], - P["light-purple"], - P["light-yellow"], - P["light-aqua"], - P["light-gray"]] - -COLORS_DARK = [ - P["dark-blue"], - P["dark-orange"], - P["dark-green"], - P["dark-red"], - P["dark-purple"], - P["dark-yellow"], - P["dark-aqua"], - P["dark-gray"]] - -COLORS_ALT = [ - P["alt-blue"], - P["alt-orange"], - P["alt-green"], - P["alt-red"], - P["alt-purple"], - P["alt-yellow"], - P["alt-aqua"], - P["alt-gray"]] - -COLORS = COLORS_LIGHT - -def set_colorscheme(variant="dark"): - global COLORS - if variant == "dark": - FIG_BG_COLOR = P["bg0"] - PLT_FG_COLOR = P["fg0"] - PLT_BG_COLOR = P["bg0"] - PLT_GRID_COLOR = P["bg2"] - LEGEND_FG_COLOR = PLT_FG_COLOR - LEGEND_BG_COLOR = P["bg1"] - LEGEND_BORDER_COLOR = P["bg2"] - COLORS = COLORS_LIGHT - else: - FIG_BG_COLOR = P["fg0"] - PLT_FG_COLOR = P["bg0"] - PLT_BG_COLOR = P["fg0"] - PLT_GRID_COLOR = P["fg2"] - LEGEND_FG_COLOR = PLT_FG_COLOR - LEGEND_BG_COLOR = P["fg1"] - LEGEND_BORDER_COLOR = P["fg2"] - COLORS = COLORS_DARK - if variant == "alt": - COLORS = COLORS_ALT - - color_rcParams = { - 'axes.edgecolor': PLT_FG_COLOR, - 'axes.facecolor': PLT_BG_COLOR, - 'axes.labelcolor': PLT_FG_COLOR, - 'axes.titlecolor': 'auto', - # 'axes.prop_cycle': cycler('color', ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']), - 'axes.prop_cycle': cycler('color', COLORS), - # 'axes3d.xaxis.panecolor': (0.95, 0.95, 0.95, 0.5), - # 'axes3d.yaxis.panecolor': (0.9, 0.9, 0.9, 0.5), - # 'axes3d.zaxis.panecolor': (0.925, 0.925, 0.925, 0.5), - # 'boxplot.boxprops.color': 'black', - # 'boxplot.capprops.color': 'black', - # 'boxplot.flierprops.color': 'black', - # 'boxplot.flierprops.markeredgecolor': 'black', - # 'boxplot.flierprops.markeredgewidth': 1.0, - # 'boxplot.flierprops.markerfacecolor': 'none', - # 'boxplot.meanprops.color': 'C2', - # 'boxplot.meanprops.markeredgecolor': 'C2', - # 'boxplot.meanprops.markerfacecolor': 'C2', - # 'boxplot.meanprops.markersize': 6.0, - # 'boxplot.medianprops.color': 'C1', - # 'boxplot.whiskerprops.color': 'black', - 'figure.edgecolor': PLT_BG_COLOR, - 'figure.facecolor': PLT_BG_COLOR, - # 'figure.figsize': [6.4, 4.8], - # 'figure.frameon': True, - # 'figure.labelsize': 'large', - 'grid.color': PLT_GRID_COLOR, - # 'hatch.color': 'black', - 'legend.edgecolor': LEGEND_BORDER_COLOR, - 'legend.facecolor': LEGEND_BG_COLOR, - 'xtick.color': PLT_FG_COLOR, - 'ytick.color': PLT_FG_COLOR, - 'xtick.labelcolor': PLT_FG_COLOR, - 'ytick.labelcolor': PLT_FG_COLOR, - # 'lines.color': 'C0', - 'text.color': PLT_FG_COLOR, - } - - for k, v in color_rcParams.items(): - plt.rcParams[k] = v - - # override single char codes - mpl.colors.get_named_colors_mapping()["b"] = COLORS[0] - mpl.colors.get_named_colors_mapping()["o"] = COLORS[1] - mpl.colors.get_named_colors_mapping()["g"] = COLORS[2] - mpl.colors.get_named_colors_mapping()["r"] = COLORS[3] - mpl.colors.get_named_colors_mapping()["m"] = COLORS[4] - mpl.colors.get_named_colors_mapping()["y"] = COLORS[5] - mpl.colors.get_named_colors_mapping()["c"] = COLORS[6] - mpl.colors.get_named_colors_mapping()["k"] = P["fg0"] - mpl.colors.get_named_colors_mapping()["w"] = P["bg0"] - diff --git a/src/mps.cpp b/src/mps.cpp index cc0664e..379188f 100644 --- a/src/mps.cpp +++ b/src/mps.cpp @@ -36,6 +36,24 @@ std::vector MPS::getChi() const { } return chis; } +template +unsigned MPS::getChiMax() const { + unsigned chiMax = this->B(0).dimension(2); + for (size_t i = 1; i < (L - 1); i++) { + if (this->B(i).dimension(2) > chiMax) { + chiMax = this->B(i).dimension(2); + } + } + return chiMax; +} +template +double MPS::getChiAverage() const { + double chiAvg = 0; + for (size_t i = 0; i < (L - 1); i++) { + chiAvg += this->B(i).dimension(2); + } + return chiAvg / (L-1); +} template T MPS::getSiteExpectationValue(const e::Tensor& op, size_t siteIdx) const { @@ -98,6 +116,7 @@ T MPS::collapse(size_t siteIndex) const { e::Tensor thetaTheta = theta.contract(theta.conjugate(), e::array{iP{0, 0}, iP{1, 1}}); for (size_t i = siteIndex+1; i < this->L; i++) { // vL vR vL* vR* + // TODO optimize by not assigning to tensor e::Tensor BB = this->B(i).contract(this->B(i).conjugate(), e::array{iP{1, 1}}); e::Tensor thetaTheta2 = thetaTheta.contract(BB, e::array{iP{0, 0}, iP{1, 2}}); std::swap(thetaTheta, thetaTheta2); @@ -139,15 +158,20 @@ split_truncate_theta(const e::Tensor& theta, size_t chiMax, dtype eps) e::Index dL = theta.dimension(1); e::Index dR = theta.dimension(2); e::Index chivR = theta.dimension(3); - - - // Reshape theta to 2D matrix for the SVD using MatX = e::Matrix; /// todo (DONE): Possible optimization: either use a Map and later adjust the matrix to tensor loops, or globally use RowMajor memory layout? // The map approach produces the same Matrix is tensor2Matrix with Variant 2 MatX theta_mat = e::Map(theta.data(), chivL * dL, dR * chivR); // MatX theta_mat = tensor2Matrix(theta); + return split_truncate_theta(theta_mat, chivL, dL, dR, chivR, chiMax, eps); +} + + +template +std::tuple, e::Tensor, e::Tensor> +split_truncate_theta(const e::MatrixX& theta, e::Index chivL, e::Index dL, e::Index dR, e::Index chivR, size_t chiMax, dtype eps) { + using MatX = e::Matrix; // std::cout << "theta:\n" << theta << "theta_mat\n": << theta_mat << std::endl; // std::cout << "theta_mat = " << std::endl; @@ -155,7 +179,7 @@ split_truncate_theta(const e::Tensor& theta, size_t chiMax, dtype eps) // perform SVD: theta = U * diag(S) * V^T constexpr int SVDOptions = e::ComputeThinU | e::ComputeThinV; - e::JacobiSVD svd(theta_mat); + e::BDCSVD svd(theta); // e::BDCSVD svd(theta_mat); const MatX& U = svd.matrixU(); const e::VectorX& S = svd.singularValues(); diff --git a/src/mps.hpp b/src/mps.hpp index a29aa2e..0d3ffd1 100644 --- a/src/mps.hpp +++ b/src/mps.hpp @@ -1,9 +1,11 @@ #pragma once #include "util.hpp" +#include #include #include #include #include +#include #include @@ -49,6 +51,10 @@ class MPS { * @brief Return the dimensions of all bonds */ std::vector getChi() const; + /// @brief Return the maximum bond dimension + unsigned getChiMax() const; + /// @brief Return the average bond dimension + double getChiAverage() const; /** * @brief Compute the expectation value of a single-site operator on a single site @@ -75,6 +81,7 @@ class MPS { STensor_t& S(size_t i) { return this->Ss.at(i); } const BTensor_t& B(size_t i) const { return this->Bs.at(i); } const STensor_t& S(size_t i) const { return this->Ss.at(i); } + // size_t getPhysicalDimension() const { return this->Bs.at(0).dimension(1); } size_t getL() const { return this->L; } void print() const { for (size_t i = 0; i < this->L; i++) { @@ -98,24 +105,49 @@ class MPS { * @returns (A, S, B) where A is the left-canonical matrix on site i (legs `vL, i, vC`), * S are the singular values and * B is the right-canonical matrix on site i+1, with legs `vC, j, vR` + * @note This functions maps the tensor and calls the Matrux overload of split_truncate_theta internally. */ template std::tuple, e::Tensor, e::Tensor> split_truncate_theta(const e::Tensor& theta, size_t chiMax, dtype eps); +template +std::tuple, e::Tensor, e::Tensor> +split_truncate_theta(const e::MatrixX& theta, e::Index chivL, e::Index dL, e::Index dR, e::Index chivR, size_t chiMax, dtype eps); + + +template +void applyPertubationOnB(std::vector>& Bs, double pertubationStrength) { + std::mt19937 generator(493852934); // Mersenne twister rng + std::uniform_real_distribution distribution(0.0, pertubationStrength); + for (auto it = Bs.begin(); it != Bs.end(); it++) { + double abs = 0; + for (unsigned i = 0; i < it->dimension(1); i++) { + (*it)(0, i, 0) += distribution(generator); + abs += std::pow((*it)(0, i, 0).real(), 2); + } + abs = std::sqrt(abs); + for (unsigned i = 0; i < it->dimension(1); i++) { + (*it)(0, i, 0) /= abs; + } + } +} + + /** - * @brief Initialize a MPS in spin-up state + * @brief Initialize a MPS using copies of the B tensor * @details - * Return a MPS spin-up MPS with length L and bond dimension 1 + * A random pertubation up to `pertubationStrength` can be applied to each occupation number in the B tensor, + * the tensors are normalized after the pertubations are applied. + * This can help break symmetries, which can make get the simulations stuck. */ template -MPS initMPS_spinUp(size_t L) { - assert(L > 0); - e::Tensor B(1, 2, 1); - B.setZero(); - B(0, 0, 0) = 1; +MPS makeMPS(const e::Tensor& B, size_t L, double pertubationStrength) { std::vector> Bs(L, B); + if (pertubationStrength > 0) { + applyPertubationOnB(Bs, pertubationStrength); + } e::Tensor S(1); S(0) = 1; std::vector> Ss(L, S); @@ -123,24 +155,43 @@ MPS initMPS_spinUp(size_t L) { } +/** + * @brief Initialize a MPS in spin-up state + * @details + * Return a MPS spin-up MPS with length L and bond dimension 1 + */ +template +MPS initMPS_spinUp(size_t L, double pertubationStrength=0) { + assert(L > 0); + e::Tensor B(1, 2, 1); + B.setZero(); + B(0, 0, 0) = 1; + return makeMPS(B, L, pertubationStrength); +} + +template +MPS initMPS_spinDown(size_t L, double pertubationStrength=0) { + assert(L > 0); + e::Tensor B(1, 2, 1); + B.setZero(); + B(0, 1, 0) = 1; + return makeMPS(B, L, pertubationStrength); +} + /** * @brief Initialize a MPS in spin-right state * @details * Return a MPS spin-right MPS with length L and bond dimension 1 */ template -MPS initMPS_spinRight(size_t L) { +MPS initMPS_spinRight(size_t L, double pertubationStrength=0) { assert(L > 0); e::Tensor B(1, 2, 1); B.setZero(); T value = std::sqrt(0.5); B(0, 0, 0) = value; B(0, 1, 0) = value; - std::vector> Bs(L, B); - e::Tensor S(1); - S(0) = 1; - std::vector> Ss(L, S); - return MPS(std::move(Bs), std::move(Ss)); + return makeMPS(B, L, pertubationStrength); } @@ -150,15 +201,11 @@ MPS initMPS_spinRight(size_t L) { * Return a MPS of dimension `nMax+1`, where each site is occupied by `nParticles`. */ template -MPS initMPS_nParticles(size_t L, size_t nMax, size_t nParticles) { +MPS initMPS_nParticles(size_t L, size_t nMax, size_t nParticles, double pertubationStrength=0) { assert(L > 0); assert(nParticles <= nMax); e::Tensor B(1, nMax+1, 1); B.setZero(); B(0, nParticles, 0) = 1; - std::vector> Bs(L, B); - e::Tensor S(1); - S(0) = 1; - std::vector> Ss(L, S); - return MPS(std::move(Bs), std::move(Ss)); + return makeMPS(B, L, pertubationStrength); } diff --git a/src/output b/src/output deleted file mode 100644 index 30aaaed..0000000 --- a/src/output +++ /dev/null @@ -1,562 +0,0 @@ -Building: main.cpp -> ../build/main.o -In file included from main.cpp:3: -mps.hpp: In instantiation of 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]': -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -mps.hpp:140:83: error: no match for 'operator>' (operand types are 'const Eigen::ArrayWrapper >' and 'std::complex') - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/bits/stl_algobase.h:64, - from /usr/include/c++/14.2.1/vector:62, - from ./pch.hpp:1, - from : -/usr/include/c++/14.2.1/bits/stl_pair.h:1023:5: note: candidate: 'template constexpr std::common_comparison_category_t(), declval<_U1&>())), decltype (std::__detail::__synth3way(declval<_T2&>(), declval<_U2&>()))> std::operator<=>(const pair<_T1, _T2>&, const pair<_U1, _U2>&)' (reversed) - 1023 | operator<=>(const pair<_T1, _T2>& __x, const pair<_U1, _U2>& __y) - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_pair.h:1023:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::pair<_T1, _T2>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/bits/stl_algobase.h:67: -/usr/include/c++/14.2.1/bits/stl_iterator.h:576:5: note: candidate: 'template requires three_way_comparable_with<_IteratorR, _IteratorL, std::partial_ordering> constexpr std::compare_three_way_result_t<_IteratorL, _IteratorR> std::operator<=>(const reverse_iterator<_IteratorL>&, const reverse_iterator<_IteratorR>&)' (reversed) - 576 | operator<=>(const reverse_iterator<_IteratorL>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:576:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::reverse_iterator<_IteratorL>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:1679:5: note: candidate: 'template requires three_way_comparable_with<_IteratorR, _IteratorL, std::partial_ordering> constexpr std::compare_three_way_result_t<_IteratorL, _IteratorR> std::operator<=>(const move_iterator<_IteratorL>&, const move_iterator<_IteratorR>&)' (reversed) - 1679 | operator<=>(const move_iterator<_IteratorL>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:1679:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::move_iterator<_IteratorL>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/bits/uses_allocator_args.h:39, - from /usr/include/c++/14.2.1/bits/memory_resource.h:41, - from /usr/include/c++/14.2.1/vector:87: -/usr/include/c++/14.2.1/tuple:2589:5: note: candidate: 'template constexpr std::common_comparison_category_t(), declval<_Ups&>()))...> std::operator<=>(const tuple<_Ts ...>&, const tuple<_Elements ...>&)' (reversed) - 2589 | operator<=>(const tuple<_Tps...>& __t, const tuple<_Ups...>& __u) - | ^~~~~~~~ -/usr/include/c++/14.2.1/tuple:2589:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::tuple<_Ts ...>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/bits/basic_string.h:47, - from /usr/include/c++/14.2.1/string:54, - from /usr/include/c++/14.2.1/bits/locale_classes.h:40, - from /usr/include/c++/14.2.1/locale:41, - from /usr/include/c++/14.2.1/format:47, - from /usr/include/c++/14.2.1/print:41, - from ./pch.hpp:2: -/usr/include/c++/14.2.1/string_view:617:5: note: candidate: 'template constexpr decltype (__char_traits_cmp_cat<_Traits>(0)) std::operator<=>(basic_string_view<_CharT, _Traits>, __type_identity_t >)' (reversed) - 617 | operator<=>(basic_string_view<_CharT, _Traits> __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/string_view:617:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'std::basic_string_view<_CharT, _Traits>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/basic_string.h:3805:5: note: candidate: 'template constexpr decltype (__char_traits_cmp_cat<_Traits>(0)) std::operator<=>(const __cxx11::basic_string<_CharT, _Traits, _Allocator>&, const _CharT*)' (reversed) - 3805 | operator<=>(const basic_string<_CharT, _Traits, _Alloc>& __lhs, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/basic_string.h:3805:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::__cxx11::basic_string<_CharT, _Traits, _Allocator>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/format:48: -/usr/include/c++/14.2.1/optional:1283:5: note: candidate: 'template requires three_way_comparable_with<_Up, _Tp, std::partial_ordering> constexpr std::compare_three_way_result_t<_IteratorL, _IteratorR> std::operator<=>(const optional<_Tp>&, const optional<_Up>&)' (reversed) - 1283 | operator<=>(const optional<_Tp>& __x, const optional<_Up>& __y) - | ^~~~~~~~ -/usr/include/c++/14.2.1/optional:1283:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::optional<_Tp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/optional:1298:5: note: candidate: 'template constexpr std::strong_ordering std::operator<=>(const optional<_Tp>&, nullopt_t)' (reversed) - 1298 | operator<=>(const optional<_Tp>& __x, nullopt_t) noexcept - | ^~~~~~~~ -/usr/include/c++/14.2.1/optional:1298:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::optional<_Tp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/optional:1435:5: note: candidate: 'template requires !(__is_optional_v<_Up>) && (three_way_comparable_with<_Up, _Tp, std::partial_ordering>) constexpr std::compare_three_way_result_t<_IteratorL, _IteratorR> std::operator<=>(const optional<_Tp>&, const _Up&)' (reversed) - 1435 | operator<=>(const optional<_Tp>& __x, const _Up& __v) - | ^~~~~~~~ -/usr/include/c++/14.2.1/optional:1435:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::optional<_Tp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/bits/shared_ptr_base.h:59, - from /usr/include/c++/14.2.1/bits/shared_ptr.h:53, - from /usr/include/c++/14.2.1/chrono:49, - from /usr/include/eigen3/unsupported/Eigen/CXX11/Tensor:38, - from ./pch.hpp:4: -/usr/include/c++/14.2.1/bits/unique_ptr.h:996:5: note: candidate: 'template requires three_way_comparable_with::pointer, typename std::unique_ptr<_Up, _Ep>::pointer, std::partial_ordering> constexpr std::compare_three_way_result_t::pointer, typename std::unique_ptr<_Up, _Ep>::pointer> std::operator<=>(const unique_ptr<_Tp, _Dp>&, const unique_ptr<_Up, _Ep>&)' (reversed) - 996 | operator<=>(const unique_ptr<_Tp, _Dp>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:996:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::unique_ptr<_Tp, _Dp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:1005:5: note: candidate: 'template requires three_way_comparable::pointer, std::partial_ordering> constexpr std::compare_three_way_result_t::pointer> std::operator<=>(const unique_ptr<_Tp, _Dp>&, nullptr_t)' (reversed) - 1005 | operator<=>(const unique_ptr<_Tp, _Dp>& __x, nullptr_t) - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:1005:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::unique_ptr<_Tp, _Dp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/shared_ptr_base.h:1809:5: note: candidate: 'template std::strong_ordering std::operator<=>(const __shared_ptr<_Tp1, _Lp>&, const __shared_ptr<_Tp2, _Lp>&)' (reversed) - 1809 | operator<=>(const __shared_ptr<_Tp, _Lp>& __a, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/shared_ptr_base.h:1809:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::__shared_ptr<_Tp1, _Lp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/shared_ptr_base.h:1815:5: note: candidate: 'template std::strong_ordering std::operator<=>(const __shared_ptr<_Tp, _Lp>&, nullptr_t)' (reversed) - 1815 | operator<=>(const __shared_ptr<_Tp, _Lp>& __a, nullptr_t) noexcept - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/shared_ptr_base.h:1815:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::__shared_ptr<_Tp, _Lp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/shared_ptr.h:566:5: note: candidate: 'template std::strong_ordering std::operator<=>(const shared_ptr<_Tp>&, const shared_ptr<_Tp>&)' (reversed) - 566 | operator<=>(const shared_ptr<_Tp>& __a, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/shared_ptr.h:566:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::shared_ptr<_Tp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/shared_ptr.h:572:5: note: candidate: 'template std::strong_ordering std::operator<=>(const shared_ptr<_Tp>&, nullptr_t)' (reversed) - 572 | operator<=>(const shared_ptr<_Tp>& __a, nullptr_t) noexcept - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/shared_ptr.h:572:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::shared_ptr<_Tp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:594:5: note: candidate: 'template requires three_way_comparable<_Iterator, std::partial_ordering> constexpr std::compare_three_way_result_t<_Iterator, _Iterator> std::operator<=>(const reverse_iterator<_IteratorL>&, const reverse_iterator<_IteratorL>&)' (rewritten) - 594 | operator<=>(const reverse_iterator<_Iterator>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:594:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::reverse_iterator<_IteratorL>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:1745:5: note: candidate: 'template requires three_way_comparable<_Iterator, std::partial_ordering> constexpr std::compare_three_way_result_t<_Iterator, _Iterator> std::operator<=>(const move_iterator<_IteratorL>&, const move_iterator<_IteratorL>&)' (rewritten) - 1745 | operator<=>(const move_iterator<_Iterator>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:1745:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::move_iterator<_IteratorL>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/vector:66: -/usr/include/c++/14.2.1/bits/stl_vector.h:2069:5: note: candidate: 'template constexpr std::__detail::__synth3way_t<_Iterator> std::operator<=>(const vector<_Tp, _Alloc>&, const vector<_Tp, _Alloc>&)' (rewritten) - 2069 | operator<=>(const vector<_Tp, _Alloc>& __x, const vector<_Tp, _Alloc>& __y) - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_vector.h:2069:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::vector<_Tp, _Alloc>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/format:43: -/usr/include/c++/14.2.1/array:310:5: note: candidate: 'template constexpr std::__detail::__synth3way_t<_Iterator> std::operator<=>(const array<_Tp, _Nm>&, const array<_Tp, _Nm>&)' (rewritten) - 310 | operator<=>(const array<_Tp, _Nm>& __a, const array<_Tp, _Nm>& __b) - | ^~~~~~~~ -/usr/include/c++/14.2.1/array:310:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::array<_Tp, _Nm>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/basic_string.h:3790:5: note: candidate: 'template constexpr decltype (__char_traits_cmp_cat<_Traits>(0)) std::operator<=>(const __cxx11::basic_string<_CharT, _Traits, _Allocator>&, const __cxx11::basic_string<_CharT, _Traits, _Allocator>&)' (rewritten) - 3790 | operator<=>(const basic_string<_CharT, _Traits, _Alloc>& __lhs, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/basic_string.h:3790:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::__cxx11::basic_string<_CharT, _Traits, _Allocator>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/format:52: -/usr/include/c++/14.2.1/variant:1280:5: note: candidate: 'template requires (three_way_comparable<_Types, std::partial_ordering> && ...) constexpr std::common_comparison_category_t::type ...> std::operator<=>(const variant<_Types ...>&, const variant<_Types ...>&)' (rewritten) - 1280 | operator<=>(const variant<_Types...>& __v, const variant<_Types...>& __w) - | ^~~~~~~~ -/usr/include/c++/14.2.1/variant:1280:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::variant<_Types ...>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/map:63, - from /usr/include/eigen3/Eigen/SparseCore:16, - from /usr/include/eigen3/Eigen/Sparse:26, - from /usr/include/eigen3/Eigen/Eigen:2, - from ./pch.hpp:5: -/usr/include/c++/14.2.1/bits/stl_map.h:1533:5: note: candidate: 'template std::__detail::__synth3way_t > std::operator<=>(const map<_Key, _Tp, _Compare, _Allocator>&, const map<_Key, _Tp, _Compare, _Allocator>&)' (rewritten) - 1533 | operator<=>(const map<_Key, _Tp, _Compare, _Alloc>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_map.h:1533:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::map<_Key, _Tp, _Compare, _Allocator>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/map:64: -/usr/include/c++/14.2.1/bits/stl_multimap.h:1155:5: note: candidate: 'template std::__detail::__synth3way_t > std::operator<=>(const multimap<_Key, _Tp, _Compare, _Allocator>&, const multimap<_Key, _Tp, _Compare, _Allocator>&)' (rewritten) - 1155 | operator<=>(const multimap<_Key, _Tp, _Compare, _Alloc>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_multimap.h:1155:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::multimap<_Key, _Tp, _Compare, _Allocator>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/list:65, - from /usr/include/eigen3/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h:15, - from /usr/include/eigen3/Eigen/IterativeLinearSolvers:47, - from /usr/include/eigen3/Eigen/Sparse:31: -/usr/include/c++/14.2.1/bits/stl_list.h:2158:5: note: candidate: 'template std::__detail::__synth3way_t<_Iterator> std::operator<=>(const __cxx11::list<_Tp, _Alloc>&, const __cxx11::list<_Tp, _Alloc>&)' (rewritten) - 2158 | operator<=>(const list<_Tp, _Alloc>& __x, const list<_Tp, _Alloc>& __y) - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_list.h:2158:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::__cxx11::list<_Tp, _Alloc>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -In file included from /usr/include/c++/14.2.1/bits/ios_base.h:46, - from /usr/include/c++/14.2.1/bits/locale_facets.h:43, - from /usr/include/c++/14.2.1/locale:42: -/usr/include/c++/14.2.1/system_error:316:3: note: candidate: 'std::strong_ordering std::operator<=>(const error_code&, const error_code&)' (rewritten) - 316 | operator<=>(const error_code& __lhs, const error_code& __rhs) noexcept - | ^~~~~~~~ -/usr/include/c++/14.2.1/system_error:316:33: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper >' to 'const std::error_code&' - 316 | operator<=>(const error_code& __lhs, const error_code& __rhs) noexcept - | ~~~~~~~~~~~~~~~~~~^~~~~ -/usr/include/c++/14.2.1/system_error:498:3: note: candidate: 'std::strong_ordering std::operator<=>(const error_condition&, const error_condition&)' (rewritten) - 498 | operator<=>(const error_condition& __lhs, - | ^~~~~~~~ -/usr/include/c++/14.2.1/system_error:498:38: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper >' to 'const std::error_condition&' - 498 | operator<=>(const error_condition& __lhs, - | ~~~~~~~~~~~~~~~~~~~~~~~^~~~~ -/usr/include/c++/14.2.1/variant:1303:3: note: candidate: 'constexpr std::strong_ordering std::operator<=>(monostate, monostate)' (rewritten) - 1303 | operator<=>(monostate, monostate) noexcept { return strong_ordering::equal; } - | ^~~~~~~~ -/usr/include/c++/14.2.1/variant:1303:15: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper >' to 'std::monostate' - 1303 | operator<=>(monostate, monostate) noexcept { return strong_ordering::equal; } - | ^~~~~~~~~ -In file included from /usr/include/eigen3/unsupported/Eigen/CXX11/Tensor:45: -/usr/include/c++/14.2.1/thread:75:3: note: candidate: 'std::strong_ordering std::operator<=>(thread::id, thread::id)' (rewritten) - 75 | operator<=>(thread::id __x, thread::id __y) noexcept - | ^~~~~~~~ -/usr/include/c++/14.2.1/thread:75:26: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper >' to 'std::thread::id' - 75 | operator<=>(thread::id __x, thread::id __y) noexcept - | ~~~~~~~~~~~^~~ -In file included from /usr/include/eigen3/Eigen/src/Core/ArrayBase.h:99, - from /usr/include/eigen3/Eigen/Core:309, - from /usr/include/eigen3/unsupported/Eigen/CXX11/Tensor:14: -/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: candidate: 'template const Eigen::CwiseBinaryOp::Scalar, Eigen::internal::cmp_LT>, const OtherDerived, const Derived> Eigen::ArrayBase::operator>(const Eigen::ArrayBase&) const [with Derived = Eigen::ArrayWrapper >]' - 204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT) - | ^~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const Eigen::ArrayBase' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: candidate: 'const Eigen::ArrayBase::RCmpLTReturnType Eigen::ArrayBase::operator>(const Scalar&) const [with Derived = Eigen::ArrayWrapper >; RCmpLTReturnType = Eigen::CwiseBinaryOp, const Eigen::CwiseNullaryOp, Eigen::Array >, const Eigen::ArrayWrapper > >; Scalar = double]' - 204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT) - | ^~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: no known conversion for argument 1 from 'std::complex' to 'const Eigen::ArrayBase > >::Scalar&' {aka 'const double&'} - 204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT) - | ^~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:551:5: note: candidate: 'template constexpr bool std::operator>(const reverse_iterator<_IteratorL>&, const reverse_iterator<_IteratorR>&) requires requires{{std::operator>::__x->base() < std::operator>::__y->base()} -> decltype(auto) [requires std::convertible_to<, bool>];}' - 551 | operator>(const reverse_iterator<_IteratorL>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:551:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::reverse_iterator<_IteratorL>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:1714:5: note: candidate: 'template constexpr bool std::operator>(const move_iterator<_IteratorL>&, const move_iterator<_IteratorR>&) requires requires{{std::operator>::__y->base() < std::operator>::__x->base()} -> decltype(auto) [requires std::convertible_to<, bool>];}' - 1714 | operator>(const move_iterator<_IteratorL>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/stl_iterator.h:1714:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::move_iterator<_IteratorL>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/optional:1258:5: note: candidate: 'template constexpr std::__optional_gt_t<_Tp, _Up> std::operator>(const optional<_Tp>&, const optional<_Up>&)' - 1258 | operator>(const optional<_Tp>& __lhs, const optional<_Up>& __rhs) - | ^~~~~~~~ -/usr/include/c++/14.2.1/optional:1258:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::optional<_Tp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/optional:1396:5: note: candidate: 'template constexpr std::__optional_gt_t<_Tp, _Up> std::operator>(const optional<_Tp>&, const _Up&)' - 1396 | operator>(const optional<_Tp>& __lhs, const _Up& __rhs) - | ^~~~~~~~ -/usr/include/c++/14.2.1/optional:1396:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::optional<_Tp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/optional:1402:5: note: candidate: 'template constexpr std::__optional_gt_t<_Up, _Tp> std::operator>(const _Up&, const optional<_Tp>&)' - 1402 | operator>(const _Up& __lhs, const optional<_Tp>& __rhs) - | ^~~~~~~~ -/usr/include/c++/14.2.1/optional:1402:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::optional<_Tp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/variant:1269:3: note: candidate: 'template constexpr bool std::operator>(const variant<_Types ...>&, const variant<_Types ...>&)' - 1269 | _VARIANT_RELATION_FUNCTION_TEMPLATE(>, greater) - | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/c++/14.2.1/variant:1269:3: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::variant<_Types ...>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:942:5: note: candidate: 'template constexpr bool std::operator>(const unique_ptr<_Tp, _Dp>&, const unique_ptr<_Up, _Ep>&)' - 942 | operator>(const unique_ptr<_Tp, _Dp>& __x, - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:942:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::unique_ptr<_Tp, _Dp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:950:5: note: candidate: 'template constexpr bool std::operator>(const unique_ptr<_Tp, _Dp>&, nullptr_t)' - 950 | operator>(const unique_ptr<_Tp, _Dp>& __x, nullptr_t) - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:950:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'const Eigen::ArrayWrapper >' is not derived from 'const std::unique_ptr<_Tp, _Dp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:960:5: note: candidate: 'template constexpr bool std::operator>(nullptr_t, const unique_ptr<_Tp, _Dp>&)' - 960 | operator>(nullptr_t, const unique_ptr<_Tp, _Dp>& __x) - | ^~~~~~~~ -/usr/include/c++/14.2.1/bits/unique_ptr.h:960:5: note: template argument deduction/substitution failed: -mps.hpp:140:83: note: 'std::complex' is not derived from 'const std::unique_ptr<_Tp, _Dp>' - 140 | int chivC = static_cast(std::min(chi_max, static_cast((S.array() > eps).count()))); - | ~~~~~~~~~~~^~~~~~ -/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: candidate: 'const Eigen::ArrayBase > >::CmpLTReturnType Eigen::operator>(const ArrayBase > >::Scalar&, const ArrayWrapper >&)' - 204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT) - | ^~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.inc:204:1: note: no known conversion for argument 1 from 'const Eigen::ArrayWrapper >' to 'const Eigen::ArrayBase > >::Scalar&' {aka 'const double&'} - 204 | EIGEN_MAKE_CWISE_COMP_R_OP(operator>, operator<, LT) - | ^~~~~~~~~~~~~~~~~~~~~~~~~~ -In file included from /usr/include/eigen3/Eigen/Core:300: -/usr/include/eigen3/Eigen/src/Core/Product.h: In instantiation of 'struct Eigen::internal::traits, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >': -/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1, 8>' - 248 | class dense_product_base : public internal::dense_xpr_base>::type {}; - | ^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1, Eigen::Dense>' - 276 | class ProductImpl : public internal::dense_product_base { - | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>' - 198 | class Product - | ^~~~~~~ -mps.hpp:159:21: required from 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]' - 159 | MatX theta2 = U * S.asDiagonal() * V; - | ~~^~~~~~~~~~~~~~~~ -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:33:90: error: no type named 'ReturnType' in 'struct Eigen::ScalarBinaryOpTraits, double, Eigen::internal::scalar_product_op, double> >' - 33 | typename traits::Scalar>::ReturnType Scalar; - | ^~~~~~ -In file included from /usr/include/eigen3/Eigen/Core:296: -/usr/include/eigen3/Eigen/src/Core/DenseBase.h: In instantiation of 'class Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >': -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:52:7: required from 'class Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >' - 52 | class MatrixBase : public DenseBase { - | ^~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1, 8>' - 248 | class dense_product_base : public internal::dense_xpr_base>::type {}; - | ^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1, Eigen::Dense>' - 276 | class ProductImpl : public internal::dense_product_base { - | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>' - 198 | class Product - | ^~~~~~~ -mps.hpp:159:21: required from 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]' - 159 | MatX theta2 = U * S.asDiagonal() * V; - | ~~^~~~~~~~~~~~~~~~ -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:72:15: error: 'coeff' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 72 | using Base::coeff; - | ^~~~~ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:73:15: error: 'coeffByOuterInner' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 73 | using Base::coeffByOuterInner; - | ^~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:81:24: error: 'operator()' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 81 | using Base::operator(); - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:82:24: error: 'operator[]' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 82 | using Base::operator[]; - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:88:15: error: 'w' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 88 | using Base::w; - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:89:15: error: 'x' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 89 | using Base::x; - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:90:15: error: 'y' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 90 | using Base::y; - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:91:15: error: 'z' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 91 | using Base::z; - | ^ -In file included from /usr/include/eigen3/Eigen/Core:297: -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h: In instantiation of 'class Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >': -/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1, 8>' - 248 | class dense_product_base : public internal::dense_xpr_base>::type {}; - | ^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1, Eigen::Dense>' - 276 | class ProductImpl : public internal::dense_product_base { - | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>' - 198 | class Product - | ^~~~~~~ -mps.hpp:159:21: required from 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]' - 159 | MatX theta2 = U * S.asDiagonal() * V; - | ~~^~~~~~~~~~~~~~~~ -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:72:15: error: 'coeff' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 72 | using Base::coeff; - | ^~~~~ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:77:15: error: 'eval' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 77 | using Base::eval; - | ^~~~ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:81:23: error: 'operator-' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 81 | using Base::operator-; - | ^ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:84:23: error: 'operator*=' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 84 | using Base::operator*=; - | ^~ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:85:23: error: 'operator/=' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1> >::Base' - 85 | using Base::operator/=; - | ^~ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h: In instantiation of 'class Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >': -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:52:7: required from 'class Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >' - 52 | class MatrixBase : public DenseBase { - | ^~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0, 8>' - 248 | class dense_product_base : public internal::dense_xpr_base>::type {}; - | ^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0, Eigen::Dense>' - 276 | class ProductImpl : public internal::dense_product_base { - | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0>' - 198 | class Product - | ^~~~~~~ -mps.hpp:159:38: required from 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]' - 159 | MatX theta2 = U * S.asDiagonal() * V; - | ~~~~~~~~~~~~~~~~~~~^~~ -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:72:15: error: 'coeff' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 72 | using Base::coeff; - | ^~~~~ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:73:15: error: 'coeffByOuterInner' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 73 | using Base::coeffByOuterInner; - | ^~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:81:24: error: 'operator()' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 81 | using Base::operator(); - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:82:24: error: 'operator[]' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 82 | using Base::operator[]; - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:88:15: error: 'w' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 88 | using Base::w; - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:89:15: error: 'x' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 89 | using Base::x; - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:90:15: error: 'y' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 90 | using Base::y; - | ^ -/usr/include/eigen3/Eigen/src/Core/DenseBase.h:91:15: error: 'z' has not been declared in 'Eigen::DenseBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 91 | using Base::z; - | ^ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h: In instantiation of 'class Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >': -/usr/include/eigen3/Eigen/src/Core/Product.h:248:7: required from 'class Eigen::internal::dense_product_base, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0, 8>' - 248 | class dense_product_base : public internal::dense_xpr_base>::type {}; - | ^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:276:7: required from 'class Eigen::ProductImpl, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0, Eigen::Dense>' - 276 | class ProductImpl : public internal::dense_product_base { - | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Product.h:198:7: required from 'class Eigen::Product, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0>' - 198 | class Product - | ^~~~~~~ -mps.hpp:159:38: required from 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]' - 159 | MatX theta2 = U * S.asDiagonal() * V; - | ~~~~~~~~~~~~~~~~~~~^~~ -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:72:15: error: 'coeff' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 72 | using Base::coeff; - | ^~~~~ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:77:15: error: 'eval' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 77 | using Base::eval; - | ^~~~ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:81:23: error: 'operator-' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 81 | using Base::operator-; - | ^ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:84:23: error: 'operator*=' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 84 | using Base::operator*=; - | ^~ -/usr/include/eigen3/Eigen/src/Core/MatrixBase.h:85:23: error: 'operator/=' has not been declared in 'Eigen::MatrixBase, -1, -1, 1, -1, -1>, Eigen::DiagonalWrapper >, 1>, Eigen::Matrix, -1, -1, 1, -1, -1>, 0> >::Base' - 85 | using Base::operator/=; - | ^~ -/usr/include/c++/14.2.1/format: In instantiation of 'class std::__format::_Checking_scanner >': -/usr/include/c++/14.2.1/format:4210:4: required from 'consteval std::basic_format_string<_CharT, _Args>::basic_format_string(const _Tp&) [with _Tp = char [4]; _CharT = char; _Args = {std::complex&}]' - 4210 | __scanner(_M_str); - | ^~~~~~~~~ -mps.hpp:162:19: required from 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]' - 162 | std::print("{} ", theta2.data()[i]); - | ~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~ -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -/usr/include/c++/14.2.1/format:4065:10: error: static assertion failed: std::formatter must be specialized for each type being formatted - 4065 | (is_default_constructible_v> && ...), - | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/c++/14.2.1/format:4065:10: note: 'std::is_default_constructible_v, char> >' evaluates to false -/usr/include/c++/14.2.1/format: In instantiation of 'constexpr void std::__format::_Checking_scanner<_CharT, _Args>::_M_parse_format_spec(std::size_t) [with _Tp = std::complex; _OtherArgs = {}; _CharT = char; _Args = {std::complex}; std::size_t = long unsigned int]': -/usr/include/c++/14.2.1/format:4082:33: required from 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]' - 4082 | _M_parse_format_spec<_Args...>(__id); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~ -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -mps.hpp:162:19: in 'constexpr' expansion of 'std::basic_format_string&>("{} ")' -/usr/include/c++/14.2.1/format:4211:19: in 'constexpr' expansion of '__scanner.std::__format::_Checking_scanner >::std::__format::_Scanner.std::__format::_Scanner::_M_scan()' -/usr/include/c++/14.2.1/format:3952:30: in 'constexpr' expansion of '((std::__format::_Scanner*)this)->std::__format::_Scanner::_M_on_replacement_field()' -/usr/include/c++/14.2.1/format:4004:15: in 'constexpr' expansion of '((std::__format::_Scanner*)this)->std::__format::_Scanner::_M_format_arg(__id)' -/usr/include/c++/14.2.1/format:4095:38: error: use of deleted function 'std::formatter<_Tp, _CharT>::formatter() [with _Tp = std::complex; _CharT = char]' - 4095 | formatter<_Tp, _CharT> __f; - | ^~~ -/usr/include/c++/14.2.1/format:178:7: note: declared here - 178 | formatter() = delete; // No std::formatter specialization for this type. - | ^~~~~~~~~ -/usr/include/c++/14.2.1/format:4095:38: note: use '-fdiagnostics-all-candidates' to display considered candidates - 4095 | formatter<_Tp, _CharT> __f; - | ^~~ -/usr/include/c++/14.2.1/format:4096:42: error: 'class std::formatter, char>' has no member named 'parse' - 4096 | this->_M_pc.advance_to(__f.parse(this->_M_pc)); - | ~~~~^~~~~ -mps.hpp: In instantiation of 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]': -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -mps.hpp:162:19: error: call to consteval function 'std::basic_format_string&>("{} ")' is not a constant expression - 162 | std::print("{} ", theta2.data()[i]); - | ~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~ -In file included from /usr/include/eigen3/Eigen/Core:175: -/usr/include/eigen3/Eigen/src/Core/AssignEvaluator.h: In instantiation of 'constexpr void Eigen::internal::call_assignment_no_alias(Dst&, const Src&, const Func&) [with Dst = Eigen::Matrix; Src = Eigen::Matrix; Func = assign_op]': -/usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:770:39: required from 'constexpr Derived& Eigen::PlainObjectBase::_set_noalias(const Eigen::DenseBase&) [with OtherDerived = Eigen::Matrix; Derived = Eigen::Matrix]' - 770 | internal::call_assignment_no_alias(this->derived(), other.derived(), - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - 771 | internal::assign_op()); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:565:17: required from 'Eigen::PlainObjectBase::PlainObjectBase(const Eigen::DenseBase&) [with OtherDerived = Eigen::Matrix; Derived = Eigen::Matrix]' - 565 | _set_noalias(other); - | ~~~~~~~~~~~~^~~~~~~ -/usr/include/eigen3/Eigen/src/Core/Matrix.h:386:108: required from 'Eigen::Matrix::Matrix(const Eigen::EigenBase&) [with OtherDerived = Eigen::Matrix; Scalar_ = double; int Rows_ = -1; int Cols_ = 1; int Options_ = 0; int MaxRows_ = -1; int MaxCols_ = 1]' - 386 | EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Matrix(const EigenBase& other) : Base(other.derived()) {} - | ^ -mps.hpp:136:24: required from 'std::tuple, Eigen::Tensor, Eigen::Tensor > split_truncate_theta(const Eigen::Tensor&, size_t, dtype) [with dtype = std::complex; unsigned int Storage = 1; size_t = long unsigned int]' - 136 | const e::VectorXd& S = svd.singularValues(); - | ^ -main.cpp:37:62: required from here - 37 | auto [A, S, B] = split_truncate_theta(t4, chi_max, eps); - | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/AssignEvaluator.h:839:3: error: static assertion failed: YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY - 839 | EIGEN_CHECK_BINARY_COMPATIBILIY(Func, typename ActualDstTypeCleaned::Scalar, typename Src::Scalar); - | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -/usr/include/eigen3/Eigen/src/Core/AssignEvaluator.h:839:3: note: 'Eigen::internal::has_ReturnType > >::value' evaluates to false -make: *** [Makefile:47: ../build/main.o] Error 1 diff --git a/src/run.cpp b/src/run.cpp deleted file mode 100644 index 3b18903..0000000 --- a/src/run.cpp +++ /dev/null @@ -1,390 +0,0 @@ -#include -#include -#ifdef EMSCRIPTEN -#include -#include -// #include -#include "mps.hpp" -#include "transverse_field_ising.hpp" -#include "bose_hubbard.hpp" -#include "tebd.hpp" -#include "util.hpp" -#include "type.hpp" -#include -#include -#include -#include -namespace ems = emscripten; - -using TfiModel = TransverseFieldIsingModel1D; -using BhModel = BoseHubbardModel1D; - -/** - * @brief A class for generating a range of values. - * @details - * Returns `nSteps` values with `start` being the first and `stop` being the last value. - * The same value is returned `countsPerStep` times before the next value is returned. - * After the last `stop` value has been returned, the cycle is restarted from the first `start` value. - */ -template -struct Range { - using value_type = T; - Range(T start, T stop, unsigned nSteps, unsigned countsPerStep=1) : start(start), stop(stop), nSteps(nSteps), stepCounter(0), countsPerStep(countsPerStep), counter(0) { - if (nSteps <= 1) { - throw std::invalid_argument("nSteps must be greater than 1"); - } - step = (stop - start) / (nSteps - 1); - }; - T getNext() { - T result = start + stepCounter * step; - if (++counter >= countsPerStep) { - counter = 0; - if (++stepCounter >= nSteps) { - stepCounter = 0; - } - } - return result; - }; - private: - T start; - T stop; - T step; - unsigned nSteps; - unsigned stepCounter; - unsigned countsPerStep; - unsigned counter; -}; - -/** - * This ugly mechanism lets the javascript pass - * any combination of Range and basic types to a function. - * Every time the function is called, the Range classes will return a new value while the rest stays the same. - */ -template class U> -inline constexpr bool is_instance_of_v = std::false_type{}; -template class U, class V> -inline constexpr bool is_instance_of_v,U> = std::true_type{}; - -template -concept isRange = is_instance_of_v, Range>; - - -// All types except Range -template -struct transform_type { - using type = T; -}; -// Range -template -struct transform_type>> { - using type = typename std::remove_reference_t::value_type; // If T is a Range, return its value_type -}; - - -static_assert(isRange>); -static_assert(isRange&>); -static_assert(!isRange); - -// Dispatcher -template -std::remove_reference_t::value_type processType(T&& obj) { - return obj.getNext(); -} -template -requires(!isRange) -T processType(T obj) { - return obj; -} - -/// Make a shared model pointer from a arguments, which may contain Range objects -template -std::shared_ptr makeModel(Args&&... args) { - std::tuple::type...> results = std::make_tuple(processType(std::forward(args))...); - return std::apply([](auto&&... args) { return std::make_shared(args...); }, results); -} - -// template -// void printArgs(Arg&& val) { -// std::cout << val << std::endl; -// } -// template -// void printArgs(Arg&& val, Args&&... values) { -// std::cout << val << ", "; -// printArgs(std::forward(values)...); -// } - -// template -// void processDebug(Args&&... args) { -// std::tuple::type...> results = std::make_tuple(processType(std::forward(args))...); -// std::apply([](auto&&... args){ printArgs(args...); }, results); -// } - -/** - */ -template -class Results { - public: - virtual ~Results() {}; - virtual postUpdateFunction_t getPostUpdateFunction(std::shared_ptr model) = 0; - virtual void postUpdateFunction(std::shared_ptr model, const MPS& psi) = 0; -}; - - -class TfiResults : public Results { - public: - using Model_t = TfiModel; - TfiResults(ems::val callback) : callback(callback) { - sigmaZ.setValues({{1, 0}, {0, -1}}); - }; - /// @todo write a second postUpdateFunction without the JS callback and return that one if callback is null - postUpdateFunction_t getPostUpdateFunction(std::shared_ptr model) override { - return [model, this](const MPS& psi, unsigned i){ this->postUpdateFunction(model, psi); }; - } - void postUpdateFunction(std::shared_ptr model, const MPS& psi) override { - auto Mval = psi.getSiteExpectationValueSum(this->sigmaZ).real(); - auto Eval = model->energy(psi).real(); - auto Sval = psi.getEntanglementEntropy(psi.getL()/2).real(); - callback(Mval, Eval, Sval); - // hand control back to browser, otherwise site hangs and callback will not be executed - emscripten_sleep(0); - } - private: - e::TensorFixedSize> sigmaZ; - ems::val callback; -}; - - -/** - * @brief Measures , E and S - */ -class BoseHubbardResults : public Results { - public: - BoseHubbardResults(ems::val callback) : callback(callback) {}; - postUpdateFunction_t getPostUpdateFunction(std::shared_ptr model) override { - return [this, model](const MPS& psi, unsigned) { this->postUpdateFunction(model, psi); }; - } - void postUpdateFunction(std::shared_ptr model, const MPS& psi) override { - e::TensorMap> particleNumberOperator(model->getLadderOperators().getN().data(), model->localDim, model->localDim); - auto Nval = psi.getSiteExpectationValueSum(particleNumberOperator).real() / model->L; - auto Eval = model->energy(psi).real(); - auto Sval = psi.getEntanglementEntropy(model->L/2).real(); - callback(Nval, Eval, Sval); - // hand control back to browser, otherwise site hangs and callback will not be executed - emscripten_sleep(0); - } - private: - ems::val callback; -}; - -// EMSCRIPTEN_BINDINGS(my_class_example) { -// class_("MyClass") -// .constructor() -// .function("incrementX", &MyClass::incrementX) -// .property("x", &MyClass::getX, &MyClass::setX) -// .property("x_readonly", &MyClass::getX) -// .class_function("getStringFromInstance", &MyClass::getStringFromInstance) -// ; -// }; - - -enum Tfi_INIT { - INIT_SPIN_UP = 0, -}; -// struct Tfi { -// Tfi(unsigned L, double _g, double _J, Tfi_INIT init) { -// throwIfLessThan(L, 2u, "L"); -// dtype g(_g); -// dtype J(_J); -// throwIfNaN(g, "g"); -// throwIfNaN(J, "J"); -// this->model = TfiModel(L, J, g); -// if (init == INIT_SPIN_UP) { -// mps = initMPSspinUp(L); -// } - -// } -// TfiModel model; -// MPS mps; -// }; - -// TODO unify -// Model1D initModel_Tfi(unsigned L, double _g, double _J) { -// dtype g(_g); -// dtype J(_J); -// throwIfNaN(g, "g"); -// throwIfNaN(J, "J"); -// throwIfLessThan(L, 3u, "L"); -// return TfiModel(L, J, g); -// } - - -enum class TEBD_TYPE { - BRICKWALL = 0, - SECOND_ORDER = 1, -}; -// EMSCRIPTEN_BINDINGS(tebd_module) { -// enum_(TEBD_BRICKWALL) -// .value("TEBD_BRICKWALL", TEBD_BRICKWALL) -// ; -// class>("MPS") -// .constructor() -// .function("getL", &MPS::getL) -// .function("getN", &MPS::getN) -// .function("getN_real", -// }; - -/** - * @brief Run TEBD on a model and measure energy, magnetization and half-chain entropy after each TEBD iteration - * @param callback: A javascript function taking 3 parameters, M, E and S after each iteration - */ -template -void runTEBD(MPS& psi, const std::shared_ptr model, unsigned nSteps, double dtReal, double dtImag, double _eps, unsigned chiMax, TEBD_TYPE type, Results& result) { - dtype eps(_eps); - dtype dt(dtReal, dtImag); - // Sanity checks - throwIfNaN(eps, "eps"); - throwIfNaN(dt, "dt"); - throwIfLessThan(nSteps, 1u, "nSteps"); - throwIfLessThan(chiMax, 1u, "chiMax"); - if (psi.getL() != model->L) { - throw std::runtime_error("psi and model have different L"); - } - std::println("Running TEBD type={}, nSteps={}, dt={}, eps={}, chiMax={}", static_cast(type), nSteps, dt, eps, chiMax); - - // collect M over time in an array - auto F = result.getPostUpdateFunction(model); - // initial values now - F(psi, 0); - - // std::println("Bonds:"); - const auto& bondsMat = model->getH_BondsMatrices(); - // std::cout << "First: \n" << *(bondsMat.front()) << std::endl; - // std::cout << "Secnd: \n" << *(bondsMat.at(1)) << std::endl; - // std::cout << "Last : \n" << *(bondsMat.back()) << std::endl; - - // const Bonds bonds = model.getH_Bonds(); - const Bonds bonds = getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt, model->localDim); - // std::println("U:"); - // int i = 0; - // for (const auto& bond : bonds) { - // printTensor(*bond, std::format("Bond[{:02}]", i)); - // i++; - // } - if (type == TEBD_TYPE::BRICKWALL) { - runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, F); - } else if (type == TEBD_TYPE::SECOND_ORDER) { - const Bonds bondsHalf = getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt/2., model->localDim); - runTEBD_secondOrder(psi, bonds, bondsHalf, chiMax, eps, nSteps, F); - } else { - throw std::invalid_argument(std::format("Got invalid TEBD_TYPE: {}", static_cast(type))); - } - std::println("Chis: {}", psi.getChi()); - std::println("Collapse: {}", psi.collapse(0)); -}; - - - -// @TODO -template -void runTEBD_phaseDiagram(MPS& psiI, unsigned nSteps, double dtReal, double dtImag, double _eps, unsigned chiMax, TEBD_TYPE type, unsigned nValues, Results& result, ModelArgs&&... args) { - dtype eps(_eps); - dtype dt(dtReal, dtImag); - size_t L = psiI.getL(); - // Sanity checks - throwIfNaN(eps, "eps"); - throwIfNaN(dt, "dt"); - throwIfLessThan(nSteps, 2u, "nSteps"); - throwIfLessThan(chiMax, 1u, "chiMax"); - std::println("Running TEBD type={}, nSteps={}, dt={}, eps={}, chiMax={}", static_cast(type), nSteps, dt, eps, chiMax); - - postUpdateFunction_t F = [](const MPS&, unsigned) { - emscripten_sleep(0); - }; - for (unsigned i = 0; i < nValues; i++) { - // TfiModel model(L, U, g); - std::shared_ptr model = makeModel(std::forward(args)...); - std::println("{:02} - Model: {}", i, *model); - auto psi = psiI; - const Bonds bonds = getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt, model->localDim); - runTEBD_brickwall(psi, bonds, chiMax, eps, nSteps, F); - result.postUpdateFunction(model, psi); - } -} - - -void test(unsigned x) { - Range range(0, 1.0, 3, 2); - - std::println("test({})", x); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - std::println("range({})", range.getNext()); - - std::println("{}", processType(0.5)); - std::println("{}", processType(range)); - - // processDebug(range, 4.3, -3, range); - // process2(range); -} - - -EMSCRIPTEN_BINDINGS(run) { - // see https://emscripten.org/docs/porting/connecting_cpp_and_javascript/embind.html#advanced-class-concepts - // Results - ems::class_>("ResultsBaseTfi"); - ems::class_>("ResultsBaseBh"); - ems::class_>>("ResultsTfi") - .constructor() - ; - ems::class_>>("ResultsBh") - .constructor() - ; - // Models - ems::class_>("Model1D") - .smart_ptr>>("Model1D") - ; - ems::class_>>("ModelTfi") - // .constructor() - .smart_ptr_constructor("ModelTfi", &std::make_shared) - ; - ems::class_>>("ModelBh") - // .constructor() - .smart_ptr_constructor("ModelBh", &std::make_shared) - ; - ems::class_>("DoubleRange") - .constructor() - ; - // Other - ems::class_>("MPS"); - ems::function("initMpsSpinUp", &initMPS_spinUp); - ems::function("initMpsSpinRight", &initMPS_spinRight); - ems::function("initMpsNParticles", &initMPS_nParticles); - // ems::function("initMPS_spinUp", &initMPS_spinUp); - ems::register_vector("stdVectorDouble"); - ems::register_vector>("stdVectorDouble2D"); - // ems::class_("TEBDResults") - // .property("M", &Results::M) - // .property("E", &Results::E) - // .property("S", &Results::S) - // ; - ems::enum_("TEBD_TYPE") - .value("BRICKWALL", TEBD_TYPE::BRICKWALL) - .value("SECOND_ORDER", TEBD_TYPE::SECOND_ORDER) - ; - ems::function("runTebdTfi", &runTEBD); - ems::function("runTebdPhaseDiagramTfi", &runTEBD_phaseDiagram, Range>); - ems::function("runTebdBh", &runTEBD); - ems::function("runTebdPhaseDiagramBh", &runTEBD_phaseDiagram, Range, unsigned>); - ems::function("test", &test); -}; - - - -#endif diff --git a/src/type.hpp b/src/type.hpp index 8227b87..f95c0ca 100644 --- a/src/type.hpp +++ b/src/type.hpp @@ -2,3 +2,5 @@ /// the datatype for all template instantiations using dtype = std::complex; +// template +// using Tensor = std::shared_ptr>; diff --git a/src/wasm/emscripten.hpp b/src/wasm/emscripten.hpp new file mode 100644 index 0000000..641b7a0 --- /dev/null +++ b/src/wasm/emscripten.hpp @@ -0,0 +1,27 @@ +#pragma once + +#ifdef EMSCRIPTEN +#include +#include +namespace ems = emscripten; +// #warning emscripten +#else +// #warning noEmscripten +#include +/** + * @brief Allow using ems::val without emscripten + */ +namespace ems { + struct val { + template + void operator()(Args...) { + std::println("Callback called"); + } + }; + +} +inline void emscripten_sleep(int time) { + std::println("emscripten_sleep({})", time); +} +#endif + diff --git a/src/wasm/range.hpp b/src/wasm/range.hpp new file mode 100644 index 0000000..cc3e817 --- /dev/null +++ b/src/wasm/range.hpp @@ -0,0 +1,88 @@ +#pragma once +#include +#include + +/** + * @brief A class for generating a range of values. + * @details + * Returns `nSteps` values with `start` being the first and `stop` being the last value. + * The same value is returned `countsPerStep` times before the next value is returned. + * After the last `stop` value has been returned, the cycle is restarted from the first `start` value. + */ +template +struct Range { + using value_type = T; + Range(T start, T stop, unsigned nSteps, unsigned countsPerStep=1) : start(start), stop(stop), nSteps(nSteps), stepCounter(0), countsPerStep(countsPerStep), counter(0) { + if (nSteps <= 1) { + throw std::invalid_argument("nSteps must be greater than 1"); + } + step = (stop - start) / (nSteps - 1); + }; + T getNext() { + T result = start + stepCounter * step; + if (++counter >= countsPerStep) { + counter = 0; + if (++stepCounter >= nSteps) { + stepCounter = 0; + } + } + return result; + }; + private: + T start; + T stop; + T step; + unsigned nSteps; + unsigned stepCounter; + unsigned countsPerStep; + unsigned counter; +}; + +/** + * This ugly mechanism lets the javascript pass + * any combination of Range and basic types to a function. + * Every time the function is called, the Range classes will return a new value while the rest stays the same. + */ +template class U> +inline constexpr bool is_instance_of_v = std::false_type{}; +template class U, class V> +inline constexpr bool is_instance_of_v,U> = std::true_type{}; + +template +concept isRange = is_instance_of_v, Range>; + + +// All types except Range +template +struct transform_type { + using type = T; +}; +// Range +template +struct transform_type>> { + using type = typename std::remove_reference_t::value_type; // If T is a Range, return its value_type +}; + + +static_assert(isRange>); +static_assert(isRange&>); +static_assert(!isRange); + +// Dispatcher +template +std::remove_reference_t::value_type processType(T&& obj) { + return obj.getNext(); +} +template +requires(!isRange) +T processType(T obj) { + return obj; +} + +/// Make a shared model pointer from a arguments, which may contain Range objects +template +std::shared_ptr makeModel(Args&&... args) { + std::tuple::type...> results = std::make_tuple(processType(std::forward(args))...); + return std::apply([](auto&&... args) { return std::make_shared(args...); }, results); +} + diff --git a/src/wasm/result.hpp b/src/wasm/result.hpp new file mode 100644 index 0000000..1faf6cb --- /dev/null +++ b/src/wasm/result.hpp @@ -0,0 +1,81 @@ +#pragma once +#include "algorithms/tebd.hpp" +#include "mps.hpp" +#include "models/transverse_field_ising.hpp" +#include "models/bose_hubbard.hpp" +#include "emscripten.hpp" +#include +#include +#include +#include +#include + +using TfiModel = TransverseFieldIsingModel1D; +using BhModel = BoseHubbardModel1D; + +template +class Results { + public: + virtual ~Results() {}; + virtual postUpdateFunction_t getPostUpdateFunction(std::shared_ptr model) = 0; + virtual bool postUpdateFunction(std::shared_ptr model, const MPS& psi) = 0; +}; + + +class TfiResults : public Results { + public: + using Model_t = TfiModel; + TfiResults(ems::val callback, double convergence) : callback(callback), convergence(convergence), previousE(std::numeric_limits::max()) { + sigmaZ.setValues({{1, 0}, {0, -1}}); + }; + /// @todo write a second postUpdateFunction without the JS callback and return that one if callback is null + postUpdateFunction_t getPostUpdateFunction(std::shared_ptr model) override { + return [this, model](const MPS& psi, unsigned i){ return this->postUpdateFunction(model, psi); }; + } + bool postUpdateFunction(std::shared_ptr model, const MPS& psi) override { + auto Mval = psi.getSiteExpectationValueSum(this->sigmaZ).real() / psi.getL(); + auto Eval = model->energy(psi).real(); + auto Sval = psi.getEntanglementEntropy(psi.getL()/2).real(); + auto chiVal = psi.getChiMax(); + callback(Mval, Eval, Sval, chiVal); + // hand control back to browser, otherwise site hangs and callback will not be executed + emscripten_sleep(0); + if (std::abs(Eval - previousE) < convergence) return true; + previousE = Eval; + return false; + } + private: + e::TensorFixedSize> sigmaZ; + ems::val callback; + double convergence; + double previousE; +}; + + +/** + * @brief Measures , E and S + */ +class BoseHubbardResults : public Results { + public: + BoseHubbardResults(ems::val callback, double convergence) : callback(callback), convergence(convergence), previousE(std::numeric_limits::max()){}; + postUpdateFunction_t getPostUpdateFunction(std::shared_ptr model) override { + return [this, model](const MPS& psi, unsigned) { return this->postUpdateFunction(model, psi); }; + } + bool postUpdateFunction(std::shared_ptr model, const MPS& psi) override { + e::TensorMap> particleNumberOperator(model->getLadderOperators().getN().data(), model->localDim, model->localDim); + auto Nval = psi.getSiteExpectationValueSum(particleNumberOperator).real() / model->L; + auto Eval = model->energy(psi).real(); + auto Sval = psi.getEntanglementEntropy(model->L/2).real(); + auto chiVal = psi.getChiMax(); + callback(Nval, Eval, Sval, chiVal); + // hand control back to browser, otherwise site hangs and callback will not be executed + emscripten_sleep(0); + if (std::abs(Eval - previousE) < convergence) return true; + previousE = Eval; + return false; + } + private: + ems::val callback; + double convergence; + double previousE; +}; diff --git a/src/wasm/run.cpp b/src/wasm/run.cpp new file mode 100644 index 0000000..06251e0 --- /dev/null +++ b/src/wasm/run.cpp @@ -0,0 +1,107 @@ +#include "wasm/run.hpp" + +#ifndef BUILD_OPTIONS_STRING +#define BUILD_OPTIONS_STRING "unknown" +#endif + +std::string getBuildOptions() { + return BUILD_OPTIONS_STRING; +} + +void test(unsigned x) { + Range range(0, 1.0, 3, 2); + + std::println("test({})", x); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + std::println("range({})", range.getNext()); + + std::println("{}", processType(0.5)); + std::println("{}", processType(range)); + + // processDebug(range, 4.3, -3, range); + // process2(range); +} + +/** + * @brief Function used to test performance in the browser + */ +int testPerformance() { + e::Tensor t4(4, 16, 16, 4); + auto data = t4.data(); + for (unsigned i = 0; i < t4.size(); i++) { + data[i] = i; + } + auto [A, S, B] = split_truncate_theta(t4, 20, 1e-10); + e::Tensor t4reconstructedButNormalized = A.contract(diagonal(S).contract(B, e::array{iP{1, 0}}), e::array{iP{2, 0}}); + return t4reconstructedButNormalized.size(); +} + +#ifdef EMSCRIPTEN +EMSCRIPTEN_BINDINGS(run) { + // see https://emscripten.org/docs/porting/connecting_cpp_and_javascript/embind.html#advanced-class-concepts + // RESULTS + ems::class_>("ResultsBaseTfi"); + ems::class_>("ResultsBaseBh"); + ems::class_>>("ResultsTfi") + .constructor() + ; + ems::class_>>("ResultsBh") + .constructor() + ; + // MODELS + ems::class_>("Model1D") + .smart_ptr>>("Model1D") + ; + ems::class_>>("ModelTfi") + // .constructor() + .smart_ptr_constructor("ModelTfi", &std::make_shared) + ; + ems::class_>>("ModelBh") + // .constructor() + .smart_ptr_constructor("ModelBh", &std::make_shared) + ; + ems::class_>("DoubleRange") + .constructor() + ; + // MPS + ems::class_>("MPS"); + ems::function("initMpsSpinUp", &initMPS_spinUp); + ems::function("initMpsSpinRight", &initMPS_spinRight); + ems::function("initMpsNParticles", &initMPS_nParticles); + // ems::function("initMPS_spinUp", &initMPS_spinUp); + ems::register_vector("stdVectorDouble"); + ems::register_vector>("stdVectorDouble2D"); + // ems::class_("TEBDResults") + // .property("M", &Results::M) + // .property("E", &Results::E) + // .property("S", &Results::S) + // ; + // TEBD + ems::enum_("TEBD_TYPE") + .value("BRICKWALL", TEBD_TYPE::BRICKWALL) + .value("SECOND_ORDER", TEBD_TYPE::SECOND_ORDER) + ; + ems::function("runTebdTfi", &runTEBD); + ems::function("runTebdPhaseDiagramTfi", &runTEBD_phaseDiagram, Range>); + ems::function("runTebdBh", &runTEBD); + ems::function("runTebdPhaseDiagramBh", &runTEBD_phaseDiagram, Range, unsigned>); + // DMRG + ems::function("runDmrgTfi", &runDMRG); + ems::function("runDmrgPhaseDiagramTfi", &runDMRG_phaseDiagram, Range>); + ems::function("runDmrgBh", &runDMRG); + ems::function("runDmrgPhaseDiagramBh", &runDMRG_phaseDiagram, Range, unsigned>); + // OTHER + ems::function("test", &test); + ems::function("getBuildOptions", &getBuildOptions); + ems::function("testPerformance", &testPerformance); +}; + +#endif diff --git a/src/wasm/run.hpp b/src/wasm/run.hpp new file mode 100644 index 0000000..f46fed2 --- /dev/null +++ b/src/wasm/run.hpp @@ -0,0 +1,166 @@ +#include +#include + +#include "wasm/result.hpp" +#include "algorithms/tebd.hpp" +#include "algorithms/dmrg.hpp" +#include "util.hpp" +#include "type.hpp" +#include "wasm/emscripten.hpp" +#include "wasm/range.hpp" + + +// template +// void printArgs(Arg&& val) { +// std::cout << val << std::endl; +// } +// template +// void printArgs(Arg&& val, Args&&... values) { +// std::cout << val << ", "; +// printArgs(std::forward(values)...); +// } + +// template +// void processDebug(Args&&... args) { +// std::tuple::type...> results = std::make_tuple(processType(std::forward(args))...); +// std::apply([](auto&&... args){ printArgs(args...); }, results); +// } + +/** + */ + +// EMSCRIPTEN_BINDINGS(my_class_example) { +// class_("MyClass") +// .constructor() +// .function("incrementX", &MyClass::incrementX) +// .property("x", &MyClass::getX, &MyClass::setX) +// .property("x_readonly", &MyClass::getX) +// .class_function("getStringFromInstance", &MyClass::getStringFromInstance) +// ; +// }; + +/** + * @brief Return a string that was set at compilation using the BUILD_OPTIONS_STRING variable + */ +std::string getBuildOptions(); + + +enum class TEBD_TYPE { + BRICKWALL = 0, + SECOND_ORDER = 1, +}; +// EMSCRIPTEN_BINDINGS(tebd_module) { +// enum_(TEBD_BRICKWALL) +// .value("TEBD_BRICKWALL", TEBD_BRICKWALL) +// ; +// class>("MPS") +// .constructor() +// .function("getL", &MPS::getL) +// .function("getN", &MPS::getN) +// .function("getN_real", +// }; + +// TEBD +template +void runTEBD(MPS& psi, const std::shared_ptr model, unsigned nSteps, double dtReal, double dtImag, double _eps, unsigned chiMax, TEBD_TYPE type, Results& result) { + dtype eps(_eps); + dtype dt(dtReal, dtImag); + // sanity checks + throwIfNaN(eps, "eps"); + throwIfNaN(dt, "dt"); + throwIfLessThan(nSteps, 1u, "nSteps"); + throwIfLessThan(chiMax, 1u, "chiMax"); + if (psi.getL() != model->L) { + throw std::runtime_error("psi and model have different L"); + } + std::println("Running TEBD type={}, nSteps={}, dt={}, eps={}, chiMax={}", static_cast(type), nSteps, dt, eps, chiMax); + + auto F = result.getPostUpdateFunction(model); + // initial values now + F(psi, 0); + + const Bonds bonds = tebd::getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt, model->localDim); + if (type == TEBD_TYPE::BRICKWALL) { + tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, F); + } else if (type == TEBD_TYPE::SECOND_ORDER) { + const Bonds bondsHalf = tebd::getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt/2., model->localDim); + tebd::runSecondOrder(psi, bonds, bondsHalf, chiMax, eps, nSteps, F); + } else { + throw std::invalid_argument(std::format("Got invalid TEBD_TYPE: {}", static_cast(type))); + } + std::println("Chis: {}", psi.getChi()); + std::println("Collapse: {}", psi.collapse(0)); +}; + + +template +void runTEBD_phaseDiagram(const MPS& psiI, unsigned nSteps, double dtReal, double dtImag, double _eps, unsigned chiMax, TEBD_TYPE type, unsigned nValues, Results& result, ModelArgs&&... args) { + dtype eps(_eps); + dtype dt(dtReal, dtImag); + // sanity checks + throwIfNaN(eps, "eps"); + throwIfNaN(dt, "dt"); + throwIfLessThan(nSteps, 2u, "nSteps"); + throwIfLessThan(chiMax, 1u, "chiMax"); + std::println("Running TEBD type={}, nSteps={}, dt={}, eps={}, chiMax={}", static_cast(type), nSteps, dt, eps, chiMax); + + postUpdateFunction_t F = [](const MPS&, unsigned) { + emscripten_sleep(0); + return false; + }; + for (unsigned i = 0; i < nValues; i++) { + std::shared_ptr model = makeModel(std::forward(args)...); + std::println("{:02} - Model: {}", i, *model); + MPS psi(psiI); + std::println("New psi: value at B[2](0, 1, 0)={}", psi.B(2)(0, 1, 0)); + const Bonds bonds = tebd::getTimeEvolutionMatrices(model->getH_BondsMatrices(), dt, model->localDim); + tebd::runBrickwall(psi, bonds, chiMax, eps, nSteps, F); + result.postUpdateFunction(model, psi); + } +} + + +// DMRG +template +void runDMRG(MPS& psi, const std::shared_ptr model, unsigned nSweeps, double _eps, unsigned chiMax, Results& result) { + dtype eps(_eps); + // sanity checks + throwIfNaN(eps, "eps"); + throwIfLessThan(nSweeps, 1u, "nSweeps"); + throwIfLessThan(chiMax, 1u, "chiMax"); + if (psi.getL() != model->L) { + throw std::runtime_error("psi and model have different L"); + } + std::println("Running DMRG nSweeps={}, eps={}, chiMax={}", nSweeps, eps, chiMax); + + auto F = result.getPostUpdateFunction(model); + // initial values now + F(psi, 0); + + dmrg::run(psi, model->getH_MPO(), chiMax, eps, nSweeps, F); + std::println("Chis: {}", psi.getChi()); + std::println("Collapse: {}", psi.collapse(0)); +}; + + +template +void runDMRG_phaseDiagram(const MPS& psiI, unsigned nSweeps, double _eps, unsigned chiMax, unsigned nValues, Results& result, ModelArgs&&... args) { + dtype eps(_eps); + // sanity checks + throwIfNaN(eps, "eps"); + throwIfLessThan(nSweeps, 1u, "nSweeps"); + throwIfLessThan(chiMax, 1u, "chiMax"); + std::println("Running DMRG nSweeps={}, eps={}, chiMax={}", nSweeps, eps, chiMax); + + postUpdateFunction_t F = [](const MPS&, unsigned) { + emscripten_sleep(0); + return false; + }; + for (unsigned i = 0; i < nValues; i++) { + std::shared_ptr model = makeModel(std::forward(args)...); + std::println("{:02} - Model: {}", i, *model); + MPS psi(psiI); + dmrg::run(psi, model->getH_MPO(), chiMax, eps, nSweeps, F); + result.postUpdateFunction(model, psi); + } +} diff --git a/src/wasm/util.hpp b/src/wasm/util.hpp new file mode 100644 index 0000000..e69de29