mapmaking

Functions that wrap useful minkasi recipes

apply_mask(input_mapset, maskset)

Apply a mask to a set of input maps.

Parameters:

Name	Type	Description	Default
input_mapset	Mapset	Maps to apply the mask to.	required
maskset	Mapset	Masks to apply; defined per map	required

Returns:

Name	Type	Description
target_mapset	Mapset	Masked version of input mapset.

Source code in witch/external/minkasi/mapmaking.py

def apply_mask(input_mapset: minkasi.maps.Mapset, maskset: minkasi.maps.Mapset):
    """
    Apply a mask to a set of input maps.

    Parameters
    ----------
    input_mapset : minkasi.maps.Mapset
        Maps to apply the mask to.
    maskset : minkasi.maps.Mapset
        Masks to apply; defined per map

    Returns
    -------
    target_mapset : minkasi.maps.Mapset
        Masked version of input mapset.

    """
    for i, imap in enumerate(input_mapset.maps):
        imap.map *= maskset.maps[i].map

    return input_mapset

get_grad_prior(todvec, mapset, gradmap, *args, **kwargs)

Make a gradient based prior from a map. This helps avoid errors due to sharp features.

Parameters:

Name	Type	Description	Default
todvec	TodVec	The TODs what we are mapmaking.	required
mapset	Mapset	The mapset to compute priors with. We assume that the first element is the map we care about.	required
gradmap	MapType	Containter to use as the gradient map.	required
*args	Unpack[tuple]	Additional arguments to pass to get_grad_mask_2d.	()
**kwargs		Keyword arguments to pass to get_grad_mask_2d.	{}

Returns:

Name	Type	Description
new_mapset	Mapset	A mapset with the original map and a cleared prior map.

Source code in witch/external/minkasi/mapmaking.py

def get_grad_prior(
    todvec: minkasi.tods.TodVec,
    mapset: minkasi.maps.Mapset,
    gradmap: minkasi.maps.MapType,
    *args: Unpack[tuple],
    **kwargs,
) -> minkasi.maps.Mapset:
    """
    Make a gradient based prior from a map.
    This helps avoid errors due to sharp features.

    Parameters
    ----------
    todvec : minkasi.tods.TodVec
        The TODs what we are mapmaking.
    mapset : minkasi.maps.Mapset
        The mapset to compute priors with.
        We assume that the first element is the map we care about.
    gradmap : minkasi.maps.MapType
        Containter to use as the gradient map.
    *args : Unpack[tuple]
        Additional arguments to pass to get_grad_mask_2d.
    **kwargs
        Keyword arguments to pass to get_grad_mask_2d.

    Returns
    -------
    new_mapset : minkasi.maps.Mapset
        A mapset with the original map and a cleared prior map.
    """
    gradmap.map[:] = minkasi.mapmaking.noise.get_grad_mask_2d(
        mapset.maps[0], todvec, *args, **kwargs
    )
    prior = minkasi.mapmaking.timestream.tsModel(todvec, minkasi.tods.cuts.CutsCompact)
    for tod in todvec.tods:
        prior.data[tod.info["fname"]] = tod.prior_from_skymap(gradmap)
        print(
            "prior on tod "
            + tod.info["fname"]
            + " length is "
            + repr(prior.data[tod.info["fname"]].map.size)
        )

    new_mapset = minkasi.maps.Mapset()
    new_mapset.add_map(mapset.maps[0])
    pp = prior.copy()
    pp.clear()
    new_mapset.add_map(pp)

    return new_mapset

make_maps(todvec, skymap, noise_class, noise_args, noise_kwargs, outdir, npass, dograd, thresh=None, return_maps=False)

Make a minkasi map with multple passes and noise reestimation. Unless you are an expert this will usually be all you need.

Parameters:

Name	Type	Description	Default
todvec	TodVec	The tods to mapmake.	required
skymap	MapType	Map to use as a template. The contents don't matter only the shape and WCS info.	required
noise_class	NoiseModelType	The noise model to use on the TODs.	required
noise_args	tuple	Arguments to pass to minkasi.tods.Tod.set_noise.	required
noise_kwargs	dict	Keyword arguments to pass to minkasi.tods.Tod.set_noise.	required
outdir	str	The output directory.	required
npass	int	The number of times to mapmake and then reestimate the noise.	required
dograd	bool	If True make a map based prior to avoid biases from sharp features.	required
thresh	int	If set, mask pixels with value less than thresh durring map-making. Default; None	None
return_maps	bool	If True, return the mapset. Default; False	False

Source code in witch/external/minkasi/mapmaking.py

def make_maps(
    todvec: minkasi.tods.TodVec,
    skymap: minkasi.maps.MapType,
    noise_class: minkasi.mapmaking.noise.NoiseModelType,
    noise_args: tuple,
    noise_kwargs: dict,
    outdir: str,
    npass: int,
    dograd: bool,
    thresh: Optional[int] = None,
    return_maps: bool = False,
):
    """
    Make a minkasi map with multple passes and noise reestimation.
    Unless you are an expert this will usually be all you need.

    Parameters
    ----------
    todvec : minkasi.tods.TodVec
        The tods to mapmake.
    skymap : minkasi.maps.MapType
        Map to use as a template.
        The contents don't matter only the shape and WCS info.
    noise_class : minkasi.mapmaking.noise.NoiseModelType
        The noise model to use on the TODs.
    noise_args : tuple
        Arguments to pass to `minkasi.tods.Tod.set_noise`.
    noise_kwargs : dict
        Keyword arguments to pass to `minkasi.tods.Tod.set_noise`.
    outdir : str
        The output directory.
    npass : int
        The number of times to mapmake and then reestimate the noise.
    dograd : bool
        If True make a map based prior to avoid biases from sharp features.
    thresh : int
        If set, mask pixels with value less than thresh durring map-making. Default; None
    return_maps : bool
        If True, return the mapset. Default; False
    """
    for tod in todvec.tods:
        if "ipix" not in tod.info:
            ipix = skymap.get_pix(tod)
            tod.info["ipix"] = ipix
        if tod.noise is None:
            tod.set_noise(
                noise_class,
                tod.info["dat_calib"],
                *noise_args,
                **noise_kwargs,
            )

    naive, hits, maskset = make_naive(
        todvec=todvec, skymap=skymap, outdir=outdir, thresh=thresh
    )

    # Take 1 over hits map
    ihits = hits.copy()
    ihits.invert()

    if thresh is not None:
        ihits.map *= maskset.maps[0].map

    # Save weights and noise maps
    _ = make_weights(todvec=todvec, skymap=skymap, outdir=outdir, maskset=maskset)

    # Setup the mapset
    # For now just include the naive map so we can use it as the initial guess.
    mapset = minkasi.maps.Mapset()
    mapset.add_map(naive)

    # run PCG to solve for a first guess
    iters = [5, 25, 100]
    mapset = solve_map(
        todvec=todvec,
        x0=mapset,
        ihits=ihits,
        prior=None,
        maxiters=26,
        save_iters=iters,
        outdir=outdir,
        desc_str="initial",
        maskset=maskset,
    )

    # Now we iteratively solve and reestimate the noise
    for niter in range(npass):
        maxiter = 26 + 25 * (niter + 1)
        reestimate_noise_from_map(todvec, mapset, noise_class, noise_args, noise_kwargs)

        # Make a gradient based prior
        if dograd:
            mapset = get_grad_prior(todvec, mapset, hits.copy(), thresh=1.8)
        # Solve
        mapset = solve_map(
            todvec,
            mapset,
            ihits,
            None,
            maxiter,
            iters,
            outdir,
            f"niter_{niter+1}",
            maskset=maskset,
        )

    minkasi.barrier()

    if return_maps:
        return mapset

make_naive(todvec, skymap, outdir, thresh=None)

Make a naive map where we just bin common mode subtracted TODs.

Parameters:

Name	Type	Description	Default
todvec	TodVec	The TODs to mapmake.	required
skymap	MapType	Map to use as footprint for outputs.	required
thresh	int	If set, mask pixels with value less than thresh durring map-making. Default; None	None

Returns:

Name	Type	Description
naive	MapType	The naive map.
hits	MapType	The hit count map. We use this as a preconditioner which helps small-scale convergence quite a bit.
maskset	Mapset	Set of masks for the maps.

Source code in witch/external/minkasi/mapmaking.py

def make_naive(
    todvec: minkasi.tods.TodVec,
    skymap: minkasi.maps.MapType,
    outdir: str,
    thresh: Optional[int] = None,
) -> tuple[minkasi.maps.MapType, minkasi.maps.MapType, minkasi.maps.Mapset]:
    """
    Make a naive map where we just bin common mode subtracted TODs.

    Parameters
    ----------
    todvec : minkasi.tods.TodVec
        The TODs to mapmake.
    skymap : minkasi.maps.MapType
        Map to use as footprint for outputs.
    thresh : int
        If set, mask pixels with value less than thresh durring map-making. Default; None

    Returns
    -------
    naive : minkasi.maps.MapType
        The naive map.
    hits : minkasi.maps.MapType
        The hit count map.
        We use this as a preconditioner which helps small-scale convergence quite a bit.
    maskset : minkasi.maps.Mapset
        Set of masks for the maps.
    """
    hits = minkasi.mapmaking.make_hits(todvec, skymap)

    # Make a naive map where we just bin the CM subbed tods
    naive = skymap.copy()
    naive.clear()
    for tod in todvec.tods:
        tmp = tod.info["dat_calib"].copy()
        u, s, v = np.linalg.svd(tmp, False)
        tmp -= np.outer(u[:, 0], s[0] * v[0, :])
        naive.tod2map(tod, tmp)
    naive.mpi_reduce()
    maskset = None
    if thresh is not None:
        mask = hits.copy()
        mask.map = 1.0 * hits.map > thresh
        maskset = minkasi.maps.Mapset()
        maskset.add_map(mask)
        hits.map *= mask.map
        naive.map *= mask.map

    naive.map[hits.map > 0] = naive.map[hits.map > 0] / hits.map[hits.map > 0]
    if minkasi.myrank == 0:
        naive.write(os.path.join(outdir, "naive.fits"))
        hits.write(os.path.join(outdir, "hits.fits"))
    naive.clear()
    return naive, hits, maskset

make_weights(todvec, skymap, outdir, maskset=None)

Make weights and noise map.

Parameters:

Name	Type	Description	Default
todvec	TodVec	The TODs to mapmake.	required
skymap	MapType	Map to use as footprint for outputs.	required
maskset	Mapset	Set of masks for the maps. Default; None	None

Returns:

Name	Type	Description
weightmap	MapType	The weights map.
noisemap	MapType	The noise map. This is just 1/sqrt(weights).

Source code in witch/external/minkasi/mapmaking.py

def make_weights(
    todvec: minkasi.tods.TodVec,
    skymap: minkasi.maps.MapType,
    outdir: str,
    maskset: Optional[minkasi.maps.Mapset] = None,
) -> tuple[minkasi.maps.MapType, minkasi.maps.MapType]:
    """
    Make weights and noise map.

    Parameters
    ----------
    todvec : minkasi.tods.TodVec
        The TODs to mapmake.
    skymap : minkasi.maps.MapType
        Map to use as footprint for outputs.
    maskset : minkasi.maps.Mapset
        Set of masks for the maps. Default; None

    Returns
    -------
    weightmap : minkasi.maps.MapType
        The weights map.
    noisemap : minkasi.maps.MapType
        The noise map.
        This is just 1/sqrt(weights).
    """
    weightmap = minkasi.mapmaking.make_hits(todvec, skymap, do_weights=True)
    mask = weightmap.map > 0
    tmp = weightmap.map.copy()
    tmp[mask] = 1.0 / np.sqrt(tmp[mask])
    noisemap = weightmap.copy()
    noisemap.map[:] = tmp

    if maskset is not None:
        noisemap.map *= maskset.maps[0].map
        weightmap.map *= maskset.maps[0].map
    if minkasi.myrank == 0:
        noisemap.write(os.path.join(outdir, "noise.fits"))
        weightmap.write(os.path.join(outdir, "weights.fits"))

    return weightmap, noisemap

reestimate_noise_from_map(todvec, mapset, noise_class, noise_args, noise_kwargs)

Use the current guess at the map to reestimate the noise:

Parameters:

Name	Type	Description	Default
todvec	TodVec	The TODs to reestimate noise for.	required
mapset	Mapset	Mapset containing the current map solution.	required
noise_class	NoiseModelType	Which noise model to use.	required
noise_args	tuple	Additional arguments to pass to minkasi.tods.Tod.set_noise.	required
noise_kwargs	dict	Additional keyword arguments to pass to minkasi.tods.Tod.set_noise.	required

Source code in witch/external/minkasi/mapmaking.py

def reestimate_noise_from_map(
    todvec: minkasi.tods.TodVec,
    mapset: minkasi.maps.Mapset,
    noise_class: minkasi.mapmaking.NoiseModelType,
    noise_args: tuple,
    noise_kwargs: dict,
):
    """
    Use the current guess at the map to reestimate the noise:

    Parameters
    ----------
    todvec : minkasi.tods.TodVec
        The TODs to reestimate noise for.
    mapset : minkasi.maps.Mapset
        Mapset containing the current map solution.
    noise_class : minkasi.mapmaking.NoiseModelType
        Which noise model to use.
    noise_args : tuple
        Additional arguments to pass to `minkasi.tods.Tod.set_noise`.
    noise_kwargs : dict
        Additional keyword arguments to pass to `minkasi.tods.Tod.set_noise`.
    """
    for tod in todvec.tods:
        mat = 0 * tod.info["dat_calib"]
        for mm in mapset.maps:
            mm.map2tod(tod, mat)
        tod.set_noise(
            noise_class,
            dat=tod.info["dat_calib"] - mat,
            *noise_args,
            **noise_kwargs,
        )

solve_map(todvec, x0, ihits, prior, maxiters, save_iters, outdir, desc_str, maskset=None)

Solve for map with PCG.

Parameters:

Name	Type	Description	Default
todvec	TodVec	The TODs what we are mapmaking.	required
x0	Mapset	The initial guess mapset.	required
ihits	MapType	The inverse hits map.	required
prior	Optional[HasPrior]	Prior to use when mapmaking, set to None to not use.	required
maxiters	int	Maximum PCG iters to use.	required
save_iters	list[int]	Which iterations to save the map at.	required
outdir	str	The output directory	required
desc_str	str	String used to determine outroot.	required
maskset	Mapset	Set of masks for the maps. Default; None	None

Returns:

Name	Type	Description
mapset	Mapset	The mapset with the solved map.

Source code in witch/external/minkasi/mapmaking.py

def solve_map(
    todvec: minkasi.tods.TodVec,
    x0: minkasi.maps.Mapset,
    ihits: minkasi.maps.MapType,
    prior: Optional[minkasi.mapmaking.pcg.HasPrior],
    maxiters: int,
    save_iters: list[int],
    outdir: str,
    desc_str: str,
    maskset: Optional[minkasi.maps.Mapset] = None,
) -> minkasi.maps.Mapset:
    """
    Solve for map with PCG.

    Parameters
    ----------
    todvec : minkasi.tods.TodVec
        The TODs what we are mapmaking.
    x0 : minkasi.maps.Mapset
        The initial guess mapset.
    ihits : minkasi.maps.MapType
        The inverse hits map.
    prior : Optional[minkasi.mapmaking.pgc.HasPrior]
        Prior to use when mapmaking, set to None to not use.
    maxiters : int
        Maximum PCG iters to use.
    save_iters : list[int]
        Which iterations to save the map at.
    outdir : str
        The output directory
    desc_str : str
        String used to determine outroot.
    maskset : minkasi.maps.Mapset
        Set of masks for the maps. Default; None

    Returns
    -------
    mapset : minkasi.maps.Mapset
        The mapset with the solved map.
    """
    # make A^T N^1 d.  TODs need to understand what to do with maps
    # but maps don't necessarily need to understand what to do with TODs,
    # hence putting make_rhs in the vector of TODs.
    # Again, make_rhs is MPI-aware, so this should do the right thing
    # if you run with many processes.
    rhs = x0.copy()
    todvec.make_rhs(rhs)

    if maskset is not None:
        rhs = apply_mask(input_mapset=rhs, maskset=maskset)

    # Preconditioner is 1/ hit count map.
    # Helps a lot for convergence.
    precon = x0.copy()
    precon.maps[0].map[:] = ihits.map[:]

    if (
        maskset is not None
    ):  # Note in the context of make_maps this is redundant but in principle
        # someone could use this function independantly so we have to enforce
        # consistent masking.
        precon = apply_mask(input_mapset=precon, maskset=maskset)

    # run PCG to solve
    # Supressing print here, probably want a verbosity setting on the minkasi side...
    with open(os.devnull, "w") as f, contextlib.redirect_stdout(f):
        mapset = minkasi.mapmaking.run_pcg_wprior(
            rhs,
            x0,
            todvec,
            prior,
            precon,
            mask=maskset,
            maxiter=maxiters,
            outroot=os.path.join(outdir, desc_str),
            save_iters=save_iters,
        )

    if minkasi.myrank == 0:
        mapset.maps[0].write(
            os.path.join(outdir, f"{desc_str}.fits")
        )  # and write out the map as a FITS file

    return mapset