From b617edb95d3c327ddac361d4de4c505e928e0086 Mon Sep 17 00:00:00 2001 From: Alexander Wietek Date: Sat, 3 Aug 2024 13:45:37 +0200 Subject: [PATCH] Documented State class, fixed bug getting imaginary part of State. --- .gitignore | 1 + docs/documentation/blocks/spinhalf.md | 16 +- docs/documentation/index.md | 6 + docs/documentation/states/state.md | 217 ++++++++++++++++++ examples/usage_examples/main.cpp | 25 +- examples/usage_examples/main.jl | 18 ++ .../time_evolution/test_time_evolution.cpp | 4 +- .../test_time_evolution_distributed.cpp | 4 +- xdiag/states/fill.cpp | 2 +- xdiag/states/product_state.cpp | 21 +- xdiag/states/product_state.hpp | 5 +- xdiag/states/state.cpp | 9 +- 12 files changed, 305 insertions(+), 23 deletions(-) create mode 100644 docs/documentation/states/state.md diff --git a/.gitignore b/.gitignore index 2c7444cc..853c1708 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +cmake_commands.sh build install julia/products diff --git a/docs/documentation/blocks/spinhalf.md b/docs/documentation/blocks/spinhalf.md index 38e25799..23a26f4a 100644 --- a/docs/documentation/blocks/spinhalf.md +++ b/docs/documentation/blocks/spinhalf.md @@ -52,8 +52,22 @@ Representation of a block in a spin $S=1/2$ Hilbert space. int64_t n_sites() const ``` +!!! method "dim" + Returns the dimension of the block. + + === "Julia" + ```julia + dim(block::Spinhalf) + ``` + + === "C++" + ```c++ + int64_t dim() const; + ``` + + !!! method "size" - Returns the size of the block, i.e. its dimension. + Returns the size of the block locally. Same as "dim" for non-distributed Blocks but different for distributed blocks. === "Julia" ```julia diff --git a/docs/documentation/index.md b/docs/documentation/index.md index e4a50e67..32a90488 100644 --- a/docs/documentation/index.md +++ b/docs/documentation/index.md @@ -32,6 +32,12 @@ title: Overview | [Coupling](operators/coupling.md) | Describes the coupling of a local operator | :simple-cplusplus: :simple-julia: | +## States +| | | | +|:-------------------------|:---------------------------------------------------|----------------------------------:| +| [State](states/state.md) | A generic state describing a quantum wave function | :simple-cplusplus: :simple-julia: | + + ## Symmetries | | | | |:----------------------------------------------------|:----------------------------------------------------|----------------------------------:| diff --git a/docs/documentation/states/state.md b/docs/documentation/states/state.md new file mode 100644 index 00000000..4bc8328a --- /dev/null +++ b/docs/documentation/states/state.md @@ -0,0 +1,217 @@ +--- +title: State +--- + +A generic state describing a quantum wave function + +**Source** [state.hpp](https://github.com/awietek/xdiag/blob/main/xdiag/states/state.hpp) + +## Constructors + +=== "Julia" + ```julia + State(block::Block; real::Bool = true, n_cols::Int64 = 1) + State(block::Block, vec::Vector{Float64}) + State(block::Block, vec::Vector{ComplexF64}) + State(block::Block, mat::Matrix{Float64}) + State(block::Block, mat::Matrix{ComplexF64}) + ``` + +=== "C++" + ```c++ + State(Block const &block, bool real = true, int64_t n_cols = 1); + + template + State(block_t const &block, arma::Col const &vector); + + template + State(block_t const &block, arma::Mat const &matrix); + ``` + +| Parameter | Description | | +|:----------|:-----------------------------------------------------------------------------------------------|---| +| block | The block of a Hilbertspace on which the state is defined | | +| real | Flag whether or not the state has real coefficients | | +| n_cols | Number of columns of the state (default 1) | | +| vector | A vector containing the coefficients of the state. Must be same size as block. | | +| matrix | A matrix containing the coefficients of the state. Number of rows must be same as block size . | | + + +## Methods + +!!! method "n_sites" + + Returns the number of sites of the block the state is defined on. + + === "Julia" + ```julia + n_sites(state::State) + ``` + + === "C++" + ```c++ + int64_t n_sites() const + ``` + +!!! method "isreal" + Returns whether the state is real. + + === "Julia" + ```julia + isreal(state::State) + ``` + + === "C++" + ```c++ + int64_t isreal() const; + ``` + +!!! method "real" + Returns whether the real part of the State. + + === "Julia" + ```julia + real(state::State) + ``` + + === "C++" + ```c++ + State real() const; + ``` + +!!! method "imag" + Returns whether the imaginary part of the State. + + === "Julia" + ```julia + imag(state::State) + ``` + + === "C++" + ```c++ + State imag() const; + ``` + +!!! method "make_complex! / make_complex" + Turns a real State into a complex State. Does nothing if the state is already complex + + === "Julia" + ```julia + make_complex!(state::State) + ``` + + === "C++" + ```c++ + void make_complex(); + ``` + + +!!! method "dim" + Returns the dimension of the block the state is defined on. + + === "Julia" + ```julia + dim(block::Spinhalf) + ``` + + === "C++" + ```c++ + int64_t dim() const; + ``` + + +!!! method "size" + Returns the size of the block the state is defined on. locally. Same as "dim" for non-distributed Blocks but different for distributed blocks. + + === "Julia" + ```julia + size(block::Spinhalf) + ``` + + === "C++" + ```c++ + int64_t size() const; + ``` + +!!! method "n_rows" + Returns number of rows of the local storage. Same as "size" + + === "Julia" + ```julia + n_rows(block::Spinhalf) + ``` + + === "C++" + ```c++ + int64_t n_rows() const; + ``` + +!!! method "n_cols" + Returns number of columns of the local storage. + + === "Julia" + ```julia + n_cols(block::Spinhalf) + ``` + + === "C++" + ```c++ + int64_t n_cols() const; + ``` + +!!! method "col" + Returns a state created from the n-th column of the storage. Whether or not the storage is copied can be specified by setting the flag "copy". + + === "Julia" + ```julia + col(state::State, n::Int64 = 1; copy::Bool = true) + ``` + + === "C++" + ```c++ + State col(int64_t n, bool copy = true) const; + ``` + +!!! method "vector/vectorC" + Returns a vector from the n-th column of the storage. In C++ use "vector"/"vectorC" to either get a real or complex vector. + + === "Julia" + ```julia + vector(state::State; n::Int64 = 1) + # no vectorC method in julia + ``` + + === "C++" + ```c++ + arma::vec vector(int64_t n = 0, bool copy = true) const; + arma::cx_vec vectorC(int64_t n = 0, bool copy = true) const; + ``` + +!!! method "matrix/matrixC" + Returns matrix representing the storage. In C++ use "matrix"/"matrixC" to either get a real or complex matrix. + + === "Julia" + ```julia + matrix(state::State) + # no matrixC method in julia + ``` + + === "C++" + ```c++ + arma::vec matrix(bool copy = true) const; + arma::cx_vec matrixC(bool copy = true) const; + ``` + + +## Usage Example + +=== "Julia" + ```c++ + --8<-- "examples/usage_examples/main.jl:state" + ``` + +=== "C++" + ```c++ + --8<-- "examples/usage_examples/main.cpp:state" + ``` + diff --git a/examples/usage_examples/main.cpp b/examples/usage_examples/main.cpp index 351f0987..c16cd23a 100644 --- a/examples/usage_examples/main.cpp +++ b/examples/usage_examples/main.cpp @@ -260,7 +260,30 @@ XDIAG_SHOW(cpl.as()); // --8<-- [end:coupling] } -// clang-format on + + +{ +// --8<-- [start:state] +auto block = Spinhalf(2); +auto psi1 = State(block, arma::vec("1.0 2.0 3.0 4.0")); +XDIAG_SHOW(psi1); +XDIAG_SHOW(psi1.vector()); +psi1.make_complex(); +XDIAG_SHOW(psi1.vectorC()); + +auto psi2 = State(block, false, 3); +XDIAG_SHOW(psi2); +XDIAG_SHOW(psi2.matrixC()); + +auto psi3 = State(block, arma::cx_vec(arma::vec("1.0 2.0 3.0 4.0"), + arma::vec("4.0 3.0 2.0 1.0"))); +XDIAG_SHOW(psi3.vectorC()); +XDIAG_SHOW(psi3.real().vector()); +XDIAG_SHOW(psi3.imag().vector()); +// --8<-- [end:state] +} + + // clang-format on } catch (Error e) { error_trace(e); diff --git a/examples/usage_examples/main.jl b/examples/usage_examples/main.jl index 8752905d..95836228 100644 --- a/examples/usage_examples/main.jl +++ b/examples/usage_examples/main.jl @@ -240,3 +240,21 @@ cpl = Coupling([1 2; -2 1]) @show convert(Matrix{Float64}, cpl) @show convert(Matrix{ComplexF64}, cpl) # --8<-- [end:coupling] + +# --8<-- [start:state] +block = Spinhalf(2) +psi1 = State(block, [1.0, 2.0, 3.0, 4.0]) +@show psi1 +display(vector(psi1)) +make_complex!(psi1) +display(vector(psi1)) + +psi2 = State(block, real=false, n_cols=3) +@show psi2 +display(matrix(psi2)) + +psi3 = State(block, [1.0+4.0im, 2.0+3.0im, 3.0+2.0im, 4.0+1.0im]) +display(vector(psi3)) +display(vector(real(psi3))) +display(vector(imag(psi3))) +# --8<-- [end:state] diff --git a/tests/algorithms/time_evolution/test_time_evolution.cpp b/tests/algorithms/time_evolution/test_time_evolution.cpp index 114a7d22..27d43476 100644 --- a/tests/algorithms/time_evolution/test_time_evolution.cpp +++ b/tests/algorithms/time_evolution/test_time_evolution.cpp @@ -221,9 +221,9 @@ TEST_CASE("tj_complex_timeevo", "[time_evolution]") try { for (int x = 0; x < L; ++x) { for (int y = 0; y < L; ++y) { if (((x + y) % 2) == 0) { - pstate << "Dn"; + pstate.push_back("Dn"); } else { - pstate << "Up"; + pstate.push_back("Up"); } } } diff --git a/tests/algorithms/time_evolution/test_time_evolution_distributed.cpp b/tests/algorithms/time_evolution/test_time_evolution_distributed.cpp index c1a88613..59793795 100644 --- a/tests/algorithms/time_evolution/test_time_evolution_distributed.cpp +++ b/tests/algorithms/time_evolution/test_time_evolution_distributed.cpp @@ -47,9 +47,9 @@ TEST_CASE("time_evolution_distributed", "[time_evolution]") try { for (int x = 0; x < L; ++x) { for (int y = 0; y < L; ++y) { if (((x + y) % 2) == 0) { - pstate << "Dn"; + pstate.push_back("Dn"); } else { - pstate << "Up"; + pstate.push_back("Up"); } } } diff --git a/xdiag/states/fill.cpp b/xdiag/states/fill.cpp index 10dba175..2c81d466 100644 --- a/xdiag/states/fill.cpp +++ b/xdiag/states/fill.cpp @@ -34,7 +34,7 @@ void fill(State &state, RandomState const &rstate, int64_t col) try { } void fill(State &state, ProductState const &pstate, int64_t col) try { - if (state.n_sites() != pstate.n_sites()) { + if (state.size() != pstate.size()) { XDIAG_THROW("State and ProductState do not have the same number of sites"); } else if (col >= state.n_cols()) { XDIAG_THROW("Column index larger than number of columns in State"); diff --git a/xdiag/states/product_state.cpp b/xdiag/states/product_state.cpp index c38d253f..f510e8cb 100644 --- a/xdiag/states/product_state.cpp +++ b/xdiag/states/product_state.cpp @@ -14,8 +14,9 @@ std::string const &ProductState::operator[](int64_t i) const { } std::string &ProductState::operator[](int64_t i) { return local_states_[i]; } -int64_t ProductState::n_sites() const { return local_states_.size(); } -void ProductState::operator<<(std::string l) { local_states_.push_back(l); } +void ProductState::push_back(std::string l) { local_states_.push_back(l); } +int64_t ProductState::size() const { return local_states_.size(); } + ProductState::iterator_t ProductState::begin() const { return local_states_.begin(); } @@ -31,7 +32,7 @@ bool ProductState::operator!=(ProductState const &rhs) const { } std::ostream &operator<<(std::ostream &out, ProductState const &state) { - for (int64_t i = state.n_sites() - 1; i >= 0; --i) { + for (int64_t i = state.size() - 1; i >= 0; --i) { out << state[i] << " "; } out << "\n"; @@ -42,7 +43,7 @@ std::string to_string(ProductState const &state, std::string format) try { return to_string_generic(state); } else if (format == "fancy") { std::stringstream ss; - for (int64_t i = state.n_sites() - 1; i >= 0; --i) { + for (int64_t i = state.size() - 1; i >= 0; --i) { std::string s = state[i]; if (s == "Up") { // const char *s = u8"\u2B61"; @@ -76,7 +77,7 @@ std::string to_string(ProductState const &state, std::string format) try { template void to_product_state_spinhalf(bit_t spins, ProductState &pstate) { - for (int64_t i = 0; i < pstate.n_sites(); ++i) { + for (int64_t i = 0; i < pstate.size(); ++i) { if (spins & 1) { pstate[i] = "Up"; } else { @@ -90,7 +91,7 @@ template void to_product_state_spinhalf(uint64_t, ProductState &); template void to_product_state_tj(bit_t ups, bit_t dns, ProductState &pstate) { - for (int64_t i = 0; i < pstate.n_sites(); ++i) { + for (int64_t i = 0; i < pstate.size(); ++i) { if (ups & 1) { pstate[i] = "Up"; } else { @@ -111,7 +112,7 @@ template void to_product_state_tj(uint64_t ups, uint64_t dns, template void to_product_state_electron(bit_t ups, bit_t dns, ProductState &pstate) { - for (int64_t i = 0; i < pstate.n_sites(); ++i) { + for (int64_t i = 0; i < pstate.size(); ++i) { if (ups & 1) { if (dns & 1) { pstate[i] = "UpDn"; @@ -137,7 +138,7 @@ template void to_product_state_electron(uint64_t ups, uint64_t dns, template bit_t to_bits_spinhalf(ProductState const &pstate) try { bit_t spins = 0; - for (int64_t s = 0; s < pstate.n_sites(); ++s) { + for (int64_t s = 0; s < pstate.size(); ++s) { std::string val = pstate[s]; if (val == "Up") { spins |= ((bit_t)1 << s); @@ -159,7 +160,7 @@ template std::pair to_bits_tj(ProductState const &pstate) try { bit_t ups = 0; bit_t dns = 0; - for (int64_t s = 0; s < pstate.n_sites(); ++s) { + for (int64_t s = 0; s < pstate.size(); ++s) { std::string val = pstate[s]; if (val == "Up") { ups |= ((bit_t)1 << s); @@ -182,7 +183,7 @@ template std::pair to_bits_electron(ProductState const &pstate) try { bit_t ups = 0; bit_t dns = 0; - for (int64_t s = 0; s < pstate.n_sites(); ++s) { + for (int64_t s = 0; s < pstate.size(); ++s) { std::string val = pstate[s]; if (val == "Up") { ups |= ((bit_t)1 << s); diff --git a/xdiag/states/product_state.hpp b/xdiag/states/product_state.hpp index a9d8d5f7..d1ace2d6 100644 --- a/xdiag/states/product_state.hpp +++ b/xdiag/states/product_state.hpp @@ -17,8 +17,9 @@ class ProductState { std::string const &operator[](int64_t i) const; std::string &operator[](int64_t i); - int64_t n_sites() const; - void operator<<(std::string l); + int64_t size() const; + void push_back(std::string l); + iterator_t begin() const; iterator_t end() const; diff --git a/xdiag/states/state.cpp b/xdiag/states/state.cpp index 35b65ab1..8196614b 100644 --- a/xdiag/states/state.cpp +++ b/xdiag/states/state.cpp @@ -145,7 +145,6 @@ State State::real() const try { if (isreal()) { return (*this); } else { - // Return a State with zero values double *ptr = storage_.data(); return std::visit( [&](auto &&block) { return State(block, ptr, n_cols_, 2); }, block_); @@ -157,9 +156,12 @@ State State::real() const try { State State::imag() const try { // TODO: DOES THIS DO WHAT IT"S SUPPOSED TO? if (isreal()) { - return (*this); - } else { return State(block_, true, n_cols_); + } else { + double *ptr = storage_.data(); + return std::visit( + [&](auto &&block) { return State(block, ptr + 1, n_cols_, 2); }, + block_); } } catch (Error const &e) { XDIAG_RETHROW(e); @@ -280,7 +282,6 @@ complex *State::colptrC(int64_t col) { return memptrC() + col * n_rows_; } - template State::State(Spinhalf const &, bool, int64_t); template State::State(tJ const &, bool, int64_t); template State::State(Electron const &, bool, int64_t);