jackknife.hpp

Go to the documentation of this file.

// Copyright (c) 2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2021 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 "./concepts.hpp"
#include "./mean_error.hpp"
 
#include <mpi/mpi.hpp>
#include <nda/nda.hpp>
 
#include <cstddef>
#include <optional>
#include <ranges>
#include <stdexcept>
#include <tuple>
#include <type_traits>
#include <utility>
 
namespace triqs::stat {

 
  namespace detail {
 
    // Implementation of the jackknife resampling.
    auto jackknife_impl(std::optional<mpi::communicator> comm, auto &&f, StatCompatibleRange auto &&...rgs) { // NOLINT (ranges need not be forwarded)
      // check that the ranges all have the same size
      auto rgs_tup = std::tie(rgs...);
      auto size    = std::ranges::size(std::get<0>(rgs_tup));
      if (((std::ranges::size(rgs) != size) or ...)) throw std::runtime_error("Error in jackknife: Ranges must have the same size.");
 
      // all_reduce the size of the ranges if a communicator is provided
      if (comm) mpi::all_reduce_in_place(size, *comm);
 
      // get the mean and the jackknife samples of each range and store them in a tuple
      auto [mean_tup, jk_tup] = [&comm, size, &rgs_tup]<std::size_t... Is>(std::index_sequence<Is...>) {
        auto mtup = std::make_tuple(mean_mpi(comm, std::get<Is>(rgs_tup))...);
        auto jtup = std::make_tuple(std::ranges::transform_view(
           std::get<Is>(rgs_tup), [m = std::get<Is>(mtup), size](auto const &x) { return (size * m - x) / (size - 1); })...);
        return std::make_pair(mtup, jtup);
      }(std::make_index_sequence<sizeof...(rgs)>{});
 
      // get a zipped view of the jackknife samples from all the ranges
      auto zipped = std::apply([](auto &&...jks) { return std::ranges::zip_view(jks...); }, jk_tup);
 
      // compute the jackknife mean and error (we need regular types to avoid dangling references)
      auto f_wrapped      = [f](auto const &...args) { return nda::make_regular(f(args...)); };
      auto [jk_f, jk_err] = mean_and_err_mpi<error_tag::jk_err>(
         comm, std::ranges::transform_view(zipped, [f_wrapped](auto const &tup) { return std::apply(f_wrapped, tup); }));
 
      // compute the naive estimate of f(...) and the bias corrected estimate
      auto naive_f = nda::make_regular(std::apply(f, mean_tup));
      auto corr_f  = nda::make_regular(size * naive_f - (size - 1) * jk_f);
 
      return std::make_tuple(corr_f, jk_err, jk_f, naive_f);
    }
 
  } // namespace detail

  template <typename F, StatCompatibleRange R, StatCompatibleRange... Rs>
    requires(not std::same_as<std::remove_cvref_t<F>, mpi::communicator>)
  auto jackknife(F &&f, R &&rg, Rs &&...rgs) { // NOLINT (ranges need not be forwarded)
    return detail::jackknife_impl(std::nullopt, std::forward<F>(f), rg, rgs...);
  }

  template <typename F, StatCompatibleRange R, StatCompatibleRange... Rs>
  auto jackknife_mpi(mpi::communicator comm, F &&f, R &&rg, Rs &&...rgs) { // NOLINT (ranges need not be forwarded)
    return detail::jackknife_impl(comm, std::forward<F>(f), rg, rgs...);
  }

 
} // namespace triqs::stat