03_noise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noise\n",
    "> Load and scale 1/f and white noise simulations\n",
    "- toc: True\n",
    "\n",
    "* They include white and 1/f noise\n",
    "* Simulated for 10 days, see note about efficiency for atmosphere\n",
    "* also have a thinning factor on the number of channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# default_exp noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging as log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "\n",
    "import healpy as hp\n",
    "import numpy as np\n",
    "from pathlib import Path\n",
    "import logging as log\n",
    "\n",
    "from s4_design_sim_tool.core import (\n",
    "    get_telescope,\n",
    "    get_observing_efficiency,\n",
    "    base_folder,\n",
    "    simulations_observing_efficiency,\n",
    "    read_instrument_model,\n",
    "    mapmaking_naming,\n",
    "    simulated_tubes,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "\n",
    "thinfp_string = \"\"\"telescope\ttube\tthinfp\n",
    "LAT\tULFPL\t1\n",
    "LAT\tLFL\t4\n",
    "LAT\tLFPL\t4\n",
    "LAT\tMFL\t16\n",
    "LAT\tHFL\t16\n",
    "LAT\tMFPL\t16\n",
    "LAT\tHFPL\t16\n",
    "SAT\tLFS\t1\n",
    "SAT\tMFLS\t4\n",
    "SAT\tMFHS\t4\n",
    "SAT\tHFS\t8\n",
    "\"\"\"\n",
    "\n",
    "from astropy.io import ascii\n",
    "\n",
    "thinfp_table = ascii.read(thinfp_string)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinfp_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exports\n",
    "\n",
    "\n",
    "def get_thinfp(channel):\n",
    "    \"\"\"Get the focal plane thinning factor for noise simulations\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    channel : str\n",
    "        CMB-S4 channel tag e.g. HFL2\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    thinfp : int\n",
    "        thinning factor\n",
    "    \"\"\"\n",
    "    return (thinfp_table[thinfp_table[\"tube\"] == channel[:-1]])[\"thinfp\"][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert get_thinfp(\"ULFPL1\") == 1\n",
    "assert get_thinfp(\"HFS2\") == 8\n",
    "assert get_thinfp(\"MFLS1\") == 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exports\n",
    "\n",
    "\n",
    "def get_tube_years(config, site, channel):\n",
    "    \"\"\"Compute the number of tube/years in the CMB-S4 configuration\n",
    "\n",
    "    config_telescopes : dict\n",
    "        CMB-S4 telescopes configuration,\n",
    "        generally loaded from a TOML file\n",
    "    site : str\n",
    "        'Pole' or 'Chile', case doesn't matter\n",
    "    channel : str\n",
    "        Channel tag, e.g. 'MFHS1'\n",
    "    \"\"\"\n",
    "    tube_years = 0\n",
    "    for telescope_name, telescope_config in config[\"telescopes\"][\n",
    "        get_telescope(channel)\n",
    "    ].items():\n",
    "        if telescope_config[\"site\"].lower() == site.lower():\n",
    "            num_tubes = telescope_config.get(channel[:-1], 0)\n",
    "            tube_years += num_tubes * telescope_config.get(\n",
    "                \"years\", config[\"experiment\"][\"total_experiment_length_years\"]\n",
    "            )\n",
    "    return tube_years"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s4 = read_instrument_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tail -30 s4_design.toml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import toml\n",
    "\n",
    "config = toml.load(\"s4_design.toml\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log.basicConfig(level=log.INFO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for site in [\"Pole\", \"Chile\"]:\n",
    "    for row in s4:\n",
    "        channel = row[\"band\"]\n",
    "        tube_years = get_tube_years(config, site, channel)\n",
    "        print(site, channel, tube_years)\n",
    "        telescope = get_telescope(channel)\n",
    "        if site == \"Chile\":\n",
    "            if telescope == \"SAT\":\n",
    "                assert tube_years == 0, \"All SAT at Pole\"\n",
    "            elif channel.startswith(\"ULFL\"):\n",
    "                assert tube_years == 0, \"No ULFL in Chile\"              \n",
    "            else:\n",
    "                if channel.startswith(\"MFL\"):\n",
    "                    assert tube_years == 54*2*7, \"54+54 MFL channels\"\n",
    "        if site == \"Pole\":\n",
    "            if telescope == \"SAT\":\n",
    "                if channel.startswith(\"LFS\"):\n",
    "                    assert tube_years == 2*7, \"2 LFS Tubes\"\n",
    "                elif channel.startswith(\"MF\"):\n",
    "                    assert tube_years == 6*7, \"2 LFS Tubes\"                  \n",
    "            else:\n",
    "                if channel.startswith(\"HFPL\"):\n",
    "                    assert tube_years == 18*7, \"18 HFL tubes\"            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exports\n",
    "\n",
    "\n",
    "def load_noise(config, site, channel, realization=0):\n",
    "    \"\"\"Load noise maps for a channel\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    config : dict\n",
    "        CMB-S4 configuration,\n",
    "        generally loaded from a TOML file\n",
    "    site : str\n",
    "        'Pole' or 'Chile', case doesn't matter\n",
    "    channel : str\n",
    "        Channel tag, e.g. 'MFHS1'\n",
    "    realization : int\n",
    "        Choose one of the available 8 realizations\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    output_map : numpy array\n",
    "        Output map with all emissions combined, uses nan for missing pixels\n",
    "    \"\"\"\n",
    "    telescope = get_telescope(channel)\n",
    "    channel_noP = channel.replace(\"P\", \"\")\n",
    "\n",
    "    map_filename = (\n",
    "        Path(f\"{realization:08d}\")\n",
    "        / f\"{site.lower()}_noise_{telescope}_{channel_noP}_{mapmaking_naming[telescope]}\"\n",
    "    )\n",
    "    log.info(f\"Base folder: {base_folder}\")\n",
    "    log.info(f\"Reading {map_filename}\")\n",
    "    output_map = hp.read_map(\n",
    "        Path(base_folder) / map_filename, (0, 1, 2), dtype=None, verbose=False\n",
    "    )\n",
    "    output_map[output_map == hp.UNSEEN] = np.nan\n",
    "    output_map[output_map == 0] = np.nan\n",
    "\n",
    "    # input map is 10 days\n",
    "    output_map *= np.sqrt(\n",
    "        10\n",
    "        * simulations_observing_efficiency[site.lower()][channel]\n",
    "        / (\n",
    "            365.25\n",
    "            * get_observing_efficiency(\n",
    "                config[\"experiment\"][\"observing_efficiency\"], site, telescope, channel\n",
    "            )\n",
    "            * get_observing_efficiency(\n",
    "                config[\"experiment\"][\"sensitivity_factor\"], site, telescope, channel\n",
    "            )\n",
    "        )\n",
    "    )\n",
    "    output_map /= np.sqrt(\n",
    "        get_tube_years(config, site, channel) / simulated_tubes[channel[:-1]]\n",
    "    )\n",
    "    # focal plane thinning factor of TOD simulations\n",
    "    output_map /= np.sqrt(get_thinfp(channel))\n",
    "    return output_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Available input maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames = !ls $base_folder/00000000/*noise*bmap*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os.path\n",
    "for f in map(os.path.basename, filenames):\n",
    "    print(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "channel = \"HFS1\"\n",
    "site = \"Pole\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "telescope = get_telescope(channel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_map = hp.read_map(\n",
    "    base_folder + \"/00000000/\" + \\\n",
    "    f\"{site.lower()}_noise_{telescope}_{channel}_telescope_all_time_all_filtered.fits.gz\"\n",
    ", (0,1,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_map[input_map == 0] = hp.UNSEEN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_map = load_noise(config, site, channel, realization=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_map[np.isnan(output_map)] = hp.UNSEEN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare mask\n",
    "np.testing.assert_allclose(input_map == hp.UNSEEN, output_map == hp.UNSEEN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations_observing_efficiency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observing_efficiency = get_observing_efficiency(\n",
    "                config[\"experiment\"][\"observing_efficiency\"], site, telescope, channel\n",
    "            ) * get_observing_efficiency(\n",
    "                config[\"experiment\"][\"sensitivity_factor\"], site, telescope, channel\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observing_efficiency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(input_map[input_map != hp.UNSEEN] * \\\n",
    "    np.sqrt(10 * simulations_observing_efficiency[site.lower()][channel] / 8 / (7 * 365.25 * observing_efficiency) / 4),\n",
    "    output_map[output_map != hp.UNSEEN], rtol=1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(output_map[0], min=-1e-6, max=1e-6, unit=\"K\", title=\"Noise I\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(output_map[1], min=-1e-5, max=1e-5, unit=\"K\", title=\"Noise Q\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(output_map[2], min=-1e-5, max=1e-5, unit=\"K\", title=\"Noise U\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "condajupynersc",
   "language": "python",
   "name": "condajupynersc"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}