Exoplanet_TSO_Pipeline_Median.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext autoreload\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from astroML.plotting          import hist\n",
    "from astropy.io                import fits\n",
    "from astropy.modeling          import models, fitting\n",
    "from datetime                  import datetime\n",
    "from image_registration        import cross_correlation_shifts\n",
    "from glob                      import glob\n",
    "from functools                 import partial\n",
    "from matplotlib.ticker         import MaxNLocator\n",
    "from matplotlib                import style\n",
    "from os                        import listdir\n",
    "# from least_asymmetry.asym      import actr, moments, fitgaussian\n",
    "from multiprocessing           import cpu_count, Pool\n",
    "from numpy                     import min as npmin, max as npmax, zeros, arange, sum, float, isnan, hstack\n",
    "from numpy                     import int32 as npint, round as npround, nansum as sum, nanstd as std\n",
    "from os                        import environ, path, mkdir\n",
    "from pandas                    import DataFrame, read_csv, read_pickle, scatter_matrix\n",
    "from photutils                 import CircularAperture, CircularAnnulus, aperture_photometry, findstars\n",
    "from pylab                     import ion, gcf, sort, linspace, indices, median, mean, std, empty, figure, transpose, ceil\n",
    "from pylab                     import concatenate, pi, sqrt, ones, diag, inf, rcParams, isnan, isfinite, array, nanmax\n",
    "from pylab                     import figure, plot, imshow, scatter, legend\n",
    "from seaborn                   import *\n",
    "from scipy.special             import erf\n",
    "from scipy                     import stats\n",
    "from sklearn.cluster           import DBSCAN\n",
    "from sklearn.externals         import joblib\n",
    "from sklearn.preprocessing     import StandardScaler\n",
    "from socket                    import gethostname\n",
    "from statsmodels.robust        import scale\n",
    "from statsmodels.nonparametric import kde\n",
    "from sys                       import exit\n",
    "from time                      import time, localtime\n",
    "from tqdm                      import tqdm_notebook\n",
    "\n",
    "from numpy                     import zeros, nanmedian as median, nanmean as mean, nan\n",
    "from sys                       import exit\n",
    "from sklearn.externals         import joblib\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "startFull = time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Master Class for Exoplanet Time Series Observation Photometry**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from wanderer import wanderer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clipOutlier2D(arr2D, nSig=10):\n",
    "    arr2D     = arr2D.copy()\n",
    "    medArr2D  = median(arr2D,axis=0)\n",
    "    sclArr2D  = np.sqrt(((scale.mad(arr2D)**2.).sum()))\n",
    "    outliers  = abs(arr2D - medArr2D) >  nSig*sclArr2D\n",
    "    inliers   = abs(arr2D - medArr2D) <= nSig*sclArr2D\n",
    "    arr2D[outliers] = median(arr2D[inliers],axis=0)\n",
    "    return arr2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rcParams['figure.dpi'] = 150\n",
    "rcParams['image.interpolation'] = 'None'\n",
    "rcParams['image.cmap']          = 'Blues_r'\n",
    "rcParams['axes.grid']           = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, Spitzer data is expected to be store in the directory structure:\n",
    "\n",
    "`$HOME/PLANET_DIRECTORY/PLANETNAME/data/raw/AORDIR/CHANNEL/bcd/`\n",
    "\n",
    "EXAMPLE:\n",
    "\n",
    "1. On a Linux machine\n",
    "2. With user `tempuser`,\n",
    "3. And all Spitzer data is store in `Research/Planets`\n",
    "4. The planet named `Happy-5b`\n",
    "5. Observed during AOR r11235813\n",
    "6. In CH2 (4.5 microns)\n",
    "\n",
    "The `loadfitsdir` should read as: `/home/tempuser/Research/Planets/HATP26/data/raw/r11235813/ch2/bcd/`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os import environ\n",
    "\n",
    "planetName      = 'HATP26'\n",
    "planetDirectory = '/Research/Planets/'\n",
    "\n",
    "channel = 'ch2'\n",
    "# channel = 'ch2/'\n",
    "\n",
    "dataSub = 'bcd/'\n",
    "\n",
    "dataDir     = environ['HOME'] + planetDirectory + planetName + '/data/raw/' + channel + '/big/'\n",
    "\n",
    "AORs = []\n",
    "for dirNow in glob(dataDir + '/*'):\n",
    "    AORs.append(dirNow.split('/')[-1])\n",
    "\n",
    "fileExt = '*bcd.fits'\n",
    "uncsExt = '*bunc.fits'\n",
    "\n",
    "print(dataDir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(AORs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iAOR        = 0 # Change this to 0,1,2,3,4,5,6,7 to \"cycle\" over the 8 AORs in the base directory /home/tempuser/Research/Planets/HATP26/data/raw/\n",
    "AORNow      = AORs[iAOR] # This will flag an error if no AORs were found in the directory `dataDir`\n",
    "loadfitsdir = dataDir + AORNow + '/' + channel + '/' + dataSub\n",
    "print(loadfitsdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nCores = cpu_count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fitsFilenames = glob(loadfitsdir + fileExt);print(len(fitsFilenames))\n",
    "uncsFilenames = glob(loadfitsdir + uncsExt);print(len(uncsFilenames))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "header_test = fits.getheader(fitsFilenames[0])\n",
    "print('AORLABEL:\\t{}\\nNum Fits Files:\\t{}\\nNum Unc Files:\\t{}'.format\\\n",
    "          (header_test['AORLABEL'], len(fitsFilenames), len(uncsFilenames)))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "fitsFilenames"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "uncsFilenames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Wanderer Class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Necessary Constants Spitzer\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppm             = 1e6\n",
    "y,x             = 0,1\n",
    "\n",
    "yguess, xguess  = 15., 15.   # Specific to Spitzer circa 2010 and beyond\n",
    "filetype        = 'bcd.fits' # Specific to Spitzer Basic Calibrated Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start a New Instance with Median for the Metric\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "method = 'median'\n",
    "\n",
    "print('Initialize an instance of `wanderer` as `example_wanderer_median`\\n')\n",
    "example_wanderer_median = wanderer(fitsFileDir=loadfitsdir, filetype=filetype, telescope='Spitzer', \n",
    "                                            yguess=yguess, xguess=xguess, method=method, nCores=nCores)\n",
    "\n",
    "example_wanderer_median.AOR        = AORNow\n",
    "example_wanderer_median.planetName = planetName\n",
    "example_wanderer_median.channel    = channel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Load Data From Fits Files in ' + loadfitsdir + '\\n')\n",
    "example_wanderer_median.spitzer_load_fits_file(outputUnits='electrons')#(outputUnits='muJ_per_Pixel')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Double check for NaNs**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_wanderer_median.imageCube[np.where(isnan(example_wanderer_median.imageCube))] = \\\n",
    "                                                    np.nanmedian(example_wanderer_median.imageCube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Identifier Strong Outliers**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Find, flag, and NaN the \"Bad Pixels\" Outliers' + '\\n')\n",
    "example_wanderer_median.find_bad_pixels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flux Weighted Centroiding -- just to say we tried it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Fit for All Centers: Flux Weighted, Gaussian Fitting, Gaussian Moments, Least Asymmetry' + '\\n')\n",
    "example_wanderer_median.fit_flux_weighted_centering()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaussian centroid fitting -- the most widely used version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "start = time()\n",
    "example_wanderer_median.mp_lmfit_gaussian_centering(subArraySize=6, recheckMethod=None, median_crop=False)\n",
    "print('Operation took {} seconds with {} cores'.format(time()-start, example_wanderer_median.nCores))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the sigma-clipping outliers for plotting purpose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nSig       = 10.1\n",
    "medY       = median(example_wanderer_median.centering_GaussianFit.T[y])\n",
    "medX       = median(example_wanderer_median.centering_GaussianFit.T[x])\n",
    "stdY       = std(example_wanderer_median.centering_GaussianFit.T[y])\n",
    "stdX       = std(example_wanderer_median.centering_GaussianFit.T[x])\n",
    "\n",
    "ySig = 4\n",
    "xSig = 4\n",
    "outliers   = (((example_wanderer_median.centering_GaussianFit.T[y] - medY)/(ySig*stdY))**2 + \\\n",
    "              ((example_wanderer_median.centering_GaussianFit.T[x] - medX)/(xSig*stdX))**2) > 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the inliers (blue) vs outliers (not blue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = figure().add_subplot(111)\n",
    "cx, cy = example_wanderer_median.centering_GaussianFit.T[x],example_wanderer_median.centering_GaussianFit.T[y]\n",
    "ax.plot(cx,cy,'.',ms=1)\n",
    "ax.plot(cx[outliers],cy[outliers],'.',ms=1)\n",
    "# ax.plot(median(cx), median(cy),'ro',ms=1)\n",
    "ax.set_xlim(medX-nSig*stdX,medX+nSig*stdX)\n",
    "ax.set_ylim(medY-nSig*stdY,medY+nSig*stdY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use advanced clustering algorithms (DBSCAN) to determine the inliers vs outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dbs     = DBSCAN(n_jobs=-1, eps=0.2, leaf_size=10)\n",
    "dbsPred = dbs.fit_predict(example_wanderer_median.centering_GaussianFit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dbs_options = [k for k in range(-1,100) if (dbsPred==k).sum()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the full extent of the data to show that DBSCAN was able to identify the inliers correctly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(6,6))\n",
    "ax  = fig.add_subplot(111)\n",
    "\n",
    "medGaussCenters   = median(example_wanderer_median.centering_GaussianFit,axis=0)\n",
    "sclGaussCenters   = scale.mad(example_wanderer_median.centering_GaussianFit)\n",
    "sclGaussCenterAvg = np.sqrt(((sclGaussCenters**2.).sum()))\n",
    "\n",
    "yctrs = example_wanderer_median.centering_GaussianFit.T[y]\n",
    "xctrs = example_wanderer_median.centering_GaussianFit.T[x]\n",
    "\n",
    "nSigmas         = 5\n",
    "for nSig in linspace(1,10,10):\n",
    "    CircularAperture(medGaussCenters[::-1],nSig*sclGaussCenterAvg).plot(ax=ax)\n",
    "\n",
    "for dbsOpt in dbs_options:\n",
    "    ax.plot(xctrs[dbsPred==dbsOpt], yctrs[dbsPred==dbsOpt],'.',zorder=0, ms=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure that there are only a handful (<< 1%) of outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "npix = 3\n",
    "\n",
    "stillOutliers = np.where(abs(example_wanderer_median.centering_GaussianFit - medGaussCenters) > 4*sclGaussCenterAvg)[0]\n",
    "print(len(stillOutliers))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select the \"class\" dbsClean == 0 for the `inliers`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dbsClean  = 0\n",
    "dbsKeep   = (dbsPred == dbsClean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nCores = example_wanderer_median.nCores\n",
    "start = time()\n",
    "example_wanderer_median.mp_measure_background_annular_mask()\n",
    "print('AnnularBG took {} seconds with {} cores'.format(time() - start, nCores))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the background to make sure that the (to be subtracted) flux is stable overtime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = figure(figsize=(20,10))\n",
    "ax  = fig.add_subplot(111)\n",
    "ax.plot(example_wanderer_median.timeCube, example_wanderer_median.background_CircleMask,'.',alpha=0.2)\n",
    "ax.plot(example_wanderer_median.timeCube, example_wanderer_median.background_Annulus,'.',alpha=0.2)\n",
    "ax.plot(example_wanderer_median.timeCube, example_wanderer_median.background_MedianMask,'.',alpha=0.2)\n",
    "ax.plot(example_wanderer_median.timeCube, example_wanderer_median.background_KDEUniv,'.',alpha=0.2)\n",
    "ax.axvline(example_wanderer_median.timeCube.min()-.01+0.02)\n",
    "ax.set_ylim(-25,100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the `effective widths` of each image to use later as the \"beta pixels\" and \"optimal apertures\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_wanderer_median.measure_effective_width()\n",
    "print(example_wanderer_median.effective_widths.mean(), sqrt(example_wanderer_median.effective_widths).mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Pipeline took {} seconds thus far'.format(time() - startFull))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the time series with static aperture radii only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print('Iterating over Background Techniques, Centering Techniques, Aperture Radii' + '\\n')\n",
    "centering_choices  = ['Gaussian_Fit']#, 'Gaussian_Mom', 'FluxWeighted']#, 'LeastAsymmetry']\n",
    "background_choices = ['AnnularMask']#example_wanderer_median.background_df.columns\n",
    "staticRads         = np.arange(1, 6,0.5)#[1.0 ]# aperRads = np.arange(1, 6,0.5)\n",
    "\n",
    "for staticRad in tqdm_notebook(staticRads, total=len(staticRads), desc='Static'):\n",
    "    example_wanderer_median.mp_compute_flux_over_time_varRad(staticRad, varRad=None, centering_choices[0], background_choices[0], useTheForce=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create Beta Variable Radius"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_wanderer_median.mp_compute_flux_over_time_betaRad()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Entire Pipeline took {} seconds'.format(time() - startFull))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Advanced clustering algorithms `DBSCAN` to compute the outliers of the flux distribution. This is sensitive the structure in the data (i.e. transit vs outlier), which is not always true with sigma-clipping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "example_wanderer_median.mp_DBScan_Flux_All()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that the majority of data is an `inlier`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inlier_master = array(list(example_wanderer_median.inliers_Phots.values())).mean(axis=0) == 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "((~inlier_master).sum() / inlier_master.size)*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the PLD components -- normalized and store the PLD vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_wanderer_median.extract_PLD_components()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Advanced clustering algorithms `DBSCAN` to compute the outliers of the PLD distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_wanderer_median.mp_DBScan_PLD_All()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save all of your progress per AOR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print('Saving `example_wanderer_median` to a set of pickles for various Image Cubes and the Storage Dictionary')\n",
    "\n",
    "savefiledir         = environ['HOME']+'/Research/Planets/'+planetName+'/ExtractedData/' + channel \n",
    "saveFileNameHeader  = planetName+'_'+ AORNow +'_Median'\n",
    "saveFileType        = '.joblib.save'\n",
    "\n",
    "if not path.exists(environ['HOME']+'/Research/Planets/'+planetName+'/ExtractedData/'):\n",
    "    mkdir(environ['HOME']+'/Research/Planets/'+planetName+'/ExtractedData/')\n",
    "\n",
    "if not path.exists(savefiledir):\n",
    "    print('Creating ' + savefiledir)\n",
    "    mkdir(savefiledir)\n",
    "\n",
    "print()\n",
    "print('Saving to ' + savefiledir + saveFileNameHeader + saveFileType)\n",
    "print()\n",
    "\n",
    "example_wanderer_median.save_data_to_save_files(savefiledir=savefiledir, \\\n",
    "                                                saveFileNameHeader=saveFileNameHeader, \\\n",
    "                                                saveFileType=saveFileType)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the RMS in the raw data as a function of the apeture radius"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_cycle = rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "ax = figure().add_subplot(111)\n",
    "for key in example_wanderer_median.flux_TSO_df.keys():\n",
    "    aperRad = float(key.split('_')[-2])\n",
    "    ax.scatter(aperRad, scale.mad(np.diff(example_wanderer_median.flux_TSO_df[key])), color=color_cycle[0])\n",
    "\n",
    "ax.set_xlabel('Aperture Radius')\n",
    "ax.set_ylabel('MAD( Diff ( Flux ) )')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Entire Pipeline took {} seconds'.format(time() - startFull))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert Wanderer output to Skywalker input"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I made up this loop and did not test it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keys for `skywalker` input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_key = 'phots'\n",
    "time_key = 'times'\n",
    "flux_err_key = 'noise'\n",
    "eff_width_key = 'npix'\n",
    "pld_coeff_key = 'pld'\n",
    "ycenter_key = 'ycenters'\n",
    "xcenter_key = 'xcenters'\n",
    "ywidth_key = 'ywidths'\n",
    "xwidth_key = 'xwidths'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Things that **DON'T** change with respect to aperture radii"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeCube = example_wanderer_median.timeCube\n",
    "phots_array = example_wanderer_median.flux_TSO_df.values\n",
    "PLDFeatures = example_wanderer_median.PLD_components.T\n",
    "\n",
    "try:\n",
    "    inliers_Phots = example_wanderer_median.inliers_Phots.values()\n",
    "except:\n",
    "    inliers_Phots = np.ones(photsLocal.shape)\n",
    "\n",
    "try:\n",
    "    inliers_PLD = example_wanderer_median.inliers_PLD.values()\n",
    "except:\n",
    "    inliers_PLD = np.ones(PLDFeatureLocal.shape)\n",
    "\n",
    "inliersMaster = array(list(inliers_Phots)).all(axis=0) # Need to Switch `axis=0` for Qatar-2\n",
    "inliersMaster = inliersMaster * inliers_PLD.all(axis=1)\n",
    "\n",
    "nSig = 6 # vary this as desired for 3D sigma clipping double check\n",
    "\n",
    "ypos, xpos = clipOutlier2D(transpose([example_wanderer_median.centering_GaussianFit.T[y][inliersMaster], \\\n",
    "                                           example_wanderer_median.centering_GaussianFit.T[x][inliersMaster]])).T\n",
    "\n",
    "npix = sqrt(example_wanderer_median.effective_widths[inliersMaster])\n",
    "time_c = timeCube[inliersMaster]\n",
    "ywidths_c, xwidths_c = example_wanderer_median.widths_GaussianFit[inliersMaster].T\n",
    "pld_comp_c = example_wanderer_median.PLD_components.T # this is new to Carlos's notebook instance\n",
    "pld_output_c = np.array(list([time_c]) + list(pld_comp_c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Things that **DO** change with respect to aperture radii"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for phot_select, key_flux_now in tqdm(enumerate(example_wanderer_median.flux_TSO_df.keys())):\n",
    "    if key_flux_now[-3:] == 0.0: # only do static radii\n",
    "        flux_c = example_wanderer_median.flux_TSO_df[key_flux_now].values[inliersMaster]\n",
    "        noise_c = example_wanderer_median.noise_TSO_df[key_flux_now].values[inliersMaster]\n",
    "        \n",
    "        output_dict = {time_key: time_c, \n",
    "                       flux_key: flux_c, \n",
    "                       flux_err_key: noise_c, \n",
    "                       eff_width_key: npix_c, \n",
    "                       xcenter_key: xpos_c, \n",
    "                       ycenter_key: ypos_c, \n",
    "                       xwidth_key: xwidth_c, \n",
    "                       ywidth_key: ywidth_c, \n",
    "                       pld_coeff_key: pld_comp_c}\n",
    "        \n",
    "        # This creates 1 joblib output file for one static aperture radius -- need to be cycled from above: change `staticRad = '2.5'` to new radius\n",
    "        joblib.dump(output_dict, '{}_full_output_for_skywalker_pipeline_{}_{}_{}.joblib.save'.format(planet_dir_name, channel, staticRad, varRad))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The following code is a copy/paste from a different notebook of mine.\n",
    "\n",
    "This is the code I used to make the for loop above  \n",
    "If the for loop does not work, try / check this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeCube = example_wanderer_median.timeCube\n",
    "phots_array = example_wanderer_median.flux_TSO_df.values\n",
    "PLDFeatures = example_wanderer_median.PLD_components.T\n",
    "\n",
    "try:\n",
    "    inliers_Phots = example_wanderer_median.inliers_Phots.values()\n",
    "except:\n",
    "    inliers_Phots = np.ones(photsLocal.shape)\n",
    "\n",
    "try:\n",
    "    inliers_PLD = example_wanderer_median.inliers_PLD.values()\n",
    "except:\n",
    "    inliers_PLD = np.ones(PLDFeatureLocal.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gaussian_Fit_AnnularMask_rad_2.5_0.0\n",
    "\n",
    "staticRad = '2.5' # Need to cycle over all possible values here: [1.0, 1.5, 2.0, ..., 5.5]\n",
    "varRad = '0.0'\n",
    "key_flux_now = 'Gaussian_Fit_AnnularMask_rad_'+staticRad+'_'+varRad\n",
    "phot_select = np.where(example_wanderer_median.flux_TSO_df.keys() == key_flux_now)[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "inliersMaster = array(list(inliers_Phots)).all(axis=0) # Need to Switch `axis=0` for Qatar-2\n",
    "inliersMaster = inliersMaster * inliers_PLD.all(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "nSig = 6 # vary this as desired\n",
    "\n",
    "if inliersMaster.all():\n",
    "    # If inliersMaster keeps ALL values, then double check with 3D inlier flagging\n",
    "    print('Working on AOR {}'.format(AORNow))\n",
    "    cy_now, cx_now        = example_wanderer_median.centering_GaussianFit.T\n",
    "    phots_now             = phots_array[:,phot_select]\n",
    "    \n",
    "    phots_clipped         = clipOutlier2D(phots_now, nSig=nSig)\n",
    "    cy_clipped, cx_clipped= clipOutlier2D(transpose([cy_now, cx_now]),nSig=nSig).T\n",
    "    arr2D_clipped         = transpose([phots_clipped, cy_clipped, cx_clipped])\n",
    "    \n",
    "    # 3D inlier selection\n",
    "    inliersMaster = (phots_clipped == phots_now)*(cy_clipped==cy_now)*(cx_clipped==cx_now)\n",
    "else:\n",
    "    print(\"this box is just to double check -- keeping all inlier flags from above\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "ypos, xpos = clipOutlier2D(transpose([example_wanderer_median.centering_GaussianFit.T[y][inliersMaster], \\\n",
    "                                           example_wanderer_median.centering_GaussianFit.T[x][inliersMaster]])).T\n",
    "\n",
    "npix = sqrt(example_wanderer_median.effective_widths[inliersMaster])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_c = phots_array[:, phot_select][inliersMaster]\n",
    "\n",
    "# noise_c = np.sqrt(flux_c) # Photon limit\n",
    "noise_c = example_wanderer_median.noise_TSO_df[key_flux_now].values[inliersMaster]\n",
    "\n",
    "time_c = timeCube[inliersMaster]\n",
    "\n",
    "ywidths_c, xwidths_c = example_wanderer_median.widths_GaussianFit[inliersMaster].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I am guessing that this will work.\n",
    "# I'm keeping the commented line because that's what I used before\n",
    "# pld_comp_c = wanderer.extract_PLD_components(example_wanderer_median.imageCube, order=1)\n",
    "\n",
    "pld_comp_c = example_wanderer_median.PLD_components.T # this is new to Carlos's notebook instance\n",
    "pld_output_c = np.array(list([time_c]) + list(pld_comp_c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['qatar2_full_output_for_pipeline_ch2_2.5_0.0.joblib.save']"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flux_key = 'phots'\n",
    "time_key = 'times'\n",
    "flux_err_key = 'noise'\n",
    "eff_width_key = 'npix'\n",
    "pld_coeff_key = 'pld'\n",
    "ycenter_key = 'ycenters'\n",
    "xcenter_key = 'xcenters'\n",
    "ywidth_key = 'ywidths'\n",
    "xwidth_key = 'xwidths'\n",
    "\n",
    "output_dict = {time_key: time_c, \n",
    "               flux_key: flux_c, \n",
    "               flux_err_key: noise_c, \n",
    "               eff_width_key: npix_c, \n",
    "               xcenter_key: xpos_c, \n",
    "               ycenter_key: ypos_c, \n",
    "               xwidth_key: xwidth_c, \n",
    "               ywidth_key: ywidth_c, \n",
    "               pld_coeff_key: pld_comp_c}\n",
    "\n",
    "# This creates 1 joblib output file for one static aperture radius -- need to be cycled from above: change `staticRad = '2.5'` to new radius\n",
    "joblib.dump(output_dict, '{}_full_output_for_skywalker_pipeline_{}_{}_{}.joblib.save'.format(planet_dir_name, channel, staticRad, varRad))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "astro",
   "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.10.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "cec8ed39d91a2311e942fb88686183bc3c16ab31df95d16c9cc802d6a62ec85b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}