{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Galactic models comparison"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from healpy.newvisufunc import projview\n",
"from healpy.rotator import Rotator\n",
"import healpy as hp\n",
"import traceback\n",
"\n",
"import pandas as pd\n",
"import numpy as np\n",
"from scipy.integrate import quad\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.ticker import MultipleLocator\n",
"import itertools\n",
"\n",
"import plotly.graph_objects as go\n",
"import plotly.express as px\n",
"\n",
"import os\n",
"import sys\n",
"from radiocalibrationtoolkit import *"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# This ensures Plotly output works in multiple places:\n",
"# plotly_mimetype: VS Code notebook UI\n",
"# notebook: \"Jupyter: Export to HTML\" command in VS Code\n",
"# See https://plotly.com/python/renderers/#multiple-renderers\n",
"import plotly.io as pio\n",
"pio.renderers.default = \"plotly_mimetype+notebook\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# some global plot settings\n",
"plt.rcParams[\"axes.titleweight\"] = \"bold\"\n",
"plt.rcParams[\"axes.labelweight\"] = \"bold\"\n",
"plt.rcParams[\"font.weight\"] = \"bold\"\n",
"plt.rcParams['font.size'] = 16\n",
"plt.rcParams['legend.fontsize']= 14\n",
"\n",
"plt.rcParams['xtick.major.width']= 2\n",
"plt.rcParams['ytick.major.width']= 2\n",
"\n",
"plt.rcParams['xtick.major.size']= 5\n",
"plt.rcParams['ytick.major.size']= 5\n",
"\n",
"plt.rcParams['xtick.labelsize']= 14\n",
"plt.rcParams['ytick.labelsize']= 14"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Prepare objects\n",
"lfmap = LFmap()\n",
"lfss = LowFrequencySkyModel(freq_unit=\"MHz\")\n",
"gsm2008 = GlobalSkyModel(freq_unit=\"MHz\")\n",
"gsm2016 = GlobalSkyModel2016(freq_unit=\"MHz\")\n",
"haslam = HaslamSkyModel(freq_unit=\"MHz\", spectral_index=-2.53)\n",
"ssm = SSM()\n",
"gmoss = GMOSS()\n",
"ulsa_fdi = ULSA(index_type='freq_dependent_index')\n",
"ulsa_ci = ULSA(index_type='constant_index')\n",
"ulsa_dpi = ULSA(index_type='direction_dependent_index')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"frequency_MHz = 45\n",
"\n",
"lfmap_map = lfmap.generate(frequency_MHz)\n",
"lfss_map = lfss.generate(frequency_MHz)\n",
"gsm2008_map = gsm2008.generate(frequency_MHz)\n",
"gsm2016_map = gsm2016.generate(frequency_MHz)\n",
"haslam_map = haslam.generate(frequency_MHz)\n",
"ssm_map = ssm.generate(frequency_MHz)\n",
"gmoss_map = gmoss.generate(frequency_MHz)\n",
"ulsa_fdi_map = ulsa_fdi.generate(frequency_MHz)\n",
"ulsa_ci_map = ulsa_ci.generate(frequency_MHz)\n",
"ulsa_dpi_map = ulsa_dpi.generate(frequency_MHz)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Check size of the arrays and NSIDE number and angular resolution\n",
"\n",
"map_list = [lfmap_map, lfss_map, gsm2008_map, gsm2016_map, haslam_map, ssm_map, gmoss_map, ulsa_fdi_map, ulsa_ci_map, ulsa_dpi_map]\n",
"[print_map_properties(m) for m in map_list]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# convert to same NSIDE, NSIDE is always some power of two\n",
"new_nside = 64\n",
"lfmap_map_N = hp.ma(hp.pixelfunc.ud_grade(lfmap_map, new_nside))\n",
"lfss_map_N = hp.ma(hp.pixelfunc.ud_grade(lfss_map, new_nside))\n",
"gsm2008_map_N = hp.ma(hp.pixelfunc.ud_grade(gsm2008_map, new_nside))\n",
"gsm2016_map_N = hp.ma(hp.pixelfunc.ud_grade(gsm2016_map, new_nside))\n",
"ssm_map_N = hp.ma(hp.pixelfunc.ud_grade(ssm_map, new_nside))\n",
"haslam_map_N = hp.ma(hp.pixelfunc.ud_grade(haslam_map, new_nside))\n",
"gmoss_map_N = hp.ma(hp.pixelfunc.ud_grade(gmoss_map, new_nside))\n",
"ulsa_fdi_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_fdi_map, new_nside))\n",
"ulsa_ci_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_ci_map, new_nside))\n",
"ulsa_dpi_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_dpi_map, new_nside))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"map_dict = {\n",
" \"LFSS\": lfss_map_N,\n",
" \"GSM08\": gsm2008_map_N,\n",
" \"GSM16\": gsm2016_map_N,\n",
" \"Haslam\": haslam_map_N,\n",
" \"LFmap\": lfmap_map_N,\n",
" \"SSM\": ssm_map_N,\n",
" \"GMOSS\": gmoss_map_N,\n",
" \"ULSA\": ulsa_fdi_map_N,\n",
" # \"ULSA2\": ulsa_ci_map_N,\n",
" # \"ULSA3\": ulsa_dpi_map_N,\n",
"}\n",
"\n",
"ulsa_map_dict = {\n",
" \"ULSA FDI\": ulsa_fdi_map_N,\n",
" \"ULSA CI\": ulsa_ci_map_N,\n",
" \"ULSA DPI\": ulsa_dpi_map_N,\n",
"}\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Print the maps in the dictionary\n",
"[(key, map_dict[key]) for key in map_dict]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Galactic coordinates"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Prints all maps\n",
"cmap = \"jet\"\n",
"projection_type = \"mollweide\"\n",
"\n",
"for key in map_dict:\n",
"\n",
" projview(\n",
" map_dict[key],\n",
" norm=\"log\",\n",
" coord=[\"G\"],\n",
" graticule=True,\n",
" graticule_labels=True,\n",
" unit=\"Temperature ln[K]\",\n",
" xlabel=\"RA\",\n",
" ylabel=\"DEC\",\n",
" cb_orientation=\"vertical\",\n",
" min=3500,\n",
" max=35000,\n",
" latitude_grid_spacing=30,\n",
" projection_type=projection_type,\n",
" title=key,\n",
" xtick_label_color=\"white\",\n",
" cmap=cmap,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"compare_maps(\n",
" map_dict,\n",
" main_title=\"Frequency: {} MHz; Compared on {:.2f}{} angular resolution\".format(\n",
" frequency_MHz, np.rad2deg(hp.pixelfunc.nside2resol(new_nside)), chr(176)\n",
" ),\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"compare_maps(\n",
" ulsa_map_dict,\n",
" main_title=\"Frequency: {} MHz; Compared on {:.2f}{} angular resolution\".format(\n",
" frequency_MHz, np.rad2deg(hp.pixelfunc.nside2resol(new_nside)), chr(176)\n",
" ),\n",
" figsize=(6, 4.5)\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compare at all frequencies"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f_range = range(30, 85, 5)\n",
"for frequency_MHz in f_range:\n",
" # LFmap is by default in Celestial coordinates, to make things simpler, we rotate it to Galactic coordinates at the creation stage.\n",
" lfmap_map = hp.rotator.Rotator.rotate_map_pixel(\n",
" Rotator(coord=[\"C\", \"G\"]), lfmap.generate(frequency_MHz)\n",
" )\n",
" lfmap_map = lfmap.generate(frequency_MHz)\n",
" lfss_map = lfss.generate(frequency_MHz)\n",
" gsm2008_map = gsm2008.generate(frequency_MHz)\n",
" gsm2016_map = gsm2016.generate(frequency_MHz)\n",
" haslam_map = haslam.generate(frequency_MHz)\n",
" ssm_map = ssm.generate(frequency_MHz)\n",
" gmoss_map = gmoss.generate(frequency_MHz)\n",
" ulsa_fdi_map = ulsa_fdi.generate(frequency_MHz)\n",
" ulsa_ci_map = ulsa_ci.generate(frequency_MHz)\n",
" ulsa_dpi_map = ulsa_dpi.generate(frequency_MHz)\n",
"\n",
" # convert to same NSIDE\n",
" new_nside = 64\n",
" lfmap_map_N = hp.ma(hp.pixelfunc.ud_grade(lfmap_map, new_nside))\n",
" lfss_map_N = hp.ma(hp.pixelfunc.ud_grade(lfss_map, new_nside))\n",
" gsm2008_map_N = hp.ma(hp.pixelfunc.ud_grade(gsm2008_map, new_nside))\n",
" gsm2016_map_N = hp.ma(hp.pixelfunc.ud_grade(gsm2016_map, new_nside))\n",
" ssm_map_N = hp.ma(hp.pixelfunc.ud_grade(ssm_map, new_nside))\n",
" haslam_map_N = hp.ma(hp.pixelfunc.ud_grade(haslam_map, new_nside))\n",
" gmoss_map_N = hp.ma(hp.pixelfunc.ud_grade(gmoss_map, new_nside))\n",
" ulsa_fdi_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_fdi_map, new_nside))\n",
" ulsa_ci_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_ci_map, new_nside))\n",
" ulsa_dpi_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_dpi_map, new_nside))\n",
"\n",
" map_temp_dict = {\n",
" \"LFSS\": lfss_map_N,\n",
" \"GSM08\": gsm2008_map_N,\n",
" \"GSM16\": gsm2016_map_N,\n",
" \"Haslam\": haslam_map_N,\n",
" \"LFmap\": lfmap_map_N,\n",
" \"SSM\": ssm_map_N,\n",
" \"GMOSS\": gmoss_map_N,\n",
" \"ULSA\": ulsa_fdi_map_N,\n",
" # \"ULSA2\": ulsa_ci_map_N,\n",
" # \"ULSA3\": ulsa_dpi_map_N,\n",
" }\n",
"\n",
"\n",
" compare_dict = compare_maps(\n",
" map_temp_dict,\n",
" show_plot_comparison=False,\n",
" verbose=False\n",
" )\n",
"\n",
" if frequency_MHz == list(f_range)[0]:\n",
" compare_DF = pd.DataFrame(compare_dict).T\n",
" compare_DF.columns = [list(f_range)[0]]\n",
" else:\n",
" compare_DF[frequency_MHz] = pd.DataFrame(compare_dict).T.values"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = px.imshow(\n",
" compare_DF,\n",
" width=500,\n",
" height=600,\n",
" color_continuous_scale=\"jet\",\n",
" zmin=-0.3,\n",
" zmax=0.3,\n",
" labels={\"x\": \"frequency [MHz]\"},\n",
")\n",
"fig.update_layout(\n",
" margin=dict(l=20, r=120, t=20, b=20),\n",
" xaxis={\"nticks\": 5},\n",
" font=dict(family=\"Arial Black\", size=15, color=\"black\"),\n",
")\n",
"fig.update_yaxes(tickprefix=\"\", ticksuffix=\"\")\n",
"fig.update_xaxes(\n",
" tickprefix=\"\", ticksuffix=\"\", tick0=30, dtick=10, tickfont=dict(size=18)\n",
")\n",
"\n",
"fig.update_coloraxes(\n",
" colorbar=dict(\n",
" tickprefix=\"\",\n",
" ticksuffix=\"\",\n",
" tickfont=dict(size=18),\n",
" )\n",
")\n",
"\n",
"fig.add_annotation(\n",
" go.layout.Annotation(\n",
" text='RDIS',\n",
" x=1.4,\n",
" y=0.5,\n",
" xref='paper',\n",
" yref='paper',\n",
" showarrow=False,\n",
" font=dict(size=18),\n",
" textangle=-90,\n",
" align='center',\n",
" valign='middle',\n",
" )\n",
")\n",
"\n",
"fig.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"ax.set_title('all \\n')\n",
"\n",
"# temp = compare_DF.std(axis=0)\n",
"labels = compare_DF.columns.astype(float)\n",
"ax.errorbar(\n",
" labels,\n",
" compare_DF.mean(axis=0).values,\n",
" yerr=compare_DF.std(axis=0),\n",
" marker=\".\",\n",
" linestyle=\"\",\n",
" markersize=10,\n",
")\n",
"\n",
"ax.axes.axhline(0, ls=\"--\", c=\"r\")\n",
"ax.set_xlabel(\"frequency [MHz]\")\n",
"ax.set_ylabel(\"$_{\\mathrm{models}}$\")\n",
"ax.set_ylim(-0.15, 0.15)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.set_title('without LFSS\\n')\n",
"\n",
"ax.errorbar(\n",
" labels,\n",
" compare_DF.iloc[7:,:].mean(axis=0).values,\n",
" yerr=compare_DF.std(axis=0),\n",
" marker=\".\",\n",
" linestyle=\"\",\n",
" markersize=10,\n",
")\n",
"\n",
"ax.axes.axhline(0, ls=\"--\", c=\"r\")\n",
"ax.set_xlabel(\"frequency [MHz]\")\n",
"ax.set_ylabel(\"$_{\\mathrm{models}}$\")\n",
"ax.set_ylim(-0.15, 0.15)\n",
"\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.set_title('without LFSS and GSM08\\n')\n",
"\n",
"ax.errorbar(\n",
" labels,\n",
" compare_DF.iloc[13:,:].mean(axis=0).values,\n",
" yerr=compare_DF.std(axis=0),\n",
" marker=\".\",\n",
" linestyle=\"\",\n",
" markersize=10,\n",
")\n",
"\n",
"ax.axes.axhline(0, ls=\"--\", c=\"r\")\n",
"ax.set_xlabel(\"frequency [MHz]\")\n",
"ax.set_ylabel(\"$_{\\mathrm{models}}$\")\n",
"ax.set_ylim(-0.15, 0.15)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"compare_DF.iloc[12:,:]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"df = compare_DF.copy(deep=True)\n",
"df = compare_DF.iloc[::-1,:]\n",
"fig, ax = plt.subplots(figsize=(6, 8))\n",
"labels = df.index.values\n",
"err = df.std(axis=1)\n",
"\n",
"pos = range(labels.size)\n",
"ax.errorbar(\n",
" df.mean(axis=1), pos, xerr=err, marker=\".\", linestyle=\"\", markersize=10\n",
")\n",
"ax.set_yticks(pos, labels=labels)\n",
"ax.set_xlabel(\"$_{\\\\nu}$\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"input_data1 = compare_DF.values.flatten()\n",
"input_data2 = compare_DF.iloc[7:, :].values.flatten()\n",
"\n",
"bins = np.linspace(-0.5, 0.5, 20)\n",
"fig, ax = plt.subplots()\n",
"fig.suptitle('')\n",
"ax.hist(input_data1, bins=bins, density=True, alpha=0.5)\n",
"ax.hist(input_data2, bins=bins, density=True, alpha=0.5)\n",
"\n",
"ax.set_xlabel(\"RDIS\")\n",
"ax.set_ylabel(\"entries\")\n",
"\n",
"fig, ax = plt.subplots()\n",
"# fig.suptitle('Blue: all models, Red: without LFSS')\n",
"ax.set_xlabel(\"RDIS\")\n",
"apply_KDE(input_data1, bounds=[bins[0], bins[-1]], data_in_relative_values=True)\n",
"apply_KDE(input_data2, bounds=[bins[0], bins[-1]], data_in_relative_values=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"display(input_data1.size)\n",
"display(input_data2.size)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Local obsever"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"lst = 18\n",
"latitude = -35.206667\n",
"frequency_MHz = 45\n",
"main_title = \"Frequency: {} MHz; Compared on {:.2f}{} angular resolution; LST:{}; Latitude:{:.2f}\".format(\n",
" frequency_MHz, np.rad2deg(hp.pixelfunc.nside2resol(new_nside)), chr(176), lst, latitude\n",
")\n",
"\n",
"# Local coordinates at LST time \"LSTtime\" at latitude \"latitude\"\n",
"rotation_parameters = [(180 + (lst * 15)) % 360, -(latitude - 90)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# some plot settings\n",
"cmap = \"jet\"\n",
"projection_type = \"mollweide\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Prints all maps\n",
"for key in map_dict:\n",
"\n",
" projview(\n",
" map_dict[key],\n",
" norm=\"log\",\n",
" coord=['G','C'],\n",
" rot=rotation_parameters,\n",
" graticule=True,\n",
" graticule_labels=True,\n",
" unit=\"Temperature ln[K]\",\n",
" xlabel=\"azimuth\",\n",
" ylabel=\"altitude\",\n",
" cb_orientation=\"vertical\",\n",
" min=3500,\n",
" max=35000,\n",
" latitude_grid_spacing=30,\n",
" projection_type=projection_type,\n",
" title=key,\n",
" xtick_label_color=\"white\",\n",
" cmap=cmap,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create mask, at local coordinates we see only a hemisphere\n",
"mask = create_local_mask(new_nside, rotation_parameters)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Test\n",
"m = map_dict['LFmap'].copy()\n",
"m.mask = mask\n",
"\n",
"projview(\n",
" m,\n",
" norm=\"log\",\n",
" coord=[\"G\", 'C'],\n",
" rot=rotation_parameters,\n",
" graticule=True,\n",
" graticule_labels=True,\n",
" unit=\"Temperature ln[K]\",\n",
" xlabel=\"azimuth\",\n",
" ylabel=\"altitude\",\n",
" cb_orientation=\"vertical\",\n",
" min=3500,\n",
" max=35000,\n",
" latitude_grid_spacing=30,\n",
" xtick_label_color=\"white\",\n",
" title=main_title,\n",
" cmap=cmap,\n",
")\n",
"\n",
"projview(\n",
" m,\n",
" norm=\"log\",\n",
" # coord=[\"G\"],\n",
" # rot=False,\n",
" graticule=True,\n",
" graticule_labels=True,\n",
" unit=\"Temperature ln[K]\",\n",
" xlabel=\"RA\",\n",
" ylabel=\"DEC\",\n",
" cb_orientation=\"vertical\",\n",
" min=3500,\n",
" max=35000,\n",
" latitude_grid_spacing=30,\n",
" xtick_label_color=\"white\",\n",
" title=main_title,\n",
" cmap=cmap,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### In local coordinates"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"compare_maps(\n",
" map_dict, rotation_parameters=rotation_parameters, coord=['G','C'], main_title=main_title, mask=mask\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### In galactic coordinates (masked)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"compare_maps(map_dict, rotation_parameters=False, main_title=main_title, mask=mask)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Complex comparison"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f_range = range(30, 85, 5)\n",
"lst_range = range(0, 24, 2)\n",
"# f_range = range(30, 45, 5)\n",
"# lst_range = range(0, 6, 2)\n",
"compare_DF = pd.DataFrame(columns=['label', 'frequency [MHz]', 'LST', 'ratio'])\n",
"\n",
"for frequency_MHz in f_range:\n",
" lfmap_map = hp.rotator.Rotator.rotate_map_pixel(\n",
" Rotator(coord=[\"C\", \"G\"]), lfmap.generate(frequency_MHz)\n",
" )\n",
" lfmap_map = lfmap.generate(frequency_MHz)\n",
" lfss_map = lfss.generate(frequency_MHz)\n",
" gsm2008_map = gsm2008.generate(frequency_MHz)\n",
" gsm2016_map = gsm2016.generate(frequency_MHz)\n",
" haslam_map = haslam.generate(frequency_MHz)\n",
" ssm_map = ssm.generate(frequency_MHz)\n",
" gmoss_map = gmoss.generate(frequency_MHz)\n",
" ulsa_fdi_map = ulsa_fdi.generate(frequency_MHz)\n",
" ulsa_ci_map = ulsa_ci.generate(frequency_MHz)\n",
" ulsa_dpi_map = ulsa_dpi.generate(frequency_MHz)\n",
"\n",
" # convert to same NSIDE\n",
" new_nside = 64\n",
" lfmap_map_N = hp.ma(hp.pixelfunc.ud_grade(lfmap_map, new_nside))\n",
" lfss_map_N = hp.ma(hp.pixelfunc.ud_grade(lfss_map, new_nside))\n",
" gsm2008_map_N = hp.ma(hp.pixelfunc.ud_grade(gsm2008_map, new_nside))\n",
" gsm2016_map_N = hp.ma(hp.pixelfunc.ud_grade(gsm2016_map, new_nside))\n",
" ssm_map_N = hp.ma(hp.pixelfunc.ud_grade(ssm_map, new_nside))\n",
" haslam_map_N = hp.ma(hp.pixelfunc.ud_grade(haslam_map, new_nside))\n",
" gmoss_map_N = hp.ma(hp.pixelfunc.ud_grade(gmoss_map, new_nside))\n",
" ulsa_fdi_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_fdi_map, new_nside))\n",
" ulsa_ci_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_ci_map, new_nside))\n",
" ulsa_dpi_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_dpi_map, new_nside))\n",
"\n",
" map_temp_dict = {\n",
" \"LFSS\": lfss_map_N,\n",
" \"GSM08\": gsm2008_map_N,\n",
" \"GSM16\": gsm2016_map_N,\n",
" \"Haslam\": haslam_map_N,\n",
" \"LFmap\": lfmap_map_N,\n",
" \"SSM\": ssm_map_N,\n",
" \"GMOSS\": gmoss_map_N,\n",
" \"ULSA\": ulsa_fdi_map_N,\n",
" # \"ULSA2\": ulsa_ci_map_N,\n",
" # \"ULSA3\": ulsa_dpi_map_N,\n",
" }\n",
" \n",
" for lst in lst_range:\n",
" rotation_parameters = [(180 + (lst * 15)) % 360, -(latitude - 90)]\n",
" mask = create_local_mask(new_nside, rotation_parameters)\n",
"\n",
" compare_dict = compare_maps(\n",
" map_temp_dict,\n",
" rotation_parameters=False,\n",
" mask=mask,\n",
" show_plot_comparison=False,\n",
" verbose=False\n",
" )\n",
"\n",
" for key in compare_dict:\n",
" compare_DF.loc[compare_DF.index.size] = [key, frequency_MHz, lst, compare_dict[key][0]]\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"compare_DF\n",
"df = compare_DF.sort_values(by=['label', 'frequency [MHz]', 'LST'])\n",
"models_ratio_list = list(dict.fromkeys(df.label.values))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"default_color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
"# delete red\n",
"r = default_color_cycle[3]\n",
"del default_color_cycle[3]\n",
"b='k'\n",
"# b = default_color_cycle[0]\n",
"# del default_color_cycle[0]\n",
"\n",
"fig, ax = plt.subplots(1,2, figsize=(10,4))\n",
"\n",
"for i, frequency in enumerate(f_range):\n",
" for models_ratio in models_ratio_list:\n",
" temp = df[(df['label']==models_ratio) & (df['frequency [MHz]']==frequency)]\n",
" if 'LFSS' in models_ratio:\n",
" ax[0].plot(temp.LST.values, temp.ratio.values, c=r, alpha=1)#, zorder=len(f_range)+i)\n",
" elif 'GSM08' in models_ratio:\n",
" ax[0].plot(temp.LST.values, temp.ratio.values, alpha=0.5)#, zorder=len(f_range)+i)\n",
" else:\n",
" ax[0].plot(temp.LST.values, temp.ratio.values, alpha=0.5)\n",
"\n",
"ax[0].set_xlabel('LST')\n",
"ax[0].set_ylabel('RDIS')\n",
"\n",
"for lst in lst_range:\n",
" for models_ratio in models_ratio_list:\n",
" temp = df[(df['label']==models_ratio) & (df['LST']==lst)]\n",
" if 'LFSS' in models_ratio:\n",
" ax[1].plot(temp['frequency [MHz]'].values, temp.ratio.values, c=r, alpha=1)#, zorder=len(f_range)+i)\n",
" elif 'GSM08' in models_ratio:\n",
" ax[1].plot(temp['frequency [MHz]'].values, temp.ratio.values, alpha=0.5)#, zorder=len(f_range)+i)\n",
" else:\n",
" ax[1].plot(temp['frequency [MHz]'].values, temp.ratio.values, alpha=0.5)\n",
"\n",
"ax[1].set_xlabel('frequency [MHz]')\n",
"ax[1].set_ylabel('RDIS')\n",
"fig.subplots_adjust(wspace=0.3)\n",
"\n",
"ax[0].xaxis.set_major_locator(MultipleLocator(4))\n",
"ax[1].xaxis.set_major_locator(MultipleLocator(10))\n",
"\n",
"# Reset to default color cycle\n",
"plt.rcParams['axes.prop_cycle'] = plt.rcParamsDefault['axes.prop_cycle']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"center_DF = pd.DataFrame()\n",
"center_DF.index = lst_range\n",
"left_DF = center_DF.copy(deep=True)\n",
"right_DF = center_DF.copy(deep=True)\n",
"\n",
"bins=np.linspace(-0.5, 0.5, 20)\n",
"\n",
"for frequency in f_range:\n",
" center = []\n",
" left = []\n",
" right = []\n",
" for lst in lst_range:\n",
" data = df[(df['frequency [MHz]'] == frequency) & (df['LST'] == lst)]['ratio'].values\n",
"\n",
" xax_kde = np.linspace(bins[0], bins[-1], 2000)[:, np.newaxis]\n",
" kde = KernelDensity(kernel=\"gaussian\", bandwidth=0.05).fit(data[:, np.newaxis])\n",
" log_dens = kde.score_samples(xax_kde)\n",
"\n",
" xax_ = xax_kde.flatten()\n",
"\n",
" center_i = find_index_on_CDF(log_dens, xax_kde, 0.5)\n",
" left_i = find_index_on_CDF(log_dens, xax_kde, 0.5-0.341)\n",
" right_i = find_index_on_CDF(log_dens, xax_kde, 0.5+0.341)\n",
"\n",
" center.append(xax_[center_i])\n",
" left.append(xax_[left_i+1])\n",
" right.append(xax_[right_i])\n",
" center_DF[frequency] = center\n",
" left_DF[frequency] = left\n",
" right_DF[frequency] = right\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"font = \"Arial Black\"\n",
"\n",
"scene = dict(\n",
" xaxis=dict(\n",
" title=\"frequency [MHz]\",\n",
" color=\"black\",\n",
" tickvals=list(range(len(f_range))),\n",
" ticktext=center_DF.columns.values,\n",
" ),\n",
" yaxis=dict(\n",
" title=\"LST [hour]\",\n",
" color=\"black\",\n",
" tickvals=list(range(len(lst_range))),\n",
" ticktext=center_DF.index.values.tolist(),\n",
" ),\n",
" zaxis=dict(title=\"\", range=(-0.1, 0.05), dtick=0.05),\n",
" aspectratio=dict(x=0.75, y=0.75, z=0.75), # Adjust the aspect ratio as needed\n",
" camera=dict(eye=dict(x=1.2, y=1.2, z=0.6)), # Adjust the camera position as needed\n",
")\n",
"\n",
"# layout = go.Layout(scene=scene)\n",
"\n",
"fig = go.Figure(\n",
" data=[\n",
" # go.Surface(z=right_DF),\n",
" go.Surface(z=center_DF, showscale=True, opacity=0.9)\n",
" ],\n",
" # go.Surface(z=left_DF, showscale=False, opacity=0.9)],\n",
")\n",
"\n",
"fig.update_layout(\n",
" scene=scene,\n",
" margin=dict(r=50, b=10, l=10, t=10),\n",
" height=600,\n",
" autosize=False,\n",
" font=dict(family=font, size=18, color=\"black\"),\n",
")\n",
"\n",
"fig.update_traces(\n",
" colorbar=dict(\n",
" tickprefix=\"\",\n",
" ticksuffix=\"\",\n",
" title=\"\",\n",
" titleside=\"right\",\n",
" tickfont=dict(size=18),\n",
" len=0.75, # Adjust the length of the colorbar\n",
" )\n",
")\n",
"\n",
"fig.show()\n",
"################################################################\n",
"################################################################\n",
"################################################################\n",
"\n",
"fig = go.Figure(\n",
" data=[\n",
" go.Surface(z=(center_DF - right_DF).abs(), opacity=0.5),\n",
" # go.Surface(z=center_DF, showscale=False, opacity=0.9)],\n",
" go.Surface(z=(center_DF - left_DF).abs(), showscale=False, opacity=0.5),\n",
" ],\n",
")\n",
"\n",
"scene[\"zaxis\"] = {\"title\": \" errors\", \"dtick\": 0.05}\n",
"fig.update_layout(\n",
" scene=scene,\n",
" margin=dict(r=20, b=10, l=10, t=10),\n",
" height=600,\n",
" autosize=False,\n",
" font=dict(family=font, size=18, color=\"black\"),\n",
")\n",
"\n",
"fig.update_traces(\n",
" colorbar=dict(\n",
" tickprefix=\"\",\n",
" ticksuffix=\"\",\n",
" title=\" errors\",\n",
" titleside=\"right\",\n",
" tickfont=dict(size=18),\n",
" len=0.75, # Adjust the length of the colorbar\n",
" )\n",
")\n",
"\n",
"fig.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"input_data = df.ratio.values\n",
"# input_data = df[(df['frequency [MHz]'] >45) & (df['frequency [MHz]'] <999) & ((df['LST'] >=0) | (df['LST'] <=777))].ratio.values\n",
"input_data = df[(df['frequency [MHz]'] >0) & (df['frequency [MHz]'] <99)].ratio.values\n",
"\n",
"########################################\n",
"fig, ax = plt.subplots()\n",
"\n",
"bins=np.linspace(-0.5, 0.5, 20)\n",
"hist_data = ax.hist(input_data, bins=bins)\n",
"\n",
"xax_kde = np.linspace(bins[0], bins[-1], 2000)[:, np.newaxis]\n",
"kde = KernelDensity(kernel=\"gaussian\", bandwidth=0.05).fit(input_data[:, np.newaxis])\n",
"log_dens = kde.score_samples(xax_kde)\n",
"\n",
"ax.set_xlabel('RDIS')\n",
"ax.set_ylabel('entries')\n",
"########################################\n",
"\n",
"fig, ax = plt.subplots()\n",
"\n",
"xax_ = xax_kde.flatten()\n",
"ax.plot(xax_, np.exp(log_dens))\n",
"\n",
"center_i = find_index_on_CDF(log_dens, xax_kde, 0.5)\n",
"left_i = find_index_on_CDF(log_dens, xax_kde, 0.5-0.341)\n",
"right_i = find_index_on_CDF(log_dens, xax_kde, 0.5+0.341)\n",
"\n",
"print('Test1: p0.5={}'.format(np.sum(np.exp(log_dens[:center_i+1])*np.diff(xax_)[0]) - 0.5))\n",
"print('Test2: p0.5={}'.format(np.sum(np.exp(log_dens[center_i+1:])*np.diff(xax_)[0]) - 0.5))\n",
"print('Test3: p0.5+/-0.341={}'.format(np.sum(np.exp(log_dens[left_i+1:right_i+1])*np.diff(xax_)[0]) - 2*0.341))\n",
"print(xax_[center_i], xax_[left_i+1], xax_[right_i])\n",
"\n",
"ax.fill_between(\n",
" x= xax_, \n",
" y1= np.exp(log_dens), \n",
" where= (xax_ >= xax_[left_i+1]) & (xax_ < xax_[right_i+1]),\n",
" color= \"b\",\n",
" alpha= 0.2)\n",
"\n",
"ax.axes.axvline(xax_[center_i]) \n",
"ax.set_xlabel('RDIS')\n",
"ax.set_ylabel('PDF')\n",
"\n",
"textstr = '\\n'.join((\n",
" r'$\\mu$={:.2f}'.format(xax_[center_i]),\n",
" r'$\\sigma_-$={:.2f}'.format(xax_[left_i+1]),\n",
" r'$\\sigma_+$={:.2f}'.format(xax_[right_i+1])))\n",
"\n",
"props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)\n",
"\n",
"# place a text box in upper left in axes coords\n",
"ax.text(0.05, 0.95, textstr, transform=ax.transAxes, fontsize=14,\n",
" verticalalignment='top', bbox=props)\n",
"\n",
"############################\n",
"fig, ax = plt.subplots()\n",
"\n",
"ax.plot(xax_, np.cumsum(np.exp(log_dens)*np.diff(xax_)[0]))\n",
"\n",
"ax.set_xlabel('RDIS')\n",
"ax.set_ylabel('CDF')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"input_data1 = df.ratio.values\n",
"input_data2 = df[[False if 'LFSS' in s else True for s in df.label.values]].ratio.values\n",
"\n",
"fig, ax = plt.subplots()\n",
"fig.suptitle('')\n",
"ax.hist(input_data1, bins=bins, density=True, alpha=0.5)\n",
"ax.hist(input_data2, bins=bins, density=True, alpha=0.5)\n",
"\n",
"ax.set_xlabel(\"RDIS\")\n",
"ax.set_ylabel(\"entries\")\n",
"\n",
"fig, ax = plt.subplots()\n",
"# fig.suptitle('Blue: all models, Red: without LFSS')\n",
"ax.set_xlabel('RDIS')\n",
"apply_KDE(input_data1, bounds=[bins[0], bins[-1]], data_in_relative_values=True)\n",
"apply_KDE(input_data2, bounds=[bins[0], bins[-1]], data_in_relative_values=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.10.6"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
},
"vscode": {
"interpreter": {
"hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}