types.hpp

// Copyright (c) 2017--present, The Simons Foundation
// This file is part of TRIQS/ctint and is licensed under the terms of GPLv3 or later.
// SPDX-License-Identifier: GPL-3.0-or-later
// See LICENSE in the root of this distribution for details.
 
 
#pragma once
 
#include <triqs/gfs.hpp>
#include <triqs/mesh.hpp>
#include <triqs/arrays/block_matrix.hpp>
#include <triqs/operators/many_body_operator.hpp>
#include <triqs/hilbert_space/fundamental_operator_set.hpp>
#include <triqs/operators/util/extractors.hpp>
 
#include <itertools/itertools.hpp>
#include <mpi/mpi.hpp>
 
#include <iostream>
#include <string>
#include <utility>
#include <variant>
 
namespace triqs_ctint {
 
  using namespace std::complex_literals; // Complex Unity 1i
  using namespace triqs;
  using namespace triqs::gfs;
  using namespace nda;
  using namespace triqs::operators;
  using namespace triqs::operators::utils;
  using namespace triqs::hilbert_space;
  using namespace triqs::utility;
  using namespace h5;
 
  using namespace itertools;

  enum class Chan_t { PP, PH, XPH };

  using alpha_t = array<double, 4>;

  using triqs::gfs::gf_struct_t;

#ifdef GTAU_IS_COMPLEX
  using g_tau_t = block_gf<imtime, matrix_valued>;
#else
  using g_tau_t = block_gf<imtime, matrix_real_valued>;
#endif

  using g_tau_cv_t = g_tau_t::const_view_type;

  using g_tau_v_t = g_tau_t::view_type;

  using g_tau_scalar_t = g_tau_t::g_t::scalar_t;

  using g_iw_t    = block_gf<mesh::dlr_imfreq, matrix_valued>;
  using g_dlr_iw_t = g_iw_t; // explicit alias used by M_iw measurement

  using g_dlr_iw_cv_t = g_iw_t::const_view_type;
  using g_dlr_iw_v_t  = g_iw_t::view_type;

  using g_reg_iw_t    = block_gf<imfreq, matrix_valued>;
  using g_reg_iw_cv_t = g_reg_iw_t::const_view_type;

#if defined GTAU_IS_COMPLEX || defined INTERACTION_IS_COMPLEX
  using M_tau_target_t = matrix_valued;
#else
  using M_tau_target_t = matrix_real_valued;
#endif
  using M_tau_scalar_t = M_tau_target_t::scalar_t;

#ifdef INTERACTION_IS_COMPLEX
  using U_scalar_t = dcomplex;
#else
  using U_scalar_t = double;
#endif

  using mc_weight_t = decltype(U_scalar_t{} * g_tau_scalar_t{});

  using block_matrix_t = std::vector<matrix<M_tau_scalar_t>>;

  using block_matrix_v_t = std::vector<matrix_view<M_tau_scalar_t>>;

  using chi2_iw_t = block2_gf<imfreq, tensor_valued<4>>;

  using chi2_tau_t = block2_gf<imtime, tensor_valued<4>>;

  using chi2_tau_v_t = chi2_tau_t::view_type;

  using chi2_tau_cv_t = chi2_tau_t::const_view_type;

  using chi3_iw_t = block2_gf<prod<imfreq, imfreq>, tensor_valued<4>>;

  using chi3_iw_v_t = chi3_iw_t::view_type;

  using chi3_iw_cv_t = chi3_iw_t::const_view_type;

  using chi3_tau_t = block2_gf<prod<imtime, imtime>, tensor_valued<4>>;

  using chi3_tau_v_t = chi3_tau_t::view_type;

  using chi3_tau_cv_t = chi3_tau_t::const_view_type;

  using chi4_iw_t = block2_gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4>>;

  using chi4_tau_t = block2_gf<prod<imtime, imtime, imtime>, tensor_valued<4>>;
 
  // Declare some placeholders for the rest of the code. Use anonymous namespace for proper linkage
  // in this code, all variables with trailing _ are placeholders by convention.
  namespace {
    using nda::clef::placeholder;
    const placeholder<0> i_;
    const placeholder<1> j_;
    const placeholder<2> k_;
    const placeholder<3> l_;
    const placeholder<4> iw_;
    const placeholder<5> iw1_;
    const placeholder<6> iw2_;
    const placeholder<7> iw3_;
    const placeholder<8> iw4_;
    const placeholder<9> t_;
    const placeholder<10> t1_;
    const placeholder<11> t2_;
    const placeholder<12> t3_;
    const placeholder<13> bl_;
    const placeholder<14> bl1_;
    const placeholder<15> bl2_;
    const placeholder<16> iW_;
    const placeholder<17> iwp_;
  } // anonymous namespace
 
} // namespace triqs_ctint
 
namespace triqs::gfs {

  //FIXME Implement with is_array trait
  template <typename Value_t, int Rank> double max_norm(array_const_view<Value_t, Rank> const &arr) {
    auto max_itr = std::max_element(arr.begin(), arr.end(), [](auto a, auto b) { return std::abs(a) < std::abs(b); });
    return std::abs(*max_itr);
  }
  //FIXME array_const_view should be constructable from array_view
  template <typename Value_t, int Rank> double max_norm(array_view<Value_t, Rank> const &arr) { return max_norm(make_const_view(arr)); }

  template <typename Gf> std::enable_if_t<is_gf_v<Gf>, double> max_norm(Gf const &G) { return max_norm(G.data()); }
 
  template <typename M, typename Target = tensor_valued<4>> block2_gf<M, Target> make_block2_gf(M const &m, gf_struct_t const &gf_struct) {
 
    std::vector<std::vector<gf<M, Target>>> gf_vecvec;
    std::vector<std::string> block_names;
 
    for (auto const &[bl1, bl1_size] : gf_struct) {
      block_names.push_back(bl1);
      std::vector<gf<M, Target>> gf_vec;
      for (auto const &[bl2, bl2_size] : gf_struct) {
        if constexpr (Target::rank == 4)
          gf_vec.emplace_back(m, make_shape(bl1_size, bl1_size, bl2_size, bl2_size));
        else
          gf_vec.emplace_back(m, make_shape(bl1_size, bl2_size));
      }
      gf_vecvec.emplace_back(std::move(gf_vec));
    }
 
    return make_block2_gf(block_names, block_names, std::move(gf_vecvec));
  }
 
  template <typename M1, typename M2, typename Target = tensor_valued<4>>
  block2_gf<M1, Target> make_block2_gf(M1 const &m, block2_gf_const_view<M2, Target> g_in) {
 
    std::vector<std::vector<gf<M1, Target>>> gf_vecvec;
 
    int n_blocks0 = g_in.block_names()[0].size();
    int n_blocks1 = g_in.block_names()[1].size();
 
    for (int i : range(n_blocks0)) {
      std::vector<gf<M1, Target>> gf_vec;
      for (int j : range(n_blocks1)) { gf_vec.emplace_back(m, g_in(i, j).target_shape()); }
      gf_vecvec.emplace_back(std::move(gf_vec));
    }
 
    return make_block2_gf(g_in.block_names()[0], g_in.block_names()[1], std::move(gf_vecvec));
  }
 
  template <typename Scalar_t> std::vector<matrix<Scalar_t>> make_block_vector(gf_struct_t const &gf_struct) {
 
    std::vector<matrix<Scalar_t>> res;
    for (auto const &[bl, bl_size] : gf_struct) { res.emplace_back(nda::zeros<Scalar_t>(bl_size, bl_size)); }
    return res;
  }
} // namespace triqs::gfs
 
namespace triqs::operators {
  inline bool is_densdens_interact(monomial_t m) { return m.size() == 4 and m[0].indices == m[3].indices and m[1].indices == m[2].indices; }

  inline std::pair<int, int> get_int_indices(canonical_ops_t const &op, hilbert_space::gf_struct_t const &gf_struct) {
 
    // The Fundamental operator-set allows for easy check of index validity
    hilbert_space::fundamental_operator_set fs(gf_struct);
    if (!fs.has_indices(op.indices)) TRIQS_RUNTIME_ERROR << " Index of c/c^+ operator not compatible with Green Function structure ";
 
    // Get block-name with apply visitor, lambda(0) is called to determine return type ...
    std::string op_bl_name = visit([](auto idx) { return std::to_string(idx); }, op.indices[0]);
    long nonbl_int_idx     = std::get<long>(op.indices[1]);
 
    // Capture positions in block and nonblock list
    for (auto [bl_int_idx, bl] : itertools::enumerate(gf_struct)) {
      auto const &[bl_name, bl_size] = bl;
      if (bl_name == op_bl_name and 0 <= nonbl_int_idx and nonbl_int_idx < bl_size) { return std::make_pair(bl_int_idx, nonbl_int_idx); }
    }
    TRIQS_RUNTIME_ERROR << "Error: Failed to retrieve integer indices for operator";
  }
 
  // Function that takes a bosonic operator Op = Sum_i a_i c^+(bi, ui) c(bi, vi)
  // and returns a vector<tuple> with v[i] = (b_i, a_i, (c^+(bi, ui), c(bi, vi)))
  inline auto get_terms(many_body_operator const &A, hilbert_space::gf_struct_t const &gf_struct) {
    std::vector<std::tuple<std::complex<double>, std::pair<int, int>, std::pair<int, int>>> terms;
    for (auto const &term : A) {
      auto const &m = term.monomial;
      if (m.size() != 2 or !m[0].dagger or m[1].dagger)
        TRIQS_RUNTIME_ERROR << " Monomial in bosonic operator of chiAB measurement not of the proper form c^+ c \n";
      auto [bl1, i] = get_int_indices(m[0], gf_struct);
      auto [bl2, j] = get_int_indices(m[1], gf_struct);
      auto bl_pair  = std::make_pair(bl1, bl2);
      auto idx_pair = std::make_pair(i, j);
      terms.emplace_back(term.coef, bl_pair, idx_pair);
    }
    return terms;
  }
 
} // namespace triqs::operators
 
// Useful macros
 
#define STR(x) #x
#define STRINGIZE(x) STR(x)
 
#ifdef DEBUG_CTINT
#define TRIQS_EXCEPTION_SHOW_CPP_TRACE
#define __FILENAME__ (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
#define DEBUG(X) std::cerr << AS_STRING(X) << " = " << X << "      at " << __FILENAME__ << ':' << __LINE__ << std::endl
#define BREAK(X)                                                                                                                                     \
  std::cerr << X << " ... " << std::endl;                                                                                                            \
  getchar()
#else
#define DEBUG(X)
#define BREAK(X)
#endif