fix: NaMaster linear binning lmax mismatch and add fine pseudo-Cl rule

cailmdaley · claude · cailmdaley · commit 089a40d4647e · 2025-12-16T22:39:37.000+01:00
- Fix linear binning to respect actual lmax (2*nside-1) instead of using from_nside_linear which assumes 3*nside-1. This was causing ValueError when bins lmax (3071) didn't match fields lmax (2047). - Add fine_pseudo_cl Snakemake rule for generating finely-binned pseudo-Cls with configurable ell_step for accurate C_ℓ → COSEBIS conversion. - Add generate_fine_pseudo_cl.py script for the new rule. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
diff --git a/cosmo_inference/notebooks/2D_bmodes_paper_workflow/rules/bmodes.smk b/cosmo_inference/notebooks/2D_bmodes_paper_workflow/rules/bmodes.smk
@@ -151,3 +151,56 @@ localrules:
     bmodes_all,
     bmodes_pdf,
     # Note: plot_eb and plot_cosebis removed from localrules - they now run on compute nodes
+
+
+# Output directory for cosmo_val results
+COSMO_VAL_OUTPUT = "/n17data/cdaley/unions/pure_eb/code/sp_validation/notebooks/cosmo_val/output"
+CAT_CONFIG = "/n17data/cdaley/unions/pure_eb/code/sp_validation/notebooks/cosmo_val/cat_config.yaml"
+
+
+rule fine_pseudo_cl:
+    """Generate finely-binned pseudo-Cls for accurate C_ℓ → COSEBIS conversion.
+
+    Uses linear binning with configurable ell_step (default 1 for single-ell bins).
+    """
+    output:
+        pseudo_cl=f"{COSMO_VAL_OUTPUT}/pseudo_cl_{{version}}_ellstep={{ell_step}}.fits",
+        pseudo_cl_cov=f"{COSMO_VAL_OUTPUT}/pseudo_cl_cov_{{version}}_ellstep={{ell_step}}.fits",
+    params:
+        version="{version}",
+        ell_step=lambda w: int(w.ell_step),
+        cat_config=CAT_CONFIG,
+        nside=1024,
+        npatch=1,
+    resources:
+        mem_mb=32000,
+        runtime=120,  # minutes
+    threads: 48
+    script:
+        "../scripts/generate_fine_pseudo_cl.py"
+
+
+rule fine_pseudo_cl_fiducial:
+    """Generate finely-binned pseudo-Cls for fiducial version with ell_step=1."""
+    input:
+        rules.fine_pseudo_cl.output.pseudo_cl.format(
+            version=config["fiducial"]["version"],
+            ell_step=1,
+        ),
+        rules.fine_pseudo_cl.output.pseudo_cl_cov.format(
+            version=config["fiducial"]["version"],
+            ell_step=1,
+        ),
+
+
+rule fine_pseudo_cl_all:
+    """Generate finely-binned pseudo-Cls for all versions."""
+    input:
+        [
+            rules.fine_pseudo_cl.output.pseudo_cl.format(version=ver, ell_step=1)
+            for ver in bmodes_versions()
+        ],
+        [
+            rules.fine_pseudo_cl.output.pseudo_cl_cov.format(version=ver, ell_step=1)
+            for ver in bmodes_versions()
+        ],
diff --git a/cosmo_inference/notebooks/2D_bmodes_paper_workflow/rules/epistemics_v2.smk b/cosmo_inference/notebooks/2D_bmodes_paper_workflow/rules/epistemics_v2.smk
@@ -199,7 +199,7 @@ rule cl_version_comparison:
             f"{CONFIG_DIR}/cl_fiducial.md",
         ],
         config=f"{CONFIG_DIR}/config.yaml",
-        cl_fiducial_evidence=f"{CLAIMS_DIR}/cl_fiducial/evidence.json",
+        cl_fiducial_evidence=rules.cl_data_vector.output.evidence,
         pseudo_cl=[
             f"{SASHA_COSMO_VAL_OUTPUT}/pseudo_cl_{ver}.fits"
             for ver in config["versions"]
@@ -231,7 +231,7 @@ rule pure_eb_pte_matrix:
             f"{CONFIG_DIR}/2d_plots.md",
         ],
         config=f"{CONFIG_DIR}/config.yaml",
-        pure_eb_evidence=f"{CLAIMS_DIR}/pure_eb_data_vector/evidence.json",
+        pure_eb_evidence=rules.pure_eb_data_vector.output.evidence,
         # Pre-computed PTE matrices (from calculate_pure_eb_ptes.py)
         pte_data=f"results/paper_plots/intermediate/{config['fiducial']['version']}_pure_eb_ptes.npz",
     output:
@@ -424,12 +424,13 @@ rule harmonic_config_cosebis_comparison:
     input:
         spec=f"{CONFIG_DIR}/harmonic_config_cosebis_comparison.md",
         config=f"{CONFIG_DIR}/config.yaml",
+        # Finely-binned pseudo-Cls (nellbins=1000, linear Δℓ=16) for accurate COSEBIS conversion
         pseudo_cls=[
-            f"{SASHA_COSMO_VAL_OUTPUT}/pseudo_cl_{ver}.fits"
+            f"{COSMO_VAL_OUTPUT}/pseudo_cl_{ver}_nellbins=1000.fits"
             for ver in config["versions"]
         ],
         pseudo_cls_cov=[
-            f"{SASHA_COSMO_VAL_OUTPUT}/pseudo_cl_cov_{ver}.fits"
+            f"{COSMO_VAL_OUTPUT}/pseudo_cl_cov_{ver}_nellbins=1000.fits"
             for ver in config["versions"]
         ],
         cov_integration=[
@@ -474,9 +475,9 @@ rule xi_cosmology_paper:
     """
     input:
         spec=f"{CONFIG_DIR}/xi_cosmology_paper.md",
-        cosebis_evidence=f"{CLAIMS_DIR}/cosebis_version_comparison/evidence.json",
-        pure_eb_evidence=f"{CLAIMS_DIR}/pure_eb_data_vector/evidence.json",
-        covariance_evidence=f"{CLAIMS_DIR}/covariance_blind_consistency/evidence.json",
+        cosebis_evidence=rules.cosebis_version_comparison.output.evidence,
+        pure_eb_evidence=rules.pure_eb_data_vector.output.evidence,
+        covariance_evidence=rules.covariance_blind_consistency.output.evidence,
     output:
         macros="docs/unions_release/unions_2d_shear_xi/claims_macros.tex",
         evidence=f"{CLAIMS_DIR}/xi_cosmology_paper/evidence.json",
@@ -489,11 +490,11 @@ rule xi_cosmology_paper:
 rule paper_macros:
     """Generate LaTeX macros and tables for B-modes paper (Daley et al.)."""
     input:
-        cosebis_evidence=f"{CLAIMS_DIR}/cosebis_version_comparison/evidence.json",
-        pure_eb_evidence=f"{CLAIMS_DIR}/pure_eb_data_vector/evidence.json",
+        cosebis_evidence=rules.cosebis_version_comparison.output.evidence,
+        pure_eb_evidence=rules.pure_eb_data_vector.output.evidence,
         # PTE composite evidence for table generation
-        config_space_pte=f"{CLAIMS_DIR}/config_space_pte_matrices/evidence.json",
-        harmonic_space_pte=f"{CLAIMS_DIR}/harmonic_space_pte_matrices/evidence.json",
+        config_space_pte=rules.config_space_pte_matrices.output.evidence,
+        harmonic_space_pte=rules.harmonic_space_pte_matrices.output.evidence,
     output:
         bmodes="docs/unions_release/unions_bmodes/claims_macros.tex",
         pte_table_results="docs/unions_release/unions_bmodes/pte_table_results.tex",
@@ -536,17 +537,35 @@ rule all_method_specs:
 rule bmodes_paper_spec:
     """Generate evidence.json for bmodes_paper spec.
 
-    Dependencies parsed from markdown 'Depends on:' line.
+    Dependencies include both evidence files and figure outputs to ensure
+    dashboard regenerates all stale plots.
     """
     input:
         spec=f"{CONFIG_DIR}/bmodes_paper.md",
-        # Explicit upstream evidence (from Depends on: line in spec)
-        pure_eb_covariance=f"{CLAIMS_DIR}/pure_eb_covariance/evidence.json",
-        pure_eb_data_vector=f"{CLAIMS_DIR}/pure_eb_data_vector/evidence.json",
-        cosebis_version_comparison=f"{CLAIMS_DIR}/cosebis_version_comparison/evidence.json",
-        cl_fiducial=f"{CLAIMS_DIR}/cl_fiducial/evidence.json",
-        config_space_pte=f"{CLAIMS_DIR}/config_space_pte_matrices/evidence.json",
-        harmonic_space_pte=f"{CLAIMS_DIR}/harmonic_space_pte_matrices/evidence.json",
+        # Upstream evidence (using rules.X.output for single source of truth)
+        pure_eb_covariance=rules.pure_eb_covariance.output.evidence,
+        pure_eb_data_vector=rules.pure_eb_data_vector.output.evidence,
+        cosebis_version_comparison=rules.cosebis_version_comparison.output.evidence,
+        cl_fiducial=rules.cl_data_vector.output.evidence,
+        config_space_pte=rules.config_space_pte_matrices.output.evidence,
+        harmonic_space_pte=rules.harmonic_space_pte_matrices.output.evidence,
+        # Figure dependencies from bmodes.smk (ensures dashboard triggers regeneration)
+        plot_eb_xis=expand(
+            "/n17data/cdaley/unions/pure_eb/code/sp_validation/notebooks/cosmo_val/output/{version}_eb_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_minsepint={min_sep_int}_maxsepint={max_sep_int}_nbinsint={nbins_int}_npatch={npatch}_varmethod=semi-analytic_xis.png",
+            version=config["versions"],
+            min_sep=config["fiducial"]["min_sep"],
+            max_sep=config["fiducial"]["max_sep"],
+            nbins=config["fiducial"]["nbins"],
+            min_sep_int=config["fiducial"]["min_sep_int"],
+            max_sep_int=config["fiducial"]["max_sep_int"],
+            nbins_int=config["fiducial"]["nbins_int"],
+            npatch=config["fiducial"]["npatch"],
+        ),
+        plot_cosebis_modes=expand(
+            "/n17data/cdaley/unions/pure_eb/code/sp_validation/notebooks/cosmo_val/output/{version}_cosebis_minsep=" + str(config['fiducial']['min_sep_int']) + "_maxsep=" + str(config['fiducial']['max_sep_int']) + "_nbins=" + str(config['fiducial']['nbins_int']) + "_npatch={npatch}_varmethod=analytic_nmodes=" + str(config['fiducial']['nmodes']) + "_scalecut=" + str(config['fiducial']['fiducial_min_scale']) + "-" + str(config['fiducial']['fiducial_max_scale']) + "_cosebis.png",
+            version=config["versions"],
+            npatch=config["fiducial"]["npatch"],
+        ),
     output:
         evidence=f"{CLAIMS_DIR}/bmodes_paper/evidence.json",
     script:
diff --git a/cosmo_inference/notebooks/2D_bmodes_paper_workflow/scripts/generate_fine_pseudo_cl.py b/cosmo_inference/notebooks/2D_bmodes_paper_workflow/scripts/generate_fine_pseudo_cl.py
@@ -0,0 +1,109 @@
+"""Generate finely-binned pseudo-Cls for COSEBIS comparison.
+
+Uses linear binning with configurable ell_step for accurate C_ℓ → COSEBIS conversion.
+Designed to be called by Snakemake or standalone.
+"""
+
+import os
+import sys
+
+from astropy.io import fits
+import numpy as np
+
+from sp_validation.cosmo_val import CosmologyValidation
+
+
+def generate_pseudo_cl(
+    version: str,
+    ell_step: int,
+    output_cl: str,
+    output_cov: str,
+    cat_config: str,
+    nside: int = 1024,
+    npatch: int = 1,
+):
+    """Generate pseudo-Cl and covariance with linear binning.
+
+    Parameters
+    ----------
+    version : str
+        Catalog version (e.g., "SP_v1.4.6_leak_corr")
+    ell_step : int
+        Width of ell bins (1 for single-ell bins)
+    output_cl : str
+        Output path for pseudo-Cl FITS file
+    output_cov : str
+        Output path for covariance FITS file
+    cat_config : str
+        Path to catalog configuration YAML
+    nside : int
+        HEALPix nside for map-based estimation
+    npatch : int
+        Number of jackknife patches
+    """
+    output_dir = os.path.dirname(output_cl)
+    os.makedirs(output_dir, exist_ok=True)
+
+    print(f"\n{'='*60}")
+    print(f"Generating pseudo-Cl for {version} with ell_step={ell_step}")
+    print(f"{'='*60}\n")
+
+    # Create CosmologyValidation with fine ell binning
+    cv = CosmologyValidation(
+        versions=[version],
+        catalog_config=cat_config,
+        output_dir=output_dir,
+        binning="linear",
+        ell_step=ell_step,
+        nside=nside,
+        cell_method="map",
+        nrandom_cell=10,
+        # Standard params (not used for pseudo-Cl but required)
+        npatch=npatch,
+        theta_min=1.0,
+        theta_max=250.0,
+        nbins=20,
+    )
+
+    # Calculate pseudo-Cls
+    cv.calculate_pseudo_cl()
+
+    # Rename to final output
+    src_cl = os.path.join(output_dir, f"pseudo_cl_{version}.fits")
+    if os.path.exists(src_cl) and src_cl != output_cl:
+        with fits.open(src_cl) as hdul:
+            data = hdul["PSEUDO_CELL"].data
+            n_ell = len(data["ELL"])
+            print(f"Generated pseudo-Cl with {n_ell} ell bins")
+            print(f"ell range: [{data['ELL'].min():.1f}, {data['ELL'].max():.1f}]")
+        os.rename(src_cl, output_cl)
+        print(f"Saved to: {output_cl}")
+
+    # Calculate covariance
+    print(f"\nCalculating covariance...")
+    cv.calculate_pseudo_cl_eb_cov()
+
+    src_cov = os.path.join(output_dir, f"pseudo_cl_cov_{version}.fits")
+    if os.path.exists(src_cov) and src_cov != output_cov:
+        os.rename(src_cov, output_cov)
+        print(f"Saved covariance to: {output_cov}")
+
+
+if __name__ == "__main__":
+    # Check if running under Snakemake
+    try:
+        snakemake  # noqa: F821
+    except NameError:
+        # Standalone usage
+        print("Usage: Run via Snakemake rule 'fine_pseudo_cl'")
+        sys.exit(1)
+
+    generate_pseudo_cl(
+        version=snakemake.params.version,  # noqa: F821
+        ell_step=snakemake.params.ell_step,  # noqa: F821
+        output_cl=snakemake.output.pseudo_cl,  # noqa: F821
+        output_cov=snakemake.output.pseudo_cl_cov,  # noqa: F821
+        cat_config=snakemake.params.cat_config,  # noqa: F821
+        nside=snakemake.params.nside,  # noqa: F821
+        npatch=snakemake.params.npatch,  # noqa: F821
+    )
diff --git a/src/sp_validation/cosmo_val.py b/src/sp_validation/cosmo_val.py
@@ -2637,11 +2637,14 @@ def calculate_pseudo_cl_eb_cov(self):
                 lmax = 2*self.nside
                 b_lmax = lmax - 1
 
+                ells = np.arange(lmin, lmax+1)
+
                 if self.binning == 'linear':
-                    b = nmt.NmtBin.from_nside_linear(self.nside, self.ell_step)
+                    # Linear bands of width ell_step, respecting actual lmax
+                    bpws = (ells - lmin) // self.ell_step
+                    bpws = np.minimum(bpws, bpws[-1])  # Ensure last bin captures all
+                    b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax)
                 elif self.binning == 'powspace':
-                    ells = np.arange(lmin, lmax+1)
-
                     start = np.power(lmin, self.power)
                     end = np.power(lmax, self.power)
                     bins_ell = np.power(np.linspace(start, end, self.n_ell_bins+1), 1/self.power)
@@ -3032,11 +3035,14 @@ def get_pseudo_cls_map(self, map, wsp=None):
         lmax = 2*self.nside
         b_lmax = lmax - 1
 
+        ells = np.arange(lmin, lmax+1)
+
         if self.binning == 'linear':
-            b = nmt.NmtBin.from_nside_linear(self.nside, self.ell_step)
+            # Linear bands of width ell_step, respecting actual lmax
+            bpws = (ells - lmin) // self.ell_step
+            bpws = np.minimum(bpws, bpws[-1])  # Ensure last bin captures all
+            b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax)
         elif self.binning == 'powspace':
-            ells = np.arange(lmin, lmax+1)
-
             start = np.power(lmin, self.power)
             end = np.power(lmax, self.power)
             bins_ell = np.power(np.linspace(start, end, self.n_ell_bins+1), 1/self.power)
@@ -3053,7 +3059,7 @@ def get_pseudo_cls_map(self, map, wsp=None):
         factor = -1 if self.pol_factor else 1
 
         f_all = nmt.NmtField(mask=(map!=0), maps=[map.real, factor*map.imag], lmax=b_lmax)
-        
+
         if wsp is None:
             wsp = nmt.NmtWorkspace.from_fields(f_all, f_all, b)
         
@@ -3071,11 +3077,14 @@ def get_pseudo_cls_catalog(self, catalog, params, wsp=None):
         lmax = 2*self.nside
         b_lmax = lmax - 1
 
+        ells = np.arange(lmin, lmax+1)
+
         if self.binning == 'linear':
-            b = nmt.NmtBin.from_nside_linear(self.nside, self.ell_step)
+            # Linear bands of width ell_step, respecting actual lmax
+            bpws = (ells - lmin) // self.ell_step
+            bpws = np.minimum(bpws, bpws[-1])  # Ensure last bin captures all
+            b = nmt.NmtBin(ells=ells, bpws=bpws, lmax=b_lmax)
         elif self.binning == 'powspace':
-            ells = np.arange(lmin, lmax+1)
-
             start = np.power(lmin, self.power)
             end = np.power(lmax, self.power)
             bins_ell = np.power(np.linspace(start, end, self.n_ell_bins+1), 1/self.power)
