reorg

docs/source/tutorials/mock-observations.ipynb

@@ -47,7 +47,7 @@
    "metadata": {},
    "source": [
     "## Set up the Simulation\n",
-    "The \"Simulation\" class offers the flexibility to incorporate various inputs tailored to atmospheric simulations. For instance, it allows the creation of an observatory like 'AtLAST,' situated at the APEX site. Similarly, if the intention is to simulate Mustang-2 observations but, for instance, under less favorable weather conditions, adjusting the 'quantiles' parameter within the Simulation keywords can achieve this.\n",
+    "The \"Simulation\" class offers the flexibility to incorporate various inputs tailored to atmospheric simulations. For instance, it allows the creation of an observatory like 'AtLAST,' situated at the APEX site. Similarly, if the intention is to simulate Mustang-2 observations but, for instance, under less favorable weather conditions, adjusting the 'weather_quantiles' parameter within the Simulation keywords can achieve this.\n",
     "\n",
     "Furthermore, the Simulation class accommodates the exploration of how different scanning strategies affect the filtering of large scales. Standard strategies, such as the back-and-forth or Daisy scan patterns, have been implemented for this purpose. This feature enables the study of how the chosen scanning approach influences the overall large-scale filtering.\n",
     "\n",
@@ -93,7 +93,7 @@
     "\n",
     "    # Additional inputs:\n",
     "    # ----------------------\n",
-    "    quantiles    = {'column_water_vapor' : 0.5},    # Weather conditions specific for that site\n",
+    "    weather_quantiles    = {'column_water_vapor' : 0.5},    # Weather conditions specific for that site\n",
     "    map_units    = 'Jy/pixel',                      # Kelvin Rayleigh Jeans (K, defeault) or Jy/pixel \n",
     "    map_inbright = -2e-9,                           # Linearly scale the map to have this peak value.\n",
     "    map_res      = 0.1 / 1000,                      # degree, overwrites header information\n",

maria/__init__.py

@@ -1,37 +1,13 @@
 # Ave, María, grátia plena, Dóminus tecum
-
 from ._version import get_versions
 __version__ = get_versions()["version"]
 
 del get_versions
 
-import os
-import numpy as np
-import scipy as sp
-
-import astropy as ap
-
-import pandas as pd
-import h5py
-import glob
-import re
-import json
-import time
-import copy
-
-import weathergen
-from tqdm import tqdm
-
-import warnings
-import healpy as hp
-
-from matplotlib import pyplot as plt
-from astropy.io import fits
+from . import site
 
 from .array import get_array, get_array_config
 from .pointing import get_pointing, get_pointing_config
 from .site import get_site, get_site_config
 
-from .sim import Simulation
-
-here, this_filename = os.path.split(__file__)
+from .sim import Simulation

maria/array.py

@@ -4,26 +4,27 @@
 
 import matplotlib.pyplot as plt
 
+import glob
 import os
 
 from . import utils
 
 from collections.abc import Mapping
 
-
 here, this_filename = os.path.split(__file__)
 
-# -- Specific packages --
-ARRAY_CONFIGS = utils.read_yaml(f'{here}/configs/arrays.yml')
-
-ARRAYS = list((ARRAY_CONFIGS.keys()))
-
+ARRAY_CONFIGS = utils.io.read_yaml(f"{here}/configs/arrays.yml")
+ARRAY_PARAMS = set()
+for key, config in ARRAY_CONFIGS.items():
+    ARRAY_PARAMS |= set(config.keys())
 
 class UnsupportedArrayError(Exception):
     def __init__(self, invalid_array):
+        array_df = pd.DataFrame(columns=["description"])
+        for key, config in ARRAY_CONFIGS.items():
+            array_df.loc[key] = config["description"]
         super().__init__(f"The array \'{invalid_array}\' is not in the database of default arrays. "
-        f"Default arrays are:\n\n{sorted(list(ARRAY_CONFIGS.keys()))}")
-
+        f"Default arrays are:\n\n{array_df.sort_index()}")
 
 def get_array_config(array_name, **kwargs):
     if not array_name in ARRAY_CONFIGS.keys():
@@ -33,35 +34,18 @@ def get_array_config(array_name, **kwargs):
         ARRAY_CONFIG[k] = v
     return ARRAY_CONFIG
 
-
 def get_array(array_name, **kwargs):
     return Array(**get_array_config(array_name, **kwargs))
 
-
 def get_array_from_fits(array_name, **kwargs):
     return Array(**get_array_config(array_name, **kwargs))
 
 
 class Array:
     def __init__(self, **kwargs):
 
-        DEFAULT_ARRAY_CONFIG = {
-        "detectors": [
-            [150e9, 10e9, 100], 
-        ],
-        "geometry": "hex", 
-        "field_of_view": 1.3, 
-        "primary_size": 50,   
-        "band_grouping": "randomized",   
-        "az_bounds": [0, 360],  
-        "el_bounds": [20, 90],
-        "max_az_vel": 3,
-        "max_el_vel": 2,
-        "max_az_acc": 1,
-        "max_el_acc": 0.25
-        }
-        for k, v in kwargs.items():
-            setattr(self, k, v)
+        for key, default_value in ARRAY_CONFIGS["default"].items():
+            setattr(self, key, kwargs.get(key, default_value))
 
         detectors = kwargs.get("detectors", "")
 
@@ -97,12 +81,6 @@ def __init__(self, **kwargs):
         else:
             raise ValueError("Supplied arg 'detectors' must be either a mapping or a dataframe!")
 
-        #self.band = np.array([f"f{int(nu/10**(3*int(np.log10(nu)/3))):03}" for nu in self.band_center])
-
-
-
-
-
         # compute detector offset
         self.hull = sp.spatial.ConvexHull(self.offset)
 

maria/atmosphere.py → maria/atmosphere/__init__.py

@@ -7,7 +7,7 @@
 import warnings
 from importlib import resources
 import time as ttime
-from . import utils
+from .. import utils
 import weathergen
 from os import path
 from datetime import datetime
@@ -17,28 +17,18 @@
 
 here, this_filename = os.path.split(__file__)
 
-from . import base
-
-class AtmosphericSpectrum:
-    def __init__(self, filepath):
-        """
-        A class to hold spectra as attributes
-        """
-        with h5py.File(filepath, "r") as f:
-            self.nu = f["nu"][:]  # frequency axis of the spectrum, in GHz
-            self.tcwv = f["tcwv"][:]  # total column water vapor, in mm
-            self.elev = f["elev"][:]  # elevation, in degrees
-            self.t_rj = f["t_rj"][:]  # Rayleigh-Jeans temperature, in Kelvin
+from .. import base
 
 
 class BaseAtmosphericSimulation(base.BaseSimulation):
     """
     The base class for modeling atmospheric fluctuations.
 
-    A model needs to have the functionality to generate spectra for any pointing data we supply it with.
+    The methods to simulate e.g. line-of-sight water and temeperature profiles should be implemented by
+    classes which inherit from this one. 
     """
-    def __init__(self, array, pointing, site):
-        super().__init__(array, pointing, site)
+    def __init__(self, array, pointing, site, **kwargs):
+        super().__init__(array, pointing, site, **kwargs)
 
         self.AZ, self.EL = utils.xy_to_lonlat(
             self.array.sky_x[:, None],
@@ -49,7 +39,6 @@ def __init__(self, array, pointing, site):
 
         utils.validate_pointing(self.AZ, self.EL)
 
-        self.spectrum = AtmosphericSpectrum(filepath=f"{here}/spectra/{self.site.region}.h5")
         self.weather = weathergen.Weather(
             region=self.site.region,
             seasonal=self.site.seasonal,
@@ -71,14 +60,14 @@ def _run(self, units='K_RJ'):
 
                 band_mask = self.array.dets.band == uband
 
-                passband  = (np.abs(self.spectrum.nu - self.array.dets.band_center[band_mask].mean()) < 0.5 * self.array.dets.band_width[band_mask].mean()).astype(float)
+                passband  = (np.abs(self.site.spectrum.nu - self.array.dets.band_center[band_mask].mean()) < 0.5 * self.array.dets.band_width[band_mask].mean()).astype(float)
                 passband /= passband.sum()
 
-                band_T_RJ_interpolator = sp.interpolate.RegularGridInterpolator((self.spectrum.elev, 
-                                                                                self.spectrum.tcwv),
-                                                                                (self.spectrum.t_rj * passband).sum(axis=-1))
+                band_T_RJ_interpolator = sp.interpolate.RegularGridInterpolator((self.site.spectrum.side_pwv, 
+                                                                                 self.site.spectrum.side_elevation),
+                                                                                (self.site.spectrum.trj * passband).sum(axis=-1))
 
-                self.data[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((self.integrated_water_vapor[band_mask], np.degrees(self.EL[band_mask])))
 
         if units == 'F_RJ': # Fahrenheit Rayleigh-Jeans 🇺🇸🇺🇸🇺🇸
 
@@ -103,7 +92,7 @@ class LinearAngularSimulation(BaseAtmosphericSimulation):
     """
 
     def __init__(self, array, pointing, site, config=DEFAULT_LA_CONFIG, verbose=False, **kwargs):
-        super().__init__(array, pointing, site)
+        super().__init__(array, pointing, site, **kwargs)
 
         self.config = config
         for key, val in config.items():
@@ -119,13 +108,13 @@ def __init__(self, array, pointing, site, config=DEFAULT_LA_CONFIG, verbose=Fals
 
         # returns a beam waists and angular waists for each frequency for each layer depth
         self.waists = self.array.get_beam_waist(
-            self.layer_depths[:, None], self.array.primary_size, self.spectrum.nu
+            self.layer_depths[:, None], self.array.primary_size, self.array.dets.band_center.values
         )
         self.angular_waists = self.waists / self.layer_depths[:, None]
 
         self.min_ang_res = self.angular_waists / self.min_beam_res
 
-        self.weather.generate(time=self.pointing.time, fixed_quantiles=self.site.quantiles)
+        self.weather.generate(time=self.pointing.time, fixed_quantiles=self.site.weather_quantiles)
 
         self.heights = self.site.altitude + self.layer_depths[:, None] * np.sin(self.pointing.el)[None, :]
 
@@ -302,11 +291,11 @@ def __init__(self, array, pointing, site, config=DEFAULT_LA_CONFIG, verbose=Fals
 
                 alpha = 1e-4
 
-                C00 = utils._approximate_normalized_matern(R00, self.outer_scale / depth, 5 / 6) + alpha * np.eye(
+                C00 = utils.approximate_matern(R00, self.outer_scale / depth, 5 / 6) + alpha * np.eye(
                     len(X0)
                 )
-                C01 = utils._approximate_normalized_matern(R01, self.outer_scale / depth, 5 / 6)
-                C11 = utils._approximate_normalized_matern(R11, self.outer_scale / depth, 5 / 6) + alpha * np.eye(
+                C01 = utils.approximate_matern(R01, self.outer_scale / depth, 5 / 6)
+                C11 = utils.approximate_matern(R11, self.outer_scale / depth, 5 / 6) + alpha * np.eye(
                     len(X1)
                 )
 
@@ -664,11 +653,11 @@ def initialize(self):
         alpha = 1e-3
 
         # this is all very computationally expensive stuff (n^3 is rough!)
-        self.Cii = utils._approximate_normalized_matern(
+        self.Cii = utils.approximate_matern(
             Rii, r0=self.outer_scale, nu=1 / 3, n_test_points=4096
         ) + 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.Cij = utils.approximate_matern(Rij, r0=self.outer_scale, nu=1 / 3, n_test_points=4096)
+        self.Cjj = utils.approximate_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.B, _ = sp.linalg.lapack.dpotrf(self.Cii - np.matmul(self.A, self.Cij.T))