.. _state_vectors:

State vectors
=============

.. _state_vector:

``StateVector`` concept
-----------------------

.. default-domain:: cpp

.. namespace:: libcommute

Let us say we want to make type ``SV`` a *libcommute*-compatible state
vector type so that :ref:`linear operators <loperator>` can act on instances
of ``SV``. For this we have to make ``SV`` model the ``StateVector``
concept.
In a nutshell, a ``StateVector`` type is a one-dimensional array of numbers
(quantum amplitudes) allowing integer indexing and implementing a certain
interface. Elements of the array do not have to be stored contiguously.
Acceptable index values must be at least 64-bit wide unsigned integers since
*libcommute* uses the following type for basis state indexing.

.. type:: sv_index_type = std::uint64_t

  *Defined in <libcommute/loperator/state_vector.hpp>*

  Type of basis state indices.

The table below shows the interface (a set of free functions and a
metafunction) that needs be implemented for an object ``sv`` of type
``SV``.

*libcommute* provides an implementation of the ``StateVector`` concept for
``std::vector`` (see *<libcommute/loperator/state_vector.hpp>*).

.. list-table::
  :header-rows: 1

  * - Function/metafunction
    - Description
    - Implementation for ``std::vector<T>``
  * - ``element_type<SV>::type``
    - Type of the elements.
    - ``T``

  * - ``get_element(sv, n)``
    - Return the ``n``-th element of ``sv``.
    - ``sv[n]``

  * - ``update_add_element(sv, n, value)``
    - Add a value of some type ``U`` to the ``n``-th element of
      ``sv``.
    - ``sv[n] += value`` or ``sv[n] = sv[n] + value``

      The compound-assignment from type ``U`` will be used
      whenever ``sv``'s elements support it. Otherwise, the implementation
      will fall back to the simple addition.

  * - ``set_zeros(sv)``
    - Fill ``sv`` with zeros.
    - ``std::fill(sv.begin(), sv.end(), zero)``.

      The zero value is created by ``make_const(0)`` as described in
      ":ref:`custom_scalar_type`".

  * - ``zeros_like(sv)``
    - Return an object of the same type and size as ``sv`` but filled with
      zeros.
    - Creates a new object as ``std::vector<T>(sv.size(), zero)``.

  * - ``foreach(sv, f)``
    - Apply a function-like object ``f`` to all basis state index/non-zero
      element pairs ``(n, a)`` in ``sv``.
    - In a for-loop, calls ``f(n, a)`` for all non-zero elements ``a``
      as detected by ``is_zero()`` (see ":ref:`custom_scalar_type`").

Inclusion of *<libcommute/loperator/state_vector_eigen3.hpp>* makes some
`Eigen 3 <https://eigen.tuxfamily.org/>`_ types (`column vectors`_,
`vector blocks`_,
`column-like matrix blocks`_ and one-dimensional `Eigen::Map views`_)
compatible with the ``StateVector`` concept as well.

.. _column vectors:
  https://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html
  #TutorialMatrixVectors

.. _vector blocks:
  https://eigen.tuxfamily.org/dox/classEigen_1_1VectorBlock.html

.. _column-like matrix blocks:
  https://eigen.tuxfamily.org/dox/group__TutorialBlockOperations.html
  #TutorialBlockOperationsSyntaxColumnRows

.. _Eigen::Map views:
  https://eigen.tuxfamily.org/dox/classEigen_1_1Map.html


.. _sparse_state_vector:

Sparse state vector
-------------------

:class:`sparse_state_vector` is a state vector that saves memory by storing only
the non-zero elements. It is essentially a wrapper around
:class:`std::unordered_map` modelling the ``StateVector`` concept. Here, we show
only the part of its interface not covered by ``StateVector``.

.. class:: template<typename ScalarType> sparse_state_vector

  *Defined in <libcommute/loperator/sparse_state_vector.hpp>*

  State vector with a sparse storage of elements (quantum amplitudes).
  :type:`ScalarType` is the type of the elements.

  .. function::   sparse_state_vector() = delete
                  sparse_state_vector(sv_index_type size)

    Construct a zero (empty) sparse vector with a given :var:`size` --
    dimension of the corresponding Hilbert space.

  .. function:: sv_index_type size() const

    Size (dimension) of the vector.

  .. function:: ScalarType & operator[](sv_index_type n)

    Access the :var:`n`-th element. If it is zero (missing from the storage),
    then a new value-initialized element will be inserted and a reference to
    it will be returned.

    .. warning::

      Improper use of this method may result in zero elements being stored in
      the unordered map. Only the non-zero values should be assigned to the
      references returned by it.

  .. function:: sv_index_type n_nonzeros() const

    Get the number of non-zero (stored) elements.

  .. function:: void prune()

    Remove all zero elements (as defined by
    :ref:`scalar_traits\<ScalarType\>::is_zero() <custom_scalar_type>`)
    from the unordered map.

  .. function:: template<typename UnaryPredicate> void prune(UnaryPredicate&& p)

    Remove unordered map elements (amplitudes) for which predicate :var:`p`
    returns ``true``.

.. _compressed_state_view:

View of a compressed state vector
---------------------------------

Another way to save memory, in particular in the situation when the
:ref:`Hilbert space is sparse <sparse_hilbert_space>`, is by using
:class:`compressed_state_view`. A simple array-like container
(e.g. ``std::vector``) compatible with a sparse Hilbert space must have the
size :func:`hilbert_space::vec_size()`, which exceeds the actual dimension of
the space,
:func:`hilbert_space::dim() <sv_index_type hilbert_space::dim() const>`.
The extra required memory is the price to pay for improved performance. This
price, however, can be very steep for the bigger and sparser (i.e. including
many elementary spaces of non-power-of-two dimensions) Hilbert spaces. For
example, a vector storing a state of a system of :math:`N=10` spins
:math:`S=1` needs to be :math:`4^N` elements long, while only :math:`3^N` of its
elements represent amplitudes of the physical basis states. As a result, the
memory overhead in this case is huge, :math:`1 - (3/4)^N \approx 94\%`.

One can solve this issue by allocating an array of the size :math:`d =`
:func:`hilbert_space::dim() <sv_index_type hilbert_space::dim() const>`
(only the physical amplitudes are actually stored) and using a
:class:`compressed_state_view` object to adapt it.
:class:`compressed_state_view` models the ``StateVector`` concept and maps
indices of the physical basis states onto a continuous integer range
:math:`[0; d-1]` (this operation consumes a bit of extra time). The mapped
indices are then used to address elements of the underlying array.

.. class:: template<typename StateVector, bool Ref = true> compressed_state_view

  *Defined in <libcommute/loperator/compressed_state_view.hpp>*

  View of a :type:`StateVector` object that translates indices of basis states
  from a (possibly) :ref:`sparse Hilbert space <sparse_hilbert_space>` onto
  a continuous range. It is, therefore, sufficient to store only the
  :func:`hilbert_space::dim() <sv_index_type hilbert_space::dim() const>`
  physical amplitudes in the underlying state vector object, which can be
  much more memory-efficient than storing all :func:`hilbert_space::vec_size()`
  elements.

  :type:`StateVector` - type of the underlying state vector object. Defining a
  read-only view (such that prohibits :expr:`update_add_element()` operations)
  requires using a ``const``-qualified type here. For example, one can use
  ``StateVector = std::vector<double>`` for a read-write view, and
  ``StateVector = const std::vector<double>`` for a read-only view.

  :type:`Ref` - by default, :type:`compressed_state_view`
  stores a reference to the underlying state vector. Setting this option to
  ``false`` will result in a copy being created and stored instead. This feature
  can be useful when the underlying type is already a view-like object similar
  to ``Eigen::Map``.

  .. function:: template <typename SV, typename HSType> \
                         compressed_state_view(SV&& sv, HSType const& hs)

    Construct a view of the state vector :var:`sv`, defined in the
    (sparse) Hilbert space :var:`hs`.

  .. function:: sv_index_type map_index(sv_index_type index) const

    Translate a basis state :var:`index` from the Hilbert space to the
    continuous range ``[0; hs.dim()-1]``. If the Hilbert space is sparse,
    and :var:`index` does not correspond to a physical basis state, then
    the result is undefined.

*<libcommute/loperator/compressed_state_view.hpp>* defines two supplemental
factory functions for the :class:`compressed_state_view` objects.

.. function:: template <typename StateVector, typename HSType> \
              auto make_comp_state_view(StateVector&& sv, HSType const& hs)
              template <typename StateVector, typename HSType> \
              auto make_const_comp_state_view(StateVector&& sv, \
                                              HSType const& hs)

  Make and return a read/write or constant view of a compressed state vector
  :var:`sv` from the Hilbert space :var:`hs`.

.. _mapped_basis_view:

Mapped basis view
-----------------

:class:`mapped_basis_view` is another utility type modelling the ``StateVector``
concept. It is a view of a state vector, which translates basis state
index arguments of :func:`get_element()` and :func:`update_add_element()`
according to a predefined map :type:`sv_index_type` -> :type:`sv_index_type`.
The element access functions throw :type:`std::out_of_range` if their index
argument is missing from the map.

:class:`mapped_basis_view` can be used in situations
where a :ref:`linear operator <loperator>` acts in a small subspace of
a full Hilbert space, and it is desirable to store vector components only within
that subspace. Such a situation naturally emerges when working with
:ref:`invariant subspaces of operators <space_partition>`.

.. class:: template<typename StateVector, bool Ref = true> mapped_basis_view

  *Defined in <libcommute/loperator/mapped_basis_view.hpp>*

  View of a :type:`StateVector` object that translates basis state indices
  according to a certain mapping.

  :type:`StateVector` - type of the underlying state vector object. Defining a
  read-only view (such that prohibits :expr:`update_add_element()` operations)
  requires using a ``const``-qualified type here. For example, one can use
  ``StateVector = std::vector<double>`` for a read-write view, and
  ``StateVector = const std::vector<double>`` for a read-only view.

  .. _mapped_basis_view_Ref:

  :type:`Ref` - by default, :type:`mapped_basis_view`
  stores a reference to the underlying state vector. Setting this option to
  ``false`` will result in a copy being created and stored instead. This feature
  can be useful when the underlying type is already a view-like object similar
  to ``Eigen::Map``.

The mapped basis views should always be constructed by means of a special
factory class :class:`basis_mapper` and its methods
:func:`basis_mapper:: make_view()`/:func:`basis_mapper::make_const_view()`.

.. class:: basis_mapper

  *Defined in <libcommute/loperator/mapped_basis_view.hpp>*

  Factory class for :class:`mapped_basis_view`.

  .. rubric:: Constructors

  .. function:: basis_mapper(std::vector<sv_index_type> const& \
                             basis_state_indices)

    Build a mapping from a list of basis states :var:`basis_state_indices`
    to their positions within the list.

    .. code-block:: cpp

      std::vector<sv_index_type> basis_indices{3, 5, 6};
      basis_mapper mapper(basis_indices);

      // Views created by 'mapper' will translate basis state indices
      // according to
      // 0 -> std::out_of_range
      // 1 -> std::out_of_range
      // 2 -> std::out_of_range
      // 3 -> 0
      // 4 -> std::out_of_range
      // 5 -> 1
      // 6 -> 2
      // 7 -> std::out_of_range
      // ...

  .. function:: template<typename HSType, \
                         typename LOpScalarType, \
                         int... LOpAlgebraIDs> \
                basis_mapper(loperator<LOpScalarType,LOpAlgebraIDs...>const& O,\
                             HSType const& hs)

    Build a mapping from a set of all basis states contributing to
    :math:`\hat O|0\rangle`.

    Operator :var:`O` acts in the Hilbert space :var:`hs`.
    :math:`|0\rangle` is the basis state with index 0 ('vacuum' state in
    the case of fermions and bosons).
    Mapped values are assigned continuously starting from 0 without any specific
    order.

  .. function:: template<typename HSType, \
                         typename LOpScalarType, \
                         int... LOpAlgebraIDs> \
                basis_mapper( \
                std::vector<loperator<LOpScalarType, LOpAlgebraIDs...>> \
                  const& O_list, \
                HSType const& hs, unsigned int N)

    Given a list of operators
    :math:`\{\hat O_1, \hat O_2, \hat O_3, \ldots, \hat O_M\}`, build a mapping
    from all basis states contributing to all states
    :math:`\hat O_1^{n_1} \hat O_2^{n_2} \ldots \hat O_M^{n_M} |0\rangle`,
    where :math:`n_m \geq 0` and :math:`\sum_{m=1}^M n_M = N`.

    Operators in :var:`O_list` act in the Hilbert space :var:`hs`.
    :math:`|0\rangle` is the basis state with index 0 ('vacuum' state in
    the case of fermions and bosons).
    Mapped values are assigned continuously starting from 0 without any specific
    order.

    This constructor is useful to create a mapping from a fixed-particle-number
    subspace of a fermionic/bosonic Hilbert space.

  .. rubric:: :class:`mapped_basis_view` factory functions

  .. function:: template<typename StateVector> \
                mapped_basis_view<StateVector> \
                make_view(StateVector && sv) const
                template<typename StateVector> \
                mapped_basis_view<StateVector const> \
                make_const_view(StateVector && sv) const

    Make a read/write or constant view of :var:`sv`.
    Constant views will not be accepted by :func:`update_add_element()`.
    If :var:`sv` is not an lvalue reference, the resulting view will
    :ref:`hold a copy <mapped_basis_view_Ref>` of :var:`sv`.

    .. warning::

      To reduce memory footprint, :class:`mapped_basis_view` objects store
      a reference to the basis index map owned by their parent
      :class:`basis_mapper` object. For this reason, the views should never
      outlive the mapper.

  .. rubric:: Other methods

  .. function:: sv_index_type size() const

    Number of elements in the index map.

  .. function:: std::unordered_map<sv_index_type, sv_index_type> \
                const& map() const

    Direct access to the underlying index map.

  .. function:: std::unordered_map<sv_index_type, sv_index_type> \
                inverse_map() const

    Build and return an inverse index map. Depending on map's size, building
    the inverse can be an expensive operation. Calling this method on a
    non-invertible map is undefined behavior.

.. _n_fermion_sector_view:

N-fermion sector views
----------------------

There are two more specialised flavours of the basis mapping views called
:math:`N`-fermion sector views and :math:`N`-fermion multisector views. They
can come in handy when working with particle-number preserving models with
fermions. If a model involves :math:`M` fermionic degrees of freedom, then
storing a basis state index map for :class:`mapped_basis_view` requires
exponentially much memory, :math:`O(2^M)`.
It is possible to alleviate the memory consumption problem by employing a
(somewhat slower) algorithm that ranks bit patterns in the binary representation
of a basis state index. A computed rank is then used to index into the
:math:`N`-fermion (multi)sector.

.. warning::

  :class:`n_fermion_sector_view` and :class:`n_fermion_multisector_view` are
  incompatible with :ref:`sparse Hilbert spaces <sparse_hilbert_space>`.

.. class:: template<typename StateVector, \
                    bool Ref = true, \
                    typename RankingAlgorithm = combination_ranking> \
                    n_fermion_sector_view

  *Defined in <libcommute/loperator/n_fermion_sector_view.hpp>*

  View of a :type:`StateVector` object that translates basis state indices from
  a full :ref:`Hilbert space <hilbert_space>` to its subspace (sector) with a
  fixed total occupation of fermionic degrees of freedom :math:`N`. The full
  Hilbert space does not have to be purely fermionic.

  :type:`StateVector` - type of the underlying state vector object. Defining a
  read-only view (such that prohibits :expr:`update_add_element()` operations)
  requires using a ``const``-qualified type here. For example, one can use
  ``StateVector = std::vector<double>`` for a read-write view, and
  ``StateVector = const std::vector<double>`` for a read-only view.

  .. _n_fermion_sector_view_Ref:

  :type:`Ref` - by default, :class:`n_fermion_sector_view`
  stores a reference to the underlying state vector. Setting this option to
  ``false`` will result in a copy being created and stored instead. This feature
  can be useful when the underlying type is already a view-like object similar
  to ``Eigen::Map``.

  :type:`RankingAlgorithm` - one of the types implementing
  :ref:`bit pattern ranking <ranking_algorithms>`.

  .. function:: template <typename SV, typename HSType> \
                         n_fermion_sector_view(SV&& sv, \
                         HSType const& hs, unsigned int N)

    Construct a view of the state vector :var:`sv`, defined in the
    :var:`N`-fermion sector of the full Hilbert space :var:`hs`.

  .. function:: sv_index_type map_index(sv_index_type index) const

    Translate a basis state :var:`index` from the full Hilbert space to the
    sector.

.. struct:: template <typename HSType> sector_descriptor

  Description of an :math:`N`-fermion sector defined over a subset of fermionic
  degrees of freedom.

  :type:`HSType` - type of the full Hilbert space this sector belongs to.

  .. member:: std::set<typename HSType::index_types> indices

    Set of indices corresponding to the relevant fermionic degrees of freedom.

  .. member:: unsigned int N

    Total occupation of the sector.

.. class:: template<typename StateVector, \
                    bool Ref = true, \
                    typename RankingAlgorithm = combination_ranking> \
           n_fermion_multisector_view

  *Defined in <libcommute/loperator/n_fermion_sector_view.hpp>*

  View of a :type:`StateVector` object that translates basis state indices from
  a full :ref:`Hilbert space <hilbert_space>` to an :math:`N`-fermion
  multisector. A multisector is a set of all basis states, which have
  :math:`N_1` particles within a subset of fermionic modes :math:`\{S_1\}`,
  :math:`N_2` particles within another subset :math:`\{S_2\}` and so on. There
  can be any number of individual pairs :math:`(\{S_i\}, N_i)` (sectors
  contributing to the multisector) as long as all subsets :math:`\{S_i\}` are
  disjoint. The full Hilbert space does not have to be purely fermionic.

  It is advised to use :class:`n_fermion_sector_view` instead, if there is only
  one contributing sector that also spans all fermionic degrees of freedom.

  :type:`StateVector` - type of the underlying state vector object. Defining a
  read-only view (such that prohibits :expr:`update_add_element()` operations)
  requires using a ``const``-qualified type here. For example, one can use
  ``StateVector = std::vector<double>`` for a read-write view, and
  ``StateVector = const std::vector<double>`` for a read-only view.

  .. _n_fermion_multisector_view_Ref:

  :type:`Ref` - by default, :class:`n_fermion_multisector_view`
  stores a reference to the underlying state vector. Setting this option to
  ``false`` will result in a copy being created and stored instead. This feature
  can be useful when the underlying type is already a view-like object similar
  to ``Eigen::Map``.

  :type:`RankingAlgorithm` - one of the types implementing
  :ref:`bit pattern ranking <ranking_algorithms>`.

  .. function:: template <typename SV, typename HSType> \
                n_fermion_multisector_view(SV&& sv, HSType const& hs, \
                std::vector<sector_descriptor<HSType>> const& sectors)

    Construct a view of the state vector :var:`sv`, defined in the
    :math:`N`-fermion multisector of the full Hilbert space :var:`hs`.
    The multisector is defined via a list of contributing :var:`sectors`
    (list of :math:`(\{S_i\}, N_i)` pairs).

  .. function:: sv_index_type map_index(sv_index_type index) const

    Translate a basis state :var:`index` from the full Hilbert space to the
    multisector.

Besides :class:`n_fermion_sector_view` and :class:`n_fermion_multisector_view`,
*<libcommute/loperator/n_fermion_sector_view.hpp>* defines a few supplemental
utility functions that help working with (multi)sectors.

.. function:: template <typename StateVector, typename HSType> \
              auto make_nfs_view(StateVector&& sv, HSType const& hs, \
              unsigned int N)
              template <typename StateVector, typename HSType> \
              auto make_const_nfs_view(StateVector&& sv, HSType const& hs, \
              unsigned int N)

  Make and return a read/write or constant :var:`N`-fermion sector view of
  :var:`sv` within the full Hilbert space :var:`hs`. If :var:`sv` is not an
  lvalue reference, the resulting view will
  :ref:`hold a copy <n_fermion_sector_view_Ref>` of :var:`sv`.
  A returned view uses :class:`combination_ranking` as its bit pattern ranking
  algorithm.

.. function:: template <typename StateVector, typename HSType> \
              auto make_nfms_view(StateVector&& sv, HSType const& hs, \
              std::vector<sector_descriptor<HSType>> const& sectors)
              template <typename StateVector, typename HSType> \
              auto make_const_nfms_view(StateVector&& sv, HSType const& hs, \
              std::vector<sector_descriptor<HSType>> const& sectors)

  Make and return a read/write or constant :math:`N`-fermion multisector view of
  :var:`sv` within the full Hilbert space :var:`hs`. The multisector is
  defined via a list of contributing :var:`sectors` (list of
  :math:`(\{S_i\}, N_i)` pairs). If :var:`sv` is not an lvalue reference,
  the resulting view will
  :ref:`hold a copy <n_fermion_sector_view_Ref>` of :var:`sv`.
  A returned view uses :class:`combination_ranking` as its bit pattern ranking
  algorithm.

.. function:: template <typename HSType> sv_index_type \
              n_fermion_sector_size(HSType const& hs, unsigned int N)

  Size of the :var:`N`-fermion sector within the full Hilbert space :var:`hs`.

.. function:: template <typename HSType> sv_index_type \
              n_fermion_multisector_size(HSType const& hs, \
              std::vector<sector_descriptor<HSType>> const& sectors)

  Size of the :math:`N`-fermion multisector within the full Hilbert space
  :var:`hs`. The multisector is defined via a list of contributing
  :var:`sectors` (list of :math:`(\{S_i\}, N_i)` pairs).

.. function:: template <typename HSType> std::vector<sv_index_type> \
              n_fermion_sector_basis_states(HSType const& hs, unsigned int N)

  Build and return a list of basis state indices forming the :var:`N`-fermion
  sector within the full Hilbert space :var:`hs`. The order of the indices in
  the list is consistent with the results of
  :func:`n_fermion_sector_view::map_index()`.

  .. code-block:: cpp

    auto basis_states = n_fermion_sector_basis_states(hs, N);
    auto view = n_fermion_sector_view(st, hs, N);

    for(sv_index_type n = 0; n < basis_states.size(); ++n) {
      view.map_index(basis_states[n]) == n; // true for all n
    }

.. function:: template <typename HSType> std::vector<sv_index_type> \
              n_fermion_multisector_basis_states(HSType const& hs, \
              std::vector<sector_descriptor<HSType>> const& sectors)

  Build and return a list of basis state indices forming an :math:`N`-fermion
  multisector within the full Hilbert space :var:`hs`. The multisector is
  defined via a list of contributing :var:`sectors` (list of
  :math:`(\{S_i\}, N_i)` pairs). The order of the indices in the list is
  consistent with the results of
  :func:`n_fermion_multisector_view::map_index()`.

  .. code-block:: cpp

    auto basis_states = n_fermion_multisector_basis_states(hs, sectors);
    auto view = n_fermion_multisector_view(st, hs, sectors);

    for(sv_index_type n = 0; n < basis_states.size(); ++n) {
      view.map_index(basis_states[n]) == n; // true for all n
    }

.. _ranking_algorithms:

The following classes implement the three ranking algorithms described in
[WH22]_. They precompute and store a certain amount of data in order to speed up
calculations.

*
  .. class:: combination_ranking

    *Defined in <libcommute/loperator/n_fermion_sector_view.hpp>*

    The ranking algorithm based on the combinatorial number system.
    The :math:`N`-fermion (multi)sector view types use this algorithm by
    default. The storage space required by this class scales as
    :math:`O(M \min(N, M - N))`, where :math:`M` is the total number of
    the fermionic degrees of freedom.

*
  .. class:: template <unsigned int R> staggered_ranking

    *Defined in <libcommute/loperator/n_fermion_sector_view.hpp>*

    The improved combinatorial ranking with staggered lookup and a chunk size of
    :var:`R` bits.
    The storage space required by this class scales as
    :math:`O\left(2^R (M-R+2)(\frac{M}{2R}+1)\right)`, where :math:`M` is the
    total number of the fermionic degrees of freedom.

*
  .. class:: template <unsigned int R> trie_ranking

    *Defined in <libcommute/loperator/n_fermion_sector_view.hpp>*

    The trie-based ranking algorithm with a chunk size of :var:`R` bits.
    The storage space required by this class is roughly proportional to the size
    of the (multi)sector.

.. [WH22] "Trie-based ranking of quantum many-body states",
   M. Wallerberger and K. Held,
   Phys. Rev. Research **4**, 033238 (2022),
   https://doi.org/10.1103/PhysRevResearch.4.033238