lsst · czwa · Apr 13, 2023 · Apr 5, 2023 · Apr 5, 2023 · Apr 11, 2023
diff --git a/.gitignore b/.gitignore
@@ -21,4 +21,4 @@ version.py
 bin/
 .pytest_cache
 pytest_session.txt
-examples/.ipynb_checkpoints
+notebooks/.ipynb_checkpoints
diff --git a/examples/cpVerifyBfk.ipynb b/examples/cpVerifyBfk.ipynb
diff --git a/examples/cpVerifyBias.ipynb b/examples/cpVerifyBias.ipynb
diff --git a/examples/cpVerifyCrosstalk.ipynb b/examples/cpVerifyCrosstalk.ipynb
diff --git a/examples/cpVerifyDark.ipynb b/examples/cpVerifyDark.ipynb
diff --git a/examples/cpVerifyDefects.ipynb b/examples/cpVerifyDefects.ipynb
diff --git a/examples/cpVerifyFlat.ipynb b/examples/cpVerifyFlat.ipynb
diff --git a/examples/cpVerifyGainFromFlatPair.ipynb b/examples/cpVerifyGainFromFlatPair.ipynb
diff --git a/examples/cpVerifyLinearity.ipynb b/examples/cpVerifyLinearity.ipynb
diff --git a/examples/cpVerifyPtc.ipynb b/examples/cpVerifyPtc.ipynb
diff --git a/examples/cpVerifyPtcV2.ipynb b/examples/cpVerifyPtcV2.ipynb
diff --git a/examples/cpVerifyPtcWithLinearity.ipynb b/examples/cpVerifyPtcWithLinearity.ipynb
diff --git a/notebooks/cpVerifyBfk.ipynb b/notebooks/cpVerifyBfk.ipynb
@@ -0,0 +1,392 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8aa75eb3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import scipy as sp\n",
+    "\n",
+    "import lsst.daf.butler as dB\n",
+    "import lsst.cp.verify.notebooks.utils as utils\n",
+    "import lsst.afw.display as afwDisplay\n",
+    "plt.rcParams['figure.figsize'] = (8, 8)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "84dccd11",
+   "metadata": {
+    "tags": [
+     "parameters"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# This cell contains parameters that can be automatically set via the papermill package.\n",
+    "# Examples:\n",
+    "#  Update parameters in input.ipynb, writing output.ipynb, but do not execute:\n",
+    "#   papermill --prepare-only -p calibType newBias -p cameraName LSSTCam <input> <output>\n",
+    "#  Disable interactive cells in input.ipynb, execute it, and write output.ipynb.\n",
+    "#   papermill -p interactive False <input> <output>\n",
+    "interactive = True\n",
+    "\n",
+    "# Which repository to use.\n",
+    "repository = '/repo/main/'\n",
+    "\n",
+    "# Which calibration type to analyse.\n",
+    "calibType = 'bfk'\n",
+    "detectorNumber = 0\n",
+    "exposureNumber = 2021021700347\n",
+    "\n",
+    "# Which camera the calibration is for.\n",
+    "cameraName = 'LATISS'\n",
+    "\n",
+    "# Which display to use.\n",
+    "displayBackend = 'astrowidgets'\n",
+    "\n",
+    "# Which collection the calibration was constructed in.\n",
+    "genCollection = 'u/czw/DM-28920/bfkGen.20210928-exp'\n",
+    "\n",
+    "# Which collection the PTC data was constructed in.\n",
+    "ptcCollection = 'u/czw/DM-28920/ptcGen.20210723Ya'\n",
+    "\n",
+    "# Which collection containing the verification outputs.\n",
+    "correctedCollection = 'u/czw/DM-30172/bfkV.20211102Xa'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "01855373",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get butler and camera\n",
+    "butler = dB.Butler(repository, collections=[genCollection, ptcCollection,\n",
+    "                                            correctedCollection])\n",
+    "camera = butler.get('camera', instrument=cameraName)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e535f6b6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get the calibrations.\n",
+    "bfk = butler.get(calibType, instrument=cameraName, detector=detectorNumber,\n",
+    "                 exposure=exposureNumber,\n",
+    "                 collections=genCollection)\n",
+    "ptc = butler.get('ptc', instrument=cameraName, detector=detectorNumber,\n",
+    "                 exposure=exposureNumber,\n",
+    "                 collections=ptcCollection)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "486dabee",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "runStats = butler.get('verifyBfkStats', instrument=cameraName, \n",
+    "                      collections=correctedCollection)\n",
+    "runSuccess = runStats.pop('SUCCESS')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a695fc25",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Display summary table of tests and failure counts.\n",
+    "utils.failureTable(runStats)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0071e01a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot raw cross-correlation values\n",
+    "inputSelection = 1\n",
+    "iterator = 0\n",
+    "fig, axes = plt.subplots(nrows=4, ncols=4,\n",
+    "                         sharex=True, sharey=False, figsize=[8.0, 8.0], \n",
+    "                         constrained_layout=True)\n",
+    "for ampName in sorted(set(ptc.ampNames)):\n",
+    "    covariance = ptc.covariances[ampName][inputSelection]\n",
+    "    axes[iterator // 4, iterator % 4].set_aspect(\"equal\")\n",
+    "    axes[iterator // 4, iterator % 4].imshow(np.log10(np.abs(covariance)))\n",
+    "    axes[iterator // 4, iterator % 4].set_title(f\"Amp {ampName} {covariance[0][0]:.2e}\")\n",
+    "\n",
+    "    iterator += 1\n",
+    "fig.suptitle(f\"Measured covariance {ptc.inputExpIdPairs[ampName][inputSelection]}\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "09b866f1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot mean cross-correlation values\n",
+    "iterator = 0\n",
+    "fig, axes = plt.subplots(nrows=4, ncols=4,\n",
+    "                         sharex=True, sharey=False, figsize=[8.0, 8.0], \n",
+    "                         constrained_layout=True)\n",
+    "for ampName in sorted(set(ptc.ampNames)):\n",
+    "    axes[iterator // 4, iterator % 4].set_aspect(\"equal\")\n",
+    "    axes[iterator // 4, iterator % 4].imshow(np.log10(np.abs(bfk.meanXcorrs[ampName])))\n",
+    "    axes[iterator // 4, iterator % 4].set_title(f\"Amp {ampName} {bfk.meanXcorrs[ampName][8, 8]:.2e}\")\n",
+    "\n",
+    "    iterator += 1\n",
+    "fig.suptitle(\"Mean cross-correlations\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2ba69bdf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot cuts across mean cross-correlations\n",
+    "iterator = 0\n",
+    "fig, axes = plt.subplots(nrows=4, ncols=4,\n",
+    "                         sharex=True, sharey=False, figsize=[8.0, 8.0], \n",
+    "                         constrained_layout=True)\n",
+    "for ampName in sorted(set(ptc.ampNames)):\n",
+    "    axes[iterator // 4, iterator % 4].set_title(f\"Amp {ampName} {bfk.meanXcorrs[ampName][8, 8]:.2e}\")\n",
+    "\n",
+    "    axes[iterator // 4, iterator % 4].step(range(0, 17), bfk.meanXcorrs[ampName][:, 8],\n",
+    "                                           where='mid', label='Y-axis')\n",
+    "    axes[iterator // 4, iterator % 4].step(range(0, 17), bfk.meanXcorrs[ampName][8, :],\n",
+    "                                           where='mid', label='X-axis')\n",
+    "    axes[iterator // 4, iterator % 4].set_xlabel(\"Location (pxl)\")\n",
+    "    axes[iterator // 4, iterator % 4].set_ylabel(\"Mean cross-corr\")\n",
+    "    axes[iterator // 4, iterator % 4].legend()\n",
+    "    iterator += 1\n",
+    "fig.suptitle(\"Central cuts of mean cross-correlation\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "da8d8e91",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot average kernel <-- This is what gets used in IsrTask\n",
+    "detector = camera[0]\n",
+    "detName = detector.getName()\n",
+    "bfk.makeDetectorKernelFromAmpwiseKernels(detName)\n",
+    "averageKernel = bfk.detKernels[detName]\n",
+    "plt.imshow(averageKernel)\n",
+    "plt.title(\"Detector averaged kernel\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b05a9c4a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot central cuts of detector averaged kernel.\n",
+    "plt.step(range(0, 17), bfk.detKernels[detName][:, 8], \n",
+    "         where='mid', label='Y-axis')\n",
+    "plt.step(range(0, 17), bfk.detKernels[detName][8, :], \n",
+    "         where='mid', label='X-axis')\n",
+    "plt.legend()\n",
+    "plt.xlabel(\"Location (pxl)\")\n",
+    "plt.ylabel(\"Kernel value\")\n",
+    "plt.title(\"Central cuts of detector averaged kernel\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3b15ae8d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot amp kernels\n",
+    "iterator = 0\n",
+    "fig, axes = plt.subplots(nrows=4, ncols=4,\n",
+    "                         sharex=True, sharey=False, figsize=[8.0, 8.0], \n",
+    "                         constrained_layout=True)\n",
+    "for ampName in sorted(set(ptc.ampNames)):\n",
+    "    axes[iterator // 4, iterator % 4].set_aspect(\"equal\")\n",
+    "    axes[iterator // 4, iterator % 4].imshow(bfk.ampKernels[ampName])\n",
+    "    axes[iterator // 4, iterator % 4].set_title(f\"Amp {ampName} {bfk.ampKernels[ampName][8, 8]:.2e}\")\n",
+    "\n",
+    "    iterator += 1\n",
+    "fig.suptitle(\"Amplfier level kernels\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "636d339d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot cuts across the amp kernels\n",
+    "iterator = 0\n",
+    "fig, axes = plt.subplots(nrows=4, ncols=4,\n",
+    "                         sharex=True, sharey=False, figsize=[8.0, 8.0], \n",
+    "                         constrained_layout=True)\n",
+    "for ampName in sorted(set(ptc.ampNames)):\n",
+    "    axes[iterator // 4, iterator % 4].set_title(f\"Amp {ampName} {bfk.ampKernels[ampName][8, 8]:.2e}\")\n",
+    "\n",
+    "    axes[iterator // 4, iterator % 4].step(range(0, 17), bfk.ampKernels[ampName][:, 8], \n",
+    "                                           where='mid', label='Y-axis')\n",
+    "    axes[iterator // 4, iterator % 4].step(range(0, 17), bfk.ampKernels[ampName][8, :], \n",
+    "                                           where='mid', label='X-axis')\n",
+    "    axes[iterator // 4, iterator % 4].set_xlabel(\"Location (pxl)\")\n",
+    "    axes[iterator // 4, iterator % 4].set_ylabel(\"Mean cross-corr\")\n",
+    "    axes[iterator // 4, iterator % 4].legend()\n",
+    "\n",
+    "    iterator += 1\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "70c079cd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get object information and sizes.\n",
+    "Fall = []\n",
+    "Tall = []\n",
+    "TuAll = []\n",
+    "for detector in camera:\n",
+    "    detId = detector.getId()\n",
+    "    for exposureId, stats in runStats.items():\n",
+    "        icSrc = butler.get('icSrc', collections=correctedCollection,\n",
+    "                           visit=exposureId, detector=detId)\n",
+    "        icSrcUn = butler.get('verifyUncorrBfkSrc', collections=correctedCollection, \n",
+    "                             visit=exposureId, detector=detId)\n",
+    "        if len(icSrc) != len(icSrcUn):\n",
+    "            continue\n",
+    "        T = icSrc['base_SdssShape_xx'] + icSrc['base_SdssShape_yy']\n",
+    "        Tp = icSrc['base_SdssShape_psf_xx'] + icSrc['base_SdssShape_psf_yy']\n",
+    "        F = (icSrc['base_SdssShape_instFlux'])\n",
+    "        Tu = icSrcUn['base_SdssShape_xx'] + icSrcUn['base_SdssShape_yy']\n",
+    "        Tup = icSrcUn['base_SdssShape_psf_xx'] + icSrcUn['base_SdssShape_psf_yy']\n",
+    "        plt.xlabel(\"InstFlux\")\n",
+    "        plt.ylabel(\"Uncorrected size - corrected size\")\n",
+    "        plt.ylim(-10, 10)\n",
+    "        plt.scatter(F, 100*(Tu - T), label=exposureId, s=1.5)\n",
+    "        Fall.extend(F)\n",
+    "        Tall.extend(T)\n",
+    "        TuAll.extend(Tu)\n",
+    "    plt.legend()\n",
+    "    plt.xlim(0, 1e6)\n",
+    "    plt.title(detector.getName())\n",
+    "    plt.grid()\n",
+    "    plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e3e3eb0a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.optimize import least_squares\n",
+    "def f(x, a, b, c):\n",
+    "    return a * np.exp(b * (x + c))\n",
+    "\n",
+    "def residual(p, x, y):\n",
+    "    return y - f(x, *p)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "384d8777",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot functional fits.\n",
+    "XX = np.arange(-20, -11, 0.1)\n",
+    "\n",
+    "for exposureId in runStats.keys():\n",
+    "    for detector in camera:\n",
+    "        detStats = butler.get(\"verifyBfkDetStats\", instrument=cameraName, \n",
+    "                              detector=detector.getId(),\n",
+    "                              exposure=exposureId, visit=exposureId)\n",
+    "        print(detStats['CATALOG']['BRIGHT_SLOPE'])\n",
+    "        M = detStats['CATALOG']['MAGNITUDES']\n",
+    "        D = detStats['CATALOG']['SIZE_DIFF']\n",
+    "        mask = np.isfinite(M) * np.isfinite(D)\n",
+    "        \n",
+    "        p0 = [0.0, 0.0, 0.0]\n",
+    "        ZZ = least_squares(residual, p0, args=(np.array(M)[mask], np.array(D)[mask]), \n",
+    "                           loss='cauchy')\n",
+    "        pf = np.polyfit(np.array(M)[mask], np.array(D)[mask], 2)\n",
+    "        \n",
+    "        plt.scatter(detStats['CATALOG']['MAGNITUDES'], detStats['CATALOG']['SIZE_DIFF'])\n",
+    "        plt.plot(XX, np.polyval(pf, XX), label=f\"polynomial (N=2) {pf}\")\n",
+    "        plt.plot(XX, f(XX, *ZZ.x), label=f\"exponential {ZZ.x}\")\n",
+    "        plt.ylim(-1, 5)\n",
+    "        plt.legend()\n",
+    "\n",
+    "        plt.title(exposureId)\n",
+    "        plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7398f2a0",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "LSST",
+   "language": "python",
+   "name": "lsst"
+  },
+  "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.10.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}