Error creating spin-2 field with purify_b=True and specified lmax

[cbischoff]

Using pymaster version 2.2.2

It seems to be necessary to set the lmax parameter of my NmtField objects equal to the lmax of my NmtBin object, which was created by specifying bin edges. But setting lmax doesn't seem to work for a spin-2 field with purify_b=True.

nside = 128
apod = np.ones(hp.nside2npix(nside))
q = np.random.randn(hp.nside2npix(nside))
u = np.random.randn(hp.nside2npix(nside))
spin2 = nmt.NmtField(apod, [q, u], purify_b=True) # No error
spin2 = nmt.NmtField(apod, [q, u], lmax=200) # No error
spin2 = nmt.NmtField(apod, [q, u], purify_b=True, lmax=200) # This gives an error with following traceback

----> 1 spin2 = nmt.NmtField(apod, [q, u], purify_b=True, lmax=200)

File ~/.venv/cmb/lib/python3.12/site-packages/pymaster/field.py:254, in NmtField.__init__(self, mask, maps, spin, templates, beam, purify_e, purify_b, n_iter, n_iter_mask, tol_pinv, wcs, lmax, lmax_mask, masked_on_input, lite)
    252 task = [self.pure_e, self.pure_b]
    253 alm_mask = self.get_mask_alms()
--> 254 self.alm, maps = self._purify(self.mask, alm_mask, maps_unmasked,
    255                               n_iter=n_iter, task=task)
    256 if w_temp and (not self.lite):
    257     alm_temp = np.array([self._purify(self.mask, alm_mask, t,
    258                                       n_iter=n_iter,
    259                                       task=task)[0]
    260                          for t in temp_unmasked])

File ~/.venv/cmb/lib/python3.12/site-packages/pymaster/field.py:398, in NmtField._purify(self, mask, alm_mask, maps_u, n_iter, task, return_maps)
    396 for ipol, purify in enumerate(task):
    397     if purify:
--> 398         alms[ipol] += hp.almxfl(palm[ipol],
    399                                 fl[:self.ainfo.lmax+1],
    400                                 mmax=self.ainfo.mmax)
    402 # 3. Spin-2 mask bit
    403 # Compute spin-0 mask
    404 # The extra minus sign is because of the scalar SHT below
    405 # (E-mode definition for spin=0).
    406 fl[2:] = -np.sqrt((ls[2:]+2.0)*(ls[2:]-1.0))

ValueError: operands could not be broadcast together with shapes (73920,) (20301,) (73920,)

I guess one way to work around this would be to add a junk bin to my binning scheme that extends up to 3*nside-1?

[damonge]

Sorry for the delay. Thanks for picking this one up! The reason for the failure is that purification requires operating on both the field and its mask, so you need to specify both lmax and lmax_mask.

I realise this may be inconvenient, and it is not clear to me that having to couple these two parameters is the right thing in this case. E.g. I could imagine that, for some masks, you may want to go to higher lmax when purifying, even if it means padding the field with zeros at higher ell while purifying...

In short, if the workaround above (setting both lmax=200 and lmax_mask=200) does not work for you (e.g. returns biased spectra because the mask is not sufficiently well resolved in harmonic space), do let me know.

[cbischoff]

This solved the problem. Thanks!

But now I have a related question: When creating a workspace, it seems to be necessary for the lmax of the fields to match the lmax of my binning scheme. Why is this necessary? The bandpower window functions returned by the workspace get_bandpower_windows() stops at this lmax, but I would generally expect that the window function support would extend beyond the nominal maximum ell of the highest bin. Again, this seems to be something that I could work around by adding a junk bin (or bins) that extends higher in ell than I really need. But maybe I am missing something?

[damonge]

Yes, perhaps this could be better documented, but bins in NaMaster are intended to define all the angular scales at which the power spectra will be calculated (even the ones you may not care about).

As you say, unless this is very inconvenient, adding a few garbage bins at the end should do the trick. This shouldn't affect the code speed, since the slowest step is either estimating the coupling matrix or purifying, both of which need to be carried out over all ells on which the window functions are calculated (which are more than you'd like to include in the bins, if I understand correctly).

[cbischoff]

Great. That makes sense.

While I have you here, I have one more question (getting further and further from the original GitHub issue): Suppose I have a workspace for the cross-spectrum between two spin-2 fields. The output of ws.get_bandpower_windows() has shape (4,nbin,4,lmax+1). If I look at the [0,i,3,:] slice, is that the window function for true sky EE leaking into my BB bandpower? Or is it true BB leaking into my EE bandpower? I think that the latter is correct, because if I turn on purify_b=True, then the [0,i,3,:] slice is non-zero but the [3,i,0,:] slice is identically zero (so I think this it the EE->BB leakage function).

This isn't clear from the function documentation (and I also looked up the reference equations from the NaMaster paper).

Thanks!

[damonge]

Yes, that's right. Essentially, the measured bandpowers and the true C_ell are related via:

# below, X and Y are indices going over {EE, EB, BE, BB}
bpw_measured[X, b1] = sum_{Y, ell2} MCM[X, b1, Y, ell2] * cl_true[Y, ell2]

Therefore MCM[0, b1, 3, ell2] is the amount of BB power in ell=ell2 that leaks into the b1-th EE bandpower.

[damonge]

(indeed the equations from the paper that are linked in the docs are not particularly enlightening about this point, so perhaps I should include the example above in the docs -- thanks for bringing this up)

[cbischoff]

Thanks for all the help! I'll close this issue now...

[damonge]

Just reopening and will close by fixing the docs.