space_partition.hpp

Go to the documentation of this file.

// Copyright (c) 2015-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2015-2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2023 Simons Foundation
// Copyright (c) 2015-2016 Igor Krivenko
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You may obtain a copy of the License at
//     https://www.gnu.org/licenses/gpl-3.0.txt
//
// Authors: Maxime Charlebois, Michel Ferrero, Igor Krivenko, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "./state.hpp"
#include "../utility/numeric_ops.hpp"
 
#include <boost/pending/disjoint_sets.hpp>
 
#include <cstdint>
#include <map>
#include <set>
#include <utility>
 
namespace triqs::hilbert_space {


  template <typename S, typename OP> class space_partition {
    public:
    using idx_t = uint32_t;

    using state_t = S;

    using operator_t = OP;

    using amplitude_t = typename state_t::value_type;

    using block_mapping_t = std::set<std::pair<idx_t, idx_t>>;

    using matrix_element_map_t = std::map<std::pair<idx_t, idx_t>, amplitude_t>;

    space_partition(state_t const &psi, operator_t const &H, bool store_matrix_elements = true, operator_t const &delta = operator_t())
       : tmp_state(make_zero_state(psi)), subspaces(psi.size()) {
      auto size = tmp_state.size();
 
      // Iteration over all initial basis states
      for (idx_t i = 0; i < size; ++i) {
        tmp_state(i)        = amplitude_t(1);
        state_t final_state = H(tmp_state);
 
        auto mapping = [&](idx_t f, amplitude_t amplitude) {
          using triqs::utility::is_zero;
          if (is_zero(amplitude)) return;
          auto i_subspace = subspaces.find_set(i);
          auto f_subspace = subspaces.find_set(f);
          if (i_subspace != f_subspace) subspaces.link(i_subspace, f_subspace);
 
          if (store_matrix_elements) matrix_elements[std::make_pair(i, f)] = amplitude;
        };
 
        // Iterate over non-zero final amplitudes
        foreach (final_state, mapping);
 
        // redo for additionnal Hyb
        if (not delta.is_empty()) {
          final_state = delta(tmp_state);
          foreach (final_state, mapping);
        }
 
        tmp_state(i) = amplitude_t(0.);
      }
 
      _update_index();
    }

    space_partition(space_partition const &) = default;

    auto merge_subspaces(operator_t const &cd, operator_t const &c, bool store_matrix_elements = true) {
 
      matrix_element_map_t Cd_elements, C_elements;
      std::multimap<idx_t, idx_t> Cd_connections, C_connections;
 
      auto size = tmp_state.size();
 
      // Fill connection multimaps
      for (idx_t i = 0; i < size; ++i) {
        tmp_state(i)    = amplitude_t(1);
        auto i_subspace = subspaces.find_set(i);
 
        auto fill_conn = [this, i, i_subspace, store_matrix_elements](operator_t const &op, std::multimap<idx_t, idx_t> &conn,
                                                                      matrix_element_map_t &elem) {
          state_t final_state = op(tmp_state);
          // Iterate over non-zero final amplitudes
          foreach (final_state, [&](idx_t f, amplitude_t amplitude) {
            using triqs::utility::is_zero;
            if (is_zero(amplitude)) return;
            auto f_subspace = subspaces.find_set(f);
            conn.insert({i_subspace, f_subspace});
            if (store_matrix_elements) elem[{i, f}] = amplitude;
          });
        };
 
        fill_conn(cd, Cd_connections, Cd_elements);
        fill_conn(c, C_connections, C_elements);
 
        tmp_state(i) = amplitude_t(0.);
      }
 
      // 'Zigzag' traversal algorithm
      while (!Cd_connections.empty()) {
 
        // Take one C^+ - connection
        // C^+|lower_subspace> = |upper_subspace>
        auto const [lower_subspace, upper_subspace] = *std::begin(Cd_connections);
 
        // - Reveals all subspaces reachable from lower_subspace by application of
        //   a 'zigzag' product C^+ C C^+ C C^+ ... of any length.
        // - Removes all visited connections from Cd_connections/C_connections.
        // - Merges lower_subspace with all subspaces generated from lower_subspace by application of (C C^+)^(2*n).
        // - Merges upper_subspace with all subspaces generated from upper_subspace by application of (C^+ C)^(2*n).
        std::function<void(idx_t, bool)> zigzag_traversal = [this, lower_subspace, upper_subspace, &Cd_connections, &C_connections,
                                                             &zigzag_traversal](idx_t i_subspace, // find all connections starting from i_subspace
                                                                                bool upwards      // if true, C^+ connection, otherwise C connection
                                                            ) {
          std::multimap<idx_t, idx_t>::iterator it;
          while ((it = (upwards ? Cd_connections : C_connections).find(i_subspace)) != (upwards ? Cd_connections : C_connections).end()) {
 
            auto f_subspace = it->second;
            (upwards ? Cd_connections : C_connections).erase(it);
 
            if (upwards)
              subspaces.link(f_subspace, upper_subspace);
            else
              subspaces.link(f_subspace, lower_subspace);
 
            // Recursively apply to all found f_subspace's with a 'flipped' direction
            zigzag_traversal(f_subspace, !upwards);
          }
        };
 
        // Apply to all C^+ connections starting from lower_subspace
        zigzag_traversal(lower_subspace, true);
      }
 
      _update_index();
 
      return std::make_pair(Cd_elements, C_elements);
    }

    idx_t n_subspaces() const { return representative_to_index.size(); }

    template <typename F> friend void foreach (space_partition &sp, F f) {
      for (idx_t n = 0; n < sp.tmp_state.size(); ++n) f(n, sp.lookup_basis_state(n));
    }

    idx_t lookup_basis_state(idx_t f) { return representative_to_index[subspaces.find_set(f)]; }

    matrix_element_map_t const &get_matrix_elements() const { return matrix_elements; }

    block_mapping_t find_mappings(operator_t const &op, bool diagonal_only = false) {
 
      block_mapping_t mapping;
 
      // Iteration over all initial basis states
      for (idx_t i = 0; i < tmp_state.size(); ++i) {
        state_t initial_state = tmp_state;
        initial_state(i)      = amplitude_t(1);
        auto i_subspace       = subspaces.find_set(i);
 
        state_t final_state = op(initial_state);
 
        // Iterate over non-zero final amplitudes
        foreach (final_state, [&](idx_t f, amplitude_t amplitude) {
          using triqs::utility::is_zero;
          if (is_zero(amplitude)) return;
          auto f_subspace = subspaces.find_set(f);
          if ((!diagonal_only) || i_subspace == f_subspace)
            mapping.insert(std::make_pair(representative_to_index[i_subspace], representative_to_index[f_subspace]));
        });
      }
 
      return mapping;
    }
 
    private:
    void _update_index() {
      auto p = subspaces.parents();
      subspaces.compress_sets(p.begin(), p.end());  // parents are representatives
      subspaces.normalize_sets(p.begin(), p.end()); // the representative has the smallest index in the set
 
      // Update representative_to_index
      representative_to_index.clear();
      for (idx_t n = 0; n < tmp_state.size(); ++n) {
        representative_to_index.insert(std::make_pair(subspaces.find_set(n), representative_to_index.size()));
      }
    }
 
    private:
    mutable state_t tmp_state;
    boost::disjoint_sets_with_storage<> subspaces;
    matrix_element_map_t matrix_elements;
    std::map<idx_t, idx_t> representative_to_index;
  };

 
} // namespace triqs::hilbert_space