issue about data vector computation in the Fourier space

[junzhou502]

below is part of the function compute_X_N_masked() in generic_interface.hpp. it uses functions like C_ss_tomo_limber to compute the datavector in fourier space, however it seems faster and more accurate to use C_ss_tomo_limber_noninterp, because the first one will initiate Cells over 300 different ells and obtain angular power spectrum by interpolation, while the second one compute the Cell directly for about 20 ells user requested.

template <int N, int M> 
void compute_X_N_masked(arma::Col<double>& dv, const int start)
{
  static constexpr std::string_view errbegins = "Begins Execution"sv;
  static constexpr std::string_view errends = "Ends Execution"sv;
  static constexpr std::string_view fname = "compute_X_N_masked"sv;
  using spdlog::debug; using spdlog::info;
  debug("{}<{},{}>: {}", fname, N, M, errbegins);
  static_assert(0 == N || 1 == N, "N must be 0 (real) or 1 (fourier)");
  IP& survey = IP::get_instance();
  arma::Col<int>::fixed<2> Nlen = {Ntable.Ntheta, like.Ncl};

  if constexpr (0 == M) {
    if (1 == like.shear_shear) {    
      for (int nz=0; nz<tomo.shear_Npowerspectra; nz++) {
        const int z1 = Z1(nz);
        const int z2 = Z2(nz);
        for (int i=0; i<Nlen[N]; i++) {
          int index = start + Nlen[N]*nz + i;
          if constexpr (N == 0) {
            if (survey.get_mask(index)) {
              dv(index) = xi_pm_tomo(1, i, z1, z2, 1);
            }  
            index += Nlen[N]*tomo.shear_Npowerspectra;
            if (survey.get_mask(index)) {
              dv(index) = xi_pm_tomo(-1, i, z1, z2, 1);
            } 
          }
          else {
            if (survey.get_mask(index) && (like.ell[i]<like.lmax_shear)) {
!              dv(index) = C_ss_tomo_limber(like.ell[i], z1, z2, 1);
            }
          }
        }
      }
      add_calib_and_set_mask_X_N<N,M>(dv, start);
    }
  }
  else if constexpr (1 == M) {
    if (1 == like.shear_pos) {
      for (int nz=0; nz<tomo.ggl_Npowerspectra; nz++) {
        const int zl = ZL(nz);
        const int zs = ZS(nz);
        for (int i=0; i<Nlen[N]; i++) {
          const int index = start + Nlen[N]*nz + i;
          if (survey.get_mask(index)) {
            if constexpr (0 == N)
              dv(index) = w_gammat_tomo(i,zl,zs,like.adopt_limber_gs);
            else
!              dv(index) = C_gs_tomo_limber(like.ell[i], zl, zs);
          }
        }
      }
      add_calib_and_set_mask_X_N<N,M>(dv, start);
    }
  }
  else if constexpr (2 == M) {
    if (1 == like.pos_pos) {
      for (int nz=0; nz<tomo.clustering_Npowerspectra; nz++) {
        for (int i=0; i<Nlen[N]; i++) {
          const int index = start + Nlen[N]*nz + i;
          if (survey.get_mask(index)) {
            if constexpr (N == 0) {  
              dv(index) = w_gg_tomo(i, nz, nz, like.adopt_limber_gg);
            }
            else {
!              dv(index) = C_gg_tomo_limber(like.ell[i], nz, nz);
            }
          }
        }
      }
      add_calib_and_set_mask_X_N<N,M>(dv, start);
    }
  }

[vivianmiranda]

Yes absolutely - we need to change the fourier to use the non-interpolator version :)