Maximum likelihood mapmaking¶
[1]:
import maria
from maria.io import fetch
input_map = maria.map.load(fetch("maps/cluster2.fits"), nu=150e9)
input_map.data *= 2e1
input_map[..., 256:-256, 256:-256].to("K_RJ").plot(cmap="cmb")
print(input_map)
2025-11-28 16:48:32.425 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster2.fits
Downloading: 100%|████████████████| 4.20M/4.20M [00:00<00:00, 57.1MB/s]
ProjectionMap:
shape(nu, y, x): (1, 1024, 1024)
stokes: naive
nu: [150.] GHz
t: naive
z: naive
quantity: spectral_flux_density_per_pixel
units: Jy/pixel
min: -2.311e-04
max: -5.838e-08
rms: 2.941e-05
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=900,
max_chunk_duration=900,
sample_rate=50,
scan_options={"radius": input_map.width.deg / 2})
plans[0].plot()
print(plans)
PlanList(1 plans, 900 s):
start_time duration target(ra,dec) center(az,el)
chunk
0 2025-11-28 16:48:33.810 +00:00 900 s (260°, -10°) (25.97°, 75.74°)
[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.2,
"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 678 11.91’ 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 678 11.91’ 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, 900 s):
│ start_time duration target(ra,dec) center(az,el)
│ chunk
│ 0 2025-11-28 16:48:33.810 +00:00 900 s (260°, -10°) (25.97°, 75.74°)
├ '2d'
└ ProjectionMap:
shape(stokes, nu, t, y, x): (1, 1, 1, 1024, 1024)
stokes: I
nu: [150.] GHz
t: [1.76434851e+09]
z: naive
quantity: spectral_flux_density_per_pixel
units: Jy/pixel
min: -2.311e-04
max: -5.838e-08
rms: 2.941e-05
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()
2025-11-28 16:48:42.347 INFO: Simulating observation 1 of 1
Constructing atmosphere: 100%|████████████████| 8/8 [00:13<00:00, 1.72s/it]
Generating turbulence: 100%|████████████████| 8/8 [00:01<00:00, 7.88it/s]
Sampling turbulence: 100%|████████████████| 8/8 [00:05<00:00, 1.35it/s]
Computing atmospheric emission: 100%|████████████████| 1/1 [00:01<00:00, 1.51s/it, band=f150]
Sampling map: 100%|████████████████| 1/1 [00:13<00:00, 13.79s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|████████████████| 1/1 [00:01<00:00, 1.30s/it, band=f150]
2025-11-28 16:49:39.837 INFO: Simulated observation 1 of 1 in 57.48 s
[6]:
from maria.mappers import MaximumLikelihoodMapper
ml_mapper = MaximumLikelihoodMapper(tods=tods,
width=0.75 * input_map.width.deg,
height=0.75 * input_map.height.deg,
resolution=10 * input_map.resolution.deg,
units="mK_RJ")
print(f"{ml_mapper.loss() = }")
2025-11-28 16:49:49.335 INFO: Inferring center {'ra': '17ʰ20ᵐ0.18ˢ', 'dec': '-10°00’10.52”'} for mapper.
2025-11-28 16:49:49.337 INFO: Inferring mapper stokes parameters 'I' for mapper.
Preprocessing TODs: 100%|████████████████| 1/1 [00:01<00:00, 1.35s/it]
Mapping: 100%|██████████| 1/1 [00:02<00:00, 2.88s/it, tod=1/1]
Computing noise model: 100%|██████████| 1/1 [00:04<00:00, 4.75s/it, tod=1/1]
ml_mapper.loss() = tensor(0.0557, 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, 76, 76)
stokes: I
nu: [150.] GHz
t: [1.76434896e+09]
z: naive
quantity: rayleigh_jeans_temperature
units: mK_RJ
min: -4.102e-01
max: 1.357e-01
rms: 6.770e-02
center:
ra: 17ʰ20ᵐ0.18ˢ
dec: -10°00’10.52”
size(y, x): (0.7422°, 0.7422°)
resolution(y, x): (0.5859’, 0.5859’)
beam(maj, min, rot): [[[[0.34992376 0.34992376 0. ]]]]’
memory: 92.42 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 [01:51<00:00, 3.48s/it, loss=3.680e-03]
epoch 2/4: 100%|████████████████| 32/32 [01:51<00:00, 3.48s/it, loss=1.065e-03]
epoch 3/4: 100%|████████████████| 32/32 [01:51<00:00, 3.48s/it, loss=5.771e-04]
epoch 4/4: 100%|████████████████| 32/32 [01:51<00:00, 3.47s/it, loss=3.561e-04]
which gives us an improved map
[9]:
ml_mapper.map.plot(cmap="cmb")
that will improve more as it continues to fit.