Merge pull request #26 from thomaswmorris/tod-in-base-sim

Inherit TODs from base sim

.github/workflows/docs.yml

@@ -6,8 +6,7 @@ on:
 
 jobs:
   build_docs:
-    # pull requests are a duplicate of a branch push if they are from within
-    # the same repo. Skip these
+    if: github.event_name != 'pull_request' || github.event.pull_request.head.repo.full_name != github.repository
     runs-on: ubuntu-latest
     strategy:
       matrix:

.github/workflows/testing.yml

@@ -1,13 +1,20 @@
-name: Tests
+name: Unit Tests
 
 on:
   push:
   pull_request:
+  schedule:
+    - cron: '00 12 * * *'  # daily at 4AM
 
 jobs:
   run_tests:
-    name: Test maria simulations
+    if: github.event_name != 'pull_request' || github.event.pull_request.head.repo.full_name != github.repository
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ["3.10"]
+      fail-fast: false
+
     defaults:
       run:
         shell: bash -l {0}

docs/source/tutorials/getting-started.ipynb

@@ -10,15 +10,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "id": "ee8a9126",
+   "execution_count": 1,
+   "id": "16efcfb6",
    "metadata": {},
    "outputs": [],
    "source": [
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
+    "from maria import Simulation\n",
     "\n",
-    "import maria "
+    "sim = Simulation(\n",
+    "    array     = 'MUSTANG-2',       # Array type\n",
+    "    pointing  = 'DAISY-2deg',      # Scanning strategy \n",
+    "    site      = 'GBT',             # Site\n",
+    "    atm_model = 'linear_angular',  # The atmospheric model, set to None if you want a noiseless observation.\n",
+    ")   \n"
    ]
   },
   {

docs/source/tutorials/ingredients.ipynb

@@ -0,0 +1,162 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "6970217f",
+   "metadata": {},
+   "source": [
+    "### Getting started"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ee8a9126",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "import maria "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6dfa512c",
+   "metadata": {},
+   "source": [
+    "In order to construct a model, we need some ingredients.\n",
+    "\n",
+    "The first ingredient is the array, which entirely defines the telescope. The Array object contains all the information about detector offsets, baselines, frequencies, and more, as well as apertures and observational constraints on the instrument.\n",
+    "\n",
+    "We can load one of `maria`'s predefined arrays ([MUSTANG-2](https://greenbankobservatory.org/science/gbt-observers/mustang-2/)) "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8e3b0f75",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m2 = maria.get_array(\"MUSTANG-2\")\n",
+    "m2.plot_dets()\n",
+    "m2.description"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ee6e0f60",
+   "metadata": {},
+   "source": [
+    "The second ingredient is the pointing, which determines where and how the telescope moves, and how fast it samples. \n",
+    "\n",
+    "We can loading a default scan with a two-degree radius from `maria`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "662ebf40",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "daisy_scan = maria.get_pointing(\"DAISY-2deg\")\n",
+    "daisy_scan.description"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "332041a8-2564-4fe7-951b-78c330b7134e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "daisy_scan.scan_radius"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "542b3f89",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure(figsize=(4,4))\n",
+    "\n",
+    "plt.plot(np.degrees(daisy_scan.ra), \n",
+    "         np.degrees(daisy_scan.dec),\n",
+    "         lw=5e-1)\n",
+    "\n",
+    "plt.xlabel(r'RA (degrees)')\n",
+    "plt.ylabel(r'Dec. (degrees)')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "352f68d6",
+   "metadata": {},
+   "source": [
+    "The last ingredient is the site, which determines the translation between local coordinates and sky coordinates. It also determines the weather parameters that go into the atmospheric simulation. \n",
+    "\n",
+    "Loading the Green Bank Telescope site from `maria`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cdb1db9d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "green_bank = maria.get_site(\"GBT\")\n",
+    "green_bank.description"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "81d5997d",
+   "metadata": {},
+   "source": [
+    "We can combine these objects by passing them to one of `maria`'s models. On construction, the simulation is initialized:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "05d4cc35",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sim = maria.Simulation(array=m2, pointing=daisy_scan, site=green_bank)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "f88e5bf5f0be7bd9faa88242bcf07e5ffa7d513506f360fd3aab8d7c0fbd2667"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

maria/atmosphere.py

@@ -60,12 +60,12 @@ def __init__(self, array, pointing, site):
     def simulate_integrated_water_vapor(self):
         raise NotImplementedError('Atmospheric simulations are not implemented in the base class!')
 
-    def run(self, units='K_RJ'):
+    def _run(self, units='K_RJ'):
 
         if units == 'K_RJ': # Kelvin Rayleigh-Jeans
 
             self.simulate_integrated_water_vapor() 
-            self.temperature = np.empty((self.array.n_dets, self.pointing.n_time))
+            self.data = np.empty((self.array.n_dets, self.pointing.n_time))
 
             for iub, uband in enumerate(self.array.ubands):
 
@@ -78,12 +78,14 @@ def run(self, units='K_RJ'):
                                                                                 self.spectrum.tcwv),
                                                                                 (self.spectrum.t_rj * passband).sum(axis=-1))
 
-                self.temperature[band_mask] = band_T_RJ_interpolator((np.degrees(self.EL[band_mask]), self.integrated_water_vapor[band_mask]))
+                self.data[band_mask] = band_T_RJ_interpolator((np.degrees(self.EL[band_mask]), self.integrated_water_vapor[band_mask]))
 
         if units == 'F_RJ': # Fahrenheit Rayleigh-Jeans 🇺🇸🇺🇸🇺🇸
 
             self.simulate_temperature(self, units='K_RJ')
-            self.temperature = 1.8 * (self.temperature - 273.15) + 32
+            self.data = 1.8 * (self.data - 273.15) + 32
+
+
 
 
 DEFAULT_LA_CONFIG = {
@@ -311,7 +313,7 @@ def __init__(self, array, pointing, site, config=DEFAULT_LA_CONFIG, verbose=Fals
                 self.C00.append(C00)
                 self.C01.append(C01)
                 self.C11.append(C11)
-                self.A.append(np.matmul(C01, utils._fast_psd_inverse(C00)))
+                self.A.append(np.matmul(C01, utils.fast_psd_inverse(C00)))
                 self.B.append(sp.linalg.lapack.dpotrf(C11 - np.matmul(self.A[-1], C01.T))[0])
 
                 prog.update(1)
@@ -667,7 +669,7 @@ def initialize(self):
         ) + alpha**2 * np.eye(self.n_cells)
         self.Cij = utils._approximate_normalized_matern(Rij, r0=self.outer_scale, nu=1 / 3, n_test_points=4096)
         self.Cjj = utils._approximate_normalized_matern(Rjj, r0=self.outer_scale, nu=1 / 3, n_test_points=4096)
-        self.A = np.matmul(self.Cij, utils._fast_psd_inverse(self.Cjj))
+        self.A = np.matmul(self.Cij, utils.fast_psd_inverse(self.Cjj))
         self.B, _ = sp.linalg.lapack.dpotrf(self.Cii - np.matmul(self.A, self.Cij.T))
 
         self.initialized = True

maria/base.py

@@ -12,9 +12,14 @@
 from os import path
 from datetime import datetime
 
-here, this_filename = os.path.split(__file__)
+
 
 from . import utils
+from .tod import TOD
+
+from astropy.io import fits
+
+here, this_filename = os.path.split(__file__)
 
 class BaseSimulation:
     """
@@ -47,7 +52,43 @@ def __init__(self,
                 out_frame="az_el",
             )
 
-
 
+    def _run(self):
+
+        raise NotImplementedError()
+
+
+    def run(self):
+
+        self._run()
+
+        tod = TOD()
+
+        tod.data = self.data # this should be set in the _run() method
+
+        tod.time = self.pointing.time
+        tod.az   = self.pointing.az
+        tod.el   = self.pointing.el
+        tod.ra   = self.pointing.ra
+        tod.dec  = self.pointing.dec
+        tod.cntr = self.pointing.scan_center
 
+        if hasattr(self, "map_sim"):
+            if self.map_sim is not None:
+                tod.unit = self.map_sim.input_map.units
+                tod.header = self.map_sim.input_map.header
+            else:
+                tod.unit = 'K'
+                tod.header = fits.header.Header()
+
+
+        tod.dets = self.array.dets
+
+        tod.meta = {'latitude': self.site.latitude,
+                    'longitude': self.site.longitude,
+                    'altitude': self.site.altitude}
+
+        return tod
+
+
 

maria/noise.py

@@ -11,7 +11,7 @@ def __init__(self, array, pointing, site, **kwargs):
 
         self.white_noise_level = kwargs.get("white_noise", 1)
 
-    def run(self):
+    def _run(self):
 
-        self.temperature = self.white_noise_level * np.random.standard_normal(size=(self.pointing.nt, self.pointing.nt))
+        self.data = self.white_noise_level * np.random.standard_normal(size=(self.pointing.nt, self.pointing.nt))