Tomographic 2pcf calculation with treecorr#297
Conversation
|
Great to see this PR Lisa. Something I was considering doing after the harmonic space and rho-/tau-statistics so thanks a lot for taking this task over. I mostly have global implementation design comment that we should agree on and that mostly come from my work on harmonic space (#221) and rho-/tau-statistics (#295). Global design comment:
I will comment directly in the code for implementation related comments. |
| "var_method": "jackknife" if int(npatch) > 1 else "shot", | ||
| } | ||
|
|
||
| pol_factor = self.pol_factor |
There was a problem hiding this comment.
The polarisation factor is necessary to match NaMaster convention. I am not sure this is necessary for real space. Have you check it? I think the previous version was not using it.
In addition, if we need to flip sign, I think treecorr has built-in arguments for that we should probably use rather than multiplying manually.
There was a problem hiding this comment.
I see, ok I had thought it was included to add flexibility to account for the case where e2 flips sign. If we agree upstream that e2 doesn't require a sign flip thenI can remove it!
| out_fname = self._output_path( | ||
| f"{ver}_xi_minsep={treecorr_config['min_sep']}_maxsep={treecorr_config['max_sep']}_nbins={treecorr_config['nbins']}_npatch={npatch}.txt" | ||
| ) | ||
| if tomo_bin1 is None and tomo_bin2 is None: |
There was a problem hiding this comment.
In the global design logic, we can avoid having a if ... else statement for the tomography with some duplication of code. See e.g. here for an implementation design I adopted. We should homogeneise either to this choice or another one.
| # LG TO-DO: No longer writing out text file, change to sacc_io method | ||
| # gg.write(out_fname, write_patch_results=True, write_cov=True) | ||
|
|
||
| if not os.path.exists(patch_file_bin1) or not os.path.exists( |
There was a problem hiding this comment.
You left a comment about this elsewhere. I think we can be smart about the patches and use the same across all the tomographic bins. I don't see a reason why we would need different patches for the different tomo bins. From the first produced catalog it must be possible to extract the patch_centers and use it as input to any other treecorr.Catalog produced afterwards; hence the need to have a function that computes all the tomographic bins at once.
| # Add correlation object to class | ||
| if not hasattr(self, "cat_ggs"): | ||
| self.cat_ggs = {} | ||
| self.cat_ggs[ver] = gg |
There was a problem hiding this comment.
The saving dictionnary should now also specify the bins for the tomographic case.
| else: | ||
| print(f"Writing Map2 to output file {out_fname_map2} ") | ||
| gg.writeMapSq(out_fname_map2, R=theta_map, m2_uform="Schneider") | ||
| self._map2[ver][f"{tomo_bin1}_{tomo_bin2}"] = { |
There was a problem hiding this comment.
Sure, I'll follow the convention of f"tomo_bin_{bin_key1}_tomo_bin_{bin_key2}"?
There was a problem hiding this comment.
If it's fine for everyone. I am happy to take any suggestion
|
Thank you for your comments! Yes I see what you mean, I was contemplating between doing it this way (having the tomographic bins as input) vs giving it a compute_tomography as a boolean, then computing all the tomographic 2PCFs at one go. One reason why I chose to do it this way was that I saw other scripts (such as |
|
I think that what does not undermine other functions is to keep this function as a child function called by parents functions from |
cailmdaley
left a comment
cailmdaley
left a comment
There was a problem hiding this comment.
The core catalogue-splitting logic in calculate_2pcf looks right, and thanks for keeping the non-tomo default path intact. Most of my comments are about the interim I/O state pending sacc_io — all fine to defer, with one guard-rail — plus two small real bugs and a couple of convention alignments with #295/#298. And +1 to the design you and Sacha converged on (gating on a compute_tomography boolean inside CosmologyValidation); happy for the plotting/wiring bits to land in that future PR.
— Claude, on behalf of Cail
| print(f"Evaluating {len(generated_cuts)} scale cut combinations") | ||
|
|
||
| # Call b_modes function with scale cuts list | ||
| results = calculate_cosebis( |
There was a problem hiding this comment.
This rebinding clobbers the accumulator dict initialized above (results = {bin_key: None …}), and then line 161's results[bin_key] = results self-references — even the plain non-tomo single-scale-cut case returns a mangled dict. A differently-named inner variable (res = calculate_cosebis(...); results[bin_key] = res) should fix it — same in the else branch below.
| len(self.results[ver].dat_shear), dtype=bool | ||
| ) | ||
| else: | ||
| mask_bin1 = self.results[ver].dat_shear["tom_bin_id"] == bin1 |
There was a problem hiding this comment.
| ggs = {f"tomo_bin_{b1}_tomo_bin_{b2}": None for b1, b2 in tomo_bin_pairs} | ||
|
|
||
| # LG TO-DO: Change to sacc_io method | ||
| out_fname = self._output_path( |
There was a problem hiding this comment.
Since the writes are deferred to sacc_io, suggest also disabling this os.path.exists read branch until then: the filename has no bin pair in it, so a stale file from an old run would be read into every bin pair — silently identical data across all auto/cross combinations. (When sacc_io lands, encoding the pair à la #295's basename() solves it.) Same applies to the aperture-mass cache below, where the read branch additionally stores under a different key format/value type than the compute and plot code expect.
| if compute_tomography: | ||
| tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver) | ||
| self.print_magenta( | ||
| f"Computing tomographic ξ± for {ver} with {len(tomo_bin_pairs)} bins." |
There was a problem hiding this comment.
tiny: if a version lacks tomography info, _get_tomo_bins returns (None, None), so this len() raises TypeError before the intended, informative ValueError — moving the guard above the print fixes it (same in calculate_aperture_mass_dispersion).
| cov_path=cov_path, | ||
| scale_cuts=scale_cuts, | ||
| evaluate_all_scale_cuts=evaluate_all_scale_cuts, | ||
| compute_tomography=False, # LG: Hardcoded to False for now |
There was a problem hiding this comment.
Fine to leave for the future PR — just noting plot_cosebis also still assumes the old un-binned return shape of calculate_cosebis (results are now one level deeper, even non-tomographically), so COSEBIs plotting is broken until that pass; worth a TODO so it's not a surprise.
|
@LisaGoh just a note to point you to this function I coded. It allows to plot xi_+/xi_- like plot as long as you provide a function that says what to plot in each tomo bin axes. I was thiking this could be use for any plot involving 2pcfs of any sort so don't hesitate to use it/update it if it's useful and avoids code duplication. |
3 participants
This PR resolves #296. I've added functionalities in
calculate_2pcf()inreal_space.pyto calculate the 2PCFs given two tomographic bin IDs as input. Ditto forxip_xim.py. The functions default to the non-tomographic case if the bins are not defined.I've commented out all the functions that do the writing out of the 2PCFs, in preparation for the sacc_io output format instead.