diff --git a/Makefile b/Makefile index ab6b1944..2566a62c 100644 --- a/Makefile +++ b/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) diff --git a/RascalC/Makefile b/RascalC/Makefile index bcc54ec2..2afdde0c 100644 --- a/RascalC/Makefile +++ b/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) diff --git a/RascalC/__init__.py b/RascalC/__init__.py index 93654904..07ab1226 100644 --- a/RascalC/__init__.py +++ b/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 \ No newline at end of file +from .raw_covariance_matrices import cat_raw_covariance_matrices + +from .interface_3pcf import run_cov_3pcf \ No newline at end of file diff --git a/RascalC/convergence_check_extra.py b/RascalC/convergence_check_extra.py index 90d4dad6..896ec110 100644 --- a/RascalC/convergence_check_extra.py +++ b/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) \ No newline at end of file + return convergence_check_extra(f, n_samples, print_function) diff --git a/RascalC/correction_function_3pcf.py b/RascalC/correction_function_3pcf.py index 57a79126..9aa6c731 100644 --- a/RascalC/correction_function_3pcf.py +++ b/RascalC/correction_function_3pcf.py @@ -1,110 +1,229 @@ -### Function to fit a model (piecewise polynomial) to the 3PCF survey correction function, defined as the ratio between idealistic and true RRR pair counts for a single survey. - -## NB: Input RRR counts should be normalized by the cube of the sum of random weights here. ## is this useful/necessary? -## NB: Assume mu is in [-1,1] limit here - - import os import numpy as np -import scipy.spatial as ss +import numpy.typing as npt from typing import Callable +from scipy.special import legendre +from .correction_function import compute_V_n_w_bar, load_randoms + + +def compute_inv_phi_periodic_3pcf(n: int, n_multipoles: int) -> npt.NDArray[np.float64]: + "Compute the inverse 3PCF survey correction function multipoles for the periodic box geometry." + ## Output periodic survey correction function + phi_inv_mult = np.zeros([n, n, n_multipoles]) + + ## Set to correct periodic survey values + phi_inv_mult[:, :, 0] = 1 + + return phi_inv_mult + + +def check_triple_counts_positive(leg_triple: npt.NDArray[np.float64], lenient_samebins: bool = False, print_function: Callable[[str], None] = print) -> None: + "Check for negative counts, which should be problematic" + n_chi = 2001 + triple_counts = np.zeros(list(leg_triple.shape[:-1]) + [n_chi]) + chi_values = np.linspace(-1, 1, n_chi) + for ell in range(leg_triple.shape[-1]): + triple_counts += leg_triple[:, :, ell][:, :, None] * legendre(ell)(chi_values)[None, None, :] + problem_indices = np.argwhere(triple_counts <= 0) + if len(problem_indices) > 0: + rbins, chi_counts = np.unique(problem_indices[:, :2], axis=0, return_counts=True) + for ((rbin1, rbin2), chi_count) in zip(rbins, chi_counts): + if rbin1 > rbin2: continue # this case can be skipped by symmetry + title = "WARNING" + if lenient_samebins and rbin1 == rbin2: title = "INFO" # the problem for same-bin pairs is less critical (and seems more likely), for different-bin pairs it is more critical + print_function(f"{title}: counts are not positive for radial bin pair {rbin1}, {rbin2} for {chi_count} chi values of {n_chi} checked") + + +def check_inv_phi_values(phi_inv_mult: npt.NDArray[np.float64], print_function: Callable[[str], None] = print) -> None: + """ + Check that the mean of the monopole is neither too small nor too large. + Its values should be close to 1, which is the value for the ideal survey. But it is not clear how far they can deviate for a realistic survey with complicated geometry, varying density and weights. Moreover, the values also depend on the volume estimate, which may be off by a factor of a few. + So, this check is given a wide margin, but it may still be useful to catch some obvious problems with the RRR counts normalization. + """ + print_function(f"INFO: mean±std of the monopole of the inverse survey correction function is {np.mean(phi_inv_mult[:, :, 0])}±{np.std(phi_inv_mult[:, :, 0], ddof=1)}, expected to be not very far from 1. NB: this diagnostic depends on an empirical volume estimate with ConvexHull, which may be off (overestimated) by a factor of a few for realistic survey geometries.") + if np.mean(phi_inv_mult[:, :, 0]) < 1e-3: + raise ValueError("Survey correction function seems too small - are the RRR counts normalized correctly?") + if np.mean(phi_inv_mult[:, :, 0]) > 1e3: + raise ValueError("Survey correction function seems too large - are the RRR counts normalized correctly?") + + +def compute_inv_phi_aperiodic_3pcf(n: int, m: int, n_multipoles: int, r_bins: npt.NDArray[np.float64], RRRbin_over_3Vn3v3: npt.NDArray[np.float64], print_function: Callable[[str], None] = print) -> npt.NDArray[np.float64]: + "Compute the inverse 3PCF survey correction function multipoles for the realistic survey geometry." + + chi_all = np.linspace(-1,1,m+1) + chi_cen = 0.5*(chi_all[1:]+chi_all[:-1]) + + ## reshape RRR counts and add symmetries + # 3Vn3v3 is just part of the normalization for the (inverse) survey correction function, see Equation 4.10 of https://arxiv.org/pdf/1910.04764 + RRRbin_over_3Vn3v3 = RRRbin_over_3Vn3v3.reshape(n, n, m) + RRRbin_over_3Vn3v3 = (RRRbin_over_3Vn3v3 + RRRbin_over_3Vn3v3.swapaxes(0, 1)) / 2 # triple_counts code accumulates each triple to the three possible pairs of radial bins, but only to one of the two possible orderings of the pair, with twice the weight. so this symmetrization should give the counts where each triple contributes to the 6 bin triplets it can according to the RascalC convention (see Section 4.1 of https://arxiv.org/pdf/1910.04764, Equations 4.2-4.4). in the end, this is exactly symmetric under the interchange of the two radial bins, as it should be but wasn't before this step + + ## Now construct Legendre moments + RRRleg_over_3Vn3v3 = np.zeros([n, n, n_multipoles]) + for ell in range(n_multipoles): + # we've absorbed a factor of delta_chi into RRRbin... here, because RRRbin is RRR^{ab}_c ≈ RRR^{ab}(chi_central) delta_chi in terms of Section 4.2 of https://arxiv.org/pdf/1910.04764 + # we want the Legendre polynomial decomposition of the inverse survey correction function 1/Phi(r_a, r_b, chi), which is proportional to RRR^ab(chi) according to Equation 4.10 of https://arxiv.org/pdf/1910.04764 + # so what we want is proportinal to RRR^{ab}_ell = (2 ell + 1)/2 int_{-1}^{1} L_ell(chi) RRR^{ab}(chi) dchi ≈ (2 ell + 1)/2 sum_c L_ell(chi_central_c) RRR^{ab}(chi_central_c) delta_chi ≈ (2 ell + 1)/2 sum_c L_ell(chi_central_c) RRR^{ab}_c + # absorption of delta_chi makes sense, but the lack of division by 2 is not completely understood (whereas it is empirically validated on the ENCORE side, see below) + RRRleg_over_3Vn3v3[:, :, ell] += (2.*ell+1.) * np.sum(legendre(ell)(chi_cen)[None, None, :] * RRRbin_over_3Vn3v3, axis=-1) + + # as a precaution, check for negative counts + check_triple_counts_positive(RRRleg_over_3Vn3v3, n_multipoles, print_function=print_function) + + vol_r = 4 * np.pi / 3 * (r_bins[:, 1] **3 - r_bins[:, 0] ** 3) + + ## Construct multipoles of inverse Phi + phi_inv_mult = RRRleg_over_3Vn3v3 / (vol_r[:, None, None] * vol_r[None, :, None]) # here we restore the bin volume factors from the numerator in Equation 4.10 of https://arxiv.org/pdf/1910.04764. so we finally get to the Legendre polynomial decomposition coefficients of 1/Phi(r_a, r_b, chi) + + # Check that values seem reasonable + check_inv_phi_values(phi_inv_mult, print_function=print_function) + + return phi_inv_mult + +def compute_3pcf_correction_function(randoms_pos: npt.NDArray[np.float64], randoms_weights: npt.NDArray[np.float64], binfile: str, outdir: str, periodic: bool, RRR_file: str | None = None, print_function: Callable[[str], None] = print) -> str: + """ + Function to compute the multipole decomposition of the 3PCF inverse survey correction function. + The 3PCF survey correction function is defined as the ratio between idealistic and true RRR pair counts for a single survey. + + NB: Input RRR counts should be normalized by the cube of the sum of random weights here. + NB: Assume chi is in [-1,1] limit here - it is the cosine of the angle in the triangle, and its sign is important (unlike for mu in the 2PCF context) + """ + + n_multipoles = 7 # matches the value hard-coded in the C++ code -def compute_3pcf_correction_function(random_filename: str, binfile: str, outdir: str, periodic: bool, RRR_file: str | None = None, print_function: Callable[[str], None] = print) -> None: if periodic: - print_function("\nAssuming periodic boundary conditions - so Phi(r,mu) = 1 everywhere") + print_function("Assuming periodic boundary conditions - so Phi(r_a, r_b, chi) = 1 everywhere") elif RRR_file is None: raise TypeError("RRR counts file is required if aperiodic") - ## Load galaxies - print_function("\nLoading randoms") - randoms_data = np.loadtxt(random_filename) - random_pos = randoms_data[:, :3] - random_w = randoms_data[:, 3] - - N_gal = len(randoms_data) - w_bar = np.mean(random_w) - ## Find survey volume via ConvexHull in Scipy - hull = ss.ConvexHull(random_pos) - print_function('\nSurvey volume is approximately: %.2f (Gpc/h)^3'%(hull.volume/1e9)) - V=hull.volume # in (Mpc/h)^3 - - ## Galaxy number density - n_bar = N_gal/V + V, n_bar, w_bar = compute_V_n_w_bar(randoms_pos, randoms_weights) # Load in binning files r_bins = np.loadtxt(binfile) n=len(r_bins) ## Define normalization constant - norm = 6. * V * n_bar**3 * w_bar**3 # I don't think there is an exactly right answer once number density or weights vary across the survey + norm = 3. * V * n_bar**3 * w_bar**3 # this nearly corresponds to the definitions in the paper (the numerator in Equation 4.10 of https://arxiv.org/pdf/1910.04764, barring the bin volumes). It doesn't match exactly, but it is easier to compute in practice, and I don't think there is a completely right answer once number density or weights vary across the survey + # this is also equal to 3 * np.sum(randoms_weights)**3 / V**2; this form is more useful for thinking about scaling with random weight normalization or the volume estimate + # the value of the norm is actually important in the periodic case - it scales the final output + # in the aperiodic case, the final output does not depend on this norm value (see below), but it is used for a sanity check on the values (albeit with a wide margin) print_function("Normalizing output survey correction by %.2e"%norm) if periodic: - - ## Output periodic survey correction function - phi_inv_mult = np.zeros([n,n,7]); - - ## Set to correct periodic survey values - phi_inv_mult[:,:,0]=1. + phi_inv_mult = compute_inv_phi_periodic_3pcf(n, n_multipoles) else: - from scipy.special import legendre - - ## Load triple counts and renormalize - tmp_triple_counts = np.loadtxt(RRR_file)*np.sum(random_w)**3 - + triple_counts = np.loadtxt(RRR_file)*np.sum(randoms_weights)**3 + # Compute number of angular bins in data-set - m = (len(tmp_triple_counts)//n)//n - if len(tmp_triple_counts) % m != 0: raise ValueError("Incorrect RRR format") - + m = (len(triple_counts)//n)//n + if len(triple_counts) % m != 0: raise ValueError("Incorrect RRR format") - mu_all = np.linspace(-1,1,m+1) - mu_cen = 0.5*(mu_all[1:]+mu_all[:-1]) - - RRR_true = np.zeros([n,n,m]) - - ## load in RRR counts (and add symmetries) - for i in range(len(tmp_triple_counts)): - RRR_true[(i//m)//n,(i//m)%n,i%m] += tmp_triple_counts[i]*0.5 - RRR_true[(i//m)%n,(i//m)//n,i%m] += tmp_triple_counts[i]*0.5 - - ## Now construct Legendre moments - leg_triple = np.zeros([n,n,7]) - for a in range(n): - for b in range(n): - for ell in range(7): - # (NB: we've absorbed a factor of delta_mu into RRR_true here) - leg_triple[a,b,ell]+=np.sum(legendre(ell)(mu_cen)*RRR_true[a,b,:])*(2.*ell+1.) - - vol_r = lambda b: 4.*np.pi/3.*(r_bins[b,1]**3.-r_bins[b,0]**3.) - - ## Construct inverse multipoles of Phi - phi_inv_mult = np.zeros([n,n,7]) - for b1 in range(n): - for b2 in range(n): - phi_inv_mult[b1,b2,:] = leg_triple[b1,b2,:] / (.5 * norm * vol_r(b1) * vol_r(b2)) - - ## Check all seems reasonable - if np.mean(phi_inv_mult[:,:,0])<1e-3: - print_function(phi_inv_mult[:,:,0]) - raise ValueError("Survey correction function seems too small - are the RRR counts normalized correctly?") - if np.mean(phi_inv_mult[:,:,0])>1e3: - raise ValueError("Survey correction function seems too large - are the RRR counts normalized correctly?") + phi_inv_mult = compute_inv_phi_aperiodic_3pcf(n, m, n_multipoles, r_bins, triple_counts / norm, print_function=print_function) if periodic: outfile = os.path.join(outdir, 'BinCorrectionFactor3PCF_n%d_periodic.txt'%(n)) else: outfile = os.path.join(outdir, 'BinCorrectionFactor3PCF_n%d_m%d.txt'%(n,m)) + + np.savetxt(outfile, phi_inv_mult.reshape(n*n, n_multipoles) * norm * 2) + # saved object is not the survey correction function (or its inverse) as defined in the paper (Equation 4.10 of https://arxiv.org/pdf/1910.04764). the reason is probably that we don't want to deal with the survey volume and weights in the C++ code. rather, the saved object is the Legendre polynomial decomposition of 6 * V * n_bar^3 * w_bar^3 / Phi(r1, r2, mu), the inverse of which appears in the covariance estimators (Equation 5.13 of https://arxiv.org/pdf/1910.04764 combined with Equation 4.16 defining/explaining the K factors) as the correction_factor in modules/integrals_3pcf.h. not sure why 2 from the K factor couldn't be absorbed to cancel the 2 in norm * 2 (it is applied in the integrals in the C++ code instead, see modules/integrals_3pcf.h), but probably easier to leave it be + # the bin volume and (2*p+1) factors are re-applied inside the C++ code at the end of the computation with the normalize() function in modules/integrals_3pcf.h + print_function("Saved correction factors to %s\n"%outfile) + + return outfile + + +def compute_3pcf_correction_function_from_files(random_filename: str, binfile: str, outdir: str, periodic: bool, RRR_file: str | None = None, print_function: Callable[[str], None] = print) -> str: + print_function("Loading randoms") + return compute_3pcf_correction_function(*load_randoms(random_filename), binfile, outdir, periodic, RRR_file, print_function = print_function) + + +def compute_3pcf_correction_function_from_encore(randoms_pos: npt.NDArray[np.float64], randoms_weights: npt.NDArray[np.float64], binfile: str, outdir: str, triple_counts: npt.NDArray[np.float64], print_function: Callable[[str], None] = print) -> str: + """ + Function to compute the multipole decomposition of the 3PCF inverse survey correction function from ENCORE triple counts. + The 3PCF survey correction function is defined as the ratio between idealistic and true RRR pair counts for a single survey. + + NB: Input RRR counts are not normalized here, and are already in multipole format. + Caveat: ENCORE only computes the triple counts for pairs of different radial bins, whereas RascalC expects them for all pairs of radial bins, crucially including the pairs of identical bins. Here we try to fill the missing data for those identical-bin pairs using the neighboring bin pairs. These should only affect the covariance rows and columns corresponding to the identical-bin pairs, which should be removed from the covariance for use with ENCORE 3PCF measurements in the end. So, those missing bin pairs should not matter in the end, but it is nicer not to have complete nonsence in the intermediate products. + """ + + n_multipoles = 7 # matches the value hard-coded in the C++ code + + ells = np.arange(len(triple_counts)) # rows in triple_counts correspond to multipoles, columns - to radial bin pairs; the first column might just contain these ells + + if np.array_equal(triple_counts[:, 0], ells): triple_counts = triple_counts[:, 1:] # remove the first column if it is all ells + + # Load in binning files + r_bins = np.loadtxt(binfile) + n = len(r_bins) + + if triple_counts.shape[1] != n*(n-1)//2: raise ValueError("The shape of RRR_counts is inconsistent with the radial bins provided") + # check this after removing the ells column if present + # the columns correspond to radial bins + + # conversion from ENCORE to RascalC (triple_counts) convention + triple_counts = 6 * triple_counts # the factor of 6 comes from the fact that ENCORE counts each triplet only once, whereas RascalC counts each triplet 6 times (once for each permutation of the three points, see the end of Section 4.2 below Equation 4.16 in https://arxiv.org/pdf/1910.04764). both should encounter every triplet in every relevant permutations on top of this. (don't use *= at first to avoid overriding the variable passed to the function) + triple_counts *= 8 * np.pi**2 # the factor of 8 pi^2 comes from the fact that ENCORE seems to weigh each triplet by its basis functions (given by Equation 17 in https://arxiv.org/pdf/2105.08722). it is not the same as projecting onto the basis function, because the square norm of the basis function, its integral over chi (cosine of the angle in the triangle) from -1 to 1, is not 1 but rather 1/(8 pi^2). so, to get the coefficient in the basis decomposition, we need to multiply by 8 pi^2. (this issue probably only applies to counts from ENCORE. its 3PCF are obtained by division by the RRR monopole, which should cancel the constant normalization factor. so in the other two places in the code, where we convert the 3PCF between conventions, the basis conversion factor below is sufficient) + # convert ell coefficients from ENCORE basis to simple Legendre polynomials used in RascalC (according to Equation 4.11 in https://arxiv.org/pdf/1910.04764) + triple_counts *= ((-1)**ells * np.sqrt(2 * ells + 1) / (4 * np.pi))[:, None] # add the second dimension, corresponding to the radial bins, to avoid indexing errors. should be fine to use *= now because triple_counts is a new local variable + # 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 + triple_counts *= 2 # this factor of 2 is not yet completely understood, but it corresponds to the lack of division of 2 that seems to be needed in the decomposition onto Legendre polynomials in compute_inv_phi_aperiodic_3pcf() above, and Karim validated empirically that the resulting final covariance is consistent with the ENCORE mock sample covariance at <~20% level, so it is evidently not off by a factor of 2 (or rather 4) + + # ensure the number of multipoles in triple_counts is right to avoid indexing errors + if len(triple_counts) < n_multipoles: # this seems more likely + print_function(f"INFO: ENCORE triple counts have {len(triple_counts)} multipoles, fewer than {n_multipoles} used for the survey correction function, extending by zeros") + triple_counts = np.vstack([triple_counts, np.zeros([n_multipoles - len(triple_counts), triple_counts.shape[1]])]) + elif len(triple_counts) > n_multipoles: # this seems less likely + print_function(f"INFO: ENCORE triple counts have {len(triple_counts)} multipoles, more than {n_multipoles} used for the survey correction function, discarding the higher multipoles") + triple_counts = triple_counts[:n_multipoles] + + bin_indices = np.arange(n) + bin_index1 = np.repeat(bin_indices, n-1-bin_indices) + bin_index2 = np.concatenate([bin_indices[i+1:] for i in range(n)]) + # bin_index1 and bin_index2 cover all the bin pairs under the condition bin_index1 < bin_index2, the order follows the ENCORE format + # they could be read from first two non-comment rows of the ENCORE file, but this seems unnecessary + + leg_triple = np.zeros([n, n, n_multipoles]) + leg_triple[bin_index1, bin_index2] = triple_counts.T # fill above the diagonal; transposition puts radial bin pair index first and multipole index last + leg_triple[bin_index2, bin_index1] = triple_counts.T # fill below the diagonal symmetrically + # leave zeros at the diagonal for now + + vol_r = 4 * np.pi / 3 * (r_bins[:, 1] ** 3 - r_bins[:, 0] ** 3) # volume of radial/separation bins as 1D array + + V, n_bar, w_bar = compute_V_n_w_bar(randoms_pos, randoms_weights) + + ## Define normalization constant + norm = 3. * V * n_bar**3 * w_bar**3 # this nearly corresponds to the definitions in the paper (the numerator in Equation 4.10 of https://arxiv.org/pdf/1910.04764, barring the bin volumes). It doesn't match exactly, but it is easier to compute in practice, and I don't think there is a completely right answer once number density or weights vary across the survey + # this is also equal to 3 * np.sum(randoms_weights)**3 / V**2; this form is more useful for thinking about scaling with random weight normalization or the volume estimate + # the value of the norm is actually important in the periodic case - it scales the final output + # in the aperiodic case, the final output does not depend on this norm value (see below), but it is used for a sanity check on the values (albeit with a wide margin) + + ## Construct multipoles of inverse Phi + phi_inv_mult = leg_triple / (norm * vol_r[:, None, None] * vol_r[None, :, None]) # here we also add the bin volume factors from the numerator in Equation 4.10 of https://arxiv.org/pdf/1910.04764 + + # fill the middle diagonal elements, which have been zeros + # seems better to do in phi_inv_mult because its values change less + bin_indices_middle = bin_indices[1:-1] + phi_inv_mult[bin_indices_middle, bin_indices_middle] = (phi_inv_mult[bin_indices_middle+1, bin_indices_middle] + phi_inv_mult[bin_indices_middle-1, bin_indices_middle]) / 2 # average the neighboring elements along the column. the neighboring elements along the row are the same due to symmetry + + # fill the edge/corner diagonal elements, which are more tricky and have been zeros + phi_inv_mult[0, 0] = (2 * (2 * phi_inv_mult[1, 0] - phi_inv_mult[2, 0]) + (2 * phi_inv_mult[1, 1] - phi_inv_mult[2, 2])) / 3 + phi_inv_mult[-1, -1] = (2 * (2 * phi_inv_mult[-2, -1] - phi_inv_mult[-3, -1]) + (2 * phi_inv_mult[-2, -2] - phi_inv_mult[-3, -3])) / 3 + + # check for negative triple counts (or rather the correction function as it is passed to RascalC), which should be problematic + check_triple_counts_positive(phi_inv_mult * norm, lenient_samebins=True, print_function=print_function) + + # Check that the values seem reasonable + check_inv_phi_values(phi_inv_mult, print_function=print_function) - with open(outfile,"w+") as out: - for b1 in range(n): - for b2 in range(n): - for ell in range(7): - out.write("%.8e"%(phi_inv_mult[b1,b2,ell]*norm)) - if ell<6: - out.write("\t") - if ell==6: - out.write("\n") - print_function("\nSaved (normalized) output to %s\n"%outfile) \ No newline at end of file + outfile = os.path.join(outdir, 'BinCorrectionFactor3PCF_n%d.txt' % n) + + np.savetxt(outfile, phi_inv_mult.reshape(n*n, n_multipoles) * norm * 2) + # saved object is not the survey correction function (or its inverse) as defined in the paper (Equation 4.10 of https://arxiv.org/pdf/1910.04764). the reason is probably that we don't want to deal with the survey volume and weights in the C++ code. rather, the saved object is the Legendre polynomial decomposition of 6 * V * n_bar^3 * w_bar^3 / Phi(r1, r2, mu), the inverse of which appears in the covariance estimators (Equation 5.13 of https://arxiv.org/pdf/1910.04764 combined with Equation 4.16 defining/explaining the K factors) as the correction_factor in modules/integrals_3pcf.h. not sure why 2 from the K factor couldn't be absorbed to cancel the 2 in norm * 2 (it is applied in the integrals in the C++ code instead, see modules/integrals_3pcf.h), but probably easier to leave it be + # the bin volume and (2*p+1) factors are re-applied inside the C++ code at the end of the computation with the normalize() function in modules/integrals_3pcf.h + print_function("Saved correction factors to %s\n"%outfile) + + return outfile \ No newline at end of file diff --git a/RascalC/cov_utils_3pcf.py b/RascalC/cov_utils_3pcf.py new file mode 100644 index 00000000..38230f63 --- /dev/null +++ b/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)) \ No newline at end of file diff --git a/RascalC/interface_3pcf.py b/RascalC/interface_3pcf.py new file mode 100644 index 00000000..b72461a9 --- /dev/null +++ b/RascalC/interface_3pcf.py @@ -0,0 +1,477 @@ +""" +Python wrapper of ``RascalC`` for 3-point correlation function covariance, being interfaced with the ``ENCORE`` `code for efficient NPCF (including 3PCF) estimation `. + +Please bear with the long description; you can pay less attention to settings labeled optional in the beginning. +""" + +import pycorr +import lsstypes +import numpy as np +import os +from datetime import datetime +import numpy.typing as npt +from typing import Iterable, Literal +from warnings import warn +from .xi.utils import write_xi_file +from .write_binning_file import write_binning_file +from .pycorr_utils.input_xi import get_input_xi_from_pycorr +from .lsstypes_utils.input_xi import get_input_xi_from_lsstypes +from .allcounts_utils import get_s_edges_from_allcounts, get_mu_edges_from_allcounts, get_s_avg_from_allcounts, get_mu_avg_from_allcounts, allcount_switch_function, fix_and_wrap_allcounts +from .correction_function_3pcf import compute_3pcf_correction_function, compute_3pcf_correction_function_from_encore +from .convergence_check_extra import convergence_check_extra +from .utils import rmdir_if_exists_and_empty, suffixes_tracer_all, indices_corr_all, suffixes_corr_all +from .post_process_3pcf import post_process_3pcf + + +def run_cov_3pcf(mode: Literal["legendre_accumulated"], + s_edges: npt.NDArray[np.float64], max_l: int, + nthread: int, N2: int, N3: int, N4: int, N5: int, N6: int, n_loops: int, loops_per_sample: int, + out_dir: str, tmp_dir: str, + randoms_positions1: npt.NDArray[np.float64], randoms_weights1: npt.NDArray[np.float64], + xi_table_11: pycorr.twopoint_estimator.BaseTwoPointEstimator | lsstypes.Count2Correlation | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]] | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]], + no_data_galaxies1: float, effective_no_def: bool = False, + RRR_counts: npt.NDArray[np.float64] | None = None, normalized_RRR: bool = False, + n_mu_bins: int = 120, + position_type: Literal["rdd", "xyz", "pos"] = "pos", + xi_cut_s: float = 250, xi_refinement_iterations: int = 10, + boxsize: float | None = None, + shot_noise_rescaling1: float = 1, + skip_s_bins: int | tuple[int, int] = 0, skip_l: int = 0, + exclude_samebins: bool = True, exclude_odd_l: bool = False, + sampling_grid_size: int = 301, coordinate_scaling: float = 1, seed: int | None = None, + start_integral_index: Literal[1, 2] | None = None, last_integral_index: Literal[1, 2] | None = None, + verbose: bool = False) -> dict[str, npt.NDArray[np.float64]]: + r""" + Run the 3-point correlation function covariance integration. + Only supports single tracer in "accumulated" Legendre mode without jackknives. + + Parameters + ---------- + mode : string + Choice of binning setup, only one is supported now: + + - ``"legendre_accumulated"``: compute covariance of the 3-point correlation function Legendre multipoles in pairs separation (s) bins accumulated directly, without first doing :math:`\mu`-binned counts. + + s_edges : Numpy array of floats + Separation/radial bin edges. + + max_l : integer + Max Legendre multipole index (required in Legendre ``mode``\s). + Can be odd. + + boxsize : None or float + Periodic box side (one number — so far, only cubic boxes are supported). + All the coordinates need to be between 0 and ``boxsize``. + If None (default) or 0, assumed aperiodic. + + position_type : string, default="pos" + Type of input positions, one of: + + - "rdd": RA, Dec (both in degrees), distance; shape (3, N) + - "xyz": Cartesian positions, shape (3, N) + - "pos": Cartesian positions, shape (N, 3). + + randoms_positions1 : array of floats, shaped according to `position_type` + Coordinates of random points for the first tracer. + + randoms_weights1 : array of floats of length N_randoms + Weights of random points for the first tracer. + If they were made negative for ENCORE, make them positive again before passing to RascalC. + + no_data_galaxies1 : float + Number of first tracer data (not random!) points for the covariance rescaling. + See the definition nuances under ``effective_no_def`` below and set the value of that parameter accordingly. + + effective_no_def: boolean + (Optional) whether to use the effective number of random particles (:math:`\sum(w_i)^2/\sum(w_i^2)`, taking into account weights but invariant to overall weight rescaling) instead of a simple count for the default (Gaussian) shot noise. Should match the way the ``no_data_galaxies1`` value was computed externally. Default False (simple count). + This definition does not affect the final results obtained with shot noise rescaling based on a mock-based sample covariance as a reference, but it seems to provide a slightly better default (Gaussian) shot noise for non-uniform weights. + + RRR_counts : Numpy array of floats, or None + (Optional) RRR (random triplet) counts in ENCORE format. You should be able to load them via `np.genfromtxt("....r_3pcf.txt", skip_header=8)`. Removing the first/ell column is not necessary, the code should be able to do this. + If provided and normalized_RRR=False (default), need to be computed with the same set of randoms, positions and weights, without splits or weight rescaling. If normalized_RRR=True, you need to normalize the RRR counts by the cube of the sum of weights of the randoms used to compute them (survey volume still needs to be the same, but the number of randoms and scaling of their weights can change between the ENCORE RRR and RascalC computations). Otherwise, the covariance normalization will be off. + If not provided and the data is not in a periodic box, triple counts will be estimated with importance sampling (expect longer runtime). + In case of periodic box, the RRR counts are not needed because they are trivial. + + normalized_RRR: boolean + (Optional) whether the RRR counts are normalized by the cube of the sum of weights. Default False (not normalized). If True, the code will renormalize the random weights for consistency. + + n_mu_bins : integer + (Optional) number of angular (mu) bins for the RRR (random triplet) counts computation. Default 120. + If the RRR counts are provided in ENCORE format, this parameter has no effect. + + xi_table_11 : :class:`pycorr.TwoPointEstimator` or :class:`lsstypes.Count2Correlation`, or sequence (tuple or list) of 3 elements: ``(s_values, mu_values, xi_values)``, or sequence (tuple or list) of 4 elements: ``(s_values, mu_values, xi_values, s_edges)`` + Table of first tracer auto-correlation function in separation (s) and :math:`\mu` bins. + The code will use it for interpolation in the covariance matrix integrals. + Important: if the given correlation function is an average in :math:`(s, \mu)` bins, the separation bin edges need to be provided (and the :math:`\mu` bins are assumed to be linear) for rescaling procedure which ensures that the interpolation results averaged over :math:`(s, \mu)` bins returns the given correlation function. In case of ``pycorr.TwoPointEstimator`` or :class:`lsstypes.Count2Correlation`, the edges will be recovered automatically. Unwrapped estimators (:math:`\mu` from -1 to 1) are preferred, because symmetry allows to fix some issues. + In the sequence format: + + - ``s_values`` must be a 1D array of reference separation (s) values for the table, of length N; + - ``mu_values`` must be a 1D array of reference :math:`\mu` values (covering the range from 0 to 1) for the table, of length M; + - ``xi_values`` must be an array of correlation function values at those :math:`(s, \mu)` values of shape (N, M); + - ``s_edges``, if given, must be a 1D array of separation bin edges of length N+1. The bins must come close to zero separation (say start at ``s <= 0.01``). + + The sequence containing 3 elements should be used for theoretical models evaluated at a grid of s, mu values. + The 4-element format should be used for bin-averaged estimates (unless they are in a :class:`pycorr.TwoPointEstimator` or :class:`lsstypes.Count2Correlation`). + + xi_cut_s : float + (Optional) separation value beyond which the correlation function is assumed to be zero for the covariance matrix integrals. Default: 250. + Between the maximum separation from ``xi_table``\s and ``xi_cut_s``, the correlation function is extrapolated as :math:`\propto s^{-4}`. + + xi_refinement_iterations : integer + (Optional) number of iterations in the correlation function refinement procedure for interpolation inside the code, ensuring that the bin-averaged interpolated values match the binned correlation function estimates. Default: 10. + Important: the refinement procedure is disabled completely regardless of this setting if the ``xi_table``\s are sequences of 3 elements, since they are presumed to be a theoretical model evaluated at a grid of s, mu values and not bin averages. + + nthread : integer + Number of (OpenMP) threads to use. + Can not utilize more threads than ``n_loops``. + + - IMPORTANT: AVOID multi-threading in the Python process calling this function (e.g. at NERSC this would mean not setting ``OMP_*`` and other ``*_THREADS`` environment variables; the code should be able to set them by itself). Otherwise the code may run effectively single-threaded. If you need other multi-threaded calculations, run them separately or spawn sub-processes. + + N2 : integer + Number of secondary points to sample per each primary random point. + Setting too low (below 5) is not recommended. + + N3 : integer + Number of tertiary points to sample per each secondary point. + Setting too low (below 5) is not recommended. + + N4 : integer + Number of quaternary points to sample per each tertiary point. + Setting too low (below 5) is not recommended. + + N5 : integer + Number of quinary points to sample per each quaternary point. + Setting too low (below 5) is not recommended. + + N6 : integer + Number of senary points to sample per each quinary point. + Setting too low (below 5) is not recommended. + + n_loops : integer + Number of integration loops. + For optimal balancing and minimal idle time, should be a few times (at least twice) ``nthread`` and exactly divisible by it. + The runtime roughly scales as the number of quads per the number of threads, :math:`\mathcal{O}`\(``N_randoms * N2 * N3 * N4 * n_loops / nthread``). + For reference, on NERSC Perlmutter CPU half-node the code processed about 27 millions (``2.7e7``) quads per second per thread (using 64 threads on half a node) as of December 2024. (In Legendre projected mode, which is probably the slowest, with ``N2 = 5``, ``N3 = 10``, ``N4 = 20``.) + In single-tracer mode, the number of quads is ``N_randoms * N2 * N3 * N4 * n_loops``. + In two-tracer mode, the number of quads is ``(5 * N_randoms1 + 2 * N_randoms2) * N2 * N3 * N4 * n_loops``. + + loops_per_sample : integer + Number of loops to merge into one output sample. + Must divide ``max_loops``. + Recommended to keep the number of samples = ``n_loops / loops_per_sample`` roughly between 10 and 30. + + out_dir : string + Directory for important outputs. + Moderate disk space required (up to a few hundred megabytes), but increases with covariance matrix size and number of samples (see above). + + tmp_dir : string + Directory for temporary files. Can be deleted after running the code, and is cleaned up after the normal execution. + More disk space required - needs to store all the input arrays in the current implementation. + + shot_noise_rescaling1 : float + (Optional) shot-noise rescaling value for the first tracer if known beforehand. Default 1 (no rescaling). + + skip_s_bins : integer or tuple of two integers + (Optional) removal of separations bins at the post-processing stage. + First (or the only) number sets the number of radial/separation bins to skip from the beginning (lowest-separation bins tend to converge worse and probably will not be precise due to the limitations of the formalism). + Second number (if provided) sets the number of radial/separation bins to skip from the end. + By default, no bins are skipped at the post-processing stage. + + skip_l : integer + (Only for the Legendre modes; optional) number of highest (even) multipoles to skip at the post-processing stage. (Higher multipole moments of the correlation function tend to converge worse.) Default 0 (no skipping). + + exclude_samebins : bool + Whether to exclude the same-bin pairs. Default (True) corresponds to ENCORE convention. + + exclude_odd_l : boolean + (Optional) If True, the covariance will exclude the odd multipoles; note that then they will also not count in ``skip_l``. By default (False value), odd multipoles are kept and counted in ``skip_l``. + + seed : integer or None + (Optional) If given as an integer, sets the base RNG (random number generator) seed, allowing to reproduce the results with the same input data and settings (except the number of threads, which can be varied). + Individual subsamples (and accordingly the intrinsic precision estimates and convergence test results) may differ (slightly) because they are accumulated/written in order of loop completion which may depend on external factors at runtime, but the final integrals should be the same. + If None (default), the initialization will be random. The randomly generated seed value can be found afterwards in the log file and/or output after ``the base RNG seed is``, but using the same seed might not reproduce the runs without a preset seed using the code before commit fd2d2c41 (12 June 2025). + Note that False in Python is equivalent to 0, which is a legitimate RNG seed, and therefore falls under the integer case, not like None. (True is equivalent to 1.) + + start_integral_index : integer (1 or 2) or None + (Optional) If given as an integer, chooses from which of 7 integrals (numbered 1, 2) to start. + This parameter is intended to help complete the timed-out runs by skipping the integral computed in an unfinished run. + Using this setting will disable post-processing for safety, because not all integrals may be present in the output directory. Carefully merge the output directories with all the integrals (if applicable) and invoke the post-processing. + + last_integral_index : integer (1 or 2) or None + (Optional) If given as an integer, chooses at which of 2 integrals (numbered 1, 2) to stop. + This parameter can be used to fit a run into a time limit if the full 2 integrals are expected to take longer. + Using this setting will disable post-processing for safety, because not all integrals may be present in the output directory. Carefully merge the output directories with all the integrals (if applicable) and invoke the post-processing. + + sampling_grid_size : integer + (Optional) first guess for the sampling grid size. + The code should be able to find a suitable number automatically. + + verbose : bool + (Optional) report each 5% of each loop's progress by printing. Default False (off). + This can be a lot of output, only use when the number of loops is small. + + coordinate_scaling : float + (Optional) scaling factor for all the Cartesian coordinates. Default 1 (no rescaling). + This option is supported by the C++ code, but its use cases are not very clear. + Zero or negative value is reset to ``boxsize``, rescaling an unit cube to full periodicity. + + Returns + ------- + post_processing_results : dict[str, npt.NDArray[np.float64]] + Post-processing results as a dictionary with string keys and Numpy array values. All this information is also saved in a ``Rescaled_Covariance_Matrices*.npz`` file in the output directory. + Selected common keys are: ``"full_theory_covariance"`` for the final covariance matrix and ``"shot_noise_rescaling"`` for the shot-noise rescaling value(s). + There will also be a ``Raw_Covariance_Matices*.npz`` file in the output directory (as long as the C++ code has run without errors), which can be post-processed separately in a different way using e.g. :func:`RascalC.post_process_auto`. + For convenience, in the output dictionary only, ``"filename"`` contains the exact name of this file, and ``"path"`` contains path to it (also obtainable by :func:`os.path.join`-ing ``out_dir`` with the filename) + """ + + if mode not in ("legendre_accumulated",): raise ValueError("Given mode not supported") + + if not isinstance(max_l, int): raise TypeError("Max ell must be an integer") + if max_l < 0: raise ValueError("Max ell must not be negative") + + # Set some other flags + periodic = bool(boxsize) # False for None (default) and 0 + + if periodic and boxsize < 0: raise ValueError("Periodic box size must be positive") + + max_integrals = 2 + if start_integral_index is not None: + if not isinstance(start_integral_index, int): raise TypeError("start_integral_index must be an integer") + if not (1 <= start_integral_index <= max_integrals): raise ValueError(f"start_integral_index must take a value between 1 and {max_integrals}") + if last_integral_index is not None: + if not isinstance(last_integral_index, int): raise TypeError("last_integral_index must be an integer") + if not (1 <= last_integral_index <= max_integrals): raise ValueError(f"last_integral_index must take a value between 1 and {max_integrals}") + if start_integral_index is not None and last_integral_index < start_integral_index: raise ValueError("last_integral_index must be greater or equal than start_integral_index") + + if n_loops % loops_per_sample != 0: raise ValueError("The sample collapsing factor must divide the number of loops") + + if not isinstance(RRR_counts, (np.ndarray, type(None))): raise TypeError("RRR_counts must be a Numpy array or None") + + if not isinstance(seed, (int, type(None))): raise TypeError("Seed must be an integer or None") + + if no_data_galaxies1 <= 0: raise ValueError("Number of data galaxies (no_data_galaxies1) must be positive") + + ntracers = 1 + ncorr = ntracers * (ntracers + 1) // 2 + suffixes_tracer = suffixes_tracer_all[:ntracers] + indices_corr = indices_corr_all[:ncorr] + suffixes_corr = suffixes_corr_all[:ncorr] + + n_r_bins = len(s_edges) - 1 + + # this condition should be updated for 3PCF, but it is not top priority - for 2PCF, it is only relevant in small boxes, relative to s_max and xi_cut_s + # if periodic and 2 * (max(s_edges) + xi_cut_s) > boxsize: + # warn("Some of the interparticle distances may not be correctly periodically wrapped because of the small box period, so some 6-point configurations may be missed in error. To avoid this, keep the sum of s_max (maximum separation in the covariance bins) and the xi cutoff scale smaller than half of the box size.") + + # set the technical filenames + input_filenames = [os.path.join(tmp_dir, str(t) + ".txt") for t in range(ntracers)] + cornames = [os.path.join(out_dir, f"xi/xi_{index}.dat") for index in indices_corr] + + # make sure the dirs exist + # os.makedirs() will become confused if the path elements to create include pardir (eg. “..” on UNIX systems). + # os.path.abspath should prevent this + out_dir_safe = os.path.abspath(out_dir) + os.makedirs(out_dir_safe, exist_ok = True) + os.makedirs(os.path.join(out_dir_safe, "xi"), exist_ok = True) + os.makedirs(os.path.join(out_dir_safe, "weights"), exist_ok = True) + + # before creating the temporary directory, find which of its parent directories do not exist (yet), to remove them once they are not needed, but only if they are empty + tmp_max_iterations = 100 + tmp_dirs_to_clean_up = [] + tmp_path = os.path.abspath(tmp_dir) # start walking up the path from the temporary directory; this should also remove the slash at the end of tmp_dir if it was there + for _ in range(tmp_max_iterations): # this could be a "while not" loop, but I wanted to prevent a possibility of an endless/unreasonably long loop + if os.path.exists(tmp_path): break # terminate the loop when we found a pre-existing directory + tmp_dirs_to_clean_up.append(tmp_path) # if we reach this, this directory does not exist, add it to the cleanup list + tmp_path = os.path.dirname(tmp_path) # should always jump to the parent directory + else: # this is when the loop executes all tmp_max_iterations iterations without encountering the break statement, i.e. potentially becomes infinite even though it shouldn't + tmp_dirs_to_clean_up = [os.path.abspath(tmp_dir)] # reset the cleanup list to the temporary directory only (the old default behavior) + warn(f"The check of the newly created temporary directory path components failed to finish in {tmp_max_iterations} steps. This is not critical, but should not happen. Please inform the RascalC maintainer of this occurrence.") + os.makedirs(os.path.abspath(tmp_dir), exist_ok = True) # finally, create the temporary dir and all its missing parents + + # Create a log file in output directory + logfilename = "log.txt" + logfile = os.path.join(out_dir, logfilename) + + def print_log(s: object) -> None: + os.system(f"echo \"{s}\" >> {logfile}") + + def print_and_log(s: object) -> None: + print(s) + print_log(s) + + print_and_log(f"Mode: {mode}") + print_and_log(f"Periodic box: {periodic}") + if periodic: print_and_log(f"Box side: {boxsize}") + print_and_log(datetime.now()) + + ndata = [no_data_galaxies1] + + print_and_log(f"Number(s) of data galaxies: {ndata}") + + # write the xi file(s); need to set the 2PCF binning (even if only technical) and decide whether to rescale the 2PCF in the C++ code + all_xi = (xi_table_11,) + xi_s_edges = None + xi_n_mu_bins = None + refine_xi = False + for c, xi in enumerate(all_xi[:ncorr]): + if c > 0: + if type(xi) != type(all_xi[0]): raise TypeError(f"xi_table_{indices_corr[c]} must have the same type as xi_table_11") + if xi.shape != all_xi[0].shape: raise ValueError(f"xi_table_{indices_corr[c]} must have the same shape as xi_table_11") + if isinstance(xi, (pycorr.twopoint_estimator.BaseTwoPointEstimator, lsstypes.Count2Correlation)): + refine_xi = True + xi = fix_and_wrap_allcounts(xi) # try to fix and wrap to |mu| if needed + if c == 0: + xi_n_mu_bins = len(get_mu_edges_from_allcounts(xi)) - 1 + xi_s_edges = get_s_edges_from_allcounts(xi) + elif not np.allclose(xi_s_edges, get_s_edges_from_allcounts(xi)): raise ValueError("Different binning for different correlation functions not supported") + if not np.allclose(get_mu_edges_from_allcounts(xi), np.linspace(0, 1, xi_n_mu_bins + 1)): raise ValueError(f"xi_table_{indices_corr[c]} mu binning is not consistent with linear between 0 and 1 (after wrapping)") + write_xi_file(cornames[c], get_s_avg_from_allcounts(xi), get_mu_avg_from_allcounts(xi), allcount_switch_function(xi, get_input_xi_from_pycorr, get_input_xi_from_lsstypes)) + elif isinstance(xi, Iterable): + if len(xi) == 4: # the last element is the edges + refine_xi = True + if c == 0: xi_s_edges = xi[-1] + elif not np.allclose(xi_s_edges, xi[-1]): raise ValueError("Different binning for different correlation functions not supported") + xi = xi[:-1] + if len(xi) != 3: raise ValueError(f"xi_table {indices_corr[c]} must have 3 or 4 elements if a tuple/list") + r_vals, mu_vals, xi_vals = xi + if len(xi_vals) != len(r_vals): raise ValueError(f"xi_values {indices_corr[c]} must have the same number of rows as r_values") + if len(xi_vals[0]) != len(mu_vals): raise ValueError(f"xi_values {indices_corr[c]} must have the same number of columns as mu_values") + if c == 0: + xi_n_mu_bins = len(mu_vals) + if not refine_xi: + xi_s_edges = (r_vals[:-1] + r_vals[1:]) / 2 # middle values as midpoints of r_vals to be safe + xi_s_edges = np.concatenate([[1e-4], xi_s_edges, [2 * r_vals[-1] - xi_s_edges[-1]]]) # set the lowest edge near 0 and the highest beyond the last point of r_vals + write_xi_file(cornames[c], r_vals, mu_vals, xi_vals) + else: raise TypeError(f"Xi table {indices_corr[c]} must be either a pycorr.TwoPointEstimator, a lsstypes.Count2Correlation, or a tuple/list") + xi_refinement_iterations *= refine_xi # True is 1; False is 0 => 0 iterations => no refinement + + # write the randoms file(s) + randoms_positions = [randoms_positions1] + randoms_weights = [randoms_weights1] + for t, input_filename in enumerate(input_filenames): + randoms_properties = pycorr.twopoint_counter._format_positions(randoms_positions[t], mode = "smu", position_type = position_type, dtype = np.float64) # list of x, y, z coordinate arrays; weights (and jackknife region numbers if any) will be appended + randoms_positions[t] = np.array(randoms_properties) # save the formatted positions as an array for correction function computation + nrandoms = len(randoms_properties[0]) + if randoms_weights[t].ndim != 1: raise ValueError(f"Weights of randoms {t+1} not contained in a 1D array") + if len(randoms_weights[t]) != nrandoms: raise ValueError(f"Number of weights for randoms {t+1} mismatches the number of positions") + if normalized_RRR: randoms_weights[t] /= np.sum(randoms_weights[t]) + randoms_properties.append(randoms_weights[t]) + np.savetxt(input_filename, np.column_stack(randoms_properties)) + del randoms_properties # free memory + + # write the binning files + binfile = os.path.join(out_dir, "radial_binning_cov.csv") + write_binning_file(binfile, s_edges) + binfile_cf = os.path.join(out_dir, "radial_binning_corr.csv") + write_binning_file(binfile_cf, xi_s_edges) + + # deal with RRR counts + if RRR_counts is None: + if periodic: RRR_filename = None # RRR counts not needed + else: # need to run triple_counts + # Select the executable name + exec_name = "bin/triple.s_mu" + "_periodic" * periodic + "_verbose" * verbose + # the above must be true relative to the script location + # below we should make it absolute, i.e. right regardless of the working directory + exec_path = os.path.join(os.path.realpath(os.path.dirname(__file__)), exec_name) + + # form the command line + command = "env OMP_PROC_BIND=spread OMP_PLACES=threads " # set OMP environment variables, should not be set before + command += f"{exec_path} -output {out_dir}/ -nside {sampling_grid_size} -rescale {coordinate_scaling} -nthread {nthread} -maxloops {n_loops} -N2 {N2} -N3 {N3} -xicut {xi_cut_s} -binfile {binfile} -binfile_cf {binfile_cf} -mbin_cf {xi_n_mu_bins} -cf_loops {xi_refinement_iterations}" # here are universally acceptable parameters + command += "".join([f" -in{suffixes_tracer[t]} {input_filenames[t]}" for t in range(ntracers)]) # provide all the random filenames + # don't ask to delete the random files here to be able to reuse them for the main computation later + command += "".join([f" -norm{suffixes_tracer[t]} {ndata[t]}" for t in range(ntracers)]) # provide all ndata for normalization + command += "".join([f" -cor{suffixes_corr[c]} {cornames[c]}" for c in range(ncorr)]) # provide all correlation functions + command += f" -mbin {n_mu_bins}" + if periodic: # append periodic flag and box size + command += f" -perbox -boxsize {boxsize}" + + # deal with the seed + if seed is not None: # need to pass to the C++ code and make sure it can be received properly. 0 (False) is not equivalent to None in this case + seed &= 2**32 - 1 # this bitwise AND truncates the seed into a 32-bit unsigned (positive) integer (definitely a subset of unsigned long) + command += f" -seed {seed}" + + # run the triple_counts code + print_and_log(datetime.now()) + print_and_log(f"Launching the triple_counts C++ code with command: {command}") + status = os.system(f"bash -c 'set -o pipefail; stdbuf -oL -eL {command} 2>&1 | tee -a {logfile}'") + # tee prints what it gets to stdout AND saves to file + # stdbuf -oL -eL should solve the output delays due to buffering without hurting the performance too much + # without pipefail, the exit_code would be of tee, not reflecting main command failures + # feed the command to bash because on Ubuntu it was executed in sh (dash) where pipefail is not supported + + # check the run status + exit_code = os.waitstatus_to_exitcode(status) # assumes we are in Unix-based OS; on Windows status is the exit code + if exit_code: raise RuntimeError(f"The triple_counts C++ code terminated with an error: exit code {exit_code}") + print_and_log("The triple_counts C++ code finished succesfully") + + RRR_filename = f"{out_dir}/RRR_counts_n{n_r_bins}_m{n_mu_bins}_full.txt" + + inv_phi_filename = compute_3pcf_correction_function(randoms_positions[0], randoms_weights[0], binfile, out_dir, periodic, RRR_filename, print_function=print_and_log) + else: # need to convert RRR counts from ENCORE/CADENZA format + inv_phi_filename = compute_3pcf_correction_function_from_encore(randoms_positions[0], randoms_weights[0], binfile, out_dir, RRR_counts, print_function=print_and_log) + del randoms_positions, randoms_weights # free memory + + # Select the executable name + exec_name = "bin/cov.3pcf_" + mode + "_periodic" * periodic + "_verbose" * verbose + # the above must be true relative to the script location + # below we should make it absolute, i.e. right regardless of the working directory + exec_path = os.path.join(os.path.realpath(os.path.dirname(__file__)), exec_name) + + # form the command line + command = "env OMP_PROC_BIND=spread OMP_PLACES=threads " # set OMP environment variables, should not be set before + command += f"{exec_path} -output {out_dir}/ -nside {sampling_grid_size} -rescale {coordinate_scaling} -nthread {nthread} -maxloops {n_loops} -loopspersample {loops_per_sample} -N2 {N2} -N3 {N3} -N4 {N4} -N5 {N5} -N6 {N6} -xicut {xi_cut_s} -binfile {binfile} -binfile_cf {binfile_cf} -mbin_cf {xi_n_mu_bins} -cf_loops {xi_refinement_iterations}" # here are universally acceptable parameters + command += "".join([f" -in{suffixes_tracer[t]} {input_filenames[t]}" for t in range(ntracers)]) # provide all the random filenames + command += " -delete_in" # ask the code to delete the random particle input files after reading them for easy cleanup + command += "".join([f" -norm{suffixes_tracer[t]} {ndata[t]}" for t in range(ntracers)]) # provide all ndata for normalization + if effective_no_def: command += " -effective_norm" # ask the code to use the effective number definition on randoms for shot noise normalization instead of the simple counting + command += "".join([f" -cor{suffixes_corr[c]} {cornames[c]}" for c in range(ncorr)]) # provide all correlation functions + command += f" -max_l {max_l}" + command += f" -phi_file {inv_phi_filename}" + if periodic: # append periodic flag and box size + command += f" -perbox -boxsize {boxsize}" + if start_integral_index: command += f" -start_integral_index {start_integral_index}" + if last_integral_index: command += f" -last_integral_index {last_integral_index}" + + # deal with the seed + if seed is not None: # need to pass to the C++ code and make sure it can be received properly. 0 (False) is not equivalent to None in this case + seed &= 2**32 - 1 # this bitwise AND truncates the seed into a 32-bit unsigned (positive) integer (definitely a subset of unsigned long) + command += f" -seed {seed}" + + # run the main code + print_and_log(datetime.now()) + print_and_log(f"Launching the C++ code with command: {command}") + status = os.system(f"bash -c 'set -o pipefail; stdbuf -oL -eL {command} 2>&1 | tee -a {logfile}'") + # tee prints what it gets to stdout AND saves to file + # stdbuf -oL -eL should solve the output delays due to buffering without hurting the performance too much + # without pipefail, the exit_code would be of tee, not reflecting main command failures + # feed the command to bash because on Ubuntu it was executed in sh (dash) where pipefail is not supported + + # clean up + for tmp_path in tmp_dirs_to_clean_up: rmdir_if_exists_and_empty(tmp_path) # safely remove the temporary directory and all its parents that did not exist before, but only if they are empty + + # check the run status + exit_code = os.waitstatus_to_exitcode(status) # assumes we are in Unix-based OS; on Windows status is the exit code + if exit_code: raise RuntimeError(f"The C++ code terminated with an error: exit code {exit_code}") + print_and_log("The C++ code finished succesfully") + + # post-processing + print_and_log(datetime.now()) + if (start_integral_index is not None and start_integral_index > 1) or (last_integral_index is not None and last_integral_index < max_integrals): + print_and_log("Skipping post-processing, because not all integrals may be present in the output directory.") + print_and_log("Carefully merge the output directories with all the integrals (if applicable) and invoke the post-processing with e.g. RascalC.post_process_auto()") + return + print_and_log("Starting post-processing") + results = post_process_3pcf(out_dir, n_r_bins, max_l, out_dir, shot_noise_rescaling1, skip_s_bins, skip_l, exclude_samebins=exclude_samebins, exclude_odd_l=exclude_odd_l, print_function=print_and_log) + # TODO: + # add mock post-processing + + print_and_log("Finished post-processing") + print_and_log(datetime.now()) + + print_and_log("Performing an extra convergence check") + convergence_check_extra(results, print_function = print_and_log) + + print_and_log("Finished.") + print_and_log(datetime.now()) + return results \ No newline at end of file diff --git a/RascalC/post_process/auto.py b/RascalC/post_process/auto.py index 22b9a28d..ebf48a0f 100644 --- a/RascalC/post_process/auto.py +++ b/RascalC/post_process/auto.py @@ -219,7 +219,7 @@ def post_process_auto(file_root: str, if len(raw_filenames) > 1: warn(f"Found multiple `max_l` options, chosen {max_l=}") else: prefix = f"n{n_r_bins}_l" - raw_matrices = collect_raw_covariance_matrices(file_root, dry_run=True, check_finished=check_finished, two_tracers=two_tracers, print_function=print_function) # safer to skip the actual collection at this stage + raw_matrices = collect_raw_covariance_matrices(file_root, dry_run=True, check_finished=check_finished, two_tracers=two_tracers, print_function=print_function)[0] # safer to skip the actual collection at this stage matched_labels = [label[len(prefix):] for label in raw_matrices.keys() if label.startswith(prefix)] if not matched_labels: if not dry_run: raise ValueError("No Legendre results matched by the number of radial bins.") diff --git a/RascalC/post_process/default.py b/RascalC/post_process/default.py index 96a96033..ab124a9e 100644 --- a/RascalC/post_process/default.py +++ b/RascalC/post_process/default.py @@ -30,7 +30,7 @@ def post_process_default(file_root: str, n: int, m: int, outdir: str, alpha: flo # Check matrix convergence check_eigval_convergence(c2, c4, alpha, print_function=print_function) - # Compute full covariance matrices and precision + # Compute full covariance matrix full_cov = add_cov_terms_single(c2, c3, c4, alpha) # Check positive definiteness diff --git a/RascalC/post_process/default_mocks.py b/RascalC/post_process/default_mocks.py index aaedd8c3..09752241 100644 --- a/RascalC/post_process/default_mocks.py +++ b/RascalC/post_process/default_mocks.py @@ -42,7 +42,7 @@ def post_process_default_mocks(mock_cov_file: str, file_root: str, n: int, m: in # Check matrix convergence for the optimal alpha: if it is <1, the eigenvalue criterion should be strengthened if eigval_ok and alpha_best < 1: check_eigval_convergence(c2f, c4f, alpha_best, print_function=print_function) - # Compute full covariance matrices and precision + # Compute full covariance matrix full_cov = add_cov_terms_single(c2f, c3f, c4f, alpha_best) # Check positive definiteness diff --git a/RascalC/post_process/legendre.py b/RascalC/post_process/legendre.py index f34bdbbd..0dbdeae8 100644 --- a/RascalC/post_process/legendre.py +++ b/RascalC/post_process/legendre.py @@ -30,7 +30,7 @@ def post_process_legendre(file_root: str, n: int, max_l: int, outdir: str, alpha # Check matrix convergence check_eigval_convergence(c2, c4, alpha, print_function=print_function) - # Compute full covariance matrices and precision + # Compute full covariance matrix full_cov = add_cov_terms_single(c2, c3, c4, alpha) # Check positive definiteness diff --git a/RascalC/post_process/legendre_mocks.py b/RascalC/post_process/legendre_mocks.py index 1b65442d..082afc4f 100644 --- a/RascalC/post_process/legendre_mocks.py +++ b/RascalC/post_process/legendre_mocks.py @@ -43,7 +43,7 @@ def post_process_legendre_mocks(mock_cov_file: str, file_root: str, n: int, max_ # Check matrix convergence for the optimal alpha: if it is <1, the eigenvalue criterion should be strengthened if eigval_ok and alpha_best < 1: check_eigval_convergence(c2f, c4f, alpha_best, print_function=print_function) - # Compute full covariance matrices and precision + # Compute full covariance matrix full_cov = add_cov_terms_single(c2f, c3f, c4f, alpha_best) # Check positive definiteness diff --git a/RascalC/post_process/utils.py b/RascalC/post_process/utils.py index 5d2155f0..301fedef 100644 --- a/RascalC/post_process/utils.py +++ b/RascalC/post_process/utils.py @@ -60,6 +60,13 @@ def cov_filter_legendre_pycorr(n: int, max_l: int, skip_r_bins: int = 0, skip_l: return np.ix_(indices_1d, indices_1d) +def apply_cov_filter(A: np.typing.NDArray[np.float64], cov_filter: np.typing.NDArray[np.int_]) -> np.typing.NDArray[np.float64]: + "Apply the 2D cov filter, accounting for the possibility for an extra first dimension for subsamples" + if A.ndim == 2: return A[cov_filter] + if A.ndim == 3: return A[np.arange(A.shape[0])[:, None, None], cov_filter[0][None, :, :], cov_filter[1][None, :, :]] # this correctly selects according to the cov_filter along the last 2 axes + raise ValueError("Matrix dimension must be 2 or 3") + + def load_matrices_single(input_data: dict[str], cov_filter: npt.NDArray[np.int_], tracer: Literal[1, 2] = 1, full: bool = True, jack: bool = False) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]]: """Load the single-tracer covariance matrix terms. Allows to choose the tracer index (1-based) to load.""" joint = "j" * jack + "_" @@ -67,19 +74,10 @@ def load_matrices_single(input_data: dict[str], cov_filter: npt.NDArray[np.int_] c2 = input_data["c2" + joint + suffix] c3 = input_data["c3" + joint + str(tracer) + "," + suffix] c4 = input_data["c4" + joint + str(tracer) * 2 + "," + suffix] - matrices = (c2, c3, c4) + return tuple([symmetrized(apply_cov_filter(a, cov_filter)) for a in (c2, c3, c4)]) - def finalize_matrix(a: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: - return symmetrized(a[cov_filter]) - - if full: # 2D matrices, filter can be applied directly - c234 = [finalize_matrix(a) for a in matrices] - else: # 3D matrices, need to loop over the first index first - c234 = [np.array(list(map(finalize_matrix, a))) for a in matrices] - return tuple(c234) - -def check_eigval_convergence(c2: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], alpha: float = 1, kind: str = "", warn_function: Callable = warn, print_function: Callable = blank_function) -> bool: +def check_eigval_convergence(cN: npt.NDArray[np.float64], c2N: npt.NDArray[np.float64], alpha: float = 1, kind: str = "", Npcf: Literal[2, 3] = 2, warn_function: Callable = warn, print_function: Callable = blank_function) -> bool: """ Perform the eigenvalue convergence test on the covariance matrix terms. Warn about violated condition(s) using the `warn_function` (default: `warnings.warn`). @@ -90,17 +88,17 @@ def check_eigval_convergence(c2: npt.NDArray[np.float64], c4: npt.NDArray[np.flo """ result = True if kind and not kind.endswith(" "): kind += " " - min_eig_c2 = min(np.linalg.eigvalsh(c2)) - min_eig_c4 = min(np.linalg.eigvalsh(c4)) - print_function(f"{kind}4-point covariance matrix convergence check: min eigenvalue of C2 = {min_eig_c2:.2e}, min eigenvalue of C4 = {min_eig_c4:.2e}") - inv_sqrt_c2 = np.linalg.inv(sqrtm(c2)) - min_eig_comb = min(np.linalg.eigvalsh(inv_sqrt_c2.dot(c4).dot(inv_sqrt_c2))) - print_function(f"{kind}4-point covariance matrix convergence check: min eigenvalue of C2^{{-1/2}} C4 C2^{{-1/2}} = {min_eig_comb:.2f}") + min_eig_cN = min(np.linalg.eigvalsh(cN)) + min_eig_c2N = min(np.linalg.eigvalsh(c2N)) + print_function(f"{kind}{2*Npcf}-point covariance matrix convergence check: min eigenvalue of C{Npcf} = {min_eig_cN:.2e}, min eigenvalue of C{2*Npcf} = {min_eig_c2N:.2e}") + inv_sqrt_cN = np.linalg.inv(sqrtm(cN)) + min_eig_comb = min(np.linalg.eigvalsh(inv_sqrt_cN.dot(c2N).dot(inv_sqrt_cN))) + print_function(f"{kind}{2*Npcf}-point covariance matrix convergence check: min eigenvalue of C{Npcf}^{{-1/2}} C{2*Npcf} C{Npcf}^{{-1/2}} = {min_eig_comb:.2f}") if min_eig_comb <= -alpha**2: - warn_function(f"{kind}4-point covariance matrix has not converged properly via the weaker eigenvalue test for shot-noise rescaling >= {alpha:.2f}. Min eigenvalue of C2^{{-1/2}} C4 C2^{{-1/2}} = {min_eig_comb:.2f}, should be > {-alpha**2:.2f}") + warn_function(f"{kind}{2*Npcf}-point covariance matrix has not converged properly via the weaker eigenvalue test for shot-noise rescaling >= {alpha:.2f}. Min eigenvalue of C{Npcf}^{{-1/2}} C{2*Npcf} C{Npcf}^{{-1/2}} = {min_eig_comb:.2f}, should be > {-alpha**2:.2f}") result = False - if min_eig_c4 < - min_eig_c2 * alpha**2: - warn_function(f"{kind}4-point covariance matrix has not converged properly via the stronger eigenvalue test for shot-noise rescaling >= {alpha:.2f}. Min eigenvalue of C2 = {min_eig_c2:.2e}, min eigenvalue of C4 = {min_eig_c4:.2e}") + if min_eig_c2N < - min_eig_cN * alpha**2: + warn_function(f"{kind}{2*Npcf}-point covariance matrix has not converged properly via the stronger eigenvalue test for shot-noise rescaling >= {alpha:.2f}. Min eigenvalue of C{Npcf} = {min_eig_cN:.2e}, min eigenvalue of C{2*Npcf} = {min_eig_c2N:.2e}") result = False return result @@ -144,7 +142,7 @@ def compute_N_eff_D(full_D_est: npt.NDArray[np.float64], print_function = blank_ return N_eff_D -def Psi(alpha: float, c2: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c2s: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64]): +def Psi(alpha: float, c2: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c2s: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: """Compute precision matrix from covariance matrix, removing quadratic order bias terms.""" c_tot = add_cov_terms_single(c2, c3, c4, alpha) partial_covs = add_cov_terms_single(c2s, c3s, c4s, alpha) @@ -152,7 +150,7 @@ def Psi(alpha: float, c2: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], return Psi -def neg_log_L1(alpha: float, target_cov: npt.NDArray[np.float64], c2: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c2s: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64]): +def neg_log_L1(alpha: float, target_cov: npt.NDArray[np.float64], c2: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c2s: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64]) -> float: """Return negative log L1 likelihood between theory and target (data jackknife or mock sample) covariance matrices. log L1 is the Kullback-Leibler divergence with constant terms (including log(det(target_cov))) removed. As a result, the `target_cov` can be a singular matrix. @@ -166,7 +164,7 @@ def neg_log_L1(alpha: float, target_cov: npt.NDArray[np.float64], c2: npt.NDArra return np.trace(np.matmul(Psi_alpha, target_cov)) - logdet[1] -def fit_shot_noise_rescaling(target_cov: npt.NDArray[np.float64], c2: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c2s: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64]): +def fit_shot_noise_rescaling(target_cov: npt.NDArray[np.float64], c2: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c2s: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64]) -> float: """Fit the covariance matrix model to `target_cov` to find the optimal shot-noise rescaling. `target_cov` can be a singular matrix.""" alpha_best = fmin(neg_log_L1, 1., args = (target_cov, c2, c3, c4, c2s, c3s, c4s)) @@ -190,18 +188,12 @@ def load_matrices_multi(input_data: dict[str], cov_filter: npt.NDArray[np.int_], n3 = np.zeros([ntracers] * 3) n4 = np.zeros([ntracers] * 4) - def finalize_matrix(a: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: - return a[cov_filter] - # accumulate the values for matrix_name, matrices in input_data.items(): matrix_name_split = matrix_name.split("_") if len(matrix_name_split) != 2 + full: continue # should skip full if not loading full, and skip subsamples if loading full if full and matrix_name_split[-1] != "full": continue # double-check for safety - if full: # 2D matrix, can apply cov filter directly - matrices = finalize_matrix(matrices) - else: # 3D matrix, should loops over the first index first - matrices = np.array(list(map(finalize_matrix, matrices))) + matrices = apply_cov_filter(matrices, cov_filter) matrix_name = matrix_name_split[0] index = matrix_name_split[1] if matrix_name == "c2" + suffix_jack: diff --git a/RascalC/post_process_3pcf/__init__.py b/RascalC/post_process_3pcf/__init__.py index d4aedce4..e260bb01 100644 --- a/RascalC/post_process_3pcf/__init__.py +++ b/RascalC/post_process_3pcf/__init__.py @@ -1 +1,2 @@ -from .default import post_process_3pcf \ No newline at end of file +from .default import post_process_3pcf +from .mocks import post_process_3pcf_mocks \ No newline at end of file diff --git a/RascalC/post_process_3pcf/default.py b/RascalC/post_process_3pcf/default.py index b318c45c..704f4dec 100644 --- a/RascalC/post_process_3pcf/default.py +++ b/RascalC/post_process_3pcf/default.py @@ -1,116 +1,123 @@ """ -Function to post-process the single-field 3PCF Legendre binned integrals computed by the C++ code. +Function to post-process the single-field 3PCF Legendre binned integrals computed by the C++ code with a given shot-noise rescaling parameter value, alpha. We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. """ import numpy as np import os -from tqdm import trange -from typing import Callable +from ..post_process.utils import check_eigval_convergence, check_positive_definiteness, compute_D_precision_matrix, compute_N_eff_D +from ..raw_covariance_matrices import load_raw_covariances_3pcf_legendre +from .utils import cov_filter_3pcf_legendre, load_matrices, add_cov_terms +from typing import Callable, Iterable -def post_process_3pcf(file_root: str, n: int, max_l: int, n_samples: int, outdir: str, alpha: float = 1, print_function: Callable[[str], None] = print, dry_run: bool = False) -> dict[str]: +def post_process_3pcf(file_root: str, n: int, max_l: int, outdir: str | None = None, alpha: float = 1, overall_scaling: float = 1, skip_r_bins: int | tuple[int, int] = 0, skip_l: int = 0, n_samples: None | int | Iterable[int] | Iterable[bool] = None, exclude_samebins: bool = True, exclude_odd_l: bool = False, check_finished: bool = True, print_function: Callable[[str], None] = print, dry_run: bool = False) -> dict[str]: + r""" + 3PCF post-processing for Legendre (accumulated) mode for a given shot-noise rescaling parameter value, alpha. + + Now it should be safe to run this post-processing while the main RascalC computation is still running, as long as you do not put multiple runs into one output directory. + This is achieved by a default heuristic check for the normal finishing of the main RascalC computation. + With this, inspecting the output of an aborted or timed-out run is harder, by default it will be considered unfinished. But if you are sure that the main computation is not running, you can disable the check via ``check_finished=False``. Doing this once should be sufficient, repeated post-processing attempts should no longer detect the run as unfinished. + + Parameters + ---------- + file_root : string + Path to the RascalC (:func:`RascalC.run_cov_3pcf` or command-line) output directory. + + n : integer + The number of radial bins used in the RascalC run (before applying ``skip_r_bins`` if it is provided). + + max_l : integer + The maximum ell (Legendre moment index) used in the RascalC run (before applying ``skip_l`` if it is provided). + + outdir : string or None + (Optional) path to the directory in which the post-processing results should be saved. If None (default), is set to ``file_root``. Empty string means the current working directory. + We advise to use different output directories for different post-processing options. + + alpha : float + Fixed shot-noise rescaling value to use. In principle optional, but the default value of 1 may not be particularly good. + + overall_scaling : float + Fixed overall scaling value to use. In principle optional, but the default value of 1 may not be particularly good. + + skip_r_bins : integer or tuple of two integers + (Optional) removal of some radial bins. + First (or the only) number sets the number of radial/separation bins to skip from the beginning. + Second number (if provided) sets the number of radial/separation bins to skip from the end. + By default, no bins are skipped. + + skip_l : integer + (Optional) number of higher multipoles to skip (from the end, counting all multipoles by default and only even multipoles if exlude_odd_l is True). + + n_samples : None, integer, array/list/tuple/etc of integers or boolean values + (Optional) selection of RascalC subsamples (independent realizations of Monte-Carlo integrals). + + - If None, use all (default). + - If an integer, use the given number of samples from the beginning. + - If an array/list/tuple/etc of integers, it will be used as a NumPy index array. + - If an array/list/tuple/etc of boolean, it will be used as a NumPy boolean array mask. + + exclude_samebins : boolean + (Optional) If False, the covariance will include the pairs of the same radial bins. + The default behavior (for the True value) is to exclude them for compatibility with ENCORE. + In either case, the post-processed covariances only include each pair of different radial bins in one ordering, ``bin1 < bin2``; the raw covariances also include ``bin1 > bin2`` pairs. + + exclude_odd_l : boolean + (Optional) If True, the covariance will exclude the odd multipoles; note that then they will also not count in ``skip_l``. By default (False value), odd multipoles are kept and counted in ``skip_l``. + + print_function : Callable + (Optional) custom function to use for printing. Default is ``print``. + + dry_run: boolean + (Optional) If True, this will not run actual post-processing, only determine the filename and path (see below). + + Returns + ------- + post_processing_results : dict[str] + Post-processing results as a dictionary with string keys and Numpy array values (mostly). All this information is also saved in a ``Rescaled_Covariance_Matrices*.npz`` file in the ``out_dir`` (in ``file_root`` if the former is not provided). + Selected common keys are: ``"full_theory_covariance"`` for the final covariance matrix and ``"shot_noise_rescaling"`` for the shot-noise rescaling value(s). + For convenience, in the output dictionary only, ``"filename"`` contains the name of the file where the results were saved (which can be inconvenient to predict), and ``"path"`` contains its path (also obtainable by :func:`os.path.join`-ing ``out_dir`` with the filename) + """ + # Set default output directory if not set + if outdir is None: outdir = file_root + output_name = os.path.join(outdir, 'Rescaled_Covariance_Matrices_3PCF_n%d_l%d.npz' % (n, max_l)) name_dict = dict(path=output_name, filename=os.path.basename(output_name)) if dry_run: return name_dict - m = max_l+1 + cov_filter = cov_filter_3pcf_legendre(n, max_l, skip_r_bins, skip_l, exclude_samebins, exclude_odd_l) + + input_file = load_raw_covariances_3pcf_legendre(file_root, n, max_l, n_samples, check_finished=check_finished, print_function=print_function) + # Create output directory if not os.path.exists(outdir): os.makedirs(outdir) - def symmetrize(mat): - """ Add in symmetries to matrices """ - out_mat = np.zeros_like(mat) - for i in range(len(mat)//m): - a = i//n - b = i%n - for j in range(len(mat)//m): - c = j//n - d = j%n - # Add to all relevant bins - these_mat = mat[i*m:(i+1)*m,j*m:(j+1)*m]*0.25 - out_mat[(a*n+b)*m:(a*n+b)*m+m,(c*n+d)*m:(c*n+d)*m+m]+=these_mat - out_mat[(b*n+a)*m:(b*n+a)*m+m,(c*n+d)*m:(c*n+d)*m+m]+=these_mat - out_mat[(b*n+a)*m:(b*n+a)*m+m,(d*n+c)*m:(d*n+c)*m+m]+=these_mat - out_mat[(a*n+b)*m:(a*n+b)*m+m,(d*n+c)*m:(d*n+c)*m+m]+=these_mat - return 0.5*(out_mat+out_mat.T) - - def load_matrices(index): - """Load intermediate or full covariance matrices""" - cov_root = os.path.join(file_root, '3PCFCovMatricesAll/') - c3_0 = np.loadtxt(cov_root+'c3_n%d_l%d_0_%s.txt'%(n,max_l,index)) - c4_0 = np.loadtxt(cov_root+'c4_n%d_l%d_0_%s.txt'%(n,max_l,index)) - c5_0 = np.loadtxt(cov_root+'c5_n%d_l%d_0_%s.txt'%(n,max_l,index)) - c6_0 = np.loadtxt(cov_root+'c6_n%d_l%d_0_%s.txt'%(n,max_l,index)) - c3_1 = np.loadtxt(cov_root+'c3_n%d_l%d_1_%s.txt'%(n,max_l,index)) - c4_1 = np.loadtxt(cov_root+'c4_n%d_l%d_1_%s.txt'%(n,max_l,index)) - c5_1 = np.loadtxt(cov_root+'c5_n%d_l%d_1_%s.txt'%(n,max_l,index)) - c6_1 = np.loadtxt(cov_root+'c6_n%d_l%d_1_%s.txt'%(n,max_l,index)) - - c3 = c3_0+c3_1 - c4 = c4_0+c4_1 - c5 = c5_0+c5_1 - c6 = c6_1 # do not include c6_0 term here - - # Now symmetrize and return matrices - return symmetrize(c3),symmetrize(c4),symmetrize(c5),symmetrize(c6) - # Load in full theoretical matrices print_function("Loading best estimate of covariance matrix") - c3,c4,c5,c6=load_matrices('full') + c3, c4, c5, c6 = load_matrices(input_file, n, max_l, cov_filter, full=True) + + # Check matrix convergence by analogy with 2PCF, may be less helpful + check_eigval_convergence(c3, c6, alpha, Npcf=3, print_function=print_function) + + # Compute full covariance matrix + full_cov = add_cov_terms(c3, c4, c5, c6, alpha) * overall_scaling - # Compute full covariance matrices and precision - full_cov = c6+c5*alpha+c4*alpha**2.+c3*alpha**3. - n_bins = len(c6) + # Check positive definiteness + check_positive_definiteness(full_cov) # Compute full precision matrix print_function("Computing the full precision matrix estimate:") # Load in partial theoretical matrices - c3s, c4s, c5s, c6s = [], [], [], [] - for i in trange(n_samples, desc="Loading full subsamples"): - c3t, c4t, c5t, c6t = load_matrices(i) - c3s.append(c3t) - c4s.append(c4t) - c5s.append(c5t) - c6s.append(c6t) - c3s, c4s, c5s, c6s = [np.array(a) for a in (c3s, c4s, c5s, c6s)] - - partial_cov = alpha**3 * c3s + alpha**2 * c4s + alpha * c5s + c6s - sum_partial_cov = np.sum(partial_cov, axis=0) - - tmp=0. - - for _ in range(1): - for i in range(n_samples): - c_excl_i = (sum_partial_cov - partial_cov[i]) / (n_samples - 1) - try: - tmp+=np.matmul(np.linalg.inv(c_excl_i), partial_cov[i]) - except np.linalg.linalg.LinAlgError: - print_function("Could not invert submatrix, so setting overall precision to zero. Matrix is not fully converged") - tmp = np.inf - if tmp==np.inf: - # Couldn't estimate precision so break loop - full_prec = np.zeros_like(full_cov) - full_D_est = np.zeros_like(full_cov) - N_eff_D = 0. - continue - - full_D_est=(n_samples-1.)/n_samples * (-1.*np.eye(n_bins) + tmp/n_samples) - full_prec = np.matmul(np.eye(n_bins)-full_D_est,np.linalg.inv(full_cov)) - print_function("Full precision matrix estimate computed") - - # Now compute effective N: - slogdetD=np.linalg.slogdet(full_D_est) - D_value = slogdetD[0]*np.exp(slogdetD[1]/n_bins) - if slogdetD[0]<0: - print_function("N_eff is negative! Setting to zero") - N_eff_D = 0. - else: - N_eff_D = (n_bins+1.)/D_value+1. - print_function("Total N_eff Estimate: %.4e"%N_eff_D) - - output_dict = dict(full_theory_covariance=full_cov, shot_noise_rescaling=alpha, + c3s, c4s, c5s, c6s = load_matrices(input_file, n, max_l, cov_filter, full=False) + partial_cov = add_cov_terms(c3s, c4s, c5s, c6s, alpha) * overall_scaling + full_D_est, full_prec = compute_D_precision_matrix(partial_cov, full_cov) + print_function("Full precision matrix estimate computed") + + # Now compute effective N: + N_eff_D = compute_N_eff_D(full_D_est, print_function) + + output_dict = dict(full_theory_covariance=full_cov, shot_noise_rescaling=alpha, overall_scaling=overall_scaling, full_theory_precision=full_prec, N_eff=N_eff_D, full_theory_D_matrix=full_D_est, individual_theory_covariances=partial_cov) diff --git a/RascalC/post_process_3pcf/mocks.py b/RascalC/post_process_3pcf/mocks.py new file mode 100644 index 00000000..954a42e7 --- /dev/null +++ b/RascalC/post_process_3pcf/mocks.py @@ -0,0 +1,209 @@ +""" +Function to post-process the single-field 3PCF Legendre binned integrals computed by the C++ code, obtaining the shot-noise rescaling parameter, alpha, from a mock derived covariance matrix. +We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. +""" + +import numpy as np +import numpy.typing as npt +from scipy.optimize import fmin +import os +from ..utils import format_skip_r_bins +from ..post_process.utils import check_eigval_convergence, check_positive_definiteness, compute_D_precision_matrix, compute_N_eff_D +from ..raw_covariance_matrices import load_raw_covariances_3pcf_legendre +from .utils import cov_filter_3pcf_legendre, load_matrices, add_cov_terms +from typing import Callable, Iterable + + +def Psi(alpha: float, c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c5: npt.NDArray[np.float64], c6: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64], c5s: npt.NDArray[np.float64], c6s: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + """Compute precision matrix from covariance matrix, removing quadratic order bias terms.""" + c_tot = add_cov_terms(c3, c4, c5, c6, alpha) + partial_covs = add_cov_terms(c3s, c4s, c5s, c6s, alpha) + _, Psi = compute_D_precision_matrix(partial_covs, c_tot) + return Psi + + +def neg_log_L1(alpha: float, overall_scaling: float, target_cov: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c5: npt.NDArray[np.float64], c6: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64], c5s: npt.NDArray[np.float64], c6s: npt.NDArray[np.float64]) -> float: + """Return negative log L1 likelihood between 3PCF theory and target (data jackknife or mock sample) covariance matrices. + log L1 is the Kullback-Leibler divergence with constant terms (including log(det(target_cov))) removed. + As a result, the `target_cov` can be a singular matrix. + This function does not allow negative shot-noise rescaling `alpha` or non-positive overall scaling `overall_scaling` by returning infinity.""" + if alpha < 0 or overall_scaling <= 0: return np.inf # negative shot-noise rescaling causes problems and does not make sense + Psi_alpha = Psi(alpha, c3, c4, c5, c6, c3s, c4s, c5s, c6s) / overall_scaling + logdet = np.linalg.slogdet(Psi_alpha) + if logdet[0] < 0: + # Remove any dodgy inversions + return np.inf + return np.trace(np.matmul(Psi_alpha, target_cov)) - logdet[1] + + +def fit_shot_noise_rescaling(target_cov: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c5: npt.NDArray[np.float64], c6: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64], c5s: npt.NDArray[np.float64], c6s: npt.NDArray[np.float64]) -> float: + """Fit the 3PCF covariance matrix model to `target_cov` to find the optimal shot-noise rescaling. + `target_cov` can be a singular matrix.""" + alpha_best = fmin(neg_log_L1, x0=1, args=(1, target_cov, c3, c4, c5, c6, c3s, c4s, c5s, c6s)) + return alpha_best[0] + + +def fit_shot_noise_and_overall_rescaling(target_cov: npt.NDArray[np.float64], c3: npt.NDArray[np.float64], c4: npt.NDArray[np.float64], c5: npt.NDArray[np.float64], c6: npt.NDArray[np.float64], c3s: npt.NDArray[np.float64], c4s: npt.NDArray[np.float64], c5s: npt.NDArray[np.float64], c6s: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + """Fit the 3PCF covariance matrix model to `target_cov` to find the optimal shot-noise rescaling and overall scaling. + `target_cov` can be a singular matrix.""" + return fmin(lambda x: neg_log_L1(x[0], x[1], target_cov, c3, c4, c5, c6, c3s, c4s, c5s, c6s), x0=(1, 1)) + + + +def post_process_3pcf_mocks(mock_cov_file: str, file_root: str, n: int, max_l: int, outdir: str | None = None, overall_scaling: float | None = None, skip_r_bins: int | tuple[int, int] = 0, skip_l: int = 0, n_samples: None | int | Iterable[int] | Iterable[bool] = None, exclude_samebins: bool = True, exclude_odd_l: bool = False, check_finished: bool = True, max_l_mock: int | None = None, print_function: Callable[[str], None] = print, dry_run: bool = False) -> dict[str]: + r""" + 3PCF post-processing for Legendre (accumulated) mode, obtaining the shot-noise rescaling parameter, alpha, from a mock-derived covariance matrix. + + Now it should be safe to run this post-processing while the main RascalC computation is still running, as long as you do not put multiple runs into one output directory. + This is achieved by a default heuristic check for the normal finishing of the main RascalC computation. + With this, inspecting the output of an aborted or timed-out run is harder, by default it will be considered unfinished. But if you are sure that the main computation is not running, you can disable the check via ``check_finished=False``. Doing this once should be sufficient, repeated post-processing attempts should no longer detect the run as unfinished. + + Parameters + ---------- + mock_cov_file : string + Path to the text file containing the mock sample covariance matrix (from ENCORE 3PCF measurements), which will be used to determine the shot-noise rescaling parameter, alpha. The covariance matrix should have the same radial bins as the theoretical matrices computed by RascalC, have the repeated radial bin pairs removed and same-bin pairs removed according to the ``exclude_samebins`` option (typically, yes). It may have a different (higher) maximum ell than the computed ``RascalC`` covariance matrices (in such a case, specify it with max_l_mock). The ``skip_r_bins`` and ``skip_l`` filters will be applied after the mock covariance matrix is loaded. The covariance matrix should be in a text format that can be loaded by ``numpy.loadtxt``. + + file_root : string + Path to the RascalC (:func:`RascalC.run_cov_3pcf` or command-line) output directory. + + n : integer + The number of radial bins used in the RascalC run (before applying ``skip_r_bins`` if it is provided). + + max_l : integer + The maximum ell (Legendre moment index) used in the RascalC run (before applying ``skip_l`` if it is provided). + + outdir : string or None + (Optional) path to the directory in which the post-processing results should be saved. If None (default), is set to ``file_root``. Empty string means the current working directory. + We advise to use different output directories for different post-processing options. + + overall_scaling : float or None + (Optional) If a specific value is provided, it will be used as the fixed overall scaling factor. If None (default), the scaling will be determined by the best fit to the mock covariance matrix (together with the shot-noise rescaling parameter, alpha). + + skip_r_bins : integer or tuple of two integers + (Optional) removal of some radial bins. + First (or the only) number sets the number of radial/separation bins to skip from the beginning. + Second number (if provided) sets the number of radial/separation bins to skip from the end. + By default, no bins are skipped. + + skip_l : integer + (Optional) number of higher multipoles to skip (from the end, counting all multipoles by default and only even multipoles if exlude_odd_l is True). + + n_samples : None, integer, array/list/tuple/etc of integers or boolean values + (Optional) selection of RascalC subsamples (independent realizations of Monte-Carlo integrals). + + - If None, use all (default). + - If an integer, use the given number of samples from the beginning. + - If an array/list/tuple/etc of integers, it will be used as a NumPy index array. + - If an array/list/tuple/etc of boolean, it will be used as a NumPy boolean array mask. + + exclude_samebins : boolean + (Optional) If False, the covariance will include the pairs of the same radial bins. This probably will not work with the mock sample covariance file produced from ENCORE. + The default behavior (for the True value) is to exclude them for compatibility with ENCORE. + In either case, the post-processed covariances only include each pair of different radial bins in one ordering, ``bin1 < bin2``; the raw covariances also include ``bin1 > bin2`` pairs. + + exclude_odd_l : boolean + (Optional) If True, the covariance will exclude the odd multipoles; note that then they will also not count in ``skip_l``. By default (False value), odd multipoles are kept and counted in ``skip_l``. + + max_l_mock : integer or None + (Optional) The maximum ell (Legendre moment index) used in the mock covariance matrix. If None (default), it is assumed to be the same as `max_l` for theoretical covariance matrices. This option is provided in case the mock covariance matrix has a different (higher) maximum ell than the theoretical covariance matrices computed by RascalC. + + print_function : Callable + (Optional) custom function to use for printing. Default is ``print``. + + dry_run: boolean + (Optional) If True, this will not run actual post-processing, only determine the filename and path (see below). + + Returns + ------- + post_processing_results : dict[str] + Post-processing results as a dictionary with string keys and Numpy array values (mostly). All this information is also saved in a ``Rescaled_Covariance_Matrices*.npz`` file in the ``out_dir`` (in ``file_root`` if the former is not provided). + Selected common keys are: ``"full_theory_covariance"`` for the final covariance matrix and ``"shot_noise_rescaling"`` for the shot-noise rescaling value(s). + For convenience, in the output dictionary only, ``"filename"`` contains the name of the file where the results were saved (which can be inconvenient to predict), and ``"path"`` contains its path (also obtainable by :func:`os.path.join`-ing ``out_dir`` with the filename) + """ + # Set default output directory if not set + if outdir is None: outdir = file_root + + output_name = os.path.join(outdir, 'Rescaled_Covariance_Matrices_3PCF_Mocks_n%d_l%d.npz' % (n, max_l)) + name_dict = dict(path=output_name, filename=os.path.basename(output_name)) + if dry_run: return name_dict + + # Load the mock covariance matrix and convert it to the same convention as RascalC theoretical covariance matrices + mock_cov = np.loadtxt(mock_cov_file) # load external mock covariance matrix from text file, would be based on ENCORE + # prepare r bin filtering + skip_r_bins_start, skip_r_bins_end = format_skip_r_bins(skip_r_bins) + r_bin_indices = np.arange(n) + r_bin_index1 = np.repeat(r_bin_indices, n-exclude_samebins-r_bin_indices) + r_bin_index2 = np.concatenate([r_bin_indices[i+exclude_samebins:] for i in r_bin_indices]) + n_r_pairs_orig = len(r_bin_index1) + # r_bin_index1 and r_bin_index2 cover all the bin pairs under the condition r_bin_index1 < r_bin_index2, the order follows the ENCORE format + r_filter = (r_bin_index1 >= skip_r_bins_start) & (r_bin_index1 < n - skip_r_bins_end) & (r_bin_index2 >= skip_r_bins_start) & (r_bin_index2 < n - skip_r_bins_end) # filter for the bin pairs to keep based on the skip_r_bins option + n_r_pairs = r_filter.sum() + # prepare ell indexing and scaling factor accounting for the different basis + ells = np.arange(0, max_l + 1, 1+exclude_odd_l) + if skip_l > 0: ells = ells[:-skip_l] # without the condition, wouldn't work right for skip_l=0 + 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 + # scale and transpose the covariance + if max_l_mock is None: max_l_mock = max_l + mock_cov = mock_cov.reshape(max_l_mock + 1, n_r_pairs_orig, max_l_mock + 1, n_r_pairs_orig) # reshape the covariance from 2D to 4D, the ENCORE ordering is [l, r_bin_pair] for both rows and columns + mock_cov = mock_cov[ells][:, :, ells] # apply the multipole selection + mock_cov = mock_cov.transpose(1, 0, 3, 2) # change ordering to [r_bin_pair, l] for both rows and columns as in RascalC + mock_cov *= ell_factor[:, None, None] * ell_factor[None, None, :] # apply the factor, NumPy broadcasting matches trailing dimensions. need to double-check if it is not division; there might also be a factor of 2 or something similar + mock_cov = mock_cov[r_filter][:, :, r_filter] # apply the r bin pair filter + mock_cov = mock_cov.reshape(n_r_pairs * len(ells), n_r_pairs * len(ells)) # reshape back to 2D, should be ready + + cov_filter = cov_filter_3pcf_legendre(n, max_l, skip_r_bins, skip_l, exclude_samebins, exclude_odd_l) + + input_file = load_raw_covariances_3pcf_legendre(file_root, n, max_l, n_samples, check_finished=check_finished, print_function=print_function) + + # Create output directory + if not os.path.exists(outdir): + os.makedirs(outdir) + + # Load in full theoretical matrices + print_function("Loading best estimate of covariance matrix") + c3, c4, c5, c6 = load_matrices(input_file, n, max_l, cov_filter, full=True) + + # Check matrix convergence by analogy with 2PCF, may be less helpful + eigval_ok = check_eigval_convergence(c3, c6, Npcf=3, print_function=print_function) + + # Load in partial theoretical matrices + c3s, c4s, c5s, c6s = load_matrices(input_file, n, max_l, cov_filter, full=False) + + # Now optimize for shot-noise rescaling parameter alpha, and possibly overall scaling parameter + if overall_scaling is None: + print_function("Optimizing for the shot-noise rescaling and overall scaling parameters") + alpha_best, overall_scaling = fit_shot_noise_and_overall_rescaling(mock_cov, c3, c4, c5, c6, c3s, c4s, c5s, c6s) + print_function("Optimization complete - optimal shot-noise rescaling parameter is %.6f and optimal overall scaling parameter is %.6f" % (alpha_best, overall_scaling)) + else: + print_function("Using the provided overall scaling parameter: %.6f" % overall_scaling) + print_function("Optimizing for the shot-noise rescaling parameter") + alpha_best = fit_shot_noise_rescaling(mock_cov, c3, c4, c5, c6, c3s, c4s, c5s, c6s) + print_function("Optimization complete - optimal rescaling parameter is %.6f" % alpha_best) + + # Check matrix convergence for the optimal alpha: if it is <1, the eigenvalue criterion should be strengthened + if eigval_ok and alpha_best < 1: check_eigval_convergence(c3, c6, alpha_best, Npcf=3, print_function=print_function) + + # Compute full covariance matrix + full_cov = add_cov_terms(c3, c4, c5, c6, alpha_best) * overall_scaling + + # Check positive definiteness + check_positive_definiteness(full_cov) + + # Compute full precision matrix + print_function("Computing the full precision matrix estimate:") + partial_cov = add_cov_terms(c3s, c4s, c5s, c6s, alpha_best) * overall_scaling + full_D_est, full_prec = compute_D_precision_matrix(partial_cov, full_cov) + print_function("Full precision matrix estimate computed") + + # Now compute effective N: + N_eff_D = compute_N_eff_D(full_D_est, print_function) + + output_dict = dict(full_theory_covariance=full_cov, shot_noise_rescaling=alpha_best, overall_scaling=overall_scaling, + full_theory_precision=full_prec, N_eff=N_eff_D, + full_theory_D_matrix=full_D_est, individual_theory_covariances=partial_cov, + mock_covariance=mock_cov) + + np.savez_compressed(output_name, **output_dict) + print_function("Saved output covariance matrices as %s"%output_name) + + return output_dict | name_dict \ No newline at end of file diff --git a/RascalC/post_process_3pcf/utils.py b/RascalC/post_process_3pcf/utils.py new file mode 100644 index 00000000..15377c44 --- /dev/null +++ b/RascalC/post_process_3pcf/utils.py @@ -0,0 +1,48 @@ +import numpy as np +from ..utils import format_skip_r_bins +from ..post_process.utils import apply_cov_filter + + +def cov_filter_3pcf_legendre(n: int, max_l: int, skip_r_bins: int | tuple[int, int] = 0, skip_l: int = 0, exclude_samebins: bool = True, exclude_odd_l: bool = False): + """Produce a 2D indexing array for 3PCF Legendre covariance matrices in RascalC convention (not yet fully compatible with ENCORE bin ordering).""" + skip_r_bins_start, skip_r_bins_end = format_skip_r_bins(skip_r_bins) + n_l = max_l + 1 + l_indices = np.arange(0, n_l, 1+exclude_odd_l) + if skip_l > 0: l_indices = l_indices[:-skip_l] # without the condition, wouldn't work right for skip_l=0 + r_indices = np.arange(skip_r_bins_start, n - skip_r_bins_end) + r_indices1, r_indices2 = [a.ravel() for a in np.meshgrid(r_indices, r_indices, indexing='ij')] # flattened array indices + r_filter = (r_indices1 <= r_indices2 - exclude_samebins) # strictly less for exclude_samebin=True, less or equal otherwise + r_indices1, r_indices2 = r_indices1[r_filter], r_indices2[r_filter] # apply filter to both index arrays + indices_l_r = (n_l * (n * r_indices1 + r_indices2))[:, None] + l_indices[None, :] + # indices_l_r = (n_l * (n * r_indices1 + r_indices2))[None, :] + l_indices[:, None] # could switch the l vs bin pair ordering right here easily but decided not to + indices_1d = indices_l_r.ravel() + return np.ix_(indices_1d, indices_1d) + + +def symmetrized_3pcf(A: np.typing.NDArray[np.float64], n: int, max_l: int) -> np.typing.NDArray[np.float64]: + "Symmetrize a 3PCF covariance matrix (2D array), or an array of 3PCF covariance matrices (3D array)" + if len(A.shape) not in (2, 3): raise ValueError("Dimension of the input array must be 2 or 3") + m = max_l + 1 + if not np.array_equal(A.shape[-2:], [n * n * m] * 2): raise ValueError("Unexpected shape in the last 2 dimensions") + leading_dims = list(A.shape[:-2]) # list containing a leading dimension for the array of covariance matrices, and empty for a single covariance matrix + A = A.reshape(leading_dims + [n, n, m] * 2) # last 6 axes will be [r1, r2, l12, r3, r4, l34] + A = (A + A.swapaxes(-2, -3)) / 2 # symmetrize wrt swaps of r3 and r4 + A = (A + A.swapaxes(-5, -6)) / 2 # symmetrize wrt swaps of r1 and r2 + A = A.reshape(leading_dims + [n * n * m] * 2) # back to the original shape + return (A + A.swapaxes(-1, -2)) / 2 # finally, symmetrize wrt full covariance matrix bin swap + + +def load_matrices(input_data: dict[str], n: int, max_l: int, cov_filter: np.typing.NDArray[np.int_], full: bool = True) -> tuple[np.typing.NDArray[np.float64], np.typing.NDArray[np.float64], np.typing.NDArray[np.float64], np.typing.NDArray[np.float64]]: + """Load the 3PCF single-tracer covariance matrix terms.""" + matrices = [] + for npoints in range(3, 7): + these_matrices = [input_data[f"c{npoints}_{index}" + "_full" * full] for index in range(npoints == 6, 2)] # exclude c6_0 term, because it should be small but is also hard to compute (see Section 5.2.3 and Appendix A of https://arxiv.org/abs/1910.04764) + this_matrix = these_matrices[0] + if npoints != 6: this_matrix += these_matrices[1] + matrices.append(apply_cov_filter(symmetrized_3pcf(this_matrix, n, max_l), cov_filter)) # symmetrize before filtering, because filtering removes repeating bin pairs + return tuple(matrices) + + +def add_cov_terms(c3: np.typing.NDArray[np.float64], c4: np.typing.NDArray[np.float64], c5: np.typing.NDArray[np.float64], c6: np.typing.NDArray[np.float64], alpha: float = 1) -> np.typing.NDArray[np.float64]: + """Add the 3PCF single-tracer covariance matrix terms with a given shot-noise rescaling value.""" + return c6 + c5 * alpha + c4 * alpha**2 + c3 * alpha**3 \ No newline at end of file diff --git a/RascalC/raw_covariance_matrices.py b/RascalC/raw_covariance_matrices.py index 23c6f785..fc29c1e4 100644 --- a/RascalC/raw_covariance_matrices.py +++ b/RascalC/raw_covariance_matrices.py @@ -1,4 +1,4 @@ -## Script to collect the raw covariance matrices from the output directory of the C++ code +"Functions to collect the raw covariance matrices from the output directory of the C++ code into a numpy file, load them, and catenate raw covariance matrices from one directory or multiple directories into another directory." import numpy as np import numpy.typing as npt @@ -15,16 +15,17 @@ def convert_suffix(suffix: str) -> int | str: return suffix -def organize_filename(filename: str, output_groups: dict, jack: bool = False) -> None: +def organize_filename(filename: str, output_groups: dict, jack: bool = False, threepcf: bool = False) -> None: "Interpret the filename as saved by the C++ code and put it into output_groups array. Jack is for jackknife" filename_no_ext = os.path.basename(filename) # remove the directory name filename_no_ext = ".".join(filename_no_ext.split(".")[:-1]) # remove the extension filename_parts = filename_no_ext.split("_") # split the rest by underscore term_name = filename_parts[0] # cN, RRN or EEN - if term_name not in ("c2", "c3", "c4", "RR1", "RR2", "EE1", "EE2"): return # do not process other arrays + term_names = [f"c{i}" for i in range(3, 7)] if threepcf else ("c2", "c3", "c4", "RR1", "RR2", "EE1", "EE2") + if term_name not in term_names: return # do not process other arrays if jack and term_name.startswith("c"): term_name += "j" # add "j" to cN if doing jack output_group_name = "_".join(filename_parts[1:-2]) # nN and mM or lL joined back - indices = filename_parts[-2] # tracer numbers + indices = filename_parts[-2] # tracer numbers for 2PCF and integral index for 3PCF try: suffix = convert_suffix(filename_parts[-1]) # "full" or subsample number except ValueError: @@ -40,14 +41,14 @@ def organize_filename(filename: str, output_groups: dict, jack: bool = False) -> output_groups[output_group_name][matrix_name][suffix] = filename -def save_safe(output_dir: str, output_group_name: str, output_dictionary: dict[str], print_function: Callable[[str], None] = print) -> None: +def save_safe(output_dir: str, output_group_name: str, output_dictionary: dict[str], threepcf: bool = False, print_function: Callable[[str], None] = print) -> None: "Save the dictionary of Numpy arrays into directory avoiding name clashes" - output_filename = os.path.join(output_dir, f"Raw_Covariance_Matrices_{output_group_name}.npz") + output_filename = os.path.join(output_dir, "Raw_Covariance_Matrices" + "_3PCF" * threepcf + f"_{output_group_name}.npz") if os.path.exists(output_filename): print_function(f"The default output filename for group {output_group_name}, {output_filename}, already exists. Will try to find a replacement.") i = 1 while True: - output_filename = os.path.join(output_dir, f"Raw_Covariance_Matrices_{output_group_name}.{i}.npz") + output_filename = os.path.join(output_dir, "Raw_Covariance_Matrices" + "_3PCF" * threepcf + f"_{output_group_name}.{i}.npz") if not os.path.exists(output_filename): print_function(f"Found unused name {output_filename}, will save there.") break @@ -57,15 +58,17 @@ def save_safe(output_dir: str, output_group_name: str, output_dictionary: dict[s np.savez_compressed(output_filename, **output_dictionary) -def collect_raw_covariance_matrices(cov_dir: str, dry_run: bool = False, cleanup: bool = True, check_finished: bool = True, two_tracers: bool | None = None, print_function: Callable[[str], None] = print) -> dict[str, dict[str, npt.NDArray[np.float64]]]: +def collect_raw_covariance_matrices(cov_dir: str, dry_run: bool = False, cleanup: bool = True, check_finished: bool = True, two_tracers: bool | None = None, print_function: Callable[[str], None] = print) -> list[dict[str, dict[str, npt.NDArray[np.float64]]]]: """ Collect the covariance matrices from text files written by the C++ code and organize them into a Numpy .npz file. + Returns a list of two dictionaries, for 2PCF and 3PCF separately, with the structure {output_group_name: {matrix_name: array}}, with {matrix_name: array} replicating the structure of the .npz file. With dry_run enabled, this function will only return the output group names via the dictionary, and not perform the actual collection. With cleanup enabled (default), deletes the text files after collection. With check_finished enabled (default), performs a heuristic check whether the run seems finished by looking for the presence of the full covariance matrices. If they are not found, the collection does not proceed and a warning is issued. This is to prevent disrupting ongoing runs by collecting and/or deleting the text files. If you want to check convergence of timed-out run(s) that did not produce the full matrices, you can disable this check, but be careful not to disrupt ongoing runs. With check_finished enabled, the code additionally checks presence of all types of matrices for two tracers. Pass two_tracers = True to enable this check explicitly, or two_tracers = False to disable it. By default (two_tracers = None), the code will try to guess whether the run is for two tracers by looking for the presence of the xi_22.dat file, which is only produced for two tracers by :func:`RascalC.run_cov`. """ + # 2PCF cov_dir_all = os.path.join(cov_dir, 'CovMatricesAll/') cov_dir_jack = os.path.join(cov_dir, 'CovMatricesJack/') @@ -73,93 +76,111 @@ def collect_raw_covariance_matrices(cov_dir: str, dry_run: bool = False, cleanup # load the full matrices for input_filename in glob(cov_dir_all + "*.txt"): - organize_filename(input_filename, output_groups, jack=False) + organize_filename(input_filename, output_groups, jack=False, threepcf=False) # load the jack matrices if present for input_filename in glob(cov_dir_jack + "*.txt"): - organize_filename(input_filename, output_groups, jack=True) + organize_filename(input_filename, output_groups, jack=True, threepcf=False) - print_function(f"Detected {len(output_groups)} output group(s) in {cov_dir}") + print_function(f"Detected {len(output_groups)} 2PCF output group(s) in {cov_dir}") - if dry_run: return {output_group_name: {} for output_group_name in output_groups} # return the group names via the dictionary, skipping any actual collection + # 3PCF - should be done separately, because the output group name might match 2PCF and some arrays would be overwritten + cov_dir_3pcf_all = os.path.join(cov_dir, '3PCFCovMatricesAll/') - return_dictionary = {} + output_groups_3pcf = {} + + # load the full 3PCF matrices + for input_filename in glob(cov_dir_3pcf_all + "*.txt"): + organize_filename(input_filename, output_groups_3pcf, jack=False, threepcf=True) + + print_function(f"Detected {len(output_groups_3pcf)} 3PCF output group(s) in {cov_dir}") + + if dry_run: return [{output_group_name: {} for output_group_name in output_groups}, {output_group_name: {} for output_group_name in output_groups_3pcf}] # return the group names via the list of dictionaries, skipping any actual collection + + return_dictionaries = [{}, {}] # for 2PCF and 3PCF separately - for output_group_name, output_group in output_groups.items(): - print_function(f"Processing output group {output_group_name}") - - # check if the run seems finished by looking for the presence of the full matrices, if check_finished is enabled - if check_finished: - unfinished_names = [matrix_name for matrix_name, matrix_filenames_dictionary in output_group.items() if "full" not in matrix_filenames_dictionary] - if unfinished_names: - print_function(f"WARNING: {unfinished_names} full matrices not found for output group {output_group_name}. This may indicate that the C++ code is still running. Will skip this output group to avoid disrupting an ongoing run by collecting and/or deleting the text files. If you are sure that the run(s) in this directory finished (e.g., it/they timed out without producing the full matrices), you can disable this check by setting check_finished to False, but be careful not to disrupt ongoing runs.") - continue - # additional check that all the required types of matrices are present - if two_tracers is None: two_tracers = os.path.isfile(os.path.join(cov_dir, f"xi/xi_22.dat")) # simple heuristic if not provided explicitly, copied from post_process_auto - expected_matrix_names = ['c2_' + i + j for i in "12" for j in "12"] + ['c3_' + i + ',' + index2 for i in "12" for index2 in ("11", "12", "22")] + ['c4_' + indices for indices in ("11,11", "11,22", "12,11", "12,12", "12,21", "21,22", "22,22")] if two_tracers else ['c2_11', 'c3_1,11', 'c4_11,11'] # not including jackknife, because it is hard to know if we expect jackknife at this point - not_found_names = [matrix_name for matrix_name in expected_matrix_names if matrix_name not in output_group] - if not_found_names: - print_function(f"WARNING: {not_found_names} matrices not found in output group {output_group_name}. This may indicate that the C++ code is still running. Will skip this output group to avoid disrupting an ongoing run by collecting and/or deleting the text files. If you are sure that the run(s) in this directory finished (e.g., it/they timed out without producing the full matrices), you can disable this check by setting check_finished to False, but be careful not to disrupt ongoing runs.") - continue - - # check that the different matrices have the same subsamples by index - subsample_indices = [{suffix for suffix in matrix_filenames_dictionary.keys() if isinstance(suffix, int)} for matrix_filenames_dictionary in output_group.values()] - common_subsample_indices = subsample_indices[0].intersection(*subsample_indices) - - if any(this_subsample_indices != common_subsample_indices for this_subsample_indices in subsample_indices): - print_function(f"WARNING: Some matrices in output group {output_group_name} have different subsample indices. " + ("This may indicate that the C++ code is still running. Will skip this output group to avoid disrupting an ongoing run by collecting and/or deleting the text files. If you are sure that the run(s) in this directory finished (e.g., it/they timed out without producing the full matrices), you can disable this check by setting check_finished to False, but be careful not to disrupt ongoing runs." if check_finished else f"Will cut to the common subsample indices, leaving {len(common_subsample_indices)} subsamples.")) - if check_finished: continue # skip the output group to avoid disrupting an ongoing run - # otherwise, cut all to the common subsample indices - output_group = {matrix_name: {suffix: input_filename for suffix, input_filename in matrix_filenames_dictionary.items() if suffix in common_subsample_indices or isinstance(suffix, str)} for matrix_name, matrix_filenames_dictionary in output_group.items()} - - # now create and fill the dictionary to be saved in the numpy file - output_dictionary = {} - for matrix_name, matrix_filenames_dictionary in output_group.items(): - output_dictionary[matrix_name] = {} - for suffix, input_filename in matrix_filenames_dictionary.items(): - matrix = np.loadtxt(input_filename) - if matrix_name.startswith("c2") and matrix.ndim == 1: matrix = np.diag(matrix) # convert 1D c2 to a 2D diagonal matrix - if isinstance(suffix, str): # special treatment for string suffixes (at the moment, only "full") - output_dictionary[matrix_name + "_" + suffix] = matrix # this creates a separate array to be saved, which seems more convenient - elif isinstance(suffix, int): # integer suffixes are the subsample indices - output_dictionary[matrix_name][suffix] = matrix # save the subsample matrices in a sub-dictionary, which will be converted to a numpy array later - else: raise TypeError(f"Unexpected suffix type {type(suffix)} for {matrix_name} in {output_group_name}; expected int or str") - - # now all the remaining suffixes must be integers so can be sorted easily - output_dictionary[matrix_name] = np.array([output_dictionary[matrix_name][i_subsample] for i_subsample in sorted(output_dictionary[matrix_name].keys())]) - # this transformed the dictionary to numpy array, ordered by increasing subsample index - - # calculate the full as the mean of the subsamples - full_matrix_computed = np.mean(output_dictionary[matrix_name], axis = 0) - full_matrix_name = matrix_name + "_full" - if full_matrix_name in output_dictionary: - if not np.allclose(output_dictionary[full_matrix_name], full_matrix_computed): - print_function(f"WARNING: For {matrix_name} matrix, the loaded full is different from the average of subsamples. The latter will be saved.") - matrix_filenames_dictionary.pop("full") # remove the filename since it will be technically unused - output_dictionary[full_matrix_name] = full_matrix_computed # only override the full matrix if it is different from the computed one - - return_dictionary[output_group_name] = output_dictionary - - save_safe(cov_dir, output_group_name, output_dictionary, print_function) - - # now that the file is saved (not any earlier to be sure), can remove all the text files - # the list contains only the files that had their contents loaded and saved - if cleanup: - for matrix_filenames_dictionary in output_group.values(): - for input_filename in matrix_filenames_dictionary.values(): - os.remove(input_filename) - - print_function(f"Finished with output group {output_group_name}") + for threepcf, these_output_groups in enumerate([output_groups, output_groups_3pcf]): + for output_group_name, output_group in these_output_groups.items(): + print_function(f"Processing {2+threepcf}PCF output group {output_group_name}") + + # check if the run seems finished by looking for the presence of the full matrices, if check_finished is enabled + if check_finished: + unfinished_names = [matrix_name for matrix_name, matrix_filenames_dictionary in output_group.items() if "full" not in matrix_filenames_dictionary] + if unfinished_names: + print_function(f"WARNING: {unfinished_names} full matrices not found for {2+threepcf}PCF output group {output_group_name}. This may indicate that the C++ code is still running. Will skip this output group to avoid disrupting an ongoing run by collecting and/or deleting the text files. If you are sure that the run(s) in this directory finished (e.g., it/they timed out without producing the full matrices), you can disable this check by setting check_finished to False, but be careful not to disrupt ongoing runs.") + continue + # additional checks that all types of matrices are present + if threepcf: # different names for 3PCF, and only single-tracer for now + expected_matrix_names = [f'c{n}_{i}' for n in range(3, 6) for i in range(2)] + ['c6_1'] # exclude c6_0 term, because it should be small but is also hard to compute (see Section 5.2.3 and Appendix A of https://arxiv.org/abs/1910.04764) + else: # set 2PCF names, excluding jackknife + if two_tracers is None: two_tracers = os.path.isfile(os.path.join(cov_dir, f"xi/xi_22.dat")) # simple heuristic if not provided explicitly, copied from post_process_auto + expected_matrix_names = ['c2_' + i + j for i in "12" for j in "12"] + ['c3_' + i + ',' + index2 for i in "12" for index2 in ("11", "12", "22")] + ['c4_' + indices for indices in ("11,11", "11,22", "12,11", "12,12", "12,21", "21,22", "22,22")] if two_tracers else ['c2_11', 'c3_1,11', 'c4_11,11'] + not_found_names = [matrix_name for matrix_name in expected_matrix_names if matrix_name not in output_group] + if not_found_names: + print_function(f"WARNING: {not_found_names} matrices not found in {2+threepcf}PCF output group {output_group_name}. This may indicate that the C++ code is still running. Will skip this output group to avoid disrupting an ongoing run by collecting and/or deleting the text files. If you are sure that the run(s) in this directory finished (e.g., it/they timed out without producing the full matrices), you can disable this check by setting check_finished to False, but be careful not to disrupt ongoing runs.") + continue + + # check that the different matrices have the same subsamples by index + subsample_indices = [{suffix for suffix in matrix_filenames_dictionary.keys() if isinstance(suffix, int)} for matrix_filenames_dictionary in output_group.values()] + common_subsample_indices = subsample_indices[0].intersection(*subsample_indices) + + if any(this_subsample_indices != common_subsample_indices for this_subsample_indices in subsample_indices): + print_function(f"WARNING: Some matrices in output group {output_group_name} have different subsample indices. " + ("This may indicate that the C++ code is still running. Will skip this output group to avoid disrupting an ongoing run by collecting and/or deleting the text files. If you are sure that the run(s) in this directory finished (e.g., it/they timed out without producing the full matrices), you can disable this check by setting check_finished to False, but be careful not to disrupt ongoing runs." if check_finished else f"Will cut to the common subsample indices, leaving {len(common_subsample_indices)} subsamples.")) + if check_finished: continue # skip the output group to avoid disrupting an ongoing run + # otherwise, cut all to the common subsample indices + output_group = {matrix_name: {suffix: input_filename for suffix, input_filename in matrix_filenames_dictionary.items() if suffix in common_subsample_indices or isinstance(suffix, str)} for matrix_name, matrix_filenames_dictionary in output_group.items()} + + # now create and fill the dictionary to be saved in the numpy file + output_dictionary = {} + for matrix_name, matrix_filenames_dictionary in output_group.items(): + output_dictionary[matrix_name] = {} + for suffix, input_filename in matrix_filenames_dictionary.items(): + matrix = np.loadtxt(input_filename) + if matrix_name.startswith("c2") and matrix.ndim == 1: matrix = np.diag(matrix) # convert 1D c2 to a 2D diagonal matrix + if isinstance(suffix, str): # special treatment for string suffixes (at the moment, only "full") + output_dictionary[matrix_name + "_" + suffix] = matrix # this creates a separate array to be saved, which seems more convenient + elif isinstance(suffix, int): # integer suffixes are the subsample indices + output_dictionary[matrix_name][suffix] = matrix # save the subsample matrices in a sub-dictionary, which will be converted to a numpy array later + else: raise TypeError(f"Unexpected suffix type {type(suffix)} for {matrix_name} in {output_group_name}; expected int or str") + + # now all the remaining suffixes must be integers so can be sorted easily + output_dictionary[matrix_name] = np.array([output_dictionary[matrix_name][i_subsample] for i_subsample in sorted(output_dictionary[matrix_name].keys())]) + # this transformed the dictionary to numpy array, ordered by increasing subsample index + + # calculate the full as the mean of the subsamples + full_matrix_computed = np.mean(output_dictionary[matrix_name], axis = 0) + full_matrix_name = matrix_name + "_full" + if full_matrix_name in output_dictionary: + if not np.allclose(output_dictionary[full_matrix_name], full_matrix_computed): + print_function(f"For {matrix_name} matrix from {2+threepcf}PCF ouput group {output_group_name}, the loaded full is different from the average of subsamples. The latter will be saved.") + matrix_filenames_dictionary.pop("full") # remove the filename since it will be technically unused + output_dictionary[full_matrix_name] = full_matrix_computed + + return_dictionaries[threepcf][output_group_name] = output_dictionary + + save_safe(cov_dir, output_group_name, output_dictionary, threepcf=threepcf, print_function=print_function) + + # now that the file is saved (not any earlier to be sure), can remove all the text files + # the list contains only the files that had their contents loaded and saved + if cleanup: + for matrix_filenames_dictionary in output_group.values(): + for input_filename in matrix_filenames_dictionary.values(): + os.remove(input_filename) + + print_function(f"Finished with {2+threepcf}PCF output group {output_group_name}") # remove subdirectories too if they are empty - if cleanup and len(output_groups) > 0: # also don't want to delete directories if no output groups were found. This can happen when the run is ongoing at an early stage and the text files are not yet produced, removing a directory will cause an error when the C++ code eventually tries to write them - rmdir_if_exists_and_empty(cov_dir_all) - rmdir_if_exists_and_empty(cov_dir_jack) + if cleanup: + if len(output_groups) > 0: # also don't want to delete 2PCF directories if no 2PCF output groups were found. This can happen when the run is ongoing at an early stage and the text files are not yet produced, removing a directory will cause an error when the C++ code eventually tries to write them + rmdir_if_exists_and_empty(cov_dir_all) + rmdir_if_exists_and_empty(cov_dir_jack) + if len(output_groups_3pcf) > 0: # similarly, don't want to delete the 3PCF directory if no 3PCF output groups were found + rmdir_if_exists_and_empty(cov_dir_3pcf_all) - return return_dictionary + return return_dictionaries -def load_raw_covariances(file_root: str, label: str, n_samples: None | int | Iterable[int] | Iterable[bool] = None, check_finished: bool = True, two_tracers: bool | None = None, print_function: Callable[[str], None] = print) -> dict[str]: +def load_raw_covariances(file_root: str, label: str, threepcf: bool = False, n_samples: None | int | Iterable[int] | Iterable[bool] = None, check_finished: bool = True, two_tracers: bool | None = None, print_function: Callable[[str], None] = print) -> dict[str]: """ Load the raw covariance matrices as a dictionary. Uses the Numpy file if it exists, otherwise tried to run the collection function. @@ -167,6 +188,8 @@ def load_raw_covariances(file_root: str, label: str, n_samples: None | int | Ite label specifies the number of radial bins and either the number of angular (mu) bins (nN_mM) or the maximum (even) multipole index (nN_lL). + threepcf chooses between 2PCF (False, default) and 3PCF (True). + n_samples allows to select subsamples flexibly: - if None (default) returns all samples; - if a positive integer, returns as many samples from the beginning; @@ -176,13 +199,13 @@ def load_raw_covariances(file_root: str, label: str, n_samples: None | int | Ite With check_finished enabled (default), before collecting raw covariance matrices (if needed), performs a heuristic check whether the run seems finished by looking for the presence of the full covariance matrices. If they are not found, the collection does not proceed and a warning is issued. This is to prevent disrupting ongoing runs by collecting and/or deleting the text files. If you want to check convergence of timed-out run(s) that did not produce the full matrices, you can disable this check, but be careful not to disrupt ongoing runs. With check_finished enabled, the code additionally checks presence of all types of matrices for two tracers. Pass two_tracers = True to enable this check explicitly, or two_tracers = False to disable it. By default (two_tracers = None), the code will try to guess whether the run is for two tracers by looking for the presence of the xi_22.dat file, which is only produced for two tracers by :func:`RascalC.run_cov`. """ - input_filename = os.path.join(file_root, f"Raw_Covariance_Matrices_{label}.npz") + input_filename = os.path.join(file_root, "Raw_Covariance_Matrices" + "_3PCF" * threepcf + f"_{label}.npz") if os.path.isfile(input_filename): raw_cov = dict(np.load(input_filename)) else: print_function(f"Collecting the raw covariance matrices from {file_root}") - result = collect_raw_covariance_matrices(file_root, check_finished=check_finished, two_tracers=two_tracers, print_function=print_function) + result = collect_raw_covariance_matrices(file_root, check_finished=check_finished, two_tracers=two_tracers, print_function=print_function)[threepcf] if label not in result: - raise ValueError(f"Raw covariance matrices for {label} not produced. Check the n and m/max_l values.") + raise ValueError(f"Raw {2+threepcf}PCF covariance matrices for {label} not produced. Check the n and m/max_l values.") raw_cov = result[label] if n_samples is None: return raw_cov # return the full set elif isinstance(n_samples, int): @@ -218,7 +241,7 @@ def load_raw_covariances_smu(file_root: str, n: int, m: int, n_samples: None | i With check_finished enabled, the code additionally checks presence of all types of matrices for two tracers. Pass two_tracers = True to enable this check explicitly, or two_tracers = False to disable it. By default (two_tracers = None), the code will try to guess whether the run is for two tracers by looking for the presence of the xi_22.dat file, which is only produced for two tracers by :func:`RascalC.run_cov`. """ label = f"n{n}_m{m}" - return load_raw_covariances(file_root, label, n_samples, check_finished, two_tracers, print_function) + return load_raw_covariances(file_root, label, threepcf=False, n_samples=n_samples, check_finished=check_finished, two_tracers=two_tracers, print_function=print_function) def load_raw_covariances_legendre(file_root: str, n: int, max_l: int, n_samples: None | int | Iterable[int] | Iterable[bool] = None, check_finished: bool = True, two_tracers: bool | None = None, print_function = print) -> dict[str]: @@ -238,10 +261,29 @@ def load_raw_covariances_legendre(file_root: str, n: int, max_l: int, n_samples: With check_finished enabled, the code additionally checks presence of all types of matrices for two tracers. Pass two_tracers = True to enable this check explicitly, or two_tracers = False to disable it. By default (two_tracers = None), the code will try to guess whether the run is for two tracers by looking for the presence of the xi_22.dat file, which is only produced for two tracers by :func:`RascalC.run_cov`. """ label = f"n{n}_l{max_l}" - return load_raw_covariances(file_root, label, n_samples, check_finished, two_tracers, print_function) + return load_raw_covariances(file_root, label, threepcf=False, n_samples=n_samples, check_finished=check_finished, two_tracers=two_tracers, print_function=print_function) + + +def load_raw_covariances_3pcf_legendre(file_root: str, n: int, max_l: int, n_samples: None | int | Iterable[int] | Iterable[bool] = None, check_finished: bool = True, print_function = print) -> dict[str]: + """ + Load the raw covariance matrices from the 3PCF Legendre mode as a dictionary. Uses the Numpy file if it exists, otherwise tried to run the collection function. + + file_root is the directory to look in. + n and m are numbers of radial and angular (mu) bins respectively. + + n_samples allows to select subsamples flexibly: + - if None (default) returns all samples; + - if a positive integer, returns as many samples from the beginning; + - if a sequence of integers, returns subsamples with indices from this sequence; + - sequence of boolean values is interpreted as a boolean mask for subsamples. + + With check_finished enabled (default), before collecting raw covariance matrices (if needed), performs a heuristic check whether the run seems finished by looking for the presence of the full covariance matrices. If they are not found, the collection does not proceed and a warning is issued. This is to prevent disrupting ongoing runs by collecting and/or deleting the text files. If you want to check convergence of timed-out run(s) that did not produce the full matrices, you can disable this check, but be careful not to disrupt ongoing runs. + """ + label = f"n{n}_l{max_l}" + return load_raw_covariances(file_root, label, threepcf=True, n_samples=n_samples, check_finished=check_finished, print_function=print_function) -def cat_raw_covariance_matrices(n: int, mstr: str, input_roots: list[str], ns_samples: list[None | int | Iterable[int] | Iterable[bool]], output_root: str, collapse_factor: int = 1, check_finished: bool = True, two_tracers: bool | None = None, print_function: Callable[[str], None] = print) -> dict[str]: +def cat_raw_covariance_matrices(n: int, mstr: str, input_roots: list[str], ns_samples: list[None | int | Iterable[int] | Iterable[bool]], output_root: str, collapse_factor: int = 1, threepcf: bool = False, check_finished: bool = True, two_tracers: bool | None = None, print_function: Callable[[str], None] = print) -> dict[str]: """ Catenate the raw covariance matrices from one or more input directories, assuming they are from similar runs: the same number of input random points, ``N2``, ``N3``, ``N4`` and ``loops_per_sample``. @@ -270,6 +312,9 @@ def cat_raw_covariance_matrices(n: int, mstr: str, input_roots: list[str], ns_sa collapse_factor : positive integer (Optional) allows to reduce the number of subsamples by averaging over a given number of adjacent subsamples. Default value is 1, meaning no reduction. + + threepcf : boolean + (Optional) whether the runs are for 3PCF. Default is False, meaning 2PCF. check_finished : boolean (Optional) whether to check if the runs are finished before collecting raw covariance matrices (if that is needed). Default is True. @@ -277,7 +322,7 @@ def cat_raw_covariance_matrices(n: int, mstr: str, input_roots: list[str], ns_sa two_tracers : boolean | None (Optional) whether the runs are for two tracers (for an additional check for finished runs). If None (default), the code will try to guess by looking for the presence of the xi_22.dat file. - print_function : Callable + print_function : Callable[[str], None] (Optional) custom function to use for printing. Default is ``print``. Returns @@ -292,7 +337,7 @@ def cat_raw_covariance_matrices(n: int, mstr: str, input_roots: list[str], ns_sa label = f"n{n}_{mstr}" result = {} for index, (input_root, n_samples) in enumerate(zip(input_roots, ns_samples)): - input_file = load_raw_covariances(input_root, label, n_samples, check_finished, two_tracers, print_function) + input_file = load_raw_covariances(input_root, label, threepcf=threepcf, n_samples=n_samples, check_finished=check_finished, two_tracers=two_tracers, print_function=print_function) # ignore full arrays input_file = {key: value for (key, value) in input_file.items() if not key.endswith("_full")} # check that the keys are the same, unless the result is brand new diff --git a/docs/index.rst b/docs/index.rst index 462a27bf..60cee154 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -21,12 +21,13 @@ For general usage, we provide :doc:`library/fundamentals` and :doc:`library/tuto .. toctree:: :maxdepth: 2 - :caption: Python library (2PCF) + :caption: Python library library/installation library/fundamentals library/tutorials-examples library/advanced-functions + library/3pcf .. toctree:: :maxdepth: 2 diff --git a/docs/legacy/main-code.rst b/docs/legacy/main-code.rst index ba062d3b..adf13ca5 100644 --- a/docs/legacy/main-code.rst +++ b/docs/legacy/main-code.rst @@ -128,8 +128,8 @@ Input parameters for the RascalC code may be specified by passing options on the - ``-boxsize`` (*boxsize*): Sets the size of the computational domain if the periodic box is enabled, or the random particles are created in RascalC. (Default: 2000) - ``-effective_norm``: Flag indicating to use the effective number of random particles (:math:`\sum(w_i)^2/\sum(w_i^2)`, taking into account weights but invariant to overall weight rescaling) instead of simple counting (default). Should match the way the -norm value was computed for the data galaxies externally. - ``-seed`` (*seed*): Random number generator seed. If given, allows to reproduce the results with the same settings, except the number of threads. -- ``-start_integral_index`` (*start_integral_index*): If given, chooses from which of 7 integrals (numbered 1 through 7) to start in the multi-tracer regime. Default 1 (to cover all integrals). Intended to help complete the timed-out multi-tracer runs by skipping the integrals computed in an unfinished run. -- ``-last_integral_index`` (*last_integral_index*): If given, chooses at which of 7 integrals (numbered 1 through 7) to stop in the multi-tracer regime. Default 7 (to cover all integrals). Can be used to fit into the time limit if the full 7 integrals are expected to take longer. +- ``-start_integral_index`` (*start_integral_index*): If given, chooses from which of 7 integrals (numbered 1 through 7) to start in the 2PCF multi-tracer regime, or 2 integrals for 3PCF. Default 1 (to cover all integrals). Intended to help complete the timed-out runs by skipping the integrals computed in an unfinished run. +- ``-last_integral_index`` (*last_integral_index*): If given, chooses at which of 7 integrals (numbered 1 through 7) to stop in the 2PCF multi-tracer regime, or 2 integrals for 3PCF. Default 7 for 2PCF multi-tracer and 2 for 3PCF (to cover all integrals). Can be used to fit into the time limit if the full set of integrals is expected to take longer. - ``-mumin`` (*mumin*): Minimum :math:`\mu` binning to use in the analysis. (Default: 0, or -1 in 3PCF mode) - ``-mumax`` (*mumax*): Maximum :math:`\mu` binning to use in the analysis. (Default: 1) - ``-cf_loops`` (*cf_loops*): Number of iterations over which to refine the correlation functions. diff --git a/docs/library/3pcf.rst b/docs/library/3pcf.rst new file mode 100644 index 00000000..65de5c2e --- /dev/null +++ b/docs/library/3pcf.rst @@ -0,0 +1,31 @@ +3PCF covariances (under development) +==================================== + +The main computation wrapper +---------------------------- + +This Python wrapper for ``RascalC`` 3PCF is interfaced with ``ENCORE`` for 3PCF ``pycorr`` for 2PCF. +Many of the arguments are intentionally similar to ``pycorr.TwoPointCorrelationFunction`` `high-level interface `_. + +Please bear with the long description; you can pay less attention to settings labeled optional in the beginning. + +.. autofunction:: RascalC.run_cov_3pcf + + +Post-processing +--------------- + +The post-processing routine is invoked automatically at the end of the main wrapper function (:func:`RascalC.run_cov_3pcf`), so in many circumstances you may not need to run it separately. +However, knowing the post-processing routine is useful for timed-out runs, switching the mode, testing different cuts and/or output combinations in cases of insufficient convergence, etc. + +After post-processing, you probably want to run the extra convergence check (see :mod:`RascalC.convergence_check_extra`). + +.. autofunction:: RascalC.post_process_3pcf.post_process_3pcf + + +Loading and exporting the final covariance matrices +--------------------------------------------------- + +.. automodule:: RascalC.cov_utils_3pcf + :members: + :member-order: bysource \ No newline at end of file diff --git a/docs/library/fundamentals.rst b/docs/library/fundamentals.rst index 9860827d..579615b0 100644 --- a/docs/library/fundamentals.rst +++ b/docs/library/fundamentals.rst @@ -1,5 +1,5 @@ -Fundamentals -============ +Fundamentals (2PCF covariances) +=============================== Introduction ------------ diff --git a/grid_covariance.cpp b/grid_covariance.cpp index c40de0df..ae3ed3f3 100644 --- a/grid_covariance.cpp +++ b/grid_covariance.cpp @@ -281,32 +281,33 @@ int main(int argc, char *argv[]) { rescale.refine_wrapper(&par, all_grid, all_cf, all_rd, no_functions); #ifdef THREE_PCF - // Compute threePCF integrals - compute_integral(&all_grid[0],&par,&all_cf[0],&all_rd[0],&all_survey[0],1); // final digit is iteration number - compute_integral(&all_grid[0],&par,&all_cf[0],&all_rd[0],&all_survey[0],2); + // Compute 3PCF integrals + compute_integral(&all_grid[0], &par, &all_cf[0], &all_rd[0], &all_survey[0], 0); // final digit is the iteration number. For 3PCF, 0 and nonzero (1) encode different sampling strategies as described in the paper (https://arxiv.org/abs/1910.04764, Section 5.2.2) + compute_integral(&all_grid[0], &par, &all_cf[0], &all_rd[0], &all_survey[0], 1); + // multi-tracer for 3PCF is not implemented #elif (defined LEGENDRE || defined POWER) // Compute integrals - compute_integral(all_grid,&par,all_cf,all_rd,all_survey,1,1,1,1,1); // final digit is iteration number - - if(par.multi_tracers==true){ - compute_integral(all_grid,&par,all_cf,all_rd,all_survey,1,2,1,1,2); - compute_integral(all_grid,&par,all_cf,all_rd,all_survey,1,2,2,1,3); - compute_integral(all_grid,&par,all_cf,all_rd,all_survey,1,2,1,2,4); - compute_integral(all_grid,&par,all_cf,all_rd,all_survey,1,1,2,2,5); - compute_integral(all_grid,&par,all_cf,all_rd,all_survey,2,1,2,2,6); - compute_integral(all_grid,&par,all_cf,all_rd,all_survey,2,2,2,2,7); + compute_integral(all_grid, &par, all_cf, all_rd, all_survey, 1, 1, 1, 1, 1); // first 4 digits encode the choice of 4 tracers involved, final digit is the iteration number (mostly used for informational printouts). Need 1 combination/iteration for 1 tracer, and 7 iterations for 2 tracers, can be extended for more + + if (par.multi_tracers == true) { + compute_integral(all_grid, &par, all_cf, all_rd, all_survey, 1, 2, 1, 1, 2); + compute_integral(all_grid, &par, all_cf, all_rd, all_survey, 1, 2, 2, 1, 3); + compute_integral(all_grid, &par, all_cf, all_rd, all_survey, 1, 2, 1, 2, 4); // tracers 1, 2, 1, 2 should in principle be identical to 1, 2, 2, 1, but the implementation only considers first-third and second-fourth pairings for the 4-point term C4 (not including the first-fourth and second-third pairings) which makes them different. In other cases (for 2 tracers) this does not matter + compute_integral(all_grid, &par, all_cf, all_rd, all_survey, 1, 1, 2, 2, 5); + compute_integral(all_grid, &par, all_cf, all_rd, all_survey, 2, 1, 2, 2, 6); + compute_integral(all_grid, &par, all_cf, all_rd, all_survey, 2, 2, 2, 2, 7); } #else // Compute integrals - compute_integral(all_grid,&par,all_weights,all_cf,all_rd,1,1,1,1,1); // final digit is iteration number - - if(par.multi_tracers==true){ - compute_integral(all_grid,&par,all_weights,all_cf,all_rd,1,2,1,1,2); - compute_integral(all_grid,&par,all_weights,all_cf,all_rd,1,2,2,1,3); - compute_integral(all_grid,&par,all_weights,all_cf,all_rd,1,2,1,2,4); - compute_integral(all_grid,&par,all_weights,all_cf,all_rd,1,1,2,2,5); - compute_integral(all_grid,&par,all_weights,all_cf,all_rd,2,1,2,2,6); - compute_integral(all_grid,&par,all_weights,all_cf,all_rd,2,2,2,2,7); + compute_integral(all_grid, &par, all_weights, all_cf, all_rd, 1, 1, 1, 1, 1); // first 4 digits encode the choice of 4 tracers involved, final digit is the iteration number (mostly used for informational printouts). Need 1 combination/iteration for 1 tracer, and 7 iterations for 2 tracers, can be extended for more + + if (par.multi_tracers == true) { + compute_integral(all_grid, &par, all_weights, all_cf, all_rd, 1, 2, 1, 1, 2); + compute_integral(all_grid, &par, all_weights, all_cf, all_rd, 1, 2, 2, 1, 3); + compute_integral(all_grid, &par, all_weights, all_cf, all_rd, 1, 2, 1, 2, 4); // tracers 1, 2, 1, 2 should in principle be identical to 1, 2, 2, 1, but the implementation only considers first-third and second-fourth pairings for the 4-point term C4 (not including the first-fourth and second-third pairings) which makes them different. In other cases (for 2 tracers) this does not matter + compute_integral(all_grid, &par, all_weights, all_cf, all_rd, 1, 1, 2, 2, 5); + compute_integral(all_grid, &par, all_weights, all_cf, all_rd, 2, 1, 2, 2, 6); + compute_integral(all_grid, &par, all_weights, all_cf, all_rd, 2, 2, 2, 2, 7); } #endif diff --git a/modules/compute_integral_3pcf.h b/modules/compute_integral_3pcf.h index a82c3db9..1fe52280 100644 --- a/modules/compute_integral_3pcf.h +++ b/modules/compute_integral_3pcf.h @@ -52,14 +52,19 @@ class compute_integral{ // MAIN FUNCTION TO COMPUTE INTEGRALS int tot_iter=2; // total number of integral sets to compute + assert((iter_no >= 0) && (iter_no < tot_iter)); // ensure that the iteration number/index is as it should be + + if ((iter_no+1 < par->start_integral_index) || (iter_no+1 > par->last_integral_index)) { // here, iter_no starts from 0 for compatibility with other pieces of 3PCF code, whereas integral_indices start from 1 for compatibility with 2PCF + printf("# Skipping integral computation %d of %d; only running %d through %d (both inclusive).\n", iter_no+1, tot_iter, par->start_integral_index, par->last_integral_index); + return; + } nbin = par->nbin; // number of radial bins mbin = par->mbin; // number of Legendre bins + STimer* LoopTimes = new STimer[par->max_loops]; STimer initial, TotalTime; // Time initialization - initial.Start(); - - int convergence_counter=0, printtime=0;// counter to stop loop early if convergence is reached. + initial.Start(); //-----------INITIALIZE OPENMP + CLASSES---------- unsigned long seed_step = (unsigned long)(std::numeric_limits::max()) / (unsigned long)(par->max_loops); // restrict the seeds to 32 bits to avoid possible collisions, as most GSL random generators only accept 32-bit seeds and truncate the higher bits. unsigned long is at least 32 bits but can be more. @@ -67,11 +72,17 @@ class compute_integral{ if (par->random_seed) { std::random_device urandom("/dev/urandom"); std::uniform_int_distribution dist(0, seed_step-1); // distribution of integers with these limits, inclusive - seed_shift = dist(urandom); // for each thread, it will be seed = seed_step * n_loops + seed_shift, where n_loops goes from 0 to max_loops-1 => no overflow and guaranteed unique for each thread + seed_shift = dist(urandom); // for each iteration, it will be seed = seed_step * n_loops + seed_shift, where n_loops goes from 0 to max_loops-1 => no overflow and guaranteed unique for each thread + par->seed = seed_shift; // save the random seed + par->random_seed = false; // avoid random seed resetting in the future — different seeds for xi rescaling and the main computation are much harder to reproduce } + printf("# INFO: the base RNG seed is %lu, incremented by %lu in each integration loop (matching these should guarantee the same results in a completed computation with the same input catalogs and correlation functions).\n", seed_shift, seed_step); gsl_rng_env_setup(); // initialize gsl rng + int completed_loops = 0; // counter of completed loops, since may not finish in index order + uint64 used_triples_per_sample = 0, used_quads_per_sample = 0, used_quints_per_sample = 0, used_hexes_per_sample = 0; // triple, quad, quint and hex counts for the normalization of the current output sample Integrals sumint(par, cf, survey_corr); // total integral + Integrals outint(par, cf, survey_corr); // current output integral uint64 tot_triples=0, tot_quads=0, tot_quints=0, tot_hexes=0; // global number of particle sets used (including those rejected for being in the wrong bins) uint64 cell_attempt3=0,cell_attempt4=0,cell_attempt5=0,cell_attempt6=0; // number of k,l,m,n cells attempted @@ -87,15 +98,15 @@ class compute_integral{ TotalTime.Start(); // Start timer #ifdef OPENMP - #pragma omp parallel firstprivate(seed_step,seed_shift,par,printtime,grid,cf) shared(sumint,TotalTime,gsl_rng_default,rd) reduction(+:convergence_counter,cell_attempt3,cell_attempt4,cell_attempt5,cell_attempt6,used_cell3,used_cell4,used_cell5,used_cell6,tot_triples,tot_quads,tot_quints,tot_hexes) + #pragma omp parallel firstprivate(seed_step,seed_shift,par,grid,cf) shared(sumint,TotalTime,LoopTimes,gsl_rng_default,rd) reduction(+:cell_attempt3,cell_attempt4,cell_attempt5,cell_attempt6,used_cell3,used_cell4,used_cell5,used_cell6,tot_triples,tot_quads,tot_quints,tot_hexes) { // start parallel loop // Decide which thread we are in int thread = omp_get_thread_num(); assert(omp_get_num_threads()<=par->nthread); - if (thread==0) printf("# Starting integral computation %d of %d on %d threads.\n", iter_no, tot_iter, omp_get_num_threads()); + if (thread==0) printf("# Starting integral computation %d of %d on %d threads.\n", iter_no+1, tot_iter, omp_get_num_threads()); #else int thread = 0; - printf("# Starting integral computation %d of %d single threaded.\n",iter_no,tot_iter); + printf("# Starting integral computation %d of %d single threaded.\n", iter_no + 1, tot_iter); { // start loop #endif @@ -110,7 +121,9 @@ class compute_integral{ Float *correction_ijk, *legendre_ijk, *xi_pass, *xi_pass2, xi_pass3=0, *w_ijk, *w_ijkl, *w_ijklm; Float norm_jk=0, norm_kl=0, norm_lm=0; // arrays to store xi and weight values int *bins_ijk, bin_jk=0, bin_kl=0, bin_lm=0; // i-j separation bin +#ifdef PRINTPERCENTS Float percent_counter; +#endif int x, prim_id_1D; integer3 delta2, delta3, delta4, delta5, delta6, prim_id, sec_id, thi_id, fou_id, fif_id, six_id; Float3 cell_sep2,cell_sep3,cell_sep4,cell_sep5,cell_sep6; @@ -140,15 +153,11 @@ class compute_integral{ #pragma omp for schedule(dynamic) #endif for (int n_loops = 0; n_loopsmax_loops; n_loops++){ +#ifdef PRINTPERCENTS percent_counter=0.; - loc_used_triples=0; loc_used_quads=0; loc_used_quints=0; loc_used_hexes=0; - - // End loops early if convergence has been acheived - if (convergence_counter==10){ - if (printtime==0) printf("1 percent convergence acheived in C6 10 times, exiting.\n"); - printtime++; - continue; - } +#endif + loc_used_triples=0; loc_used_quads=0; loc_used_quints=0; loc_used_hexes=0; + LoopTimes[n_loops].Start(); // Set/reset the RNG seed based on loop number instead of thread number to reproduce results with different number of threads but other parameters kept the same. Individual subsamples may differ because they are accumulated/written in order of loop completion which may depend on external factors at runtime, but the final integrals should be the same. gsl_rng_set(locrng, seed_step * (unsigned long)n_loops + seed_shift); // the second number, seed, will not overflow and will be unique for each loop in one run. Here, seed_shift is a random number between 0 and seed_step-1, inclusively. @@ -157,11 +166,13 @@ class compute_integral{ // LOOP OVER ALL FILLED I CELLS for (int n1=0; n1nf;n1++){ +#ifdef PRINTPERCENTS // Print time left if((float(n1)/float(grid->nf)*100)>=percent_counter){ - printf("Integral %d of %d, run %d of %d on thread %d: Using cell %d of %d - %.0f percent complete\n",iter_no+1,tot_iter,1+n_loops/par->nthread, int(ceil(float(par->max_loops)/(float)par->nthread)),thread, n1+1,grid->nf,percent_counter); + printf("Integral %d of %d, iteration %d of %d on thread %d: Using cell %d of %d - %.0f percent complete\n", iter_no+1, tot_iter, n_loops, par->max_loops, thread, n1+1, grid->nf, percent_counter); percent_counter+=5.; } +#endif // Pick first particle prim_id_1D = grid->filled[n1]; // 1d ID for cell i @@ -241,7 +252,7 @@ class compute_integral{ // LOOP OVER N5 M CELLS for (int n5=0; n5N5; n5++){ - cell_attempt5+=1; // new fourth cell attempted + cell_attempt5+=1; // new fifth cell attempted // Draw fifth cell from l or i cell delta5 = rd->random_xidraw(locrng,&p5); @@ -258,16 +269,18 @@ class compute_integral{ if(x==1) continue; if((pid_m==pid_l)||(pid_m==pid_k)||(pid_m==pid_j)) continue; - used_cell5+=1; // new fourth cell used + used_cell5+=1; // new fifth cell used p5*=p4/(double)filn; - // Now compute the four-point integral - locint.fifth(prim_list, prim_ids, pln, particle_j, particle_k, particle_l, particle_m, pid_j, pid_k, pid_l, pid_m, p5, w_ijkl, bins_ijk, correction_ijk, legendre_ijk, xi_pass, xi_pass2, norm_kl, bin_kl, w_ijklm, xi_pass3, norm_lm, bin_lm, iter_no); + // Now compute the five-point integral + locint.fifth(prim_list, prim_ids, pln, particle_j, particle_k, particle_l, particle_m, pid_j, pid_k, pid_l, pid_m, p5, w_ijkl, bins_ijk, correction_ijk, legendre_ijk, xi_pass, xi_pass2, norm_kl, bin_kl, w_ijklm, xi_pass3, norm_lm, bin_lm, iter_no); + + if (iter_no == 0) continue; // skip the first 6-point term, because it should be small but is also hard to compute (see Section 5.2.3 and Appendix A of https://arxiv.org/abs/1910.04764) // LOOP OVER N6 N CELLS for (int n6=0; n6N6; n6++){ - cell_attempt6+=1; // new fourth cell attempted + cell_attempt6+=1; // new sixth cell attempted // Draw sixth cell from m cell weighted by xi(r) delta6 = rd->random_xidraw(locrng,&p6); @@ -288,7 +301,7 @@ class compute_integral{ p6*=p5/(double)siln; - // Now compute the four-point integral + // Now compute the six-point integral locint.sixth(prim_list, prim_ids, pln, particle_j, particle_k, particle_l, particle_m, particle_n, pid_j, pid_k, pid_l, pid_m, pid_n, p6, w_ijklm, bins_ijk, correction_ijk, legendre_ijk, xi_pass, xi_pass2, xi_pass3, norm_lm, bin_lm, iter_no); } } @@ -307,35 +320,49 @@ class compute_integral{ #pragma omp critical // only one processor can access at once #endif { - if ((n_loops+1)%par->nthread==0){ // Print every nthread loops + int subsample_index = completed_loops / par->loops_per_sample; // index of output subsample for this loop + completed_loops++; // increment completed loops counter, since they may be done not according to n_loops order + printf("Integral %d of %d, iteration %d of %d on thread %d completed (%d/%d)\n", iter_no+1, tot_iter, n_loops, par->max_loops, thread, completed_loops, par->max_loops); + if (completed_loops % par->nthread == 0) { // Print every nthread completed loops TotalTime.Stop(); // interrupt timing to access .Elapsed() int current_runtime = TotalTime.Elapsed(); - int remaining_time = current_runtime/((n_loops+1)/par->nthread)*(par->max_loops/par->nthread-(n_loops+1)/par->nthread); // estimated remaining time - fprintf(stderr,"\nFinished integral loop %d of %d after %d s. Estimated time left: %2.2d:%2.2d:%2.2d hms, i.e. %d s.\n",n_loops+1,par->max_loops, current_runtime,remaining_time/3600,remaining_time/60%60, remaining_time%60,remaining_time); - + int remaining_time = current_runtime*(par->max_loops - completed_loops)/completed_loops; // estimated remaining time + fprintf(stderr, "\nFinished %d integral loops of %d after %d s. Estimated time left: %2.2d:%2.2d:%2.2d hms, i.e. %d s.\n", completed_loops, par->max_loops, current_runtime, remaining_time/3600, remaining_time/60%60, remaining_time%60, remaining_time); + TotalTime.Start(); // Restart the timer + } + int accumulated_loops = subsample_index * par->loops_per_sample; // number of (completed) integral loops accumulated in sumint (as can be seen later, it is updated every loops_per_sample completed loops) + if ((subsample_index > 0) && (completed_loops == accumulated_loops + 1)) { // the condition when the Frobenius difference after adding one loop is most straightforward to compute sensibly, because sumint is updated only every loops_per_sample completed loops. Also works for loops_per_sample=1 unlike the possible alternative condition, completed_loops % loops_per_sample == 1 Float frob_C3, frob_C4, frob_C5, frob_C6; - sumint.frobenius_difference_sum(&locint,n_loops, frob_C3, frob_C4, frob_C5, frob_C6); - if(frob_C6<0.01) convergence_counter++; - if (n_loops!=0){ - fprintf(stderr,"Frobenius percent difference after loop %d is %.3f (C3), %.3f (C4), %.3f (C5), %.3f (C6) \n",n_loops,frob_C3, frob_C4, frob_C5, frob_C6); - } + sumint.frobenius_difference_sum(&locint, accumulated_loops, frob_C3, frob_C4, frob_C5, frob_C6); // computes the Frobenius relative differences (in percents) after adding one integral loop result (locint) to the accumulation of accumulated_loops loops in sumint; the method signature is different with and without jackknives + fprintf(stderr, "Frobenius percent difference between %d and %d loops is %.3f (C3), %.3f (C4), %.3f (C5), %.3f (C6) \n", accumulated_loops, completed_loops, frob_C3, frob_C4, frob_C5, frob_C6); } - + // Sum up integrals - sumint.sum_ints(&locint); - - // Save output after each loop - std::string output_string = string_format("%d", n_loops); - - locint.normalize(grid->norm,(Float)loc_used_triples, (Float)loc_used_quads, (Float)loc_used_quints, (Float)loc_used_hexes); - locint.save_integrals(output_string.c_str(), 0, iter_no); - locint.sum_total_counts(cnt3, cnt4, cnt5, cnt6); + outint.sum_ints(&locint); // subsample total + + // Sum up pairs, triples, quads for the group + used_triples_per_sample += loc_used_triples; + used_quads_per_sample += loc_used_quads; + used_quints_per_sample += loc_used_quints; + used_hexes_per_sample += loc_used_hexes; + + // Save output if the group is done + if (completed_loops % par->loops_per_sample == 0) { + sumint.sum_ints(&outint); // add to grand total + std::string output_string = string_format("%d", subsample_index); + outint.normalize(grid->norm, (Float)used_triples_per_sample, (Float)used_quads_per_sample, (Float)used_quints_per_sample, (Float)used_hexes_per_sample); + outint.save_integrals(output_string.c_str(), 0, iter_no); + // Reset the current output sample variables + outint.reset(); + used_triples_per_sample = used_quads_per_sample = used_quints_per_sample = used_hexes_per_sample = 0; + } + locint.sum_total_counts(cnt3, cnt4, cnt5, cnt6); locint.reset(); } - + LoopTimes[n_loops].Stop(); } // end cycle loop // Free up allocated memory at end of process @@ -365,7 +392,11 @@ class compute_integral{ cnt6/=(9.*mbin*mbin); int runtime = TotalTime.Elapsed(); - printf("\n\nINTEGRAL %d OF %d COMPLETE\n",iter_no+1,tot_iter); + printf("\n\nINTEGRAL %d OF %d COMPLETE\n",iter_no+1,tot_iter); + for (int n_loop = 0; n_loop < par->max_loops; n_loop++) { + int loop_runtime = LoopTimes[n_loop].Elapsed(); + fprintf(stderr, "Loop %d time: %d s, i.e. %2.2d:%2.2d:%2.2d hms\n", n_loop, loop_runtime, loop_runtime/3600, loop_runtime/60%60, loop_runtime%60); + } fprintf(stderr, "\nTotal process time for %.2e sets of cells and %.2e hexes of particles: %d s, i.e. %2.2d:%2.2d:%2.2d hms\n", double(used_cell6),double(tot_hexes),runtime, runtime/3600,runtime/60%60,runtime%60); printf("We tried %.2e triples, %.2e quads, %.2e quints and %.2e hexes of cells.\n",double(cell_attempt3),double(cell_attempt4),double(cell_attempt5),double(cell_attempt6)); printf("Of these, we accepted %.2e triples, %.2e quads, %.2e quints and %.2e hexes of cells.\n",double(used_cell3),double(used_cell4),double(used_cell5),double(used_cell6)); diff --git a/modules/correlation_function.h b/modules/correlation_function.h index 5b0b5728..42d46a40 100644 --- a/modules/correlation_function.h +++ b/modules/correlation_function.h @@ -116,7 +116,8 @@ class CorrelationFunction{ void readData(const char *filename,double **x,double **y, double **z,int *np,int *mp){ - char line[100000]; + const size_t line_size = 1000000; + char line[line_size]; char * pch; int n=0,m=0; double x0; @@ -127,7 +128,7 @@ class CorrelationFunction{ } //Count lines and columns - while (fgets(line,100000,fp)!=NULL) { + while (fgets(line, line_size, fp)!=NULL) { if (line[0]=='#') continue; if (line[0]=='\n') continue; if(n==0){ @@ -176,7 +177,7 @@ class CorrelationFunction{ } //Read content of lines and columns - while (fgets(line,100000,fp)!=NULL) { + while (fgets(line, line_size, fp)!=NULL) { if (line[0]=='#') continue; if (line[0]=='\n') continue; if(lnnr==0) readLine(line, x, &nx); diff --git a/modules/integrals_3pcf.h b/modules/integrals_3pcf.h index 903a2496..b5e75053 100644 --- a/modules/integrals_3pcf.h +++ b/modules/integrals_3pcf.h @@ -38,7 +38,7 @@ class Integrals{ mbin = par->mbin; // number of Legendre bins, should be simply max_l+1 out_file = par->out_file; // output directory max_l = par->max_l; // maximum Legendre polynomial to output - max_leg = std::max(max_l, sc->n_param); // maximum Legendre multipole to use + max_leg = std::max(max_l, sc->n_param-1); // maximum Legendre multipole to use, selected among the max ell for the covariance and the max ell for the survey correction function decomposition, the latter being 1 less than the number of parameters (multipoles) there int ec=0; array_len = nbin*nbin*mbin; @@ -339,17 +339,17 @@ class Integrals{ // Define triangle sides independent of i triangle_bins(pj.pos,pk.pos,pl.pos,norm_jkl,ang_jkl,norm_jk,2); + // Save lengths and bins for later + norm_kl = norm_jkl[0]; + // Define first xi function if(index==0){ - los_tmp = compute_los(pk.pos,pl.pos,norm_jkl[0]); - tmp_xi1 = cf->xi(norm_jkl[0],los_tmp); //xi_kl + los_tmp = compute_los(pk.pos, pl.pos, norm_kl); + tmp_xi1 = cf->xi(norm_kl, los_tmp); //xi_kl xi_pass2[0] = tmp_xi1; } else tmp_xi1 = xi_pass[0]; //xi_jk - // Save lengths and bins for later - norm_kl = norm_jkl[0]; - // Compute radial bins for this triangle int x = all_bins(norm_jkl,bins_jkl,bin_jk,2); @@ -462,13 +462,13 @@ class Integrals{ // Define first xi function if(index==0){ - los_tmp = compute_los(pl.pos,pm.pos,norm_klm[0]); - tmp_xi1 = cf->xi(norm_klm[0],los_tmp); // xi_lm + los_tmp = compute_los(pl.pos, pm.pos, norm_lm); + tmp_xi1 = cf->xi(norm_lm, los_tmp); // xi_lm } else{ - Float norm_tmp; - cleanup_l(pj.pos,pm.pos,norm_tmp,los_tmp); - tmp_xi1 = cf->xi(norm_tmp,los_tmp); // xi_jm + Float norm_jm; + cleanup_l(pj.pos, pm.pos, norm_jm, los_tmp); + tmp_xi1 = cf->xi(norm_jm, los_tmp); // xi_jm xi_pass3 = tmp_xi1; // save for next integrator } @@ -767,15 +767,15 @@ class Integrals{ } } - void save_counts(uint64 triple_counts,uint64 quad_counts, uint64 quint_counts, uint64 hex_counts, int index){ + void save_counts(uint64 triple_counts, uint64 quad_counts, uint64 quint_counts, uint64 hex_counts, int index){ // Print the counts for each integral (used for combining the estimates outside of C++) // This is the number of counts used in each loop [always the same] std::string counts_file = string_format("%s3PCFCovMatricesAll/total_counts_n%d_l%d_%d.txt", out_file, nbin, max_l, index); FILE * CountsFile = fopen(counts_file.c_str(), "w"); - fprintf(CountsFile,"%llu\n",triple_counts); - fprintf(CountsFile,"%llu\n",quad_counts); - fprintf(CountsFile,"%llu\n",quint_counts); - fprintf(CountsFile,"%llu\n",hex_counts); + fprintf(CountsFile, "%llu\n", triple_counts); + fprintf(CountsFile, "%llu\n", quad_counts); + fprintf(CountsFile, "%llu\n", quint_counts); + if (index != 0) fprintf(CountsFile, "%llu\n", hex_counts); // only print hex counts for second loop since we don't compute hexes in first loop fflush(NULL); fclose(CountsFile); @@ -790,25 +790,25 @@ class Integrals{ std::string c3name = string_format("%s3PCFCovMatricesAll/c3_n%d_l%d_%d_%s.txt", out_file, nbin, max_l, index, suffix); std::string c4name = string_format("%s3PCFCovMatricesAll/c4_n%d_l%d_%d_%s.txt", out_file, nbin, max_l, index, suffix); - std::string c5name = string_format("%s3PCFCovMatricesAll/c5_n%d_l%d_%d_%s.txt", out_file,nbin, max_l,index, suffix); - std::string c6name = string_format("%s3PCFCovMatricesAll/c6_n%d_l%d_%d_%s.txt", out_file,nbin, max_l,index,suffix); + std::string c5name = string_format("%s3PCFCovMatricesAll/c5_n%d_l%d_%d_%s.txt", out_file, nbin, max_l, index, suffix); + std::string c6name = string_format("%s3PCFCovMatricesAll/c6_n%d_l%d_%d_%s.txt", out_file, nbin, max_l, index, suffix); FILE * C3File = fopen(c3name.c_str(), "w"); // for c3 part of integral FILE * C4File = fopen(c4name.c_str(), "w"); // for c4 part of integral FILE * C5File = fopen(c5name.c_str(), "w"); // for c5 part of integral - FILE * C6File = fopen(c6name.c_str(), "w"); // for c6 part of integral - - for(int i=0;inbin*par->mbin; // define number of bins in total - char line[1000000]; + const size_t line_size = 1000000; + char line[line_size]; int index; int ec2=0; @@ -113,7 +114,7 @@ class JK_weights{ fprintf(stderr,"\nReading RR bin count file '%s'\n",RR_file); index = 0; - while (fgets(line,5000,fp2)!=NULL){ + while (fgets(line, line_size, fp2)!=NULL){ // Select required lines in file if (line[0]=='#') continue; if (line[0]=='\n') continue; @@ -141,7 +142,7 @@ class JK_weights{ fprintf(stderr,"\nReading jackknife file '%s'\n",jk_file); // Count lines to construct the correct size - while (fgets(line,1000000,fp)!=NULL){ + while (fgets(line, line_size, fp)!=NULL){ if (line[0]=='#') continue; // comment line if (line[0]=='\n') continue; n_JK_filled++; @@ -159,8 +160,8 @@ class JK_weights{ int line_count=0; // line counter int counter; // counts which element in line - // Read in values to file - while (fgets(line,1000000,fp)!=NULL) { + // Read in values from file + while (fgets(line, line_size, fp)!=NULL) { // Select required lines in file if (line[0]=='#') continue; if (line[0]=='\n') continue; diff --git a/modules/legendre_mix_utilities.h b/modules/legendre_mix_utilities.h index 70a54314..c30b8bb6 100644 --- a/modules/legendre_mix_utilities.h +++ b/modules/legendre_mix_utilities.h @@ -33,7 +33,8 @@ class MuBinLegendreFactors{ n_l = max_l/2+1; // number of multipoles, even only // READ IN FILE - char line[1000000]; + const size_t line_size = 1000000; + char line[line_size]; int line_no = 0; FILE *fp; @@ -46,7 +47,7 @@ class MuBinLegendreFactors{ fprintf(stderr, "\nMu bin Legendre factors file '%s'\n", filename); // Count lines to construct the correct size - while (fgets(line, 1000000, fp) != NULL) { + while (fgets(line, line_size, fp) != NULL) { if (line[0] == '#') continue; // comment line if (line[0] == '\n') continue; line_no++; @@ -62,7 +63,7 @@ class MuBinLegendreFactors{ int index = 0; // indexes array // Read in values from file - while (fgets(line, 1000000, fp) != NULL) { + while (fgets(line, line_size, fp) != NULL) { // Select required lines in file if (line[0] == '#') continue; diff --git a/modules/legendre_utilities.h b/modules/legendre_utilities.h index 5721c8c3..7d8d2cd7 100644 --- a/modules/legendre_utilities.h +++ b/modules/legendre_utilities.h @@ -92,7 +92,7 @@ class SurveyCorrection{ public: #ifdef THREE_PCF Float* inv_phi_multipoles = nullptr; // houses multipole moments of the inverse correction function - int n_param = 7; + int n_param = 7; // number of multipole coefficients in the Legendre polynomial decomposition of the inverse correction function #else Float* phi_coeffs = nullptr; // houses polynomial coefficients for the correction function int n_param = 7; // number of polynomial coefficients in fitting model (only 5 are independent) @@ -140,8 +140,8 @@ class SurveyCorrection{ #ifdef THREE_PCF SurveyCorrection(Parameters *par){ - // This initializes the function and reads in the relevant polynomial coefficients for each radial bin. - // NB: coefficients are indexed as INDEX = RADIAL_BIN_1*N_COEFF^2. + RADIAL_BIN_2*N_COEFF + COEFF_ID where N_COEFF is the total number of coefficients for each model; here 7. + // This initializes the function and reads in the relevant polynomial coefficients for each radial bin. + // NB: coefficients are indexed as INDEX = RADIAL_BIN_1*N_COEFF*NBIN + RADIAL_BIN_2*N_COEFF + COEFF_ID where N_COEFF is the total number of coefficients for each model; here 7. #else SurveyCorrection(Parameters *par, int index1, int index2){ // This initializes the function and reads in the relevant polynomial coefficients for each radial bin @@ -151,7 +151,8 @@ class SurveyCorrection{ nbin = par->nbin; // number of radial bins // READ IN FILE - char line[1000000], *phi_file; + const size_t line_size = 1000000; + char line[line_size], *phi_file; int line_no = 0; FILE *fp; @@ -171,7 +172,7 @@ class SurveyCorrection{ fprintf(stderr,"\nReading survey correction function coefficient file '%s'\n",phi_file); // Count lines to construct the correct size - while (fgets(line,1000000,fp)!=NULL){ + while (fgets(line, line_size, fp)!=NULL){ if (line[0]=='#') continue; // comment line if (line[0]=='\n') continue; line_no++; @@ -191,8 +192,8 @@ class SurveyCorrection{ int line_count=0; // line counter int index=0; // indexes array - // Read in values to file - while (fgets(line,1000000,fp)!=NULL) { + // Read in values from file + while (fgets(line, line_size, fp)!=NULL) { // Select required lines in file if (line[0]=='#') continue; diff --git a/modules/parameters.h b/modules/parameters.h index ee4dc3df..9e595c9d 100644 --- a/modules/parameters.h +++ b/modules/parameters.h @@ -104,10 +104,12 @@ class Parameters{ // Maximum number of iterations to compute the C_ab integrals over int max_loops = 40; +#ifndef TRIPLE // Number of loops to output into each subsample/file int loops_per_sample = 1; // Number of output subsamples/files int no_subsamples = 120; +#endif // Number of random cells to draw at each stage int N2 = 20; // number of j cells per i cell @@ -135,10 +137,16 @@ class Parameters{ // Number of galaxies in second dataset Float nofznorm2=3398430; // - // parameters allowing to run only some of the 7 integrals +#ifndef TRIPLE + // parameters allowing to run only some of the 7 integrals for multi-tracer 2PCF or 2 for single-tracer 3PCF int start_integral_index = 1; // the first integral to be computed +#ifdef THREE_PCF + int last_integral_index = 2; // the last integral to be computed +#else int last_integral_index = 7; // the last integral to be computed +#endif // defaults compute all integrals +#endif //---------- (r,mu) MULTI-FIELD PARAMETERS ------------------------------ @@ -191,7 +199,7 @@ class Parameters{ //-------- OTHER PARAMETERS ---------------------------------------------- // The minimum mu of the smallest bin. -#ifdef THREE_PCF +#if (defined THREE_PCF || defined TRIPLE) Float mumin = -1.0; #else Float mumin = 0.0; @@ -267,7 +275,9 @@ class Parameters{ rect_boxsize = {boxsize, boxsize, boxsize}; } else if (!strcmp(argv[i],"-maxloops")) max_loops = atoi(argv[++i]); +#ifndef TRIPLE else if (!strcmp(argv[i],"-loopspersample")) loops_per_sample = atoi(argv[++i]); +#endif else if (!strcmp(argv[i],"-rescale")) rescale = atof(argv[++i]); else if (!strcmp(argv[i],"-mumax")) mumax = atof(argv[++i]); else if (!strcmp(argv[i],"-mumin")) mumin = atof(argv[++i]); @@ -296,8 +306,10 @@ class Parameters{ else if (!strcmp(argv[i],"-binfile_cf")) radial_bin_file_cf=argv[++i]; else if (!strcmp(argv[i],"-N2")) N2=atof(argv[++i]); else if (!strcmp(argv[i],"-N3")) N3=atof(argv[++i]); +#ifndef TRIPLE else if (!strcmp(argv[i],"-N4")) N4=atof(argv[++i]); -#if (defined LEGENDRE || defined LEGENDRE_MIX) +#endif +#if (defined LEGENDRE || defined LEGENDRE_MIX || defined THREE_PCF || defined POWER) else if (!strcmp(argv[i],"-max_l")) max_l=atoi(argv[++i]); #endif #ifdef JACKKNIFE @@ -312,7 +324,6 @@ class Parameters{ else if (!strcmp(argv[i],"-phi_file12")) phi_file12=argv[++i]; else if (!strcmp(argv[i],"-phi_file2")) phi_file2=argv[++i]; #elif defined POWER - else if (!strcmp(argv[i],"-max_l")) max_l=atoi(argv[++i]); else if (!strcmp(argv[i],"-inv_phi_file")) inv_phi_file=argv[++i]; else if (!strcmp(argv[i],"-inv_phi_file12")) inv_phi_file12=argv[++i]; else if (!strcmp(argv[i],"-inv_phi_file2")) inv_phi_file2=argv[++i]; @@ -321,16 +332,17 @@ class Parameters{ else if (!strcmp(argv[i],"-power_norm12")) power_norm12 = atof(argv[++i]); else if (!strcmp(argv[i],"-power_norm2")) power_norm2 = atof(argv[++i]); #elif defined THREE_PCF - else if (!strcmp(argv[i],"-max_l")) max_l=atoi(argv[++i]); else if (!strcmp(argv[i],"-phi_file")) phi_file=argv[++i]; else if (!strcmp(argv[i],"-N5")) N5=atof(argv[++i]); else if (!strcmp(argv[i],"-N6")) N6=atof(argv[++i]); #endif #if (!defined LEGENDRE && !defined POWER && !defined THREE_PCF) else if (!strcmp(argv[i],"-mbin")) mbin = atoi(argv[++i]); +#ifndef TRIPLE else if (!strcmp(argv[i],"-RRbin")) RR_bin_file=argv[++i]; else if (!strcmp(argv[i],"-RRbin12")) RR_bin_file12=argv[++i]; else if (!strcmp(argv[i],"-RRbin2")) RR_bin_file2=argv[++i]; +#endif #endif else if (!strcmp(argv[i],"-perbox")) perbox = 1; else if (!strcmp(argv[i],"-np")) { @@ -344,8 +356,10 @@ class Parameters{ make_random=1; } else if (!strcmp(argv[i],"-seed")) { random_seed = false; sscanf(argv[++i], "%lu", &seed); } +#ifndef TRIPLE else if (!strcmp(argv[i],"-start_integral_index")) start_integral_index = atoi(argv[++i]); else if (!strcmp(argv[i],"-last_integral_index")) last_integral_index = atoi(argv[++i]); +#endif else if (!strcmp(argv[i],"-def")) { fname = NULL; } else { fprintf(stderr, "Don't recognize %s\n", argv[i]); @@ -375,10 +389,12 @@ class Parameters{ assert(nside%2!=0); // The probability integrator needs an odd grid size assert(nofznorm>0); // need some galaxies! +#ifndef TRIPLE assert(max_loops % loops_per_sample == 0); // group size need to divide the number of loops no_subsamples = max_loops / loops_per_sample; +#endif #ifndef THREE_PCF - //assert(mumin>=0); // We take the absolte value of mu + //assert(mumin>=0); // We take the absolute value of mu #endif assert(mumax<=1); // mu > 1 makes no sense @@ -535,13 +551,19 @@ class Parameters{ #endif // check integral selection parameters +#ifndef TRIPLE assert(start_integral_index >= 1); assert(last_integral_index >= 1); - if (!multi_tracers) last_integral_index = 1; - assert(start_integral_index <= last_integral_index); +#ifdef THREE_PCF + int no_integrals = 2; // there are 2 integrals for single-tracer 3PCF +#else int no_integrals = multi_tracers ? 7 : 1; + if (!multi_tracers) last_integral_index = 1; // only 1 integral for single-tracer 2PCF +#endif + assert(start_integral_index <= last_integral_index); assert(start_integral_index <= no_integrals); assert(last_integral_index <= no_integrals); +#endif create_directory(); @@ -573,7 +595,9 @@ class Parameters{ #endif printf("Number of galaxies = %6.5e\n",nofznorm); printf("Maximum number of integration loops = %d\n",max_loops); +#ifndef TRIPLE printf("Number of output subsamples = %d\n", no_subsamples); +#endif printf("Output directory: '%s'\n",out_file); } @@ -595,7 +619,11 @@ class Parameters{ private: void usage() { - fprintf(stderr, "\nUsage for grid_covariance:\n\n"); +#ifdef TRIPLE + fprintf(stderr, "\nUsage for triple (triple_counts):\n\n"); +#else + fprintf(stderr, "\nUsage for cov (grid_covariance):\n\n"); +#endif fprintf(stderr, " -def: This allows one to accept the defaults without giving other entries.\n"); fprintf(stderr, " -in : The input random particle file for particle-set 1 (space-separated x,y,z,w).\n"); fprintf(stderr, " -binfile : File containing the desired radial bins\n"); @@ -604,10 +632,12 @@ class Parameters{ fprintf(stderr, " -norm : Number of galaxies in the first tracer set.\n"); fprintf(stderr, " -effective_norm: Flag indicating to use the effective number of random particles (sum(w)^2/sum(w^2), taking into account weights but invariant to overall weight rescaling) instead of simple counting (default). Should match the way the -norm value was computed for the data galaxies externally.\n"); #ifdef JACKKNIFE - fprintf(stderr, " -jackknife : File containing the {1,1} jackknife weights (normally computed from Corrfunc)\n"); + fprintf(stderr, " -jackknife : File containing the {1,1} jackknife weights\n"); #endif #if (!defined LEGENDRE && !defined POWER && !defined THREE_PCF) - fprintf(stderr, " -RRbin : File containing the {1,1} jackknife RR bin counts (computed from Corrfunc)\n"); +#ifndef TRIPLE + fprintf(stderr, " -RRbin : File containing the {1,1} RR bin counts\n"); +#endif fprintf(stderr, " -mbin : The number of mu bins (spaced linearly).\n"); #endif fprintf(stderr, " -mbin_cf : The number of mu bins in the correlation function (spaced linearly).\n"); @@ -624,7 +654,7 @@ class Parameters{ fprintf(stderr, " -cor2 : (Optional) File location of input xi_2 correlation function file.\n"); fprintf(stderr, " -norm2 : (Optional) Number of galaxies in the survey for the second tracer set.\n"); -#if (defined LEGENDRE || defined LEGENDRE_MIX) +#if (defined LEGENDRE || defined LEGENDRE_MIX || defined THREE_PCF || defined POWER) fprintf(stderr, " -max_l : Maximum legendre multipole (must be even)\n"); #endif #ifdef LEGENDRE_MIX @@ -646,9 +676,9 @@ class Parameters{ fprintf(stderr, " -power_norm12: Power spectrum normalization = V*<(nw)^2> = Sum(nw^2) for field 1 x 2\n"); fprintf(stderr, " -power_norm2: Power spectrum normalization = V*<(nw)^2> = Sum(nw^2) for field 2\n"); #endif -#if (!defined LEGENDRE && !defined POWER && !defined THREE_PCF) - fprintf(stderr, " -RRbin12 : (Optional) File containing the {1,2} jackknife RR bin counts (computed from Corrfunc)\n"); - fprintf(stderr, " -RRbin2 : (Optional) File containing the {2,2} jackknife RR bin counts (computed from Corrfunc)\n"); +#if (!defined LEGENDRE && !defined POWER && !defined THREE_PCF && !defined TRIPLE) + fprintf(stderr, " -RRbin12 : (Optional) File containing the {1,2} RR bin counts\n"); + fprintf(stderr, " -RRbin2 : (Optional) File containing the {2,2} RR bin counts\n"); fprintf(stderr, "\n"); #endif #ifdef JACKKNIFE @@ -656,10 +686,14 @@ class Parameters{ fprintf(stderr, " -jackknife2 : (Optional) File containing the {2,2} jackknife weights (normally computed from Corrfunc)\n"); #endif fprintf(stderr, " -maxloops : Maximum number of integral loops\n"); +#ifndef TRIPLE fprintf(stderr, " -loopspersample : Number of loops to collapse into each subsample. Default 1.\n"); +#endif fprintf(stderr, " -N2 : Number of secondary particles to choose per primary particle\n"); fprintf(stderr, " -N3 : Number of tertiary particles to choose per secondary particle\n"); +#ifndef TRIPLE fprintf(stderr, " -N4 : Number of quaternary particles to choose per tertiary particle\n"); +#endif fprintf(stderr, "\n"); #ifdef THREE_PCF fprintf(stderr, " -N5 : Number of fifth particles to choose per quaternary particle\n"); @@ -687,10 +721,17 @@ class Parameters{ fprintf(stderr, " -seed : If given, allows to reproduce the results with the same settings, except the number of threads.\n"); fprintf(stderr, " Individual subsamples may differ because they are accumulated/written in order of loop completion which may depend on external factors at runtime, but the final integrals should be the same.\n"); fprintf(stderr, " Needs to be a non-negative integer fitting into 32 bits (max 4294967295). If not given (default), the seed will be created randomly.\n"); +#ifdef THREE_PCF + fprintf(stderr, " -start_integral_index : If given, chooses from which of 2 integrals (numbered 1 and 2) to stop for 3PCF. Default 1 (to cover all integrals).\n"); + fprintf(stderr, " Intended to help complete the timed-out multi-tracer runs by skipping the integrals computed in an unfinished run.\n"); + fprintf(stderr, " -last_integral_index : If given, chooses at which of 2 integrals (numbered 1 and 2) to stop for 3PCF. Default 2 (to cover all integrals).\n"); + fprintf(stderr, " Can be used to fit into the time limit if the full 2 integrals are expected to take longer.\n"); +#elif (!defined TRIPLE) fprintf(stderr, " -start_integral_index : If given, chooses from which of 7 integrals (numbered 1 through 7) to start in the multi-tracer regime. Default 1 (to cover all integrals).\n"); fprintf(stderr, " Intended to help complete the timed-out multi-tracer runs by skipping the integrals computed in an unfinished run.\n"); fprintf(stderr, " -last_integral_index : If given, chooses at which of 7 integrals (numbered 1 through 7) to stop in the multi-tracer regime. Default 7 (to cover all integrals).\n"); fprintf(stderr, " Can be used to fit into the time limit if the full 7 integrals are expected to take longer.\n"); +#endif fprintf(stderr, "\n"); fprintf(stderr, "\n"); @@ -711,6 +752,7 @@ class Parameters{ if (mkdir(out_file,S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH)==0){ printf("\nCreating output directory\n"); } +#ifndef TRIPLE #ifdef THREE_PCF std::string cname = string_format("%s3PCFCovMatricesAll/", out_file); #elif defined POWER @@ -734,6 +776,7 @@ class Parameters{ printf("\nCreation of directory %s failed\n", cjname.c_str()); exit(1); } +#endif #endif } diff --git a/power_spectra/power_mod/survey_correction_legendre.h b/power_spectra/power_mod/survey_correction_legendre.h index eb90188a..631e2dc2 100644 --- a/power_spectra/power_mod/survey_correction_legendre.h +++ b/power_spectra/power_mod/survey_correction_legendre.h @@ -74,7 +74,8 @@ class SurveyCorrection{ // NB: coefficients are indexed as INDEX = MULTIPOLE_INDEX*N_COEFF. + COEFF_ID where N_COEFF is the total number of coefficients for each multipole; here 3. // READ IN FILE - char line[1000000], *phi_file; + const size_t line_size = 1000000; + char line[line_size], *phi_file; int line_no = 0; FILE *fp; @@ -91,7 +92,7 @@ class SurveyCorrection{ printf("\nReading survey correction function coefficient file '%s'\n",phi_file); // Count lines to construct the correct size - while (fgets(line,1000000,fp)!=NULL){ + while (fgets(line, line_size, fp)!=NULL){ if (line[0]=='#') continue; // comment line if (line[0]=='\n') continue; line_no++; @@ -107,8 +108,8 @@ class SurveyCorrection{ int line_count=0; // line counter int index=0; // indexes array - // Read in values to file - while (fgets(line,1000000,fp)!=NULL) { + // Read in values from file + while (fgets(line, line_size, fp)!=NULL) { // Select required lines in file if (line[0]=='#') continue; diff --git a/pyproject.toml b/pyproject.toml index 35fb2dfe..a6b0466e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -51,4 +51,4 @@ path = "RascalC/__init__.py" out_dir = "./RascalC/bin" work_dir = "./RascalC" commands = ["make"] -artifacts = ["RascalC/bin/cov.*"] +artifacts = ["RascalC/bin/cov.*", "RascalC/bin/triple.*"] diff --git a/scripts/cat_raw_covariance_matrices.py b/scripts/cat_raw_covariance_matrices.py index 5c7d255b..f86375df 100644 --- a/scripts/cat_raw_covariance_matrices.py +++ b/scripts/cat_raw_covariance_matrices.py @@ -1,24 +1,21 @@ -## Script to catenate (take a subset of subsamples or resample) raw covariance matrix files produced by collect_raw_covariance_matrices.py +## Script to catenate (take a subset of subsamples or resample) raw covariance matrix files ## More specifically, copy part of partial results to other directory and recompute totals by averaging ## Determines single-field vs multi-field and jackknife automatically ## Do not use if subsamples have different numbers of pairs/triples/quadruplets -import sys +import argparse -# PARAMETERS -if len(sys.argv)<6: # if too few - print("Usage: python cat_raw_covariance_matrices.py {N_R_BINS} {mN_MU_BINS/lMAX_L} {COVARIANCE_INPUT_DIR1} {N_SUBSAMPLES_TO_USE1} [{COVARIANCE_INPUT_DIR2} {N_SUBSAMPLES_TO_USE2} ...] [{COLLAPSE_FACTOR}] {COVARIANCE_OUTPUT_DIR}") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to catenate (take a subset of subsamples or resample) raw covariance matrix files or directories. More specifically, copy part of partial results to other directory and recompute totals by averaging. Determines single-field vs multi-field and jackknife automatically. Do not use if subsamples have different numbers of pairs/triples/quadruplets.") +parser.add_argument("n_r_bins", type=int, help="number of radial/separation bins") +parser.add_argument("n_mu_bins_or_max_l_specifier", type=int, help="string specifying the number of angular/mu bins or max_l (Legendre), the format is m{n_mu_bins} or l{max_l} accordingly") +parser.add_argument("input_covariance_dirs", type=str, nargs='+', help="one directory or multiple directories to read the covariance matrix subsamples from") +parser.add_argument("--ns_subsamples", type=int, nargs='*', default=None, help="number of subsamples to use, one number for each covariance directory") +parser.add_argument("--collapse_factor", type=int, default=1, help="reduce the number of samples by this factor by averaging batches of this size (default 1 = no reduction)") +parser.add_argument("output_covariance_dir", type=str, help="directory to write the resulting covariance matrices and subsamples") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.raw_covariance_matrices import cat_raw_covariance_matrices -n = int(sys.argv[1]) -mstr = str(sys.argv[2]) -input_roots = [str(s) for s in sys.argv[3:-1:2]] -ns_samples = [int(s) for s in sys.argv[4:-1:2]] -output_root = str(sys.argv[-1]) -collapse_factor = int(input_roots.pop()) if len(input_roots) > len(ns_samples) else 1 # recover the collapse factor if present - -cat_raw_covariance_matrices(n, mstr, input_roots, ns_samples, output_root, collapse_factor) \ No newline at end of file +cat_raw_covariance_matrices(args.n, args.mstr, args.input_covariance_dirs, args.ns_subsamples or [None] * len(args.input_covariance_dirs), args.output_covariance_dir, args.collapse_factor) \ No newline at end of file diff --git a/scripts/collect_raw_covariance_matrices.py b/scripts/collect_raw_covariance_matrices.py index ec944a41..1043199a 100644 --- a/scripts/collect_raw_covariance_matrices.py +++ b/scripts/collect_raw_covariance_matrices.py @@ -1,16 +1,13 @@ ## Script to collect the raw covariance matrices from the output directory of the C++ code -import sys +import argparse -# PARAMETERS -if len(sys.argv) != 2: # if too few - print("Usage: python collect_raw_covariance_matrices.py {COVARIANCE_DIR}") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to collect the raw covariance matrices from the output directory of the C++ code. It converts multiple text files to a single binary (Numpy) file.") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.raw_covariance_matrices import collect_raw_covariance_matrices -cov_dir = str(sys.argv[1]) - -collect_raw_covariance_matrices(cov_dir) \ No newline at end of file +collect_raw_covariance_matrices(args.covariance_dir) \ No newline at end of file diff --git a/scripts/combine_covs.py b/scripts/combine_covs.py index ad7280b7..43344d8f 100644 --- a/scripts/combine_covs.py +++ b/scripts/combine_covs.py @@ -1,25 +1,20 @@ "This reads two sets of RascalC results and two cosmodesi/pycorr .npy files to combine two covs following NS/GCcomb procedure. Covariance of N and S 2PCF is neglected." -import sys +import argparse -## PARAMETERS -if len(sys.argv) not in (9, 11): - print("Usage: python combine_covs.py {RASCALC_RESULTS1} {RASCALC_RESULTS2} {PYCORR_FILE1} {PYCORR_FILE2} {R_STEP} {N_MU_BINS} {SKIP_R_BINS} {OUTPUT_COV_FILE} [{OUTPUT_COV_FILE1} {OUTPUT_COV_FILE2}].") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads two sets of RascalC results (in radial/separation and angular/mu bins) and two cosmodesi/pycorr .npy files to combine two covs following NS/GCcomb procedure. Covariance of N and S 2PCF is neglected.") +parser.add_argument("rascalc_files", type=str, nargs=2, help="RascalC post-processed .npy filenames for the two regions") +parser.add_argument("pycorr_files", type=str, nargs=2, help="pycorr .npy filenames for the two regions") +parser.add_argument("r_step", type=float, help="desired width of radial bins") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins (0 means to preserve the original number)") +parser.add_argument("skip_r_bins", type=int, help="number of last radial/separation bins to discard") +parser.add_argument("output_cov_file", type=str, help="output text file for the covariance in the combined region") +parser.add_argument("output_cov_file1", type=str, default=None, nargs='?', help="optional output text file for the covariance in the first region") +parser.add_argument("output_cov_file2", type=str, default=None, nargs='?', help="optional output text file for the covariance in the second region") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.comb.combine_covs import combine_covs -rascalc_results1 = str(sys.argv[1]) -rascalc_results2 = str(sys.argv[2]) -pycorr_file1 = str(sys.argv[3]) -pycorr_file2 = str(sys.argv[4]) -r_step = float(sys.argv[5]) -n_mu_bins = int(sys.argv[6]) -skip_r_bins = int(sys.argv[7]) -output_cov_file = str(sys.argv[8]) -output_cov_file1 = get_arg_safe(9, str, None) -output_cov_file2 = get_arg_safe(10, str, None) - -combine_covs(rascalc_results1, rascalc_results2, pycorr_file1, pycorr_file2, output_cov_file, n_mu_bins, r_step, skip_r_bins, output_cov_file1, output_cov_file2) \ No newline at end of file +combine_covs(args.rascalc_files[0], args.rascalc_files[1], args.pycorr_files[0], args.pycorr_files[1], args.output_cov_file, args.n_mu_bins or None, args.r_step, args.skip_r_bins, args.output_cov_file1, args.output_cov_file2) \ No newline at end of file diff --git a/scripts/combine_covs_legendre.py b/scripts/combine_covs_legendre.py index de160622..5cf240f6 100644 --- a/scripts/combine_covs_legendre.py +++ b/scripts/combine_covs_legendre.py @@ -1,25 +1,20 @@ "This reads two sets of RascalC results and two cosmodesi/pycorr .npy files to combine two covs following NS/GCcomb procedure in Legendre mode. Covariance of N and S 2PCF is neglected." -import sys +import argparse -## PARAMETERS -if len(sys.argv) not in (9, 11): - print("Usage: python combine_covs_legendre.py {RASCALC_RESULTS1} {RASCALC_RESULTS2} {PYCORR_FILE1} {PYCORR_FILE2} {R_STEP} {MAX_L} {SKIP_R_BINS} {OUTPUT_COV_FILE} [{OUTPUT_COV_FILE1} {OUTPUT_COV_FILE2}].") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads two sets of RascalC results and two cosmodesi/pycorr .npy files to combine two covs following NS/GCcomb procedure in Legendre mode. Covariance of N and S 2PCF is neglected.") +parser.add_argument("rascalc_files", type=str, nargs=2, help="RascalC post-processed .npy filenames for the two regions") +parser.add_argument("pycorr_files", type=str, nargs=2, help="pycorr .npy filenames for the two regions") +parser.add_argument("r_step", type=float, default=1, help="desired width of radial bins (default 1)") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("skip_r_bins", type=int, default=0, help="number of last radial/separation bins to discard") +parser.add_argument("output_cov_file", type=str, help="output text file for the covariance in the combined region") +parser.add_argument("output_cov_file1", type=str, default=None, nargs='?', help="optional output text file for the covariance in the first region") +parser.add_argument("output_cov_file2", type=str, default=None, nargs='?', help="optional output text file for the covariance in the second region") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.comb.combine_covs_legendre import combine_covs_legendre -rascalc_results1 = str(sys.argv[1]) -rascalc_results2 = str(sys.argv[2]) -pycorr_file1 = str(sys.argv[3]) -pycorr_file2 = str(sys.argv[4]) -r_step = float(sys.argv[5]) -max_l = int(sys.argv[6]) -skip_r_bins = int(sys.argv[7]) -output_cov_file = str(sys.argv[8]) -output_cov_file1 = get_arg_safe(9, str, None) -output_cov_file2 = get_arg_safe(10, str, None) - -combine_covs_legendre(rascalc_results1, rascalc_results2, pycorr_file1, pycorr_file2, output_cov_file, max_l, r_step, skip_r_bins, output_cov_file1, output_cov_file2) \ No newline at end of file +combine_covs_legendre(args.rascalc_files[0], args.rascalc_files[1], args.pycorr_files[0], args.pycorr_files[1], args.output_cov_file, args.max_l, args.r_step, args.skip_r_bins, args.output_cov_file1, args.output_cov_file2) \ No newline at end of file diff --git a/scripts/combine_covs_legendre_multi.py b/scripts/combine_covs_legendre_multi.py index 82d879c1..4b40229f 100644 --- a/scripts/combine_covs_legendre_multi.py +++ b/scripts/combine_covs_legendre_multi.py @@ -1,25 +1,20 @@ "This reads two sets of RascalC results and two triplets of cosmodesi/pycorr .npy files to combine two full 2-tracer covs following NS/GCcomb procedure for 2 tracers in Legendre mode. Covariance of N(GC) and S(GC) 2PCF is neglected." -import sys +import argparse -## PARAMETERS -if len(sys.argv) not in (13, 15): - print("Usage: python combine_covs_legendre_multi.py {RASCALC_RESULTS1} {RASCALC_RESULTS2} {PYCORR_FILE1_11} {PYCORR_FILE2_11} {PYCORR_FILE1_12} {PYCORR_FILE2_12} {PYCORR_FILE1_22} {PYCORR_FILE2_22} {N_R_BINS} {MAX_L} {SKIP_R_BINS} {OUTPUT_COV_FILE} [{OUTPUT_COV_FILE1} {OUTPUT_COV_FILE2}].") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads two sets of RascalC results and two cosmodesi/pycorr .npy files to combine two covs following NS/GCcomb procedure for 2 tracers in Legendre mode. Covariance of N and S 2PCF is neglected.") +parser.add_argument("rascalc_files", type=str, nargs=2, help="RascalC post-processed .npy filenames for the two regions") +parser.add_argument("pycorr_files", type=str, nargs=6, help="6 pycorr .npy filenames for the two regions and 3 distinct correlation functions in the following order: first tracer auto-counts for both regions, then cross-counts between the two tracers for both regions, and finally the second tracer auto-counts for both regions") +parser.add_argument("r_step", type=float, help="desired width of radial bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("skip_r_bins", type=int, help="number of last radial/separation bins to discard") +parser.add_argument("output_cov_file", type=str, help="output text file for the covariance in the combined region") +parser.add_argument("output_cov_file1", type=str, default=None, nargs='?', help="optional output text file for the covariance in the first region") +parser.add_argument("output_cov_file2", type=str, default=None, nargs='?', help="optional output text file for the covariance in the second region") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.comb.combine_covs_legendre_multi import combine_covs_legendre_multi -rascalc_results1 = str(sys.argv[1]) -rascalc_results2 = str(sys.argv[2]) -pycorr_files1 = [str(sys.argv[i]) for i in (3, 5, 7)] -pycorr_files2 = [str(sys.argv[i]) for i in (4, 6, 8)] -r_step = float(sys.argv[9]) -max_l = int(sys.argv[10]) -skip_r_bins = int(sys.argv[11]) -output_cov_file = str(sys.argv[12]) -output_cov_file1 = get_arg_safe(13, str, None) -output_cov_file2 = get_arg_safe(14, str, None) - -combine_covs_legendre_multi(rascalc_results1, rascalc_results2, pycorr_files1, pycorr_files2, output_cov_file, max_l, r_step, skip_r_bins, output_cov_file1, output_cov_file2) \ No newline at end of file +combine_covs_legendre_multi(args.rascalc_files[0], args.rascalc_files[1], [args.pycorr_files[i] for i in range(0, 6, 2)], [args.pycorr_files[i] for i in range(1, 6, 2)], args.output_cov_file, args.max_l, args.r_step, args.skip_r_bins, args.output_cov_file1, args.output_cov_file2) \ No newline at end of file diff --git a/scripts/combine_covs_legendre_multi_to_cat.py b/scripts/combine_covs_legendre_multi_to_cat.py index 4f934281..f0074587 100644 --- a/scripts/combine_covs_legendre_multi_to_cat.py +++ b/scripts/combine_covs_legendre_multi_to_cat.py @@ -1,27 +1,22 @@ -"This reads two sets of RascalC results and two triplets of cosmodesi/pycorr .npy files to combine two full 2-tracer covs following NS/GCcomb procedure for 2 tracers in Legendre mode and convert them . Covariance of N(GC) and S(GC) 2PCF is neglected." +"This reads two sets of RascalC results and two triplets of cosmodesi/pycorr .npy files to combine two full 2-tracer covs following NS/GCcomb procedure for 2 tracers in Legendre mode and convert them to the covariance for the single combined tracer (obtained by additionally weighting each tracer by its bias and concatenating the catalogs). Covariance of N(GC) and S(GC) 2PCF is neglected." -import sys +import argparse -## PARAMETERS -if len(sys.argv) not in (13, 15, 17): - print("Usage: python combine_covs_multi_to_cat.py {RASCALC_RESULTS1} {RASCALC_RESULTS2} {PYCORR_FILE1_11} {PYCORR_FILE2_11} {PYCORR_FILE1_12} {PYCORR_FILE2_12} {PYCORR_FILE1_22} {PYCORR_FILE2_22} {R_STEP} {MAX_L} {R_BINS_SKIP} {OUTPUT_COV_FILE} [{BIAS1} {BIAS2} [{OUTPUT_COV_FILE1} {OUTPUT_COV_FILE2}]].") - sys.exit(1) +parser = argparse.ArgumentParser(description="This reads two sets of RascalC results and two triplets of cosmodesi/pycorr .npy files to combine two full 2-tracer covs following NS/GCcomb procedure for 2 tracers in Legendre mode and convert them to the covariance for the single combined tracer (obtained by additionally weighting each tracer by its bias and concatenating the catalogs). Covariance of N(GC) and S(GC) 2PCF is neglected.") +parser.add_argument("rascalc_files", type=str, nargs=2, help="RascalC post-processed .npy filenames for the two regions") +parser.add_argument("pycorr_files", type=str, nargs=6, help="6 pycorr .npy filenames for the two regions and 3 distinct correlation functions in the following order: first tracer auto-counts for both regions, then cross-counts between the two tracers for both regions, and finally the second tracer auto-counts for both regions") +parser.add_argument("r_step", type=float, help="desired width of radial bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("skip_r_bins", type=int, help="number of last radial/separation bins to discard") +parser.add_argument("output_cov_file", type=str, help="output text file for the combined tracer covariance in the combined region") +parser.add_argument("bias1", type=float, default=1, nargs='?', help="linear bias for the first tracer (default 1)") +parser.add_argument("bias2", type=float, default=1, nargs='?', help="linear bias for the second tracer (default 1)") +parser.add_argument("output_cov_file1", type=str, default=None, nargs='?', help="optional output text file for the combined tracer covariance in the first region") +parser.add_argument("output_cov_file2", type=str, default=None, nargs='?', help="optional output text file for the combined tracer covariance in the second region") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.comb.combine_covs_legendre_multi_to_cat import combine_covs_legendre_multi_to_cat -rascalc_results1 = str(sys.argv[1]) -rascalc_results2 = str(sys.argv[2]) -pycorr_files1 = [str(sys.argv[i]) for i in (3, 5, 7)] -pycorr_files2 = [str(sys.argv[i]) for i in (4, 6, 8)] -r_step = float(sys.argv[9]) -max_l = int(sys.argv[10]) -skip_r_bins = int(sys.argv[11]) -output_cov_file = str(sys.argv[12]) -bias1 = get_arg_safe(13, float, 1) -bias2 = get_arg_safe(14, float, 1) -output_cov_file1 = get_arg_safe(15, str, None) -output_cov_file2 = get_arg_safe(16, str, None) - -combine_covs_legendre_multi_to_cat(rascalc_results1, rascalc_results2, pycorr_files1, pycorr_files2, output_cov_file, max_l, r_step, skip_r_bins, bias1, bias2, output_cov_file1, output_cov_file2) \ No newline at end of file +combine_covs_legendre_multi_to_cat(args.rascalc_files[0], args.rascalc_files[1], [args.pycorr_files[i] for i in range(0, 6, 2)], [args.pycorr_files[i] for i in range(1, 6, 2)], args.output_cov_file, args.max_l, args.r_step, args.skip_r_bins, args.bias1, args.bias2, args.output_cov_file1, args.output_cov_file2) \ No newline at end of file diff --git a/scripts/combine_covs_multi.py b/scripts/combine_covs_multi.py index 2f5089cd..c3968ccb 100644 --- a/scripts/combine_covs_multi.py +++ b/scripts/combine_covs_multi.py @@ -1,25 +1,20 @@ "This reads two sets of RascalC results and two triplets of cosmodesi/pycorr .npy files to combine two covs following NS/GCcomb procedure for 2 tracers. Covariance of N and S 2PCF is neglected." -import sys +import argparse -## PARAMETERS -if len(sys.argv) not in (13, 15): - print("Usage: python combine_covs_multi.py {RASCALC_RESULTS1} {RASCALC_RESULTS2} {PYCORR_FILE1_11} {PYCORR_FILE2_11} {PYCORR_FILE1_12} {PYCORR_FILE2_12} {PYCORR_FILE1_22} {PYCORR_FILE2_22} {R_STEP} {N_MU_BINS} {SKIP_R_BINS} {OUTPUT_COV_FILE} [{OUTPUT_COV_FILE1} {OUTPUT_COV_FILE2}].") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads two sets of RascalC results (in radial/separation and angular/mu bins) and two cosmodesi/pycorr .npy files to combine two covs following NS/GCcomb procedure for 2 tracers. Covariance of N and S 2PCF is neglected.") +parser.add_argument("rascalc_files", type=str, nargs=2, help="RascalC post-processed .npy filenames for the two regions") +parser.add_argument("pycorr_files", type=str, nargs=6, help="6 pycorr .npy filenames for the two regions and 3 distinct correlation functions in the following order: first tracer auto-counts for both regions, then cross-counts between the two tracers for both regions, and finally the second tracer auto-counts for both regions") +parser.add_argument("r_step", type=float, help="desired width of radial bins") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins (0 means to preserve the original number)") +parser.add_argument("skip_r_bins", type=int, help="number of last radial/separation bins to discard") +parser.add_argument("output_cov_file", type=str, help="output text file for the covariance in the combined region") +parser.add_argument("output_cov_file1", type=str, default=None, nargs='?', help="optional output text file for the covariance in the first region") +parser.add_argument("output_cov_file2", type=str, default=None, nargs='?', help="optional output text file for the covariance in the second region") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.comb.combine_covs_multi import combine_covs_multi -rascalc_results1 = str(sys.argv[1]) -rascalc_results2 = str(sys.argv[2]) -pycorr_files1 = [str(sys.argv[i]) for i in (3, 5, 7)] -pycorr_files2 = [str(sys.argv[i]) for i in (4, 6, 8)] -r_step = float(sys.argv[9]) -n_mu_bins = int(sys.argv[10]) -skip_r_bins = int(sys.argv[11]) -output_cov_file = str(sys.argv[12]) -output_cov_file1 = get_arg_safe(13, str, None) -output_cov_file2 = get_arg_safe(14, str, None) - -combine_covs_multi(rascalc_results1, rascalc_results2, pycorr_files1, pycorr_files2, output_cov_file, n_mu_bins, r_step, skip_r_bins, output_cov_file1, output_cov_file2) \ No newline at end of file +combine_covs_multi(args.rascalc_files[0], args.rascalc_files[1], [args.pycorr_files[i] for i in range(0, 6, 2)], [args.pycorr_files[i] for i in range(1, 6, 2)], args.output_cov_file, args.n_mu_bins or None, args.r_step, args.skip_r_bins, args.output_cov_file1, args.output_cov_file2) \ No newline at end of file diff --git a/scripts/compute_correction_function.py b/scripts/compute_correction_function.py index 77a832a3..bd8ffa21 100644 --- a/scripts/compute_correction_function.py +++ b/scripts/compute_correction_function.py @@ -1,23 +1,20 @@ -### Function to fit a model to the survey correction function, defined as the ratio between model and true RR pair counts for a single survey. This fits a piecewise polynomial model to the data. +### Script to fit a model to the survey correction function, defined as the ratio between model and true RR pair counts for a single survey. This fits a piecewise polynomial model to the data. ## NB: Input RR counts should NOT be normalized by summed galaxy weights here. ## NB: Assume mu is in [0,1] limit here -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (5, 6): - print("Usage: python compute_correction_function.py {RANDOM_PARTICLE_FILE} {BIN_FILE} {OUTPUT_DIR} {PERIODIC} [{RR_COUNTS}] ") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to fit a model to the survey correction function, defined as the ratio between model and true RR pair counts for a single survey.") +parser.add_argument("random_particle_file", type=str, help="file containing the random points' coordinates and weights") +parser.add_argument("r_bin_file", type=str, help="file containing the radial/separation bin boundaries in its rows") +parser.add_argument("output_dir", type=str, help="directory to write the resulting correction function coefficients") +parser.add_argument("periodic", type=bool, help="periodic boundary conditions flag") +parser.add_argument("RR_counts_file", type=str, default=None, nargs='?', help="file containing the RR counts in radial/separation and angular bins. necessary with aperiodic (realistic survey) geometry. not used with periodic geometry.") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.correction_function import compute_correction_function_from_file -random_file = str(sys.argv[1]) -binfile = str(sys.argv[2]) -outdir = str(sys.argv[3]) -periodic = int(sys.argv[4]) -RR_file = get_arg_safe(5) - -compute_correction_function_from_file(random_file, binfile, outdir, periodic, RR_file) \ No newline at end of file +compute_correction_function_from_file(args.random_particle_file, args.r_bin_file, args.output_dir, args.periodic, args.RR_counts_file) \ No newline at end of file diff --git a/scripts/compute_correction_function_multi.py b/scripts/compute_correction_function_multi.py index c6e1d86f..5cc5da81 100644 --- a/scripts/compute_correction_function_multi.py +++ b/scripts/compute_correction_function_multi.py @@ -1,26 +1,20 @@ -### Function to fit a model to the survey correction function, defined as the ratio between model and true RR pair counts for a two surveys. This fits a piecewise polynomial model to the data. +### Script to fit a model to the survey correction function, defined as the ratio between model and true RR pair counts for a two-tracer survey. This fits a piecewise polynomial model to the data. ## NB: Input RR counts should NOT be normalized by summed galaxy weights here. ## NB: Assume mu is in [0,1] limit here -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (6, 9): - print("Usage: python compute_correction_function_multi.py {RANDOM_PARTICLE_FILE_1} {RANDOM_PARTICLE_FILE_2} {BIN_FILE} {OUTPUT_DIR} {PERIODIC} [{RR_COUNTS_11} {RR_COUNTS_12} {RR_COUNTS_22}]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to fit a model to the survey correction function, defined as the ratio between model and true RR pair counts for a two-tracer survey.") +parser.add_argument("random_particle_files", type=str, nargs=2, help="2 files containing the random points' coordinates and weights (one file for each tracer)") +parser.add_argument("r_bin_file", type=str, help="file containing the radial/separation bin boundaries in its rows") +parser.add_argument("output_dir", type=str, help="directory to write the resulting correction function coefficients") +parser.add_argument("periodic", type=bool, help="periodic boundary conditions flag") +parser.add_argument("RR_counts_files", type=str, default=[None]*3, nargs='*', help="3 files containing the RR counts in radial/separation and angular bins: first auto-counts for the first tracer, then the cross-counts between the tracers, and finally auto-counts for the second tracer. necessary with aperiodic (realistic survey) geometry. not used with periodic geometry.") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.correction_function import compute_correction_function_multi_from_files -random_file = str(sys.argv[1]) -random_file2 = str(sys.argv[2]) -binfile = str(sys.argv[3]) -outdir = str(sys.argv[4]) -periodic = int(sys.argv[5]) -RR_file = get_arg_safe(6) -RR_file12 = get_arg_safe(7) -RR_file2 = get_arg_safe(8) - -compute_correction_function_multi_from_files(random_file, random_file2, binfile, outdir, periodic, RR_file, RR_file12, RR_file2) \ No newline at end of file +compute_correction_function_multi_from_files(args.random_particle_files[0], args.random_particle_files[1], args.r_bin_file, args.output_dir, args.periodic, args.RR_counts_files[0], args.RR_counts_files[1], args.RR_counts_files[2]) \ No newline at end of file diff --git a/scripts/convergence_check_extra.py b/scripts/convergence_check_extra.py index 5eabc663..23142da6 100644 --- a/scripts/convergence_check_extra.py +++ b/scripts/convergence_check_extra.py @@ -2,18 +2,15 @@ ## More specifically, divide integral subsamples into halves and check similarity of their average results ## Should work in any case - default, jackknife, Legendre, multi-tracer - as it utilizes universal data from RascalC file -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (2, 3): - print("Usage: python convergence_check_extra.py {RASCALC_RESULTS_FILE} [{N_SUBSAMPLES_TO_USE}]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to perform an extra convergence check on full integrals; more specifically, divide integral subsamples into halves and check similarity of their average results. Should work in any case - default, jackknife, Legendre, multi-tracer - as it utilizes universal data from a RascalC .npy file.") +parser.add_argument("rascalc_results", type=str, help="RascalC .npy filename") +parser.add_argument("n_subsamples", type=int, default=None, nargs='?', help="number of covariance subsamples to use (optional; by default all are used)") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.convergence_check_extra import convergence_check_extra_file -rascalc_results = str(sys.argv[1]) -n_samples = get_arg_safe(2, int) - -convergence_check_extra_file(rascalc_results, n_samples, print_function = print) \ No newline at end of file +convergence_check_extra_file(args.rascalc_results, args.n_subsamples, print_function = print) \ No newline at end of file diff --git a/scripts/convert_counts_from_pycorr.py b/scripts/convert_counts_from_pycorr.py index c5fb1f21..03ddd58e 100644 --- a/scripts/convert_counts_from_pycorr.py +++ b/scripts/convert_counts_from_pycorr.py @@ -1,25 +1,21 @@ "This reads a cosmodesi/pycorr .npy file and generates binned pair counts text file for RascalC to use" -import sys import numpy as np +import argparse -## PARAMETERS -if len(sys.argv) not in (3, 4, 5, 6, 7, 8): - print("Usage: python convert_counts_from_pycorr.py {INPUT_NPY_FILE} {OUTPUT_PAIRCOUNTS_TEXT_FILE} [{R_STEP} [{N_MU} [{COUNTS_FACTOR} [{SPLIT_ABOVE}] [{R_MAX}]]]]]. COUNTS_FACTOR=0 is a special value to use normalized counts.") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads cosmodesi/pycorr .npy file(s) and generates a binned pair counts text file for RascalC to use") +parser.add_argument("pycorr_file", type=str, help="pycorr .npy filename") +parser.add_argument("output_file", type=str, help="output text file for the binned RR pair counts") +parser.add_argument("r_step", type=float, default=1, nargs='?', help="desired width of radial bins (default 1)") +parser.add_argument("n_mu_bins", type=int, default=0, nargs='?', help="number of angular (mu) bins (default 0, which means to preserve the original number)") +parser.add_argument("counts_factor", type=float, default=0, nargs='?', help="basically number of randoms used for these counts, used to convert from total to 1 catalog count estimate. 0 is a special value to use normalized counts, which is easiest and also the default") +parser.add_argument("split_above", type=float, default=np.inf, nargs='?', help="divide weighted RR counts by counts_factor squared below this and by counts_factor above. default is infinity") +parser.add_argument("r_max", type=float, default=np.inf, nargs='?', help="maximum radius (cutoff), the default is infinity (no cutoff)") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.pycorr_utils.counts import convert_counts_from_pycorr_files -infile_name = str(sys.argv[1]) -outfile_name = str(sys.argv[2]) -r_step = get_arg_safe(3, float, 1) -n_mu = get_arg_safe(4, int, None) -counts_factor = get_arg_safe(5, float, None) # basically number of randoms used for these counts, used to convert from total to 1 catalog count estimate -split_above = get_arg_safe(6, float, np.inf) # divide weighted RR counts by counts_factor**2 below this and by counts_factor above -r_max = get_arg_safe(7, float, np.inf) # if given, limit used r_bins at that - -if n_mu == 0: n_mu = None -if counts_factor == 0: counts_factor = None -convert_counts_from_pycorr_files(infile_name, outfile_name, n_mu, r_step, r_max, counts_factor, split_above) \ No newline at end of file +# x or y should be a shortcut for x if x else y (or more explicitly x if x != 0 else y) +convert_counts_from_pycorr_files(args.pycorr_file, args.output_file, args.n_mu_bins or None, args.r_step, args.r_max, args.counts_factor or None, args.split_above) \ No newline at end of file diff --git a/scripts/convert_cov.py b/scripts/convert_cov.py index 43ba3186..41147980 100644 --- a/scripts/convert_cov.py +++ b/scripts/convert_cov.py @@ -1,17 +1,14 @@ -"Simple convenience script that reads RascalC results and saves full cov to text file, in addition checking eigenvalues of bias matrix" +"Simple convenience script that reads RascalC results and saves full cov to text file, in addition checking eigenvalues of the bias matrix" -import sys +import argparse -## PARAMETERS -if len(sys.argv) != 3: - print("Usage: python convert_cov.py {RASCALC_RESULTS_FILE} {OUTPUT_COV_FILE}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="Simple convenience script that reads RascalC results and saves full cov to text file, in addition checking eigenvalues of the bias matrix") +parser.add_argument("rascalc_file", type=str, help="RascalC post-processed .npy filename") +parser.add_argument("output_cov_file", type=str, help="output text file for the covariance") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.cov_utils import export_cov -rascalc_results = str(sys.argv[1]) -output_cov_file = str(sys.argv[2]) - -export_cov(rascalc_results, output_cov_file) \ No newline at end of file +export_cov(args.rascalc_file, args.output_cov_file) \ No newline at end of file diff --git a/scripts/convert_cov_legendre.py b/scripts/convert_cov_legendre.py index 4a1e440e..3509a8a4 100644 --- a/scripts/convert_cov_legendre.py +++ b/scripts/convert_cov_legendre.py @@ -1,18 +1,15 @@ "Slightly more complicated (than convert_cov.py) convenience script that reads RascalC Legendre mode results and saves full cov to text file, changing the indexing from [r, l] to [l, r], in addition checking eigenvalues of bias matrix" -import sys +import argparse -## PARAMETERS -if len(sys.argv) != 4: - print("Usage: python convert_cov_legendre.py {RASCALC_RESULTS_FILE} {MAX_L} {OUTPUT_COV_FILE}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="Simple convenience script that reads RascalC results and saves full cov to text file, changing the single-tracer Legendre-mode indexing from [r, l] to [l, r], in addition checking eigenvalues of bias matrix") +parser.add_argument("rascalc_file", type=str, help="RascalC post-processed .npy filename") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("output_cov_file", type=str, help="output text file for the covariance") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.cov_utils import export_cov_legendre -rascalc_results = str(sys.argv[1]) -max_l = int(sys.argv[2]) -output_cov_file = str(sys.argv[3]) - -export_cov_legendre(rascalc_results, max_l, output_cov_file) +export_cov_legendre(args.rascalc_file, args.max_l, args.output_cov_file) diff --git a/scripts/convert_cov_legendre_multi.py b/scripts/convert_cov_legendre_multi.py index f6970353..623b330c 100644 --- a/scripts/convert_cov_legendre_multi.py +++ b/scripts/convert_cov_legendre_multi.py @@ -1,18 +1,15 @@ "Slightly more complicated (than convert_cov.py) convenience script that reads RascalC 2-tracer Legendre mode results and saves full cov to text file, changing the indexing from [t, r, l] to [t, l, r], in addition checking eigenvalues of bias matrix" -import sys +import argparse -## PARAMETERS -if len(sys.argv) != 4: - print("Usage: python convert_cov_legendre_multi.py {RASCALC_RESULTS_FILE} {MAX_L} {OUTPUT_COV_FILE}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="Simple convenience script that reads RascalC results and saves full cov to text file, changing the 2-tracer Legendre-mode indexing from [t, r, l] to [t, l, r], in addition checking eigenvalues of bias matrix") +parser.add_argument("rascalc_file", type=str, help="RascalC post-processed .npy filename") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("output_cov_file", type=str, help="output text file for the covariance") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.cov_utils import export_cov_legendre_multi -rascalc_results = str(sys.argv[1]) -max_l = int(sys.argv[2]) -output_cov_file = str(sys.argv[3]) - -export_cov_legendre_multi(rascalc_results, max_l, output_cov_file) +export_cov_legendre_multi(args.rascalc_file, args.max_l, args.output_cov_file) diff --git a/scripts/convert_cov_legendre_multi_to_cat.py b/scripts/convert_cov_legendre_multi_to_cat.py index 18358562..342b2da3 100644 --- a/scripts/convert_cov_legendre_multi_to_cat.py +++ b/scripts/convert_cov_legendre_multi_to_cat.py @@ -1,24 +1,20 @@ "This reads a .npy file of RascalC Legendre results (or txt covariance converted previously) and a triplet of cosmodesi/pycorr .npy files to produce a covariance for a catalog of these two tracers concatenated." -import sys +import argparse -## PARAMETERS -if len(sys.argv) not in (9, 11): - print("Usage: python combine_cov_multi_to_cat.py {RASCALC_RESULTS} {PYCORR_FILE_11} {PYCORR_FILE_12} {PYCORR_FILE_22} {R_STEP} {MAX_L} {SKIP_R_BINS} {OUTPUT_COV_FILE} [{BIAS1} {BIAS2}].") - sys.exit(1) +parser = argparse.ArgumentParser(description="This reads a .npy file of RascalC Legendre results (or txt covariance converted previously) and a triplet of cosmodesi/pycorr .npy files to produce a covariance for a catalog of these two tracers concatenated (each tracer upweighted by its linear bias).") +parser.add_argument("rascalc_file", type=str, help="RascalC post-processed .npy filename") +parser.add_argument("pycorr_files", type=str, nargs=3, help="3 pycorr .npy filenames for the distinct correlation functions in the following order: first tracer auto-counts, then cross-counts between the two tracers, and finally the second tracer auto-counts") +parser.add_argument("r_step", type=float, help="desired width of radial bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("skip_r_bins", type=int, help="number of last radial/separation bins to discard") +parser.add_argument("output_cov_file", type=str, help="output text file for the combined tracer covariance") +parser.add_argument("bias1", type=float, default=1, nargs='?', help="linear bias for the first tracer (default 1)") +parser.add_argument("bias2", type=float, default=1, nargs='?', help="linear bias for the second tracer (default 1)") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.comb.convert_cov_legendre_multi_to_cat import convert_cov_legendre_multi_to_cat -rascalc_results = str(sys.argv[1]) -pycorr_files = [str(sys.argv[i]) for i in range(2, 5)] -r_step = float(sys.argv[5]) -max_l = int(sys.argv[6]) -skip_r_bins = int(sys.argv[7]) -output_cov_file = str(sys.argv[8]) -from .utils import get_arg_safe -bias1 = get_arg_safe(9, float, 1) -bias2 = get_arg_safe(10, float, 1) - -convert_cov_legendre_multi_to_cat(rascalc_results, pycorr_files, output_cov_file, max_l, r_step, skip_r_bins, bias1, bias2) \ No newline at end of file +convert_cov_legendre_multi_to_cat(args.rascalc_file, args.pycorr_files, args.output_cov_file, args.max_l, args.r_step, args.skip_r_bins, args.bias1, args.bias2) \ No newline at end of file diff --git a/scripts/convert_to_xyz.py b/scripts/convert_to_xyz.py index 51469b47..357804a1 100644 --- a/scripts/convert_to_xyz.py +++ b/scripts/convert_to_xyz.py @@ -15,30 +15,24 @@ """ -import sys - -if len(sys.argv) not in (3, 4, 5, 6, 7, 8, 9): - print("Usage: python convert_to_xyz.py {INFILE} {OUTFILE} [{OMEGA_M} {OMEGA_K} {W_DARK_ENERGY} [{USE_FKP_WEIGHTS or P0,NZ_name} [{MASK} [{USE_WEIGHTS}]]]]") - sys.exit(1) - -from utils import adjust_path, get_arg_safe +import argparse + +parser = argparse.ArgumentParser(description="Convenience script to convert an input (Ra,Dec,w) FITS or txt file to comoving (x,y,z) coordinates saved as a text file for use with the main C++ code. Output file format has (x,y,z,w) coordinates in Mpc/h units") +parser.add_argument("input_file", type=str, help="input text or FITS file name (text must be rdzw, i.e. RA, DEC, redshift and weight columns)") +parser.add_argument("output_file", type=str, help="output text file name (xyzw, i.e. 3D Cartesian coordinates and weight columns)") +parser.add_argument("Omega_m", type=float, default=0.31, nargs="?", help="relative matter density (default 0.31)") +parser.add_argument("Omega_K", type=float, default=0, nargs="?", help="relative curvature density (default 0)") +parser.add_argument("w_DE", type=float, default=-1, nargs="?", help="dark energy equation of state parameter (default -1)") +parser.add_argument("FKP_weights", type=str, default="0", nargs="?", help="whether to use FKP weights column (default False/0). also accepts format P0,NZ_name. only applies to (DESI) FITS files") +parser.add_argument("mask", type=int, default=0, nargs="?", help="sets bins that all must be set in STATUS for the particle to be selected (default 0). only applies to (DESI) FITS files") +parser.add_argument("use_weights", type=bool, default=True, nargs="?", help="whether to use WEIGHTS column. if not, set unit weights. default True/1") +args = parser.parse_args() + +from utils import adjust_path adjust_path() from RascalC.pre_process.convert_to_xyz import convert_to_xyz_files -from RascalC.pre_process.utils import parse_FKP_arg, my_str_to_bool - -# Load file names -input_file = str(sys.argv[1]) -output_file = str(sys.argv[2]) - -# Read in optional cosmology parameters -Omega_m = get_arg_safe(3, float, 0.31) -Omega_k = get_arg_safe(4, float, 0) -w_dark_energy = get_arg_safe(5, float, -1) -# defaults from the BOSS DR12 2016 clustering paper assuming LCDM +from RascalC.pre_process.utils import parse_FKP_arg -# The next only applies to (DESI) FITS files -FKP_weights = parse_FKP_arg(get_arg_safe(6, str, "False")) # determine whether to use FKP weights -mask = get_arg_safe(7, int, 0) # default is 0 - no mask filtering -use_weights = my_str_to_bool(get_arg_safe(8, str, "True")) # use weights by default +FKP_weights = parse_FKP_arg(args.FKP_weights) # FKP weights needs to be parsed additionally for backward compatibility (could be redesigned better) -convert_to_xyz_files(input_file, output_file, Omega_m, Omega_k, w_dark_energy, FKP_weights, mask, use_weights) \ No newline at end of file +convert_to_xyz_files(args.input_file, args.output_file, args.Omega_m, args.Omega_K, args.w_DE, FKP_weights, args.mask, args.use_weights) \ No newline at end of file diff --git a/scripts/convert_xi_from_pycorr.py b/scripts/convert_xi_from_pycorr.py index 3f64a4c5..2214a340 100644 --- a/scripts/convert_xi_from_pycorr.py +++ b/scripts/convert_xi_from_pycorr.py @@ -1,19 +1,16 @@ "This reads cosmodesi/pycorr .npy file(s) and generates input xi text file for RascalC to use" -import sys +import argparse -## PARAMETERS -if len(sys.argv) < 5: - print("Usage: python convert_xi_from_pycorr.py {INPUT_NPY_FILE1} [{INPUT_NPY_FILE2} ...] {OUTPUT_XI_DAT_FILE} {R_STEP} {N_MU}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads cosmodesi/pycorr .npy file(s) and generates input xi text file for RascalC to use") +parser.add_argument("pycorr_files", type=str, nargs='+', help="pycorr .npy filename(s). if multiple are provided, the correlation function will be averaged") +parser.add_argument("output_xi_file", type=str, help="output correlation function (xi) text file") +parser.add_argument("r_step", type=float, help="desired width of radial bins") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.pycorr_utils.input_xi import convert_xi_from_pycorr_files -infile_names = sys.argv[1:-3] -outfile_name = str(sys.argv[-3]) -r_step = float(sys.argv[-2]) -n_mu = int(sys.argv[-1]) - -convert_xi_from_pycorr_files(infile_names, outfile_name, n_mu, r_step) \ No newline at end of file +convert_xi_from_pycorr_files(args.pycorr_files, args.outfile_name, args.n_mu_bins, args.r_step) \ No newline at end of file diff --git a/scripts/convert_xi_jack_from_pycorr.py b/scripts/convert_xi_jack_from_pycorr.py index 1dfa46d3..713d8fcf 100644 --- a/scripts/convert_xi_jack_from_pycorr.py +++ b/scripts/convert_xi_jack_from_pycorr.py @@ -1,26 +1,24 @@ "This reads a cosmodesi/pycorr .npy file and generates jackknife xi, weight, paircounts, and total paircounts text files for RascalC to use" import numpy as np -import sys +import argparse -## PARAMETERS -if len(sys.argv) not in (6, 7, 8, 9, 10, 11): - print("Usage: python convert_xi_jack_from_pycorr.py {INPUT_NPY_FILE} {OUTPUT_XI_JACK_FILE} {OUTPUT_JACKKNIFE_WEIGHTS_FILE} {OUTPUT_JACKKNIFE_PAIRCOUNTS_FILE} {OUTPUT_BINNED_PAIRCOUNTS_FILE} [{R_STEP} [{N_MU} [{COUNTS_FACTOR} [{SPLIT_ABOVE} [{R_MAX}]]]]]. COUNTS_FACTOR=0 is a special value to use normalized counts.") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads cosmodesi/pycorr .npy file(s) and generates a binned pair counts text file for RascalC to use") +parser.add_argument("pycorr_file", type=str, help="pycorr .npy filename") +parser.add_argument("output_xi_jack_file", type=str, help="output text file for the jackknife 2PCF (xi)") +parser.add_argument("output_jack_weights_file", type=str, help="output text file for the jackknife region weights") +parser.add_argument("output_jack_RR_file", type=str, help="output text file for the jackknife binned RR pair counts") +parser.add_argument("output_RR_file", type=str, help="output text file for the total binned RR pair counts") +parser.add_argument("r_step", type=float, default=1, nargs='?', help="desired width of radial bins (default 1)") +parser.add_argument("n_mu_bins", type=int, default=0, nargs='?', help="number of angular (mu) bins (default 0, which means to preserve the original number)") +parser.add_argument("counts_factor", type=float, default=0, nargs='?', help="basically number of randoms used for these counts, used to convert from total to 1 catalog count estimate. 0 is a special value to use normalized counts, which is easiest and also the default") +parser.add_argument("split_above", type=float, default=np.inf, nargs='?', help="divide weighted RR counts by counts_factor squared below this and by counts_factor above. default is infinity") +parser.add_argument("r_max", type=float, default=np.inf, nargs='?', help="maximum radius (cutoff), the default is infinity (no cutoff)") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.pycorr_utils.jack import convert_jack_xi_weights_counts_from_pycorr_files -infile_name = str(sys.argv[1]) -xi_jack_name = str(sys.argv[2]) -jackweights_name = str(sys.argv[3]) -jackpairs_name = str(sys.argv[4]) -binpairs_name = str(sys.argv[5]) -r_step = get_arg_safe(6, float, 1) -n_mu = get_arg_safe(7, int, None) -counts_factor = get_arg_safe(8, float, None) # basically number of randoms used for these counts, used to convert from total to 1 catalog count estimate -split_above = get_arg_safe(9, float, np.inf) # divide weighted RR counts by counts_factor**2 below this and by counts_factor above -r_max = get_arg_safe(10, float, np.inf) # if given, limit used r_bins at that - -convert_jack_xi_weights_counts_from_pycorr_files(infile_name, xi_jack_name, jackweights_name, jackpairs_name, binpairs_name, n_mu, r_step, r_max, counts_factor, split_above) \ No newline at end of file +# x or y should be a shortcut for x if x else y (or more explicitly x if x != 0 else y) +convert_jack_xi_weights_counts_from_pycorr_files(args.pycorr_file, args.output_xi_jack_file, args.output_jack_weights_file, args.output_jack_RR_file, args.output_RR_file, args.n_mu_bins or None, args.r_step, args.r_max, args.counts_factor or None, args.split_above) \ No newline at end of file diff --git a/scripts/convert_xi_to_multipoles.py b/scripts/convert_xi_to_multipoles.py index 95a01590..58cde7d1 100644 --- a/scripts/convert_xi_to_multipoles.py +++ b/scripts/convert_xi_to_multipoles.py @@ -2,19 +2,16 @@ ### This computes all even multipoles up to a specified maximum ell, approximating the integral by a sum. ### The output form is a text file with the first column specifying the r-bin, the second giving xi_0(r), the third with xi_2(r) etc. -import sys +import argparse -## PARAMETERS -if len(sys.argv) != 4: - print("Usage: python convert_xi_to_multipoles.py {INFILE} {MAX_L} {OUTFILE}") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to convert a measured 2PCF in xi(r,mu) format to Legendre multipoles, i.e. xi_ell(r). This computes all even multipoles up to a specified maximum ell, approximating the integral by a sum. The output form is a text file with the first column specifying the r-bin, the second giving xi_0(r), the third with xi_2(r) etc.") +parser.add_argument("input_file", type=str, help="input xi text file") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("output_file", type=str, help="output xi text file") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.xi.convert_to_multipoles import convert_xi_to_multipoles_files -infile = str(sys.argv[1]) -max_l = int(sys.argv[2]) -outfile = str(sys.argv[3]) - -convert_xi_to_multipoles_files(infile, max_l, outfile) \ No newline at end of file +convert_xi_to_multipoles_files(args.input_file, args.max_l, args.output_file) \ No newline at end of file diff --git a/scripts/create_jackknives_pycorr.py b/scripts/create_jackknives_pycorr.py index 87df1672..3b3c355f 100644 --- a/scripts/create_jackknives_pycorr.py +++ b/scripts/create_jackknives_pycorr.py @@ -1,19 +1,16 @@ -""" This function will assign a jackknife region to each input particle, in a way compatible with pycorr used in DESI data processing. Data is saved as a 5 column text file.""" +""" This script will assign a jackknife region to each input particle, in a way compatible with pycorr used in DESI data processing. Data is saved as a 5 column text file.""" -import sys +import argparse -## PARAMETERS -if len(sys.argv) != 5: - print("Usage: python create_jackknives_pycorr.py {REF_RDZ_FILE} {INPUT_XYZW_FILE} {OUTPUT_XYZWJ_FILE} {NJACK}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script will assign a jackknife region to each input particle, in a way compatible with pycorr used in DESI data processing. Data is saved as a 5 column text file.") +parser.add_argument("reference_file", type=str, help="reference (galaxies/data) text file name (rdz, i.e. RA, DEC, redshift columns; can have additional columns but they will not be used)") +parser.add_argument("input_file", type=str, help="input (randoms) text file name (xyzw, i.e. 3D Cartesian coordinates and weight columns)") +parser.add_argument("output_file", type=str, help="output (randoms) text file name (xyzwj, i.e. 3D Cartesian coordinates, weight and jackknife region number columns)") +parser.add_argument("n_jack", type=int, help="the number of jackknife regions to create") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.pre_process.create_jackknives_pycorr import create_jackknives_pycorr_files -reffile_name = str(sys.argv[1]) -infile_name = str(sys.argv[2]) -outfile_name = str(sys.argv[3]) -njack = int(sys.argv[4]) - -create_jackknives_pycorr_files(reffile_name, infile_name, outfile_name, njack) \ No newline at end of file +create_jackknives_pycorr_files(args.reference_file, args.input_file, args.output_file, args.n_jack) \ No newline at end of file diff --git a/scripts/mu_bin_legendre_factors.py b/scripts/mu_bin_legendre_factors.py index 89e28abb..a610192d 100644 --- a/scripts/mu_bin_legendre_factors.py +++ b/scripts/mu_bin_legendre_factors.py @@ -1,18 +1,15 @@ "Simple script to produce the file of mu bin Legendre factors for the C++ code in LEGENDRE_MIX mode." -import sys +import argparse -## PARAMETERS -if len(sys.argv) != 4: - print("Usage: python mu_bin_legendre_factors.py {N_MU_BINS} {MAX_L} {OUTPUT_DIR}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="Simple script to produce the file of mu bin Legendre factors for the C++ code in LEGENDRE_MIX mode.") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("output_dir", type=str, help="directory to save the resulting file") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.mu_bin_legendre_factors import write_mu_bin_legendre_factors -n_mu_bins = int(sys.argv[1]) -max_l = int(sys.argv[2]) -output_dir = str(sys.argv[3]) - -write_mu_bin_legendre_factors(n_mu_bins, max_l, output_dir) \ No newline at end of file +write_mu_bin_legendre_factors(args.n_mu_bins, args.max_l, args.output_dir) \ No newline at end of file diff --git a/scripts/normalize_weights.py b/scripts/normalize_weights.py index e45a4a80..4bb463bc 100644 --- a/scripts/normalize_weights.py +++ b/scripts/normalize_weights.py @@ -6,19 +6,15 @@ """ -import sys +import argparse -# Check number of parameters -if len(sys.argv) != 3: - print("Usage: python normalize_weights.py {INFILE} {OUTFILE}") - sys.exit(1) +parser = argparse.ArgumentParser(description="Convenience script to normalize weights in an input (x,y,z,w) or (x,y,z,w,j) text file and save the result in another text file for use with further scripts.") +parser.add_argument("input_file", type=str) +parser.add_argument("output_file", type=str) +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.pre_process.normalize_weights import normalize_weights_files -# Load file names -input_file = str(sys.argv[1]) -output_file = str(sys.argv[2]) - -normalize_weights_files(input_file, output_file) \ No newline at end of file +normalize_weights_files(args.input_file, args.output_file) \ No newline at end of file diff --git a/scripts/post_process_default.py b/scripts/post_process_default.py index 7a210d4a..ef2928c3 100644 --- a/scripts/post_process_default.py +++ b/scripts/post_process_default.py @@ -1,22 +1,19 @@ ## Script to post-process the single-field integrals computed by the C++ code. ## We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (5, 6, 7): - print("Usage: python post_process_default.py {COVARIANCE_DIR} {N_R_BINS} {N_MU_BINS} {OUTPUT_DIR} [{SHOT_NOISE_RESCALING}]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the single-field integrals computed by the C++ code (in radial and angular/mu bins).") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_r_bins", type=int, help="number of radial/separation bins") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("shot_noise_rescaling", type=float, default=1, nargs='?', help="shot-noise rescaling parameter value to use") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_default -file_root = str(sys.argv[1]) -n = int(sys.argv[2]) -m = int(sys.argv[3]) -outdir = str(sys.argv[4]) -alpha = get_arg_safe(5, float, 1) -skip_r_bins = get_arg_safe(6, int, 0) - -post_process_default(file_root, n, m, outdir, alpha, skip_r_bins) \ No newline at end of file +post_process_default(args.covariance_dir, args.n_r_bins, args.n_mu_bins, args.output_dir, args.shot_noise_rescaling, args.skip_r_bins) \ No newline at end of file diff --git a/scripts/post_process_default_mocks.py b/scripts/post_process_default_mocks.py index 407050f3..43dc63fa 100644 --- a/scripts/post_process_default_mocks.py +++ b/scripts/post_process_default_mocks.py @@ -1,22 +1,19 @@ ## Script to post-process the single-field integrals computed by the C++ code. This computes the shot-noise rescaling parameter, alpha, from a mock derived covariance matrix. ## We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (6, 7): - print("Usage: python post_process_default_mocks.py {MOCK_COV_FILE} {COVARIANCE_DIR} {N_R_BINS} {N_MU_BINS} {OUTPUT_DIR} [{SKIP_R_BINS}]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the single-field integrals computed by the C++ code (in radial and angular/mu bins), optimizing the shot-noise rescaling parameter to the mock-derived covariance matrix.") +parser.add_argument("mock_cov_file", type=str, help="name of the text file with the mock sample covariance matrix") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_r_bins", type=int, help="number of radial/separation bins") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_default_mocks - -mock_cov_file = str(sys.argv[1]) -file_root = str(sys.argv[2]) -n = int(sys.argv[3]) -m = int(sys.argv[4]) -outdir = str(sys.argv[5]) -skip_r_bins = get_arg_safe(6, int, 0) -post_process_default_mocks(mock_cov_file, file_root, n, m, outdir, skip_r_bins) +post_process_default_mocks(args.mock_cov_file, args.covariance_dir, args.n_r_bins, args.n_mu_bins, args.output_dir, args.skip_r_bins) diff --git a/scripts/post_process_default_mocks_multi.py b/scripts/post_process_default_mocks_multi.py index 834b34ce..94f660f2 100644 --- a/scripts/post_process_default_mocks_multi.py +++ b/scripts/post_process_default_mocks_multi.py @@ -1,23 +1,19 @@ ## Script to post-process the multi-field integrals computed by the C++ code. This computes two shot-noise rescaling parameters, alphas, from a mock derived covariance matrix. ## We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (7, 8): - print("Usage: python post_process_default_mocks_multi.py {MOCK_COV_FILE} {COVARIANCE_DIR} {N_R_BINS} {N_MU_BINS} {N_SUBSAMPLES} {OUTPUT_DIR} [{SKIP_R_BINS}]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the multi-field integrals computed by the C++ code (in radial and angular/mu bins), optimizing the two shot-noise rescaling parameters to the mock-derived covariance matrix.") +parser.add_argument("mock_cov_file", type=str, help="name of the text file with the mock sample covariance matrix") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_r_bins", type=int, help="number of radial/separation bins") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_default_mocks_multi - -mock_cov_file = str(sys.argv[1]) -file_root = str(sys.argv[2]) -n = int(sys.argv[3]) -m = int(sys.argv[4]) -n_samples = int(sys.argv[5]) -outdir = str(sys.argv[6]) -skip_r_bins = get_arg_safe(7, int, 0) # convert from radial to total number of bins right away -post_process_default_mocks_multi(mock_cov_file, file_root, n, m, n_samples, outdir, skip_r_bins) +post_process_default_mocks_multi(args.mock_cov_file, args.covariance_dir, args.n_r_bins, args.n_mu_bins, args.output_dir, args.skip_r_bins) diff --git a/scripts/post_process_default_multi.py b/scripts/post_process_default_multi.py index 387c5e7e..e38d144b 100644 --- a/scripts/post_process_default_multi.py +++ b/scripts/post_process_default_multi.py @@ -1,23 +1,20 @@ ## Script to post-process the multi-field integrals computed by the C++ code. ## We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (5, 7, 8): - print("Usage: python post_process_default_multi.py {COVARIANCE_DIR} {N_R_BINS} {N_MU_BINS} {OUTPUT_DIR} [{SHOT_NOISE_RESCALING_1} {SHOT_NOISE_RESCALING_2} [{SKIP_R_BINS}]]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the multi-field integrals computed by the C++ code (in radial and angular/mu bins).") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_r_bins", type=int, help="number of radial/separation bins") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("shot_noise_rescaling1", type=float, default=1, nargs='?', help="shot-noise rescaling parameter value to use for the first tracer") +parser.add_argument("shot_noise_rescaling2", type=float, default=1, nargs='?', help="shot-noise rescaling parameter value to use for the second tracer") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_default_multi -file_root = str(sys.argv[1]) -n = int(sys.argv[2]) -m = int(sys.argv[3]) -outdir = str(sys.argv[4]) -alpha_1 = get_arg_safe(5, float, 1) -alpha_2 = get_arg_safe(6, float, 1) -skip_r_bins = get_arg_safe(7, int, 0) - -post_process_default_multi(file_root, n, m, outdir, alpha_1, alpha_2, skip_r_bins) +post_process_default_multi(args.covariance_dir, args.n_r_bins, args.n_mu_bins, args.output_dir, args.shot_noise_rescaling1, args.shot_noise_rescaling2, args.skip_r_bins) diff --git a/scripts/post_process_jackknife.py b/scripts/post_process_jackknife.py index 62fc7d3e..0ba0651f 100644 --- a/scripts/post_process_jackknife.py +++ b/scripts/post_process_jackknife.py @@ -1,22 +1,19 @@ ## Script to post-process the single-field integrals computed by the C++ code. This computes the shot-noise rescaling parameter, alpha, from a data derived covariance matrix. ## We output the data and theory jackknife covariance matrices, in addition to full theory covariance matrices and (quadratic-bias corrected) precision matrices. The effective number of samples, N_eff, is also computed. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (5, 6): - print("Usage: python post_process_jackknife.py {XI_JACKKNIFE_FILE} {WEIGHTS_DIR} {COVARIANCE_DIR} {N_MU_BINS} {OUTPUT_DIR} [{SKIP_R_BINS}]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the single-field integrals computed by the C++ code (in radial and angular/mu bins), optimizing the shot-noise rescaling parameter to a data-derived (jackknife) covariance matrix.") +parser.add_argument("xi_jack_file", type=str, help="name of the text file with the jackknife correlation function estimates (in proper RascalC format)") +parser.add_argument("weight_dir", type=str, help="directory containing the jackknife region weights and RR counts") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_jackknife -jackknife_file = str(sys.argv[1]) -weight_dir = str(sys.argv[2]) -file_root = str(sys.argv[3]) -m = int(sys.argv[4]) -outdir = str(sys.argv[5]) -skip_r_bins = get_arg_safe(6, int, 0) - -post_process_jackknife(jackknife_file, weight_dir, file_root, m, outdir, skip_r_bins) \ No newline at end of file +post_process_jackknife(args.xi_jack_file, args.weight_dir, args.file_root, args.n_mu_bins, args.output_dir, args.skip_r_bins) \ No newline at end of file diff --git a/scripts/post_process_jackknife_multi.py b/scripts/post_process_jackknife_multi.py index 6b710ced..374aed18 100644 --- a/scripts/post_process_jackknife_multi.py +++ b/scripts/post_process_jackknife_multi.py @@ -2,24 +2,19 @@ ## We output the data and theory jackknife covariance matrices, in addition to full theory covariance matrices and (quadratic-bias corrected) precision matrices. ## The effective number of samples, N_eff, is also computed. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (8, 9): - print("Usage: python post_process_jackknife_multi.py {XI_JACKKNIFE_FILE_11} {XI_JACKKNIFE_FILE_12} {XI_JACKKNIFE_FILE_22} {WEIGHTS_DIR} {COVARIANCE_DIR} {N_MU_BINS} {OUTPUT_DIR} [{SKIP_R_BINS}]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the single-field integrals computed by the C++ code (in radial and angular/mu bins), optimizing the shot-noise rescaling parameter to a data-derived (jackknife) covariance matrix.") +parser.add_argument("xi_jack_files", type=str, nargs=3, help="name of the 3 text files with the jackknife correlation function estimates (in proper RascalC format) in the following order: the first tracer auto-correlation, the cross-correlation of the two tracers, the second tracer auto-correlation") +parser.add_argument("weight_dir", type=str, help="directory containing the jackknife region weights and RR counts") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_jackknife_multi -jackknife_file_11 = str(sys.argv[1]) -jackknife_file_12 = str(sys.argv[2]) -jackknife_file_22 = str(sys.argv[3]) -weight_dir = str(sys.argv[4]) -file_root = str(sys.argv[5]) -m = int(sys.argv[6]) -outdir = str(sys.argv[7]) -skip_r_bins = get_arg_safe(8, int, 0) - -post_process_jackknife_multi(jackknife_file_11, jackknife_file_12, jackknife_file_22, weight_dir, file_root, m, outdir, skip_r_bins) +post_process_jackknife_multi(args.xi_jack_files[0], args.xi_jack_files[1], args.xi_jack_files[2], args.weight_dir, args.covariance_dir, args.n_mu_bins, args.output_dir, args.skip_r_bins) diff --git a/scripts/post_process_legendre.py b/scripts/post_process_legendre.py index db3a1d8c..4b6623e7 100644 --- a/scripts/post_process_legendre.py +++ b/scripts/post_process_legendre.py @@ -1,23 +1,20 @@ ## Script to post-process the single-field Legendre binned integrals computed by the C++ code. ## We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (5, 6, 7, 8): - print("Usage: python post_process_legendre.py {COVARIANCE_DIR} {N_R_BINS} {MAX_L} {OUTPUT_DIR} [{SHOT_NOISE_RESCALING} [{SKIP_R_BINS} [{SKIP_L}]]]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the single-field Legendre binned integrals computed by the C++ code") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_r_bins", type=int, help="number of radial/separation bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("shot_noise_rescaling", type=float, default=1, nargs='?', help="shot-noise rescaling parameter value to use") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +parser.add_argument("skip_l", type=int, default=0, nargs='?', help="number of last Legendre multipoles to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_legendre -file_root = str(sys.argv[1]) -n = int(sys.argv[2]) -max_l = int(sys.argv[3]) -outdir = str(sys.argv[4]) -alpha = get_arg_safe(5, float, 1) -skip_r_bins = get_arg_safe(6, int, 0) -skip_l = get_arg_safe(7, int, 0) - -post_process_legendre(file_root, n, max_l, outdir, alpha, skip_r_bins, skip_l) \ No newline at end of file +post_process_legendre(args.covariance_dir, args.n_r_bins, args.max_l, args.output_dir, args.shot_noise_rescaling, args.skip_r_bins, args.skip_l) \ No newline at end of file diff --git a/scripts/post_process_legendre_mix_jackknife.py b/scripts/post_process_legendre_mix_jackknife.py index cc529670..7f9f58c8 100644 --- a/scripts/post_process_legendre_mix_jackknife.py +++ b/scripts/post_process_legendre_mix_jackknife.py @@ -1,24 +1,21 @@ ## Script to post-process the single-field integrals computed by the C++ code in mixed Legendre (LEGENDRE_MIX) mode. This computes the shot-noise rescaling parameter, alpha, from a data derived covariance matrix. ## We output the data and theory jackknife covariance matrices, in addition to full theory covariance matrices and (quadratic-bias corrected) precision matrices. The effective number of samples, N_eff, is also computed. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (7, 8, 9): - print("Usage: python post_process_legendre_mix_jackknife.py {XI_JACKKNIFE_FILE} {WEIGHTS_DIR} {COVARIANCE_DIR} {N_MU_BINS} {MAX_L} {OUTPUT_DIR} [{SKIP_R_BINS} [{SKIP_L}]]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the single-field Legendre binned integrals computed by the C++ code, optimizing the shot-noise rescaling parameter to a data-derived (jackknife) covariance matrix.") +parser.add_argument("xi_jack_file", type=str, help="name of the text file with the jackknife correlation function estimates (in proper RascalC format)") +parser.add_argument("weight_dir", type=str, help="directory containing the jackknife region weights and RR counts") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +parser.add_argument("skip_l", type=int, default=0, nargs='?', help="number of last Legendre multipoles to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_legendre_mix_jackknife - -jackknife_file = str(sys.argv[1]) -weight_dir = str(sys.argv[2]) -file_root = str(sys.argv[3]) -m = int(sys.argv[4]) -max_l = int(sys.argv[5]) -outdir = str(sys.argv[6]) -skip_r_bins = get_arg_safe(7, int, 0) -skip_l = get_arg_safe(8, int, 0) -post_process_legendre_mix_jackknife(jackknife_file, weight_dir, file_root, m, max_l, outdir, skip_r_bins, skip_l) \ No newline at end of file +post_process_legendre_mix_jackknife(args.xi_jack_file, args.weight_dir, args.file_root, args.n_mu_bins, args.max_l, args.output_dir, args.skip_r_bins, args.skip_l) \ No newline at end of file diff --git a/scripts/post_process_legendre_mocks.py b/scripts/post_process_legendre_mocks.py index d6fa1f4d..70aaf74c 100644 --- a/scripts/post_process_legendre_mocks.py +++ b/scripts/post_process_legendre_mocks.py @@ -1,23 +1,20 @@ ## Script to post-process the single-field Legendre binned integrals computed by the C++ code. This computes the shot-noise rescaling parameter, alpha, from a mock derived covariance matrix. ## We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (6, 7, 8): - print("Usage: python post_process_legendre_mocks.py {MOCK_COV_FILE} {COVARIANCE_DIR} {N_R_BINS} {MAX_L} {OUTPUT_DIR} [{SKIP_R_BINS} [{SKIP_L}]]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the single-field Legendre binned integrals computed by the C++ code, optimizing the shot-noise rescaling parameter to the mock-derived covariance matrix.") +parser.add_argument("mock_cov_file", type=str, help="name of the text file with the mock sample covariance matrix") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_r_bins", type=int, help="number of radial/separation bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +parser.add_argument("skip_l", type=int, default=0, nargs='?', help="number of last Legendre multipoles to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_legendre_mocks - -mock_cov_file = str(sys.argv[1]) -file_root = str(sys.argv[2]) -n = int(sys.argv[3]) -max_l = int(sys.argv[4]) -outdir = str(sys.argv[5]) -skip_r_bins = get_arg_safe(6, int, 0) -skip_l = get_arg_safe(7, int, 0) -post_process_legendre_mocks(mock_cov_file, file_root, n, max_l, outdir, skip_r_bins, skip_l) \ No newline at end of file +post_process_legendre_mocks(args.mock_cov_file, args.covariance_dir, args.n_r_bins, args.max_l, args.output_dir, args.skip_r_bins, args.skip_l) \ No newline at end of file diff --git a/scripts/post_process_legendre_multi.py b/scripts/post_process_legendre_multi.py index dc3adfaa..0d1fd404 100644 --- a/scripts/post_process_legendre_multi.py +++ b/scripts/post_process_legendre_multi.py @@ -1,24 +1,21 @@ ## Script to post-process the multi-field Legendre binned integrals computed by the C++ code, given a shot-noise rescaling parameter alpha. ## We output the theoretical covariance matrices, (quadratic-bias corrected) precision matrices and the effective number of samples, N_eff. -import sys +import argparse -# PARAMETERS -if len(sys.argv) not in (5, 7, 8, 9): - print("Usage: python post_process_legendre_multi.py {COVARIANCE_DIR} {N_R_BINS} {MAX_L} {N_SUBSAMPLES} {OUTPUT_DIR} [{SHOT_NOISE_RESCALING_1} {SHOT_NOISE_RESCALING_2}]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to post-process the multi-field Legendre binned integrals computed by the C++ code") +parser.add_argument("covariance_dir", type=str, help="directory containing the covariance matrix subdirectory") +parser.add_argument("n_r_bins", type=int, help="number of radial/separation bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("output_dir", type=str, help="directory to write the post-processing results") +parser.add_argument("shot_noise_rescaling1", type=float, default=1, nargs='?', help="shot-noise rescaling parameter value to use for the first tracer") +parser.add_argument("shot_noise_rescaling2", type=float, default=1, nargs='?', help="shot-noise rescaling parameter value to use for the second tracer") +parser.add_argument("skip_r_bins", type=int, default=0, nargs='?', help="number of last radial/separation bins to discard") +parser.add_argument("skip_l", type=int, default=0, nargs='?', help="number of last Legendre multipoles to discard") +args = parser.parse_args() -from utils import adjust_path, get_arg_safe +from utils import adjust_path adjust_path() from RascalC.post_process import post_process_legendre_multi -file_root = str(sys.argv[1]) -n = int(sys.argv[2]) -max_l = int(sys.argv[3]) -outdir = str(sys.argv[4]) -alpha_1 = get_arg_safe(5, float, 1) -alpha_2 = get_arg_safe(6, float, 1) -skip_r_bins = get_arg_safe(7, int, 0) -skip_l = get_arg_safe(8, int, 0) - -post_process_legendre_multi(file_root, n, max_l, outdir, alpha_1, alpha_2, skip_r_bins, skip_l) +post_process_legendre_multi(args.covariance_dir, args.n_r_bins, args.max_l, args.output_dir, args.shot_noise_rescaling1, args.shot_noise_rescaling2, args.skip_r_bins, args.skip_l) diff --git a/scripts/print_shot_noise_rescaling.py b/scripts/print_shot_noise_rescaling.py index 49a3017a..32c66b41 100644 --- a/scripts/print_shot_noise_rescaling.py +++ b/scripts/print_shot_noise_rescaling.py @@ -1,14 +1,13 @@ ## Simple convenience script to print shot noise rescaling from one or many RascalC files -import sys +import argparse -# PARAMETERS -if len(sys.argv) < 2: - print("Usage: python print_shot_noise_rescaling.py {RASCALC_RESULTS_1} [{RASCALC_RESULTS_2} ...]") - sys.exit(1) +parser = argparse.ArgumentParser(description="Simple convenience script to print shot noise rescaling from one or many RascalC files") +parser.add_argument("rascalc_results", type=str, nargs="+", help="RascalC .npy filename(s)") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.get_shot_noise_rescaling import get_shot_noise_rescalings -get_shot_noise_rescalings(sys.argv[1:], print_function = print) \ No newline at end of file +get_shot_noise_rescalings(args.rascalc_results, print_function = print) \ No newline at end of file diff --git a/scripts/redshift_cut.py b/scripts/redshift_cut.py index ada1d663..753bee92 100644 --- a/scripts/redshift_cut.py +++ b/scripts/redshift_cut.py @@ -13,30 +13,23 @@ """ -import sys - -# Check number of parameters -if len(sys.argv) not in (5, 6, 7, 8): - print("Usage: python redshift_cut.py {INFILE} {OUTFILE} {Z_MIN} {Z_MAX} [{USE_FKP_WEIGHTS or P0,NZ_name} [{MASK} [{USE_WEIGHTS}]]]") - sys.exit(1) - -from utils import adjust_path, get_arg_safe +import argparse + +parser = argparse.ArgumentParser(description="Convenience script to perform redshift cut on an input (Ra,Dec,z,w) FITS or txt file and save the result in a text file for use with further scripts. Output file format has (Ra,Dec,z,w) coordinates") +parser.add_argument("input_file", type=str, help="input text or FITS file name (text must be rdzw, i.e. RA, DEC, redshift and weight columns)") +parser.add_argument("output_file", type=str, help="output text file name (rdzw, i.e. RA, DEC, redshift and weight columns)") +parser.add_argument("z_min", type=float) +parser.add_argument("z_max", type=float) +parser.add_argument("FKP_weights", type=str, default=0, nargs="?", help="whether to use FKP weights column (default False/0). also accepts format P0,NZ_name. only applies to (DESI) FITS files") +parser.add_argument("mask", type=int, default="0", nargs="?", help="sets bins that all must be set in STATUS for the particle to be selected (default 0). only applies to (DESI) FITS files") +parser.add_argument("use_weights", type=bool, default=True, nargs="?", help="whether to use WEIGHTS column. if not, set unit weights. default True/1") +args = parser.parse_args() + +from utils import adjust_path adjust_path() from RascalC.pre_process.redshift_cut import redshift_cut_files -from RascalC.pre_process.utils import parse_FKP_arg, my_str_to_bool - -# Load file names -input_file = str(sys.argv[1]) -output_file = str(sys.argv[2]) -print("\nUsing input file %s in Ra,Dec,z coordinates\n"%input_file) - -# Load min and max redshifts -z_min = float(sys.argv[3]) -z_max = float(sys.argv[4]) +from RascalC.pre_process.utils import parse_FKP_arg -# The next only applies to (DESI) FITS files -FKP_weights = parse_FKP_arg(get_arg_safe(5, str, "False")) # determine whether to use FKP weights -mask = get_arg_safe(6, int, 0) # default is 0 - no mask filtering -use_weights = my_str_to_bool(get_arg_safe(7, str, "True")) # use weights by default +FKP_weights = parse_FKP_arg(args.FKP_weights) # FKP weights needs to be parsed additionally for backward compatibility (could be redesigned better) -redshift_cut_files(input_file, output_file, z_min, z_max, FKP_weights, mask, use_weights) \ No newline at end of file +redshift_cut_files(args.input_file, args.output_file, args.z_min, args.z_max, FKP_weights, args.mask, args.use_weights) \ No newline at end of file diff --git a/scripts/sample_cov_from_pycorr.py b/scripts/sample_cov_from_pycorr.py index 1fa6ee06..d053220f 100644 --- a/scripts/sample_cov_from_pycorr.py +++ b/scripts/sample_cov_from_pycorr.py @@ -1,27 +1,27 @@ "This reads cosmodesi/pycorr .npy file(s) and generates sample covariance of xi(s,mu) in text format" -import sys +import argparse -## PARAMETERS -if len(sys.argv) < 8: - print("Usage: python sample_cov_from_pycorr.py {INPUT_NPY_FILE1} {INPUT_NPY_FILE2} [{INPUT_NPY_FILE3} ...] {OUTPUT_COV_FILE} {R_STEP} {N_MU} {R_MAX} {N_CORRELATIONS}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads cosmodesi/pycorr .npy file(s) and generates sample covariance of xi(s,mu) in text format") +parser.add_argument("pycorr_files", type=str, nargs='+', help="pycorr .npy filenames. in multi-tracer case, you should give all the realizations for the first tracer's auto-correlation, then the same realizations in the same order for the cross-correlation between tracers 1 and 2, and then for the auto-correlation of the second tracer (this can be continued further, but RascalC can not compute semi-analytic covariances for more than 2 tracers at a time.)") +parser.add_argument("output_file", type=str, help="output text file for the sample covariance matrix") +parser.add_argument("r_step", type=float, help="desired width of radial bins") +parser.add_argument("n_mu_bins", type=int, help="number of angular (mu) bins") +parser.add_argument("r_max", type=float, help="maximum radius (cutoff)") +parser.add_argument("n_correlations", type=int, help="number of correlation functions (1 for single tracer, typically 3 for two tracers)") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.pycorr_utils.sample_cov import sample_cov_from_pycorr_files -infile_names = sys.argv[1:-5] -outfile_name = str(sys.argv[-5]) -r_step = float(sys.argv[-4]) -n_mu = int(sys.argv[-3]) -r_max = float(sys.argv[-2]) -n_corr = int(sys.argv[-1]) +n_corr: int = args.n_correlations +pycorr_files: list[str] = args.pycorr_files assert n_corr >= 1, "Need to have at least one correlation" -n_files = len(infile_names) +n_files = len(pycorr_files) assert n_files % n_corr == 0, "Need to have the same number of files for all correlations" n_samples = n_files // n_corr -infile_names = [infile_names[n_samples * i, n_samples * (i+1)] for i in range(n_corr)] +pycorr_files = [pycorr_files[n_samples * i, n_samples * (i+1)] for i in range(n_corr)] -sample_cov_from_pycorr_files(infile_names, outfile_name, n_mu, r_step, r_max) \ No newline at end of file +sample_cov_from_pycorr_files(pycorr_files, args.output_file, args.n_mu_bins, args.r_step, args.r_max) \ No newline at end of file diff --git a/scripts/sample_cov_multipoles_from_pycorr.py b/scripts/sample_cov_multipoles_from_pycorr.py index c4a8a650..276cc06a 100644 --- a/scripts/sample_cov_multipoles_from_pycorr.py +++ b/scripts/sample_cov_multipoles_from_pycorr.py @@ -1,27 +1,27 @@ "This reads cosmodesi/pycorr .npy file(s) and generates sample covariance of xi_l(s) in text format" -import sys +import argparse -## PARAMETERS -if len(sys.argv) < 8: - print("Usage: python sample_cov_multipoles_from_pycorr.py {INPUT_NPY_FILE1} {INPUT_NPY_FILE2} [{INPUT_NPY_FILE3} ...] {OUTPUT_COV_FILE} {R_STEP} {MAX_L} {R_MAX} {N_CORRELATIONS}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script reads cosmodesi/pycorr .npy file(s) and generates sample covariance of xi_l(s) in text format") +parser.add_argument("pycorr_files", type=str, nargs='+', help="pycorr .npy filenames. in multi-tracer case, you should give all the realizations for the first tracer's auto-correlation, then the same realizations in the same order for the cross-correlation between tracers 1 and 2, and then for the auto-correlation of the second tracer (this can be continued further, but RascalC can not compute semi-analytic covariances for more than 2 tracers at a time.)") +parser.add_argument("output_file", type=str, help="output text file for the sample covariance matrix") +parser.add_argument("r_step", type=float, help="desired width of radial bins") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("r_max", type=float, help="maximum radius (cutoff)") +parser.add_argument("n_corr", type=int, help="number of correlation functions (1 for single tracer, typically 3 for two tracers)") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.pycorr_utils.sample_cov_multipoles import sample_cov_multipoles_from_pycorr_files -infile_names = sys.argv[1:-5] -outfile_name = str(sys.argv[-5]) -r_step = float(sys.argv[-4]) -max_l = int(sys.argv[-3]) -r_max = float(sys.argv[-2]) -n_corr = int(sys.argv[-1]) +n_corr: int = args.n_correlations +pycorr_files: list[str] = args.pycorr_files assert n_corr >= 1, "Need to have at least one correlation" -n_files = len(infile_names) +n_files = len(pycorr_files) assert n_files % n_corr == 0, "Need to have the same number of files for all correlations" n_samples = n_files // n_corr -infile_names = [infile_names[n_samples * i, n_samples * (i+1)] for i in range(n_corr)] +pycorr_files = [pycorr_files[n_samples * i, n_samples * (i+1)] for i in range(n_corr)] -sample_cov_multipoles_from_pycorr_files(infile_names, outfile_name, max_l, r_step, r_max) \ No newline at end of file +sample_cov_multipoles_from_pycorr_files(pycorr_files, args.output_file, args.max_l, args.r_step, args.r_max) \ No newline at end of file diff --git a/scripts/smoothen_xi.py b/scripts/smoothen_xi.py index 77b7a8b6..64e80d67 100644 --- a/scripts/smoothen_xi.py +++ b/scripts/smoothen_xi.py @@ -3,21 +3,18 @@ ### Then r^2*xi_ell(r) are smoothened with Savitzky-Golay filter ### Finally, they are transformed back to xi(r,mu) format. -import sys +import argparse -## PARAMETERS -if len(sys.argv) != 6: - print("Usage: python smoothen_xi.py {INFILE} {MAX_L} {RADIAL_WINDOW_LENGTH} {RADIAL_POLYORDER} {OUTFILE}") - sys.exit(1) +parser = argparse.ArgumentParser(description="Script to smoothen a measured 2PCF in xi(r,mu) format. This computes all even multipoles xi_ell(r) up to a specified maximum ell, approximating the integral by a sum. Then r^2*xi_ell(r) are smoothened with Savitzky-Golay filter. Finally, they are transformed back to xi(r,mu) format.") +parser.add_argument("input_file", type=str, help="input xi text file") +parser.add_argument("max_l", type=int, help="maximum multipole index") +parser.add_argument("radial_window_length", type=int, help="width for the Savitzky-Golay filter") +parser.add_argument("radial_poly_order", type=int, help="polynomial order for the Savitzky-Golay filter") +parser.add_argument("output_file", type=str, help="output xi text file") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.xi.smoothen import smoothen_xi_files -infile = str(sys.argv[1]) -max_l = int(sys.argv[2]) -window_length = int(sys.argv[3]) -polyorder = int(sys.argv[4]) -outfile = str(sys.argv[5]) - -smoothen_xi_files(infile, max_l, window_length, polyorder, outfile) \ No newline at end of file +smoothen_xi_files(args.input_file, args.max_l, args.radial_window_length, args.radial_pol_yorder, args.output_file) \ No newline at end of file diff --git a/scripts/take_subset_of_particles.py b/scripts/take_subset_of_particles.py index 0df1e1dd..b9f194ab 100644 --- a/scripts/take_subset_of_particles.py +++ b/scripts/take_subset_of_particles.py @@ -1,17 +1,15 @@ """This script takes an input txt file listing the (x,y,z,w,j) coordinates of random particles and selects N random particles from this, which it writes to file.""" -import sys +import argparse -if len(sys.argv) != 4: - print("Usage: python take_subset_of_particles.py {INPUT_FILE} {OUTPUT_FILE} {N_PARTICLES}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="This script takes an input txt file listing the (x,y,z,w,j) coordinates of random particles, selects N random particles and writes them to the output file.") +parser.add_argument("input_file", type=str) +parser.add_argument("output_file", type=str) +parser.add_argument("N", type=int, help="number of particles to randomly select") +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.pre_process.take_subset_of_particles import take_subset_of_particles -infile_name = str(sys.argv[1]) -outfile_name = str(sys.argv[2]) -N = int(sys.argv[3]) - -take_subset_of_particles(infile_name, outfile_name, N) \ No newline at end of file +take_subset_of_particles(args.input_file, args.output_file, args.N) \ No newline at end of file diff --git a/scripts/utils/__init__.py b/scripts/utils/__init__.py index 8ef80952..a8a5d068 100644 --- a/scripts/utils/__init__.py +++ b/scripts/utils/__init__.py @@ -3,9 +3,3 @@ def adjust_path(): # insert the parent directory, override the one installed in the system sys.path.insert(0, os.path.realpath(os.path.join(os.path.dirname(__file__), "../.."))) - - -def get_arg_safe(index: int, type = str, default: object = None) -> object: - # get argument by index from sys.argv and convert it to the requested type if there are enough elements there - # otherwise return the default value - return type(sys.argv[index]) if len(sys.argv) > index else default diff --git a/scripts/write_binning_file_hybrid.py b/scripts/write_binning_file_hybrid.py index eeebefd7..db3a93f8 100644 --- a/scripts/write_binning_file_hybrid.py +++ b/scripts/write_binning_file_hybrid.py @@ -1,20 +1,18 @@ -## Utility function to create the Corrfunc radial-binning file for a given number of radial bins in a certain range of comoving radii. +## Utility script to create the Corrfunc radial binning file for a given number of logarithmic and then linear radial bins in a certain range of comoving radii. -import sys +import argparse -if len(sys.argv) < 7: - print("Please specify input parameters in the form {N_LOG_BINS} {N_LIN_BINS} {MIN_R} {CUTOFF_R} {MAX_R} {OUTPUT_FILE}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="Utility script to create the Corrfunc radial binning file with logarithmic bins at smaller radii and linear bins at larger.") +parser.add_argument("n_log_r_bins", type=int, help="number of logarithmic radial bins") +parser.add_argument("n_lin_r_bins", type=int, help="number of linear radial bins") +parser.add_argument("r_min", type=float) +parser.add_argument("r_cut", type=float, help="boundary between the logarithmic and linear bins") +parser.add_argument("r_max", type=float) +parser.add_argument("output_file", type=str) +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.write_binning_file import write_binning_file_hybrid -n_log_bins = int(sys.argv[1]) -n_lin_bins = int(sys.argv[2]) -r_min = float(sys.argv[3]) -r_cut = float(sys.argv[4]) -r_max = float(sys.argv[5]) -out_file = str(sys.argv[6]) - -write_binning_file_hybrid(out_file, r_min, r_cut, r_max, n_log_bins, n_lin_bins) \ No newline at end of file +write_binning_file_hybrid(args.output_file, args.r_min, args.r_cut, args.r_max, args.n_log_r_bins, args.n_lin_r_bins) \ No newline at end of file diff --git a/scripts/write_binning_file_linear.py b/scripts/write_binning_file_linear.py index 307ab81e..d44adfad 100644 --- a/scripts/write_binning_file_linear.py +++ b/scripts/write_binning_file_linear.py @@ -1,18 +1,16 @@ -## Utility function to create the Corrfunc radial-binning file for a given number of radial bins in a certain range of comoving radii. +## Utility script to create the Corrfunc radial binning file for a given number of linear radial bins in a certain range of comoving radii. -import sys +import argparse -if len(sys.argv) < 5: - print("Please specify input parameters in the form {N_BINS} {MIN_R} {MAX_R} {OUTPUT_FILE}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="Utility script to create the Corrfunc radial binning file with linear bins.") +parser.add_argument("n_r_bins", type=int, help="number of radial bins") +parser.add_argument("r_min", type=float) +parser.add_argument("r_max", type=float) +parser.add_argument("output_file", type=str) +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.write_binning_file import write_binning_file_linear -nrbins = int(sys.argv[1]) -r_min = float(sys.argv[2]) -r_max = float(sys.argv[3]) -out_file = str(sys.argv[4]) - -write_binning_file_linear(out_file, r_min, r_max, nrbins) \ No newline at end of file +write_binning_file_linear(args.output_file, args.r_min, args.r_max, args.n_r_bins) \ No newline at end of file diff --git a/scripts/write_binning_file_log.py b/scripts/write_binning_file_log.py index bb2be208..fd911270 100644 --- a/scripts/write_binning_file_log.py +++ b/scripts/write_binning_file_log.py @@ -1,18 +1,16 @@ -## Utility function to create the Corrfunc radial-binning file for a given number of radial bins in a certain range of comoving radii. +## Utility script to create the Corrfunc radial binning file for a given number of logarithmic radial bins in a certain range of comoving radii. -import sys +import argparse -if len(sys.argv) < 5: - print("Please specify input parameters in the form {N_LOG_BINS} {MIN_R} {MAX_R} {OUTPUT_FILE}.") - sys.exit(1) +parser = argparse.ArgumentParser(description="Utility script to create the Corrfunc radial binning file with logarithmic bins.") +parser.add_argument("n_r_bins", type=int, help="number of radial bins") +parser.add_argument("r_min", type=float) +parser.add_argument("r_max", type=float) +parser.add_argument("output_file", type=str) +args = parser.parse_args() from utils import adjust_path adjust_path() from RascalC.write_binning_file import write_binning_file_log -nrbins = int(sys.argv[1]) -r_min = float(sys.argv[2]) -r_max = float(sys.argv[3]) -out_file = str(sys.argv[4]) - -write_binning_file_log(out_file, r_min, r_max, nrbins) \ No newline at end of file +write_binning_file_log(args.output_file, args.r_min, args.r_max, args.n_r_bins) \ No newline at end of file diff --git a/triple_counts/Makefile b/triple_counts/Makefile index a1af9197..776a228b 100644 --- a/triple_counts/Makefile +++ b/triple_counts/Makefile @@ -4,14 +4,10 @@ CC = gcc CFLAGS = -O3 -Wall -MMD CXX = g++ CXXFLAGS = -O3 -Wall -MMD -std=c++11 -ffast-math $(shell pkg-config --cflags gsl) -CXXFLAGS += -DOPENMP -DTHREE_PCF -DPRINTPERCENTS +CXXFLAGS += -DOPENMP -DTRIPLE -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) -#-DLEGENDRE_MIX # also compute 2PCF covariances in Legendre bins, but in other, "mixed" mode, corresponding to projection of s,µ bins into multipoles -# without either of the two Legendre flags above, the covariance is computed in s,µ bins -#-DJACKKNIFE # use this to compute (r,mu)-space 2PCF covariances and jackknife covariances. Incompatible with -DLEGENDRE but works with -DLEGENDRE_MIX -#-DTHREE_PCF # use this to compute 3PCF autocovariances +#-DTRIPLE # need this to have mumin=-1 and not 0 to be consistent with the 3PCF covariance #-DPRINTPERCENTS # use this to print percentage of progress in each loop. This can be a lot of output LFLAGS = $(shell pkg-config --libs gsl) # common part diff --git a/triple_counts/triple_counts.cpp b/triple_counts/triple_counts.cpp index e710422f..77419e28 100644 --- a/triple_counts/triple_counts.cpp +++ b/triple_counts/triple_counts.cpp @@ -309,10 +309,14 @@ class compute_triples{ if (par->random_seed) { std::random_device urandom("/dev/urandom"); std::uniform_int_distribution dist(0, seed_step-1); // distribution of integers with these limits, inclusive - seed_shift = dist(urandom); // for each thread, it will be seed = seed_step * n_loops + seed_shift, where n_loops goes from 0 to max_loops-1 => no overflow and guaranteed unique for each thread + seed_shift = dist(urandom); // for each iteration, it will be seed = seed_step * n_loops + seed_shift, where n_loops goes from 0 to max_loops-1 => no overflow and guaranteed unique for each thread + par->seed = seed_shift; // save the random seed + par->random_seed = false; // avoid random seed resetting in the future — different seeds for xi rescaling and the main computation are much harder to reproduce } + printf("# INFO: the base RNG seed is %lu, incremented by %lu in each integration loop (matching these should guarantee the same results in a completed computation with the same input catalogs and correlation functions).\n", seed_shift, seed_step); gsl_rng_env_setup(); // initialize gsl rng + int completed_loops = 0; // counter of completed loops, since may not finish in index order uint64 tot_pairs=0, tot_triples=0; // global number of particle pairs/triples/quads used (including those rejected for being in the wrong bins) uint64 cell_attempt2=0,cell_attempt3=0; // number of j,k,l cells attempted @@ -351,7 +355,9 @@ class compute_triples{ int* prim_ids; // list of particle IDs in primary cell double p2,p3; // probabilities int mnp = grid->maxnp; // max number of particles in a grid cell +#ifdef PRINTPERCENTS Float percent_counter; +#endif int x, prim_id_1D; integer3 delta2, delta3, prim_id, sec_id, thi_id; Float3 cell_sep2,cell_sep3; @@ -373,7 +379,9 @@ class compute_triples{ #pragma omp for schedule(dynamic) #endif for (int n_loops = 0; n_loopsmax_loops; n_loops++){ +#ifdef PRINTPERCENTS percent_counter=0.; +#endif loc_used_pairs=0; loc_used_triples=0; // Set/reset the RNG seed based on loop number instead of thread number to reproduce results with different number of threads but other parameters kept the same. Individual subsamples may differ because they are accumulated/written in order of loop completion which may depend on external factors at runtime, but the final integrals should be the same. @@ -383,11 +391,13 @@ class compute_triples{ // LOOP OVER ALL FILLED I CELLS for (int n1=0; n1nf;n1++){ +#ifdef PRINTPERCENTS // Print time left if((float(n1)/float(grid->nf)*100)>=percent_counter){ - printf("Run %d of %d on thread %d: Using cell %d of %d - %.0f percent complete\n",1+n_loops/par->nthread, int(ceil(float(par->max_loops)/(float)par->nthread)),thread, n1+1,grid->nf,percent_counter); + printf("Iteration %d of %d on thread %d: Using cell %d of %d - %.0f percent complete\n", n_loops, par->max_loops, thread, n1+1, grid->nf, percent_counter); percent_counter+=5.; } +#endif // Pick first particle prim_id_1D = grid-> filled[n1]; // 1d ID for cell i @@ -445,11 +455,13 @@ class compute_triples{ #pragma omp critical // only one processor can access at once #endif { - if ((n_loops+1)%par->nthread==0){ // Print every nthread loops + completed_loops++; // increment completed loops counter, since they may be done not according to n_loops order + printf("Iteration %d of %d on thread %d completed (%d/%d)\n", n_loops, par->max_loops, thread, completed_loops, par->max_loops); + if (completed_loops % par->nthread == 0) { // Print every nthread completed loops TotalTime.Stop(); // interrupt timing to access .Elapsed() int current_runtime = TotalTime.Elapsed(); - int remaining_time = current_runtime/((n_loops+1)/par->nthread)*(par->max_loops/par->nthread-(n_loops+1)/par->nthread); // estimated remaining time - fprintf(stderr,"\nFinished integral loop %d of %d after %d s. Estimated time left: %2.2d:%2.2d:%2.2d hms, i.e. %d s.\n",n_loops+1,par->max_loops, current_runtime,remaining_time/3600,remaining_time/60%60, remaining_time%60,remaining_time); + int remaining_time = current_runtime*(par->max_loops - completed_loops)/completed_loops; // estimated remaining time + fprintf(stderr, "\nFinished %d integral loops of %d after %d s. Estimated time left: %2.2d:%2.2d:%2.2d hms, i.e. %d s.\n", completed_loops, par->max_loops, current_runtime, remaining_time/3600, remaining_time/60%60, remaining_time%60, remaining_time); TotalTime.Start(); // Restart the timer }