<commit_before>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from rhlUtils import BBox, CCD, Image, imshow\n",
    "\n",
    "%matplotlib notebook \n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "#%matplotlib qt\n",
    "#%gui qt\n",
    "\n",
    "dataDir = \"/Users/rhl/TeX/Talks/DSFP/2018-01/Exercises/Detectors\"\n",
    "mag0 = 33   # Magnitude of an object with 1 detected photon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have provided three images containing stars, taken with 3 different CCDs, in \"stars_X.npz\" where X = 0, 1, 2.\n",
    "The true magnitudes of all the stars are also known, taken from a reliable reference catalogue.\n",
    "\n",
    "OK, I'll be honest; these are simulations.  The stars' centers are all at centres of pixels\n",
    "make your lives easier; there is no spatial structure in the PSF; and there is also no noise added of any kind.  You may buy me a drink later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To read the image data and the calibs say something like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.load(os.path.join(dataDir, \"stars_0.npz\"))\n",
    "image, calibs = data[\"image\"], data[\"calibs\"]\n",
    "\n",
    "image0 = image.copy()    # Keep a copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to make sure that we're all on the same page, here's code to display the image and stars (imshow is a simple utility imported from rhlUtils -- feel free to use plt.imshow if you'd rather)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(1)\n",
    "plt.clf()\n",
    "imshow(image, vmin=0, vmax=1000)\n",
    "plt.title(\"Data\")\n",
    "\n",
    "plt.plot(calibs[:, 0], calibs[:, 1], '+')   # calibs[:, 2] contains the object's magnitude (not flux)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time for you to do some work.  Write some code to estimate a PSF model by simply averaging all the objects, giving each equal weight.  You should add the option to use only some subset of the stars (e.g. the faintest 25%).\n",
    "\n",
    "Your model should simply be an image of the PSF, normalised to have a maximum value of 1.0\n",
    "\n",
    "You could use the calibs object to find all the stars, but in the real world there are stars in the data that are\n",
    "not in the catalogue so you'll have to write a simple object finder (i.e. don't use the calibs!).  It's sufficient to find pixels that are larger than any neighbours to the left, right, top, or bottom; no, this isn't quite how you'd do it in the real world but it's not far off.\n",
    "\n",
    "I told you that there were three data sets, taken in different places on the sky and with CCDs with different properties.  Choose one set (e.g. \"stars_0.npz\") to carry out the following activities, then go back and look at the other two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = image0.copy()\n",
    "\n",
    "#...\n",
    "\n",
    "plt.figure(2)\n",
    "plt.clf()\n",
    "imshow(psfIm, vmin=0, vmax=1.1)\n",
    "plt.title(\"PSF model\")\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, now use your PSF model to create an image of the residuals created by subtracting the scaled PSF from the stars.\n",
    "\n",
    "How does it look?  Do you see what's going on?  Remember, I said that the PSF isn't a function of position, but CCDs aren't perfect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = image0.copy()\n",
    "\n",
    "# ...\n",
    "\n",
    "plt.figure(3)\n",
    "plt.clf()\n",
    "imshow(image, vmin=image.min(), vmax=image.max()) # , vmin=0, vmax=100)\n",
    "plt.title(\"Residuals\")\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A powerful diagnostic is to measure the fluxes of the objects in different ways, and compare with the truth (i.e. the\n",
    "magnitudes given in the calibs data structure).\n",
    "\n",
    "For the purpose of this exercise you should use:\n",
    "- aperture fluxes (the sum of all the pixels in a square centred on the star is good enough)\n",
    "- psf fluxes (the flux in the best-fit PSF model, $\\phi$.  The amplitude is simply $A = \\frac{\\sum \\phi I}{\\sum \\phi^2}$, and the flux is $A \\sum \\phi$)\n",
    "\n",
    "Convert them to magnitudes using the known zero-point mag0 (the magnitude corresponding to one count), and make suitable plots.\n",
    "\n",
    "What do you think is going on?  Are there other measurements that you can make on the data to test your hypotheses?  Are there other observations that you'd like (I might be able to make them for you)."
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}<commit_msg>Add note about downloading the image data<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before beginning this exercise you must download some data files, which can be retrieved from here: \n",
    "\n",
    "https://northwestern.box.com/s/rsb5wjb7dztg2128flzgex4fmq19havu\n",
    "\n",
    "Be sure to move the corresponding files to the directory where you are running this notebook. Furthermore, you will need to provide the path to these data as the variable `dataDir`, below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from rhlUtils import BBox, CCD, Image, imshow\n",
    "\n",
    "%matplotlib notebook \n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "#%matplotlib qt\n",
    "#%gui qt\n",
    "\n",
    "dataDir = # complete\n",
    "mag0 = 33   # Magnitude of an object with 1 detected photon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have provided three images containing stars, taken with 3 different CCDs, in \"stars_X.npz\" where X = 0, 1, 2.\n",
    "The true magnitudes of all the stars are also known, taken from a reliable reference catalogue.\n",
    "\n",
    "OK, I'll be honest; these are simulations.  The stars' centers are all at centres of pixels\n",
    "make your lives easier; there is no spatial structure in the PSF; and there is also no noise added of any kind.  You may buy me a drink later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To read the image data and the calibs say something like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = np.load(os.path.join(dataDir, \"stars_0.npz\"))\n",
    "image, calibs = data[\"image\"], data[\"calibs\"]\n",
    "\n",
    "image0 = image.copy()    # Keep a copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to make sure that we're all on the same page, here's code to display the image and stars (imshow is a simple utility imported from rhlUtils -- feel free to use plt.imshow if you'd rather)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.figure(1)\n",
    "plt.clf()\n",
    "imshow(image, vmin=0, vmax=1000)\n",
    "plt.title(\"Data\")\n",
    "\n",
    "plt.plot(calibs[:, 0], calibs[:, 1], '+')   # calibs[:, 2] contains the object's magnitude (not flux)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time for you to do some work.  Write some code to estimate a PSF model by simply averaging all the objects, giving each equal weight.  You should add the option to use only some subset of the stars (e.g. the faintest 25%).\n",
    "\n",
    "Your model should simply be an image of the PSF, normalised to have a maximum value of 1.0\n",
    "\n",
    "You could use the calibs object to find all the stars, but in the real world there are stars in the data that are\n",
    "not in the catalogue so you'll have to write a simple object finder (i.e. don't use the calibs!).  It's sufficient to find pixels that are larger than any neighbours to the left, right, top, or bottom; no, this isn't quite how you'd do it in the real world but it's not far off.\n",
    "\n",
    "I told you that there were three data sets, taken in different places on the sky and with CCDs with different properties.  Choose one set (e.g. \"stars_0.npz\") to carry out the following activities, then go back and look at the other two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "image = image0.copy()\n",
    "\n",
    "#...\n",
    "\n",
    "plt.figure(2)\n",
    "plt.clf()\n",
    "imshow(psfIm, vmin=0, vmax=1.1)\n",
    "plt.title(\"PSF model\")\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, now use your PSF model to create an image of the residuals created by subtracting the scaled PSF from the stars.\n",
    "\n",
    "How does it look?  Do you see what's going on?  Remember, I said that the PSF isn't a function of position, but CCDs aren't perfect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "image = image0.copy()\n",
    "\n",
    "# ...\n",
    "\n",
    "plt.figure(3)\n",
    "plt.clf()\n",
    "imshow(image, vmin=image.min(), vmax=image.max()) # , vmin=0, vmax=100)\n",
    "plt.title(\"Residuals\")\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A powerful diagnostic is to measure the fluxes of the objects in different ways, and compare with the truth (i.e. the\n",
    "magnitudes given in the calibs data structure).\n",
    "\n",
    "For the purpose of this exercise you should use:\n",
    "- aperture fluxes (the sum of all the pixels in a square centred on the star is good enough)\n",
    "- psf fluxes (the flux in the best-fit PSF model, $\\phi$.  The amplitude is simply $A = \\frac{\\sum \\phi I}{\\sum \\phi^2}$, and the flux is $A \\sum \\phi$)\n",
    "\n",
    "Convert them to magnitudes using the known zero-point mag0 (the magnitude corresponding to one count), and make suitable plots.\n",
    "\n",
    "What do you think is going on?  Are there other measurements that you can make on the data to test your hypotheses?  Are there other observations that you'd like (I might be able to make them for you)."
   ]
  }
 ],
 "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
{'fraction_non_alphanumeric': 0.1639101456430511, 'fraction_numerical': 0.008804410433637786, 'mean_word_length': 2.7862928348909657, 'pattern_counts': {'":': 130, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '641712', 'n_tokens_mistral': 4353, 'n_tokens_neox': 4152, 'n_words': 1460}