Maximum likelihood mapmaking¶
[1]:
import maria
from maria.io import fetch
input_map = maria.map.load(fetch("maps/cluster1.fits"), nu=150e9).to("uK_RJ")
input_map.data *= 2e1
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+02
max: -4.984e-01
rms: 5.775e+01
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. 0. 0.]] rad
memory: 16.78 MB
[2]:
from maria import Planner
planner = Planner(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 / 2})
plans[0].plot()
print(plans)
PlanList(1 plans, 3600 s):
start_time duration target(ra,dec) center(az,el)
chunk
0 2026-01-22 12:42:33.862 +00:00 3600 s (260°, -10°) (32.31°, 74.95°)
[3]:
# import maria
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.15,
"beam_spacing": 1.25,
"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 380 8.878’ 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 0.3499’
[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 380 8.878’ 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 0.3499’
├ Site:
│ region: chajnantor
│ timezone: America/Santiago
│ location:
│ longitude: 67°47’16.08” W
│ latitude: 22°57’30.96” S
│ altitude: 5190 m
│ seasonal: True
│ diurnal: True
├ PlanList(1 plans, 3600 s):
│ start_time duration target(ra,dec) center(az,el)
│ chunk
│ 0 2026-01-22 12:42:33.862 +00:00 3600 s (260°, -10°) (32.31°, 74.95°)
├ '2d'
└ ProjectionMap:
shape(stokes, nu, t, y, x): (1, 1, 1, 1024, 1024)
stokes: I
nu: [150.] GHz
t: [1.76901915e+09]
z: naive
quantity: rayleigh_jeans_temperature
units: uK_RJ
min: -3.830e+02
max: -4.984e-01
rms: 5.775e+01
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. 0. 0.]]]] rad
memory: 16.78 MB
[5]:
tods = sim.run()
tods[0].plot()
2026-01-21 18:12:42.959 INFO: Simulating observation 1 of 1
Constructing atmosphere: 100%|████████████████| 8/8 [00:31<00:00, 3.92s/it]
Generating turbulence: 100%|████████████████| 8/8 [00:02<00:00, 3.72it/s]
Sampling turbulence: 100%|████████████████| 8/8 [00:06<00:00, 1.30it/s]
Computing atmospheric emission: 100%|████████████████| 1/1 [00:01<00:00, 1.34s/it, band=f150]
Sampling map: 100%|████████████████| 1/1 [00:17<00:00, 17.12s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|████████████████| 1/1 [00:01<00:00, 1.28s/it, band=f150]
2026-01-21 18:14:03.624 INFO: Simulated observation 1 of 1 in 80.65 s
[6]:
from maria.mappers import MaximumLikelihoodMapper
ml_mapper = MaximumLikelihoodMapper(tods=tods,
width=0.8 * input_map.width.deg,
height=0.8 * input_map.height.deg,
resolution=10 * input_map.resolution.deg,
units="mK_RJ")
print(f"{ml_mapper.loss() = }")
2026-01-21 18:14:14.355 INFO: Inferring center {'ra': '17ʰ19ᵐ59.97ˢ', 'dec': '-9°59’57.75”'} for mapper.
2026-01-21 18:14:14.357 INFO: Inferring mapper stokes parameters 'I' for mapper.
Preprocessing TODs: 100%|████████████████| 1/1 [00:01<00:00, 1.41s/it]
Mapping: 100%|██████████| 1/1 [00:03<00:00, 3.26s/it, tod=1/1]
Computing noise model: 100%|██████████| 1/1 [00:05<00:00, 5.37s/it, tod=1/1]
ml_mapper.loss() = tensor(0.0137, grad_fn=<AddBackward0>)
The initial map is a “guess” constructing by heavily filtering the input TODs:
[7]:
print(ml_mapper.map)
ml_mapper.map.plot(cmap="cmb")
ProjectionMap:
shape(stokes, nu, t, y, x): (1, 1, 1, 81, 81)
stokes: I
nu: [150.] GHz
t: [1.76908755e+09]
z: naive
quantity: rayleigh_jeans_temperature
units: mK_RJ
min: -1.117e-01
max: 8.590e-02
rms: 2.086e-02
center:
ra: 17ʰ19ᵐ59.97ˢ
dec: -9°59’57.75”
size(y, x): (0.791°, 0.791°)
resolution(y, x): (0.5859’, 0.5859’)
beam(maj, min, rot): [[[[0.34992376 0.34992376 0. ]]]]’
memory: 105 kB
To fit the map we run
[8]:
ml_mapper.fit(epochs=4, steps_per_epoch=32, lr=2e-1)
epoch 1/4: 100%|████████████████| 32/32 [02:14<00:00, 4.20s/it, loss=1.050e-03]
epoch 2/4: 100%|████████████████| 32/32 [02:14<00:00, 4.19s/it, loss=3.014e-04]
epoch 3/4: 100%|████████████████| 32/32 [02:14<00:00, 4.19s/it, loss=1.675e-04]
epoch 4/4: 100%|████████████████| 32/32 [02:14<00:00, 4.20s/it, loss=1.208e-04]
which gives us an improved map
[9]:
ml_mapper.map.plot(cmap="cmb")
that will improve more as it continues to fit.