Units are wrong

maria/_version.py

@@ -13,5 +13,5 @@
 __version_tuple__: VERSION_TUPLE
 version_tuple: VERSION_TUPLE
 
-__version__ = version = "1.0.0"
-__version_tuple__ = version_tuple = (1, 0, 0)
+__version__ = version = "1.0.1.dev18"
+__version_tuple__ = version_tuple = (1, 0, 1, "dev18")

maria/functions.py

@@ -34,8 +34,8 @@ def normalized_matern(r, nu):
     return (
         2 ** (1 - nu)
         / sp.special.gamma(nu)
-        * sp.special.kv(nu, r + 1e-16)
-        * (r + 1e-16) ** nu
+        * sp.special.kv(nu, np.sqrt(2 * nu) * r + 1e-16)
+        * (np.sqrt(2 * nu) * r + 1e-16) ** nu
     )
 
 

maria/instrument/bands/configs/mustang2.yml → maria/instrument/bands/configs/m2.yml

@@ -9,7 +9,7 @@ f093:
 
 
 f093-custom:
-  passband: data/passbands/mustang2/f093.csv
+  passband: data/passbands/m2/f093.csv
   time_constant: 0
   NEP: 3.e-5 # pW / sqrt(s)
   knee: 1.0 # Hz

maria/instrument/bands/configs/toltec.yml

@@ -0,0 +1,26 @@
+f150:
+  center: 150
+  width: 40
+  shape: top_hat
+  time_constant: 0
+  NEP: 3.e-5 # pW / sqrt(s)
+  knee: 1.0 # Hz
+  efficiency: 0.5
+
+f220:
+  center: 220
+  width: 40
+  shape: top_hat
+  time_constant: 0
+  NEP: 3.e-5 # pW / sqrt(s)
+  knee: 1.0 # Hz
+  efficiency: 0.5
+
+f270:
+  center: 270
+  width: 50
+  shape: gaussian
+  time_constant: 0
+  NEP: 3.e-5 # pW / sqrt(s)
+  knee: 1.0 # Hz
+  efficiency: 0.5

maria/instrument/configs.yml

@@ -81,8 +81,35 @@ MUSTANG-2:
     n: 217
     field_of_view: 0.07
     array_shape: hex
-    bands: [mustang2/f093]
+    bands: [m2/f093]
     bath_temp: 100.e-3 # in K
   field_of_view: 0.07 # in degrees
   primary_size: 100 # in meters
   documentation: https://www.atlast.uio.no/
+
+
+TolTEC:
+  aliases: ["toltec"]
+  description: TolTEC
+  subarrays:
+    array-1:
+      n: 586
+      array_shape: hex
+      bands: [toltec/f150]
+      polarization: random
+      bath_temp: 100.e-3 # in K
+    array-2:
+      n: 1266
+      array_shape: hex
+      bands: [toltec/f220]
+      polarization: random
+      bath_temp: 100.e-3 # in K
+    array-3:
+      n: 2006
+      array_shape: hex
+      bands: [toltec/f270]
+      polarization: random
+      bath_temp: 100.e-3 # in K
+  field_of_view: 0.07 # in degrees
+  primary_size: 50 # in meters
+  documentation: http://toltec.astro.umass.edu/

maria/map/__init__.py

@@ -77,7 +77,7 @@ def _sample_maps(self):
             self.data["map"][band_mask] += map_power
 
 
-def from_fits(
+def read_fits(
     filename: str,
     index: int = 0,
     **map_kwargs,

maria/map/map.py

@@ -34,7 +34,8 @@ class Map:
 
     def __init__(
         self,
-        data: float,
+        name: str = None,
+        data: float = np.zeros((1, 1)),
         weight: float = None,
         width: float = None,
         resolution: float = None,
@@ -47,8 +48,13 @@ def __init__(
         if units not in UNITS_CONFIG:
             raise ValueError(f"'units' must be one of {list(UNITS_CONFIG.keys())}.")
 
+        self.name = (
+            np.atleast_1d(name)
+            if name is not None
+            else [f"{int(nu)} GHz" for nu in np.atleast_1d(frequency)]
+        )
         self.data = data if data.ndim > 2 else data[None]
-        self.weight = np.ones(self.data.shape) if weight is None else weight
+        self.weight = weight if weight is not None else np.ones(self.data.shape)
         self.center = tuple(np.radians(center)) if degrees else center
         self.frequency = np.atleast_1d(frequency)
         self.frame = frame
@@ -68,8 +74,8 @@ def __init__(
                 f"frequency dimension of the supplied map ({self.n_f})."
             )
 
-        if self.units == "Jy/pixel":
-            self.to("K_RJ")
+        # if self.units == "Jy/pixel":
+        # self.to("K_RJ", inplace=True)
 
         self.header = ap.io.fits.header.Header()
 
@@ -126,18 +132,18 @@ def y_side(self):
         return y - y.mean()
 
     def to(self, units, inplace=False):
-        data = self.data.copy()
+        data = np.zeros(self.data.shape)
 
         for i, nu in enumerate(self.frequency):
             if units == self.units:
                 data[i] = self.data[i]
 
-            if units == "K_RJ":
-                data[i] = self.data[i] * KbrightToJyPix(
+            elif units == "K_RJ":
+                data[i] = self.data[i] / KbrightToJyPix(
                     nu * 1e9, np.degrees(self.resolution)
                 )
             elif units == "Jy/pixel":
-                data[i] = self.data[i] / KbrightToJyPix(
+                data[i] = self.data[i] * KbrightToJyPix(
                     nu * 1e9, np.degrees(self.resolution)
                 )
             else:

maria/map/mappers.py

@@ -5,7 +5,6 @@
 import scipy as sp
 from tqdm import tqdm
 
-from .. import utils
 from ..tod import TOD
 from .map import Map
 
@@ -72,16 +71,44 @@ def _run(self):
         raise ValueError("Not implemented!")
 
     def run(self):
-        self.maps = {}
+        self.map_data = {}
 
         for band in self.bands:
-            self.maps[band] = self._run(band)
+            self.map_data[band] = self._run(band)
+
+        map_data = np.zeros((len(self.map_data), self.n_y, self.n_x))
+        map_weight = np.zeros((len(self.map_data), self.n_y, self.n_x))
+        map_names = []
+        map_freqs = []
+
+        for i, (band_name, band_map_data) in enumerate(self.map_data.items()):
+            map_names.append(band_name)
+            map_freqs.append(band_map_data["nom_freq"])
+
+            if "gaussian_filter" in self.map_postprocessing.keys():
+                sigma = self.map_postprocessing["gaussian_filter"]["sigma"]
+
+                band_map_numer = sp.ndimage.gaussian_filter(
+                    band_map_data["sum"], sigma=sigma
+                )
+                band_map_denom = sp.ndimage.gaussian_filter(
+                    band_map_data["weight"], sigma=sigma
+                )
+
+                map_data[i] = band_map_numer / band_map_denom
+                map_weight[i] = band_map_denom
+
+            if "median_filter" in self.map_postprocessing.keys():
+                map_data[i] = sp.ndimage.median_filter(
+                    map_data[i], size=self.map_postprocessing["median_filter"]["size"]
+                )
 
         return Map(
-            data=np.concatenate([m.data for m in self.maps.values()], axis=0),
-            weight=np.concatenate([m.weight for m in self.maps.values()], axis=0),
+            data=map_data,
+            name=map_names,
+            weight=map_weight,
+            frequency=map_freqs,
             resolution=self.resolution,
-            frequency=[float(m.frequency) for m in self.maps.values()],
             center=self.center,
             degrees=False,
         )
@@ -122,115 +149,46 @@ def _run(self, band):
         The actual mapper for the BinMapper.
         """
 
-        self.band_data = {}
-
-        self.raw_map_sums = np.zeros((self.n_x, self.n_y))
-        self.raw_map_cnts = np.zeros((self.n_x, self.n_y))
-
-        self.map_data = np.zeros((self.n_y, self.n_x))
-        self.map_weight = np.zeros((self.n_y, self.n_x))
+        band_map_data = {
+            "sum": np.zeros((self.n_y, self.n_x)),
+            "weight": np.zeros((self.n_y, self.n_x)),
+        }
 
         tods_pbar = tqdm(
             self.tods, desc=f"Running mapper ({band})", disable=not self.verbose
-        )
+        )  # noqa
 
         for tod in tods_pbar:
-            dx, dy = tod.coords.offsets(frame=self.frame, center=self.center)
-
-            band_mask = tod.dets.band_name == band
-
-            band_center = tod.dets.loc[band_mask, "band_center"].mean()
-
-            D = tod.cal[band_mask, None] * tod.data[band_mask]
-
-            # windowing
-            W = np.ones(D.shape[0])[:, None] * sp.signal.windows.tukey(
-                D.shape[-1], alpha=0.1
-            )
-
-            WD = W * sp.signal.detrend(D, axis=-1)
-            del D
-
-            if "highpass" in self.tod_postprocessing.keys():
-                WD = utils.signal.highpass(
-                    WD,
-                    fc=self.tod_postprocessing["highpass"]["f"],
-                    fs=tod.fs,
-                    order=self.tod_postprocessing["highpass"].get("order", 1),
-                    method="bessel",
-                )
+            band_tod = tod.subset(band=band)
 
-            if "remove_modes" in self.tod_postprocessing.keys():
-                n_modes_to_remove = self.tod_postprocessing["remove_modes"]["n"]
+            W, D = band_tod.process(**self.tod_postprocessing)
+            D *= band_tod.dets.cal.values[..., None]  # convert to KRJ
 
-                U, V = utils.signal.decompose(
-                    WD, downsample_rate=np.maximum(int(tod.fs / 16), 1), mode="uv"
-                )
-                WD = U[:, n_modes_to_remove:] @ V[n_modes_to_remove:]
+            dx, dy = band_tod.coords.offsets(frame=self.frame, center=self.center)
 
-            if "despline" in self.tod_postprocessing.keys():
-                B = utils.signal.get_bspline_basis(
-                    tod.time.compute(),
-                    spacing=self.tod_postprocessing["despline"].get("knot_spacing", 10),
-                    order=self.tod_postprocessing["despline"].get("spline_order", 3),
-                )
+            nu = band_tod.dets.band_center.mean()
 
-                A = np.linalg.inv(B @ B.T) @ B @ WD.T
-                WD -= A.T @ B
+            del band_tod
 
             map_sum = sp.stats.binned_statistic_2d(
-                dx[band_mask].ravel(),
-                dy[band_mask].ravel(),
-                WD.ravel(),
+                dx.ravel(),
+                dy.ravel(),
+                D.ravel(),
                 bins=(self.x_bins, self.y_bins),
                 statistic="sum",
             )[0]
 
-            map_cnt = sp.stats.binned_statistic_2d(
-                dx[band_mask].ravel(),
-                dy[band_mask].ravel(),
+            map_weight = sp.stats.binned_statistic_2d(
+                dx.ravel(),
+                dy.ravel(),
                 W.ravel(),
                 bins=(self.x_bins, self.y_bins),
                 statistic="sum",
             )[0]
 
-            self.DATA = WD
-
-            self.raw_map_sums += map_sum
-            self.raw_map_cnts += map_cnt
-
-            # self.band_data["band_center"] = tod.dets.band_center.mean()
-            # self.band_data["band_width"] = 30
-
-            band_map_numer = self.raw_map_sums.copy()
-            band_map_denom = self.raw_map_cnts.copy()
-
-            if "gaussian_filter" in self.map_postprocessing.keys():
-                sigma = self.map_postprocessing["gaussian_filter"]["sigma"]
-
-                band_map_numer = sp.ndimage.gaussian_filter(band_map_numer, sigma=sigma)
-                band_map_denom = sp.ndimage.gaussian_filter(band_map_denom, sigma=sigma)
-
-            band_map_data = band_map_numer / band_map_denom
+            band_map_data["sum"] += map_sum
+            band_map_data["weight"] += map_weight
 
-            mask = band_map_denom > 0
+        band_map_data["nom_freq"] = nu
 
-            if "median_filter" in self.map_postprocessing.keys():
-                band_map_data = sp.ndimage.median_filter(
-                    band_map_data, size=self.map_postprocessing["median_filter"]["size"]
-                )
-
-            self.map_data = np.where(mask, band_map_numer, np.nan) / np.where(
-                mask, band_map_denom, 1
-            )
-
-            self.map_weight = band_map_denom
-
-        return Map(
-            data=self.map_data,
-            weight=self.map_weight,
-            resolution=self.resolution,
-            frequency=band_center,
-            center=self.center,
-            degrees=False,
-        )
+        return band_map_data

maria/sim/base.py

@@ -154,7 +154,9 @@ def run(self, dtype=np.float32):
             coords=self.coords,
         )
 
-        tod.cal = 1 / self.instrument.dets.dP_dTRJ  # takes the data to TRJ
+        tod.dets.loc[:, "cal"] = (
+            1 / self.instrument.dets.dP_dTRJ
+        )  # takes the data to TRJ
 
         # tod.metadata = pd.Series({"pwv": self.atmosphere.weather.pwv,
         #                         "region": self.site.region})

maria/tests/data/instruments.yaml

@@ -83,7 +83,7 @@ MUSTANG-2:
     n: 217
     field_of_view: 0.07
     array_shape: hex
-    bands: [mustang2/f093]
+    bands: [m2/f093]
     bath_temp: 100.e-3 # in K
   field_of_view: 0.07 # in degrees
   primary_size: 100 # in meters

maria/tests/test_pipeline.py

@@ -24,7 +24,7 @@ def test_map_sim():
 
     map_filename = fetch("maps/cluster.fits")
 
-    input_map = maria.map.from_fits(filename=map_filename, width=0.1, center=(150, 10))
+    input_map = maria.map.read_fits(filename=map_filename, width=0.1, center=(150, 10))
 
     input_map.data *= 1e3