legendre.hpp

Go to the documentation of this file.

// Copyright (c) 2013-2016 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2013-2016 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2019-2023 Simons Foundation
//
// 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: Philipp Dumitrescu, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "./mesh_iterator.hpp"
#include "./utils.hpp"
#include "../utility/legendre.hpp"
#include "../utility/macros.hpp"
 
#include <fmt/format.h>
#include <h5/h5.hpp>
#include <itertools/itertools.hpp>
#include <nda/nda.hpp>
 
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
 
namespace triqs::mesh {


  class C2PY_RENAME(MeshLegendre) legendre {
    public:
    using index_t = long;

    using data_index_t = long;

    class C2PY_IGNORE mesh_point_t {
      public:
      using mesh_t = legendre;

      mesh_point_t() = default;

      mesh_point_t(long n, long d, uint64_t mhash) : index_(n), data_index_(d), mesh_hash_(mhash) {}

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

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

      [[nodiscard]] uint64_t mesh_hash() const noexcept { return mesh_hash_; }
 
      private:
      long index_         = 0;
      long data_index_    = 0;
      uint64_t mesh_hash_ = 0;
    };

    legendre() = default;

    C2PY_DEPRECATED_PARAMETER_NAME(S : statistic, n_max : max_n)
    legendre(double beta, statistic_enum statistic, long max_n) : beta_(beta), stat_(statistic), N_(max_n), mesh_hash_(hash(beta, statistic, max_n)) {
      EXPECTS(beta_ > 0);
      EXPECTS(N_ >= 0);
    }

    bool operator==(legendre const &) const = default;

    [[nodiscard]] bool is_index_valid(index_t n) const noexcept { return 0 <= n and n < N_; }

    [[nodiscard]] data_index_t to_data_index(index_t n) const noexcept {
      EXPECTS(is_index_valid(n));
      return n;
    }

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

    [[nodiscard]] mesh_point_t operator[](long d) const { return {to_index(d), d, mesh_hash_}; }

    [[nodiscard]] mesh_point_t operator()(long n) const { return {n, to_data_index(n), mesh_hash_}; }

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

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

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

    [[nodiscard]] long size() const { return N_; }

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

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

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

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

    friend std::ostream &operator<<(std::ostream &sout, legendre const &m) {
      auto stat_cstr = (m.stat_ == Boson ? "Boson" : "Fermion");
      return sout << fmt::format("Legendre mesh with beta = {}, statistics = {}, N = {}", m.beta_, stat_cstr, m.size());
    }

    void serialize(auto &ar) const { ar & beta_ & stat_ & N_ & mesh_hash_; }

    void deserialize(auto &ar) { ar & beta_ & stat_ & N_ & mesh_hash_; }

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

    friend void h5_write(h5::group g, std::string const &name, legendre 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, "max_n", m.N_);
    }

    friend void h5_read(h5::group g, std::string const &name, legendre &m) {
      h5::group gr = g.open_group(name);
      h5::assert_hdf5_format(gr, m, true); // NOLINT (downcasting to base class)
 
      // for backward compatibility
      if (gr.has_key("domain")) { gr = gr.open_group("domain"); }
 
      auto beta      = h5::read<double>(gr, "beta");
      auto statistic = (h5::read<std::string>(gr, "statistic") == "F" ? Fermion : Boson);
 
      // for backward compatibility
      long max_n{};
      if (not h5::try_read(gr, "max_n", max_n)) max_n = static_cast<long>(h5::read<size_t>(gr, "n_max"));
 
      m = legendre(beta, statistic, max_n);
    }
 
    private:
    double beta_         = 1.0;
    statistic_enum stat_ = Fermion;
    long N_              = 0;
    uint64_t mesh_hash_  = 0;
  };

  inline auto evaluate(legendre const &m, auto const &f, double tau) {
    EXPECTS(m.size() > 0 and tau >= 0 and tau <= m.beta());
    utility::legendre_generator gen{};
    gen.reset(2 * tau / m.beta() - 1);
    return detail::sum_to_regular(nda::range(m.size()), [&](auto n) { return f(n) * std::sqrt(2 * n + 1) * gen.next() / m.beta(); });
  }

 
} // namespace triqs::mesh