histograms.hpp

Go to the documentation of this file.

// Copyright (c) 2016-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2016-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: Philipp Dumitrescu, Igor Krivenko, Olivier Parcollet, Hugo U. R. Strand, Nils Wentzell

 
#pragma once
 
#include "../arrays.hpp"
#include "../utility/exceptions.hpp"
 
#include <h5/h5.hpp>
#include <nda/mpi.hpp>
#include <nda/nda.hpp>
 
#include <cstddef>
#include <cstdint>
#include <iosfwd>
#include <string>
#include <utility>
 
namespace triqs::stat {


  class histogram {
    public:
    histogram() = default;

    histogram(int a, int b);

    histogram(double a, double b, std::size_t nbins);

    histogram &operator<<(double x);

    void clear() {
      n_lost_pts_ = 0;
      n_data_pts_ = 0;
      data_()     = 0.0;
    }

    auto mesh_point(int n) const { return a_ + n * binsize_; }

    auto size() const { return data_.size(); }

    C2PY_PROPERTY_GET(limits) auto limits() const { return std::pair{a_, b_}; }

    C2PY_PROPERTY_GET(data) auto const &data() const { return data_; }

    C2PY_PROPERTY_GET(n_data_pts) auto n_data_pts() const { return n_data_pts_; }

    C2PY_PROPERTY_GET(n_lost_pts) auto n_lost_pts() const { return n_lost_pts_; }

    friend histogram operator+(histogram h1, histogram const &h2);

    bool operator==(histogram const &h) const = default;

    C2PY_IGNORE friend void mpi_broadcast(histogram &h, mpi::communicator c = {}, int root = 0) {
      mpi::broadcast(h.a_, c, root);
      mpi::broadcast(h.b_, c, root);
      mpi::broadcast(h.data_, c, root);
      mpi::broadcast(h.n_data_pts_, c, root);
      mpi::broadcast(h.n_lost_pts_, c, root);
      if (c.rank() != root) h.initialize();
    }

    C2PY_IGNORE friend histogram mpi_reduce(histogram const &h, mpi::communicator c = {}, int root = 0, bool all = false, MPI_Op op = MPI_SUM) {
      TRIQS_ASSERT(op == MPI_SUM);
      histogram h2(h.a_, h.b_, h.size());
      mpi::reduce_into(h.data_, h2.data_, c, root, all, MPI_SUM);
      h2.n_data_pts_ = mpi::reduce(h.n_data_pts_, c, root, all, MPI_SUM);
      h2.n_lost_pts_ = mpi::reduce(h.n_lost_pts_, c, root, all, MPI_SUM);
      return h2;
    }

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

    friend void h5_write(h5::group g, std::string const &name, histogram const &h);

    friend void h5_read(h5::group g, std::string const &name, histogram &h);

    friend std::ostream &operator<<(std::ostream &os, histogram const &h);

    void serialize(auto &ar) const { ar & a_ & b_ & n_data_pts_ & n_lost_pts_ & data_ & binsize_; }

    void deserialize(auto &ar) { ar & a_ & b_ & n_data_pts_ & n_lost_pts_ & data_ & binsize_; }

    inline friend histogram pdf(histogram const &h) {
      auto pdf = h;
      pdf.data_ /= static_cast<double>(h.n_data_pts());
      return pdf;
    }

    inline friend histogram cdf(histogram const &h) {
      auto cdf = h;
      for (int i = 1; i < h.size(); ++i) cdf.data_[i] += cdf.data_[i - 1];
      cdf.data_ /= static_cast<double>(h.n_data_pts());
      return cdf;
    }
 
    private:
    // Initialize the histogram by checking the interval and setting the bin size.
    void initialize();
 
    private:
    double a_{0.0};
    double b_{0.0};
    double binsize_{0.0};
    std::uint64_t n_data_pts_{0};
    std::uint64_t n_lost_pts_{0};
    nda::vector<double> data_{};
  };

 
} // namespace triqs::stat