How to use the healpy.ud_grade function in healpy

if npix > 12:
        # Determine which pixels comprise multi-pixel tiles.
        ipix = np.flatnonzero(
            (m[0::4] == m[1::4]) &
            (m[0::4] == m[2::4]) &
            (m[0::4] == m[3::4]))

        if len(ipix):
            ipix = (4 * ipix +
                np.expand_dims(np.arange(4, dtype=np.intp), 1)).T.ravel()

            nside = hp.npix2nside(npix)

            # Downsample.
            m_lores = hp.ud_grade(
                m, nside // 2, order_in='NESTED', order_out='NESTED')

            # Interpolate recursively.
            _interpolate_level(m_lores)

            # Record interpolated multi-pixel tiles.
            m[ipix] = hp.get_interp_val(
                m_lores, *hp.pix2ang(nside, ipix, nest=True), nest=True)

mask_apo[id0]=0.
    mask_apo[id1]=1.
    mask_apo[idb]=x-np.sin(2*np.pi*x)/(2*np.pi)
    return mask_apo


nside_out=64
mask=getmaskapoana(nside_out,20.,0.4,dec0=90.)
dl,dw_q,dw_u=hp.synfast([(cltt+nltt)[:3*nside_out],
                         (clee+nlee)[:3*nside_out],
                         (clbb+nlbb)[:3*nside_out],
                         (clte+nlte)[:3*nside_out]],
                        nside_out,new=True,verbose=False,pol=True)
sl=hp.ud_grade(hp.read_map("../../sandbox_validation/data/cont_lss_dust_ns64.fits",
                           field=0,verbose=False),nside_out=nside_out);
sw_q=hp.ud_grade(hp.read_map("../../sandbox_validation/data/cont_wl_psf_ns64.fits",
                             field=0,verbose=False),nside_out=nside_out);
sw_u=hp.ud_grade(hp.read_map("../../sandbox_validation/data/cont_wl_psf_ns64.fits",
                             field=1,verbose=False),nside_out=nside_out);
hp.write_map("msk.fits",mask,overwrite=True)
hp.write_map("mps.fits",[dl,dw_q,dw_u],overwrite=True)
hp.write_map("tmp.fits",[sl,sw_q,sw_u],overwrite=True)

b=nmt.NmtBin(nside_out,nlb=16)
leff=b.get_effective_ells()
lfull=np.arange(3*nside_out)

#No contaminants
prefix='bm_nc_np'
f0=nmt.NmtField(mask,[dl])
f2=nmt.NmtField(mask,[dw_q,dw_u])
w00=nmt.NmtWorkspace(); w00.compute_coupling_matrix(f0,f0,b);

mask_apo[idb]=x-np.sin(2*np.pi*x)/(2*np.pi)
    return mask_apo


nside_out=64
mask=getmaskapoana(nside_out,20.,0.4,dec0=90.)
dl,dw_q,dw_u=hp.synfast([(cltt+nltt)[:3*nside_out],
                         (clee+nlee)[:3*nside_out],
                         (clbb+nlbb)[:3*nside_out],
                         (clte+nlte)[:3*nside_out]],
                        nside_out,new=True,verbose=False,pol=True)
sl=hp.ud_grade(hp.read_map("../../sandbox_validation/data/cont_lss_dust_ns64.fits",
                           field=0,verbose=False),nside_out=nside_out);
sw_q=hp.ud_grade(hp.read_map("../../sandbox_validation/data/cont_wl_psf_ns64.fits",
                             field=0,verbose=False),nside_out=nside_out);
sw_u=hp.ud_grade(hp.read_map("../../sandbox_validation/data/cont_wl_psf_ns64.fits",
                             field=1,verbose=False),nside_out=nside_out);
hp.write_map("msk.fits",mask,overwrite=True)
hp.write_map("mps.fits",[dl,dw_q,dw_u],overwrite=True)
hp.write_map("tmp.fits",[sl,sw_q,sw_u],overwrite=True)

b=nmt.NmtBin(nside_out,nlb=16)
leff=b.get_effective_ells()
lfull=np.arange(3*nside_out)

#No contaminants
prefix='bm_nc_np'
f0=nmt.NmtField(mask,[dl])
f2=nmt.NmtField(mask,[dw_q,dw_u])
w00=nmt.NmtWorkspace(); w00.compute_coupling_matrix(f0,f0,b);
cw00=nmt.NmtCovarianceWorkspace(); cw00.compute_coupling_coefficients(f0,f0);
cw00.write_to(prefix+'_cw00.dat')

plt.plot(larr,clc,label='Star contaminant');
        plt.plot(larr,cld,label='True $\\delta_g$');
        plt.plot(larr,clt,label='True + contaminant');
        plt.xlabel('$\\ell$',fontsize=18);
        plt.ylabel('$C_\\ell$',fontsize=18);
        plt.legend(); plt.loglog(); plt.xlim([2,3*o.nside_out])

        hp.mollview(dmap*msk,title='True $\\delta_g$')
        hp.mollview(starmap_b*msk,title='Star contaminant')
        hp.mollview((starmap_b+dmap)*msk,title='True + contaminant')
        plt.show()
    
    hp.write_map("cont_lss_star_ns%d.fits"%o.nside_out,starmap_b,overwrite=True)

if not os.path.isfile("cont_lss_dust_ns%d.fits"%o.nside_out) :
    dust=hp.ud_grade(hp.read_map("lambda_sfd_ebv.fits",verbose=False),nside_out=o.nside_out)
    mean=np.amin(dust*msk)
    dust=mean-dust
    
    cl=hp.anafast(dust*msk)/fsky
    ratio=np.sqrt(0.2*np.sum(cltt[50:200])/np.sum(cl[50:200])) #20% contamination in the power spectrum at ell~100
    dust=dust*ratio

    if o.plot_stuff :
        cl=hp.anafast(dust*msk)/fsky
        plt.figure();
        plt.plot(larr,cl,label='Dust contaminant');
        plt.plot(l,cltt,label='True $\\delta_g$');
        plt.plot(larr,cl+cltt[:3*o.nside_out],label='True + contaminant')
        plt.xlabel('$\\ell$',fontsize=18);
        plt.ylabel('$C_\\ell$',fontsize=18);
        plt.legend(); plt.loglog(); plt.xlim([2,3*o.nside_out])

parser.add_option('--plot', dest='plot_stuff', default=False, action='store_true',
                  help='Set if you want to produce plots')
parser.add_option('--nside', dest='nside_out', default=512, type=int,
                  help='Resolution parameter')
(o, args) = parser.parse_args()

l,cltt,clee,clbb,clte,nltt,nlee,nlbb,nlte=np.loadtxt("cls_lss.txt",unpack=True)
msk=hp.ud_grade(hp.read_map("mask_lss.fits",verbose=False),nside_out=o.nside_out); fsky=np.mean(msk)

larr=np.arange(3*o.nside_out)

if not os.path.isfile("cont_lss_nvar_ns%d.fits"%o.nside_out) :
    #First, variations in \bar{n} due to stars
    starmap=hp.read_map("star_template.fits",verbose=False)
    starmap[starmap<=0]=1E-15
    starmap_low=hp.smoothing(hp.ud_grade(hp.ud_grade(starmap,nside_out=64),nside_out=o.nside_out),fwhm=3.5*np.pi/180,verbose=False) #Smooth out the star map to 3.5deg FWHM
    r_starmap=starmap_low/np.amax(starmap_low[msk>0]) #Normalize by maximum value
    shotnoise_factor=(r_starmap-4)/(2*r_starmap-4) #Compute shot noise factor. This interpolates between 1 (no stars) and 1.5 (lots of stars)

    #Now, impose additional depth variations on scales of 3.5 degrees
    cl_extravar=np.exp(-larr*(larr+1)*(3.5*np.pi/180/2.355)**2)
    norm=0.67E-3/np.sum(cl_extravar*(2*larr+1)/4/np.pi); cl_extravar*=norm; #This yields <10% variations
    depth_factor=1+hp.synfast(cl_extravar,o.nside_out,new=True,verbose=False);

    #This total map corresponds to the relative variation in the shot-noise variance
    snvmap=shotnoise_factor*depth_factor*msk

    if o.plot_stuff :
        hp.mollview(snvmap,min=0.9,title='Relative local shot noise variance');
        plt.show()

    hp.write_map("cont_lss_nvar_ns%d.fits"%o.nside_out,snvmap,overwrite=True)

def match_hp_resolution(in_map, nside_out, UNSEEN2nan=True):
    """Utility to convert healpix map resolution if needed and change hp.UNSEEN values to
    np.nan.

    Parameters
    ----------
    in_map : np.array
        A valie healpix map
    nside_out : int
        The desired resolution to convert in_map to
    UNSEEN2nan : bool (True)
        If True, convert any hp.UNSEEN values to np.nan
    """
    current_nside = hp.npix2nside(np.size(in_map))
    if current_nside != nside_out:
        out_map = hp.ud_grade(in_map, nside_out=nside_out)
    else:
        out_map = in_map
    if UNSEEN2nan:
        out_map[np.where(out_map == hp.UNSEEN)] = np.nan
    return out_map

# Late imports

from lalinference.io import fits
from lalinference.bayestar import postprocess
import healpy as hp
import numpy as np
import json


# Read input file
prob, _ = fits.read_sky_map(opts.input.name, nest=True)

# Resample if requested
if opts.nside is not None and opts.interpolate in ('nearest', 'nested'):
    prob = hp.ud_grade(prob, opts.nside, order_in='NESTED', power=-2)
elif opts.nside is not None and opts.interpolate == 'bilinear':
    prob = postprocess.smooth_ud_grade(prob, opts.nside, nest=True)
if opts.interpolate == 'nested':
    prob = postprocess.interpolate_nested(prob, nest=True)

# Find credible levels
i = np.flipud(np.argsort(prob))
cumsum = np.cumsum(prob[i])
cls = np.empty_like(prob)
cls[i] = cumsum * 100

# Generate contours
paths = list(postprocess.contour(
    cls, opts.contour, nest=True, degrees=True, simplify=opts.simplify))

json.dump({

#Now, impose additional depth variations on scales of 3.5 degrees
    cl_extravar=np.exp(-larr*(larr+1)*(3.5*np.pi/180/2.355)**2)
    norm=0.67E-3/np.sum(cl_extravar*(2*larr+1)/4/np.pi); cl_extravar*=norm; #This yields <10% variations
    depth_factor=1+hp.synfast(cl_extravar,o.nside_out,new=True,verbose=False);

    #This total map corresponds to the relative variation in the shot-noise variance
    snvmap=shotnoise_factor*depth_factor*msk

    if o.plot_stuff :
        hp.mollview(snvmap,min=0.9,title='Relative local shot noise variance');
        plt.show()

    hp.write_map("cont_lss_nvar_ns%d.fits"%o.nside_out,snvmap,overwrite=True)

if not os.path.isfile("cont_lss_star_ns%d.fits"%o.nside_out) :
    starmap=hp.ud_grade(hp.read_map("star_template.fits",verbose=False),nside_out=o.nside_out)
    starmap[starmap<=0]=1E-15
    starmap=-np.log(starmap) #Contaminant will be proportional to log(n_star)
    mean=np.sum(starmap*msk)/np.sum(msk)

    clcont=hp.anafast((starmap-mean)*msk)/fsky
    nflat=np.mean(clcont[max(2,o.nside_out//3-50):o.nside_out//3+50])*np.ones_like(clcont); #Add extra fluctuations beyond ell~nside/3 with flat power spectrum
    dstar=hp.synfast(nflat,o.nside_out,new=True,verbose=False)
    starmap_b=np.maximum(starmap+dstar,0)
    clcont_b=hp.anafast((starmap_b-mean)*msk)/fsky
    ratio=-np.sqrt(0.03*np.sum(cltt[max(2,o.nside_out//3-50):o.nside_out//3+50])/np.sum(clcont_b[max(2,o.nside_out//3-50):o.nside_out//3+50])) #3% contamination at ell~nside/3
    starmap_b=(starmap_b-mean)*ratio

    dmap=hp.synfast(cltt[:3*o.nside_out],o.nside_out,new=True,verbose=False);

    if o.plot_stuff :
        cld=hp.anafast(dmap*msk)/fsky