Maximum likelihood mapmaking

[1]:
import maria
from maria.io import fetch

input_map = maria.map.get("maps/cluster1.fits", nu=150e9).to("uK_RJ")
input_map.data *= 2e2

input_map.plot(cmap="cmb", center_zero=True)
print(input_map)
ProjectionMap:
  data(1, 1024, 1024):
    min: -3.830e+03
    max: -4.984e+00
    units: uK_RJ
    quantity: rayleigh_jeans_temperature
  nu(1):
    values: [150.] GHz
  y(1024):
    height: 1°
    res: 3.516”
  x(1024):
    width: 1°
    res: 3.516”
  frame: ra/dec
  center:
    ra: 17ʰ20ᵐ0.00ˢ
    dec: -10°00’0.00”
  beam(maj, min, psi): (0 rad, 0 rad, 0 rad)
  memory: 8.389 MB
../_images/tutorials_maximum-likelihood-mapper_1_1.png
[2]:
from maria.instrument import Band

f150 = Band(
    center=150e9,
    width=30e9,
    NET_RJ=50e-6,
    knee=1e1,
    gain_error=2e-2)

array = {"field_of_view": 0.5,
         "beam_spacing": 1.8,
         "primary_size": 12,
         "packing": "sunflower",
         "shape": "circle",
         "bands": [f150]}

instrument = maria.get_instrument(array=array)

print(instrument)
instrument.plot()
Instrument(1 array)
├ arrays:
│            n     FOV baseline   bands polarized
│  array1  372  30.01’      0 m  [f150]     False
│
└ bands:
      name   center   width    η         NEP      NET_RJ        NET_CMB    FWHM
   0  f150  150 GHz  30 GHz  0.5  10.21 aW√s  50 uK_RJ√s  86.7 uK_CMB√s  43.74”
../_images/tutorials_maximum-likelihood-mapper_2_1.png
[3]:
import numpy as np
from maria import Planner

planner = Planner(start_time="2026-03-16T12:00:00",
                  target=input_map,
                  site="cerro_toco",
                  constraints={"el": (70, 90)})

plans = planner.generate_plans(
                               total_duration=1200,
                               max_chunk_duration=1200,
                               sample_rate=25,
                               scan_options={
                                   "radius": 0.5 * input_map.width.deg,
                                   "miss_factor": 0.25,
                                   "miss_freq": 0.5,
                               })

plans[0].plot()
print(plans)
PlanList(1 plans, 1200 s):
                           start_time duration target(ra,dec)     center(az,el)
chunk
0      2026-03-17 09:10:00.000 +00:00   1200 s   (260°, -10°)  (46.43°, 71.95°)
../_images/tutorials_maximum-likelihood-mapper_3_1.png
[4]:
sim = maria.Simulation(
    instrument,
    plans=plans,
    site="cerro_toco",
    map=input_map,
    atmosphere="2d",
    atmosphere_kwargs={"weather": {"pwv": 1}},
)

print(sim)
Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │            n     FOV baseline   bands polarized
│ │  array1  372  30.01’      0 m  [f150]     False
│ │
│ └ bands:
│       name   center   width    η         NEP      NET_RJ        NET_CMB    FWHM
│    0  f150  150 GHz  30 GHz  0.5  10.21 aW√s  50 uK_RJ√s  86.7 uK_CMB√s  43.74”
├ Site:
│   region: chajnantor
│   timezone: America/Santiago
│   location:
│     longitude: 67°47’16.08” W
│     latitude:  22°57’30.96” S
│     altitude: 5.19 km
│   seasonal: True
│   diurnal: True
├ PlanList(1 plans, 1200 s):
│                            start_time duration target(ra,dec)     center(az,el)
│ chunk
│ 0      2026-03-17 09:10:00.000 +00:00   1200 s   (260°, -10°)  (46.43°, 71.95°)
├ '2d'
└ ProjectionMap:
    data(1, 1, 1, 1024, 1024):
      min: -3.830e+03
      max: -4.984e+00
      units: uK_RJ
      quantity: rayleigh_jeans_temperature
    stokes(1):
      components: I
    nu(1):
      values: [150.] GHz
    t(1):
      values: [1.77867254e+09] s
    y(1024):
      height: 1°
      res: 3.516”
    x(1024):
      width: 1°
      res: 3.516”
    frame: ra/dec
    center:
      ra: 17ʰ20ᵐ0.00ˢ
      dec: -10°00’0.00”
    beam(maj, min, psi): (0 rad, 0 rad, 0 rad)
    memory: 8.389 MB
[5]:
tods = sim.run()
tods[0].plot()
2026-05-13 11:42:31.626 INFO: Simulating observation 1 of 1
Constructing atmosphere: 100%|██████████| 8/8 [00:37<00:00,  4.70s/it]
Generating turbulence: 100%|██████████| 8/8 [00:01<00:00,  5.37it/s]
Sampling turbulence: 100%|██████████| 8/8 [00:06<00:00,  1.32it/s]
Computing atmospheric emission: 100%|██████████| 1/1 [00:01<00:00,  1.33s/it, band=f150]
Sampling map: 100%|██████████| 1/1 [00:06<00:00,  6.69s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|██████████| 1/1 [00:00<00:00,  1.89it/s, band=f150]
2026-05-13 11:43:33.854 INFO: Simulated observation 1 of 1 in 62.21 s
../_images/tutorials_maximum-likelihood-mapper_5_1.png

We can map the TOD with the MaximumLikelihoodMapper. Here the k parameter controls how many detector-detector eigenmodes to use in the noise model.

[6]:
from maria.mappers.ml_mapper import *

ml_mapper = MaximumLikelihoodMapper(tods=tods,
                                    resolution=4 * input_map.resolution,
                                    tod_preprocessing={
                                        "remove_spline": {"knot_spacing": 30, "remove_el_gradient_order": 3},
                                    },
                                    bilinear=True,
                                    k=3)
2026-05-13 11:43:39.385 INFO: Inferring center {'ra': '17ʰ20ᵐ0.03ˢ', 'dec': '-9°59’55.45”'} for mapper.
2026-05-13 11:43:39.401 INFO: Inferring mapper width 1.466° for mapper from observation patch.
2026-05-13 11:43:39.401 INFO: Inferring mapper height 1.466° to match supplied width.
2026-05-13 11:43:39.405 INFO: Inferring mapper stokes parameters 'I' for mapper.
Preprocessing TODs: 100%|██████████| 1/1 [00:04<00:00,  4.70s/it]
Computing pointing matrices: 100%|██████████| 1/1 [00:06<00:00,  6.43s/it]

The initial map is a “guess” constructing by lightly filtering and binning the input TODs:

[7]:
print(ml_mapper.map)
ml_mapper.map.plot(cmap="cmb", contrast=1e-4)
ProjectionMap:
  data(1, 1, 1, 375, 375):
    min: -2.363e-03
    max: 1.440e-03
    units: K_RJ
    quantity: rayleigh_jeans_temperature
  stokes(1):
    components: I
  nu(1):
    values: [150.] GHz
  t(1):
    values: [1.7737392e+09] s
  y(375):
    height: 1.465°
    res: 14.06”
  x(375):
    width: 1.465°
    res: 14.06”
  frame: ra/dec
  center:
    ra: 17ʰ20ᵐ0.03ˢ
    dec: -9°59’55.45”
  beam(maj, min, psi): (43.74”, 43.74”, 0 rad)
  memory: 1.125 MB
../_images/tutorials_maximum-likelihood-mapper_9_1.png
[8]:
ml_mapper.fit(epochs=8,
              steps_per_epoch=25,
              plot=True,
              plot_kwargs={"cmap": "cmb",
                           "contrast": 1e-4})
Updating noise model: 100%|██████████| 1/1 [00:05<00:00,  5.35s/it, tod=1/1]
Fitting map (epoch 1/8): 100%|██████████| 25/25 [02:50<00:00,  6.83s/it, mll=10.038, pgrad=0.710, step=2.43e-06]
../_images/tutorials_maximum-likelihood-mapper_10_1.png
Updating noise model: 100%|██████████| 1/1 [00:05<00:00,  5.36s/it, tod=1/1]
Fitting map (epoch 2/8): 100%|██████████| 25/25 [02:47<00:00,  6.72s/it, mll=9.303, pgrad=0.996, step=1.18e-05]
../_images/tutorials_maximum-likelihood-mapper_10_3.png
Updating noise model: 100%|██████████| 1/1 [00:05<00:00,  5.36s/it, tod=1/1]
Fitting map (epoch 3/8): 100%|██████████| 25/25 [02:50<00:00,  6.83s/it, mll=8.946, pgrad=0.998, step=1.04e-05]
../_images/tutorials_maximum-likelihood-mapper_10_5.png
Updating noise model: 100%|██████████| 1/1 [00:05<00:00,  5.35s/it, tod=1/1]
Fitting map (epoch 4/8): 100%|██████████| 25/25 [02:56<00:00,  7.08s/it, mll=8.683, pgrad=0.208, step=5.95e-07]
../_images/tutorials_maximum-likelihood-mapper_10_7.png
Updating noise model: 100%|██████████| 1/1 [00:05<00:00,  5.36s/it, tod=1/1]
Fitting map (epoch 5/8): 100%|██████████| 25/25 [02:50<00:00,  6.84s/it, mll=8.549, pgrad=-0.390, step=8.05e-07]
../_images/tutorials_maximum-likelihood-mapper_10_9.png
Updating noise model: 100%|██████████| 1/1 [00:05<00:00,  5.36s/it, tod=1/1]
Fitting map (epoch 6/8): 100%|██████████| 25/25 [02:57<00:00,  7.11s/it, mll=8.430, pgrad=1.000, step=1.62e-06]
../_images/tutorials_maximum-likelihood-mapper_10_11.png
Updating noise model: 100%|██████████| 1/1 [00:05<00:00,  5.34s/it, tod=1/1]
Fitting map (epoch 7/8): 100%|██████████| 25/25 [02:57<00:00,  7.08s/it, mll=8.300, pgrad=0.999, step=3.64e-07]
../_images/tutorials_maximum-likelihood-mapper_10_13.png
Updating noise model: 100%|██████████| 1/1 [00:05<00:00,  5.34s/it, tod=1/1]
Fitting map (epoch 8/8): 100%|██████████| 25/25 [02:50<00:00,  6.82s/it, mll=8.246, pgrad=0.988, step=7.80e-06]
../_images/tutorials_maximum-likelihood-mapper_10_15.png

that will improve as it continues to fit.