many_body_operator.hpp

Go to the documentation of this file.

// Copyright (c) 2014-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2014-2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2023 Simons Foundation
// Copyright (c) 2014-2017 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: Michel Ferrero, Alexander Hampel, Igor Krivenko, Olivier Parcollet, Hugo U. R. Strand, Nils Wentzell

 
#pragma once
 
#include "../hilbert_space/fundamental_operator_set.hpp"
#include "../utility/dressed_iterator.hpp"
#include "../utility/real_or_complex.hpp"
#include "../utility/numeric_ops.hpp"
#include "../utility/variant_extensions.hpp"
 
#include <h5/h5.hpp>
 
#include <algorithm>
#include <cmath>
#include <compare>
#include <complex>
#include <cstddef>
#include <functional>
#include <map>
#include <ostream>
#include <string>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>
 
namespace triqs::operators {

 
  // Elevate triqs::utility::real_or_complex to the `triqs::operators` namespace.
  using utility::real_or_complex;

  using indices_t = hilbert_space::fundamental_operator_set::indices_t;
 
  // Forward declarations.
  template <typename T> class many_body_operator_generic;
  template <typename T> many_body_operator_generic<T> real(many_body_operator_generic<T> const &op);
  template <typename T> many_body_operator_generic<T> imag(many_body_operator_generic<T> const &op);

  using many_body_operator = many_body_operator_generic<real_or_complex>;

  using many_body_operator_real = many_body_operator_generic<double>;

  using many_body_operator_complex = many_body_operator_generic<std::complex<double>>;

  struct canonical_ops_t {
    bool dagger{false};

    indices_t indices{};

    auto operator<=>(canonical_ops_t const &b) const {
      if (dagger != b.dagger) return !dagger <=> !b.dagger; // c+ < c
 
      // Compare indices element-wise, ordering by type index first (long < string < double < array)
      auto compare_indices = [](indices_t const &a, indices_t const &b) -> std::strong_ordering {
        for (size_t i = 0; i < std::min(a.size(), b.size()); ++i) {
          // Compare by type index first
          if (a[i].index() != b[i].index()) return a[i].index() <=> b[i].index();
          // Same type: compare values
          auto cmp = std::visit(
             [](auto const &x, auto const &y) -> std::strong_ordering {
               if constexpr (std::is_same_v<std::decay_t<decltype(x)>, std::decay_t<decltype(y)>>) {
                 if (x == y) return std::strong_ordering::equal;
                 return x < y ? std::strong_ordering::less : std::strong_ordering::greater;
               }
               return std::strong_ordering::equal; // unreachable due to index check
             },
             a[i], b[i]);
          if (cmp != std::strong_ordering::equal) return cmp;
        }
        return a.size() <=> b.size();
      };
 
      return dagger ? compare_indices(indices, b.indices) : compare_indices(b.indices, indices);
    }

    bool operator==(canonical_ops_t const &b) const { return (*this <=> b) == 0; }

    void serialize(auto &ar) const { ar & dagger & indices; }

    void deserialize(auto &ar) { ar & dagger & indices; }
  };

  std::ostream &operator<<(std::ostream &os, canonical_ops_t const &op);

  using monomial_t = std::vector<canonical_ops_t>;

  bool operator<(monomial_t const &m1, monomial_t const &m2);

  std::ostream &operator<<(std::ostream &os, monomial_t const &m);

  template <typename T> class many_body_operator_generic {
    public:
    using monomials_map_t = std::map<monomial_t, T>;

    using scalar_t = T;

    many_body_operator_generic() = default;

    template <typename S> many_body_operator_generic(many_body_operator_generic<S> const &p) {
      static_assert(std::is_constructible_v<scalar_t, S>);
      *this = p;
    }

    explicit many_body_operator_generic(scalar_t const &x) {
      using triqs::utility::is_zero;
      if (!is_zero(x)) monomials_.insert({{}, x});
    }

    many_body_operator_generic(scalar_t const &x, monomial_t monomial) {
      using triqs::utility::is_zero;
      if (!is_zero(x)) monomials_.emplace(std::move(monomial), x);
    }
 
    struct _cdress;
    C2PY_IGNORE many_body_operator_generic(_cdress const &term) { normalize_and_insert(term.monomial, term.coef, monomials_); }

    template <typename S> many_body_operator_generic &operator=(many_body_operator_generic<S> const &p) {
      static_assert(std::is_constructible_v<scalar_t, S>);
      monomials_.clear();
      for (auto const &y : p.get_monomials()) monomials_.insert(std::make_pair(monomial_t{y.first}, scalar_t(y.second)));
      return *this;
    }

    [[nodiscard]] monomials_map_t const &get_monomials() const { return monomials_; }

    [[nodiscard]] hilbert_space::fundamental_operator_set make_fundamental_operator_set() const {
      hilbert_space::fundamental_operator_set fops;
      for (auto const &m : monomials_)
        for (auto const &c_cdag_op : m.first) fops.insert_from_indices_t(c_cdag_op.indices);
      return fops;
    }

    static many_body_operator_generic make_canonical(bool is_dag, indices_t indices) {
      many_body_operator_generic res;
      auto m = monomial_t{canonical_ops_t{.dagger = is_dag, .indices = indices}};
      res.monomials_.insert({m, scalar_t(1)});
      return res;
    }
 
    // We use utility::dressed_iterator to dress iterators.
    // _cdress is a simple struct of refs that exposes (monomial, coef) to the iterator.
    struct C2PY_IGNORE _cdress {
      monomial_t const &monomial; // NOLINT
      scalar_t coef;
      _cdress(typename monomials_map_t::const_iterator _it) : monomial(_it->first), coef(_it->second) {}
      operator std::pair<std::vector<std::pair<bool, indices_t>>, scalar_t>() const {
        std::vector<std::pair<bool, indices_t>> tmp_monomial;
        tmp_monomial.reserve(monomial.size());
        for (auto cop : monomial) tmp_monomial.emplace_back(cop.dagger, cop.indices);
        return {tmp_monomial, coef};
      }
    };

    using const_iterator = utility::dressed_iterator<typename monomials_map_t::const_iterator, _cdress>;

    [[nodiscard]] const_iterator begin() const noexcept { return monomials_.begin(); }

    [[nodiscard]] const_iterator end() const noexcept { return monomials_.end(); }

    [[nodiscard]] const_iterator cbegin() const noexcept { return monomials_.cbegin(); }

    [[nodiscard]] const_iterator cend() const noexcept { return monomials_.cend(); }

    [[nodiscard]] bool is_almost_zero(double precision = 1e-10) const {
      auto term_is_zero = [precision](auto const &term) {
        using std::abs;
        return triqs::utility::is_zero(abs(term.coef), precision);
      };
      return std::all_of(this->begin(), this->end(), term_is_zero);
    }

    [[nodiscard]] bool is_zero() const { return monomials_.empty(); }

    many_body_operator_generic operator-() const {
      auto res = *this;
      for (auto &m : res.monomials_) m.second = -m.second;
      return res;
    }

    many_body_operator_generic &operator+=(scalar_t a) {
      using triqs::utility::is_zero;
      if (is_zero(a)) return *this;
      bool is_new_monomial{};
      typename monomials_map_t::iterator it;
      std::tie(it, is_new_monomial) = monomials_.insert(std::make_pair(monomial_t(0), a));
      if (!is_new_monomial) {
        it->second += a;
        erase_zero_monomial(monomials_, it);
      }
      return *this;
    }

    many_body_operator_generic &operator-=(scalar_t a) { return operator+=(-a); }

    many_body_operator_generic &operator*=(scalar_t a) {
      using triqs::utility::is_zero;
      if (is_zero(a)) {
        monomials_.clear();
      } else {
        for (auto &m : monomials_) m.second *= a;
      }
      return *this;
    }

    many_body_operator_generic &operator/=(scalar_t a) { return operator*=(scalar_t(1) / a); }

    friend many_body_operator_generic operator+(many_body_operator_generic op, scalar_t a) { return op += a; }

    friend many_body_operator_generic operator+(scalar_t a, many_body_operator_generic op) { return op += a; }

    friend many_body_operator_generic operator-(many_body_operator_generic op, scalar_t a) { return op -= a; }

    friend many_body_operator_generic operator-(scalar_t a, many_body_operator_generic const &op) { return -op + a; }

    friend many_body_operator_generic operator*(many_body_operator_generic op, scalar_t a) { return op *= a; }

    friend many_body_operator_generic operator*(scalar_t a, many_body_operator_generic op) { return op *= a; }

    friend many_body_operator_generic operator/(many_body_operator_generic op, scalar_t a) { return op /= a; }

    many_body_operator_generic &operator+=(many_body_operator_generic const &op) {
      bool is_new_monomial{};
      typename monomials_map_t::iterator it;
      for (auto const &m : op.monomials_) {
        std::tie(it, is_new_monomial) = monomials_.insert(m);
        if (!is_new_monomial) {
          it->second += m.second;
          erase_zero_monomial(monomials_, it);
        }
      }
      return *this;
    }

    many_body_operator_generic &operator-=(many_body_operator_generic const &op) {
      bool is_new_monomial{};
      typename monomials_map_t::iterator it;
      for (auto const &m : op.monomials_) {
        std::tie(it, is_new_monomial) = monomials_.insert(std::make_pair(m.first, -m.second));
        if (!is_new_monomial) {
          it->second -= m.second;
          erase_zero_monomial(monomials_, it);
        }
      }
      return *this;
    }

    many_body_operator_generic &operator*=(many_body_operator_generic const &op) {
      monomials_map_t tmp_map; // product will be stored here
      for (auto const &m : monomials_)
        for (auto const &op_m : op.monomials_) {
          monomial_t product_m;
          product_m.reserve(m.first.size() + op_m.first.size());
          for (auto const &tmp_op : m.first) product_m.push_back(tmp_op);
          for (auto const &tmp_op : op_m.first) product_m.push_back(tmp_op);
          normalize_and_insert(product_m, m.second * op_m.second, tmp_map);
        }
      std::swap(monomials_, tmp_map);
      return *this;
    }

    friend many_body_operator_generic operator+(many_body_operator_generic lhs, many_body_operator_generic const &rhs) { return lhs += rhs; }

    friend many_body_operator_generic operator-(many_body_operator_generic lhs, many_body_operator_generic const &rhs) { return lhs -= rhs; }

    friend many_body_operator_generic operator*(many_body_operator_generic lhs, many_body_operator_generic const &rhs) { return lhs *= rhs; }

    bool operator==(many_body_operator_generic const &op) const { return (*this - op).is_zero(); }
 
    private:
    // Implementation details of dagger.
    static canonical_ops_t _dagger(canonical_ops_t const &cop) { return {!cop.dagger, cop.indices}; }
 
    static monomial_t _dagger(monomial_t const &m) {
      monomial_t res;
      for (auto it = m.rbegin(); it != m.rend(); ++it) res.push_back(_dagger(*it));
      return res;
    }
 
    public:
    friend many_body_operator_generic dagger(many_body_operator_generic const &op) {
      many_body_operator_generic res;
      using triqs::utility::conj;
      for (auto const &x : op) res.monomials_.insert({_dagger(x.monomial), conj(x.coef)});
      return res;
    }

    template <typename F> friend many_body_operator_generic transform(many_body_operator_generic const &op, F &&f) { // NOLINT
      many_body_operator_generic res;
      using triqs::utility::is_zero;
      for (auto const &x : op) {
        auto c = f(x.monomial, x.coef);
        if (!is_zero(c)) res.monomials_.insert({x.monomial, c});
      }
      return res;
    }

    [[nodiscard]] C2PY_PROPERTY_GET(real) many_body_operator_generic real() const { return operators::real(*this); }

    [[nodiscard]] C2PY_PROPERTY_GET(imag) many_body_operator_generic imag() const { return operators::imag(*this); }

    friend std::ostream &operator<<(std::ostream &os, many_body_operator_generic const &op) {
      if (op.monomials_.size() != 0) {
        bool print_plus = false;
        for (auto const &m : op.monomials_) {
          os << (print_plus ? " + " : "") << m.second;
          os << m.first;
          print_plus = true;
        }
      } else
        os << "0";
      return os;
    }

    void serialize(auto &ar) const { ar & monomials_; }

    void deserialize(auto &ar) { ar & monomials_; }

    [[nodiscard]] static std::string hdf5_format() { return "Operator"; }

    friend void h5_write(h5::group g, std::string const &name, many_body_operator const &op, hilbert_space::fundamental_operator_set const &fops);

    friend void h5_write(h5::group g, std::string const &name, many_body_operator_generic const &op) {
      h5_write(g, name, op, op.make_fundamental_operator_set());
    }

    friend void h5_read(h5::group g, std::string const &name, many_body_operator &op, hilbert_space::fundamental_operator_set &fops);

    friend void h5_read(h5::group g, std::string const &name, many_body_operator_generic &op) {
      hilbert_space::fundamental_operator_set fops;
      many_body_operator op_real_cplx;
      h5_read(g, name, op_real_cplx, fops);
      op = std::move(op_real_cplx);
    }
 
    private:
    // Normalize a monomial and insert into a map
    static void normalize_and_insert(monomial_t m, scalar_t coeff, monomials_map_t &target) {
      // The normalization is done by employing a simple bubble sort algorithms.
      // Apart from sorting elements this function keeps track of the sign and
      // recursively calls itself if a permutation of two operators produces a new
      // monomial
      if (m.size() >= 2) {
        bool is_swapped{};
        do { // NOLINT
          is_swapped = false;
          for (int n = 1; n < m.size(); ++n) {
            canonical_ops_t &prev_index = m[n - 1];
            canonical_ops_t &cur_index  = m[n];
            if (prev_index == cur_index) return; // The monomial is effectively zero
            if (prev_index > cur_index) {
              // Are we swapping C and C^+ with the same indices?
              // A bit ugly ...
              canonical_ops_t cur_index_flipped_type(cur_index);
              cur_index_flipped_type.dagger = !cur_index_flipped_type.dagger;
              if (prev_index == cur_index_flipped_type) {
                monomial_t new_m;
                new_m.reserve(m.size() - 2);
                std::copy(m.begin(), m.begin() + n - 1, std::back_inserter(new_m));
                std::copy(m.begin() + n + 1, m.end(), std::back_inserter(new_m));
                normalize_and_insert(new_m, coeff, target);
              }
              coeff = -coeff;
              std::swap(prev_index, cur_index);
              is_swapped = true;
            }
          }
        } while (is_swapped);
      }
 
      // Insert the result
      bool is_new_monomial{};
      typename monomials_map_t::iterator it;
      std::tie(it, is_new_monomial) = target.insert(std::make_pair(std::move(m), coeff));
      if (!is_new_monomial) {
        it->second += coeff;
        erase_zero_monomial(target, it);
      }
    }
 
    // Erase a monomial with a close-to-zero coefficient.
    static void erase_zero_monomial(monomials_map_t &m, typename monomials_map_t::iterator &it) {
      using triqs::utility::is_zero;
      if (is_zero(it->second)) m.erase(it);
    }
 
    private:
    monomials_map_t monomials_;
  };

  template <typename T1, typename T2>
  void assert_operators_are_close(many_body_operator_generic<T1> const &op1, many_body_operator_generic<T2> const &op2, double precision) {
    if (!(op1 - op2).is_almost_zero(precision)) TRIQS_RUNTIME_ERROR << "Error in operators::assert_operators_are_close: Terms are different";
  }

  template <typename T> many_body_operator_generic<T> real(many_body_operator_generic<T> const &op) {
    return transform(op, [](monomial_t const &, T a_i) {
      using triqs::utility::real;
      return real(a_i);
    });
  }

  template <typename T> many_body_operator_generic<T> imag(many_body_operator_generic<T> const &op) {
    return transform(op, [](monomial_t const &, T a_i) {
      using triqs::utility::imag;
      return imag(a_i);
    });
  }

  template <typename T> bool is_op_hermitian(many_body_operator_generic<T> const &op, double tolerance = 0.0) {
    return (dagger(op) - op).is_almost_zero(tolerance);
  }

  template <typename T = real_or_complex, typename... IndexTypes> many_body_operator_generic<T> c(IndexTypes... indices) {
    return many_body_operator_generic<T>::make_canonical(false, indices_t{indices...});
  }

  template <typename T = real_or_complex, typename... IndexTypes> many_body_operator_generic<T> c_dag(IndexTypes... indices) {
    return many_body_operator_generic<T>::make_canonical(true, indices_t{indices...});
  }

  template <typename T = real_or_complex, typename... IndexTypes> many_body_operator_generic<T> n(IndexTypes... indices) {
    return c_dag<T>(indices...) * c<T>(indices...);
  }

 
} // namespace triqs::operators
 
// Specialization of `std::hash` for triqs::operators::canonical_ops_t combining the `dagger` flag with each element of `indices`.
template <> struct std::hash<triqs::operators::canonical_ops_t> {
  std::size_t operator()(triqs::operators::canonical_ops_t const &c) const noexcept {
    std::size_t h = std::hash<bool>{}(c.dagger);
    for (auto const &idx : c.indices) {
      std::visit(triqs::utility::overloaded{[&](std::array<long, 3> const &a) {
                                              for (long x : a) h += std::hash<long>{}(x);
                                            },
                                            [&](auto const &v) { h += std::hash<std::decay_t<decltype(v)>>{}(v); }},
                 idx);
    }
    return h;
  }
};