{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reduction\n",
    "\n",
    "This notebook will demonstrate basic OPTICAM data reduction using `phoptic`.\n",
    "\n",
    "## Performance Tip\n",
    "\n",
    "Before starting, it's worth noting that `phoptic` uses the `multiprocessing` Python module to parallelize data reduction across multiple CPU cores. For the best scaling, it's therefore a good idea to disable automatic parallelization in any underlying libraries (like [`numpy`](https://numpy.org)) using some combination of the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:20:52.721743Z",
     "iopub.status.busy": "2026-01-23T11:20:52.721659Z",
     "iopub.status.idle": "2026-01-23T11:20:52.723353Z",
     "shell.execute_reply": "2026-01-23T11:20:52.723175Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# limit underlying math libraries to a single thread for better multiprocessing performance\n",
    "os.environ['OMP_NUM_THREADS'] = '1'  # OpenMP\n",
    "# os.environ['OPENBLAS_NUM_THREADS'] = '1'  # OpenBLAS\n",
    "# os.environ['MKL_NUM_THREADS'] = '1'  # Intel Math Kernel Library\n",
    "# os.environ['VECLIB_MAXIMUM_THREADS'] = '1'  # Apple Accelerate vector library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The relevant libraries will depend on your build and the specific system you're using; in my case, I only need to include `os.environ['OMP_NUM_THREADS'] = '1'`.\n",
    "\n",
    "**Warning**: if you do not disable automatic parallelization as shown above, [cosmic ray removal](#defining-a-reducer) performance will be extremely degraded and/or may *hang*.\n",
    "\n",
    "## Generating Data\n",
    "\n",
    "Before we can begin reducing data, we need some data to reduce. `phoptic` provides a routine for generating some dummy observations via the `generate_observations()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:20:52.724102Z",
     "iopub.status.busy": "2026-01-23T11:20:52.724032Z",
     "iopub.status.idle": "2026-01-23T11:21:01.172280Z",
     "shell.execute_reply": "2026-01-23T11:21:01.172050Z"
    }
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import phoptic\n",
    "\n",
    "\n",
    "out_dir = Path('out')\n",
    "\n",
    "phoptic.generate_observations(\n",
    "    out_directory=out_dir / 'data',  # path to the directory where the generated data will be saved\n",
    "    circular_aperture=False,  # disable circular aperture shadow\n",
    "    n_images=50,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These dummy observations come in three filters: $g$, $r$, and $i$, and will be used in many of the guided tutorials. We can also see that there is a variable source at (131, 115), which we will use as our source of interest. We're also given the fractional RMS amplitude of the variability, its frequency, the phase lags between the different filters. These values can be used as a sanity check.\n",
    "\n",
    "There are a number of parameters that can be tweaked when calling `generate_observations()`, such as the number of images, the binning of the images, and whether a circular aperture shadow should be applied to the images. For this example, I've only generated a small number of images for simplicity.\n",
    "\n",
    "## Defining an instrument\n",
    "\n",
    "In `phoptic`, `Instrument`s are used to interface with data products. An `Instrument` defines the keywords and necessary routines required for `phoptic` to interpret the headers of images from a particular instrument. Since `phoptic` was developed primarily for use with OPTICAM-MX, the data generated by the `generate_observations()` function are compatible with the `OPTICAM_MX` instrument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:01.173207Z",
     "iopub.status.busy": "2026-01-23T11:21:01.173052Z",
     "iopub.status.idle": "2026-01-23T11:21:01.174634Z",
     "shell.execute_reply": "2026-01-23T11:21:01.174448Z"
    }
   },
   "outputs": [],
   "source": [
    "instrument = phoptic.OPTICAM_MX()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of `phoptic`'s reduction routines require an `Instrument` instance, and we have tried to make it straightforward for users to define custom `Instrument` instances for their instruments (see the [instruments tutorial](instruments.ipynb) for more details).\n",
    "\n",
    "## Pre-reduction checks\n",
    "\n",
    "Before we reduce any data, it's a good idea to run `phoptic.scan_data()`, especially if you're unfamiliar with the data. This scans the image headers and summarises the relevant information for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:01.175333Z",
     "iopub.status.busy": "2026-01-23T11:21:01.175257Z",
     "iopub.status.idle": "2026-01-23T11:21:01.942636Z",
     "shell.execute_reply": "2026-01-23T11:21:01.942375Z"
    }
   },
   "outputs": [],
   "source": [
    "phoptic.scan_data(\n",
    "    out_directory=out_dir / 'reduced',\n",
    "    data_directory=out_dir / 'data',\n",
    "    instrument=instrument,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the images use the 4x4 binning mode, we have $g$-, $r$-, and $i$-band images, and there are 20 images for each filter. These checks are useful to make sure that all images use the same binning, for example, which will result in an error. If an error is raised, `phoptic` will explain what caused the error such that it should be straightforward to resolve. A design choice of `phoptic` is that it will only ever read from the data directory and never *write* to it. As such, errors will need to be resolved manually. Since there were no issues with these data, however, we can proceed with reduction by defining a `Reducer` instance.\n",
    "\n",
    "## Defining a `Reducer`\n",
    "\n",
    "When defining a `Reducer`, there are two required parameters: `data_directory` and `out_directory`. `data_directory` should be a path to a directory that contains a series of FITS files (or a single multi-extension FITS file). `out_directory` determines the directory to which any output files are written. If `out_directory` does not exist, `phoptic` will attempt to create it. We can also apply flat-field corrections to our images by passing a `FlatFieldCorrector`, though we will not do this here (see the [corrections tutorial](applying_corrections.ipynb) for details on applying flat-field corrections). Another useful parameter is `remove_cosmic_rays`, which can be used to clip cosmic rays from your images. If `remove_cosmic_rays=True`, cosmic rays are clipped using [astroscrappy](https://astroscrappy.readthedocs.io/en/latest/index.html)'s [L.A.Cosmic](https://ui.adsabs.harvard.edu/abs/2001PASP..113.1420V/exportcitation) implementation. By default, `remove_cosmic_rays` is set to `False`.\n",
    "\n",
    "Let's initialise a `Reducer` using its default paramters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:01.943720Z",
     "iopub.status.busy": "2026-01-23T11:21:01.943621Z",
     "iopub.status.idle": "2026-01-23T11:21:04.418360Z",
     "shell.execute_reply": "2026-01-23T11:21:04.417996Z"
    }
   },
   "outputs": [],
   "source": [
    "reducer = phoptic.Reducer(\n",
    "    data_directory=out_dir / 'data',  # path to the data\n",
    "    out_directory=out_dir / 'reduced',  # path to where the output will be saved\n",
    "    instrument=instrument,  # the instrument (defaults to OPTICAM_MX())\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After creating a `Reducer` instance, `phoptic.scan_data()` is ran to ensure there are no issues with the data; as mentioned above, this may result in errors that will need to be resolved.\n",
    "\n",
    "We can also see that our chosen instrument, `OPTICAM_MX` passed some internal checks. See the [custom instruments](instruments.ipynb) tutorial for guidance on how to create custom instruments for use with `phoptic`, and how to fix potential errors that might arise when doing so.\n",
    "\n",
    "## Create Source Catalogs\n",
    "\n",
    "The next step is to create source catalogs for each camera. Creating source catalogs requires aligning each camera's images to track sources over time and ensure consistent labelling.\n",
    "\n",
    "The way in which images are aligned can be customised. By default, images will be aligned using `transform_type='affine'` which uses the [astroalign](https://astroalign.quatrope.org/en/latest/) package. Alternatively, passing `transform_type='translation'` will compute simple (x, y) translations between images. The number of reference sources can also be changed by passing the desired value to `n_alignment_sources`. When using `transform_type='translation'`, `n_alignment_sources` must be $\\geq 1$ while `transform_type='affine'` requires `n_alignment_sources` be $\\geq 3$.\n",
    "\n",
    "To create the source catalogs, we must call the `create_catalogs()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:04.419351Z",
     "iopub.status.busy": "2026-01-23T11:21:04.419248Z",
     "iopub.status.idle": "2026-01-23T11:21:08.264811Z",
     "shell.execute_reply": "2026-01-23T11:21:08.264561Z"
    }
   },
   "outputs": [],
   "source": [
    "reducer.create_catalogs(overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that six sources have been identified in each of the three sets of images. In this example, the source labelling is consistent across the catalogs because there are no field-of-view or pixel-scale differences between the simulated cameras. In practise, the source labelling will not usually be consistent across the catalogs, and so care must be taken when performing differential photometry that the same sources are being used for each filter (more on this later).\n",
    "\n",
    "It's also worth discussing `phoptic`'s labelling convention. For compatibility with multi-camera instruments, like OPTICAM-MX, `phoptic` groups images according to both the camera and filter used to take the image. In the above catalogue, you can see the catalogues are titled \"1:g\", \"2:r\", and \"3:i\". The value before the colon represents the camera, which may simply be the name of the instrument if only has one camera, while the value after the colon represents the filter. In the case of OPTICAM-MX, filters are unique to each camera, though this may not be the case for all instruments. See the [custom instruments tutorial](instruments.ipynb) for more details.\n",
    "\n",
    "In addition to the source catalogues, `create_catalogs()` also produces a plot of the observation's systematics (e.g., background, seeing, and airmass) for sanity checking.\n",
    "\n",
    "## Performing Photometry\n",
    "\n",
    "With our catalogs defined, we can now perform photometry. Performing photometry in `phoptic` requires a photometer object that inherits from `phoptic.photometers.BasePhotometer`. Currently, `phoptic` provides two photometers: `AperturePhotometer`, for performing simple aperture photometry, and `OptimalPhotometer`, for performing optimal photometry (as described in [Naylor 1998](https://ui.adsabs.harvard.edu/abs/1998MNRAS.296..339N/abstract)). Both photometers can be configured to perform \"forced photometry\", and can use either the `Reducer`'s 2D background estimator or estimate the local background around each source using an annulus. In this example, I'll show a couple of different photometry configurations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:08.265750Z",
     "iopub.status.busy": "2026-01-23T11:21:08.265655Z",
     "iopub.status.idle": "2026-01-23T11:21:08.267545Z",
     "shell.execute_reply": "2026-01-23T11:21:08.267299Z"
    }
   },
   "outputs": [],
   "source": [
    "# aperture photometer with local background estimations\n",
    "default_annulus_photometer = phoptic.AperturePhotometer(\n",
    "    local_background_estimator=phoptic.DefaultLocalBackground(),  # use the default local background estimator\n",
    ")\n",
    "\n",
    "# optimal photometer\n",
    "# implements the method described in Naylor 1998, MNRAS, 296, 339-346\n",
    "optimal_photometer = phoptic.OptimalPhotometer()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once a photometer has been defined, it can be passed to the `photometry()` method of `Catalog` to compute the raw light curves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:08.268326Z",
     "iopub.status.busy": "2026-01-23T11:21:08.268247Z",
     "iopub.status.idle": "2026-01-23T11:21:10.321070Z",
     "shell.execute_reply": "2026-01-23T11:21:10.320810Z"
    }
   },
   "outputs": [],
   "source": [
    "reducer.photometry(default_annulus_photometer)  # using the aperture photometer with local background estimations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you perform photometry on a catalog, the directory to which the light curves are saved will depend on how the photometer is configured:\n",
    "1. If the photometer is defined with `forced=True`, the directory will have a \"forced\" prefix.\n",
    "2. If a local background estimator is passed to the photometer's `local_background_estimator` parameter, then the directory will also have an \"annulus\" suffix.\n",
    "\n",
    "In the above case, we used the `AperturePhotometer` with a local background estimator and so we can see that the raw light curves have been saved to a `lcs/aperture_annulus` directory. Let's compare this to the subdirectory for light curves produced by our `optimal_photometer`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:10.322124Z",
     "iopub.status.busy": "2026-01-23T11:21:10.322014Z",
     "iopub.status.idle": "2026-01-23T11:21:15.834885Z",
     "shell.execute_reply": "2026-01-23T11:21:15.834620Z"
    }
   },
   "outputs": [],
   "source": [
    "reducer.photometry(optimal_photometer)  # using the optimal photometer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case there is no prefix nor suffix, and only the name of the photometer class is included.\n",
    "\n",
    "## Computing Relative Light Curves\n",
    "\n",
    "We now have some raw light curves using a couple of different photometry configurations. However, raw light curves contain a lot of atmospheric and systematic variability, and so we often want to compute the relative light curve between our source of interest and some comparison sources to reduce this atmospheric/systematic variability.\n",
    "\n",
    "In this example, let's say that Source 2 is our target of interest (since we know that it's varying), and use Source 1 for comparison. In this example, the choice of comparison source(s) is arbitraray (since there are no atmospheric/systematic variations in our simulated data). In practise, however, choosing suitable comparison sources is vital for obtaining quality light curves, though a discussion of this is beyond the scope of this tutorial.\n",
    "\n",
    "Let's now produce a relative light curve for Source 2 using the \"aperture_annulus\" light curves created by `default_annulus_photometer`. First, however, we need to initialise a `DifferentialPhotometer` object. When initialising a `DifferentialPhotometer` object, we need to pass the directory path to the reduced data created by `Reducer`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:15.835958Z",
     "iopub.status.busy": "2026-01-23T11:21:15.835726Z",
     "iopub.status.idle": "2026-01-23T11:21:16.122045Z",
     "shell.execute_reply": "2026-01-23T11:21:16.121810Z"
    }
   },
   "outputs": [],
   "source": [
    "dphot = phoptic.DifferentialPhotometer(\n",
    "    out_directory=out_dir / 'reduced',  # same as the catalog's out_directory\n",
    "    show_plots=True,  # show plots (useful for diagnosis and debugging)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When initialising a `DifferentialPhotometer` object, the source catalogs are output for convenience (unless `show_plots=False`). We can now create our relative light curves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:16.123724Z",
     "iopub.status.busy": "2026-01-23T11:21:16.123596Z",
     "iopub.status.idle": "2026-01-23T11:21:17.885707Z",
     "shell.execute_reply": "2026-01-23T11:21:17.885499Z"
    }
   },
   "outputs": [],
   "source": [
    "target = 2  # source of interest\n",
    "comparisons = [1]  # comparison sources\n",
    "\n",
    "annulus_analyzer = dphot.get_relative_light_curve(\n",
    "    '1:g',  # filter for which to compute the relative light curve\n",
    "    target,  # source of interest\n",
    "    comparisons,  # comparison sources\n",
    "    phot_label='aperture_annulus',  # label for the photometry results\n",
    "    prefix='test',  # prefix for the output files (e.g., the name of the target source)\n",
    "    match_other_cameras=True,  # match sources across cameras\n",
    "    show_diagnostics=True,  # show diagnostic plots (useful for diagnosis and debugging)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the relative light curves are generated, a number of diagnostic plots are rendered. These diagnostic plots can be suppressed by passing `show_diagnostics=False`. Regardless of the value of `show_diagnostics`, these diagnostic plots are saved to `out_directory/relative_light_curves/diag`. The resulting relative light curves are also shown for convenience. After relative light curves have been computed, an `Analyzer` object is returned (see the [timing methods](timing_methods.ipynb) tutorial for more details on this).\n",
    "\n",
    "Above, we used the `get_relative_light_curve()` method to automatically match sources across filters by setting `match_other_cameras=True`. However, this can misidentify sources, so care should be taken to check the correct sources are identified. In this case, we can see that the identified sources are correct, and so we don't have to manually create relative light curves for each filter. For the optimal light curves, however, I will show how to merge `Analyzer` instances using the `join()` method should you ever have to match the sources manually. First, let's define the individual `Analyzer` instances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:17.886684Z",
     "iopub.status.busy": "2026-01-23T11:21:17.886595Z",
     "iopub.status.idle": "2026-01-23T11:21:19.034045Z",
     "shell.execute_reply": "2026-01-23T11:21:19.033835Z"
    }
   },
   "outputs": [],
   "source": [
    "g_band_optimal_analyzer = dphot.get_relative_light_curve(\n",
    "    '1:g',\n",
    "    target,\n",
    "    comparisons,\n",
    "    phot_label='optimal',\n",
    "    prefix='test',\n",
    "    match_other_cameras=False,\n",
    "    show_diagnostics=False,\n",
    "    )\n",
    "\n",
    "r_band_optimal_analyzer = dphot.get_relative_light_curve(\n",
    "    '2:r',\n",
    "    target,\n",
    "    comparisons,\n",
    "    phot_label='optimal',\n",
    "    prefix='test',\n",
    "    match_other_cameras=False,\n",
    "    show_diagnostics=False,\n",
    "    )\n",
    "\n",
    "i_band_optimal_analyzer = dphot.get_relative_light_curve(\n",
    "    '3:i',\n",
    "    target,\n",
    "    comparisons,\n",
    "    phot_label='optimal',\n",
    "    prefix='test',\n",
    "    match_other_cameras=False,\n",
    "    show_diagnostics=False,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's combine them all into a single instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:19.035020Z",
     "iopub.status.busy": "2026-01-23T11:21:19.034842Z",
     "iopub.status.idle": "2026-01-23T11:21:19.078704Z",
     "shell.execute_reply": "2026-01-23T11:21:19.078497Z"
    }
   },
   "outputs": [],
   "source": [
    "# join g-band with r-band\n",
    "optimal_analyzer = g_band_optimal_analyzer.join(r_band_optimal_analyzer)\n",
    "\n",
    "# join g-band and r-band with i-band\n",
    "optimal_analyzer = optimal_analyzer.join(i_band_optimal_analyzer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check the `Analyzer` instances have merged properly, let's try plotting the light curves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:21:19.079449Z",
     "iopub.status.busy": "2026-01-23T11:21:19.079370Z",
     "iopub.status.idle": "2026-01-23T11:21:19.566579Z",
     "shell.execute_reply": "2026-01-23T11:21:19.566373Z"
    }
   },
   "outputs": [],
   "source": [
    "fig = optimal_analyzer.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, all three light curves have been combined successfully combined into a single `Analyzer` instance.\n",
    "\n",
    "That concludes the reduction tutorial for `phoptic`! If you're using `phoptic` to reduce data from OPTICAM-MX, I recommend checking out the [applying corrections](applying_corrections.ipynb) and [photometry](photometry.ipynb) tutorials next. To use `phoptic` with another instrument, check out the [custom instruments](instruments.ipynb) tutorial. To explore the quick-look timing analysis routines available in `phoptic`, see the [timing methods tutorial](timing_methods.ipynb)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "opticam7",
   "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.14.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}