{ "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 }