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")
print(input_map)
ProjectionMap:
shape(nu, y, x): (1, 1024, 1024)
stokes: naive
nu: [150.] GHz
t: naive
z: naive
quantity: rayleigh_jeans_temperature
units: uK_RJ
min: -3.830e+03
max: -4.984e+00
rms: 5.775e+02
center:
ra: 17ʰ20ᵐ0.00ˢ
dec: -10°00’0.00”
size(y, x): (1°, 1°)
resolution(y, x): (3.516”, 3.516”)
beam(maj, min, rot): (0 rad, 0 rad, 0 rad)
memory: 16.78 MB
[2]:
from maria.instrument import Band
f150 = Band(
center=150e9,
width=30e9,
NET_RJ=30e-6,
knee=1e1,
gain_error=2e-2)
array = {"field_of_view": 0.25,
"beam_spacing": 2,
"primary_size": 25,
"shape": "circle",
"bands": [f150]}
instrument = maria.get_instrument(array=array)
print(instrument)
instrument.plot()
Instrument(1 array)
├ arrays:
│ n FOV baseline bands polarized
│ array1 425 14.96’ 0 m [f150] False
│
└ bands:
name center width η NEP NET_RJ NET_CMB FWHM
0 f150 150 GHz 30 GHz 0.5 6.125 aW√s 30 uK_RJ√s 52.02 uK_CMB√s 21”
[3]:
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=3600,
max_chunk_duration=3600,
sample_rate=25,
scan_options={"radius": input_map.width.deg})
plans[0].plot()
print(plans)
PlanList(1 plans, 3600 s):
start_time duration target(ra,dec) center(az,el)
chunk
0 2026-03-17 09:11:15.000 +00:00 3600 s (260°, -10°) (31.48°, 75.07°)
[4]:
sim = maria.Simulation(
instrument,
plans=plans,
site="cerro_toco",
map=input_map,
atmosphere="2d",
atmosphere_kwargs={"weather": {"pwv": 1.0}},
)
print(sim)
Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │ n FOV baseline bands polarized
│ │ array1 425 14.96’ 0 m [f150] False
│ │
│ └ bands:
│ name center width η NEP NET_RJ NET_CMB FWHM
│ 0 f150 150 GHz 30 GHz 0.5 6.125 aW√s 30 uK_RJ√s 52.02 uK_CMB√s 21”
├ 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, 3600 s):
│ start_time duration target(ra,dec) center(az,el)
│ chunk
│ 0 2026-03-17 09:11:15.000 +00:00 3600 s (260°, -10°) (31.48°, 75.07°)
├ '2d'
└ ProjectionMap:
shape(stokes, nu, t, y, x): (1, 1, 1, 1024, 1024)
stokes: I
nu: [150.] GHz
t: [1.77403457e+09]
z: naive
quantity: rayleigh_jeans_temperature
units: uK_RJ
min: -3.830e+03
max: -4.984e+00
rms: 5.775e+02
center:
ra: 17ʰ20ᵐ0.00ˢ
dec: -10°00’0.00”
size(y, x): (1°, 1°)
resolution(y, x): (3.516”, 3.516”)
beam(maj, min, rot): (0 rad, 0 rad, 0 rad)
memory: 16.78 MB
[5]:
tods = sim.run()
tods[0].plot()
2026-03-20 19:23:06.689 INFO: Simulating observation 1 of 1
Constructing atmosphere: 100%|████████████████| 8/8 [00:50<00:00, 6.27s/it]
Generating turbulence: 100%|████████████████| 8/8 [00:04<00:00, 1.85it/s]
Sampling turbulence: 100%|████████████████| 8/8 [00:07<00:00, 1.10it/s]
Computing atmospheric emission: 100%|████████████████| 1/1 [00:01<00:00, 1.56s/it, band=f150]
Sampling map: 0%| | 0/1 [00:00<?, ?it/s, band=f150] /opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:4596: RuntimeWarning: invalid value encountered in scalar subtract
diff_b_a = b - a
Sampling map: 100%|████████████████| 1/1 [00:20<00:00, 20.00s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|████████████████| 1/1 [00:01<00:00, 1.46s/it, band=f150]
2026-03-20 19:24:54.504 INFO: Simulated observation 1 of 1 in 107.8 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 MaximumLikelihoodMapper
ml_mapper = MaximumLikelihoodMapper(tods=tods,
width=input_map.width.deg,
height=input_map.height.deg,
resolution=2 * input_map.resolution.deg,
units="mK_RJ",
k=3)
2026-03-20 19:25:04.991 INFO: Inferring center {'ra': '17ʰ19ᵐ59.74ˢ', 'dec': '-9°59’56.84”'} for mapper.
2026-03-20 19:25:04.996 INFO: Inferring mapper stokes parameters 'I' for mapper.
Preprocessing TODs: 100%|████████████████| 1/1 [00:02<00:00, 2.84s/it]
Mapping: 100%|██████████| 1/1 [00:04<00:00, 4.04s/it, tod=1/1]
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")
ProjectionMap:
shape(stokes, nu, t, y, x): (1, 1, 1, 512, 512)
stokes: I
nu: [150.] GHz
t: [1.77374047e+09]
z: naive
quantity: rayleigh_jeans_temperature
units: mK_RJ
min: -4.389e+00
max: 2.701e+00
rms: 8.171e-01
center:
ra: 17ʰ19ᵐ59.74ˢ
dec: -9°59’56.84”
size(y, x): (1°, 1°)
resolution(y, x): (7.031”, 7.031”)
beam(maj, min, rot): (21”, 21”, 0 rad)
memory: 4.194 MB
To fit the map solution, we run
[8]:
ml_mapper.fit(epochs=4, steps_per_epoch=64, lr=1e-1, plot=True, optimizer="adam")
Updating noise model: 100%|██████████| 1/1 [00:05<00:00, 5.71s/it, tod=1/1]
Fitting map (epoch 1/4): 100%|████████████████| 64/64 [04:11<00:00, 3.94s/it, loss=2.667e-02]
Updating noise model: 100%|██████████| 1/1 [00:04<00:00, 4.48s/it, tod=1/1]
Fitting map (epoch 2/4): 100%|████████████████| 64/64 [03:57<00:00, 3.71s/it, loss=8.969e-03]
Updating noise model: 100%|██████████| 1/1 [00:04<00:00, 4.20s/it, tod=1/1]
Fitting map (epoch 3/4): 100%|████████████████| 64/64 [03:57<00:00, 3.72s/it, loss=3.542e-03]
Updating noise model: 100%|██████████| 1/1 [00:04<00:00, 4.29s/it, tod=1/1]
Fitting map (epoch 4/4): 100%|████████████████| 64/64 [03:57<00:00, 3.70s/it, loss=1.094e-03]
that will improve as it continues to fit.