Compare power datasets

# this script compares dataframes producing using the same antenna pattern but with various galactic radio emission models
# first you need to produce these, first run get_power_dataset.py in loop using all the galactic radio emission models
# e.g.: `for m in LFSS GSM08 GSM16 Haslam LFmap SSM GMOSS ULSA; do python get_power_dataset.py -g $m -s Salla_EW; done`

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import plotly.express as px
from pathlib import Path
import os
from scipy.interpolate import CubicSpline

# This ensures Plotly output works in multiple places:
# plotly_mimetype: VS Code notebook UI
# notebook: "Jupyter: Export to HTML" command in VS Code
# See https://plotly.com/python/renderers/#multiple-renderers
import plotly.io as pio

pio.renderers.default = "plotly_mimetype+notebook"

# some global plot settings
plt.rcParams["axes.labelweight"] = "bold"
plt.rcParams["font.weight"] = "bold"
plt.rcParams["font.size"] = 16
plt.rcParams["legend.fontsize"] = 14

plt.rcParams["xtick.major.width"] = 2
plt.rcParams["ytick.major.width"] = 2

plt.rcParams["xtick.major.size"] = 5
plt.rcParams["ytick.major.size"] = 5

plt.rcParams["xtick.labelsize"] = 14
plt.rcParams["ytick.labelsize"] = 14

layout_settings = dict(
    xaxis=dict(title="<b>LST</b>", tickprefix="<b>", ticksuffix="</b>", dtick=2),
    yaxis=dict(
        title="<b>frequency [MHz]</b>",
        tickprefix="<b>",
        ticksuffix="</b>",
        range=(30, 80),
        tick0=0,
        dtick=10,
        autorange=False,
    ),
    font=dict(
        size=20,
        color="black",
    ),
    coloraxis=dict(
        colorbar=dict(
            tickprefix="<b>",
            ticksuffix="</b>",
            title=dict(text="<b>Power [pW]</b>", side="right")),
        cmin=0,
        cmax=24,
    )
)

dir_path = "./simulated_power_datasets/"
df_files = [
    os.path.join(dir_path, i) for i in os.listdir(dir_path) if (i.endswith(".csv") & ('Salla_EW' in i))
]

df_list = []
df_names = []
for f in df_files:
    df = pd.read_csv(f, index_col=0)
    df.columns = df.columns.astype(float)
    df_list.append(df)
    df_names.append(Path(f).stem)

concatenated_df = pd.concat(df_list, keys=df_names)

# HINTS
# print df main keys
# print(concatenated_df.index.levels[0])
# get the same cell from each dfs
# one_cell_values = np.asarray([concatenated_df.xs(key).iloc[0,0] for key in concatenated_df.index.levels[0]])

# mean across the cells of all dfs
mean_df = concatenated_df.groupby(level=1).mean()
# std
std_df = concatenated_df.groupby(level=1).std()

fig = px.imshow(mean_df.T, width=600, aspect="cube", color_continuous_scale="jet")
fig.update_layout(**layout_settings)
fig.update_layout(
    title='<b>Mean</b>',
    coloraxis=dict(
        colorbar=dict(
            title=dict(
                text="<power> [pW]"
            )
        )
    )
)

fig = px.imshow(
    (std_df / mean_df.values).T, width=600, aspect="cube", color_continuous_scale="jet", title='<b>Relative standard deviation</b>'
)
fig.update_layout(**layout_settings)
fig.update_layout(
    title='<b>Mean</b>',
    coloraxis=dict(
        colorbar=dict(
            title=dict(
            text="\u03C1",
            side="top",
            )
        ),
        cmin=0,  # Minimum color value
        cmax=0.1,  # Maximum color value
    )
)

rho = (std_df / mean_df.values).values.flatten()
print('Min: {:.2f}, Max: {:.2f}, Mean: {:.2f}'.format(np.min(rho), np.max(rho), np.mean(rho)))

Min: 0.03, Max: 0.09, Mean: 0.06

fig, ax = plt.subplots()

ax.hist(rho)
ax.set_xlabel(r"$\rho$")
ax.set_ylabel("entries")

Text(0, 0.5, 'entries')

# Frequency mean of the mean dataframe
mean_signal = mean_df.mean(axis=1).values
mean_gal_signal_func = CubicSpline(
    np.append(mean_df.index.values, mean_df.index.values[-1] + 1),
    np.append(mean_signal, mean_signal[0]),
    extrapolate="periodic",
    bc_type="periodic",
)
new_x = np.linspace(0, 24, 24 * 60)
mean_gal_signal_func(new_x)

fig = px.line(
    x=new_x,
    y=mean_gal_signal_func(new_x),
    title="Frequency mean of the mean dataframe",
    width=600,
)

fig.add_trace(
    px.scatter(
        x=[
            new_x[np.argmin(mean_gal_signal_func(new_x))],
            new_x[np.argmax(mean_gal_signal_func(new_x))],
        ],
        y=[np.min(mean_gal_signal_func(new_x)), np.max(mean_gal_signal_func(new_x))],
    ).data[0]
)

fig.update_traces(marker=dict(color="red", size=15))
fig.update_layout(**layout_settings)
fig.update_layout(
    dict(
        yaxis=dict(
            title="<b>power [pW]</b>",
            tickprefix="<b>",
            ticksuffix="</b>",
            range=(5, 10),
            # tick0=0,
            dtick=1,
            autorange=False,
        ),
    )
)

fig.show()