linear.hpp

Go to the documentation of this file.

// Copyright (c) 2013-2015 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2013-2015 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, Michel Ferrero, Henri Menke, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "../utils.hpp"
#include "../mesh_iterator.hpp"
#include "../../utility/macros.hpp"
 
#include <h5/h5.hpp>
#include <itertools/itertools.hpp>
#include <nda/nda.hpp>
 
#include <algorithm>
#include <concepts>
#include <cstdint>
#include <limits>
#include <string>
#include <type_traits>
#include <utility>
 
namespace triqs::mesh::detail {

  template <typename M, typename T>
    requires std::totally_ordered<T>
  class linear {
    public:
    using value_t = T;

    using index_t = long;

    using data_index_t = long;

    class mesh_point_t {
      public:
      using mesh_t = M;

      mesh_point_t() = default;

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

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

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

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

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

      operator value_t() const { return value_; }
 
#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;
      value_t value_      = {};
    };

    linear() = default;

    linear(value_t a, value_t b, long N) : N_(N), a_(a), b_(b) {
      EXPECTS(N_ >= 0);
      if (N_ == 1) {
        EXPECTS(a_ == b_);
        delta_     = 0;
        delta_inv_ = std::numeric_limits<double>::infinity();
      } else if (N_ > 1) {
        EXPECTS(a_ < b_);
        delta_     = (b_ - a_) / (N_ - 1);
        delta_inv_ = 1 / delta_;
      }
    }

    C2PY_IGNORE bool operator==(linear const &rhs) const = default;

    [[nodiscard]] bool is_index_valid(index_t n) const noexcept { return 0 <= n and n < N_; }
 
    private:
    // Check if a value is valid.
    [[nodiscard]] bool is_value_valid(value_t m) const noexcept { return N_ > 0 and a_ <= m and m <= b_; }
 
    public:
    [[nodiscard]] data_index_t to_data_index(index_t n) const noexcept {
      EXPECTS(is_index_valid(n));
      return n;
    }

    [[nodiscard]] C2PY_IGNORE index_t to_data_index(closest_mesh_point_t<value_t> 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<value_t> const &cmp) const noexcept {
      EXPECTS(is_value_valid(cmp.value));
      return (N_ == 1 ? 0 : static_cast<index_t>((cmp.value - a_) * delta_inv_ + 0.5));
    }

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

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

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

    [[nodiscard]] value_t to_value(index_t n) const noexcept {
      EXPECTS(is_index_valid(n));
      if (N_ == 1) return a_;
      auto const w = static_cast<double>(n) / (N_ - 1);
      return a_ * (1 - w) + b_ * w;
    }

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

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

    [[nodiscard]] C2PY_PROPERTY_GET(delta) value_t delta() const noexcept { return delta_; }

    [[nodiscard]] value_t delta_inv() const noexcept { return delta_inv_; }

    [[nodiscard]] long first_index() const { return 0; }

    [[nodiscard]] long last_index() const { return N_ - 1; }

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

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

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

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

    void serialize(auto &ar) const { ar & N_ & a_ & b_ & delta_ & delta_inv_ & mesh_hash_; }

    void deserialize(auto &ar) { ar & N_ & a_ & b_ & delta_ & delta_inv_ & mesh_hash_; }
 
    protected:
    void h5_write_impl(h5::group g, std::string const &name, const char *format) const {
      h5::group gr = g.create_group(name);
      h5::write_hdf5_format_as_string(gr, format); // NOLINT (downcasting to base class)
      h5::write(gr, "min", a_);
      h5::write(gr, "max", b_);
      h5::write(gr, "size", N_);
    }

    void h5_read_impl(h5::group g, std::string const &name, const char *exp_format) {
      h5::group gr = g.open_group(name);
      h5::assert_hdf5_format_as_string(gr, exp_format, true); // NOLINT (downcasting to base class)
      auto a = h5::read<value_t>(gr, "min");
      auto b = h5::read<value_t>(gr, "max");
      auto N = h5::read<long>(gr, "size");
      *this  = linear(a, b, N);
    }
 
    public:
    auto evaluate(auto const &f, double x) const {
      EXPECTS(is_value_valid(x) and N_ > 1);
      auto const n_dbl = (x - a_) * delta_inv_;
      auto const n     = std::clamp<index_t>(static_cast<index_t>(n_dbl), 0, N_ - 2);
      auto const w     = std::min(n_dbl - n, 1.0);
      return (1 - w) * f(n) + w * f(n + 1);
    }
 
    protected:
    long N_{0};
    value_t a_{0}, b_{0};
    value_t delta_{0}, delta_inv_{0};
    uint64_t mesh_hash_{hash(N_, a_, b_)};
  };
 
} // namespace triqs::mesh::detail