Maximum likelihood mapmaking¶
[1]:
import maria
from maria.io import fetch
input_map = maria.map.get("maps/cluster2.fits", nu=150e9)
input_map.data *= 2e2
input_map.plot(cmap="cmb", center_zero=True)
print(input_map)
2026-06-02 15:41:43.767 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster2.fits
Downloading: 100%|██████████| 4.20M/4.20M [00:00<00:00, 116MB/s]
ProjectionMap:
data(1, 1024, 1024):
min: -7.691e-04
max: -1.001e-06
units: compton_y
quantity: compton_y
nu(1):
values: [150.] 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
f150 = Band(
center=150e9,
width=30e9,
NET_RJ=10e-6,
knee=1e1)
array = {"field_of_view": 0.25,
"primary_size": 50,
"n": 100,
"shape": "circle",
"packing": "sunflower",
"bands": [f150]}
instrument = maria.get_instrument(array=array)
print(instrument)
instrument.plot()
Instrument(1 array)
├ arrays:
│ n FOV baseline bands polarized
│ array1 100 15’ 0 m [f150] False
│
└ bands:
name center width η NEP NET_RJ NET_CMB FWHM
0 f150 150 GHz 30 GHz 0.5 2.042 aW√s 10 uK_RJ√s 17.34 uK_CMB√s 10.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=1800,
max_chunk_duration=1800,
sample_rate=100,
scan_options={
"radius": input_map.width.deg,
"speed": 1,
})
plans[0].plot()
print(plans)
PlanList(1 plans, 1800 s):
start_time duration target(ra,dec) center(az,el)
chunk
0 2026-03-17 09:09:22.500 +00:00 1800 s (260°, -9.999°) (43.78°, 72.66°)
[4]:
sim = maria.Simulation(
instrument,
plans=plans,
site="cerro_toco",
map=input_map,
atmosphere="2d",
atmosphere_kwargs={"weather": {"pwv": 0.5}},
)
print(sim)
Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │ n FOV baseline bands polarized
│ │ array1 100 15’ 0 m [f150] False
│ │
│ └ bands:
│ name center width η NEP NET_RJ NET_CMB FWHM
│ 0 f150 150 GHz 30 GHz 0.5 2.042 aW√s 10 uK_RJ√s 17.34 uK_CMB√s 10.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, 1800 s):
│ start_time duration target(ra,dec) center(az,el)
│ chunk
│ 0 2026-03-17 09:09:22.500 +00:00 1800 s (260°, -9.999°) (43.78°, 72.66°)
├ '2d'
└ 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: [150.] GHz
t(1):
values: [1.7804149e+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-02 15:41:53.311 INFO: Simulating observation 1 of 1
Constructing atmosphere: 100%|██████████| 8/8 [00:05<00:00, 1.58it/s]
Generating turbulence: 100%|██████████| 8/8 [00:00<00:00, 11.09it/s]
Sampling turbulence: 100%|██████████| 8/8 [00:04<00:00, 1.87it/s]
Computing atmospheric emission: 100%|██████████| 1/1 [00:00<00:00, 1.01it/s, band=f150]
Sampling map: 100%|██████████| 1/1 [00:07<00:00, 7.43s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|██████████| 1/1 [00:01<00:00, 1.61s/it, band=f150]
2026-06-02 15:42:24.601 INFO: Simulated observation 1 of 1 in 31.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",
bilinear=False,
prior=False)
2026-06-02 15:42:29.128 INFO: Inferring resolution = 5.249” from detector FWHM
2026-06-02 15:42:31.400 INFO: Inferring center {'ra': '17ʰ20ᵐ0.42ˢ', 'dec': '-9°59’55.45”'} for mapper
2026-06-02 15:42:31.409 INFO: Inferring mapper width 2.215° for mapper from observation patch
2026-06-02 15:42:31.410 INFO: Inferring mapper height 2.215° to match supplied width
2026-06-02 15:42:31.412 INFO: Inferring stokes parameters 'I' for mapper from detector sensitivities
Preprocessing TODs: 100%|██████████| 1/1 [00:06<00:00, 6.99s/it]
Computing pointing matrices: 100%|██████████| 1/1 [00:03<00:00, 3.66s/it]
Binning the map yields only large-scale modes with no discernible signal:
To improve the solution, we fit the mapmaker by iteratively solving the mapmaking equation. The initial solution is a field of random noise:
[7]:
print(ml_mapper.map)
ml_mapper.map.plot()
ProjectionMap:
data(1, 1519, 1519):
min: -6.013e-04
max: 9.314e-04
units: K_RJ
quantity: rayleigh_jeans_temperature
nu(1):
values: [150.] GHz
eta(1519):
height: 2.213°
res: -5.249”
xi(1519):
width: 2.213°
res: 5.249”
frame: ra/dec
center:
ra: 17ʰ20ᵐ0.42ˢ
dec: -9°59’55.45”
beam(maj, min, psi): (10.5”, 10.5”, 0 rad)
memory: 18.46 MB
[8]:
ml_mapper.fit(epochs=4,
steps_per_epoch=25,
plot=True)
Updating noise model: 100%|██████████| 1/1 [00:04<00:00, 4.64s/it, tod=1/1]
Fitting map (epoch 1/4): 100%|██████████| 25/25 [01:21<00:00, 3.26s/it, mll=-1.099e+03, pgrad=-0.228, step=4.574e-09]
Updating noise model: 100%|██████████| 1/1 [00:04<00:00, 4.61s/it, tod=1/1]
Fitting map (epoch 2/4): 100%|██████████| 25/25 [01:21<00:00, 3.26s/it, mll=-1.065e+03, pgrad=0.930, step=1.236e-09]
Updating noise model: 100%|██████████| 1/1 [00:04<00:00, 4.62s/it, tod=1/1]
Fitting map (epoch 3/4): 100%|██████████| 25/25 [01:21<00:00, 3.26s/it, mll=-5.925e+02, pgrad=0.949, step=6.093e-09]
Updating noise model: 100%|██████████| 1/1 [00:04<00:00, 4.62s/it, tod=1/1]
Fitting map (epoch 4/4): 100%|██████████| 25/25 [01:21<00:00, 3.26s/it, mll=-4.661e+02, pgrad=-0.142, step=5.888e-09]
[9]:
from maria.mappers import compute_residual_map
residual_map = compute_residual_map(input_map.smooth(fwhm=instrument.dets.fwhm.mean()),
ml_mapper.map)
residual_map.plot()
