VSMDPL-halo-ellipticities-include-interlopers.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import os\n",
    "from easyquery import Query\n",
    "from fast3tree import fast3tree\n",
    "from dotenv import load_dotenv\n",
    "\n",
    "load_dotenv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VSMDPL_DATA = os.getenv('VSMDPL_DATA') # path to parquet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hosts_and_subs = pd.read_parquet(VSMDPL_DATA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hosts, subs = Query(\"upid == -1\").split(hosts_and_subs)\n",
    "\n",
    "subs = hosts_and_subs.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "box_len = 160"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grab Interlopers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# threshold redshift\n",
    "dz_thres = 2.75 #coordinate z + vz / 100 , this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hosts[\"z_vz\"] = np.remainder(hosts[\"z\"] + hosts[\"vz\"] / 100.0, box_len)\n",
    "subs[\"z_vz\"] = np.remainder(subs[\"z\"] + subs[\"vz\"] / 100.0, box_len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interloper_groups = pd.DataFrame()\n",
    "interloper_groups\n",
    "\n",
    "with fast3tree(np.vstack([subs[\"x\"].values, subs[\"y\"].values]).T) as tree:\n",
    "    tree.set_boundaries(0, box_len)\n",
    "    for host in hosts.itertuples(): \n",
    "        center = (host.x, host.y)\n",
    "        idx = tree.query_radius(center, host.rvir * 0.001, periodic=True)\n",
    "        dz = np.remainder(subs[\"z_vz\"].values[idx] - host.z_vz, box_len)\n",
    "        idx = idx[(dz < dz_thres) | (dz > box_len - dz_thres)]\n",
    "        sys_subs = subs.iloc[idx].copy()\n",
    "        sys_subs[\"interloper_host_id\"] = host.id\n",
    "        sys_subs[\"interloper_host_x\"] = host.x\n",
    "        sys_subs[\"interloper_host_y\"] = host.y\n",
    "        sys_subs[\"interloper_host_vmax_mpeak\"] = host.vmax_mpeak\n",
    "        interloper_groups = interloper_groups.append(sys_subs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interloper_groups['x_adj'] = interloper_groups['x']\n",
    "interloper_groups['x_adj'].mask(interloper_groups['x'] - interloper_groups['interloper_host_x'] > box_len/2, interloper_groups['x'] - 160, inplace=True)\n",
    "interloper_groups['x_adj'].mask(interloper_groups['x'] - interloper_groups['interloper_host_x'] < -box_len/2, interloper_groups['x'] + 160, inplace=True)\n",
    "\n",
    "interloper_groups['y_adj'] = interloper_groups['y']\n",
    "interloper_groups['y_adj'].mask(interloper_groups['y'] - interloper_groups['interloper_host_y'] > box_len/2, interloper_groups['y'] - 160, inplace=True)\n",
    "interloper_groups['y_adj'].mask(interloper_groups['y'] - interloper_groups['interloper_host_y'] < -box_len/2, interloper_groups['y'] + 160, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VSMDPL_subs_interlopers = interloper_groups.copy()\n",
    "%store VSMDPL_subs_interlopers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "groups_cut = interloper_groups[interloper_groups['vmax_mpeak'] >= 46.77].copy()\n",
    "\n",
    "groups_cut = groups_cut.groupby(['interloper_host_id']).filter(lambda x: len(x) > 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "average_num_subs = len(groups_cut)/len(groups_cut['interloper_host_id'].unique())\n",
    "print(\"Average # subhalos per system: {}\".format(average_num_subs))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Ellipticity calculations</h2>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Quadrupole moments</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ellipticity_df = groups_cut.copy()\n",
    "ellipticity_df['x_diff_sq'] = (ellipticity_df['x'] - ellipticity_df['interloper_host_x'])**2\n",
    "ellipticity_df['y_diff_sq'] = (ellipticity_df['y'] - ellipticity_df['interloper_host_y'])**2\n",
    "ellipticity_df['x_diff_y_diff'] = (ellipticity_df['y'] - ellipticity_df['interloper_host_y'])*(ellipticity_df['x'] - ellipticity_df['interloper_host_x'])\n",
    "ellipticity_df = ellipticity_df.groupby('interloper_host_id', as_index = False)[['x_diff_sq', 'y_diff_sq', 'x_diff_y_diff']].mean().rename(columns={'x_diff_sq':'Q_xx', 'y_diff_sq':'Q_yy', 'x_diff_y_diff':'Q_xy'}).copy()\n",
    "ellipticity_df"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<h4>Ellipticity components</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def e_1 (Q_xx, Q_yy, Q_xy):\n",
    "    return (( Q_xx - Q_yy )/(Q_xx + Q_yy + 2*(Q_xx*Q_yy-Q_xy**2)**(1/2)))\n",
    "def e_2 (Q_xx, Q_yy, Q_xy):\n",
    "    return (2*Q_xy/(Q_xx + Q_yy + 2*(Q_xx*Q_yy-Q_xy**2)**(1/2)))\n",
    "def e_3 (e_1, e_2):\n",
    "    return ( e_1**2 + e_2**2 )**(1/2)\n",
    "\n",
    "\n",
    "ellipticity_df['e_1'] = e_1(ellipticity_df['Q_xx'], ellipticity_df['Q_yy'], ellipticity_df['Q_xy'])\n",
    "ellipticity_df['e_2'] = e_2(ellipticity_df['Q_xx'], ellipticity_df['Q_yy'], ellipticity_df['Q_xy'])\n",
    "ellipticity_df['e_3'] = e_3(ellipticity_df['e_1'], ellipticity_df['e_2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "custom_style = {\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",
    "}\n",
    "\n",
    "sns.set_style(\"darkgrid\", custom_style)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.histplot(ellipticity_df[\"e_3\"], bins=7, element='poly', stat=\"density\", fill=False)\n",
    "ax.set(xlabel=\"$\\longleftarrow$ more elliptical         (b/a)           more spherical $\\longrightarrow$\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Store ellipticity DataFrame to use in other nbs</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ellipticity_df_VSMDPL_interlopers = ellipticity_df.copy()\n",
    "%store ellipticity_df_VSMDPL_interlopers"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Highest luminosity analysis</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add num satellites\n",
    "groups_cut['num_satellites'] = groups_cut.groupby('interloper_host_id')['interloper_host_id'].transform('count')\n",
    "highest_luminosity_subs = groups_cut.groupby('interloper_host_id', as_index = False)[['Mr', 'num_satellites']].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_luminosity_ellipticity = pd.merge(highest_luminosity_subs, ellipticity_df, how=\"left\", left_on=\"interloper_host_id\", right_on=\"interloper_host_id\")[['interloper_host_id', 'Mr', 'e_3']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_luminosity_ellipticity_VSMDPL_interlopers = highest_luminosity_ellipticity[~highest_luminosity_ellipticity['e_3'].isna()]\n",
    "%store highest_luminosity_ellipticity_VSMDPL_interlopers\n",
    "highest_luminosity_ellipticity_VSMDPL_interlopers"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Concat highest luminosity w/ log(vmax_mpeak)</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_vmax_subs = groups_cut.groupby('interloper_host_id', as_index = False)['vmax_mpeak'].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_vmax_ellipticity = pd.merge(highest_vmax_subs, ellipticity_df, how=\"left\", left_on=\"interloper_host_id\", right_on=\"interloper_host_id\")[['interloper_host_id', 'vmax_mpeak', 'e_3']]\n",
    "highest_vmax_ellipticity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_vmax_ellipticity['log_vmax_mpeak'] = np.log(highest_vmax_ellipticity['vmax_mpeak'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_vmax_ellipticity_VSMDPL_interlopers = highest_vmax_ellipticity[~highest_vmax_ellipticity['e_3'].isna()]\n",
    "%store highest_vmax_ellipticity_VSMDPL_interlopers\n",
    "highest_vmax_ellipticity_VSMDPL_interlopers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_xlabel('\"$\\longleftarrow$ more elliptical         ellipticity           more spherical $\\longrightarrow$\"')\n",
    "ax1.set_ylabel('log(vmax_mpeak)')\n",
    "sp1 = sns.scatterplot(data=highest_vmax_ellipticity_VSMDPL_interlopers, x='e_3', y='log_vmax_mpeak', ax=ax1, s=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sp1.figure.get_figure()\n",
    "fig.set_size_inches([6, 4])\n",
    "fig.savefig(\"figures/SimVmax_mpeakVsEllipticity.png\",dpi=600)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Average # satellites per vmax_mpeak</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hosts_num_vmax = groups_cut.groupby('interloper_host_id').first()[['interloper_host_vmax_mpeak', 'num_satellites']]\n",
    "hosts_num_vmax"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bin hosts by vmax_mpeak"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hosts_num_vmax_averages = hosts_num_vmax.groupby(pd.cut(hosts_num_vmax['interloper_host_vmax_mpeak'], np.linspace(135, 340, 20))).mean()\n",
    "hosts_num_vmax_averages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2, ax2 = plt.subplots()\n",
    "ax2.set_xlabel('vmax_mpeak')\n",
    "ax2.set_ylabel('average # satellites')\n",
    "sp2 = plt.plot(hosts_num_vmax_averages['interloper_host_vmax_mpeak'], hosts_num_vmax_averages['num_satellites'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig3, ax3 = plt.subplots()\n",
    "ax3.set_xlabel('vmax_mpeak')\n",
    "ax3.set_ylabel('# satellites')\n",
    "sp3 = plt.scatter(hosts_num_vmax['interloper_host_vmax_mpeak'], hosts_num_vmax['num_satellites'], s=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "spearmanr_results = sp.stats.spearmanr(hosts_num_vmax['interloper_host_vmax_mpeak'], hosts_num_vmax['num_satellites'], axis=0, nan_policy='propagate', alternative='two-sided')\n",
    "spearmanr_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig3 = sp3.figure.get_figure()\n",
    "fig3.set_size_inches([6, 4])\n",
    "fig3.savefig(\"figures/NumsubsVsVmax_mpeak_scatter.png\",dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}