Maximum likelihood mapmaking¶
[1]:
import maria
from maria.io import fetch
input_map = maria.map.get("maps/cluster2.fits", nu=90e9)
input_map.data *= 2e2
input_map.plot(cmap="cmb")
print(input_map)
2026-06-05 13:41:58.847 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster2.fits
Downloading: 100%|██████████| 4.20M/4.20M [00:00<00:00, 103MB/s]
ProjectionMap:
data(1, 1024, 1024):
min: -7.691e-04
max: -1.001e-06
units: compton_y
quantity: compton_y
nu(1):
values: [90.] GHz
eta(1024):
height: 59.94’
res: -3.516”
xi(1024):
width: 59.94’
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: 4.194 MB
[2]:
from maria.instrument import Band
f090 = Band(
center=90e9,
width=30e9,
NET_RJ=20e-6,
knee=2e1)
array = {"name": "my_custom_array",
"field_of_view": 3/60,
"primary_size": 50,
"n": 91,
"shape": "circle",
# "packing": "sunflower",
"bands": [f090]}
instrument = maria.get_instrument(array=array)
print(instrument)
instrument.plot()
Instrument(1 array)
├ arrays:
│ n field_of_view max_baseline bands polarized primary_size
│ my_custom_array 91 3’ 0 m [f090] False 50 m
│
└ bands:
name center width η NEP NET_RJ NET_CMB FWHM
0 f090 90 GHz 30 GHz 0.5 4.407 aW√s 20 uK_RJ√s 24.6 uK_CMB√s 17.5”
[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=900,
max_chunk_duration=900,
sample_rate=50,
scan_pattern="daisy",
scan_options={
"radius": 0.75 * input_map.width.deg,
"speed": 0.5,
})
plans[0].plot()
print(plans)
PlanList(1 plans, 900 s):
start_time duration target(ra,dec) center(az,el)
chunk
0 2026-03-17 09:09:22.500 +00:00 900 s (260°, -9.999°) (48.21°, 71.41°)
[4]:
sim = maria.Simulation(
instrument,
plans=plans,
site="cerro_toco",
map=input_map,
atmosphere="2d",
atmosphere_kwargs={"weather": {"pwv": 1.}},
)
print(sim)
Constructing atmosphere: 100%|██████████| 8/8 [00:01<00:00, 6.42it/s]
Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │ n field_of_view max_baseline bands polarized primary_size
│ │ my_custom_array 91 3’ 0 m [f090] False 50 m
│ │
│ └ bands:
│ name center width η NEP NET_RJ NET_CMB FWHM
│ 0 f090 90 GHz 30 GHz 0.5 4.407 aW√s 20 uK_RJ√s 24.6 uK_CMB√s 17.5”
├ 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, 900 s):
│ start_time duration target(ra,dec) center(az,el)
│ chunk
│ 0 2026-03-17 09:09:22.500 +00:00 900 s (260°, -9.999°) (48.21°, 71.41°)
├ Atmosphere(8 processes with 8 layers):
│ ├ spectrum:
│ │ region: chajnantor
│ └ weather:
│ region: chajnantor
│ altitude: 5.19 km
│ time: Mar 17 06:16:52 -03:00
│ pwv[mean, rms]: (1 mm, 30 um)
└ ProjectionMap:
data(1, 1, 1, 1024, 1024):
min: -7.691e-04
max: -1.001e-06
units: compton_y
quantity: compton_y
stokes(1):
components: ['I']
nu(1):
values: [90.] GHz
t(1):
values: [1.78066692e+09] s
eta(1024):
height: 59.94’
res: -3.516”
xi(1024):
width: 59.94’
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: 4.194 MB
[5]:
tods = sim.run()
tods[0].plot()
2026-06-05 13:42:11.373 INFO: Simulating observation 1 of 1
Generating turbulence: 100%|██████████| 8/8 [00:00<00:00, 20.85it/s]
Sampling turbulence: 100%|██████████| 8/8 [00:04<00:00, 1.81it/s]
Computing atmospheric emission: 100%|██████████| 1/1 [00:00<00:00, 1.11it/s, band=f090]
Sampling map: 100%|██████████| 1/1 [00:03<00:00, 3.58s/it, band=f090, channel=(0 Hz, inf Hz)]
Generating noise: 100%|██████████| 1/1 [00:00<00:00, 1.20it/s, band=f090]
2026-06-05 13:42:24.659 INFO: Simulated observation 1 of 1 in 13.28 s
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,
tod_preprocessing={
"remove_spline": {"knot_spacing": 60, "remove_el_gradient_order": 3},
},
init="bin",
units="compton_y",
bilinear=False,
prior=False, k=1)
2026-06-05 13:42:28.156 INFO: Inferring resolution = 8.748” from detector FWHM
2026-06-05 13:42:29.195 INFO: Inferring center {'ra': '17ʰ19ᵐ59.97ˢ', 'dec': '-9°59’57.48”'} for mapper
2026-06-05 13:42:29.208 INFO: Inferring mapper width 1.536° for mapper from observation patch
2026-06-05 13:42:29.208 INFO: Inferring mapper height 1.536° to match supplied width
2026-06-05 13:42:29.262 INFO: Inferring stokes parameters 'I' for mapper from detector sensitivities
Preprocessing TODs: 100%|██████████| 1/1 [00:02<00:00, 2.27s/it]
Computing pointing matrices: 100%|██████████| 1/1 [00:01<00:00, 1.18s/it]
The initial solution is just a binning of the data, which has some noise artifacts and is missing some power (especially at large scales):
[7]:
from maria.mappers import compute_residual_map
print(ml_mapper.map)
ml_mapper.map.plot()
residual_map = compute_residual_map(input_map, ml_mapper.map)
residual_map.plot()
ProjectionMap:
data(1, 632, 632):
min: -3.509e-04
max: 3.443e-04
units: compton_y
quantity: compton_y
nu(1):
values: [90.] GHz
eta(632):
height: 1.533°
res: -8.748”
xi(632):
width: 1.533°
res: 8.748”
frame: ra/dec
center:
ra: 17ʰ19ᵐ59.97ˢ
dec: -9°59’57.48”
beam(maj, min, psi): (17.5”, 17.5”, 0 rad)
memory: 3.195 MB
To improve the map, we build a noise model perform conjugate gradient descent to solve the mapmaking equation \(m = (P^\top N^{-1} P)^{-1} P^\top N^{-1} d\) where \(m\) is the map, \(P\) is the pointing matrix, \(N = \langle n \otimes n \rangle\) is the noise covariance, \(d = Pm + n\) is the data, and \(n\) is the noise.
[8]:
ml_mapper.fit(epochs=2, max_steps_per_epoch=400, plot=True)
Updating noise model: 100%|██████████| 1/1 [00:01<00:00, 1.37s/it, tod=1/1]
Fitting map (epoch 1/2): 230it [02:57, 1.12it/s, alpha=70.4]2026-06-05 13:45:40.154 INFO: Finished conjugate gradient descent, terminating
Fitting map (epoch 1/2): 230it [02:58, 1.29it/s, alpha=70.4]
Updating noise model: 100%|██████████| 1/1 [00:01<00:00, 1.36s/it, tod=1/1]
Fitting map (epoch 2/2): 242it [03:01, 1.11it/s, alpha=49.2]2026-06-05 13:48:44.708 INFO: Finished conjugate gradient descent, terminating
Fitting map (epoch 2/2): 242it [03:02, 1.33it/s, alpha=49.2]
[9]:
from maria.mappers import compute_residual_map
residual_map = compute_residual_map(input_map, ml_mapper.map)
residual_map.plot()
