Multipole order mismatch in get_coupling_matrix() routine

[Sylk1-9]

Hello,

When calling the new routine get_coupling_matrix(), on the B-modes purification test scrip for example, I get a mixing matrix with rows and column order that are mismatched.
I couldn't recover the initial EE-EE, EE-BB, etc.. blocks of the coupling matrix.

In addition, I compared a matrix from a workspace whose the polarisation purification option is on, with a matrix for which the purification option is off (purifyE/B = True/False). Both arrays are the same. I am surprised since I would have expected that the mixing kernel block EE-BB would be different.

Thanks,

Sylvain

[damonge]

Hi @Sylk1-9,

can you please post a code snippet or a script that reproduces this?

[Sylk1-9]

Yeah sure, here is an example from your script sample_pureB.py.
Actually, using this example, both Mll matrices seems different. So the second part of my issue is actually not an issue ;)
But their respective plot still show mismatched column, I think

`
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
import pymaster as nmt

nside=128

msk=np.ones(hp.nside2npix(nside))
th,ph=hp.pix2ang(nside,np.arange(hp.nside2npix(nside)))
ph[np.where(ph>np.pi)[0]]-=2*np.pi
msk[np.where((th<2.63) & (th>1.86) & (ph>-np.pi/4) & (ph<np.pi/4))[0]]=1.

msk_apo=nmt.mask_apodization(msk,10.0,apotype='C1')

#Select a binning scheme
b=nmt.NmtBin(nside,nlb=1)
leff=b.get_effective_ells()

l,cltt,clee,clbb,clte=np.loadtxt('cls.txt',unpack=True);

mp_t,mp_q,mp_u=hp.synfast([cltt,clee,clbb,clte],nside=nside,new=True,verbose=False)
f2_np=nmt.NmtField(msk_apo,[mp_q,mp_u],purify_e=False,purify_b=False)
f2_yp=nmt.NmtField(msk_apo,[mp_q,mp_u],purify_e=True,purify_b=True)

#We initialize two workspaces for the non-pure and pure fields:
w_np=nmt.NmtWorkspace(); w_np.compute_coupling_matrix(f2_np,f2_np,b)
w_yp=nmt.NmtWorkspace(); w_yp.compute_coupling_matrix(f2_yp,f2_yp,b)

Mll_yp = w_yp.get_coupling_matrix()
Mll_np = w_np.get_coupling_matrix()

plt.matshow(np.log10(abs(Mll_yp)))
plt.matshow(np.log10(abs(Mll_np)))
plt.show()
`
Here is a matrix plot of the Mll_yp

and a zoom

[damonge]

Hi @Sylk1-9
I think in that script the mask is essentially 1 everywhere (possibly a typo ones -> zeros). In that case I wouldn't be surprised if you you got essentially the same MCM regardless of purification.

See the attached script (modified from yours), which shows how to interpret the MCM returned by NaMaster. As you can see from the plot below, there are definitely differences between purified and unpurified!

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
import pymaster as nmt

nside=128

#Mask                                                                                                                                                                                                              
msk=np.zeros(hp.nside2npix(nside))
th,ph=hp.pix2ang(nside,np.arange(hp.nside2npix(nside)))
msk[th<np.pi/6]=1; msk=nmt.mask_apodization(msk,10.0,apotype='C1')

#Select a binning scheme                                                                                                                                                                                           
b=nmt.NmtBin(nside,nlb=int(1./np.mean(msk)))
leff=b.get_effective_ells()

f2_np=nmt.NmtField(msk,[msk,msk],purify_e=False,purify_b=False)
f2_yp=nmt.NmtField(msk,[msk,msk],purify_e=False,purify_b=True)

#We initialize two workspaces for the non-pure and pure fields:                                                                                                                                                    
w_np=nmt.NmtWorkspace(); w_np.compute_coupling_matrix(f2_np,f2_np,b)
w_yp=nmt.NmtWorkspace(); w_yp.compute_coupling_matrix(f2_yp,f2_yp,b)

Mll_yp = w_yp.get_coupling_matrix()
Mll_np = w_np.get_coupling_matrix()
Mll_yp-=np.diag(np.diag(Mll_yp)) #Subtract diagonal to make the differences clearer                                                                                                                                
Mll_np-=np.diag(np.diag(Mll_np))
Mll_yp=Mll_yp.reshape([3*nside,4,3*nside,4]) #Reshape to a more convenient form                                                                                                                                    
Mll_np=Mll_np.reshape([3*nside,4,3*nside,4])

plt.figure()
l0=50
plt.plot(np.arange(3*nside),Mll_yp[l0,3,:,3],'r-',label='With purification')
plt.plot(np.arange(3*nside),Mll_np[l0,3,:,3],'k--',label='Without purification')
plt.xlabel('$\\ell\'$',fontsize=16)
plt.ylabel('$M^{BB,\\ell\'}_{BB,\\ell=%d}-M^{BB,\\ell\'=%d}_{BB,\\ell=%d}$'%(l0,l0,l0),fontsize=16)
plt.legend(loc='upper right')
plt.xlim([0,3*nside])
plt.savefig("mcm_diff.png",bbox_inches='tight')
plt.show();

[damonge]

Let me know if that solved the issue!

[Sylk1-9]

Yes indeed, thank you for your response.
Also, I think I didn't get correctly the way the Mll matrix was sorted. Tell me if I am wrong, but when I resort the Mll matrix as follow and plot it

Mll_yp_rc = np.zeros_like(Mll_yp) 
for tag1 in [0, 1, 2, 3]:
	for tag2 in [0, 1, 2, 3]:
		for l1 in np.arange(3*nside_large):
			for l2 in np.arange(3*nside_large):
				Mll_yp_rc[tag1*(3*nside_large-1)+l1, tag2*(3*nside_large-1)+l2] = Mll_yp[4*l1 + tag1, 4*l2 + tag2]

matshow(log10(abs(Mll_yp_rc)))

I get

So the four block-rows/column entries are EE, EB, BE, and BB ? As defined in your paper.
Going left to right, starting from the first line, the first block is the EE-EE mixing, then EE-EB, EE-BE, EE-BB.
So the BB-EE mixing (leakage) is at the bottom left. This is consistant since for in my exemple, I only purify B modes, thus the BB-EE block values are lower than the EE-BB block. I think I got it.

[damonge]

Yes, that's correct. It follows the same way in which power spectra are stored in NaMaster (for each ell, all the different power spectrum combinations together).

I'll close this then, but feel free to reopen if you have more questions.