refreq_log.hpp

Go to the documentation of this file.

// Copyright (c) 2025 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: Nils Wentzell

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


  class C2PY_RENAME(MeshReFreqLog) refreq_log {
    public:
    using value_t = double;

    using index_t = long;

    using data_index_t = long;

    class C2PY_IGNORE mesh_point_t {
      public:
      using mesh_t = refreq_log;

      mesh_point_t() = default;

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

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

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

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

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

      operator double() const { return value_; } // NOLINT (implicit conversion intended)
 
      // Arithmetic operations with scalars.
#define IMPL_OP(OP)                                                                                                                                  \
  template <typename U> friend auto operator OP(mesh_point_t const &mp, U &&y) { return mp.value() OP std::forward<U>(y); }                          \
  template <typename U>                                                                                                                              \
    requires(not std::is_same_v<std::decay_t<U>, mesh_point_t>)                                                                                      \
  friend auto operator OP(U &&x, mesh_point_t const &mp) {                                                                                           \
    return std::forward<U>(x) OP mp.value();                                                                                                         \
  }
      IMPL_OP(+)
      IMPL_OP(-)
      IMPL_OP(*)
      IMPL_OP(/)
#undef IMPL_OP
 
      private:
      long index_         = 0;
      long data_index_    = 0;
      uint64_t mesh_hash_ = 0;
      double value_       = 0.0;
    };

    refreq_log() = default;

    refreq_log(double eps, double w_max, double ratio) : eps_(eps), w_max_(w_max), ratio_(ratio), mesh_hash_(hash(eps, w_max, ratio)) {
      EXPECTS(eps_ > 0);
      EXPECTS(w_max_ >= eps_);
      EXPECTS(ratio_ > 1.0);
      init_points();
    }

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

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

    [[nodiscard]] bool is_value_valid(double w) const noexcept { return !pts_.empty() and pts_.front() <= w and w <= pts_.back(); }

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

    [[nodiscard]] C2PY_IGNORE data_index_t to_data_index(closest_mesh_point_t<double> const &cmp) const noexcept {
      return to_data_index(to_index(cmp));
    }

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

    [[nodiscard]] C2PY_IGNORE index_t to_index(closest_mesh_point_t<double> const &cmp) const noexcept {
      EXPECTS(is_value_valid(cmp.value));
 
      auto itr_r = std::ranges::lower_bound(pts_, cmp.value);
      long i_r   = itr_r - pts_.begin();
 
      if (i_r == 0) { return 0; }
      if (i_r == size()) { return size() - 1; }
 
      long i_l = i_r - 1;
      if (std::abs(cmp.value - pts_[i_l]) < std::abs(cmp.value - pts_[i_r]))
        return i_l;
      else
        return i_r;
    }

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

    [[nodiscard]] C2PY_IGNORE mesh_point_t operator[](closest_mesh_point_t<double> const &cmp) const noexcept { return (*this)[to_data_index(cmp)]; }

    [[nodiscard]] mesh_point_t operator()(index_t n) const noexcept {
      EXPECTS(is_index_valid(n));
      return {n, n, mesh_hash_, pts_[n]};
    }

    [[nodiscard]] double to_value(index_t n) const noexcept {
      EXPECTS(is_index_valid(n));
      return pts_[n];
    }

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

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

    [[nodiscard]] C2PY_PROPERTY_GET(ratio) double ratio() const noexcept { return ratio_; }

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

    [[nodiscard]] long size() const noexcept { return static_cast<long>(pts_.size()); }

    [[nodiscard]] C2PY_PROPERTY_GET(points) std::vector<double> const &points() const noexcept { return pts_; }

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

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

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

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

    friend std::ostream &operator<<(std::ostream &sout, refreq_log const &m) {
      return sout << fmt::format("Logarithmic real frequency mesh with eps = {}, w_max = {}, ratio = {}, N = {}", m.eps_, m.w_max_, m.ratio_,
                                 m.size());
    }

    void serialize(auto &ar) const { ar & eps_ & w_max_ & ratio_ & mesh_hash_; }

    void deserialize(auto &ar) {
      ar & eps_ & w_max_ & ratio_ & mesh_hash_;
      if (w_max_ > 0) init_points();
    }

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

    friend void h5_write(h5::group g, std::string const &name, refreq_log const &m) {
      h5::group gr = g.create_group(name);
      h5::write_hdf5_format(gr, m); // NOLINT (downcasting to base class)
      h5::write(gr, "eps", m.eps_);
      h5::write(gr, "w_max", m.w_max_);
      h5::write(gr, "ratio", m.ratio_);
    }

    friend void h5_read(h5::group g, std::string const &name, refreq_log &m) {
      h5::group gr = g.open_group(name);
      h5::assert_hdf5_format(gr, m, true); // NOLINT (downcasting to base class)
      auto e = h5::read<double>(gr, "eps");
      auto w = h5::read<double>(gr, "w_max");
      auto r = h5::read<double>(gr, "ratio");
      m      = refreq_log(e, w, r);
    }
 
    private:
    void init_points() {
      pts_.clear();
      for (double w = w_max_; w >= eps_; w /= ratio_) {
        pts_.push_back(w);
        pts_.push_back(-w);
      }
      std::ranges::sort(pts_);
    }
 
    double eps_         = 0.0;
    double w_max_       = 0.0;
    double ratio_       = 0.0;
    uint64_t mesh_hash_ = 0;
    std::vector<double> pts_;
  };

  inline auto evaluate(refreq_log const &m, auto const &f, double w) {
    EXPECTS(m.is_value_valid(w));
 
    auto const &pts = m.points();
    auto itr_r      = std::ranges::lower_bound(pts, w);
    long i_r        = itr_r - pts.begin();
 
    if (i_r == 0) { return f(0); }
    if (i_r == m.size()) { return f(m.size() - 1); }
 
    long i_l = i_r - 1;
 
    double w_l = pts[i_l];
    double w_r = pts[i_r];
    double del = w_r - w_l;
 
    double wt_r = (w - w_l) / del;
    double wt_l = (w_r - w) / del;
 
    return f(i_l) * wt_l + f(i_r) * wt_r;
  }

 
} // namespace triqs::mesh