dlr.hpp

Go to the documentation of this file.

// Copyright (c) 2023 Simons Foundation
// Copyright (c) 2023 Hugo U.R. Strand
//
// 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, Olivier Parcollet, Hugo U. R. Strand, Nils Wentzell

 
#pragma once
 
#include "../gf/flatten.hpp"
#include "../gf/gf.hpp"
#include "../block/block_gf.hpp"
#include "../../utility/exceptions.hpp"
#include "../../utility/tuple_tools.hpp"
 
#include "../../mesh/dlr.hpp"
#include "../../mesh/dlr_imfreq.hpp"
#include "../../mesh/dlr_imtime.hpp"
#include "../../mesh/imfreq.hpp"
#include "../../mesh/imtime.hpp"
#include "../../mesh/prod.hpp"
 
#include <array>
#include <cmath>
#include <vector>
 
namespace triqs::gfs {
 
  // Elevate the DLR mesh types into the triqs::gfs namespace.
  using mesh::dlr;
  using mesh::dlr_imfreq;
  using mesh::dlr_imtime;

 
  //-------------------------------------------------------
  // Transformation of DLR Green's functions
  // ------------------------------------------------------
 
  // Apply f to mesh component N: helper for the make_gf_dlr* transforms; recurses over product-mesh components and blocks.
  template <int N, int... Ns, typename F, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto apply_to_mesh(F const &f, G const &g) {
    using M = typename G::mesh_t;
    if constexpr (sizeof...(Ns) > 0) {
      return apply_to_mesh<Ns...>(f, apply_to_mesh<N>(f, g));
    } else if constexpr (is_block_gf_v<G>) {
      return map_block_gf([&](auto const &gbl) { return apply_to_mesh<N>(f, gbl); }, g);
    } else {
      static_assert(mesh::is_product<M>, "requires product mesh");
      auto gfl      = f(flatten_gf_2d<N>(g));
      auto out_mesh = mesh::prod{triqs::tuple::replace<N>(g.mesh(), gfl.mesh())};
      auto g_out    = gf{out_mesh, g.target_shape()};
      unflatten_gf_2d<N>(g_out, gfl);
      return g_out;
    }
  }

  template <int N = 0, int... Ns, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto make_gf_dlr(G const &g) {
    using M = typename G::mesh_t;
    if constexpr (is_block_gf_v<G>) {
      return map_block_gf([&](auto const &gbl) { return make_gf_dlr<N, Ns...>(gbl); }, g);
    } else if constexpr (mesh::is_product<M>) {
      return apply_to_mesh<N, Ns...>([&](auto const &gfl) { return make_gf_dlr(gfl); }, g);
    } else {
      static_assert(N == 0, "N must be 0 for non-product meshes");
      static_assert(nda::AnyOf<M, dlr_imtime, dlr_imfreq>, "Mesh must be dlr_imtime or dlr_imfreq");
      auto result = gf{dlr{g.mesh()}, g.target_shape()};
      if constexpr (std::is_same_v<M, dlr_imtime>)
        result.data() = result.mesh().dlr_it().vals2coefs(g.data());
      else { // dlr_imfreq
        auto beta_inv = 1. / result.mesh().beta();
        result.data() = beta_inv * result.mesh().dlr_if().vals2coefs(g.data());
      }
      return result;
    }
  }

  template <int N = 0, int... Ns, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto fit_gf_dlr(G const &g, mesh::dlr const &dlr_mesh) {
    using M = typename G::mesh_t;
    if constexpr (is_block_gf_v<G>) {
      return map_block_gf([&](auto const &gbl) { return fit_gf_dlr<N, Ns...>(gbl, dlr_mesh); }, g);
    } else if constexpr (mesh::is_product<M>) {
      return apply_to_mesh<N, Ns...>([&](auto const &gfl) { return fit_gf_dlr(gfl, dlr_mesh); }, g);
    } else {
      static_assert(N == 0, "N must be 0 for non-product meshes");
      static_assert(std::is_same_v<M, mesh::imtime>, "Input mesh must be imtime");
      auto tvals    = nda::array_adapter(std::array{g.mesh().size()}, [&](auto i) { return g.mesh()[i].value() / g.mesh().beta(); });
      auto result   = gf{dlr_mesh, g.target_shape()};
      result.data() = result.mesh().dlr_it().fitvals2coefs(make_regular(tvals), g.data());
      return result;
    }
  }

  template <int N = 0, int... Ns, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto fit_gf_dlr(G const &g, double w_max, double eps, bool symmetrize = true) {
    // Build the DLR mesh once from a representative imtime mesh and delegate to the
    // mesh overload, which handles the block / product / leaf recursion.
    auto const &imt = get_mesh<N>(g);
    auto dlr_mesh   = dlr{imt.beta(), imt.statistic(), w_max, eps, symmetrize};
    return fit_gf_dlr<N, Ns...>(g, dlr_mesh);
  }

  template <int N = 0, int... Ns, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto make_gf_dlr_imtime(G const &g) {
    using M = typename G::mesh_t;
    if constexpr (is_block_gf_v<G>) {
      return map_block_gf([&](auto const &gbl) { return make_gf_dlr_imtime<N, Ns...>(gbl); }, g);
    } else if constexpr (mesh::is_product<M>) {
      return apply_to_mesh<N, Ns...>([&](auto const &gfl) { return make_gf_dlr_imtime(gfl); }, g);
    } else if constexpr (std::is_same_v<M, dlr_imfreq>) {
      return make_gf_dlr_imtime(make_gf_dlr(g));
    } else {
      static_assert(N == 0, "N must be 0 for non-product meshes");
      static_assert(std::is_same_v<M, mesh::dlr>, "Input mesh must be dlr");
      auto result   = gf{dlr_imtime{g.mesh()}, g.target_shape()};
      result.data() = g.mesh().dlr_it().coefs2vals(g.data());
      return result;
    }
  }

  template <int N = 0, int... Ns, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto make_gf_dlr_imfreq(G const &g) {
    using M = typename G::mesh_t;
    if constexpr (is_block_gf_v<G>) {
      return map_block_gf([&](auto const &gbl) { return make_gf_dlr_imfreq<N, Ns...>(gbl); }, g);
    } else if constexpr (mesh::is_product<M>) {
      return apply_to_mesh<N, Ns...>([&](auto const &gfl) { return make_gf_dlr_imfreq(gfl); }, g);
    } else if constexpr (std::is_same_v<M, dlr_imtime>) {
      return make_gf_dlr_imfreq(make_gf_dlr(g));
    } else {
      static_assert(N == 0, "N must be 0 for non-product meshes");
      static_assert(std::is_same_v<M, mesh::dlr>, "Input mesh must be dlr");
      auto result   = gf{dlr_imfreq{g.mesh()}, g.target_shape()};
      auto beta     = result.mesh().beta();
      result.data() = beta * g.mesh().dlr_if().coefs2vals(g.data());
      return result;
    }
  }
 
  namespace detail {
    template <typename G>
      requires(MemoryGf<G> or is_block_gf_v<G>)
    mesh::imfreq const &imfreq_mesh_of(G const &g) {
      if constexpr (is_block_gf_v<G>)
        return g[0].mesh();
      else
        return g.mesh();
    }
 
    // Factored out so make_gf_dlr_imfreq and find_w_max build the dlr_imfreq mesh
    // once (the cppdlr SVD is not cheap) and reuse it across blocks / search iterations.
    template <typename G>
      requires(MemoryGf<G> or is_block_gf_v<G>)
    auto sample_on_dlr_imfreq(G const &g, dlr_imfreq const &dlr_mesh) {
      if constexpr (is_block_gf_v<G>) {
        return map_block_gf([&](auto const &gbl) { return sample_on_dlr_imfreq(gbl, dlr_mesh); }, g);
      } else {
        static_assert(std::is_same_v<typename G::mesh_t, mesh::imfreq>, "Input mesh must be imfreq");
        auto const &iw_mesh = g.mesh();
        auto [iw_lo, iw_hi] = dlr_mesh.min_max_frequencies();
        if (not iw_mesh.is_index_valid(iw_lo.n) or not iw_mesh.is_index_valid(iw_hi.n))
          TRIQS_RUNTIME_ERROR << "make_gf_dlr_imfreq: input imfreq mesh does not cover the requested DLR frequency range";
        auto result = gf{dlr_mesh, g.target_shape()};
        for (auto const &mp : dlr_mesh) result[mp] = g[iw_mesh(mp.index())];
        return result;
      }
    }
  } // namespace detail

  template <int N = 0, int... Ns, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto make_gf_dlr_imfreq(G const &g, double w_max, double eps, bool symmetrize = true) {
    using M = typename G::mesh_t;
    if constexpr (mesh::is_product<M>) {
      return apply_to_mesh<N, Ns...>([&](auto const &gfl) { return make_gf_dlr_imfreq(gfl, w_max, eps, symmetrize); }, g);
    } else {
      static_assert(N == 0, "N must be 0 for non-product meshes");
      auto const &m0 = detail::imfreq_mesh_of(g);
      auto dlr_mesh  = dlr_imfreq{m0.beta(), m0.statistic(), w_max, eps, symmetrize};
      return detail::sample_on_dlr_imfreq(g, dlr_mesh);
    }
  }

  template <int = 0, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  double find_w_max(G const &g, double eps = 1e-10, bool symmetrize = true, double w_max_init = 1.0, double w_max_max = 200.0) {
    if (w_max_init > w_max_max) TRIQS_RUNTIME_ERROR << "find_w_max: w_max_init (" << w_max_init << ") > w_max_max (" << w_max_max << ")";
 
    auto const &m0 = detail::imfreq_mesh_of(g);
 
    // NOLINTNEXTLINE(bugprone-float-loop-counter): intentional geometric (x1.5) growth of the real-valued DLR cutoff
    for (double w_max = w_max_init; w_max <= w_max_max; w_max *= 1.5) {
      auto dlr_mesh       = dlr_imfreq{m0.beta(), m0.statistic(), w_max, eps, symmetrize};
      auto [iw_lo, iw_hi] = dlr_mesh.min_max_frequencies();
      if (not m0.is_index_valid(iw_lo.n) or not m0.is_index_valid(iw_hi.n)) continue;
 
      auto g_dlr_iw  = detail::sample_on_dlr_imfreq(g, dlr_mesh);
      auto g_rec     = make_gf_imfreq(make_gf_dlr(g_dlr_iw), m0.n_iw());
      auto data_err  = [](auto const &x, auto const &y) { return double(max_element(abs(x.data() - y.data()))); };
      double max_err = 0.0;
      if constexpr (is_block_gf_v<G>) {
        for (long b = 0; b < long(g.size()); ++b) max_err = std::max(max_err, data_err(g[b], g_rec[b]));
      } else {
        max_err = data_err(g, g_rec);
      }
      if (max_err < eps) return w_max;
    }
    TRIQS_RUNTIME_ERROR << "find_w_max: no w_max <= " << w_max_max << " yields round-trip error < eps = " << eps;
  }

  template <int N = 0, int... Ns, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto make_gf_imtime(G const &g, long n_tau) {
    using M = typename G::mesh_t;
    if constexpr (is_block_gf_v<G>) {
      return map_block_gf([&](auto const &gbl) { return make_gf_imtime<N, Ns...>(gbl, n_tau); }, g);
    } else if constexpr (mesh::is_product<M>) {
      return apply_to_mesh<N, Ns...>([&](auto const &gfl) { return make_gf_imtime(gfl, n_tau); }, g);
    } else if constexpr (nda::AnyOf<M, dlr_imtime, dlr_imfreq>) {
      return make_gf_imtime(make_gf_dlr(g), n_tau);
    } else { // M == dlr
      static_assert(N == 0, "N must be 0 for non-product meshes");
      static_assert(std::is_same_v<M, dlr>, "Input mesh must be dlr");
      auto result = gf{mesh::imtime{g.mesh().beta(), g.mesh().statistic(), n_tau}, g.target_shape()};
      for (auto tau : result.mesh()) result[tau] = g(tau.value());
      return result;
    }
  }

  template <int N = 0, int... Ns, typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto make_gf_imfreq(G const &g, long n_iw) {
    using M = typename G::mesh_t;
    if constexpr (is_block_gf_v<G>) {
      return map_block_gf([&](auto const &gbl) { return make_gf_imfreq<N, Ns...>(gbl, n_iw); }, g);
    } else if constexpr (mesh::is_product<M>) {
      return apply_to_mesh<N, Ns...>([&](auto const &gfl) { return make_gf_imfreq(gfl, n_iw); }, g);
    } else if constexpr (nda::AnyOf<M, dlr_imtime, dlr_imfreq>) {
      return make_gf_imfreq(make_gf_dlr(g), n_iw);
    } else { // M == dlr
      static_assert(N == 0, "N must be 0 for non-product meshes");
      static_assert(std::is_same_v<M, dlr>, "Input mesh must be dlr");
      auto result = gf{mesh::imfreq{g.mesh().beta(), g.mesh().statistic(), n_iw}, g.target_shape()};
      for (auto w : result.mesh()) result[w] = g(w.value());
      return result;
    }
  }

  template <typename G>
    requires(MemoryGf<G> or is_block_gf_v<G>)
  auto tau_L2_norm(G const &g) {
    using M = typename G::mesh_t;
    static_assert(nda::AnyOf<M, dlr, dlr_imfreq, dlr_imtime>, "Input mesh must be one of dlr, dlr_imfreq, dlr_imtime");
    if constexpr (is_block_gf_v<G>) {
      // Manual loop with a concrete element type. map_block_gf would deduce
      // std::vector<std::invoke_result_t<lambda, gf>> (the lambda is local),
      // and a function-scoped `using elem_t = ...` alias would also leak its
      // spelling -- in both cases clair-c2py fails to render the type.
      constexpr int rank = std::decay_t<G>::g_t::target_t::rank;
      if constexpr (rank == 0) {
        std::vector<double> result;
        result.reserve(g.data().size());
        for (auto const &gbl : g.data()) result.push_back(tau_L2_norm(gbl));
        return result;
      } else {
        std::vector<nda::array<double, rank>> result;
        result.reserve(g.data().size());
        for (auto const &gbl : g.data()) result.push_back(tau_L2_norm(gbl));
        return result;
      }
    } else if constexpr (nda::AnyOf<M, dlr_imtime, dlr_imfreq>) {
      return tau_L2_norm(make_gf_dlr(g));
    } else { // M == dlr
      if constexpr (G::target_rank == 0) {
        // use abs here to avoid sqrt of negative numbers, which can happen due to numerical errors
        return std::sqrt(std::abs(g.mesh().dlr_it().innerprod(g.data(), g.data())));
      } else {
        auto result = nda::array<double, G::target_rank>(g.target_shape());
        nda::for_each(g.target_shape(), [&](auto... is) {
          auto dat      = g.data()(nda::range::all, is...);
          result(is...) = std::sqrt(std::abs(g.mesh().dlr_it().innerprod(dat, dat)));
        });
 
        return result;
      }
    }
  }

 
} // namespace triqs::gfs