prod.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-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, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "./concepts.hpp"
#include "./utils.hpp"
#include "../utility/tuple_tools.hpp"
 
#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
 
// to silence some clangd warning and make gf independent of triqs::lattice. cf l 64. backward compat check.
namespace triqs::lattice {
  class brillouin_zone;
}
 
namespace triqs::mesh {

 
  // Forward declaration.
  template <Mesh... Ms>
    requires(((n_variables<Ms> == 1) and ...) and (not(std::is_reference_v<Ms> or ...)))
  class prod;
 
  // Specialize some traits for the product mesh type.
  template <Mesh... Ms> static constexpr bool is_product<prod<Ms...>> = true;
  template <Mesh... Ms> static constexpr int n_variables<prod<Ms...>> = sizeof...(Ms);

  template <Mesh... Ms> class prod_mesh_point : public std::tuple<typename Ms::mesh_point_t...> {
    public:
    using mesh_t = prod<Ms...>;

    using index_t = typename mesh_t::index_t;

    using data_index_t = typename mesh_t::data_index_t;

    using tuple_t = std::tuple<typename Ms::mesh_point_t...>;

    prod_mesh_point() = default;

    prod_mesh_point(tuple_t mpt) : tuple_t(std::move(mpt)) {}

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

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

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

    tuple_t const &as_tuple() const { return *this; }
 
    private:
    index_t index_           = std::apply([](auto &...ms) { return std::make_tuple(ms.index()...); }, as_tuple());
    data_index_t data_index_ = std::apply([](auto &...ms) { return std::make_tuple(ms.data_index()...); }, as_tuple());
    uint64_t mesh_hash_      = std::apply([](auto &...ms) { return (ms.mesh_hash() + ...); }, as_tuple());
  };

  template <Mesh... Ms>
    requires(((n_variables<Ms> == 1) and ...) and (not(std::is_reference_v<Ms> or ...)))
  class prod : public std::tuple<Ms...> {
    public:
    using index_t = std::tuple<typename Ms::index_t...>;

    using data_index_t = std::tuple<typename Ms::data_index_t...>;

    using mesh_point_t = prod_mesh_point<Ms...>;

    using m_tuple_t = std::tuple<Ms...>;

    prod() = default;

    prod(Ms const &...ms)
      requires(sizeof...(Ms) > 0)
       : m_tuple_t{ms...}, mesh_hash_((ms.mesh_hash() + ...)) {}

    template <typename... U>
    prod(std::tuple<U...> const &mt) : m_tuple_t{mt}, mesh_hash_(std::apply([](auto &...m) { return (m.mesh_hash() + ...); }, mt)) {}

    bool operator==(prod const &m) const { return as_tuple() == m.as_tuple(); }

    template <typename... Args> [[nodiscard]] bool is_index_valid(Args const &...ns) const {
      return triqs::tuple::fold([](auto &m, auto &n, bool res) { return res && (m.is_index_valid(n)); }, as_tuple(), std::tie(ns...), true);
    }

    template <typename... Args> [[nodiscard]] bool is_index_valid(index_t const &n) const {
      return triqs::tuple::fold([](auto &m, auto &n, bool res) { return res && (m.is_index_valid(n)); }, as_tuple(), n, true);
    }

    [[nodiscard]] data_index_t to_data_index(index_t const &n) const {
      auto l = [](auto const &m, auto const &n_i) { return m.to_data_index(n_i); };
      return triqs::tuple::map_on_zip(l, *this, n);
    }

    [[nodiscard]] index_t to_index(data_index_t const &d) const {
      auto l = [](auto const &m, auto const &d_i) { return m.to_index(d_i); };
      return triqs::tuple::map_on_zip(l, *this, d);
    }

    [[nodiscard]] mesh_point_t operator[](data_index_t const &d) const {
      auto l = [](auto const &m, auto const &d_i) { return m[d_i]; };
      return triqs::tuple::map_on_zip(l, *this, d);
    }

    [[nodiscard]] mesh_point_t operator()(index_t const &n) const {
      auto l = [](auto const &m, auto const &n_i) { return m(n_i); };
      return triqs::tuple::map_on_zip(l, *this, n);
    }

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

    [[nodiscard]] long size() const {
      return triqs::tuple::fold([](auto const &m, size_t n) { return n * m.size(); }, as_tuple(), 1);
    }

    [[nodiscard]] auto size_of_components() const {
      std::array<long, std::tuple_size_v<m_tuple_t>> res;
      auto l = [&res](int i, auto const &m) mutable { res[i] = m.size(); };
      triqs::tuple::for_each_enumerate(as_tuple(), l);
      return res;
    }
 
    // Get the underlying tuple of individual meshes, i.e. \f$ (M_1, \dots, M_k) \f$.
    [[nodiscard]] m_tuple_t const &components() const { return *this; }

    [[nodiscard]] m_tuple_t &components() { return *this; }

    [[nodiscard]] m_tuple_t const &as_tuple() const { return *this; }

    [[nodiscard]] m_tuple_t &as_tuple() { return *this; }
 
    private:
    // Get a lazy range that generates the mesh points.
    [[nodiscard]] auto r_() const {
      auto f             = [](auto &&...ms) { return itertools::product(ms...); };
      auto to_mesh_point = [](auto &&n_mp) { return mesh_point_t{n_mp}; };
      return itertools::transform(std::apply(f, components()), to_mesh_point);
    }
 
    public:
    [[nodiscard]] auto begin() const { return r_().begin(); }

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

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

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

    friend std::ostream &operator<<(std::ostream &sout, prod const &m) {
      sout << "Product Mesh";
      triqs::tuple::for_each(m.as_tuple(), [&sout](auto &mesh) { sout << "\n  -- " << mesh; });
      return sout;
    }

    void serialize(auto &ar) const { ar &as_tuple(); }

    void deserialize(auto &ar) { ar &as_tuple(); }

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

    friend void h5_write(h5::group g, std::string name, prod const &m) {
      h5::group gr = g.create_group(name);
      h5::write_hdf5_format(gr, m);
      auto l = [gr](int N, auto const &mc) { h5::write(gr, "MeshComponent" + std::to_string(N), mc); };
      triqs::tuple::for_each_enumerate(m.components(), l);
    }

    friend void h5_read(h5::group g, std::string name, prod &m) {
      h5::group gr = g.open_group(name);
      h5::assert_hdf5_format(gr, m, true);
      auto l = [gr](int N, auto &mc) { h5::read(gr, "MeshComponent" + std::to_string(N), mc); };
      triqs::tuple::for_each_enumerate(m.components(), l);
    }
 
    private:
    uint64_t mesh_hash_ = 0;
  };
 
  // Class template argument deduction rules (CTAD).
  template <typename M1, typename M2, typename... Ms> prod(M1, M2, Ms...) -> prod<M1, M2, Ms...>;
 
  template <typename M1, typename M2, typename... Ms>
  prod(std::tuple<M1, M2, Ms...>) -> prod<std::decay_t<M1>, std::decay_t<M2>, std::decay_t<Ms>...>;

  template <Mesh... Ls, Mesh... Rs> auto operator*(prod<Ls...> const &m, prod<Rs...> const &n) {
    return prod<Ls..., Rs...>{std::tuple_cat(m.components(), n.components())};
  }

  template <Mesh L, Mesh... Rs> auto operator*(L const &n, prod<Rs...> const &m) {
    return prod<L, Rs...>{std::tuple_cat(std::make_tuple(n), m.components())};
  }

  template <Mesh... Ls, Mesh R> auto operator*(prod<Ls...> const &m, R const &n) {
    return prod<Ls..., R>{std::tuple_cat(m.components(), std::make_tuple(n))};
  }

  template <Mesh L, Mesh R> auto operator*(L const &m, R const &n) { return prod<L, R>{m, n}; }

 
} // namespace triqs::mesh

template <typename... Ms> struct std::tuple_size<triqs::mesh::prod<Ms...>> : public std::integral_constant<size_t, sizeof...(Ms)> {};

template <typename... Ms> struct std::tuple_size<triqs::mesh::prod_mesh_point<Ms...>> : public std::integral_constant<size_t, sizeof...(Ms)> {};

template <size_t N, typename... Ms> struct std::tuple_element<N, triqs::mesh::prod<Ms...>> : public std::tuple_element<N, std::tuple<Ms...>> {};

template <size_t N, typename... Ms>
struct std::tuple_element<N, triqs::mesh::prod_mesh_point<Ms...>> : public std::tuple_element<N, std::tuple<typename Ms::mesh_point_t...>> {};