Python interface for 3PCF

Makefile

@@ -4,7 +4,7 @@ CC = gcc
 CFLAGS = -O3 -Wall -MMD
 CXX = g++
 CXXFLAGS	= -O3 -Wall -MMD -std=c++11 -ffast-math $(shell pkg-config --cflags gsl)
-CXXFLAGS	+= -DOPENMP -DLEGENDRE_MIX -DJACKKNIFE -DPRINTPERCENTS
+CXXFLAGS	+= -DOPENMP -DTHREE_PCF -DPRINTPERCENTS
 #-DOPENMP  # use this to run multi-threaded with OPENMP
 #-DPERIODIC # use this to enable periodic behavior
 #-DLEGENDRE # use this to compute 2PCF covariances in Legendre bins (original mode, corresponding to direct accumulation into multipoles from pair counts)

RascalC/Makefile

@@ -34,15 +34,16 @@ endif
 
 LD = $(CXX)
 
-# The code below compiles all the valid variants for the 2PCF covariance
+# The code below compiles all the valid variants for the 2PCF covariance, and one for 3PCF
 
-BASE_VARIANTS	= s_mu s_mu_jackknife legendre_accumulated legendre_projected legendre_projected_jackknife
+BASE_VARIANTS	= s_mu s_mu_jackknife legendre_accumulated legendre_projected legendre_projected_jackknife 3pcf_legendre_accumulated
 
 DEFINES_FOR_s_mu=
 DEFINES_FOR_s_mu_jackknife=-DJACKKNIFE
 DEFINES_FOR_legendre_accumulated=-DLEGENDRE
 DEFINES_FOR_legendre_projected=-DLEGENDRE_MIX
 DEFINES_FOR_legendre_projected_jackknife=-DLEGENDRE_MIX -DJACKKNIFE
+DEFINES_FOR_3pcf_legendre_accumulated=-DTHREE_PCF
 
 define add_periodic_variant
 DEFINES_FOR_$(1)_periodic	= $$(DEFINES_FOR_$(1)) -DPERIODIC
@@ -90,6 +91,37 @@ endef
 
 $(foreach variant,$(VARIANTS),$(eval $(call make_targets,$(variant))))
 
+# end of variants for the covariance
+
+# The code below compiles all the valid variants for the triple_counts
+
+VARIANTS2 = s_mu s_mu_verbose s_mu_periodic s_mu_periodic_verbose
+
+DEFINES2 = -DTRIPLE # triple_counts needs this to have mumin=-1 and not 0 to be consistent with the 3PCF covariance
+
+VAR_SRC2_BASE	= triple_counts
+VAR_SRC2	= ../$(VAR_SRC2_BASE)/$(VAR_SRC2_BASE).cpp
+VAR_OBJ2_BASE	= $(OBJDIR)/$(VAR_SRC2_BASE)
+
+EXEC2_BASE	= $(EXECDIR)/triple
+
+define make_targets2
+OBJ2_$(1)	= $$(VAR_OBJ2_BASE).$(1).o
+ALL_OBJS	+= $$(OBJ2_$(1))
+EXEC2_$(1)	= $$(EXEC2_BASE).$(1)
+ALL_EXECS	+= $$(EXEC2_$(1))
+
+$$(OBJ2_$(1)): $$(VAR_SRC2) | $$(OBJDIR)
+	$$(CXX) $$(CXXFLAGS) $$(DEFINES2) $$(DEFINES_FOR_$(1)) -c $$(VAR_SRC2) -o $$@
+
+$$(EXEC2_$(1)): $$(OBJ2_$(1)) $$(COMMON_OBJS) | $$(EXECDIR)
+	$$(LD) $$(OBJ2_$(1)) $$(COMMON_OBJS) $$(LFLAGS) -o $$@
+endef
+
+$(foreach variant,$(VARIANTS2),$(eval $(call make_targets2,$(variant))))
+
+# end of variants for the triple_counts
+
 .PHONY: all clean
 
 all: $(ALL_EXECS)

RascalC/__init__.py

@@ -1,5 +1,7 @@
-__version__ = "3.2"
+__version__ = "3.3.alpha1"
 
 from .interface import run_cov
 from .post_process import post_process_auto
-from .raw_covariance_matrices import cat_raw_covariance_matrices
+from .raw_covariance_matrices import cat_raw_covariance_matrices
+
+from .interface_3pcf import run_cov_3pcf

RascalC/convergence_check_extra.py

@@ -83,4 +83,4 @@ def convergence_check_extra_file(rascalc_results_filename: str, n_samples: int |
     Further optionally, use ``print_function`` to report the results.
     """
     with np.load(rascalc_results_filename) as f:
-        return convergence_check_extra(f, n_samples, print_function)
+        return convergence_check_extra(f, n_samples, print_function)

RascalC/cov_utils_3pcf.py

@@ -0,0 +1,44 @@
+"""
+This module contains convenience functions that
+
+- read RascalC results (checking eigenvalues of the inversion bias matrix. If they are much smaller than 1, it is safe to simply invert the covariance matrix. Otherwise, a correction factor is necessary.
+- convert (reorder) 3PCF covariance matrices (as is usually needed for Legendre moments)
+- and/or export (save) the full 3PCF covariance matrix to a text file.
+"""
+
+import numpy as np
+import numpy.typing as npt
+from typing import Callable
+from .get_shot_noise_rescaling import get_shot_noise_rescaling # for convenience
+from .cov_utils import get_cov_header, load_cov
+
+
+def convert_cov_3pcf_legendre(cov: npt.NDArray[np.float64], max_l: int, exclude_odd_l: bool = False, apply_scaling: bool = True) -> npt.NDArray[np.float64]:
+    """
+    Change the bin ordering of the 3PCF covariance matrix in Legendre mode for a single tracer.
+    The bin ordering in ``RascalC`` ``.npy`` files is by radial bin pairs (top-level) and then by multipoles.
+    The resulting bin ordering (for text files) is by multipoles (top-level) and then by radial bin pairs.
+    By default, the ell-dependent scaling factor is applied to convert to ENCORE basis functions from RascalC Legendre polynomial basis, this can be disabled with ``apply_scaling=False``.
+    """
+    n_l = max_l // (1 + exclude_odd_l) + 1
+    n_bins = len(cov)
+    if n_bins % n_l != 0: raise ValueError("Number of bins in the Legendre covariance must be divisible by the number of multipoles (even and odd)")
+    n_r_bin_pairs = n_bins // n_l
+    cov = cov.reshape(n_r_bin_pairs, n_l, n_r_bin_pairs, n_l) # convert to 4D from 2D with [r_bin_pair, l] ordering for both rows and columns
+    if apply_scaling: # prepare to apply the ell-dependent factor between RascalC and ENCORE basis functions
+        ells = np.arange(0, max_l + 1, 1 + exclude_odd_l)
+        ell_factor = ((-1)**ells * np.sqrt(2 * ells + 1) / (4 * np.pi)) # the ell-dependent factor between the ENCORE 3-point basis functions and Legendre polynomials given by Equation (17) in https://arxiv.org/pdf/2105.08722. RascalC uses Legendre polynomial basis according to Equation (4.11) in https://arxiv.org/pdf/1910.04764
+        cov /= ell_factor[:, None, None] * ell_factor[None, None, :] # apply the factor, NumPy broadcasting matches trailing dimensions. need to check if it is not multiplication; there might also be a factor of 2 or something similar
+    cov = cov.transpose(1, 0, 3, 2) # change ordering to [l, r_pair] for both rows and columns
+    cov = cov.reshape(n_bins, n_bins) # convert back from 4D to 2D
+    return cov
+
+
+def load_cov_3pcf_legendre(rascalc_results_file: str, max_l: int, exclude_odd_l: bool = False, apply_scaling: bool = True, print_function: Callable[[str], None] = print) -> npt.NDArray[np.float64]:
+    "Load the theoretical covariance matrix from RascalC results file, change the bin ordering and apply scaling as in :func:`convert_cov_3pcf_legendre`; intended for 3PCF Legendre single-tracer mode."
+    return convert_cov_3pcf_legendre(load_cov(rascalc_results_file, print_function), max_l, exclude_odd_l=exclude_odd_l, apply_scaling=apply_scaling)
+
+
+def export_cov_3pcf_legendre(rascalc_results_file: str, max_l: int, output_cov_file: str, exclude_odd_l: bool = False, apply_scaling: bool = True, print_function: Callable[[str], None] = print) -> None:
+    "Export the theoretical covariance matrix from RascalC results file to a text file with conversion appropriate for 3PCF single-tracer Legendre mode (see :func:`convert_cov_3pcf_legendre`)."
+    np.savetxt(output_cov_file, load_cov_3pcf_legendre(rascalc_results_file, max_l, exclude_odd_l=exclude_odd_l, apply_scaling=apply_scaling, print_function = print_function), header = get_cov_header(rascalc_results_file))