P03-LF+QF.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": [
    "nsat_missing_all = int(np.ceil(Query(~C.has_spec, C.faint_end_limit).filter(base, \"p_sat_corrected\").sum()))\n",
    "nsat_missing_limit = int(np.ceil(Query(~C.has_spec, C.r_abs_limit).filter(base, \"p_sat_corrected\").sum()))\n",
    "print(nsat_missing_all, nsat_missing_limit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sats = C.is_sat.filter(base)\n",
    "sats_maybe_all = Query(~C.has_spec, \"p_sat_corrected > 0\").filter(base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 1177 #923 #444 #177 # 114 # 47\n",
    "while True:\n",
    "    #seed += 1\n",
    "    sats_maybe = Query(\n",
    "        (lambda p: p > np.random.RandomState(seed).rand(len(p)), \"p_sat_corrected\"),  # must put in 1st place\n",
    "        ~C.has_spec, \n",
    "    ).filter(base)\n",
    "\n",
    "    if (np.abs((C.faint_end_limit.count(sats_maybe) - nsat_missing_all) < 1.5 and\n",
    "        C.r_abs_limit.count(sats_maybe)) == nsat_missing_limit and\n",
    "        (~Query(C.high_priority_sb_tight, C.gr_cut_tight, \"Mr > -15.5\")).count(sats_maybe) == 0):\n",
    "        print(seed)\n",
    "        break\n",
    "    else:\n",
    "        seed += 1\n",
    "        \n",
    "print(C.r_abs_limit.count(sats_maybe), len(sats_maybe))"
   ]
  },
  {
   "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": [
    "fq_lg = ascii.read(format=\"fast_tab\", names=[\"log_sm\", \"fq\", \"fq_ue\", \"fq_le\"], data_start=0,\n",
    "    table=\"\"\"3.49\t1.00\t1.00\t0.77\n",
    "4.49\t1.00\t1.00\t0.90\n",
    "5.49\t0.94\t0.96\t0.83\n",
    "6.49\t1.00\t1.00\t0.77\n",
    "7.49\t0.80\t0.88\t0.54\n",
    "8.49\t0.66\t0.81\t0.38\n",
    "9.18\t0.00\t0.60\t0.00\n",
    "\"\"\")\n",
    "\n",
    "fq_tinker = ascii.read(format=\"fast_tab\", names=[\"log_sm0\", \"fq\", \"log_sm1\",  \"fq_ue\", \"log_sm2\", \"fq_le\"], data_start=0,\n",
    "    table=\"\"\"9.7632\t0.0625\t9.7655\t0.0455\t9.7633\t0.0788\n",
    "9.8511\t0.1149\t9.8510\t0.0986\t9.8511\t0.1312\n",
    "9.9459\t0.1340\t9.9459\t0.1170\t9.9459\t0.1503\n",
    "10.0479\t0.1835\t10.0479\t0.1637\t10.0480\t0.2019\n",
    "10.1477\t0.2514\t10.1476\t0.2295\t10.1477\t0.2734\n",
    "10.2449\t0.2790\t10.2472\t0.2563\t10.2473\t0.3010\n",
    "10.3494\t0.3618\t10.3493\t0.3363\t10.3495\t0.3859\n",
    "10.4468\t0.4433\t10.4467\t0.4142\t10.4445\t0.4709\n",
    "10.5488\t0.4666\t10.5487\t0.4340\t10.5489\t0.4985\n",
    "10.6485\t0.5381\t10.6484\t0.5033\t10.6486\t0.5728\n",
    "10.7458\t0.5812\t10.7481\t0.5437\t10.7483\t0.6188\n",
    "10.8503\t0.6471\t10.8501\t0.5982\t10.8504\t0.6938\n",
    "10.9523\t0.7079\t10.9522\t0.6526\t10.9525\t0.7632\n",
    "11.0495\t0.7249\t11.0493\t0.6575\t11.0521\t0.7929\n",
    "11.1397\t0.7517\t11.1394\t0.6695\t11.1399\t0.8346\n",
    "11.2347\t0.8211\t11.2319\t0.6943\t11.2327\t0.9472\n",
    "\"\"\")\n",
    "fq_tinker[\"log_sm\"] = (fq_tinker[\"log_sm0\"] + fq_tinker[\"log_sm1\"] + fq_tinker[\"log_sm2\"]) / 3\n",
    "\n",
    "fq_geha = ascii.read(format=\"fast_tab\", names=[\"log_sm\", \"fq\",\"fq_e\"], data_start=0,\n",
    "    table=\"\"\"9.9\t0.084\t0.0052\n",
    "9.7\t0.031\t0.0032\n",
    "9.5\t0.017\t0.0023\n",
    "9.3\t0.007\t0.0016\n",
    "9.1\t0.002\t0.0014\n",
    "8.9\t0.000\t0.0019\n",
    "8.7\t0.000\t0.0028\n",
    "8.5\t0.000\t0.0043\n",
    "8.3\t0.000\t0.0065\n",
    "8.1\t0.000\t0.0091\n",
    "7.8\t0.000\t0.0093\n",
    "7.3\t0.000\t0.0193\n",
    "\"\"\")\n",
    "fq_geha[\"fq_ue\"] = fq_geha[\"fq\"] + fq_geha[\"fq_e\"]\n",
    "fq_geha[\"fq_le\"] = np.maximum(0, fq_geha[\"fq\"] - fq_geha[\"fq_e\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=2, figsize=(10, 4.3))\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",
    "\n",
    "\n",
    "ax_this.fill_between(fq_tinker[\"log_sm\"], fq_tinker[\"fq_le\"], fq_tinker[\"fq_ue\"], color=\"C4\", alpha=0.7, lw=0, rasterized=True)\n",
    "ax_this.fill_between(fq_geha[\"log_sm\"], fq_geha[\"fq_le\"], fq_geha[\"fq_ue\"], color=\"C4\", alpha=0.7, rasterized=True)\n",
    "\n",
    "\n",
    "ax_this.errorbar(fq_lg[\"log_sm\"], fq_lg[\"fq\"], marker='x', ls=\"\", c=\"C1\", ms=6)\n",
    "ax_this.fill_between(fq_lg[\"log_sm\"], fq_lg[\"fq_le\"], fq_lg[\"fq_ue\"], color=\"C1\", alpha=0.25, lw=0, rasterized=True)\n",
    "\n",
    "pbins = np.linspace(0, 100, 8)\n",
    "log_sm = np.concatenate([sats['log_sm']])\n",
    "bins = np.percentile(log_sm, pbins)\n",
    "x = np.percentile(log_sm, (pbins[1:]+pbins[:-1])*0.5)\n",
    "\n",
    "interloper_corr = lambda d: 0.95 - d*0.001\n",
    "\n",
    "n_q = np.histogram(Query('EW_Halpha < 2').filter(sats, 'log_sm'), bins)[0]\n",
    "sats_q = Query('EW_Halpha < 2').filter(sats, [\"log_sm\", \"RHOST_KPC\"])\n",
    "sats_maybe_q = Query().filter(sats_maybe_all, [\"log_sm\", \"RHOST_KPC\", \"p_sat_corrected\"])\n",
    "n_q_corr = np.histogram(sats_q['log_sm'], bins, weights=interloper_corr(sats_q[\"RHOST_KPC\"]))[0]\n",
    "n_q_corr += np.histogram(sats_maybe_q['log_sm'], bins, weights=sats_maybe_q[\"p_sat_corrected\"]*interloper_corr(sats_maybe_q[\"RHOST_KPC\"]))[0]\n",
    "\n",
    "n_all = np.histogram(sats['log_sm'], bins)[0]\n",
    "n_all_corr = np.histogram(sats['log_sm'], bins, weights=interloper_corr(sats[\"RHOST_KPC\"]))[0]\n",
    "n_all_corr += np.histogram(sats_maybe_all['log_sm'], bins, weights=sats_maybe_all[\"p_sat_corrected\"]*interloper_corr(sats_maybe_all[\"RHOST_KPC\"]))[0]\n",
    "\n",
    "p = n_q/n_all\n",
    "pb = (n_q+1)/(n_all+2)\n",
    "perr = np.sqrt(pb*(1-pb) / n_all)\n",
    "p1 = np.minimum(1, p + perr)\n",
    "p2 = np.maximum(0, p - perr)\n",
    "p_corr = np.maximum(0, (n_q_corr/n_all_corr) - p)\n",
    "\n",
    "ax_this.errorbar(x, p, yerr=(p-p2, p1-p), ls='', marker='o', lw=2, capsize=2, ms=5, c=\"C2\")\n",
    "ax_this.errorbar(x, p, yerr=(np.zeros_like(p), p1-p+p_corr), ls='', marker='', alpha=0.6, c=\"C2\", lw=2, rasterized=True)\n",
    "\n",
    "\n",
    "\n",
    "this_work_label = \"This work\" and \"Mao+20\"\n",
    "ax_this.set_xlim(5.3, 11)\n",
    "ax_this.text(6.75, 0.2, f\"SAGA sat.\\n({this_work_label})\", ha=\"right\", va=\"bottom\", color=\"C2\", fontsize=13, fontweight=\"bold\")\n",
    "ax_this.text(7.2, 0.91, \"MW+M31 sat. (Wetzel+15)\", ha=\"left\", va=\"bottom\", color=\"C1\", fontsize=13, fontweight=\"bold\")\n",
    "ax_this.text(10.7, 0.6, \"Field gal.\\n(Geha+12)\", ha=\"right\", va=\"bottom\", color=\"C4\", fontsize=13, fontweight=\"bold\")\n",
    "#plt.legend(loc=\"upper right\",  markerfirst=False, fontsize=13,frameon=True, handletextpad=0)\n",
    "ax_this.set_xlabel(\"$\\\\log\\;[M_*/M_\\\\odot]$\")\n",
    "ax_this.set_ylabel(\"Quenched fraction\")\n",
    "ax_this.set_ylim(-0.05, 1.05)\n",
    "ax_this.axhline(0, lw=0.5, color=\"grey\")\n",
    "ax_this.axhline(1, lw=0.5, color=\"grey\")\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.savefig('/home/yymao/Downloads/sat_lf_qf.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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}