fix dtype conversion bugs

.github/workflows/testing.yml

@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.9", "3.11"]
+        python-version: ["3.9", "3.12"]
       fail-fast: false
 
     defaults:

README.rst

@@ -1,50 +1,58 @@
 maria
 =====
 
+.. image:: https://github.com/thomaswmorris/maria/actions/workflows/testing.yml/badge.svg
+   :target: https://github.com/thomaswmorris/maria/actions/workflows/testing.yml
+
+.. image:: https://img.shields.io/pypi/v/maria.svg
+   :target: https://pypi.python.org/pypi/maria
+
 .. image:: ./docs/source/_static/cloud.gif
    :width: 256px
    :alt: StreamPlayer
 
-`Oh, maria blows the stars around / and sends the clouds
-a-flyin’ <https://youtu.be/qKxgfnoz2pk>`_
-
-``maria`` is a python-based package that simulates turbulent atmospheric
-emission using a auto-regressive Gaussian process framework, for
-applications in observational astronomy. Tutorials for installation and
-usage can be found in the `documentation <https://www.thomaswmorris.com/maria>`_.
-
-Background
-----------
-
-Atmospheric modeling is an important step in both experiment design and
-subsequent data analysis for ground-based cosmological telescopes
-observing the cosmic microwave background (CMB). The next generation of
-ground-based CMB experiments will be marked by a huge increase in data
-acquisition: telescopes like `AtLAST <https://www.atlast.uio.no>`_ and
-`CMB-S4 <https://cmb-s4.org>`_ will consist of hundreds of thousands of
-superconducting polarization-sensitive bolometers sampling the sky. This
-necessitates new methods of efficiently modeling and simulating
-atmospheric emission at small angular resolutions, with algorithms than
-can keep up with the high throughput of modern telescopes.
-
-maria simulates layers of turbulent atmospheric emission according to a
-statistical model derived from observations of the atmosphere in the
-Atacama Desert, from the `Atacama Cosmology Telescope
-(ACT) <https://lambda.gsfc.nasa.gov/product/act/>`_ and the `Atacama
-B-Mode Search (ABS) <https://lambda.gsfc.nasa.gov/product/abs/>`_. It
-uses a sparse-precision auto-regressive Gaussian process algorithm that
-allows for both fast simulation of high-resolution atmosphere, as well
-as the ability to simulate arbitrarily long periods of atmospheric
-evolution.
-
-Methodology
------------
-
-``maria`` auto-regressively simulates an multi-layeed two-dimensional
-“integrated” atmospheric model that is much cheaper to compute than a
-three-dimensional model, which can effectively describe time-evolving
-atmospheric emission. maria can thus effectively simulate correlated
-atmospheric emission for in excess of 100,000 detectors observing the
-sky concurrently, at resolutions as fine as one arcminute. The
-atmospheric model used is detailed
-`here <https://arxiv.org/abs/2111.01319>`_.
+*maria blows the stars around / and sends the clouds a-flyin’*
+
+.. raw:: html
+    <audio controls="controls">
+      <source src="./_static/they_call_the_wind_maria.mp4" type="audio/wav">
+      Your browser does not support the <code>audio</code> element.
+    </audio>
+
+``maria`` is a complete simulator of ground-based millimeter- and submillimeter-wave telescopes. Tutorials for installation and usage can be found in the `documentation <https://www.thomaswmorris.com/maria>`_.
+
+.. Background
+.. ----------
+
+.. Atmospheric modeling is an important step in both experiment design and
+.. subsequent data analysis for ground-based cosmological telescopes
+.. observing the cosmic microwave background (CMB). The next generation of
+.. ground-based CMB experiments will be marked by a huge increase in data
+.. acquisition: telescopes like `AtLAST <https://www.atlast.uio.no>`_ and
+.. `CMB-S4 <https://cmb-s4.org>`_ will consist of hundreds of thousands of
+.. superconducting polarization-sensitive bolometers sampling the sky. This
+.. necessitates new methods of efficiently modeling and simulating
+.. atmospheric emission at small angular resolutions, with algorithms than
+.. can keep up with the high throughput of modern telescopes.
+
+.. maria simulates layers of turbulent atmospheric emission according to a
+.. statistical model derived from observations of the atmosphere in the
+.. Atacama Desert, from the `Atacama Cosmology Telescope
+.. (ACT) <https://lambda.gsfc.nasa.gov/product/act/>`_ and the `Atacama
+.. B-Mode Search (ABS) <https://lambda.gsfc.nasa.gov/product/abs/>`_. It
+.. uses a sparse-precision auto-regressive Gaussian process algorithm that
+.. allows for both fast simulation of high-resolution atmosphere, as well
+.. as the ability to simulate arbitrarily long periods of atmospheric
+.. evolution.
+
+.. Methodology
+.. -----------
+
+.. ``maria`` auto-regressively simulates an multi-layeed two-dimensional
+.. “integrated” atmospheric model that is much cheaper to compute than a
+.. three-dimensional model, which can effectively describe time-evolving
+.. atmospheric emission. maria can thus effectively simulate correlated
+.. atmospheric emission for in excess of 100,000 detectors observing the
+.. sky concurrently, at resolutions as fine as one arcminute. The
+.. atmospheric model used is detailed
+.. `here <https://arxiv.org/abs/2111.01319>`_.

docs/source/conf.py

@@ -47,4 +47,4 @@
 autoapi_type = "python"
 autoapi_dirs = ["./../../maria"]
 
-nbsphinx_execute = "never"
+nbsphinx_execute = "always"

docs/source/index.rst

@@ -9,7 +9,7 @@ maria
 .. `Oh, maria blows the stars around / and sends the clouds
 .. a-flyin’ <https://youtu.be/qKxgfnoz2pk>`_
 
-*Oh, maria blows the stars around / and sends the clouds a-flyin’*
+*maria blows the stars around / and sends the clouds a-flyin’*
 
 .. raw:: html
 

maria/base/__init__.py

@@ -55,15 +55,17 @@ class BaseSimulation:
 
     def __init__(
         self,
-        instrument: Union[Instrument, str] = "default",
-        plan: Union[Plan, str] = "stare",
-        site: Union[Site, str] = "default",
+        instrument: Union[Instrument, str],
+        plan: Union[Plan, str],
+        site: Union[Site, str],
         verbose=False,
+        dtype=np.float32,
         **kwargs,
     ):
         if hasattr(self, "boresight"):
             return
 
+        self.dtype = dtype
         self.verbose = verbose
 
         parsed_sim_kwargs = parse_sim_kwargs(kwargs, master_params)
@@ -131,14 +133,22 @@ def __init__(
     def _run(self):
         raise NotImplementedError()
 
-    def run(self, dtype=np.float32):
+    def run(self):
         self.data = {}
 
         # Simulate all the junk
         self._run()
 
+        tod_data = {}
+
+        for k, data in self.data.items():
+            # scaling floats doesn't make them more accurate, unless they're huge or tiny
+            offset = data.mean(axis=-1)[..., None]
+            # scale = data.std(axis=-1)[..., None]
+            tod_data[k] = {"data": (data - offset).astype(self.dtype), "offset": offset}
+
         tod = TOD(
-            data={k: v.astype(dtype) for k, v in self.data.items()},
+            data=tod_data,
             dets=self.instrument.dets,
             coords=self.coords,
             units="pW",

maria/coords/coordinates.py

@@ -81,7 +81,7 @@ def __init__(
         earth_location: EarthLocation = DEFAULT_EARTH_LOCATION,
         frame: str = "ra_dec",
         distributed: bool = True,
-        dtype=np.float32,
+        dtype=np.float64,
     ):
         ref_time = ttime.monotonic()
 

maria/plan/__init__.py

@@ -127,8 +127,8 @@ def __post_init__(self):
             x_scan_offsets_radians, y_scan_offsets_radians
         ].T
 
-        # add jitter, should be well sub-arcsecond
-        self.scan_offsets_radians += 1e-7 * np.random.standard_normal(
+        # add 0.1 arcseconds of jitter
+        self.scan_offsets_radians += np.radians(0.1 / 3600) * np.random.standard_normal(
             size=self.scan_offsets_radians.shape
         )
 

maria/plotting/__init__.py

@@ -0,0 +1,112 @@
+from itertools import cycle
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy as sp
+from matplotlib.patches import Patch
+
+from ..utils.signal import detrend as detrend_signal
+
+
+def plot_tod(
+    tod,
+    detrend: bool = True,
+    n_freq_bins: int = 256,
+    max_dets: int = 100,
+    lw: float = 1e0,
+    fontsize: float = 8,
+):
+    fig, axes = plt.subplots(
+        ncols=2, nrows=len(tod.fields), sharex="col", figsize=(8, 4), dpi=160
+    )
+    gs = axes[0, 0].get_gridspec()
+
+    plt.subplots_adjust(top=0.99, bottom=0.01, hspace=0, wspace=0.02)
+
+    for ax in axes[:, -1]:
+        ax.remove()
+    ps_ax = fig.add_subplot(gs[:, -1])
+
+    handles = []
+    colors = cycle(plt.rcParams["axes.prop_cycle"].by_key()["color"])
+
+    for i, field in enumerate(tod.fields):
+        data = tod.get_field(field)
+
+        tod_ax = axes[i, 0]
+
+        for band_name in np.unique(tod.dets.band_name):
+            color = next(colors)
+
+            band_mask = tod.dets.band_name == band_name
+            signal = data[band_mask]
+
+            if band_mask.sum() > max_dets:
+                signal = signal[
+                    np.random.choice(np.arange(len(signal)), max_dets, replace=False)
+                ]
+
+            detrended_signal = detrend_signal(signal, order=1)
+            if detrend:
+                signal = detrended_signal
+
+            f, ps = sp.signal.periodogram(detrended_signal, fs=tod.fs, window="tukey")
+
+            f_bins = np.geomspace(f[1], f[-1], n_freq_bins)
+            f_mids = np.sqrt(f_bins[1:] * f_bins[:-1])
+
+            binned_ps = sp.stats.binned_statistic(
+                f, ps.mean(axis=0), bins=f_bins, statistic="mean"
+            )[0]
+
+            use = binned_ps > 0
+
+            ps_ax.plot(
+                f_mids[use],
+                binned_ps[use],
+                lw=lw,
+                color=color,
+            )
+            tod_ax.plot(
+                tod.time,
+                signal[0],
+                lw=lw,
+                color=color,
+            )
+
+            handles.append(Patch(color=color, label=f"{field} ({band_name})"))
+
+        tod_ax.set_xlim(tod.time.min(), tod.time.max())
+
+        # if tod.units == "K_RJ":
+        #     ylabel = f"{field} [K]"
+        # elif tod.units == "K_CMB":
+        #     ylabel = rf"{field} [K]"
+        # elif tod.units == "pW":
+        #     ylabel = f"{field} [pW]"
+
+        # tod_ax.set_ylabel(tod.units)
+
+        label = tod_ax.text(
+            0.01,
+            0.99,
+            f"{field} [{tod.units}]",
+            fontsize=fontsize,
+            ha="left",
+            va="top",
+            transform=tod_ax.transAxes,
+        )
+        label.set_bbox(dict(facecolor="white", alpha=0.8))
+
+        # if i + 1 < n_fields:
+        # tod_ax.set_xticklabels([])
+
+    tod_ax.set_xlabel("Time [s]", fontsize=fontsize)
+
+    ps_ax.yaxis.tick_right()
+    ps_ax.yaxis.set_label_position("right")
+    ps_ax.set_xlabel("Frequency [Hz]", fontsize=fontsize)
+    ps_ax.set_ylabel(f"[{tod.units}$^2$/Hz]", fontsize=fontsize)
+    ps_ax.legend(handles=handles, loc="upper right", fontsize=fontsize)
+    ps_ax.loglog()
+    ps_ax.set_xlim(f_mids.min(), f_mids.max())

maria/sim/atmosphere.py

@@ -67,7 +67,7 @@ def _compute_atmospheric_emission(self):
                 bounds_error=False,
                 fill_value="extrapolate",
             )(self.coords.time)
-        ).astype(np.float32)
+        )
 
     def _compute_atmospheric_transmission(self):
         self.atmosphere.transmission = da.zeros_like(self.atmosphere.zenith_scaled_pwv)
@@ -129,7 +129,7 @@ def _compute_atmospheric_transmission(self):
                 bounds_error=False,
                 fill_value="extrapolate",
             )(self.coords.time)
-        ).astype(np.float32)
+        )
 
         # if units == "F_RJ":  # Fahrenheit Rayleigh-Jeans 🇺🇸
         #     self._simulate_atmospheric_emission(self, units="K_RJ")
@@ -244,7 +244,7 @@ def _classic_simulate_atmospheric_emission(self, units="K_RJ"):
         if units == "K_RJ":  # Kelvin Rayleigh-Jeans
             self._simulate_atmospheric_fluctuations()
             self.data["atmosphere"] = np.empty(
-                (self.instrument.n_dets, self.plan.n_time), dtype=np.float32
+                (self.instrument.n_dets, self.plan.n_time)
             )
 
             bands = (
@@ -291,7 +291,7 @@ def _classic_simulate_atmospheric_emission(self, units="K_RJ"):
                 )
 
             self.atmospheric_transmission = np.empty(
-                (self.instrument.n_dets, self.plan.n_time), dtype=np.float32
+                (self.instrument.n_dets, self.plan.n_time)
             )
 
             # to make a new progress bar

maria/sim/noise.py

@@ -12,7 +12,7 @@ def _run(self):
 
     def _simulate_noise(self):
         self.data["noise"] = da.from_array(
-            np.zeros((self.instrument.n_dets, self.plan.n_time), dtype=np.float32)
+            np.zeros((self.instrument.n_dets, self.plan.n_time))
         )
 
         bands = tqdm(

maria/sim/simulation.py

@@ -39,7 +39,7 @@ def from_config(cls, config: dict = {}, **params):
     def __init__(
         self,
         instrument: tuple[Instrument, str] = "default",
-        plan: tuple[Plan, str] = "daisy",
+        plan: tuple[Plan, str] = "one_minute_zenith_stare",
         site: tuple[Site, str] = "hoagie_haven",
         atmosphere: tuple[Atmosphere, str] = None,
         cmb: tuple[CMB, str] = None,
@@ -129,24 +129,21 @@ def __init__(
                 raise ValueError(f"Invalid value for cmb: '{cmb}'.")
 
     def _run(self):
-        # number of bands are lost here
-        if self.noise:
-            self._simulate_noise()
-
         if hasattr(self, "atmosphere"):
-            if self.atmosphere.model == "2d-classic":
-                self._classic_simulate_atmospheric_emission()
-            else:
-                self._simulate_atmosphere()
-                self._compute_atmospheric_emission()
-                self._compute_atmospheric_transmission()
+            self._simulate_atmosphere()
+            self._compute_atmospheric_emission()
+            self._compute_atmospheric_transmission()
 
         if hasattr(self, "cmb"):
             self._simulate_cmb_emission()
 
         if hasattr(self, "map"):
             self._sample_maps()
 
+        # number of bands are lost here
+        if self.noise:
+            self._simulate_noise()
+
         # scale the noise so that there is
         if hasattr(self, "atmospheric_transmission"):
             for k in ["cmb", "map"]:

maria/tests/test_simulation.py

@@ -32,7 +32,7 @@ def test_complete_sim(instrument, site, plan_config):
     tod = sim.run()
 
     for field in ["atmosphere", "cmb"]:
-        if np.isnan(tod.data[field]).any():
+        if np.isnan(tod.get_field(field)).any():
             raise ValueError(f"There are NaNs in the '{field}' field.")
 
     tod = tod.to("K_RJ")