{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Galaxy Clusters & Velocity Appendix: Peaks in the Halo mass distribution\n",
    "\n",
    "Author: Julien Peloton [@JulienPeloton](https://github.com/JulienPeloton)  \n",
    "Last Verifed to Run: 2019-01-04  \n",
    "Estimated running time: < 5 min.\n",
    "\n",
    "This notebook, with the help of Apache Spark, inspects the peaks found in the halo mass distribution in the companion notebook.\n",
    "\n",
    "**Summary:**\n",
    "\n",
    "If we look at the distribution of halo masses in cosmoDC2, all seem OK. But if we look at the same distribution after filtering objects according to the stellar mass, we start to see peaks in the distribution. These peaks are regularly spaced in log, and their position is independent on the stellar mass cut applied.\n",
    "\n",
    "\n",
    "**Useful links:**\n",
    "\n",
    "- Main issue: https://github.com/LSSTDESC/DC2-analysis/issues/57\n",
    "- Potential update of the cosmoDC2 simulation: https://github.com/LSSTDESC/cosmodc2/issues/84\n",
    "\n",
    "**Logistics:** \n",
    "\n",
    "This notebook is intended to be run through the JupyterHub NERSC interface with the desc-pyspark kernel. The kernel is automatically installed in your environment when you use the kernel setup script:\n",
    "\n",
    "```bash\n",
    "source /global/common/software/lsst/common/miniconda/kernels/setup.sh\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to the data\n",
    "fn = \"/global/cscratch1/sd/peloton/cosmodc2/xyz_v1.1.4_mass_and_mag.parquet\"\n",
    "\n",
    "# Load the data - this a lazy operation, no data movement yet!\n",
    "df = spark.read.format(\"parquet\").load(fn)\n",
    "\n",
    "# Let's inspect the schema\n",
    "df.printSchema()\n",
    "\n",
    "# Number of objects in the catalog\n",
    "print(\"Number of rows: {}\".format(df.count()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation: selecting halos by their stellar mass\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We reject synthetic halos - see for example discussion in \n",
    "# https://github.com/LSSTDESC/cosmodc2/issues/82\n",
    "# In addition we select halo members according to their stellar mass\n",
    "stellar_masses_cut = [1e9, 1e10, 5e10, 1e11]\n",
    "full_data = []\n",
    "for stellar_mass_cut in stellar_masses_cut:\n",
    "    df_filt = df.filter(\"halo_id > 0\").filter(\"stellar_mass > {}\".format(stellar_mass_cut))\n",
    "\n",
    "    # Group data by halos and compute the number of halo members\n",
    "    df_disp = df_filt.groupBy(\"halo_id\").count()\n",
    "\n",
    "    # Add back the original DataFrame columns\n",
    "    # and select only central member for halo \n",
    "    # (unique halo mass for a halo)\n",
    "    data_joined = df_filt.join(df_disp, \"halo_id\")\\\n",
    "        .filter(\"is_central == True\")\\\n",
    "        .select(\"halo_mass\", 'count')\\\n",
    "        .dropna()\n",
    "\n",
    "    # Collect the data from the executors to the driver\n",
    "    data = data_joined.collect()\n",
    "    full_data.append(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Halo mass distribution\n",
    "\n",
    "(for the version without cut, see the main notebook)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 17})\n",
    "\n",
    "fig = plt.figure(figsize=(10, 7))\n",
    "plt.title(\"Histogram of halo masses ($n_{{gal}}/halo > 0$)\")\n",
    "\n",
    "# Plot peak locations (empirical!)\n",
    "for k in range(90):\n",
    "    plt.axvline(k*0.15, ls='--', color='k', alpha=0.5)\n",
    "\n",
    "# Plot halo mass data\n",
    "for index, stellar_mass_cut in enumerate(stellar_masses_cut):\n",
    "    mass, count = np.transpose(full_data[index])\n",
    "    print(r\"Number of entries ($M_*$ > {:.1e} $M_\\odot$): {}\".format(stellar_mass_cut, len(mass)))\n",
    "    plt.hist(np.log10(mass), bins=200, alpha=0.5, label='$M_*$ > {:.1e} $M_\\odot$'.format(stellar_mass_cut))\n",
    "    \n",
    "\n",
    "plt.xlim(10.3, 15)\n",
    "plt.xlabel(r'$\\log({\\rm M}_h \\, [{\\rm M}_\\odot])$')\n",
    "plt.yscale('log')\n",
    "plt.legend(title=\"Only central galaxies with:\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The position $p$ of the peaks is independent of the cut on the stellar mass, and it seems to follow:\n",
    "\n",
    "$$ \\log(M_{h}^{p}) = 0.15 * p $$\n",
    "\n",
    "\n",
    "Let's plot the same distribution, but selecting only halos with at least 2 members"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 17})\n",
    "\n",
    "mincount = 1\n",
    "\n",
    "fig = plt.figure(figsize=(10, 7))\n",
    "plt.title(\"Histogram of halo masses ($n_{{gal}}/halo >$  {})\".format(mincount))\n",
    "\n",
    "# Plot peak locations (empirical!)\n",
    "for k in range(90):\n",
    "    plt.axvline(k*0.15, ls='--', color='k', alpha=0.5)\n",
    "\n",
    "# Plot halo mass data\n",
    "for index, stellar_mass_cut in enumerate(stellar_masses_cut):\n",
    "    mass, count = np.transpose(full_data[index])\n",
    "    mask = count > mincount\n",
    "    print(r\"Number of entries ($M_*$ > {:.1e} $M_\\odot$): {}\".format(stellar_mass_cut, len(mass[mask])))\n",
    "    plt.hist(np.log10(mass[mask]), bins=200, alpha=0.5, label='$M_*$ > {:.1e} $M_\\odot$'.format(stellar_mass_cut))\n",
    "    \n",
    "\n",
    "plt.xlim(10.3, 15)\n",
    "plt.xlabel(r'$\\log({\\rm M}_h \\, [{\\rm M}_\\odot])$')\n",
    "plt.yscale('log')\n",
    "plt.legend(title=\"Only central galaxies with:\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The peaks are even sharper at low mass."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** Is that expected?\n",
    "\n",
    "**Answer:** The saw tooth shape is indeed expected (but probably not wanted!) from the way N-body simulations are done. A detailed discussion can be found at https://github.com/LSSTDESC/DC2-analysis/issues/57, and an action item for the simulation has been opened at https://github.com/LSSTDESC/cosmodc2/issues/84."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moving to cuts on the Magnitude\n",
    "\n",
    "It has been suggested to not use the stellar mass for selecting halos (not a direct observable), but rather a more physically motivated quantity like magnitudes/colors. In this part, we select the apparent magnitude, lensed, in i band. This might not be the best quantity to look at, but just want to play with something else than the stellar mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's select halo masses, and magnitudes (i, g, u) for \n",
    "# the central galaxy.\n",
    "data = df.filter(\"halo_id > 0\")\\\n",
    "    .filter(\"is_central == True\")\\\n",
    "    .select([\"halo_mass\", \"mag_i\", \"mag_g\", \"mag_u\"])\\\n",
    "    .sample(fraction=0.05).collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mass, mag_i, mag_g, mag_u = np.transpose(data)\n",
    "print(\"Number of rows: {}\".format(len(mass)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at the 2D histogram halo mass - magnitude in i band:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "joint_kws = dict(gridsize=200, mincnt=1)\n",
    "g = sns.jointplot(\n",
    "    np.log10(mass), \n",
    "    mag_i, \n",
    "    height=8, space=0.5,\n",
    "    kind='hex', color='k', \n",
    "    xlim=(10.5, 12), ylim=(20, 35),\n",
    "    joint_kws=joint_kws, \n",
    "    marginal_kws=dict(bins=200, rug=True))\n",
    "\n",
    "g.set_axis_labels(\n",
    "    r'$\\log({\\rm M}_h \\, [{\\rm M}_\\odot])$', \n",
    "    'Apparent magnitude, lensed, in i band (lsst)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hum, the stripes are still there... Let's have a look at 1D halo mass distribution as a function of magnitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We reject synthetic halos - see for example discussion in \n",
    "# https://github.com/LSSTDESC/cosmodc2/issues/82\n",
    "# In addition we select halo members according to their stellar mass\n",
    "stellar_masses_cut = [25.0, 22.5, 20.0]\n",
    "full_data = []\n",
    "for stellar_mass_cut in stellar_masses_cut:\n",
    "    df_filt = df.filter(\"halo_id > 0\").filter(\"mag_i < {}\".format(stellar_mass_cut))\n",
    "\n",
    "    # Group data by halos and compute the number of halo members\n",
    "    df_disp = df_filt.groupBy(\"halo_id\").count()\n",
    "\n",
    "    # Add back the original DataFrame columns\n",
    "    # and select only central member for halo \n",
    "    # (unique halo mass for a halo)\n",
    "    data_joined = df_filt.join(df_disp, \"halo_id\")\\\n",
    "        .filter(\"is_central == True\")\\\n",
    "        .select(\"halo_mass\", 'count')\\\n",
    "        .dropna()\n",
    "\n",
    "    # Collect the data from the executors to the driver\n",
    "    data = data_joined.collect()\n",
    "    full_data.append(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 17})\n",
    "\n",
    "fig = plt.figure(figsize=(10, 7))\n",
    "plt.title(\"Histogram of halo masses ($n_{{gal}}/halo > 0$)\")\n",
    "\n",
    "# Plot peak locations (empirical!)\n",
    "for k in range(90):\n",
    "    plt.axvline(k*0.15, ls='--', color='k', alpha=0.5)\n",
    "\n",
    "# Plot halo mass data\n",
    "for index, stellar_mass_cut in enumerate(stellar_masses_cut):\n",
    "    mass, count = np.transpose(full_data[index])\n",
    "    print(r\"Number of entries (mag_i_lsst < {}): {}\".format(stellar_mass_cut, len(mass)))\n",
    "    plt.hist(np.log10(mass), bins=200, alpha=0.5, label='mag_i_lsst < {}'.format(stellar_mass_cut))\n",
    "    \n",
    "\n",
    "plt.xlim(10.3, 15)\n",
    "plt.xlabel(r'$\\log({\\rm M}_h \\, [{\\rm M}_\\odot])$')\n",
    "plt.yscale('log')\n",
    "plt.legend(title=\"Only central galaxies with:\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to what is seen for cut on the stellar mass... According to Andrew Hearin in [#57](https://github.com/LSSTDESC/DC2-analysis/issues/57), this makes sense as\n",
    "\n",
    "```\n",
    "... in the cosmoDC2 model, restframe flux derives from stellar mass, \n",
    "which in turn drives from halo mass, so it is expected that this discreteness \n",
    "propagates through to these other variables.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desc-pyspark",
   "language": "python",
   "name": "desc-pyspark"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}