Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Analysis notebook for tie vs danish #9

Open
wants to merge 10 commits into
base: develop
Choose a base branch
from
Open

Analysis notebook for tie vs danish #9

Changes from 9 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
b99cfa4
Added notebook for tie vs danish.
jfcrenshaw Oct 28, 2024
973e47e
Plotting donut stamps.
jfcrenshaw Oct 28, 2024
6d4743b
Updated with new ISR runs.
jfcrenshaw Oct 28, 2024
cc077ac
More detailed obs_lsst info
jfcrenshaw Oct 28, 2024
6297f5b
sorted ids
jfcrenshaw Oct 28, 2024
dbe7ad8
New naming
jfcrenshaw Oct 28, 2024
7fcd0a8
Restructure directory.
jfcrenshaw Oct 28, 2024
4fb7992
Fixed stamp plotting bug
jfcrenshaw Oct 28, 2024
26230f2
Updated, including with same number of donuts for each detector.
jfcrenshaw Oct 29, 2024
73a844c
More quantitative statements.
jfcrenshaw Oct 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
382 changes: 382 additions & 0 deletions notebooks/WET/WET-007/WET-007_20241028_tie-vs-danish.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,382 @@
{
Copy link
Contributor

@suberlak suberlak Oct 28, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given that the ts_aos_analysis notebook does not have outputs, it may be a good idea ( beyond attaching PDF of the executed notebook in the PR, and in the ticket) to describe the result here in words (eg. for Z7 Danish result of .. and TIE of ... with estimated uncertainty... agree with manual results ....) . But it is hard to describe results without any pictures, and with fleeting persistency of repos / collections / packages, I wonder if we shouldn't have a repo where these notebooks are actually executed.. (ts_aos_analysis_exec ) ?

Perhaps rewrite to avoid tautology, eg. Note we employ means of combined Zernike estimates from each of the 9 detectors


Reply via ReviewNB

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Bryce pointed out that we can't post executed notebooks anywhere public because many of the them contain embargoed data (hence why he deleted the PDFs in the PRs). We have created #aos-notebook-review on Slack to make the review process easier, but I agree we might want somewhere to store the executed notebooks. A private Github repo could work, but the size of the history will quickly get cumbersome. Maybe instead we create a shared Google drive folder where we can mirror the directory structure of the repo and drop PDFs of executed notebooks? It's a little cumbersome, but we could make adding the PDF to the drive be the last stage of the PR process?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, the slack channel is not a bad idea, but it may end up requiring too much scrolling with "v1", "better v1", "better v1 after comments" .... If we can create a shared private google drive (or some other approved location) that would work better (in reply to review one can just overwrite).

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the place to have that is probably JIRA in the different tickets, also for traceability in the future for the project, it would be better than an ad-hoc solution that only AOS people can access

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gmegh we cannot put compiled PDFs in tickets because they frequently contain on-sky pixels

Copy link
Contributor

@suberlak suberlak Oct 29, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo "mostly gree".

And I think it strengthens the text by putting a numerical value ("noisiest, with 150nm rms"), and later about disagreements (~200 nm difference for Z5, with amplitude of ~1200 nm is ~16% , and for Z13,Z20...)


Reply via ReviewNB

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# WET-007: WEP Performance TIE vs Danish\n",
"\n",
"Owner: **John Franklin Crenshaw** <br>\n",
"Last Verified to Run: **2024-10-29** <br>\n",
"Software Version:\n",
" - `ts_wep`: **12.6.0**\n",
" - `lsst_distrib`: **w_2024_43**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test Description\n",
"\n",
"This notebook addresses [SITCOM-1149](https://rubinobs.atlassian.net/browse/SITCOM-1149) by comparing performance of TIE vs Danish algorithms on donuts from 20241027.\n",
"\n",
"Data for this analyis was created by running the following two commands on USDF:\n",
"\n",
"TIE:\n",
"```sh\n",
"pipetask run -j 9 -b /repo/embargo -i LSSTComCam/defaults -o u/crenshaw/WET-007_tie_20241027_3 -p $HOME/notebooks/comcam_commissioning/test_prep/WET-007/pipeline_WET-007_tie_maxSelect.yaml -d \"instrument='LSSTComCam' and exposure.observation_type='cwfs' and visit in (2024102700016,2024102700017)\" --rebase --register-dataset-types\n",
"```\n",
"\n",
"Danish:\n",
"```sh\n",
"pipetask run -j 9 -b /repo/embargo -i LSSTComCam/defaults -o u/crenshaw/WET-007_danish_20241027_3 -p $HOME/notebooks/comcam_commissioning/test_prep/WET-007/pipeline_WET-007_danish_maxSelect.yaml -d \"instrument='LSSTComCam' and exposure.observation_type='cwfs' and visit in (2024102700016,2024102700017)\" --rebase --register-dataset-types\n",
"```\n",
"\n",
"where the contents of `pipeline_WET-007_tie_maxSelect.yaml` are\n",
"```yaml\n",
"description: Pipeline for WET-007 using TIE\n",
"instrument: lsst.obs.lsst.LsstComCam\n",
"tasks:\n",
" generateDonutDirectDetectTask:\n",
" class: lsst.ts.wep.task.generateDonutDirectDetectTask.GenerateDonutDirectDetectTask\n",
" config:\n",
" donutSelector.useCustomMagLimit: True\n",
" donutSelector.sourceLimit: 20\n",
" cutOutDonutsScienceSensorGroupTask:\n",
" class: lsst.ts.wep.task.cutOutDonutsScienceSensorTask.CutOutDonutsScienceSensorTask\n",
" config:\n",
" python: |\n",
" from lsst.ts.wep.task.pairTask import GroupPairer\n",
" config.pairer.retarget(GroupPairer)\n",
" donutStampSize: 200\n",
" initialCutoutPadding: 40\n",
" calcZernikesTask:\n",
" class: lsst.ts.wep.task.calcZernikesTask.CalcZernikesTask\n",
" config:\n",
" estimateZernikes.maxNollIndex: 28\n",
" estimateZernikes.saveHistory: False\n",
" estimateZernikes.maskKwargs: {'doMaskBlends': False}\n",
" donutStampSelector.maxSelect: 5\n",
" isr:\n",
" class: lsst.ip.isr.IsrTaskLSST\n",
" config:\n",
" # Although we don't have to apply the amp offset corrections, we do want\n",
" # to compute them for analyzeAmpOffsetMetadata to report on as metrics.\n",
" doAmpOffset: true\n",
" ampOffset.doApplyAmpOffset: false\n",
" # Turn off slow steps in ISR\n",
" doBrighterFatter: false\n",
" # Mask saturated pixels,\n",
" # but turn off quadratic crosstalk because it's currently broken\n",
" doSaturation: True\n",
" crosstalk.doQuadraticCrosstalkCorrection: False\n",
"```\n",
"\n",
"and the contents of `pipeline_WET-007_danish_maxSelect.yaml` are\n",
"```yaml\n",
"description: Pipeline for WET-007 using Danish\n",
"instrument: lsst.obs.lsst.LsstComCam\n",
"tasks:\n",
" generateDonutDirectDetectTask:\n",
" class: lsst.ts.wep.task.generateDonutDirectDetectTask.GenerateDonutDirectDetectTask\n",
" config:\n",
" donutSelector.useCustomMagLimit: True\n",
" donutSelector.sourceLimit: 20\n",
" cutOutDonutsScienceSensorGroupTask:\n",
" class: lsst.ts.wep.task.cutOutDonutsScienceSensorTask.CutOutDonutsScienceSensorTask\n",
" config:\n",
" python: |\n",
" from lsst.ts.wep.task.pairTask import GroupPairer\n",
" config.pairer.retarget(GroupPairer)\n",
" donutStampSize: 200\n",
" initialCutoutPadding: 40\n",
" calcZernikesTask:\n",
" class: lsst.ts.wep.task.calcZernikesTask.CalcZernikesTask\n",
" config:\n",
" estimateZernikes.maxNollIndex: 28\n",
" estimateZernikes.saveHistory: False\n",
" donutStampSelector.maxSelect: 5\n",
" python: |\n",
" from lsst.ts.wep.task import EstimateZernikesDanishTask\n",
" config.estimateZernikes.retarget(EstimateZernikesDanishTask)\n",
" isr:\n",
" class: lsst.ip.isr.IsrTaskLSST\n",
" config:\n",
" # Although we don't have to apply the amp offset corrections, we do want\n",
" # to compute them for analyzeAmpOffsetMetadata to report on as metrics.\n",
" doAmpOffset: true\n",
" ampOffset.doApplyAmpOffset: false\n",
" # Turn off slow steps in ISR\n",
" doBrighterFatter: false\n",
" # Mask saturated pixels,\n",
" # but turn off quadratic crosstalk because it's currently broken\n",
" doSaturation: True\n",
" crosstalk.doQuadraticCrosstalkCorrection: False\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from lsst.daf.butler import Butler"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Change this path to appropriate butler when on-sky images arrive\n",
"path_to_aos_butler = '/repo/embargo'\n",
"butler = Butler(path_to_aos_butler)\n",
"\n",
"tie_collection = \"u/crenshaw/WET-007_tie_20241027_3\"\n",
"danish_collection = \"u/crenshaw/WET-007_danish_20241027_3\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot donut stamps\n",
"\n",
"First we will plot the donut stamps to see what they look like.\n",
"\n",
"Above each stamp I print the field angle of the donut so we can identify whether the TIE and Danish donut pairs are the same.\n",
"\n",
"#### Notes & Questions:\n",
"- It looks like we need to expand the default stamp size"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def plot_pairs(ID, collection, axes, title):\n",
" # Determine if each pair was selected \n",
" quality = butler.get(\"donutQualityTable\", dataId=ID, collections=collection)\n",
" intraSelected = quality[quality[\"DEFOCAL_TYPE\"] == \"intra\"][\"FINAL_SELECT\"]\n",
" extraSelected = quality[quality[\"DEFOCAL_TYPE\"] == \"extra\"][\"FINAL_SELECT\"]\n",
"\n",
" intraStamps = butler.get(\"donutStampsIntra\", dataId=ID, collections=collection)\n",
" intraStamps = [stamp for stamp, selected in zip(intraStamps, intraSelected) if selected]\n",
" extraStamps = butler.get(\"donutStampsExtra\", dataId=ID, collections=collection)\n",
" extraStamps = [stamp for stamp, selected in zip(extraStamps, extraSelected) if selected]\n",
" \n",
" usedList = butler.get(\"zernikes\", dataId=ID, collections=collection)[1:][\"used\"]\n",
" \n",
" for i, (used, intra, extra) in enumerate(zip(usedList, intraStamps, extraStamps)):\n",
" axes[0, i].imshow(intra.wep_im.image, origin=\"lower\")\n",
" angle = intra.wep_im.fieldAngle\n",
" axes[0, i].set_title(f\"({angle[0]:.2f}, {angle[1]:.2f})\", fontsize=8)\n",
" \n",
" axes[1, i].imshow(extra.wep_im.image, origin=\"lower\")\n",
" angle = extra.wep_im.fieldAngle\n",
" axes[1, i].set_title(f\"({angle[0]:.2f}, {angle[1]:.2f})\", fontsize=8)\n",
" if not used:\n",
" npix = intra.wep_im.image.shape[0] - 1\n",
" axes[0, i].plot([0, npix], [0, npix], c=\"r\", lw=1)\n",
" axes[0, i].plot([0, npix], [npix, 0], c=\"r\", lw=1)\n",
" axes[1, i].plot([0, npix], [0, npix], c=\"r\", lw=1)\n",
" axes[1, i].plot([0, npix], [npix, 0], c=\"r\", lw=1)\n",
" for ax in axes[:, len(usedList):].flatten():\n",
" ax.set_axis_off()\n",
" for ax in axes.flatten():\n",
" ax.set(xticks=[], yticks=[])\n",
" axes[0, 0].set_ylabel(f\"{title}\\nIntra\")\n",
" axes[1, 0].set_ylabel(f\"{title}\\nExtra\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tie_ids = list(butler.registry.queryDataIds(('exposure', 'visit', 'detector'), collections=tie_collection, datasets='donutStampsIntra'))\n",
"dan_ids = list(butler.registry.queryDataIds(('exposure', 'visit', 'detector'), collections=danish_collection, datasets='donutStampsIntra'))\n",
"tie_ids = sorted(set(tie_ids))\n",
"dan_ids = sorted(set(dan_ids))\n",
"\n",
"for t_id, d_id in zip(tie_ids, dan_ids):\n",
" fig, axes = plt.subplots(4, 5, dpi=120, constrained_layout=True)\n",
" plot_pairs(t_id, tie_collection, axes[:2], \"TIE\")\n",
" plot_pairs(d_id, danish_collection, axes[2:4], \"Danish\")\n",
" fig.suptitle(f\"Detector {t_id['detector']}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load Zernike Estimates\n",
"\n",
"Below we load the Zernikes from the two pipelines above and compare the outputs from the TIE and Danish"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Load zernike estimates\n",
"tie_ids = list(butler.registry.queryDataIds(('exposure', 'visit', 'detector'), collections=tie_collection, datasets='zernikeEstimateAvg'))\n",
"tie_zks = []\n",
"for data_id in tie_ids:\n",
" table = butler.get('zernikes', dataId=data_id, collections=tie_collection)\n",
" tie_zks.append(np.array([table[table[\"label\"] == \"average\"][f\"Z{i}\"][0].value for i in range(4, 29)]))\n",
"tie_zks = np.array(tie_zks)\n",
"\n",
"dan_ids = list(butler.registry.queryDataIds(('exposure', 'visit', 'detector'), collections=danish_collection, datasets='zernikeEstimateAvg'))\n",
"dan_zks = []\n",
"for data_id in tie_ids:\n",
" table = butler.get('zernikes', dataId=data_id, collections=danish_collection)\n",
" dan_zks.append(np.array([table[table[\"label\"] == \"average\"][f\"Z{i}\"][0].value for i in range(4, 29)]))\n",
"dan_zks = np.array(dan_zks)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, 1, dpi=150)\n",
"\n",
"ax.errorbar(np.arange(4, 29), np.mean(tie_zks, axis=0), yerr=np.std(tie_zks, axis=0), ls=\"--\", c=\"C2\", label=\"TIE\", capsize=4)\n",
"ax.errorbar(np.arange(4, 29), np.mean(dan_zks, axis=0), yerr=np.std(dan_zks, axis=0), ls=\"--\", c=\"cornflowerblue\", alpha=1, label=\"Danish\", capsize=4)\n",
"ax.scatter([4, 5, 6, 7, 8], [-1573.3, 819.4, 1290.6, -1920.3, 671.3], c=\"r\", marker=\"x\", label=\"Josh\", zorder=10, s=70)\n",
"\n",
"ax.legend()\n",
"\n",
"ax.set(\n",
" xlabel=\"Noll index\",\n",
" ylabel=\"Amplitude [nm]\",\n",
" xticks=np.arange(4, 29, 2),\n",
" xlim=(3.5, 28.5),\n",
")\n",
"\n",
"# inset Axes....\n",
"x1, x2, y1, y2 = 3.5, 8.75, -2350, +2300 # subregion of the original image\n",
"axins = ax.inset_axes(\n",
" [0.5, 0.075, 0.4, 0.35],\n",
" xlim=(x1, x2), ylim=(y1, y2), xticks=np.arange(4, 9), yticks=[-2e3, -1e3, 0, 1e3, 2e3], ylabel=\"nm\")\n",
"\n",
"axins.errorbar(np.arange(4, 29), np.mean(tie_zks, axis=0), yerr=np.std(tie_zks, axis=0), ls=\"--\", c=\"C2\", label=\"TIE\", capsize=4)\n",
"axins.errorbar(np.arange(4, 29), np.mean(dan_zks, axis=0), yerr=np.std(dan_zks, axis=0), ls=\"--\", c=\"cornflowerblue\", alpha=1, label=\"Danish\", capsize=4)\n",
"axins.scatter([4, 5, 6, 7, 8], [-1573.3, 819.4, 1290.6, -1920.3, 671.3], c=\"r\", marker=\"x\", label=\"Josh\", zorder=10, s=70)\n",
"\n",
"\n",
"ax.indicate_inset_zoom(axins, edgecolor=\"gray\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The 4 points labeled \"Josh\" were created manually by Josh Meyers, who used Danish to jointly fit a single intra- & extra-focal donut.\n",
"Points above are means of the combined Zernikes from each of the 9 detectors.\n",
"The uncertainties are the standard deviations of these combined Zernikes.\n",
"\n",
"The TIE and Danish mostly gree within their purported uncertainties.\n",
"Astigmatism (indices 5 and 6) seem the noisiest, followed by coma (7 and 8) and focus (4).\n",
"\n",
"The largest disagreements are in Zernikes 13 (oblique secondary astigmatism) and 20 (oblique pentafoil).\n",
"The amplitude of each is small though, so I don't think this is an issue."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, 1, dpi=150)\n",
"\n",
"ax.plot(np.arange(4, 29), np.abs(np.std(tie_zks, axis=0) / np.mean(tie_zks, axis=0)), c=\"C2\", label=\"TIE\")\n",
"ax.plot(np.arange(4, 29), np.abs(np.std(dan_zks, axis=0) / np.mean(dan_zks, axis=0)), c=\"cornflowerblue\", label=\"Danish\")\n",
"ax.legend()\n",
"\n",
"ax.set(\n",
" xlabel=\"Noll index\",\n",
" ylabel=\"Relative uncertainty\",\n",
" xticks=np.arange(4, 29, 2),\n",
" xlim=(3.5, 28.5),\n",
" #yscale=\"log\"\n",
" ylim=(0, 1),\n",
")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Relative uncertainties `(standard deviation / mean)` are relatively consistent at the ~30% level.\n",
"Note that the relative uncertainties of some high-order indices blow up because the estimated amplitude of these Zernikes is very small.\n",
"This is not a problem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Discussion\n",
"\n",
"Everything looks good in this analysis.\n",
"Donut selection works reasonably well, Zernike estimation is consistent between the three methods plotted.\n",
"\n",
"Note that some very-close blends are making it through the pipeline.\n",
"It seems this is happening because `GenerateDonutDirectDetectTask` thinks these are a single source.\n",
"Josh Meyers suggested to fix this by checking for a \"donut hole\" by measuring flux in different annuli."
]
},
{
"cell_type": "code",
"execution_count": null,
"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.11.9"
}
},
"nbformat": 4,
"nbformat_minor": 4
}