On this page... (hide)

  1.   1.  2D Weighted mean estimator :
  2.   2.  implementation :
  3.   3.  python code :
  4.   4.  directory :

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/