F05-Survey completeness.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": [
    "from pathlib import Path\n",
    "\n",
    "import SAGA\n",
    "from SAGA import ObjectCuts as C\n",
    "from SAGA import HostCuts as H\n",
    "from easyquery import *\n",
    "from SAGA.utils import *\n",
    "from SAGA.utils.display import show_images\n",
    "from SAGA.utils.distance import *\n",
    "from astropy.table import vstack, Table, join\n",
    "from astropy.coordinates import Distance, SkyCoord\n",
    "from SAGA.objects.build2 import add_surface_brightness\n",
    "\n",
    "import numexpr as ne\n",
    "from scipy.stats import binned_statistic, norm\n",
    "#from bottleneck import nanmax, nanmin\n",
    "from itertools import product, chain\n",
    "\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": [
    "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 = Query(C.basic_cut2).filter(base)\n",
    "base = annotate_base(base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_other = SAGA.database.FitsTable(\"/home/yymao/Documents/Research/SAGA/PaperII/data-archive/saga_base_other.fits\").read()\n",
    "base_other = Query(C.basic_cut2).filter(base_other)\n",
    "base_other = annotate_base(base_other)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1_hosts = saga.host_catalog.load(include_stats=True, query=\"paper1_complete\")[\"HOSTID\"]\n",
    "p2_hosts = saga.host_catalog.load(include_stats=True, query=\"paper2_complete\")[\"HOSTID\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbins = np.linspace(-20.25, -11.75, 18)\n",
    "#rbins = np.linspace(16.75, 21, )\n",
    "r = 0.5*(rbins[1:]+rbins[:-1])\n",
    "\n",
    "mag_label = \"Mr_\"\n",
    "\n",
    "hostid = []\n",
    "X = []\n",
    "for base_this in chain(base.group_by(\"HOSTID\").groups, base_other.group_by(\"HOSTID\").groups):\n",
    "    hostid.append(base_this[\"HOSTID\"][0])\n",
    "    for q in (\n",
    "        C.paper2_targeting_cut,\n",
    "        C.paper2_targeting_cut & C.has_spec,\n",
    "        C.paper1_targeting_cut,\n",
    "        C.paper1_targeting_cut & C.has_spec,\n",
    "        C.paper2_targeting_cut & (~C.has_spec) & \"human_selected > 0\",\n",
    "        C.paper1_targeting_cut & (~C.has_spec) & \"human_selected > 0\",\n",
    "        C.paper2_targeting_cut & C.has_spec & (~C.has_our_specs_only),\n",
    "    ):\n",
    "        X.extend(np.histogram(q.filter(base_this, mag_label), rbins)[0])\n",
    "X = np.array(X, dtype=np.float64).reshape(len(hostid), -1, len(r))\n",
    "hostid = np.array(hostid)\n",
    "print(X.shape)\n",
    "\n",
    "with np.errstate(invalid=\"ignore\"):\n",
    "    X_ratio = np.nan_to_num(X[:,1] / X[:,0], nan=1)\n",
    "    X_ratio_pre = np.nan_to_num(X[:,6] / X[:,0], nan=1)\n",
    "    X_sum_ratio = X[:,1].sum(axis=1) / X[:,0].sum(axis=1) \n",
    "    X_sum = X[:,0].sum(axis=1) - X[:,1].sum(axis=1)\n",
    "    X_sum_lsbg = X[:,5].sum(axis=1) - X[:,4].sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1_mask = np.in1d(hostid, p1_hosts)\n",
    "p2_mask = np.in1d(hostid, p2_hosts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=2, figsize=(10, 3.7))\n",
    "\n",
    "ax_this = ax[0]\n",
    "\n",
    "for q, color, label, marker in (\n",
    "    (~(p1_mask | p2_mask), 'C0', \"Incomplete systems ({})\", 'o'),\n",
    "    (p2_mask & (~p1_mask), 'C3', 'New complete systems ({})', 's'),\n",
    "    (p1_mask, 'C1', 'Paper I complete systems ({})', '^'),\n",
    "):\n",
    "    n = np.count_nonzero(q)\n",
    "    ax_this.errorbar(X_sum_ratio[q], X_sum[q]*10**(np.random.randn(n)*0.01), \n",
    "                     marker=marker, ls='', label=label.format(n), ms=3.3, color=color, \n",
    "                    # yerr=(np.zeros(n), X_sum_lsbg[q]),\n",
    "                    )\n",
    "ax_this.set_yscale(\"log\")\n",
    "ax_this.legend(handletextpad=0.1, fontsize=13, frameon=True)\n",
    "ax_this.set_ylim(1, 500)\n",
    "ax_this.axvline(0.8, color=\"grey\", ls=\"--\")\n",
    "ax_this.set_ylabel('# objects w/o redshift per host')\n",
    "ax_this.set_xlabel('Spec coverage in primary targeting region')\n",
    "\n",
    "ax_this = ax[1]\n",
    "Y0 = np.mean(X_ratio[p2_mask], axis=0)\n",
    "Yerr = np.percentile(X_ratio[p2_mask], norm.sf([1, -1])*100, axis=0)\n",
    "\n",
    "ax_this.errorbar(r, Y0, (Y0-Yerr[0], Yerr[1]-Y0), ls=\":\", marker='o', color=\"C2\", label=\"Spec coverage\", zorder=8)\n",
    "ax_this.errorbar(r, np.mean(X_ratio_pre, axis=0), ls=\"-\", marker='', color=\"C8\", label=\"Spec coverage\\nbefore SAGA\")\n",
    "\n",
    "ax_this.set_ylim(0, 1.05)\n",
    "ax_this.set_xlabel('$M_{r,o}$')\n",
    "ax_this.set_ylabel('Spec coverage (primary region)')\n",
    "ax_this.axhline(1, color='grey', lw=0.5)\n",
    "ax_this.axvspan(-12.3, -11.75, color='grey', alpha=0.7, zorder=99)\n",
    "\n",
    "ax_this.set_xlim(-20.25, -11.75)\n",
    "Y0 = np.mean(X[p2_mask,0]-X[p2_mask,1], axis=0)\n",
    "Yerr = np.percentile(X[p2_mask,0]-X[p2_mask,1], norm.sf([1, -1])*100, axis=0)\n",
    "\n",
    "ax_second = plt.twinx(ax_this)\n",
    "ax_second.bar(r, height=Y0, width=np.ediff1d(rbins), color=\"C4\", alpha=0.3, label=\"No $z$; $\\\\mathcal{R}_\\\\mathrm{sat}\\\\sim0.005$\")\n",
    "ax_second.bar(r, height=np.mean(X[p2_mask,5], axis=0), width=np.ediff1d(rbins), color=\"C4\", alpha=0.7, label=\"No $z$, dLSBG; $\\\\mathcal{R}_\\\\mathrm{sat}\\\\sim0.3$\")\n",
    "\n",
    "ax_second.errorbar(r, Y0, (Y0-Yerr[0], Yerr[1]-Y0), ls=\"\", marker='', color=\"C4\", alpha=0.6)\n",
    "\n",
    "ax_second.set_ylim(0, 10.5)\n",
    "ax_second.set_ylabel('# objects w/o redshift per host')\n",
    "\n",
    "ax_this.legend(loc=\"lower left\", bbox_to_anchor=(0, 0.25), fontsize=13, handlelength=1, title=\"For complete systems\")\n",
    "ax_second.legend(loc=\"lower left\",  bbox_to_anchor=(0, 0), fontsize=13, handlelength=1)\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.savefig(\"/home/yymao/Downloads/survey-completeness.pdf\", bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## calculate average R_sat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "basic = Query(C.r_abs_limit, ~C.has_spec, saga.host_catalog.construct_host_query(\"paper2_complete\"))\n",
    "sel1 = Query(basic, \"human_selected == 0\", C.paper2_targeting_cut)\n",
    "sel2 = Query(basic, \"human_selected > 0\")\n",
    "\n",
    "print(sel1.filter(base, \"p_sat_corrected\").mean())\n",
    "print(sel2.filter(base, \"p_sat_corrected\").mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}