Resolution degrading routine for polarization datasets, using 2D weighed average.

1.  2D Weighted mean estimator :

map m_i = [Q_i, U_i]

mean estimator for big pixel p : \chi^2 = \sum_{i\in p} ( \mu - m_i) C_i^{-1} ( \mu - m_i)^T

chi2 derivative : \partial_\mu \chi^2 = 2 \sum_{i\in p} ( \mu - m_i)^T C_i^{-1} = 0

optimal estimator : \mu = (\sum_{i\in p} C_i^{-1})^{-1} \cdot (\sum_{i\in p} C_i ^{-1} m_i)

2.  implementation :

• Bad pixel in datasets (UNSEEN healpix value), correspond to 0 in the datasets covariance matrices. Those must be taken into account when ud_grading datasets maps. We make use of the healpix ud_graded routine which average nearby pixels, while skipping UNSEEN pixels :

1. build a map of UNSEEN pixels in dataset
2. compute determinant | C |
3. put det[det==inf]=0, and det[UNSEEN]=0
4. compute Weight matrices W = C^{-1}
5. put W[UNSEEN] = 0
6. compute ud_graded weight matrices Wud = {\rm udgrade}(W)
7. put Wud[Wud == hp.UNSEEN] = 0
8. compute pWP = [Q,U].C^-1
9. put pWP[:, unseenmap] = hp.UNSEEN
10. Compute WP = ud_grade(pWP)
11. build a map of UNSEEN in WP
12. compute determinant of Wud
13. put dets[dets == 0] = np.inf
14. compute [Qud, Uud] = Wud^-1 . WP
15. compute normalized Wud = (nside_in/(1.0*nside_out))**2*Wud
16. compute determinant of Wud
17. put dets[dets == 0] = np.inf
18. compute Cud = Wud^-1

3.  python code :

Input map P = [Q, U] and C = [covQQ, covQU, covUU]

    unseenmap = (P[0] == hp.UNSEEN) + (P[1] == hp.UNSEEN)
    nside_in = hp.npix2nside(len(C[0]))
    dets = 1.0*(C[0]*C[2] - C[1]**2)
    dets[dets == 0] = np.inf
    dets[unseenmap] = np.inf

    # Compute Wud = ud(C^-1)
    W = np.array([C[2]/dets, -C[1]/dets, C[0]/dets])
    W[:, unseenmap] = hp.UNSEEN
    Wud = np.array([hp.ud_grade(W[0], nside_out, order_in=order, order_out=order),
                    hp.ud_grade(W[1], nside_out, order_in=order, order_out=order),
                    hp.ud_grade(W[2], nside_out, order_in=order, order_out=order)])
    Wud[Wud == hp.UNSEEN] = 0

    # Compute ud_grade([Q,U].C^-1)
    pWP = np.array([(C[2]*P[0] - C[1]*P[1])/dets, (C[1]*P[0] + C[0]*P[1])/dets])
    pWP[:, unseenmap] = hp.UNSEEN
    WP = np.array([hp.ud_grade(pWP[0], nside_out, order_in=order, order_out=order),
                   hp.ud_grade(pWP[1], nside_out, order_in=order, order_out=order)])
    #### WP[WP == hp.UNSEEN] = 0
    unseenmapWP = (WP[0] == hp.UNSEEN) + (WP[1] == hp.UNSEEN)

    dets = 1.0*(Wud[0]*Wud[2] - Wud[1]**2)
    dets[dets == 0] = np.inf
    Q = (Wud[2]*WP[0] - Wud[1]*WP[1])/dets
    U = (-Wud[1]*WP[0] + Wud[0]*WP[1])/dets
    Pud = np.array([Q, U])
    Pud[:, unseenmapWP] = hp.UNSEEN

    Wud = (nside_in/(1.0*nside_out))**2*Wud

    dets = 1.0*(Wud[0]*Wud[2] - Wud[1]**2)
    dets[dets == 0] = np.inf
    Cud = np.array([Wud[2], -Wud[1], Wud[0]])/dets

4.  directory :

/sps/planck/Users/vanneste/Work/Planck/maps/PR3/wudgraded/