time_pt.hpp

Go to the documentation of this file.

// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2022 Simons Foundation
// Copyright (c) 2016 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, Igor Krivenko, Olivier Parcollet, Priyanka Seth, Hugo U. R. Strand, Nils Wentzell

 
#pragma once
 
#include "./first_include.hpp"
#include "./macros.hpp"
 
#include <h5/h5.hpp>
 
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
 
namespace triqs::utility {
 
  // Forward declaration.
  struct time_segment;


  struct time_pt {
    // Friend declaration.
    friend struct time_segment;

    time_pt() = default;

    time_pt(uint64_t n_, double beta_) : n(n_), beta(beta_), val(beta_ * (static_cast<double>(n_) / static_cast<double>(Nmax))) {}

    time_pt &operator=(double v) = delete;

    auto operator<=>(time_pt const &tp) const { return n <=> tp.n; }

    bool operator==(time_pt const &tp) const { return n == tp.n; }

    inline friend time_pt operator+(time_pt const &tp1, time_pt const &tp2) {
      bool wrapped = ((Nmax - std::max(tp1.n, tp2.n)) < std::min(tp1.n, tp2.n));
      if (!wrapped)
        return {(tp1.n + tp2.n) % Nmax, tp1.beta};
      else
        return {((tp1.n + tp2.n) + 1) % Nmax, tp1.beta};
    }

    inline friend time_pt operator-(time_pt const &tp1, time_pt const &tp2) {
      uint64_t nres = (tp1.n >= tp2.n ? (tp1.n - tp2.n) % Nmax : Nmax - (tp2.n - tp1.n));
      return {nres, tp1.beta};
    }

    friend time_pt operator-(time_pt const &tp) { return {Nmax - tp.n, tp.beta}; }

    friend time_pt div_by_int(time_pt const &tp, size_t a) { return {tp.n / a, tp.beta}; }

    friend time_pt mult_by_int(time_pt const &tp, size_t a) { return {tp.n * a, tp.beta}; }

    friend size_t floor_div(time_pt const &tp1, time_pt const &tp2) { return tp1.n / tp2.n; }

    explicit operator double() const { return val; }

    friend std::ostream &operator<<(std::ostream &out, time_pt const &p) {
      return out << p.val << " [time_pt : beta = " << p.beta << " n = " << p.n << "]";
    }

    static constexpr uint64_t Nmax = std::numeric_limits<uint64_t>::max();

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

    friend void h5_write(h5::group g, std::string const &name, time_pt const &tp) {
      auto gr = g.create_group(name);
      write_hdf5_format(gr, tp); // NOLINT
      h5_write(gr, "beta", tp.beta);
      h5_write(gr, "val", tp.val);
      h5_write(gr, "n", tp.n);
    }

    friend void h5_read(h5::group g, std::string const &name, time_pt &tp) {
      auto gr = g.open_group(name);
      h5_read(gr, "beta", tp.beta);
      h5_read(gr, "val", tp.val);
      h5_read(gr, "n", tp.n);
    }
 
    private:
    uint64_t n  = 0;
    double beta = 0;
    double val  = 0;
  };

  inline double operator*(time_pt const &tp1, time_pt const &tp2) { return double(tp1) * double(tp2); }

  inline double operator/(time_pt const &tp1, time_pt const &tp2) { return double(tp1) / double(tp2); }
 
  // Arithmetic operations between time_pt and double values.
#define IMPL_OP(OP)                                                                                                                                  \
  inline double operator OP(time_pt const &x, double y) { return static_cast<double>(x) OP y; }                                                      \
  inline double operator OP(double y, time_pt const &x) { return y OP static_cast<double>(x); }
  IMPL_OP(+);
  IMPL_OP(-);
  IMPL_OP(*);
  IMPL_OP(/);
#undef IMPL_OP

  struct time_segment {
    double beta;

    time_segment(double beta_) : beta(beta_) {}

    template <typename RNG> [[nodiscard]] time_pt get_random_pt(RNG &rng, time_pt tp) const { return {rng(tp.n), beta}; }

    template <typename RNG> [[nodiscard]] time_pt get_random_pt(RNG &rng) const { return {rng(time_pt::Nmax), beta}; }

    template <typename RNG> [[nodiscard]] time_pt get_random_pt(RNG &rng, time_pt tp1, time_pt tp2) const {
      return {rng(tp2.n - tp1.n) + tp1.n, beta};
    }

    [[nodiscard]] time_pt get_upper_pt() const { return {time_pt::Nmax, beta}; }

    [[nodiscard]] time_pt get_lower_pt() const { return {0, beta}; }

    [[nodiscard]] time_pt get_epsilon() const { return {1, beta}; }

    [[nodiscard]] time_pt make_time_pt(double x) const {
      EXPECTS(0 <= x && x <= beta);
      uint64_t n = static_cast<uint64_t>(static_cast<double>(time_pt::Nmax) * std::min(1.0, (x / beta)));
      return {n, beta};
    }
  };

 
} // namespace triqs::utility