dlr.hpp

Go to the documentation of this file.

// Copyright (c) 2023 Simons Foundation
// Copyright (c) 2023 Hugo U.R. Strand
//
// 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: Alexander Hampel, Olivier Parcollet, Hugo U.R. Strand, Nils Wentzell

 
#pragma once
 
#include "./matsubara_freq.hpp"
#include "./mesh_iterator.hpp"
#include "./utils.hpp"
#include "../utility/macros.hpp"
 
#include <cppdlr/cppdlr.hpp>
#include <h5/h5.hpp>
#include <fmt/format.h>
#include <itertools/itertools.hpp>
#include <nda/nda.hpp>
 
#include <cstdint>
#include <iostream>
#include <memory>
#include <string>
#include <string_view>
#include <utility>
 
namespace triqs::mesh {
 
  namespace detail {
 
    // Struct that combines the DLR frequencies, DLR imaginary time operations and DLR Matsubara frequency operations.
    struct dlr_ops {
      nda::vector<double> freq;
      cppdlr::imtime_ops imt;
      cppdlr::imfreq_ops imf;
    };
 
  } // namespace detail
 
  // Forward declarations.
  class dlr_imtime;
  class dlr_imfreq;


  class C2PY_RENAME(MeshDLR) dlr {
    public:
    using value_t = double;

    using index_t = long;

    using data_index_t = long;

    class C2PY_IGNORE mesh_point_t {
      public:
      using mesh_t = dlr;

      mesh_point_t() = default;

      mesh_point_t(long l, long d, uint64_t mhash, double w_l) : index_(l), data_index_(d), mesh_hash_(mhash), value_(w_l) {}

      [[nodiscard]] long index() const { return index_; }

      [[nodiscard]] long data_index() const { return data_index_; }

      [[nodiscard]] double value() const { return value_; }

      [[nodiscard]] uint64_t mesh_hash() const noexcept { return mesh_hash_; }

      operator double() const { return value_; }
 
      private:
      long index_         = 0;
      long data_index_    = 0;
      uint64_t mesh_hash_ = 0;
      double value_       = {};
    };
 
    private:
    // Construct a DLR mesh with a given set of DLR frequencies.
    dlr(double beta, statistic_enum statistic, double w_max, double eps, bool symmetrize, nda::vector<double> const &dlr_freq)
       : dlr(beta, statistic, w_max, eps, symmetrize,
             detail::dlr_ops{.freq = dlr_freq,
                             .imt  = {w_max * beta, dlr_freq, symmetrize},
                             .imf  = {w_max * beta, dlr_freq, static_cast<cppdlr::statistic_t>(statistic), symmetrize}}) {}
 
    // Construct a DLR mesh with given DLR operations.
    dlr(double beta, statistic_enum statistic, double w_max, double eps, bool symmetrize, detail::dlr_ops ops)
       : beta_(beta),
         stat_(statistic),
         w_max_(w_max),
         eps_(eps),
         symmetrize_(symmetrize),
         mesh_hash_(hash(beta, statistic, w_max, eps, symmetrize, hash_bytes(ops.imf.get_ifnodes()), std::string_view{"dlr"})),
         dlr_{std::make_shared<detail::dlr_ops>(std::move(ops))} {}
 
    public:
    dlr() = default;

    dlr(double beta, statistic_enum statistic, double w_max, double eps, bool symmetrize = true)
       : dlr(beta, statistic, w_max, eps, symmetrize, cppdlr::build_dlr_rf(w_max * beta, eps, symmetrize)) {}

    template <nda::AnyOf<dlr_imtime, dlr_imfreq> M>
    explicit dlr(M const &m)
       : beta_(m.beta_),
         stat_(m.stat_),
         w_max_(m.w_max_),
         eps_(m.eps_),
         symmetrize_(m.symmetrize_),
         mesh_hash_(hash(beta_, stat_, w_max_, eps_, symmetrize_, hash_bytes(m.dlr_->imf.get_ifnodes()), std::string_view{"dlr"})),
         dlr_(m.dlr_) {}

    bool operator==(dlr const &m) const { return mesh_hash_ == m.mesh_hash_; }

    [[nodiscard]] bool is_index_valid(long l) const noexcept { return 0 <= l and l < size(); }

    [[nodiscard]] long to_data_index(long l) const noexcept {
      EXPECTS(is_index_valid(l));
      return l;
    }

    [[nodiscard]] long to_index(long d) const noexcept {
      EXPECTS(is_index_valid(d));
      return d;
    }

    [[nodiscard]] mesh_point_t operator[](long d) const { return (*this)(d); }

    [[nodiscard]] mesh_point_t operator()(long l) const { return {l, l, mesh_hash_, to_value(l)}; }

    [[nodiscard]] double to_value(long l) const noexcept {
      EXPECTS(is_index_valid(l));
      return (dlr_->freq)[l];
    }

    [[nodiscard]] C2PY_PROPERTY_GET(beta) double beta() const noexcept { return beta_; }

    [[nodiscard]] C2PY_PROPERTY_GET(statistic) statistic_enum statistic() const noexcept { return stat_; }

    [[nodiscard]] C2PY_PROPERTY_GET(w_max) double w_max() const noexcept { return w_max_; }

    [[nodiscard]] C2PY_PROPERTY_GET(eps) double eps() const noexcept { return eps_; }

    [[nodiscard]] C2PY_PROPERTY_GET(symmetrize) bool symmetrize() const noexcept { return symmetrize_; }

    [[nodiscard]] C2PY_PROPERTY_GET(dlr_freq) auto const &dlr_freq() const { return dlr_->freq; }

    [[nodiscard]] C2PY_IGNORE auto const &dlr_it() const { return dlr_->imt; }

    [[nodiscard]] C2PY_IGNORE auto const &dlr_if() const { return dlr_->imf; }

    [[nodiscard]] C2PY_PROPERTY_GET(mesh_hash) uint64_t mesh_hash() const noexcept { return mesh_hash_; }

    [[nodiscard]] long size() const noexcept { return (dlr_ ? dlr_->freq.size() : 0); }

    [[nodiscard]] auto begin() const { return mesh_iterator<dlr>{.mesh_ptr = this, .data_index = 0}; }

    [[nodiscard]] auto cbegin() const { return begin(); }

    [[nodiscard]] auto end() const { return mesh_iterator<dlr>{.mesh_ptr = this, .data_index = size()}; }

    [[nodiscard]] auto cend() const { return end(); }

    friend std::ostream &operator<<(std::ostream &sout, dlr const &m) {
      auto stat_cstr = (m.stat_ == Boson ? "Boson" : "Fermion");
      return sout << fmt::format("DLR coefficient mesh of size {} with beta = {}, statistics = {}, w_max = {}, eps = {}, symmetrized = {}", m.size(),
                                 m.beta_, stat_cstr, m.w_max_, m.eps_, m.symmetrize_);
    }

    void serialize(auto &ar) const {
      EXPECTS(dlr_);
      ar & beta_ & stat_ & w_max_ & eps_ & symmetrize_ & mesh_hash_ & dlr_->freq;
      dlr_->imt.serialize(ar);
      dlr_->imf.serialize(ar);
    }

    void deserialize(auto &ar) {
      nda::vector<double> freq;
      cppdlr::imtime_ops imt;
      cppdlr::imfreq_ops imf;
      ar & beta_ & stat_ & w_max_ & eps_ & symmetrize_ & mesh_hash_ & freq;
      imt.deserialize(ar);
      imf.deserialize(ar);
      dlr_ = std::make_shared<detail::dlr_ops>(freq, imt, imf);
    }

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

    friend void h5_write(h5::group g, std::string const &name, dlr const &m) {
      h5::group gr = g.create_group(name);
      h5::write_hdf5_format(gr, m); // NOLINT (downcasting to base class)
      h5::write(gr, "beta", m.beta_);
      h5::write(gr, "statistic", (m.stat_ == Fermion ? "F" : "B"));
      h5::write(gr, "w_max", m.w_max_);
      h5::write(gr, "eps", m.eps_);
      h5::write(gr, "symmetrize", m.symmetrize_);
      h5::write(gr, "dlr_freq", m.dlr_freq());
      h5::write(gr, "dlr_it", m.dlr_it());
      h5::write(gr, "dlr_if", m.dlr_if());
    }

    friend void h5_read(h5::group g, std::string const &name, dlr &m) {
      h5::group gr = g.open_group(name);
      h5::assert_hdf5_format(gr, m, true); // NOLINT (downcasting to base class)
      auto b       = h5::read<double>(gr, "beta");
      auto stat    = (h5::read<std::string>(gr, "statistic") == "F" ? Fermion : Boson);
      auto wmax    = h5::read<double>(gr, "w_max");
      auto epsilon = h5::read<double>(gr, "eps");
      bool sym     = false;
      h5::try_read(gr, "symmetrize", sym);
      auto freq = h5::read<nda::vector<double>>(gr, "dlr_freq");
      auto imt  = h5::read<cppdlr::imtime_ops>(gr, "dlr_it");
      auto imf  = h5::read<cppdlr::imfreq_ops>(gr, "dlr_if");
      m         = dlr(b, stat, wmax, epsilon, sym, {.freq = freq, .imt = imt, .imf = imf});
    }
 
    // Friend declarations.
    friend class dlr_imtime;
    friend class dlr_imfreq;
 
    private:
    double beta_                                = 1.0;
    statistic_enum stat_                        = Fermion;
    double w_max_                               = 0.0;
    double eps_                                 = 1e-10;
    bool symmetrize_                            = false;
    uint64_t mesh_hash_                         = 0;
    std::shared_ptr<const detail::dlr_ops> dlr_ = {};
  };

  auto evaluate(dlr const &m, auto const &f, double tau) {
    EXPECTS(m.size() > 0);
    EXPECTS(tau >= 0 and tau <= m.beta());
    return detail::sum_to_regular(nda::range(m.size()), [&](auto l) { return f(l) * cppdlr::k_it(tau / m.beta(), m.dlr_freq()[l]); });
  }

  auto evaluate(dlr const &m, auto const &f, matsubara_freq const &iw) {
    EXPECTS(m.size() > 0);
    return detail::sum_to_regular(nda::range(m.size()),
                                  [&](auto l) { return f(l) * cppdlr::k_if(iw.n, m.dlr_freq()[l], (cppdlr::statistic_t)iw.statistic) * m.beta(); });
  }

 
} // namespace triqs::mesh