{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c356285e-0524-4bc4-a747-5fce661e0434",
   "metadata": {},
   "source": [
    "## Run wcs-based matching to compare object tables to metadetect no-shear catalogs\n",
    "\n",
    "This takes as the metadetect no-shear catalog and object catalogs as inputs and runs matching using the wcs-based `WcsMatch` matcher."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5a4b6b3-b696-426d-b20f-3aa5149cb7a6",
   "metadata": {},
   "source": [
    "#### Standard imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fb9df94-02bf-4788-98b5-d289f0da8b38",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hpmcm\n",
    "import tables_io\n",
    "import glob\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d57ab36-f67a-40a2-abea-e62fc81b193a",
   "metadata": {},
   "source": [
    "#### Set up the configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17cd4820-7642-4f16-b226-8453adb170b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "DATADIR = \"test_data\"   # Input data directory\n",
    "shear_st = \"0p01\"       # Applied shear as a string\n",
    "shear = 0.01            # Decimal version of applied shear\n",
    "shear_type = \"wmom\"     # which object characterization to use \n",
    "tract = 10463           # which tract to study\n",
    "\n",
    "REF_DIR = (37.9, 7.0)         # RA, DEC in deg of center of match region\n",
    "REGION_SIZE = (0.375, 0.375)  # Size of match region in degrees\n",
    "PIXEL_SIZE = 0.5/3600.        # Size of pixels used in matching\n",
    "PIXEL_R2CUT = 4.              # Cut at distance**2 = 4 pixels\n",
    "PIXEL_MATCH_SCALE = 1         # Use pixel scale to do matching\n",
    "\n",
    "SOURCE_TABLEFILES = sorted(glob.glob(os.path.join(DATADIR, f\"shear_*_{shear_st}_cleaned_{tract}_ns.pq\")))\n",
    "SOURCE_TABLEFILES.append(os.path.join(DATADIR, f\"object_{tract}.pq\"))\n",
    "SOURCE_TABLEFILES.reverse()\n",
    "SOURCE_TABLEFILES = [SOURCE_TABLEFILES[0], SOURCE_TABLEFILES[1]]\n",
    "VISIT_IDS = np.arange(len(SOURCE_TABLEFILES))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef87000f-d932-44da-8fcf-5006a82f6232",
   "metadata": {},
   "source": [
    "#### Make the matcher, reduce the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4393270d-6675-4cb3-a516-033516304034",
   "metadata": {},
   "outputs": [],
   "source": [
    "matcher = hpmcm.WcsMatch.create(REF_DIR, REGION_SIZE, pixel_size=PIXEL_SIZE, pixel_r2_cut=PIXEL_R2CUT)\n",
    "matcher.reduceData(SOURCE_TABLEFILES, VISIT_IDS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38da6ab5-1a40-41a6-b463-8b32c4be1f14",
   "metadata": {},
   "source": [
    "#### Make a plot comparing the signal-to-noise in the two catalogs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a29be2c-6089-49d2-b772-5030d553a625",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt.hist(matcher.full_data[0].snr, bins=np.logspace(0, 4, 81), alpha=0.4, label=\"obj\")\n",
    "_ = plt.hist(matcher.full_data[1].snr, bins=np.logspace(0, 4, 81), alpha=0.4, label=\"pgauss\")\n",
    "_ = plt.xscale('log')\n",
    "_ = plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7823e1bd-c58a-49ab-9fcd-0bf9215db023",
   "metadata": {},
   "source": [
    "#### This should have made 3 x 3 cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c630f6e2-02b0-4527-a38f-5941040a8fd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "matcher.n_cell"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8de684e9-2bbd-4382-b74e-9bda78fb2dce",
   "metadata": {},
   "source": [
    "#### Run the loop over cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c731ef03-2f36-4f5e-8dd2-4520ada4213d",
   "metadata": {},
   "outputs": [],
   "source": [
    "matcher.analysisLoop()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24511bb4-35c2-4a18-b51e-1fcab9d448e9",
   "metadata": {},
   "source": [
    "#### Show a single cluster\n",
    "\n",
    "The x and y axes here are the in the cluster frame for a single cluster.\n",
    "The color scale shows the number of sources per/pixel.\n",
    "\n",
    "The `x` markers are the original source postions.   The `o` makters are the deshear positions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d094934-dcc7-4268-aa59-166860203ce0",
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = matcher.cell_dict[matcher.getCellIdx(2,2)]\n",
    "od = cell.analyze(None, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eada054c-06a2-4adc-9598-30468c52b7e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster = list(cell.cluster_dict.values())[0]\n",
    "_ = hpmcm.viz_utils.showCluster(od['image'], cluster, cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62a8f605-1b15-4d0e-8f78-3e1fe501a266",
   "metadata": {},
   "source": [
    "#### Classify the objects by match type\n",
    "\n",
    "This looks at the characteristics of the matched objects and categorizes them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d7df2d6-f60e-4392-a500-ecdd5ccfaea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj_lists = hpmcm.classify.classifyObjects(matcher, snr_cut=10)\n",
    "hpmcm.classify.printObjectTypes(obj_lists)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12912754-da82-43a0-8671-d4515efdc130",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_good = len(obj_lists['ideal'])\n",
    "bad_list = ['edge_mixed', 'edge_missing', 'edge_extra', 'orphan', 'missing', 'two_missing', 'many_missing', 'extra', 'caught']\n",
    "n_bad = np.sum([len(obj_lists[x]) for x in bad_list])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f4a90b8-5088-4077-bbda-b73d0d640b5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "effic = n_good/(n_good+n_bad)\n",
    "effic_err = np.sqrt(effic*(1-effic)/(n_good+n_bad))\n",
    "print(f\"Effic: {effic:.5} +- {effic_err:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7b6d6e5-d8e5-4ae0-a5ad-5096a9681538",
   "metadata": {},
   "source": [
    "#### Classify objcts using the object table as the reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b635ac31-3488-454a-8fdf-3ec5a3c23130",
   "metadata": {},
   "outputs": [],
   "source": [
    "odict = hpmcm.classify.matchObjectsAgainstRef(matcher, snrCut=10.)\n",
    "hpmcm.classify.printObjectMatchTypes(odict)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b579b36f-84bb-410f-bce1-bee04cab2db4",
   "metadata": {},
   "source": [
    "#### Display an object\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe6185dd-996b-43bb-8376-de8b70af7521",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = hpmcm.viz_utils.showShearObj(matcher, obj_lists['missing'][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3672eac6-66d9-4a73-9e8c-9ceb413efbb2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "deab44aa-c311-45ed-81f3-dd3155d68fbd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}