04_hitmap_wcov.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hitmaps and White noise covariance\n",
    "> Load and scale hitmaps and white noise covariance matrices\n",
    "- toc: True\n",
    "\n",
    "* Basically the same scaling of the noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# default_exp hitmap_wcov"
   ]
  },
  {
   "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",
    ")\n",
    "from s4_design_sim_tool.noise import get_thinfp, get_tube_years"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulations are executed at lower sampling frequency:\n",
    "\n",
    "* $20~Hz$ for SAT (production is $100~Hz$)\n",
    "* $200~Hz$ for LAT (production is $400~Hz$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "\n",
    "simulations_sampling_frequency_scaling = dict(SAT=5, LAT=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exports\n",
    "\n",
    "hitmaps_naming = {\n",
    "    \"SAT\": \"telescope_all_time_all_hits.fits.gz\",\n",
    "    \"LAT\": \"filtered_telescope_all_time_all_hmap.fits\",\n",
    "}\n",
    "wcov_naming = {\n",
    "    \"SAT\": \"telescope_all_time_all_wcov.fits.gz\",\n",
    "    \"LAT\": \"filtered_telescope_all_time_all_wcov.fits\",\n",
    "}\n",
    "\n",
    "\n",
    "def load_hitmap_wcov(config, site, channel, realization=0, raw_hitmap=False):\n",
    "    \"\"\"Load hitmaps and white noise covariance matrices for a channel\n",
    "\n",
    "    This loads the simulated hitmaps and white noise covariance matrices\n",
    "    and scales them properly to the experiment configuration and duration\n",
    "    as defined in the input config file.\n",
    "    Hitmaps assumes a sampling frequency of 100 Hz for SAT and 400 Hz for\n",
    "    LAT.\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",
    "    hitmap : numpy array\n",
    "        Hitmap for all channels all tubes combined\n",
    "    wcov : numpy array\n",
    "        White noise covariance matrix, rows are:\n",
    "        \"II\", \"IQ\", \"IU\", \"QQ\", \"QU\", \"UU\", units are K^2\n",
    "    \"\"\"\n",
    "\n",
    "    # it is the same scaling for hitmap and white noise covariance matrix,\n",
    "    # which is the same as noise except squared\n",
    "\n",
    "    telescope = get_telescope(channel)\n",
    "    tube = channel[:-1]\n",
    "\n",
    "    channel_noP = channel.replace(\"P\", \"\")\n",
    "\n",
    "    map_filename = (\n",
    "        Path(f\"{realization:08d}\")\n",
    "        / f\"{site.lower()}_noise_{telescope}_{channel_noP}_{hitmaps_naming[telescope]}\"\n",
    "    )\n",
    "    wcov_filename = (\n",
    "        Path(f\"{realization:08d}\")\n",
    "        / f\"{site.lower()}_noise_{telescope}_{channel_noP}_{wcov_naming[telescope]}\"\n",
    "    )\n",
    "    log.info(f\"Base folder: {base_folder}\")\n",
    "    log.info(f\"Reading {map_filename}\")\n",
    "    hitmap = hp.read_map(Path(base_folder) / map_filename, 0, dtype=None, verbose=False)\n",
    "\n",
    "    hitmap[hitmap == hp.UNSEEN] = 0\n",
    "    if raw_hitmap:\n",
    "        return hitmap\n",
    "    wcov = hp.read_map(\n",
    "        Path(base_folder) / wcov_filename, (0, 1, 2, 3, 4, 5), dtype=None, verbose=False\n",
    "    )\n",
    "    scaling = (\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",
    "    ) / (10 * simulations_observing_efficiency[site.lower()][channel])\n",
    "    # focal plane thinning factor of TOD simulations\n",
    "    scaling *= get_thinfp(channel)\n",
    "    scaling *= get_tube_years(config, site, channel) / simulated_tubes[tube]\n",
    "    hitmap = np.round(hitmap * scaling).astype(np.int64)\n",
    "    hitmap *= simulations_sampling_frequency_scaling[telescope]\n",
    "    wcov /= scaling\n",
    "    wcov[:, hitmap == 0] = hp.UNSEEN\n",
    "    return hitmap, wcov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import toml\n",
    "\n",
    "config = toml.load(\"s4_design.toml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimate observation efficiency\n",
    "\n",
    "In the Winter 2021 run, there is a significant efficiency loss in the turnaround of the telescopes.\n",
    "This is already taken into account in the overall observation efficiency so we want to calibrate it out.\n",
    "Best is to use the hitmaps themselves, and scale them up from whatever number of hits to the expected\n",
    "rate for 10 days at 100% efficiency for each channel.\n",
    "\n",
    "Still we want to make sure they are within 10% of the scheduling efficiencies measured by TOAST:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scheduling_efficiency = dict(\n",
    "    pole=dict(SAT=0.642, LAT=0.754),\n",
    "    chile=dict(SAT=0.988, LAT=1.042),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s4 = read_instrument_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tot_seconds = 10 * 24 * 60 * 60\n",
    "sampling_frequency = dict(SAT=20, LAT=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observing_efficiency = dict(Pole={}, Chile={})\n",
    "for row in s4:\n",
    "    channel = row[\"band\"]\n",
    "    sites = [\"Pole\"] if row[\"telescope\"] == \"SAT\" or (\"P\" in channel) else [\"Chile\"]\n",
    "    for site in sites:\n",
    "        print(channel, site)\n",
    "        raw_hitmap = load_hitmap_wcov(config, site, channel, realization=0, raw_hitmap=True)\n",
    "        observing_efficiency[site][channel] = raw_hitmap.sum() / (tot_seconds*sampling_frequency[row[\"telescope\"]]*row[\"detectors_per_tube\"]/ get_thinfp(channel))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_tubes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scheduling_efficiency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for site, channels in observing_efficiency.items():\n",
    "    for channel in channels:\n",
    "        observing_efficiency[site][channel] /= simulated_tubes[channel[:-1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for site, channels in observing_efficiency.items():\n",
    "    for channel, efficiency in channels.items():\n",
    "        telescope = get_telescope(channel)\n",
    "        print(\n",
    "f\"{site} {channel} Efficiency from hitmap {efficiency:.2f}, scheduling efficiency {scheduling_efficiency[site.lower()][telescope]:.2f}, ratio {efficiency/scheduling_efficiency[site.lower()][telescope]:.1%}\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observing_efficiency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls /global/cscratch1/sd/keskital/s4sim/reference_tool_round_2/out/00000000/pole_noise_LAT*hmap.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verification of the input hitmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames = !ls $base_folder/00000000/*noise*wcov*"
   ]
  },
  {
   "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": [
    "input_map = hp.read_map(\n",
    "    base_folder + \"/00000000/\" + \\\n",
    "    f\"{site.lower()}_noise_{telescope}_{channel}_telescope_all_time_all_hits.fits.gz\"\n",
    ", 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fits.open(base_folder + \"/00000000/\" + \\\n",
    "    f\"{site.lower()}_noise_{telescope}_{channel}_telescope_all_time_all_wcov.fits.gz\"\n",
    ")[1].header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_EE_wcov = hp.read_map(base_folder + \"/00000000/\" + \\\n",
    "    f\"{site.lower()}_noise_{telescope}_{channel}_telescope_all_time_all_wcov.fits.gz\"\n",
    ", (3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_tubes = sum([\n",
    "    telescope[channel[:-1]] for telescope in config[\"telescopes\"][telescope].values()\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampling_rate = input_map.sum() \n",
    "sampling_rate /= 10 * 24 * 3600 * 0.5410630787037037\n",
    "sampling_rate /= 8592 # HF detectors per tube\n",
    "sampling_rate *= 8 # focal plane thinning factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampling_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check we are within 5% of 20 Hz (used in sims for SAT, LAT used 200 Hz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(sampling_rate, 20, rtol=5/100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.units as u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_variance = input_EE_wcov * u.K**2 * (input_map  / (20 * u.Hz))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_variance = input_variance[input_map > 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_variance = input_variance[input_variance > 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(input_variance.value, bins=100, log=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_NET = np.sqrt(np.median(input_variance))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_NET.to(u.uK * u.s**(.5)) / np.sqrt(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.tests.helper import assert_quantity_allclose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "channel_NET = 609.1 * u.uK * u.s**(.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The channel NET is valid for temperature, when we solve for the $2$ polarization components we get a factor of $\\sqrt{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert_quantity_allclose(input_NET/np.sqrt(2), channel_NET, rtol=5/100) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "telescope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verify the output hitmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hitmap, wcov = load_hitmap_wcov(config, site, channel, realization=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampling_rate = hitmap.sum() / (config[\"experiment\"][\"total_experiment_length_years\"] \\\n",
    "                                * config[\"experiment\"][\"observing_efficiency\"][site][telescope][channel] \\\n",
    "                                * config[\"experiment\"][\"sensitivity_factor\"][site][telescope][channel])\n",
    "sampling_rate /= 365.25 * 24 * 3600\n",
    "sampling_rate /= 8592 # HF detectors per tube\n",
    "sampling_rate /= num_tubes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampling_rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(sampling_rate, 100, rtol=5/100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare mask\n",
    "np.testing.assert_allclose(input_map == 0, hitmap == 0)"
   ]
  },
  {
   "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": [
    "log_hitmap = np.log10(hitmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_hitmap[np.isinf(log_hitmap)] = hp.UNSEEN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(log_hitmap, unit=\"log10(samples)\", min=3, max=8.5, title=\"Hitmap\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i,label in [(3, \"EE\"), (4, \"EB\"), (5, \"BB\")]:\n",
    "    hp.mollview(wcov[i], unit=\"$K^2_{CMB}$\", min=0, max=1e-11, title=f\"{label} White noise covariance\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wcov[wcov == hp.UNSEEN] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variance = wcov[3] * u.K**2 * (hitmap / (100 * u.Hz)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variance[variance != 0].min(), variance.max(), variance[variance != 0].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_mission_NET = np.sqrt(np.median(variance[hitmap != 0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_mission_NET.to(u.uK * u.s**(.5))/np.sqrt(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert_quantity_allclose(full_mission_NET/np.sqrt(2), channel_NET, rtol=5/100) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "condajupynersc",
   "language": "python",
   "name": "condajupynersc"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}