Prominence-Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import SAGA\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import scipy as sp\n",
    "from operator import itemgetter\n",
    "import os\n",
    "from dotenv import load_dotenv\n",
    "import matplotlib.colors as mcolors\n",
    "\n",
    "from SAGA import ObjectCuts as C\n",
    "\n",
    "load_dotenv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "custom_palette = ['#377eb8', '#ff7f00', '#4daf4a','#f781bf', '#a65628', '#984ea3','#999999', '#e41a1c', '#dede00']\n",
    "\n",
    "deut_pallete = [\"#FFEBAE\", \"#ffe599\",\"#f1e594\",\"#e1e592\",\"#d1e591\",\"#c0e492\",\"#9ae39a\",\"#85e2a0\",'#6de0a8',\"#50deb1\",\"#1cd2b8\",\"#00c6bd\",\"#00b9c0\",\"#00acc1\",\"#009fbf\",\"#0091ba\",\"#0083b3\",\"#0075a9\",\"#2a679c\",\"#295f97\",\"#295792\",\"#2b4f8d\",\"#2e4787\",\"#313f80\",\"#343679\",\"#372e72\",\"#3a246a\"]\n",
    "\n",
    "custom_style = {\n",
    "    \"figure.facecolor\": \"#000000\",\n",
    "    \"axes.facecolor\": \"#000000\",\n",
    "    \"savefig.facecolor\": \"#000000\", \n",
    "    \"grid.color\": \"#000000\",\n",
    "    \"text.color\": \"#ffe599\",\n",
    "    \"axes.labelcolor\": \"#ffe599\",\n",
    "    \"axes.edgecolor\": \"#ffe599\",\n",
    "    \"xtick.color\": \"#ffe599\",\n",
    "    \"ytick.color\": \"#ffe599\",\n",
    "    \"grid.linestyle\": \"-\",\n",
    "    \"lines.solid_capstyle\": \"round\",\n",
    "    \"xtick.color\": \"#ffe599\",\n",
    "    \"xtick.major.size\": \"4\",\n",
    "    \"xtick.minor.size\": \"4\",\n",
    "    \"xtick.major.pad\": \"4\",\n",
    "    \"xtick.minor.pad\": \"4\",\n",
    "    \"xtick.direction\": \"in\",\n",
    "    \"xtick.labelsize\": \"medium\",\n",
    "    \"ytick.direction\": \"in\",\n",
    "    \"ytick.labelsize\": \"medium\"\n",
    "}\n",
    "new_custom_palette = ['#9fc5e8', \"#ffe599\", \"#29af7fff\"]\n",
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data manipulations\n",
    "\n",
    "def get_radii(coordinates):\n",
    "    radii = []\n",
    "    for coordinate_pair in coordinates:\n",
    "        radii.append(np.sqrt(coordinate_pair[0]**2+coordinate_pair[1]**2))\n",
    "    return radii\n",
    "    \n",
    "def get_ellipticities(coordinates):\n",
    "    # # First get quadrupole moments\n",
    "    coordinates_squared = np.square(coordinates)\n",
    "\n",
    "    xdiff_ydiff_col = np.prod(coordinates, axis=2)\n",
    "    xdiff_ydiff_col = xdiff_ydiff_col[:,:,np.newaxis]\n",
    "\n",
    "    quadrupole_moments_setup = np.insert(coordinates_squared, 0, [[[xdiff_ydiff_col]]], axis=2) # diffx*diffy now at position 0\n",
    "    quadrupole_moments = quadrupole_moments_setup.mean(axis=1)\n",
    "\n",
    "    # # Get ellipticity from quadrupole moments\n",
    "    # Q_xy, Q_xx, Q_yy = np.quadrupole_moments\n",
    "    quadrupole_xx_yy = quadrupole_moments[:,1:3]\n",
    "    quadrupole_xy = quadrupole_moments[:,0:1]\n",
    "\n",
    "    e_1_numerator = -np.diff(quadrupole_xx_yy, axis=1)\n",
    "    e_2_numerator = 2*quadrupole_xy\n",
    "\n",
    "    e_denominator_1 = np.sum(quadrupole_xx_yy, axis=1)[:,np.newaxis]\n",
    "    e_denominator_inside_sqrt = np.prod(quadrupole_xx_yy, axis=1)[:,np.newaxis] - quadrupole_xy**2\n",
    "    e_denominator = e_denominator_1 + 2*np.sqrt(e_denominator_inside_sqrt)\n",
    "\n",
    "    e_1 = e_1_numerator / e_denominator\n",
    "    e_2 = e_2_numerator / e_denominator\n",
    "\n",
    "    ellipticities = np.sqrt(e_1**2 + e_2**2)[:,0]\n",
    "\n",
    "    return ellipticities\n",
    "\n",
    "\n",
    "def get_randomized_ellipticities(radii, num_iterations):\n",
    "    num = len(radii)\n",
    "    radial_stdev = np.std(radii)\n",
    "\n",
    "    # Create \"isotropic\" (random from uniform - not uniform) angles\n",
    "    coords = np.random.randn(num_iterations, num, 2)  # (number of times to repeat data set generation, how many satellites, how many spatial dimensions)\n",
    "    calc_radii = np.sqrt((coords**2).sum(axis=2))\n",
    "    normalized_coords = coords/calc_radii[:,:,np.newaxis]\n",
    "\n",
    "    # Move the positions of satellites out or in radially until they reach original radial scale\n",
    "    reshaped_radii = np.reshape(radii, (1,len(radii), 1))\n",
    "    scaled_coords = normalized_coords*reshaped_radii\n",
    "\n",
    "    ellipticities = get_ellipticities(scaled_coords,)\n",
    "    \n",
    "    return ellipticities\n",
    "\n",
    "\n",
    "def get_prominence(actual, radii, num_iterations):\n",
    "    randomized = get_randomized_ellipticities(radii, num_iterations)\n",
    "    abv_actual = np.where(randomized >= actual, 1, 0)\n",
    "    p_above_actual = np.sum(abv_actual)/num_iterations\n",
    "    prominence = 1/p_above_actual\n",
    "    \n",
    "    return prominence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dataframe manipulations\n",
    "def filter_subhalo_dataframe_by_num_satellites(dataframe, cutoff):\n",
    "    filtered_dataframe = dataframe[dataframe['vmax_mpeak'] >= cutoff].copy()\n",
    "    filtered_dataframe = filtered_dataframe.groupby(['HOSTID']).filter(lambda x: len(x) > 2).copy()\n",
    "    return filtered_dataframe\n",
    "\n",
    "def add_coordinate_columns(dataframe):\n",
    "    dataframe['x_coord'] = dataframe['x_adj'] - dataframe['x_host']\n",
    "    dataframe['y_coord'] = dataframe['y_adj'] - dataframe['y_host']\n",
    "    # dataframe['z_coord'] = dataframe['z_adj'] - dataframe['z_host']\n",
    "    # dataframe['3d_coordinates'] = dataframe[['x_coord', 'y_coord', 'z_coord']].values.tolist()\n",
    "    dataframe['2d_coordinates'] = dataframe[['x_coord', 'y_coord']].values.tolist()\n",
    "\n",
    "def get_systems_dataframe(dataframe):\n",
    "    systems_df = dataframe.groupby('HOSTID')['2d_coordinates'].apply(tuple).reset_index(name='coordinates_list')\n",
    "    systems_df['coordinates_list'] = systems_df['coordinates_list'].tolist()\n",
    "    systems_df['radii'] = systems_df['coordinates_list'].apply(get_radii)\n",
    "    systems_df['Dataset'] = 'Actual'\n",
    "    return systems_df\n",
    "\n",
    "def add_ellipticity_column(dataframe):\n",
    "    dataframe['ellipticity'] = dataframe.apply(lambda x: get_ellipticities([x['coordinates_list']])[0], axis=1)\n",
    "\n",
    "def get_isotropic_data(dataframe):\n",
    "    dataframe_copy = dataframe.copy()\n",
    "    dataframe_copy['ellipticity'] = dataframe_copy.apply(lambda x: get_randomized_ellipticities(x['radii'], 1)[0], axis=1)\n",
    "    dataframe_copy['Dataset'] = 'Isotropic'\n",
    "    return dataframe_copy\n",
    "\n",
    "def add_num_subs_column(dataframe):\n",
    "    dataframe['num_subs'] = dataframe['radii'].str.len()\n",
    "    display(dataframe)\n",
    "\n",
    "def add_prominence_column(dataframe, num_iterations):\n",
    "    dataframe['prominence'] = dataframe.apply(lambda x: get_prominence(x['ellipticity'], x['radii'], num_iterations), axis=1)\n",
    "\n",
    "def get_prominence_dataframe(dataframe):\n",
    "    prominence_data = dataframe[['HOSTID', 'ellipticity', 'prominence', 'Dataset']].copy()\n",
    "    return prominence_data\n",
    "\n",
    "def plot_prominence(data, num_iterations):\n",
    "    format = {\n",
    "        \"figure.facecolor\": \"212946\",\n",
    "        \"axes.facecolor\": \"212946\",\n",
    "        \"savefig.facecolor\": \"212946\", \n",
    "        \"grid.color\": \"2A3459\",\n",
    "        \"text.color\": \"0.9\",\n",
    "        \"axes.labelcolor\": \"0.9\",\n",
    "        \"xtick.color\": \"0.9\",\n",
    "        \"ytick.color\": \"0.9\",\n",
    "        \"grid.linestyle\": \"-\",\n",
    "        \"lines.solid_capstyle\": \"round\",\n",
    "        \"figure.figsize\" : \"(11.7,8.27)\"\n",
    "    }\n",
    "\n",
    "    sns.set_style(\"darkgrid\", format)\n",
    "\n",
    "    fig, ax1 = plt.subplots()\n",
    "    sp1 = sns.ecdfplot(data=data, x=\"prominence\", hue=\"Dataset\", hue_order=['Actual', 'Isotropic'], complementary=True, ax=ax1)\n",
    "    ax1.set_xscale('log')\n",
    "    ax1.set_yscale('log')\n",
    "    ax1.set_xlim(1,num_iterations)\n",
    "    plt.show()\n",
    "\n",
    "def plot_prominence_histogram(data, num_iterations):\n",
    "    format = {\n",
    "        \"figure.facecolor\": \"212946\",\n",
    "        \"axes.facecolor\": \"212946\",\n",
    "        \"savefig.facecolor\": \"212946\", \n",
    "        \"grid.color\": \"2A3459\",\n",
    "        \"text.color\": \"0.9\",\n",
    "        \"axes.labelcolor\": \"0.9\",\n",
    "        \"xtick.color\": \"0.9\",\n",
    "        \"ytick.color\": \"0.9\",\n",
    "        \"grid.linestyle\": \"-\",\n",
    "        \"lines.solid_capstyle\": \"round\",\n",
    "        \"figure.figsize\" : \"(11.7,8.27)\"\n",
    "    }\n",
    "\n",
    "    sns.set_style(\"darkgrid\", format)\n",
    "\n",
    "    fig, ax1 = plt.subplots()\n",
    "    sp1 = sns.histplot(data=data, x=\"prominence\", hue=\"Dataset\", hue_order=['Actual', 'Isotropic'], bins=7, element=\"poly\", stat=\"density\", common_norm=False, common_bins=True, fill=False, ax=ax1, palette='gist_rainbow')\n",
    "    ax1.set_xscale('log')\n",
    "    # ax1.set_yscale('log')\n",
    "    # ax1.set_xlim(1,num_iterations)\n",
    "    plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ellipticity_and_prominence_data(dataframe, num_iterations, num_isotropic):\n",
    "    add_coordinate_columns(dataframe)\n",
    "    systems_dataframe = get_systems_dataframe(dataframe).copy()\n",
    "    add_ellipticity_column(systems_dataframe) # this is the original dataframe\n",
    "    add_prominence_column(systems_dataframe, num_iterations) # bottleneckyy\n",
    "    prominence_data = get_prominence_dataframe(systems_dataframe).copy()\n",
    "\n",
    "    for i in range (0,num_isotropic):\n",
    "        systems_dataframe_isotropic = get_isotropic_data(systems_dataframe).copy()\n",
    "        add_prominence_column(systems_dataframe_isotropic, num_iterations)  \n",
    "        prominence_data_isotropic = get_prominence_dataframe(systems_dataframe_isotropic).copy()\n",
    "        prominence_data = pd.concat([prominence_data, prominence_data_isotropic], ignore_index=True)\n",
    "\n",
    "    return [systems_dataframe, prominence_data]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%store -r saga_sats\n",
    "%store -r VSMDPL_subs_raw\n",
    "%store -r VSMDPL_subs_interlopers\n",
    "\n",
    "VSMDPL_subs_raw = VSMDPL_subs_raw.copy()\n",
    "saga_sats = saga_sats.copy()\n",
    "VSMDPL_subs_interlopers = VSMDPL_subs_interlopers.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(VSMDPL_subs_interlopers.interloper_host_id.unique()))\n",
    "print(len(VSMDPL_subs_raw.upid.unique()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HOST_PROPERTIES = os.getenv('HOST_PROPERTIES')\n",
    "saga_hosts = pd.read_parquet(HOST_PROPERTIES)\n",
    "\n",
    "iterations = 10000\n",
    "\n",
    "saga_sats.rename(columns = {'X':'x_adj', 'Y':'y_adj', 'HOST_X':'x_host', 'HOST_Y':'y_host'}, inplace = True)\n",
    "VSMDPL_subs_raw.rename(columns = {'upid':'HOSTID'}, inplace = True)\n",
    "VSMDPL_subs_interlopers.rename(columns = {'interloper_host_id':'HOSTID', 'interloper_host_x':'x_host', 'interloper_host_y':'y_host', 'interloper_host_vmax_mpeak':'host_vmax_mpeak'}, inplace = True)\n",
    "\n",
    "# Use vmax_mpeak as proxy for luminosity\n",
    "vmax_cutoff_for_avg_num = {\n",
    "    '3.5' : 61,\n",
    "      '4' : 51.6,\n",
    "    '4.5' : 46.4,\n",
    "      '5' : 43.5,\n",
    "    '5.5' : 40.98,\n",
    "    'saga': 44.09, # 4.855 sats per host on average\n",
    "    'saga_interlopers': 46.77 # 4.855 sats per host on average\n",
    "}\n",
    "\n",
    "VSMDPL_subs_multi = {}\n",
    "for cutoff in list(vmax_cutoff_for_avg_num.keys())[:5]:\n",
    "    VSMDPL_subs_multi[cutoff] = filter_subhalo_dataframe_by_num_satellites(VSMDPL_subs_raw, vmax_cutoff_for_avg_num[cutoff])\n",
    "\n",
    "VSMDPL_subs = filter_subhalo_dataframe_by_num_satellites(VSMDPL_subs_raw, vmax_cutoff_for_avg_num['saga'])\n",
    "VSMDPL_subs_interlopers = filter_subhalo_dataframe_by_num_satellites(VSMDPL_subs_interlopers, vmax_cutoff_for_avg_num['saga_interlopers'])\n",
    "\n",
    "saga_sats = saga_sats.groupby(['HOSTID']).filter(lambda x: len(x) > 2).copy()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SAGA Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n",
    "saga_data, saga_plot = get_ellipticity_and_prominence_data(saga_sats, iterations, 1)\n",
    "plot_prominence(saga_plot, iterations)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate 100x isotropic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga_data, saga_plot = get_ellipticity_and_prominence_data(saga_sats, iterations, 100)\n",
    "plot_prominence(saga_plot, iterations)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SIM Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VSMDPL_data, VSMDPL_plot = get_ellipticity_and_prominence_data(VSMDPL_subs, iterations, 1)\n",
    "plot_prominence(VSMDPL_plot, iterations)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SIM Interloper Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VSMDPL_interloper_data, VSMDPL_interloper_plot = get_ellipticity_and_prominence_data(VSMDPL_subs_interlopers, iterations, 1)\n",
    "plot_prominence(VSMDPL_interloper_plot, iterations)\n",
    "VSMDPL_interloper_data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SAGA + Sim (Halofinder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga_plot['Datasource'] = 'Saga'\n",
    "VSMDPL_plot['Datasource'] = 'VSMDPL (RS)'\n",
    "saga_sim_data = pd.concat([saga_plot, VSMDPL_plot], ignore_index=True)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "sp1 = sns.ecdfplot(data=saga_sim_data, x=\"prominence\", hue=saga_sim_data[['Datasource', 'Dataset']].apply(tuple, axis=1), complementary=True, ax=ax1)\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_xlim(1,10000)\n",
    "\n",
    "lss = [':', '--', '-.', '-']\n",
    "\n",
    "handles = sp1.legend_.legendHandles[::-1]\n",
    "\n",
    "for line, ls, handle in zip(sp1.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "\n",
    "fig.set_size_inches([6, 4])\n",
    "fig.savefig(\"figures/SAGAvVSMDPLProminence.png\",dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_normal_df = pd.DataFrame(list(np.arange(1,1000001,1)), columns=['ordered_list'])\n",
    "random_normal_df['prominence'] = random_normal_df.size/(random_normal_df.size - random_normal_df.iloc[lambda x: x.index]+1)\n",
    "random_normal_df['Datasource'] = 'Isotropic'\n",
    "random_normal_df"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SAGA + Sim (interlopers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n",
    "# random = pd.DataFrame(list(prominence_arr), columns=['prominence'])\n",
    "random_normal_df['Dataset'] = 'Actual'\n",
    "saga_plot['Datasource'] = 'Saga'\n",
    "VSMDPL_interloper_plot['Datasource'] = 'VSMDPL'\n",
    "saga_sim_data = pd.concat([saga_plot, VSMDPL_interloper_plot, random_normal_df], ignore_index=True)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "sp1 = sns.ecdfplot(data=saga_sim_data[saga_sim_data['Dataset'] == 'Actual'], x=\"prominence\", hue=saga_sim_data['Datasource'], complementary=True, ax=ax1, palette=new_custom_palette)\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_xlim(1,10000)\n",
    "ax1.set_ylim(0.0001,1)\n",
    "\n",
    "saga_lim = ax1.axhline(y=1/len(list(saga_plot.HOSTID.unique())), color=new_custom_palette[0], linestyle='-', alpha=0.5)\n",
    "VSMDPL_lim = ax1.axhline(y=1/len(list(VSMDPL_interloper_plot.HOSTID.unique())), color=new_custom_palette[1], linestyle='-', alpha=0.5)\n",
    "\n",
    "lss = [':', '-.', '-']\n",
    "\n",
    "handles = sp1.legend_.legendHandles[::-1]\n",
    "ax1.legend_.set_title('')\n",
    "ax1.set_xlabel('Prominence $\\mathcal{P}_e$')\n",
    "ax1.set_ylabel('Fraction > $\\mathcal{P}_e$')\n",
    "sp1.set_xticks([1, 10, 100, 1000, 10000])\n",
    "sp1.tick_params(which='both', direction=\"in\")\n",
    "\n",
    "for line, ls, handle in zip(sp1.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "fig.set_size_inches([5, 4])\n",
    "fig.savefig(\"figures/SAGAvVSMDPLProminence.png\",dpi=600)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sim (Halofinder) + Sim (Interlopers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_sim_data = pd.concat([VSMDPL_plot, VSMDPL_interloper_plot], ignore_index=True)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "sp1 = sns.ecdfplot(data=sim_sim_data[sim_sim_data['Dataset'] == 'Actual'], x=\"prominence\", hue=sim_sim_data['Datasource'], complementary=True, ax=ax1)\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_xlim(1,10000)\n",
    "\n",
    "lss = [':', '--', '-.', '-']\n",
    "\n",
    "handles = sp1.legend_.legendHandles[::-1]\n",
    "\n",
    "for line, ls, handle in zip(sp1.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "\n",
    "fig.set_size_inches([6, 4])\n",
    "fig.savefig(\"figures/VSMDPLvVSMDPLProminence.png\",dpi=600)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saga + Sim + sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"white\")\n",
    "\n",
    "saga_sim_sim_data = pd.concat([saga_plot, sim_sim_data], ignore_index=True)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "sp1 = sns.ecdfplot(data=saga_sim_sim_data[saga_sim_sim_data['Dataset'] == 'Actual'], x=\"prominence\", hue=saga_sim_sim_data['Datasource'], complementary=True, ax=ax1)\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_xlim(1,10000)\n",
    "ax1.set_xlabel('Prominence $\\mathcal{P}_e$')\n",
    "ax1.set_ylabel('Density')\n",
    "\n",
    "\n",
    "lss = [':', '--', '-.', '-']\n",
    "\n",
    "handles = sp1.legend_.legend_handles[::-1]\n",
    "ax1.legend_.set_title('')\n",
    "\n",
    "for line, ls, handle in zip(sp1.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "fig.set_size_inches([8, 5])\n",
    "fig.savefig(\"figures/VSMDPLvVSMDPLProminence.png\",dpi=600)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "fig2, ax2 = plt.subplots()\n",
    "sp2 = sns.histplot(data=saga_sim_sim_data[saga_sim_sim_data['Dataset'] == 'Actual'], x=\"prominence\", hue=\"Datasource\", bins=50, element=\"poly\", stat=\"percent\", common_norm=False, common_bins=True, fill=True, ax=ax2, palette='gist_rainbow', log_scale=[True,True])\n",
    "\n",
    "lss = ['--', '-', '-.']\n",
    "\n",
    "handles = ax2.legend_.legendHandles[::-1]\n",
    "ax2.legend_.set_title('')\n",
    "\n",
    "ax2.set_xlim(1,10000)\n",
    "\n",
    "for line, ls, handle in zip(sp1.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "fig2.set_size_inches([8, 5])\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ellipticity SAGA + Sim (Halofinder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n",
    "saga_data['Datasource'] = 'Saga'\n",
    "VSMDPL_data['Datasource'] = 'VSMDPL (RS)'\n",
    "saga_sim_ellipticity_data = pd.concat([saga_data, VSMDPL_data], ignore_index=True)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_xlabel('$\\longleftarrow$ more spherical         Ellipticity (2D)         more elliptical $\\longrightarrow$')\n",
    "sp1 = sns.histplot(saga_sim_ellipticity_data, x=\"ellipticity\", hue=\"Datasource\", hue_order=['Saga', 'VSMDPL (RS)'], bins=7, element=\"poly\", stat=\"density\", common_norm=False, common_bins=True, fill=False, ax=ax1, palette='gist_rainbow')\n",
    "\n",
    "lss = ['--', '-']\n",
    "\n",
    "handles = sp1.legend_.legendHandles[::-1]\n",
    "\n",
    "for line, ls, handle in zip(sp1.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "\n",
    "fig.set_size_inches([6, 4])\n",
    "fig.savefig(\"figures/SAGAvVSMDPLEllipticity.png\",dpi=600)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ellipticity SAGA + Sim (Interlopers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n",
    "\n",
    "saga_data['Datasource'] = 'Saga'\n",
    "VSMDPL_interloper_data['Datasource'] = 'VSMDPL'\n",
    "saga_sim_ellipticity_data = pd.concat([saga_data, VSMDPL_interloper_data], ignore_index=True)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_xlabel('$\\longleftarrow$ circular              Ellipticity (2D)                 linear $\\longrightarrow$')\n",
    "sp1 = sns.histplot(saga_sim_ellipticity_data, x=\"ellipticity\", hue=\"Datasource\", hue_order=['Saga'], bins=11, element=\"poly\", stat=\"density\", common_norm=False, common_bins=True, fill=False, ax=ax1, palette=['#9fc5e8', \"#ffe599\"])\n",
    "sp2 = sns.histplot(saga_sim_ellipticity_data, x=\"ellipticity\", hue=\"Datasource\", hue_order=['VSMDPL'], bins=20, element=\"poly\", stat=\"density\", common_norm=False, common_bins=True, fill=False, ax=ax1, palette=['#ffe599', \"#9fc5e8\"])\n",
    "sp1.set_xticks([0, .2, .4, .6, .8, 1])\n",
    "sp1.set_yticks([0, .4, .8, 1.2, 1.6, 2.0])\n",
    "sp1.set_yticklabels(['', .4, .8, 1.2, 1.6, 2.0])\n",
    "sp1.tick_params(which='both', direction=\"in\")\n",
    "\n",
    "lss = ['--', '-']\n",
    "\n",
    "handles = sp1.legend_.legendHandles[::-1]\n",
    "\n",
    "\n",
    "for line, ls, handle in zip(sp2.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "ax1.legend_.set_title('')\n",
    "ax1.legend(title='', loc='upper right', labels=['SAGA', 'VSMDPL'])\n",
    "\n",
    "\n",
    "fig.set_size_inches([5, 4])\n",
    "fig.savefig(\"figures/SAGAvVSMDPLEllipticity.png\",dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Saga mean ellipticity: {}\".format(saga_data.ellipticity.mean()))\n",
    "print(\"Sim mean ellipticity: {}\".format(VSMDPL_interloper_data.ellipticity.mean()))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ellipticity Sim + Sim (plus prominence for fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n",
    "## Ellipticity Plot\n",
    "VSMDPL_interloper_data['Datasource'] = 'VSMDPL (Int)'\n",
    "VSMDPL_data['Datasource'] = 'VSMDPL (RS)'\n",
    "saga_sim_ellipticity_data = pd.concat([VSMDPL_data, VSMDPL_interloper_data], ignore_index=True)\n",
    "\n",
    "fig, axs = plt.subplots(ncols=2, layout=\"constrained\")\n",
    "axs[0].set_xlabel('$ \\longleftarrow$ circular   Ellipticity (2D)   linear $\\longrightarrow$ ')\n",
    "sp1 = sns.histplot(saga_sim_ellipticity_data, x=\"ellipticity\", hue=\"Datasource\", hue_order=['VSMDPL (RS)', 'VSMDPL (Int)'], bins=7, element=\"poly\", stat=\"density\", common_norm=False, common_bins=True, fill=False, ax=axs[0], palette=['#9fc5e8', \"#ffe599\"])\n",
    " \n",
    "\n",
    "lss = ['--', '-']\n",
    "\n",
    "handles = ax1.legend_.legendHandles[::-1]\n",
    "axs[0].legend_.set_title('')\n",
    "axs[0].tick_params(which='both', direction=\"in\")\n",
    "axs[0].set_yticks([0, .25, .5, .75, 1.0, 1.25, 1.5, 1.75, 2.0])\n",
    "axs[0].set_yticklabels(['', .25, .50, .75, 1.00, 1.25, 1.50, 1.75, 2.00])\n",
    "axs[0].set_xticks([0, .2, .4, 0.6, 0.8, 1])\n",
    "\n",
    "for line, ls, handle in zip(sp1.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "\n",
    "## Prominence Plot\n",
    "sns.set_style(\"white\")\n",
    "VSMDPL_interloper_plot['Datasource'] = 'VSMDPL (Int)'\n",
    "sim_sim_data = pd.concat([VSMDPL_plot, VSMDPL_interloper_plot], ignore_index=True)\n",
    "\n",
    "sp2 = sns.ecdfplot(data=sim_sim_data[sim_sim_data['Dataset'] == 'Actual'], x=\"prominence\", hue=sim_sim_data['Datasource'], complementary=True, ax=axs[1], palette=['#9fc5e8', \"#ffe599\"])\n",
    "axs[1].set_xscale('log')\n",
    "axs[1].set_yscale('log')\n",
    "axs[1].set_xlim(1,10000)\n",
    "axs[1].set_xlabel('Prominence $\\mathcal{P}_e$')\n",
    "axs[1].set_ylabel('Proportion', labelpad=0)\n",
    "axs[1].legend_.set_title('')\n",
    "axs[1].tick_params(which='both', direction=\"in\")\n",
    "\n",
    "\n",
    "\n",
    "lss = ['--', '-']\n",
    "\n",
    "handles = sp2.legend_.legendHandles[::-1]\n",
    "\n",
    "for line, ls, handle in zip(sp2.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "\n",
    "fig.set_size_inches([7, 3])\n",
    "fig.savefig(\"figures/VSMDPLvVSMDPLProminenceAndEllipticity.png\",dpi=600)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ellipticity SAGA + Sim rockstar + Sim interlopers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VSMDPL_interloper_data['Datasource'] = 'VSMDPL (Int)'\n",
    "saga_data['Datasource'] = 'SAGA'\n",
    "saga_sim_ellipticity_data = pd.concat([saga_data, VSMDPL_data, VSMDPL_interloper_data], ignore_index=True)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_xlabel('$\\longleftarrow$ more spherical         Ellipticity (2D)         more elliptical $\\longrightarrow$')\n",
    "sp1 = sns.histplot(saga_sim_ellipticity_data, x=\"ellipticity\", hue=\"Datasource\", hue_order=['SAGA', 'VSMDPL (RS)', 'VSMDPL (Int)'], bins=7, element=\"poly\", stat=\"density\", common_norm=False, common_bins=True, fill=False, ax=ax1)\n",
    " \n",
    "\n",
    "lss = ['--', '-', '-.']\n",
    "\n",
    "handles = ax1.legend_.legendHandles[::-1]\n",
    "ax1.legend_.set_title('')\n",
    "\n",
    "for line, ls, handle in zip(sp1.lines, lss, handles):\n",
    "    line.set_linestyle(ls)\n",
    "    handle.set_ls(ls)\n",
    "\n",
    "\n",
    "fig.set_size_inches([5, 4])\n",
    "fig.savefig(\"figures/SAGAvVSMDPLvInterloperEllipticity.png\",dpi=600)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SIM Compare prominence for different avg # satellites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VSMDPL_systems_multi = {}\n",
    "VSMDPL_plot_multi = {}\n",
    "\n",
    "for avg in list(vmax_cutoff_for_avg_num.keys())[:5]:\n",
    "    VSMDPL_systems_multi[avg], VSMDPL_plot_multi[avg] = get_ellipticity_and_prominence_data(VSMDPL_subs_multi[avg], iterations, 1)\n",
    "    VSMDPL_plot_multi[avg]['Average # Satellites'] = avg\n",
    "    VSMDPL_systems_multi[avg]['Datasource'] = 'VSMDPL ' + avg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VSMDPL_multi_plot_data = pd.concat(VSMDPL_plot_multi, ignore_index=True)\n",
    "VSMDPL_multi_plot_data = VSMDPL_multi_plot_data[VSMDPL_multi_plot_data['Dataset'] == 'Actual'].copy()\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "sp1 = sns.ecdfplot(data=VSMDPL_multi_plot_data, x=\"prominence\", hue='Average # Satellites', hue_order=['3.5', '4', '4.5', '5', '5.5'], complementary=True, ax=ax1)\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_xlim(1,iterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for df in VSMDPL_systems_multi.values():\n",
    "    add_num_subs_column(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n",
    "VSMDPL_multi_data = pd.concat(VSMDPL_systems_multi, ignore_index=True)\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_xlabel('$\\longleftarrow$ circular                  Ellipticity (2D)                       linear $\\longrightarrow$')\n",
    "sp1 = sns.histplot(VSMDPL_multi_data, x=\"ellipticity\", hue=\"Datasource\", hue_order=['VSMDPL 3.5', 'VSMDPL 4', 'VSMDPL 4.5', 'VSMDPL 5', 'VSMDPL 5.5'], bins=7, element=\"poly\", stat=\"density\", common_norm=False, common_bins=True, fill=False, ax=ax1, palette='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_data = pd.concat([VSMDPL_multi_data, saga_sim_ellipticity_data], ignore_index=True)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_xlabel('$\\longleftarrow$ circular              Ellipticity (2D)                 linear $\\longrightarrow$')\n",
    "sp1 = sns.histplot(all_data, x=\"ellipticity\", hue=\"Datasource\", hue_order=['Saga', 'VSMDPL', 'VSMDPL 3.5', 'VSMDPL 4', 'VSMDPL 4.5', 'VSMDPL 5', 'VSMDPL 5.5'], bins=7, element=\"poly\", stat=\"density\", common_norm=False, common_bins=True, fill=False, ax=ax1, palette=custom_palette)\n",
    "\n",
    "\n",
    "fig.set_size_inches([7, 3])\n",
    "\n",
    "\n",
    "fig.savefig(\"figures/VSMDPLNumSatellitesEllipticity.png\",dpi=600)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### VSMDPL Interloper cuts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VSMDPL_interloper_data, VSMDPL_interloper_plot = get_ellipticity_and_prominence_data(VSMDPL_subs_interlopers, iterations, 1)\n",
    "plot_prominence(VSMDPL_interloper_plot, iterations)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ellipticity and Prominence vs Host properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "add_num_subs_column(saga_data)\n",
    "saga = saga_hosts.merge(saga_data, on=\"HOSTID\")\n",
    "saga = saga.rename(columns={\"num_subs\":\"num_sats\"})\n",
    "saga_hosts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga['num_sats'] = saga['num_sats'].astype(np.float64)\n",
    "saga.sort_values('prominence')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n",
    "sns.scatterplot(saga, x='ellipticity', y='prominence', edgecolor=None, hue=saga['num_sats'], palette='cividis', legend='full')\n",
    "plt.yscale('log')\n",
    "plt.xlim([0,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isotropic_saga_data = saga_plot[saga_plot['Dataset']=='Isotropic'].copy()\n",
    "isotropic_saga_data = isotropic_saga_data.merge(saga_data[['HOSTID', 'radii']], on='HOSTID').copy()\n",
    "add_num_subs_column(isotropic_saga_data)\n",
    "isotropic_saga_data = isotropic_saga_data.rename(columns={\"num_subs\":\"num_sats\"})\n",
    "\n",
    "\n",
    "sns.scatterplot(isotropic_saga_data, x='ellipticity', y='prominence', edgecolor=None, hue=isotropic_saga_data['num_sats'], palette='viridis', legend='full', size=1)\n",
    "plt.yscale('log')\n",
    "plt.xlim([0,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga_copy = saga[['HOSTID', 'ellipticity', 'prominence', 'Dataset', 'radii', 'num_sats']]\n",
    "raw_and_isotropic_saga_data = pd.concat([saga_copy, isotropic_saga_data])\n",
    "raw_and_isotropic_saga_data.loc[raw_and_isotropic_saga_data[\"Dataset\"] == \"Actual\", \"edgecolor\"] = 'white'\n",
    "raw_and_isotropic_saga_data.loc[raw_and_isotropic_saga_data[\"Dataset\"] == \"Isotropic\", \"edgecolor\"] = 'None'\n",
    "raw_and_isotropic_saga_data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "sns.scatterplot(raw_and_isotropic_saga_data[raw_and_isotropic_saga_data['Dataset']=='Isotropic'], x='ellipticity', y='prominence', hue='num_sats', edgecolor=None, palette='viridis', legend=None, ax=ax, size=0.1)\n",
    "sns.scatterplot(raw_and_isotropic_saga_data[raw_and_isotropic_saga_data['Dataset']=='Actual'], x='ellipticity', y='prominence', hue='num_sats', edgecolor='white', palette='viridis', legend=None, ax=ax, linewidth=1)\n",
    "plt.yscale('log')\n",
    "plt.xlim([0,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(ncols=4)\n",
    "fig.set_size_inches(30, 6, forward=True)\n",
    "sns.scatterplot(data=saga, x=\"MORPHTYPE\", y=\"ellipticity\", edgecolor=None, ax=axs[0])\n",
    "sns.scatterplot(data=saga, x=\"radius\", y=\"ellipticity\", edgecolor=None, ax=axs[1])\n",
    "sns.scatterplot(data=saga, x=\"ba\", y=\"ellipticity\", edgecolor=None, ax=axs[2])\n",
    "sns.scatterplot(data=saga, x=\"phi\", y=\"ellipticity\", edgecolor=None, ax=axs[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(ncols=4)\n",
    "fig.set_size_inches(30, 6, forward=True)\n",
    "sns.scatterplot(data=saga, x=\"sb_r\", y=\"ellipticity\", ax=axs[0])\n",
    "sns.scatterplot(data=saga, x=\"SERSIC\", y=\"ellipticity\", ax=axs[1])\n",
    "sns.scatterplot(data=saga, x=\"K_ABS\", y=\"ellipticity\", ax=axs[2])\n",
    "sns.scatterplot(data=saga, x=\"gr\", y=\"ellipticity\", ax=axs[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spearmanrs = {}\n",
    "props = [\"MORPHTYPE\", \"radius\", \"ba\", \"phi\", \"sb_r\", \"SERSIC\", \"K_ABS\", \"gr\", \"num_sats\"]\n",
    "for prop in props:\n",
    "    spearmanr_results = sp.stats.spearmanr(saga['prominence'], saga[prop], axis=0, nan_policy='propagate', alternative='two-sided')\n",
    "    spearmanrs[prop] = spearmanr_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.stats.spearmanr(np.array(saga['phi']), np.array(saga['prominence']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots()\n",
    "plt.scatter(x=np.array(saga['phi']), y=np.array(saga['prominence']))\n",
    "print(sp.stats.spearmanr(np.array(saga['phi']), np.array(saga['prominence'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(ncols=4)\n",
    "fig.set_size_inches(30, 6, forward=True)\n",
    "for i in range (0,4):\n",
    "    statistic = spearmanrs[props[i]].statistic.round(2)\n",
    "    pvalue = spearmanrs[props[i]].pvalue.round(2)\n",
    "    spearman_string = \"Coefficient: \" + str(statistic) + \" P-value: \" + str(pvalue)\n",
    "    sns.scatterplot(data=saga, x=props[i], y=\"prominence\", ax=axs[i])\n",
    "    axs[i-4].text(0.99, 0.99, spearman_string, ha='right', va='top', size=20, transform = axs[i].transAxes)\n",
    "axs[1].set_xscale('log') # just for radius\n",
    "for ax in axs:\n",
    "    ax.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(ncols=4)\n",
    "fig.set_size_inches(30, 6, forward=True)\n",
    "for i in range (4,8):\n",
    "    statistic = spearmanrs[props[i]].statistic.round(2)\n",
    "    pvalue = spearmanrs[props[i]].pvalue.round(2)\n",
    "    spearman_string = \"Coefficient: \" + str(statistic) + \" P-value: \" + str(pvalue)\n",
    "    sns.scatterplot(data=saga, x=props[i], y=\"prominence\", ax=axs[i-4])\n",
    "    axs[i-4].text(0.99, 0.99, spearman_string,\n",
    "        ha='right',\n",
    "        va='top',\n",
    "        size=20,\n",
    "        transform = axs[i-4].transAxes)\n",
    "for ax in axs:\n",
    "    ax.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.set_style('ticks')\n",
    "sns.set_style(\"ticks\", custom_style)\n",
    "\n",
    "fig, axs = plt.subplots(ncols=2, sharey=True)\n",
    "fig.set_size_inches(8, 3, forward=True)\n",
    "\n",
    "figs = [1, 6]\n",
    "fig_labels = ['$r_{host}$ (kpc)', '$K_{abs}$']\n",
    "\n",
    "for i in figs: # update w/ enumerate\n",
    "    i_index = figs.index(i)\n",
    "    statistic = spearmanrs[props[i]].statistic.round(2)\n",
    "    pvalue = spearmanrs[props[i]].pvalue.round(2)\n",
    "    spearman_string = '$\\\\rho_S$: ' + str(statistic) + \", P: \" + str(pvalue)\n",
    "    sns.scatterplot(data=saga, x=props[i], y=\"prominence\", ax=axs[i_index])\n",
    "    axs[i_index].text(0.8, 1.14, spearman_string,\n",
    "        ha='right',\n",
    "        va='top',\n",
    "        size=12,\n",
    "        transform = axs[i_index].transAxes)\n",
    "    axs[i_index].set_xlabel(fig_labels[i_index], fontsize=14)\n",
    "for ax in axs:\n",
    "    ax.set_yscale('log')\n",
    "axs[0].set_ylabel('Prominence $\\mathcal{P}_e$', fontsize=12)\n",
    "# axs[0].set_xticks([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.subplots_adjust(wspace=0, hspace=0)\n",
    "fig.savefig(\"figures/ProminenceVProperties.png\",dpi=1000)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Highest Luminosity satellite per system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_luminosity_sats = saga_sats.groupby('HOSTID', as_index = False)['Mr'].min()\n",
    "highest_luminosity_prominence = pd.merge(highest_luminosity_sats, saga, how=\"left\", left_on=\"HOSTID\", right_on=\"HOSTID\")\n",
    "highest_luminosity_prominence_SAGA = highest_luminosity_prominence[~highest_luminosity_prominence['prominence'].isna()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "spearmanr_results_lum = sp.stats.spearmanr(highest_luminosity_prominence_SAGA['prominence'], highest_luminosity_prominence_SAGA['Mr'], axis=0, nan_policy='propagate', alternative='two-sided')\n",
    "spearman_string_lum = \"Coefficient: \" + str(spearmanr_results_lum.statistic.round(2)) + \" P-value: \" + str(spearmanr_results_lum.pvalue.round(2))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "sns.scatterplot(data=highest_luminosity_prominence_SAGA, x='Mr', y='prominence', ax=ax)\n",
    "ax.set_yscale('log')\n",
    "plt.text(0.99, 0.99, spearman_string_lum, ha='right', va='top', size=15, transform = ax.transAxes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "host_property='num_sats'\n",
    "\n",
    "fig, axs = plt.subplots()\n",
    "spearmanr_results = sp.stats.spearmanr(saga['prominence'], saga[host_property], axis=0, nan_policy='propagate', alternative='less')\n",
    "display(spearmanr_results)\n",
    "\n",
    "statistic = spearmanr_results.statistic.round(2)\n",
    "pvalue = spearmanr_results.pvalue.round(2)\n",
    "spearman_string = \"Coefficient: \" + str(statistic) + \" P-value: \" + str(pvalue)\n",
    "\n",
    "sns.scatterplot(data=saga, x=host_property, y=\"prominence\", ax=axs)\n",
    "axs.set_yscale('log') # just for radius\n",
    "plt.text(2, 2, spearman_string, ha='right', va='top', size=15, transform = ax.transAxes)\n",
    "axs.set_xlabel('# satellites')\n",
    "axs.set_ylabel('$\\mathcal{P}_e$')\n",
    "axs.set_xticks(np.arange(2, 14, 1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "host_property='phi'\n",
    "saga[host_property]\n",
    "\n",
    "fig, axs = plt.subplots()\n",
    "spearmanr_results = sp.stats.spearmanr(saga['prominence'], saga[host_property], axis=0, nan_policy='propagate', alternative='less')\n",
    "display(spearmanr_results)\n",
    "\n",
    "statistic = spearmanr_results.statistic.round(2)\n",
    "pvalue = spearmanr_results.pvalue.round(2)\n",
    "spearman_string = \"Coefficient: \" + str(statistic) + \" P-value: \" + str(pvalue)\n",
    "\n",
    "sns.scatterplot(data=saga, x=host_property, y=\"prominence\", ax=axs)\n",
    "# axs.set_xscale('log') # just for radius\n",
    "axs.set_yscale('log') # just for radius\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "host_property='MORPHTYPE'\n",
    "\n",
    "fig, axs = plt.subplots()\n",
    "spearmanr_results = sp.stats.spearmanr(saga['prominence'], saga[host_property], axis=0, nan_policy='propagate', alternative='less')\n",
    "display(spearmanr_results)\n",
    "\n",
    "statistic = spearmanr_results.statistic.round(2)\n",
    "pvalue = spearmanr_results.pvalue.round(2)\n",
    "spearman_string = \"Coefficient: \" + str(statistic) + \" P-value: \" + str(pvalue)\n",
    "\n",
    "sns.scatterplot(data=saga, x=host_property, y=\"prominence\", ax=axs)\n",
    "# axs.set_xscale('log') # just for radius\n",
    "axs.set_yscale('log') # just for radius\n",
    "\n",
    "fig.set_size_inches([10, 7])\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "satellite-ellipticities",
   "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.9.16"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "593019684f7021e8fe6a4c3db3cb5af2bec7f7dac2e37296bac94ddb3e8f915d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}