F17-18-Ethan-LF_Panels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports and Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numba import jit\n",
    "\n",
    "import SAGA\n",
    "from SAGA.database import FitsTable\n",
    "from easyquery import Query\n",
    "from easyquery import QueryMaker\n",
    "from SAGA import ObjectCuts as C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "%config InlineBackend.rc = {'figure.facecolor': 'w'}\n",
    "\n",
    "plt.rcParams['font.size'] = 14\n",
    "plt.rcParams['figure.dpi'] = 200\n",
    "plt.rcParams['font.family'] = 'STIXGeneral'\n",
    "plt.rcParams['mathtext.fontset'] = 'stix'\n",
    "plt.rcParams['axes.labelsize'] = 'large'\n",
    "plt.rcParams['legend.frameon'] = False\n",
    "plt.rcParams['lines.markersize'] = 4\n",
    "plt.rcParams['xtick.minor.visible'] = True\n",
    "plt.rcParams['ytick.minor.visible'] = True\n",
    "plt.rcParams['xtick.major.size'] = 6\n",
    "plt.rcParams['xtick.minor.size'] = 4\n",
    "plt.rcParams['ytick.major.size'] = 6\n",
    "plt.rcParams['ytick.minor.size'] = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mr_bins_fine = np.linspace(-22,12,5000)\n",
    "Mr_bins = np.linspace(-22,-11,10)\n",
    "Mr_limit = -12.3\n",
    "\n",
    "r_bins = np.linspace(10,300,11)\n",
    "r_bins_fine = np.linspace(10,300,5000)\n",
    "rproj_limit = 300."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def N_gr_Mr_new(Mr,radii,prob,bins,radius_max,Mr_min):\n",
    "    \"\"\"\n",
    "    Returns a luminosity function of surviving subhalos, subject to maxmimum\n",
    "    distance and brightest magnitude.\n",
    "    \n",
    "    Args:\n",
    "        Mr (array of floats): Magnitude of each subhalo\n",
    "        radii (array of floats): Galacto-centric distance of each subhalo\n",
    "        prob (array of floats): Survival probability of each subhalo\n",
    "        bins (array of floats): Magnitude bin edges\n",
    "        radius_max (float): Maximum radius within which to consider subhalos\n",
    "        Mr_min (float): Brightest magnitude to which to consider subhalos\n",
    "    \n",
    "    Returns:\n",
    "        Cumulative number of surviving subhalos less than each magnitude\n",
    "    \"\"\"\n",
    "    bin_edges = np.hstack([-np.inf, bins])\n",
    "    idx = (radii < radius_max) & (Mr > Mr_min) & np.isfinite(Mr)\n",
    "    n = np.histogram(Mr[idx], bins=bin_edges, weights=prob[idx])[0]\n",
    "    return np.cumsum(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def N_less_r(SubsMr,SubsRadii,prob,binsr,Mr_min):\n",
    "    \"\"\"\n",
    "    Calculates the cumulative distribution function of Galactocentric radius, subject\n",
    "    to some cuts on magnitude.\n",
    "    \n",
    "    Inputs:\n",
    "        SubsMr (array): magnitudes of subhalos\n",
    "        SubsRadii (array): Galactocentric radii of subhalos\n",
    "        prob (array): distruption probabilities of subhalos\n",
    "        binsr (array): Galactocentric radius bin edges\n",
    "        Mr_min (float): Minimum subhalo magnitude to accept\n",
    "    \n",
    "    Returns:\n",
    "        Array with length `binsr.size-1`, with number of subhalos at radii less than\n",
    "        distances `binsr[1:]`.\n",
    "    \"\"\"\n",
    "    idx = (SubsMr > Mr_min) & np.isfinite(SubsMr)\n",
    "    SubsRadii = SubsRadii[idx]\n",
    "    prob = prob[idx]\n",
    "    n_in_bin = np.histogram(SubsRadii, bins=binsr, weights=prob)[0]\n",
    "    return np.cumsum(n_in_bin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def Mr_from_r(r,host_dist):\n",
    "    return r - 5.*np.log10(host_dist) - 25."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load SAGA Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga = SAGA.QuickStart()\n",
    "saga.set_default_base_version(\"paper2\")\n",
    "hosts = saga.host_catalog.load(include_stats=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "completed_hosts = saga.host_catalog.load(query=\"paper2_complete\", include_stats=\"remote\")\n",
    "observed_hosts = saga.host_catalog.load(query=\"paper2_observed\", include_stats=\"remote\")\n",
    "bright_completed_hosts = saga.host_catalog.load(query=\"paper2_bright_complete\", include_stats=\"remote\")\n",
    "\n",
    "sats_url = \"https://www.dropbox.com/sh/7qeuqkq0c591k2g/AAC0c5C7erNCp0nYOhcDv9kva/saga_sats_latest.fits?dl=1\"\n",
    "sats = FitsTable(sats_url).read()\n",
    "complete_host_query = saga.host_catalog.construct_host_query(\"paper2_complete\")\n",
    "sats_completed_hosts = complete_host_query.filter(sats)\n",
    "\n",
    "url = \"https://www.dropbox.com/sh/7qeuqkq0c591k2g/AADUIE4j1t2vXqKpaHoCY0L7a/saga_observed_latest.fits?dl=1\"\n",
    "specs_and_potential_sats = FitsTable(url).read()\n",
    "\n",
    "specs = C.has_spec.filter(specs_and_potential_sats)\n",
    "possible_sats = (~C.has_spec).filter(specs_and_potential_sats)\n",
    "\n",
    "sats_this_host = {}\n",
    "possible_sats_this_host = {}\n",
    "\n",
    "hosts = set(specs_and_potential_sats[\"HOSTID\"])\n",
    "for host in hosts:\n",
    "    sats_this_host[host] = (QueryMaker.equals(\"HOSTID\", host) & C.is_sat).filter(specs_and_potential_sats)\n",
    "    possible_sats_this_host[host] = (QueryMaker.equals(\"HOSTID\", host)&(~C.has_spec)).filter(specs_and_potential_sats)\n",
    "\n",
    "saga_colors = np.load('EON_data/color_dict.npz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "name_dict = {}\n",
    "for key in saga_colors.keys():\n",
    "    name_dict[key] = completed_hosts[completed_hosts['HOSTID']==key]['COMMON_NAME'][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga_lfs_uncorrected = {}\n",
    "saga_rs_uncorrected = {}\n",
    "\n",
    "for host in completed_hosts['HOSTID']:\n",
    "    sats = sats_completed_hosts[sats_completed_hosts['HOSTID']==host]\n",
    "    saga_lfs_uncorrected[host] = N_gr_Mr_new(sats['Mr'],sats['RHOST_KPC'],np.ones(len(sats['Mr'])),Mr_bins_fine,\n",
    "                                             rproj_limit,-25.)\n",
    "    saga_rs_uncorrected[host] = (N_less_r(sats['Mr'],sats['RHOST_KPC'],np.ones(len(sats['Mr'])),r_bins_fine,-25.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sats_num_draws = {}\n",
    "num_draws = 1000\n",
    "        \n",
    "saga_lfs = {}\n",
    "\n",
    "for host in completed_hosts['HOSTID']:\n",
    "    saga_lfs[host] = []\n",
    "    sats = sats_completed_hosts[sats_completed_hosts['HOSTID']==host]\n",
    "    Mr_possible = Mr_from_r(possible_sats_this_host[host]['r_mag'],possible_sats_this_host[host]['HOST_DIST'])\n",
    "    psat_possible = possible_sats_this_host[host]['p_sat_corrected']\n",
    "    for j in range(num_draws):\n",
    "        rands = np.random.rand(len(psat_possible))\n",
    "        Mrs = np.concatenate((sats['Mr'],Mr_possible[rands<psat_possible]))\n",
    "        rs = np.concatenate((sats['RHOST_KPC'],50.*np.ones(len(Mr_possible[rands<psat_possible]))))\n",
    "        saga_lfs[host].append(N_gr_Mr_new(Mrs,rs,np.ones(len(Mrs)),Mr_bins_fine,rproj_limit,-25.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Simulation Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('EON_data/LF_all_pred.pickle', 'rb') as handle:\n",
    "    LF_all_pred = pickle.load(handle)\n",
    "    \n",
    "with open('EON_data/LF_all_pred_std.pickle', 'rb') as handle:\n",
    "    LF_all_pred_std = pickle.load(handle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,8))\n",
    "position = range(1,len(np.unique(completed_hosts['HOSTID']))+1)\n",
    "\n",
    "k=0\n",
    "for saga_id in np.array(list(saga_colors.keys()))[::-1]:\n",
    "    plt.subplot(6,6,position[k])\n",
    "    plt.fill_between(Mr_bins,(LF_all_pred[saga_id]-LF_all_pred_std[saga_id]).clip(min=1e-2),\n",
    "                         LF_all_pred[saga_id]+LF_all_pred_std[saga_id],\n",
    "                         facecolor=\"#003660\",alpha=0.3,rasterized=True)\n",
    "    plt.fill_between(Mr_bins,(LF_all_pred[saga_id]-2.*LF_all_pred_std[saga_id]).clip(min=1e-5),\n",
    "                         LF_all_pred[saga_id]+2.*LF_all_pred_std[saga_id],\n",
    "                         facecolor=\"#04859B\",alpha=0.15,rasterized=True)\n",
    "    plt.plot(Mr_bins,LF_all_pred[saga_id],lw=2,c=\"#003660\")\n",
    "    plt.plot(Mr_bins_fine,saga_lfs_uncorrected[saga_id],c=saga_colors[saga_id],lw=2.5,label='SAGA (Complete)')\n",
    "    plt.fill_between(Mr_bins_fine,saga_lfs_uncorrected[saga_id],np.median(np.array(saga_lfs[saga_id]),axis=0),\n",
    "                     facecolor=saga_colors[saga_id],alpha=0.45,rasterized=True)\n",
    "    plt.plot(np.ones(10)*Mr_limit,np.linspace(0.1,40,10),'k--',lw=0.5)\n",
    "    plt.xlim(-21,-11)\n",
    "    plt.ylim(0.5,20)\n",
    "    plt.yscale('log')\n",
    "    plt.text(-20.5,10,name_dict[saga_id],fontsize=11)\n",
    "    plt.tick_params(axis='both', which='major', labelsize=11)\n",
    "    if k in np.array([30,31,32,33,34,35]):\n",
    "        plt.xlabel(r'$M_r$',fontsize=14)\n",
    "        plt.xticks([-20,-18,-16,-14,-12], fontsize=10.5)\n",
    "    else:\n",
    "        plt.xticks([-20,-18,-16,-14,-12],['','','','',''])\n",
    "        plt.gca().set_xticklabels([])\n",
    "    if k in np.array([0,6,12,18,24,30]):\n",
    "        plt.ylabel(r'$N(<M_r)$',fontsize=14)\n",
    "        #plt.yticks([1,10],[1,10],fontsize=11)        \n",
    "    else:\n",
    "        plt.gca().set_yticklabels([])\n",
    "    k+=1\n",
    "\n",
    "plt.tight_layout(h_pad=0.5, w_pad=0.5)\n",
    "plt.savefig('saga_lfs.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('EON_data/r_normed_pred.pickle', 'rb') as handle:\n",
    "    r_normed_pred = pickle.load(handle)\n",
    "    \n",
    "with open('EON_data/r_normed_pred_std.pickle', 'rb') as handle:\n",
    "    r_normed_pred_std = pickle.load(handle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,8))\n",
    "position = range(1,len(np.unique(completed_hosts['HOSTID']))+1)\n",
    "\n",
    "k=0\n",
    "for saga_id in np.array(list(saga_colors.keys()))[::-1]:\n",
    "    plt.subplot(6,6,position[k])\n",
    "    plt.ylim(0,1.2)\n",
    "    plt.fill_between(np.concatenate((np.array([10.]),r_bins[1:])),\n",
    "                     np.concatenate((np.array([0.]),(r_normed_pred[saga_id]-r_normed_pred_std[saga_id]))),\n",
    "                     (np.concatenate((np.array([0.]),(r_normed_pred[saga_id]+r_normed_pred_std[saga_id])))).clip(max=1.),\n",
    "                     facecolor=\"#003660\",alpha=0.3)\n",
    "    plt.fill_between(np.concatenate((np.array([10.]),r_bins[1:])),\n",
    "                     np.concatenate((np.array([0.]),(r_normed_pred[saga_id]-2.*r_normed_pred_std[saga_id]))),\n",
    "                     (np.concatenate((np.array([0.]),(r_normed_pred[saga_id]+2.*r_normed_pred_std[saga_id])))).clip(max=1.),\n",
    "                     facecolor=\"#04859B\",alpha=0.15)\n",
    "    plt.plot(np.concatenate((np.array([10.]),r_bins[1:])),np.concatenate((np.array([0.]),(r_normed_pred[saga_id]))),lw=2,c=\"#003660\")\n",
    "    idx = np.argmin(np.abs(Mr_bins_fine-Mr_limit))\n",
    "    if 1.*saga_lfs_uncorrected[saga_id][idx] > 0:\n",
    "        plt.plot(r_bins_fine[1:],\n",
    "                 (1.*saga_rs_uncorrected[saga_id]/(1.*saga_lfs_uncorrected[saga_id][idx])).clip(max=1.),\n",
    "                 c=saga_colors[saga_id],lw=2.5,label='SAGA (Complete)')\n",
    "    plt.xlim(0,300)\n",
    "    plt.axhline(1, c=\"k\", lw=0.8, ls=\":\", zorder=0)\n",
    "    plt.text(15.,1.03,name_dict[saga_id],fontsize=11)\n",
    "    if k in np.array([30,31,32,33,34,35]):\n",
    "        plt.xlabel(r'$r_{\\rm{proj}}\\ [\\rm{kpc}]$',fontsize=14)\n",
    "        plt.xticks([0,100,200,300],[0,100,200,300],fontsize=11)\n",
    "    else:\n",
    "        plt.xticks([0,100,200,300],['','','',''])\n",
    "    if k in np.array([0,6,12,18,24,30]):\n",
    "        plt.ylabel(r'$\\tilde{N}(<r_{\\rm{proj}})$',fontsize=14)\n",
    "        plt.yticks([0,1],[0,1],fontsize=11)\n",
    "    else:\n",
    "        plt.yticks([0,1],['',''])\n",
    "    k+=1\n",
    "    \n",
    "plt.tight_layout(h_pad=0.5, w_pad=0.5)\n",
    "plt.savefig('saga_radial.pdf')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}