block_space.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <uni20/symmetry/block_sector.hpp>
#include <uni20/symmetry/irregular_space.hpp>
#include <uni20/symmetry/local_space.hpp>
#include <uni20/symmetry/qnum.hpp>
#include <uni20/symmetry/space.hpp>
 
#include <algorithm>
#include <concepts>
#include <cstddef>
#include <initializer_list>
#include <numeric>
#include <ranges>
#include <span>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
 
namespace uni20
{
 
class BlockSpace {
  public:
    explicit BlockSpace(Symmetry sym, std::string label = {}) : sym_(sym), label_(std::move(label))
    {
      this->verify_symmetry();
    }
 
    BlockSpace(Symmetry sym, std::initializer_list<BlockSector> sectors, std::string label = {})
        : BlockSpace(sym, std::move(label))
    {
      this->append_sectors(sectors);
      this->canonicalize_sectors();
    }
 
    template <std::ranges::input_range Sectors>
      requires std::constructible_from<BlockSector, std::ranges::range_reference_t<Sectors>>
    BlockSpace(Symmetry sym, Sectors&& sectors, std::string label = {}) : BlockSpace(sym, std::move(label))
    {
      this->append_sectors(std::forward<Sectors>(sectors));
      this->canonicalize_sectors();
    }
 
    auto symmetry() const -> Symmetry { return sym_; }
 
    auto size() const -> std::size_t { return sectors_.size(); }
 
    auto empty() const -> bool { return sectors_.empty(); }
 
    auto total_dim() const -> std::size_t
    {
      return std::accumulate(sectors_.begin(), sectors_.end(), std::size_t{0},
                             [](std::size_t total, BlockSector const& sector) { return total + sector.dim; });
    }
 
    auto contains(QNum q) const -> bool
    {
      this->verify_symmetry(q);
      return std::find_if(sectors_.begin(), sectors_.end(), [&](BlockSector const& sector) { return sector.q == q; }) !=
             sectors_.end();
    }
 
    auto operator[](std::size_t index) const -> BlockSector const& { return sectors_[index]; }
 
    auto sectors() const -> std::span<BlockSector const> { return sectors_; }
 
    auto begin() const { return sectors_.begin(); }
 
    auto end() const { return sectors_.end(); }
 
    auto label() const -> std::string const& { return label_; }
 
    void set_label(std::string label) { label_ = std::move(label); }
 
    auto operator==(BlockSpace const& other) const -> bool = default;
 
  private:
    void verify_symmetry() const
    {
      if (!sym_.valid())
      {
        throw std::invalid_argument("BlockSpace requires a valid symmetry");
      }
    }
 
    void verify_symmetry(QNum const& q) const
    {
      if (!q.valid() || q.symmetry() != sym_)
      {
        throw std::invalid_argument("BlockSpace sector has the wrong symmetry");
      }
    }
 
    void verify_sector(BlockSector const& sector) const
    {
      this->verify_symmetry(sector.q);
      if (sector.dim == 0)
      {
        throw std::invalid_argument("BlockSpace sectors must have positive dimension");
      }
    }
 
    template <std::ranges::input_range Sectors> void append_sectors(Sectors&& sectors)
    {
      for (auto&& input : sectors)
      {
        BlockSector sector(std::forward<decltype(input)>(input));
        this->verify_sector(sector);
        sectors_.push_back(std::move(sector));
      }
    }
 
    void canonicalize_sectors()
    {
      std::sort(sectors_.begin(), sectors_.end(),
                [](BlockSector const& lhs, BlockSector const& rhs) { return lhs.q.raw_code() < rhs.q.raw_code(); });
 
      auto const duplicate =
          std::adjacent_find(sectors_.begin(), sectors_.end(),
                             [](BlockSector const& lhs, BlockSector const& rhs) { return lhs.q == rhs.q; });
      if (duplicate != sectors_.end())
      {
        throw std::invalid_argument("BlockSpace sectors must have unique quantum numbers");
      }
    }
 
    Symmetry sym_;
    std::vector<BlockSector> sectors_;
    std::string label_;
};
 
static_assert(SymmetrySpace<BlockSpace>);
 
struct QNumListRegularization
{
    QNumList original;
 
    BlockSpace regular;
 
    std::vector<std::size_t> block_index;
 
    std::vector<std::size_t> block_offset;
};
 
struct LocalSpaceRegularization
{
    LocalSpace original;
 
    BlockSpace regular;
 
    std::vector<std::size_t> block_index;
 
    std::vector<std::size_t> block_offset;
};
 
struct BlockRangeMapping
{
    std::size_t block_index;
 
    std::size_t first;
 
    std::size_t last;
 
    auto size() const -> std::size_t { return last - first; }
 
    auto operator==(BlockRangeMapping const& other) const -> bool = default;
};
 
struct IrregularSpaceRegularization
{
    IrregularSpace original;
 
    BlockSpace regular;
 
    std::vector<BlockRangeMapping> block_mapping;
};
 
namespace detail
{
 
inline auto qnum_code_less(QNum const& lhs, QNum const& rhs) -> bool { return lhs.raw_code() < rhs.raw_code(); }
 
template <typename Range> inline auto unique_qnums(Range const& range) -> std::vector<QNum>
{
  std::unordered_set<QNum> seen;
  std::vector<QNum> qnums;
  for (QNum const& q : range)
  {
    if (seen.insert(q).second)
    {
      qnums.push_back(q);
    }
  }
  std::sort(qnums.begin(), qnums.end(), qnum_code_less);
  return qnums;
}
 
} // namespace detail
 
inline auto regularize(QNumList const& list) -> QNumListRegularization
{
  auto const qnums = detail::unique_qnums(list);
 
  std::unordered_map<QNum, std::size_t> count_by_q;
  for (QNum const& q : list)
  {
    count_by_q[q] += 1;
  }
 
  std::vector<BlockSector> sectors;
  sectors.reserve(qnums.size());
  for (QNum const& q : qnums)
  {
    sectors.push_back({q, count_by_q[q]});
  }
  BlockSpace regular(list.symmetry(), sectors);
 
  std::unordered_map<QNum, std::size_t> index_by_q;
  for (std::size_t index = 0; index < regular.size(); ++index)
  {
    index_by_q.emplace(regular[index].q, index);
  }
 
  std::unordered_map<QNum, std::size_t> offset_by_q;
  std::vector<std::size_t> block_index;
  std::vector<std::size_t> block_offset;
  block_index.reserve(list.size());
  block_offset.reserve(list.size());
 
  for (QNum const& q : list)
  {
    block_index.push_back(index_by_q.at(q));
    block_offset.push_back(offset_by_q[q]);
    offset_by_q[q] += 1;
  }
 
  return {list, std::move(regular), std::move(block_index), std::move(block_offset)};
}
 
inline auto regularize(LocalSpace const& space) -> LocalSpaceRegularization
{
  auto const qnums = detail::unique_qnums(space);
 
  std::unordered_map<QNum, std::size_t> count_by_q;
  for (QNum const& q : space)
  {
    count_by_q[q] += 1;
  }
 
  std::vector<BlockSector> sectors;
  sectors.reserve(qnums.size());
  for (QNum const& q : qnums)
  {
    sectors.push_back({q, count_by_q[q]});
  }
  BlockSpace regular(space.symmetry(), sectors, space.label());
 
  std::unordered_map<QNum, std::size_t> index_by_q;
  for (std::size_t index = 0; index < regular.size(); ++index)
  {
    index_by_q.emplace(regular[index].q, index);
  }
 
  std::unordered_map<QNum, std::size_t> offset_by_q;
  std::vector<std::size_t> block_index;
  std::vector<std::size_t> block_offset;
  block_index.reserve(space.size());
  block_offset.reserve(space.size());
 
  for (QNum const& q : space)
  {
    block_index.push_back(index_by_q.at(q));
    block_offset.push_back(offset_by_q[q]);
    offset_by_q[q] += 1;
  }
 
  return {space, std::move(regular), std::move(block_index), std::move(block_offset)};
}
 
inline auto regularize(IrregularSpace const& space) -> IrregularSpaceRegularization
{
  std::vector<QNum> original_qnums;
  original_qnums.reserve(space.size());
  for (BlockSector const& sector : space)
  {
    original_qnums.push_back(sector.q);
  }
  auto const qnums = detail::unique_qnums(original_qnums);
 
  std::unordered_map<QNum, std::size_t> dim_by_q;
  for (BlockSector const& sector : space)
  {
    dim_by_q[sector.q] += sector.dim;
  }
 
  std::vector<BlockSector> sectors;
  sectors.reserve(qnums.size());
  for (QNum const& q : qnums)
  {
    sectors.push_back({q, dim_by_q[q]});
  }
  BlockSpace regular(space.symmetry(), sectors, space.label());
 
  std::unordered_map<QNum, std::size_t> index_by_q;
  for (std::size_t index = 0; index < regular.size(); ++index)
  {
    index_by_q.emplace(regular[index].q, index);
  }
 
  std::unordered_map<QNum, std::size_t> offset_by_q;
  std::vector<BlockRangeMapping> block_mapping;
  block_mapping.reserve(space.size());
  for (BlockSector const& sector : space)
  {
    auto& offset = offset_by_q[sector.q];
    block_mapping.push_back({index_by_q.at(sector.q), offset, offset + sector.dim});
    offset += sector.dim;
  }
 
  return {space, std::move(regular), std::move(block_mapping)};
}
 
inline auto to_block_space(QNumList const& list) -> BlockSpace { return regularize(list).regular; }
 
inline auto to_block_space(LocalSpace const& space) -> BlockSpace { return regularize(space).regular; }
 
inline auto to_block_space(IrregularSpace const& space) -> BlockSpace { return regularize(space).regular; }
 
} // namespace uni20