Skip to content

mapmaking

Functions that wrap useful minkasi recipes

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

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

Arguments:

todvec: The TODs what we are mapmaking.

mapset: The mapset to compute priors with.
        We assume that the first element is the map we care about.

gradmap: Containter to use as the gradient map.

*args: Additional arguments to pass to get_grad_mask_2d.

**kwargs: Kewword arguments to pass to get_grad_mask_2d.

Returns:

new_mapset: A mapset with the original map and a cleared prior map.
Source code in witch/mapmaking.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def get_grad_prior(
    todvec: minkasi.tods.TodVec,
    mapset: minkasi.maps.Mapset,
    gradmap: minkasi.maps.MapType,
    *args,
    **kwargs,
) -> tuple[minkasi.mapmaking.pcg.HasPrior, minkasi.maps.Mapset]:
    """
    Make a gradient based prior. This helps avoid errors due to sharp features.

    Arguments:

        todvec: The TODs what we are mapmaking.

        mapset: The mapset to compute priors with.
                We assume that the first element is the map we care about.

        gradmap: Containter to use as the gradient map.

        *args: Additional arguments to pass to get_grad_mask_2d.

        **kwargs: Kewword arguments to pass to get_grad_mask_2d.

    Returns:

        new_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_naive(todvec, skymap, outdir)

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

Arguments:

todvec: The TODs to mapmake.

skymap: Map to use as footprint for outputs.

Returns:

naive: The navie map.

hits: The hit count map.
We use this as a preconditioner which helps small-scale convergence quite a bit.
Source code in witch/mapmaking.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
def make_naive(
    todvec: minkasi.tods.TodVec, skymap: minkasi.maps.MapType, outdir: str
) -> tuple[minkasi.maps.MapType, minkasi.maps.MapType]:
    """
    Make a naive map where we just bin common mode subtracted TODs.

    Arguments:

        todvec: The TODs to mapmake.

        skymap: Map to use as footprint for outputs.

    Returns:

        naive: The navie map.

        hits: The hit count map.
        We use this as a preconditioner which helps small-scale convergence quite a bit.
    """
    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()
    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"))

    return naive, hits

make_weights(todvec, skymap, outdir)

Make weights and noise map.

Arguments:

todvec: The TODs to mapmake.

skymap: Map to use as footprint for outputs.

Returns:

weightmap: The weights map.

noisemap: The noise map.
          This is just 1/sqrt(weights).
Source code in witch/mapmaking.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def make_weights(
    todvec: minkasi.tods.TodVec, skymap: minkasi.maps.MapType, outdir: str
) -> tuple[minkasi.maps.MapType, minkasi.maps.MapType]:
    """
    Make weights and noise map.

    Arguments:

        todvec: The TODs to mapmake.

        skymap: Map to use as footprint for outputs.

    Returns:

        weightmap: The weights map.

        noisemap: 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 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:

Arguments:

todvec: The TODs to reestimate noise for.

mapset: Mapset containing the current map solution.

noise_class: Which noise model to use.

noise_args: Additional arguments to pass to set_noise.

noise_kwargs: Additional keyword argmuents to pass to set_noise.
Source code in witch/mapmaking.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def reestimate_noise_from_map(
    todvec: minkasi.tods.TodVec,
    mapset: minkasi.maps.Mapset,
    noise_class: minkasi.mapmaking.NoiseModelType,
    noise_args: list,
    noise_kwargs: dict,
):
    """
    Use the current guess at the map to reestimate the noise:

    Arguments:

        todvec: The TODs to reestimate noise for.

        mapset: Mapset containing the current map solution.

        noise_class: Which noise model to use.

        noise_args: Additional arguments to pass to set_noise.

        noise_kwargs: Additional keyword argmuents to pass to 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)

Solve for map with PCG.

Arguments:

todvec: The TODs what we are mapmaking.

x0: The initial guess mapset.

ihits: The inverse hits map.

prior: Prior to use when mapmaking, set to None to not use.

maxiters: Maximum PCG iters to use.

save_iters: Which iterations to save the map at.

outdir: The output directory

desc_str: String used to deterime outroot.

Returns:

mapset: The mapset with the solved map.
Source code in witch/mapmaking.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
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,
) -> minkasi.maps.Mapset:
    """
    Solve for map with PCG.

    Arguments:

        todvec: The TODs what we are mapmaking.

        x0: The initial guess mapset.

        ihits: The inverse hits map.

        prior: Prior to use when mapmaking, set to None to not use.

        maxiters: Maximum PCG iters to use.

        save_iters: Which iterations to save the map at.

        outdir: The output directory

        desc_str: String used to deterime outroot.

    Returns:

        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)

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

    # 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,
            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