F12-LF.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "import sys\n",
    "sys.path.insert(0, '/home/yymao/Dropbox/Academia/Collaborations/SAGA-local/saga')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from itertools import combinations\n",
    "from collections import Counter, defaultdict\n",
    "from astropy.table import Table, join, vstack\n",
    "from astropy.coordinates import SkyCoord, Distance\n",
    "from SAGA.utils import add_skycoord\n",
    "from easyquery import Query, QueryMaker\n",
    "from astropy.io import ascii\n",
    "\n",
    "import SAGA\n",
    "from SAGA.database import FitsTable, GoogleSheets\n",
    "from SAGA import ObjectCuts as C\n",
    "from SAGA.utils import add_skycoord, fill_values_by_query\n",
    "from SAGA.utils.distance import z2v, d2m, m2d\n",
    "from SAGA.utils.display import show_images\n",
    "from SAGA.objects.object_catalog import calc_fiducial_p_sat, calc_fiducial_p_sat_corrected\n",
    "print(SAGA.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga = SAGA.QuickStart()\n",
    "saga.set_default_base_version(\"paper2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lg = SAGA.database.FitsTable(\"../data/McConnachie2012.fit\").read()\n",
    "lg[\"Mr\"] = lg[\"vmag_lc\"] - 0.4\n",
    "lg[\"SubG\"][np.argmax(lg[\"Name\"] == \"Pegasus dIrr\")] = \"M31 \"\n",
    "lg = (Query((np.isfinite, \"vmag_lc\"), ~QueryMaker.equals(\"Name\", \"Canis Major\"), \"Mr < -10\")).filter(lg)\n",
    "\n",
    "mw_sat_mr = np.asarray(np.sort(QueryMaker.equals(\"SubG\", \"MW  \").filter(lg, \"Mr\")))\n",
    "m31_sat_mr = np.asarray(np.sort(QueryMaker.equals(\"SubG\", \"M31 \").filter(lg, \"Mr\")))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mw = SAGA.database.GoogleSheets(\"1O8tGgnHXRcAT8P78J3V2pWVb3CWE1Hvcowily79tdg4\", 0).read()\n",
    "mw_sat_mr[:4] = np.sort(mw[\"Mr_o\"])[:4]\n",
    "mw_sat_mr[4] = -12.3  # Leo I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def annotate_base(base):\n",
    "    base[\"Mr_\"] = base[\"r_mag\"] - d2m(base[\"HOST_DIST\"])\n",
    "    base[\"Mr\"] = np.where(np.isfinite(base[\"Mr\"]) & C.has_spec.mask(base), base[\"Mr\"], base[\"Mr_\"])\n",
    "    base[\"log_sm\"] = np.where(np.isfinite(base[\"log_sm\"]), base[\"log_sm\"], 1.254 + 1.0976 * base[\"gr\"] - 0.4 * base[\"Mr\"])\n",
    "    \n",
    "    base[\"human_selected\"] = 0\n",
    "    hs = saga.database[\"human_selected\"].read()\n",
    "    for i in range(1, 4):\n",
    "        base[\"human_selected\"] += np.in1d(base[\"OBJID\"], Query(f\"score == {i}\").filter(hs, \"OBJID\")).astype(np.int) * i\n",
    "\n",
    "    return base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base = SAGA.database.FitsTable(\"/home/yymao/Documents/Research/SAGA/PaperII/data-archive/saga_base_all.fits\").read()\n",
    "base = saga.host_catalog.construct_host_query(\"paper2_complete\").filter(base)\n",
    "assert len(np.unique(base[\"HOSTID\"])) == 36\n",
    "\n",
    "base = annotate_base(base)\n",
    "sats = C.is_sat.filter(base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = Query(QueryMaker.equals(\"survey\", \"NSA\"), \"SATS == 3\").filter(base, [\"Mr\", \"HOST_MK\", \"HOSTID\", \"log_sm\"])\n",
    "print(len(t))\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(10, 4))\n",
    "ax[0].scatter(t[\"HOST_MK\"], t[\"Mr\"])\n",
    "ax[0].plot([-24.6, -23], np.array([-24.6, -23])+2.5, c=\"C1\")\n",
    "ax[1].hist(t[\"log_sm\"]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hostlist = saga.host_catalog.load(include_stats=True, query=\"paper2_complete\")[\"HOSTID\", \"sats_Mr_limit\", \"sats_total\", \"COMMON_NAME\", \"K_ABS\"]\n",
    "hostlist = join(hostlist, t, \"HOSTID\", \"left\")\n",
    "hostlist[\"Mr\"] = np.where(hostlist[\"Mr\"].mask, hostlist[\"K_ABS\"] + 2.5, hostlist[\"Mr\"])\n",
    "print(\"completeness limit =\", np.median(-12.3 - (hostlist[\"Mr\"])))\n",
    "\n",
    "color_dict = np.load(\"/home/yymao/Documents/Research/SAGA/PaperII/data-archive/color_dict.npz\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lf_data = dict()\n",
    "dm_bins = np.linspace(0, 11, 15)\n",
    "\n",
    "for base_this in base.group_by(\"HOSTID\").groups:\n",
    "    host_this = base_this[\"HOSTID\"][0]\n",
    "    host_mag = hostlist[hostlist[\"HOSTID\"] == host_this][\"Mr\"][0]\n",
    "    base_this = Query(\"p_sat_corrected > 0\", C.is_sat | (~C.has_spec)).filter(base_this, [\"Mr\", \"p_sat_corrected\", \"SATS\"])\n",
    "    lf_data_this = dict()\n",
    "    base_this.sort('Mr')    \n",
    "    lf_data_this['lf_x'] = np.append(np.insert(np.vstack([base_this['Mr']]*2).T.flatten(), 0, -23.0), -9.0)\n",
    "    lf_data_this['lf_y_est'] = np.vstack([np.insert(np.cumsum(base_this['p_sat_corrected']), 0, 0.0)]*2).T.flatten()\n",
    "    lf_data_this['lf_y'] = np.vstack([np.insert(np.cumsum(C.is_sat.mask(base_this).astype(np.int)), 0, 0.0)]*2).T.flatten()\n",
    "    lf_data_this['dm'] = np.searchsorted(C.is_sat.filter(base_this, \"Mr\") - host_mag, dm_bins)\n",
    "    lf_data_this['dm_cor'] = np.insert(np.cumsum(base_this['p_sat_corrected']),0,0)[np.searchsorted(base_this['Mr'] - host_mag, dm_bins)]\n",
    "    lf_data[host_this] = lf_data_this\n",
    "    \n",
    "for host_this, mr in zip((\"MW\", \"M31\"), (mw_sat_mr, m31_sat_mr)):\n",
    "    lf_data_this = dict()\n",
    "    lf_data_this['lf_x'] = np.append(np.insert(np.vstack([mr]*2).T.flatten(), 0, -23.0), -9.0)\n",
    "    lf_data_this['lf_y'] = np.vstack([np.insert(np.arange(1, len(mr)+1), 0, 0.0)]*2).T.flatten()\n",
    "    lf_data[host_this] = lf_data_this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lf_sales = ascii.read(format=\"basic\", names=[\"dm\", \"n1\", \"n2\"], data_start=0,\n",
    "    table=\"\"\"0.14 0.00698 0.01323\n",
    "0.74 0.04821 0.10645\n",
    "1.34 0.09673 0.23615\n",
    "1.94 0.15467 0.40442\n",
    "2.54 0.23212 0.61988\n",
    "3.14 0.32869 0.87778\n",
    "3.74 0.44383 1.22344\n",
    "4.34 0.55660 1.61749\n",
    "4.94 0.67625 2.08275\"\"\")\n",
    "\n",
    "lf_lv = np.array([[ 0.5       ,  0.        ,  0.        ],\n",
    "       [ 1.5       ,  0.25      ,  0.1767767 ],\n",
    "       [ 2.5       ,  0.375     ,  0.21650635],\n",
    "       [ 3.5       ,  0.655     ,  0.28613808],\n",
    "       [ 4.5       ,  1.03      ,  0.3588175 ],\n",
    "       [ 5.5       ,  1.53      ,  0.43732139],\n",
    "       [ 6.5       ,  2.21875   ,  0.52663436],\n",
    "       [ 7.5       ,  2.865     ,  0.59843546],\n",
    "       [ 8.5       ,  3.02125   ,  0.61453743],\n",
    "       [ 9.5       ,  4.17875   ,  0.72273353],\n",
    "       [10.5       ,  5.0325    ,  0.7931346 ],\n",
    "       [11.5       ,  5.8475    ,  0.85494883]])\n",
    "\n",
    "lf_n12 = ascii.read(format=\"basic\", names=[\"dm\", \"n\", \"nl\", \"nu\"], data_start=0,\n",
    "    table=\"\"\"1.5  0.15319  0.06332 0.27467\n",
    "2.0  0.21857  0.11069 0.35585\n",
    "2.5  0.25843  0.13981 0.39589\n",
    "3.0  0.31824  0.14268 0.48504\n",
    "3.5  0.52342  0.31026 0.74681\n",
    "4.0  0.61890  0.42506 0.86090\n",
    "4.5  0.91498  0.67128 1.18541\n",
    "5.0  1.11534  0.75443 1.47465\n",
    "5.5  1.84379  1.31211 2.47517\n",
    "6.0  2.81017  1.99983 3.79168\n",
    "6.5  0.97246  0.23586 2.11470\n",
    "7.0  4.69293  2.06169 7.83712\n",
    "7.5  6.72983  3.28919 9.60195\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=2, figsize=(10, 4.5))\n",
    "\n",
    "ax_this = ax[0]\n",
    "\n",
    "for hostid, color in color_dict.items():\n",
    "    host_this = hostlist[np.argmax(hostlist[\"HOSTID\"] == hostid)]\n",
    "    lf_data_this = lf_data[hostid]\n",
    "    ax_this.plot(lf_data_this['lf_x'], lf_data_this['lf_y'], lw=1.5, c=color, label=host_this['COMMON_NAME'], alpha=0.95, rasterized=True)\n",
    "    ax_this.fill_between(lf_data_this['lf_x'], lf_data_this['lf_y'], lf_data_this['lf_y_est'], lw=0, color=color, alpha=0.25, rasterized=True)\n",
    "\n",
    "for hostid, ls in zip((\"MW\", \"M31\"), ('--', '-.')): \n",
    "    lf_data_this = lf_data[hostid]\n",
    "    ax_this.plot(lf_data_this['lf_x'], lf_data_this['lf_y'], lw=1, c=\"k\", label=hostid, ls=ls)\n",
    "    \n",
    "nrow1, nrow2 = 18, 13\n",
    "nrow3 = len(hostlist) + 2 - nrow1 - nrow2\n",
    "(lines, labels) = ax_this.get_legend_handles_labels()\n",
    "l = plt.Line2D([0], [0], linestyle='none', marker='')\n",
    "for _ in range(nrow1-nrow2):\n",
    "    lines.insert(nrow1 + nrow2, l)\n",
    "    labels.insert(nrow1 + nrow2, \"\")\n",
    "for _ in range(nrow1-nrow3):\n",
    "    lines.append(l)\n",
    "    labels.append(\"\")\n",
    "\n",
    "ax_this.legend(lines, labels, fontsize=10, ncol=3, loc='upper left', handlelength=0.9, handletextpad=0.5, columnspacing=0.6, borderpad=0.1, labelspacing=0.25)\n",
    "ax_this.axvspan(-12.3, -10, color='k', alpha=0.15, rasterized=True, zorder=99)\n",
    "ax_this.set_xlim(-22, -10)\n",
    "ax_this.set_ylim(-0.1, 13)\n",
    "ax_this.set_xlabel(r'$M_{r, o}$  [mag]')\n",
    "ax_this.set_ylabel(r'$N_{\\rm sat}(<M_{r, o})$')\n",
    "ax_this.yaxis.set_ticks_position('both')\n",
    "\n",
    "ax_this = ax[1]\n",
    "\n",
    "y = np.mean(np.vstack([lf_data_this['dm'] for lf_data_this in lf_data.values() if \"dm\" in lf_data_this]), axis=0)\n",
    "y_cor = np.mean(np.vstack([lf_data_this['dm_cor'] for lf_data_this in lf_data.values() if \"dm_cor\" in lf_data_this]), axis=0)\n",
    "nhost = 36\n",
    "\n",
    "ax_this.errorbar(dm_bins, y, yerr=(np.sqrt(y*nhost)/nhost, np.sqrt((y+y_cor)*nhost)/nhost + y_cor-y), ls=\"-\", marker=\"s\", ms=5, c=\"C2\", zorder=10, label=\"SAGA (This work)\", lw=1)\n",
    "\n",
    "ax_this.errorbar(lf_n12[\"dm\"], lf_n12[\"n\"], yerr=(lf_n12[\"n\"]-lf_n12[\"nl\"], lf_n12[\"nu\"]-lf_n12[\"n\"]),  color=\"C8\", ls=\"-.\", lw=1, label=\"Nierenberg+12\\n(COSMOS)\")\n",
    "\n",
    "ax_this.plot(lf_sales[\"dm\"], lf_sales[\"n1\"], color=\"C4\", lw=1, label=\"Sales+13 (SDSS)\")\n",
    "ax_this.plot(lf_sales[\"dm\"], lf_sales[\"n2\"], color=\"C4\", lw=1)\n",
    "ax_this.fill_between(lf_sales[\"dm\"], lf_sales[\"n1\"], lf_sales[\"n2\"], color=\"C4\", alpha=0.3)\n",
    "\n",
    "ax_this.errorbar(lf_lv[1:,0], lf_lv[1:,1]*2, yerr=lf_lv[1:,2]*2, color=\"C0\", ls=\"--\", lw=1, label=\"Carlsten+20 (LV)\")\n",
    "\n",
    "(lines, labels) = ax_this.get_legend_handles_labels()\n",
    "lines = [lines[i] for i in [1,2,0,3]]\n",
    "labels = [labels[i] for i in [1,2,0,3]]\n",
    "\n",
    "ax_this.legend(lines, labels, loc=\"lower right\", fontsize=13)\n",
    "ax_this.set_yscale(\"log\")\n",
    "ax_this.set_xlim(1, 9.5)\n",
    "ax_this.set_xlabel(r'$\\Delta m \\equiv m^\\mathrm{(sat)} - m^\\mathrm{(host)} }$  [mag]')\n",
    "ax_this.set_ylabel(r'$\\langle N_{\\rm sat}(<\\Delta m) \\rangle$')\n",
    "ax_this.set_ylim(0.01, 10)\n",
    "ax_this.set_xticks(list(range(1, 10)))\n",
    "ax_this.axvspan(8.8, 11, color='k', alpha=0.15, rasterized=True, zorder=99)\n",
    "ax_this.yaxis.set_ticks_position('both')\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.savefig('/home/yymao/Downloads/sat_lf.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}