Exoplanet_TSO_-_Combined_PLD_Example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import batman\n",
    "import corner\n",
    "# import datetime\n",
    "# import emcee3\n",
    "import emcee\n",
    "# import jdcal\n",
    "from os import environ, path\n",
    "import pywt\n",
    "# import spiderman as sp\n",
    "\n",
    "from batman    import TransitModel, TransitParams\n",
    "from spiderman import ModelParams as sp_ModelParams\n",
    "\n",
    "from functools import partial\n",
    "from pylab     import *\n",
    "from pandas    import DataFrame\n",
    "\n",
    "from scipy.spatial           import cKDTree\n",
    "from sklearn.externals       import joblib\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.decomposition   import PCA, FastICA\n",
    "\n",
    "from astropy.io import fits\n",
    "from glob import glob\n",
    "\n",
    "\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "# from scipy.optimize    import leastsq, minimize\n",
    "from scipy.interpolate import CubicSpline\n",
    "# from scipy.signal      import medfilt\n",
    "# from scipy.stats       import binned_statistic\n",
    "\n",
    "from lmfit import Parameters, Minimizer, report_errors\n",
    "\n",
    "from statsmodels.robust.scale import mad\n",
    "from sklearn.preprocessing    import scale\n",
    "\n",
    "from time import time\n",
    "\n",
    "# from sys import argv\n",
    "\n",
    "# from photutils import CircularAperture, CircularAnnulus, EllipticalAperture\n",
    "# from photutils import aperture_photometry\n",
    "\n",
    "# from statsmodels.robust import scale\n",
    "from datetime import datetime\n",
    "\n",
    "from exoparams import PlanetParams\n",
    "\n",
    "from numpy import cos, pi, abs\n",
    "\n",
    "# from sklearn.svm import SVR\n",
    "from sklearn.cluster import DBSCAN\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "from multiprocessing import cpu_count, Pool\n",
    "\n",
    "from astropy.constants import R_sun, au\n",
    "\n",
    "from wanderer import wanderer\n",
    "\n",
    "from spitzer_helper_functions import clipOutliers, bin_array, b2inc, deltaphase_eclipse\n",
    "from spitzer_helper_functions import extract_PLD_components, de_median, savePickleOut, spiderman_lmfit_model\n",
    "from spitzer_helper_functions import phase_cos_curve, phase_cos_sin_curve, inc2b, batman_lmfit_model\n",
    "# from spitzer_helper_functions import pixel_level_decorrelation_instrument_profile\n",
    "\n",
    "stdScaler = StandardScaler()\n",
    "stime = time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import rcParams\n",
    "rcParams[\"savefig.dpi\"] = 200\n",
    "rcParams[\"figure.dpi\"] = 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(\"\\n**DONE LOADING LIBRARIES AND DEFINING FUNCTIONS**\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# List of Function Definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batman_lmfit_model(period, tCenter, inc, aprs, edepth, tdepth, ecc, omega, times, u1=None, u2=None,\n",
    "                       ldtype = 'uniform', transittype=\"primary\", bm_params=None):\n",
    "    \n",
    "    if tdepth is not 0.0 or edepth is not 0.0:\n",
    "        if bm_params is None:\n",
    "            bm_params = TransitParams() # object to store transit parameters\n",
    "        \n",
    "        bm_params.per       = period  # orbital period\n",
    "        bm_params.t0        = tCenter # time of inferior conjunction\n",
    "        bm_params.a         = aprs    # semi-major axis (in units of stellar radii)\n",
    "        bm_params.fp        = edepth  # f\n",
    "        bm_params.tdepth    = tdepth  # from Fraine et al. 2014s\n",
    "        bm_params.rp        = sqrt(tdepth) # planet radius (in units of stellar radii)\n",
    "        bm_params.ecc       = ecc     # eccentricity\n",
    "        bm_params.w         = omega   # longitude of periastron (in degrees)\n",
    "        bm_params.inc       = inc     # orbital inclination (in degrees)\n",
    "        bm_params.limb_dark = ldtype  # limb darkening model # NEED TO FIX THIS\n",
    "        bm_params.u         = []      # limb darkening coefficients # NEED TO FIX THIS\n",
    "        \n",
    "        if u1 is not None and ldtype is not 'uniform':\n",
    "            bm_params.u.append(u1)\n",
    "        elif u1 is not None and ldtype is 'uniform':\n",
    "            raise ValueError('If you set `u1`, you must also set `ldtype` to either linear or quadratic')\n",
    "        if u2 is not None and ldtype is 'quadratic':\n",
    "            bm_params.u.append(u2)\n",
    "        elif u2 is not None and ldtype is not 'quadratic':\n",
    "            raise ValueError('If you set `u2`, you must also set `ldtype` quadratic')\n",
    "        \n",
    "        bm_params.delta_phase = deltaphase_eclipse(bm_params.ecc, bm_params.w) if bm_params.ecc is not 0.0 else 0.5\n",
    "        bm_params.t_secondary = bm_params.t0 + bm_params.per*bm_params.delta_phase\n",
    "        \n",
    "        m_eclipse = TransitModel(bm_params, times, transittype=transittype).light_curve(bm_params)\n",
    "    else:\n",
    "        return ones(times.size)\n",
    "    \n",
    "    return m_eclipse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spiderman_lmfit_model(period, tCenter, inc, aprs, tdepth, ecc, omega, times, u1, u2, xi, T_night, delta_T, \n",
    "                         nspiders=1000, nCores = cpu_count(),\n",
    "                         spider_params=None, n_layers=20, stellar_radius = 1.0, stellar_temperature=5500.,\n",
    "                         wave_low = 4.0, wave_hi = 5.0, brightness_model='zhang', if_eclipse=False):\n",
    "    \n",
    "    if tdepth is not 0.0:\n",
    "        # if spider_params is None:\n",
    "        #     spider_params                = sp_ModelParams(brightness_model=brightness_model)\n",
    "        #     spider_params.n_layers       = n_layers\n",
    "        #     spider_params.stellar_radius = stellar_radius\n",
    "        #     spider_params.T_s            = stellar_temperature\n",
    "        #     spider_params.l1             = wave_low\n",
    "        #     spider_params.l2             = wave_hi\n",
    "        #     spider_params.eclipse        = if_eclipse\n",
    "        \n",
    "        a_au  = aprs * spider_params.stellar_radius * R_sun.value / au.value\n",
    "        # inc   = b2inc(bImpact, aprs, ecc, omega)*180/pi\n",
    "        \n",
    "        spider_params.t0      = tCenter      # Central time of PRIMARY transit [days]\n",
    "        spider_params.per     = period       # Period [days]\n",
    "        spider_params.a_abs   = a_au         # The absolute value of the semi-major axis [AU]\n",
    "        spider_params.inc     = inc # orbital inclination (in degrees)\n",
    "        spider_params.ecc     = ecc          # Eccentricity\n",
    "        spider_params.w       = omega        # Argument of periastron\n",
    "        spider_params.rp      = sqrt(tdepth) # planet radius (in units of stellar radii)\n",
    "        spider_params.a       = aprs         # Semi-major axis scaled by stellar radius\n",
    "        spider_params.p_u1    = u1           # Planetary limb darkening parameter\n",
    "        spider_params.p_u2    = u2           # Planetary limb darkening parameter\n",
    "        \n",
    "        spider_params.xi      = xi           # Ratio of radiative to advective timescale             \n",
    "        spider_params.T_n     = T_night      # Temperature of nightside\n",
    "        spider_params.delta_T = delta_T      # Day-night temperature contrast\n",
    "        \n",
    "        nskips            = times.size // nspiders\n",
    "        times_subsampled  = times[::nskips]\n",
    "        spider_lightcurve = spider_params.lightcurve(times_subsampled)\n",
    "        spider_lightcurve = CubicSpline(times_subsampled, spider_lightcurve)\n",
    "        \n",
    "        return spider_lightcurve(times)\n",
    "    else:\n",
    "        return ones(times.size)\n",
    "    \n",
    "    raise Exception('Something Weird happened')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_qhull_one_point(point, x0, y0, np0, inds, sm_num):\n",
    "    # ind         = kdtree.query(kdtree.data[point],sm_num+1)[1][1:]\n",
    "    ind = inds[point]\n",
    "    dx  = x0[ind] - x0[point]\n",
    "    dy  = y0[ind] - y0[point]\n",
    "    \n",
    "    if np0.sum() != 0.0:\n",
    "        dnp         = np0[ind] - np0[point]\n",
    "    \n",
    "    sigx  = np.std(dx )\n",
    "    sigy  = np.std(dy )\n",
    "    \n",
    "    if dnp.sum() != 0.0:\n",
    "                \n",
    "        signp   = np.std(dnp)\n",
    "        gw_temp = np.exp(-dx**2./(2.0*sigx**2.)) * \\\n",
    "                      np.exp(-dy**2./(2.*sigy**2.))  * \\\n",
    "                      np.exp(-dnp**2./(2.*signp**2.))\n",
    "    else:\n",
    "        gw_temp = np.exp(-dx**2./(2.0*sigx**2.)) * \\\n",
    "                  np.exp(-dy**2./(2.*sigy**2.))\n",
    "    \n",
    "    gw_sum  = gw_temp.sum()\n",
    "    \n",
    "    return gw_temp/gw_sum#, ind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mp_find_nbr_qhull(xpos, ypos, npix = None, inds = None, n_nbr = 50, returnInds=False,\n",
    "                      a = 1.0, b = 0.7, c = 1.0, expansion = 1000., nCores=cpu_count()):\n",
    "    '''\n",
    "        Python Implimentation of N. Lewis method, described in Lewis etal 2012, Knutson etal 2012, Fraine etal 2013\n",
    "        \n",
    "        Taken from N. Lewis IDL code:\n",
    "            \n",
    "            Construct a 3D surface (or 2D if only using x and y) from the data\n",
    "            using the qhull.pro routine.  Save the connectivity information for\n",
    "            each data point that was used to construct the Delaunay triangles (DT)\n",
    "            that form the grid.  The connectivity information allows us to only\n",
    "            deal with a sub set of data points in determining nearest neighbors\n",
    "            that are either directly connected to the point of interest or\n",
    "            connected through a neighboring point\n",
    "        \n",
    "        Python Version:\n",
    "            J. Fraine    first edition, direct translation from IDL 12.05.12\n",
    "    '''\n",
    "    #The surface fitting performs better if the data is scattered about zero\n",
    "    if npix is not None and bool(c):\n",
    "        np0 = np.sqrt(npix)\n",
    "        np0 = (np0 - np.median(np0))/c\n",
    "        features  = np.transpose((y0, x0, np0))\n",
    "    else:\n",
    "        features  = np.transpose((y0, x0))\n",
    "        \n",
    "        if np.sum(np0) == 0.0:\n",
    "            print('SKIPPING Noise Pixel Sections of Gaussian Kernel because Noise Pixels are Zero')\n",
    "        if c == 0:\n",
    "            print('SKIPPING Noise Pixel Sections of Gaussian Kernel because c == 0')\n",
    "    \n",
    "    x0  = (xpos - np.median(xpos))/a\n",
    "    y0  = (ypos - np.median(ypos))/b\n",
    "    \n",
    "    k   = sm_num                           # This is the number of nearest neighbors you want\n",
    "    n   = x0.size                          # This is the number of data points you have\n",
    "    gw  = np.zeros((k,n),dtype=np.float64) # This is the gaussian weight for each data point determined from the nearest neighbors\n",
    "    \n",
    "    if inds is None:\n",
    "        kdtree    = cKDTree(features * expansion) #Multiplying `features` by 1000.0 avoids precision problems\n",
    "        inds      = kdtree.query(kdtree.data, n_nbr+1)[1][:,1:]\n",
    "        \n",
    "        print('WARNING: Because `inds` was not provided, we must now compute and return it here')\n",
    "        returnInds= True\n",
    "    \n",
    "    func  = partial(find_qhull_one_point, x0=x0, y0=y0, np0=np0, inds=inds, sm_num=sm_num)\n",
    "    \n",
    "    pool  = Pool(nCores)\n",
    "    \n",
    "    gw_pool = pool.starmap(func, zip(range(n)))\n",
    "    \n",
    "    pool.close()\n",
    "    pool.join()\n",
    "    \n",
    "    if returnInds:\n",
    "        return np.array(gw_pool), inds\n",
    "    else:\n",
    "        return np.array(gw_pool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.spatial      import cKDTree\n",
    "\n",
    "def find_nbr_qhull(xpos, ypos, npix, sm_num = 100, a = 1.0, b = 0.7, c = 1.0, print_space = 10000.):\n",
    "    '''\n",
    "        Python Implimentation of N. Lewis method, described in Lewis etal 2012, Knutson etal 2012, Fraine etal 2013\n",
    "        \n",
    "        Taken from N. Lewis IDL code:\n",
    "            \n",
    "            Construct a 3D surface (or 2D if only using x and y) from the data\n",
    "            using the qhull.pro routine.  Save the connectivity information for\n",
    "            each data point that was used to construct the Delaunay triangles (DT)\n",
    "            that form the grid.  The connectivity information allows us to only\n",
    "            deal with a sub set of data points in determining nearest neighbors\n",
    "            that are either directly connected to the point of interest or\n",
    "            connected through a neighboring point\n",
    "        \n",
    "        Python Version:\n",
    "            J. Fraine    first edition, direct translation from IDL 12.05.12\n",
    "    '''\n",
    "    from scipy import spatial\n",
    "    #The surface fitting performs better if the data is scattered about zero\n",
    "    \n",
    "    npix    = np.sqrt(npix)\n",
    "    \n",
    "    x0  = (xpos - np.median(xpos))/a\n",
    "    y0  = (ypos - np.median(ypos))/b\n",
    "    \n",
    "    if np.sum(npix) != 0.0 and c != 0:\n",
    "        np0 = (npix - np.median(npix))/c\n",
    "    else:\n",
    "        if np.sum(npix) == 0.0:\n",
    "            print('SKIPPING Noise Pixel Sections of Gaussian Kernel because Noise Pixels are Zero')\n",
    "        if c == 0:\n",
    "            print('SKIPPING Noise Pixel Sections of Gaussian Kernel because c == 0')\n",
    "    \n",
    "    k            = sm_num                           # This is the number of nearest neighbors you want\n",
    "    n            = x0.size                          # This is the number of data points you have\n",
    "    nearest      = np.zeros((k,n),dtype=np.int64)   # This stores the nearest neighbors for each data point\n",
    "    \n",
    "    #Multiplying by 1000.0 avoids precision problems\n",
    "    if npix.sum() != 0.0 and c != 0:\n",
    "        kdtree  = cKDTree(np.transpose((y0*1000., x0*1000., np0*1000.)))\n",
    "    else:\n",
    "        kdtree  = cKDTree(np.transpose((y0*1000., x0*1000.)))\n",
    "    \n",
    "    gw  = np.zeros((k,n),dtype=np.float64) # This is the gaussian weight for each data point determined from the nearest neighbors\n",
    "    \n",
    "    start   = time()\n",
    "    for point in tqdm_notebook(range(n),total=n):\n",
    "        ind         = kdtree.query(kdtree.data[point],sm_num+1)[1][1:]\n",
    "        dx          = x0[ind] - x0[point]\n",
    "        dy          = y0[ind] - y0[point]\n",
    "        \n",
    "        if npix.sum() != 0.0 and c != 0:\n",
    "            dnp         = np0[ind] - np0[point]\n",
    "        \n",
    "        sigx        = np.std(dx )\n",
    "        sigy        = np.std(dy )\n",
    "        if npix.sum() != 0.0 and c != 0:\n",
    "            signp       = np.std(dnp)\n",
    "        if npix.sum() != 0.0 and c != 0:\n",
    "            gw_temp     = np.exp(-dx**2./(2.0*sigx**2.)) * \\\n",
    "                          np.exp(-dy**2./(2.*sigy**2.))  * \\\n",
    "                          np.exp(-dnp**2./(2.*signp**2.))\n",
    "        else:\n",
    "            gw_temp     = np.exp(-dx**2./(2.0*sigx**2.)) * \\\n",
    "                          np.exp(-dy**2./(2.*sigy**2.))\n",
    "        \n",
    "        gw_sum      = gw_temp.sum()\n",
    "        gw[:,point] = gw_temp/gw_sum\n",
    "        \n",
    "        #if (gw_sum == 0.0) or ~np.isfinite(gw_sum):\n",
    "        #    raise Exception('(gw_sum == 0.0) or ~isfinite(gw_temp))')\n",
    "        \n",
    "        nearest[:,point]  = ind\n",
    "    \n",
    "    return gw.transpose(), nearest.transpose() # nearest  == nbr_ind.transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pixel_level_decorrelation_instrument_profile(times, input_features, \n",
    "                              pld1_l, pld2_l, pld3_l, pld4_l, pld5_l, pld6_l, pld7_l, pld8_l, pld9_l,\n",
    "                              pld1_q, pld2_q, pld3_q, pld4_q, pld5_q, pld6_q, pld7_q, pld8_q, pld9_q):# ,\n",
    "                              # intcpt=1.0, slope=0.0, crvtur=0.0):\n",
    "    \n",
    "    feature_coeffs   = [pld1_l, pld2_l, pld3_l, pld4_l, pld5_l, pld6_l, pld7_l, pld8_l, pld9_l, \\\n",
    "                        pld1_q, pld2_q, pld3_q, pld4_q, pld5_q, pld6_q, pld7_q, pld8_q, pld9_q]\n",
    "    \n",
    "    instrumental  = dot(input_features, feature_coeffs) # This is now quadratic\n",
    "    \n",
    "    return instrumental"
   ]
  },
  {
   "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(((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": [
    "def plot_MAD_AperRads(instance, minRad=None, maxRad=None, varRadFlag=True):\n",
    "    color_cycle = cycler(rcParams['axes.prop_cycle']).by_key()['color']\n",
    "    \n",
    "    quad_width= instance.quadrature_widths.values\n",
    "    vrad_dist = quad_width - np.median(quad_width)\n",
    "    vrad_dist = clipOutlier2D(vrad_dist, nSig=5)\n",
    "    vrad_dist_med = np.median(vrad_dist)\n",
    "    \n",
    "    betaColor = 7\n",
    "    quadColor = 8\n",
    "    \n",
    "    ax = figure().add_subplot(111)\n",
    "    \n",
    "    for key in instance.flux_TSO_df.keys():\n",
    "        staticRad = float(key.split('_')[-2])\n",
    "        varRad    = float(key.split('_')[-1])\n",
    "        aperRad   = staticRad + varRad*vrad_dist_med\n",
    "        colorNow  = color_cycle[int(varRad*4)]\n",
    "        \n",
    "        if 'betaRad' in key:\n",
    "            aperRad    = median(sqrt(instance.effective_widths))\n",
    "            colorNow  = color_cycle[int(betaColor)]\n",
    "        \n",
    "        if 'quadRad' in key:\n",
    "            aperRad   = 2*sqrt(2*log(2))*median(instance.quadrature_widths.values)\n",
    "            colorNow  = color_cycle[int(quadColor)]\n",
    "        \n",
    "        if minRad is not  None and maxRad is not None:\n",
    "            if aperRad > minRad and aperRad < maxRad:\n",
    "                ax.scatter(aperRad, mad(np.diff(instance.flux_TSO_df[key])), color=colorNow, zorder=int(varRad*4))\n",
    "        else:\n",
    "            ax.scatter(aperRad, mad(np.diff(instance.flux_TSO_df[key])), color=colorNow, zorder=int(varRad*4))\n",
    "    \n",
    "    for varRad in [0.,0.25, 0.5, 0.75, 1.0, 1.25, 1.5]:\n",
    "        colorNow  = color_cycle[int(varRad*4)]\n",
    "        ax.scatter([],[], color=colorNow, label=varRad)\n",
    "    \n",
    "    ax.scatter([],[],color=color_cycle[int(betaColor)], label='Beta')\n",
    "    ax.scatter([],[],color=color_cycle[int(quadColor)], label='Quad')\n",
    "    \n",
    "    ax.set_xlabel('StaticRad + Average(varRad)')\n",
    "    ax.set_ylabel('MAD( Diff ( Flux ) )')\n",
    "    ax.legend(loc=0,fontsize=10)\n",
    "    \n",
    "    ax.set_title('' if not hasattr(instance, '__name__') else instance.__name__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_model_over_reduced_lc_half_PLDsq_universal(times, phots, phots_err, input_features, in_eclipse, fit_values,\n",
    "                                                    spider_params = None, refTime=0.0, figsize=None, nbins=200, \n",
    "                                                    plotRawData=False, alpha=1.0, phase_method='fourier', \n",
    "                                                    baseline='PLD',ax = None, returnAx = False, \n",
    "                                                    model_color=0, data_color=1):\n",
    "    \n",
    "    color_cycle = cycler(rcParams['axes.prop_cycle']).by_key()['color']\n",
    "    model_color = color_cycle[model_color] if isinstance(model_color, int) else model_color\n",
    "    data_color  = color_cycle[data_color]  if isinstance(data_color , int) else data_color \n",
    "    \n",
    "    # print('model:{}\\tdata:{}'.format(model_color, data_color))\n",
    "    \n",
    "    clean_model = model_selector(fit_values, times, input_features, in_eclipse, spider_params=spider_params, \n",
    "                                 phase_method=phase_method, baseline=baseline, model_only=True)\n",
    "    \n",
    "    model_params_keys = np.sort(list(fit_values.keys()))\n",
    "    \n",
    "    feature_coeffs  = [fit_values[key].value for key in model_params_keys if 'pld' in key and '_l' in key] + \\\n",
    "                      [fit_values[key].value for key in model_params_keys if 'pld' in key and '_q' in key]\n",
    "    \n",
    "    # Check if line fit parameters exist, then allocate them; or use defaults\n",
    "    intcpt  = fit_values['intcpt'] if 'intcpt' in fit_values.keys() else 1.0\n",
    "    slope   = fit_values['slope']  if 'slope'  in fit_values.keys() else 0.0\n",
    "    crvtur  = fit_values['crvtur'] if 'crvtur' in fit_values.keys() else 0.0\n",
    "    \n",
    "    instrumental = intcpt * np.ones(times.size)\n",
    "    if slope  != 0.0:\n",
    "        instrumental += slope * (times-times.mean())\n",
    "    if crvtur != 0.0:\n",
    "        instrumental += crvtur * (times-times.mean())**2.\n",
    "    \n",
    "    if baseline == 'PLD':\n",
    "        instrumental += pixel_level_decorrelation_instrument_profile(times, input_features, \n",
    "                        *feature_coeffs)#, intcpt=intcpt, slope=slope, crvtur=crvtur)\n",
    "    \n",
    "    if baseline == 'KRData':\n",
    "        phots_now, gk, nbr_ind = input_features\n",
    "        instrumental *= np.sum((phots_now/clean_model)[nbr_ind]*gk,axis=1)\n",
    "    \n",
    "    # clean_model  = whole_model / instrumental \n",
    "    \n",
    "    binsize = phots.size // nbins\n",
    "    \n",
    "    reduced = phots/instrumental\n",
    "    diff_dm = 0.0#np.median(reduced - clean_model)\n",
    "    \n",
    "    print('Diff DM: {:0f}'.format(diff_dm*ppm))\n",
    "    \n",
    "    bin_flux    , _ = bin_array((reduced - diff_dm), binsize=binsize)\n",
    "    bin_flux_err, _ = bin_array(phots_err, binsize=binsize)\n",
    "    \n",
    "    bin_time    , _  = bin_array(times, binsize=binsize)\n",
    "    \n",
    "    bin_flux_err     = bin_flux_err/sqrt(binsize)\n",
    "    \n",
    "    if ax is None:\n",
    "        fig= figure() if figsize is None else figure(figsize=figsize)\n",
    "        ax = fig.add_subplot(111)\n",
    "    \n",
    "    if plotRawData:\n",
    "        ax.plot(times, phots/ instrumental, '.', ms=1, label='data', alpha=alpha)\n",
    "    \n",
    "    ingress   = np.where(in_eclipse)[0].min()\n",
    "    egress    = np.where(in_eclipse)[0].max()\n",
    "    IT_mean   = np.mean([clean_model[ingress:egress]])\n",
    "    gap       = IT_mean - 1.0\n",
    "    \n",
    "    print('Gap: {:0f}'.format(gap*ppm))\n",
    "    \n",
    "    ax.errorbar(bin_time-refTime, bin_flux    - gap, bin_flux_err, color=color_cycle[0], \\\n",
    "                    fmt='o', label='binned data', ms=1, zorder=0, lw=1,alpha=alpha)\n",
    "    ax.plot(times-refTime, clean_model - gap, color=color_cycle[1], \\\n",
    "                lw=2,label='best fit model', zorder=1)\n",
    "    \n",
    "    ax.legend(loc=0)\n",
    "    \n",
    "    # ax.axhline(1.0, ls='--', lw=0.5)\n",
    "    \n",
    "    if returnAx:\n",
    "        return ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_model_over_reduced_lc_full_PLDsq_universal(times_0, times_1, phots_0, phots_1, phots_err_0, phots_err_1, \n",
    "                                                    input_features_0, input_features_1, in_eclipse_0, in_eclipse_1, \n",
    "                                                    fit_values, refTime=0.0, planetName='', figsize=None, \n",
    "                                                    nbins=200, alpha=1.0, plotRawData=False, \n",
    "                                                    phase_method='fourier', baseline='PLD', \n",
    "                                                    ax = None, returnAx = False, model_color=0, data_color=1):\n",
    "    \n",
    "    none_to_inf  = lambda x, sign=1: sign*np.inf if x is None else x\n",
    "    \n",
    "    _params_names = ['period', 'tCenter', 'inc', 'aprs', 'edepth', 'tdepth', 'ecc', 'omega', 'u1', 'u2', \n",
    "                     'amp1', 'amp2', 'amp3', 'amp4', 'amp5', 'amp6', 'amp7', 'amp8', # 'mean', \n",
    "                     'pld1_l', 'pld2_l', 'pld3_l', 'pld4_l', 'pld5_l', 'pld6_l', 'pld7_l', 'pld8_l', 'pld9_l', \n",
    "                     'pld1_q', 'pld2_q', 'pld3_q', 'pld4_q', 'pld5_q', 'pld6_q', 'pld7_q', 'pld8_q', 'pld9_q']#, \n",
    "                     # 'intcpt', 'slope', 'crvtur']\n",
    "    \n",
    "    fit_values_0 = Parameters()\n",
    "    fit_values_1 = Parameters()\n",
    "    \n",
    "    for pname in _params_names:\n",
    "        if 'pld' not in pname:\n",
    "            valu = fit_values[pname].value\n",
    "            vary = fit_values[pname].vary\n",
    "            rmin = none_to_inf(fit_values[pname].min, -1)\n",
    "            rmax = none_to_inf(fit_values[pname].max,  1)\n",
    "            \n",
    "            fit_values_0.add(pname, valu, vary, rmin, rmax)\n",
    "            fit_values_1.add(pname, valu, vary, rmin, rmax)\n",
    "        else:\n",
    "            pname_0 = pname + '_1'\n",
    "            valu = fit_values[pname_0].value\n",
    "            vary = fit_values[pname_0].vary\n",
    "            rmin = none_to_inf(fit_values[pname_0].min, -1)\n",
    "            rmax = none_to_inf(fit_values[pname_0].max,  1)\n",
    "            \n",
    "            fit_values_0.add(pname, valu, vary, rmin, rmax)\n",
    "            \n",
    "            pname_1 = pname + '_2'\n",
    "            valu = fit_values[pname_1].value\n",
    "            vary = fit_values[pname_1].vary\n",
    "            rmin = none_to_inf(fit_values[pname_1].min, -1)\n",
    "            rmax = none_to_inf(fit_values[pname_1].max,  1)\n",
    "            fit_values_1.add(pname, valu, vary, rmin, rmax)\n",
    "    \n",
    "    fit_values_0['mean'] = fit_values['mean_0']\n",
    "    fit_values_1['mean'] = fit_values['mean_1']\n",
    "    \n",
    "    # nPeriods = int((times_0.mean() - fit_values_0['tCenter'])/fit_values_0['period'])\n",
    "    \n",
    "    ax_0 = plot_model_over_reduced_lc_half_PLDsq_universal(times_0, phots_0, phots_err_0,\n",
    "                                                           input_features_0, in_eclipse_0, \n",
    "                                                           fit_values_0, refTime=refTime, \n",
    "                                                           figsize=figsize, nbins=nbins, plotRawData=plotRawData,\n",
    "                                                           phase_method=phase_method, baseline=baseline,\n",
    "                                                           ax = ax, returnAx = True, alpha=alpha, \n",
    "                                                           model_color=model_color, data_color=data_color)\n",
    "    \n",
    "    ax_0.set_title(planetName + ' E-depth: {:.0f}'.format(fit_values['edepth'].value*ppm))\n",
    "    \n",
    "    return plot_model_over_reduced_lc_half_PLDsq_universal(times_1, phots_1, phots_err_1,\n",
    "                                                           input_features_1, in_eclipse_1, \n",
    "                                                           fit_values_1, refTime=refTime, alpha=alpha,  \n",
    "                                                           figsize=figsize, nbins=nbins, plotRawData=plotRawData,\n",
    "                                                           phase_method=phase_method, baseline=baseline, \n",
    "                                                           ax = ax_0, returnAx = returnAx,\n",
    "                                                           model_color=model_color, data_color=data_color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**EMCEE Setup**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_prior(params):\n",
    "    \"\"\"\n",
    "    emccee uses a uniform prior for every variable.\n",
    "    Here we create a functions which checks the bounds\n",
    "    and returns np.inf if a value is outside of its\n",
    "    allowed range. WARNING: A uniform prior may not be\n",
    "    what you want!\n",
    "    \"\"\"\n",
    "    none_to_inf  = lambda x, sign=1: sign*np.inf if x is None else x\n",
    "    lower_bounds = np.array([none_to_inf(i.min, -1) for i in params.values() if i.vary])\n",
    "    upper_bounds = np.array([none_to_inf(i.max,  1) for i in params.values() if i.vary])\n",
    "    \n",
    "    def bounds_prior(values):\n",
    "        values = np.asarray(values)\n",
    "        is_ok = np.all((lower_bounds < values) & (values < upper_bounds))\n",
    "        return 0 if is_ok else -np.inf\n",
    "    \n",
    "    return bounds_prior\n",
    "\n",
    "def create_lnliklihood(mini, sigma=None):\n",
    "    \"\"\"create a normal-likihood from the residuals\"\"\"\n",
    "    def lnprob(vals, sigma=sigma):\n",
    "        for v, p in zip(vals, [p for p in mini.params.values() if p.vary]):\n",
    "            p.value = v\n",
    "        \n",
    "        # residuals = mini.residual\n",
    "        \n",
    "        if not sigma:\n",
    "            # sigma is either the RMS estimate or it will\n",
    "            # be part of the sampling.\n",
    "            sigma = vals[-1]\n",
    "        \n",
    "        val = -0.5*np.sum(np.log(2*np.pi*sigma**2) + (partial_residuals(mini.params)/sigma)**2)\n",
    "        \n",
    "        return val\n",
    "    \n",
    "    return lnprob\n",
    "\n",
    "def starting_guess(mini, estimate_sigma=True):\n",
    "    \"\"\"\n",
    "    Use best a fit as a starting point for the samplers.\n",
    "    If no sigmas are given, it is assumed that\n",
    "    all points have the same uncertainty which will\n",
    "    be also part of the sampled parameters.\n",
    "    \"\"\"\n",
    "    vals = [i.value for i in mini.params.values() if i.vary]\n",
    "    \n",
    "    if estimate_sigma:\n",
    "        vals.append(mini.residual.std())\n",
    "    \n",
    "    return vals\n",
    "\n",
    "def create_all(mini, sigma=None):\n",
    "    \"\"\"\n",
    "    creates the log-poposterior function from a minimizer.\n",
    "    sigma should is either None or an array with the\n",
    "    1-sigma uncertainties of each residual point. If None,\n",
    "    sigma will be assumed the same for all residuals and\n",
    "    is added to the sampled parameters.\n",
    "    \"\"\"\n",
    "    sigma_given = not sigma is None\n",
    "    \n",
    "    lnprior = create_prior(mini.params)\n",
    "    lnprob  = create_lnliklihood(mini, sigma=sigma)\n",
    "    guess   = starting_guess(mini, not sigma_given)\n",
    "    \n",
    "    if sigma_given:\n",
    "        func = lambda x: lnprior(x[:]) + lnprob(x)\n",
    "    else:\n",
    "        func = lambda x: lnprior(x[:-1]) + lnprob(x)\n",
    "    \n",
    "    return func, guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_selector(model_params, times, input_features, in_eclipse, \n",
    "                   spider_params=None, phase_method='fourier', \n",
    "                   baseline = 'PLD', simple=False, model_only=False):\n",
    "    \"\"\"Identifies the phase curve + noise model to use from `method`, allocates \n",
    "        then model specific params and returns the model\n",
    "    \n",
    "    Args:\n",
    "        model_params    : set of model parmaeters  -- named to avoid global / local confusion\n",
    "        times           : temporal array to use with model\n",
    "        input_features  : basis vectors for linear fitting to noise array\n",
    "        in_eclipse      : boolean array to identify the eclipse location from initial guess\n",
    "        method          : string to identify if residuals should be taken over \n",
    "                            `fourier` or `spiderman` phase curve models\n",
    "    \n",
    "    Returns:\n",
    "        Returns the current model, which is phase_curve x transit x eclipse x normalization\n",
    "    \"\"\"\n",
    "    \n",
    "    # Grab list of pld1..9 linear, then add list of pld1..9 quadratic\n",
    "    model_params_keys      = np.sort(list(model_params.keys()))\n",
    "    \n",
    "    feature_coeffs  = [model_params[key].value for key in model_params_keys if 'pld' in key and '_l' in key] + \\\n",
    "                      [model_params[key].value for key in model_params_keys if 'pld' in key and '_q' in key]\n",
    "    \n",
    "    # Check if line fit parameters exist, then allocate them; or use defaults\n",
    "    intcpt  = model_params['intcpt'].value if 'intcpt' in model_params.keys() else 1.0\n",
    "    slope   = model_params['slope'].value  if 'slope'  in model_params.keys() else 0.0\n",
    "    crvtur  = model_params['crvtur'].value if 'crvtur' in model_params.keys() else 0.0\n",
    "    \n",
    "    # Transit Parameters\n",
    "    period  = model_params['period'].value\n",
    "    tCenter = model_params['tCenter'].value\n",
    "    inc     = model_params['inc'].value\n",
    "    aprs    = model_params['aprs'].value\n",
    "    edepth  = model_params['edepth'].value\n",
    "    tdepth  = model_params['tdepth'].value\n",
    "    ecc     = model_params['ecc'].value\n",
    "    omega   = model_params['omega'].value\n",
    "    u1      = model_params['u1'].value\n",
    "    u2      = model_params['u2'].value\n",
    "    \n",
    "    delta_phase = deltaphase_eclipse(ecc, omega) if ecc is not 0.0 else 0.5\n",
    "    t_secondary = tCenter + period*delta_phase\n",
    "    \n",
    "    # Create phase curve model -- including PLD noise parameters\n",
    "    if phase_method == 'spiderman':\n",
    "        # Siderman Phase Curve Parameters\n",
    "        xi      = model_params['xi'].value\n",
    "        T_night = model_params['T_night'].value\n",
    "        delta_T = model_params['delta_T'].value\n",
    "        mean    = model_params['mean'].value\n",
    "        \n",
    "        spider_model  = partial(spiderman_lmfit_model, period=period, tCenter=tCenter, inc=inc, aprs=aprs, \n",
    "                                           tdepth=tdepth, ecc=ecc, omega=omega, u1=u1, u2=u2, \n",
    "                                           spider_params=spider_params, \n",
    "                                           xi=xi, T_night=T_night, delta_T=delta_T, if_eclipse=False)\n",
    "        \n",
    "        phase_curve   = spider_model(times=times) + mean\n",
    "    \n",
    "    if phase_method == 'fourier':\n",
    "        mean = model_params['mean'].value\n",
    "        amp1 = model_params['amp1'].value\n",
    "        amp2 = model_params['amp2'].value\n",
    "        amp3 = model_params['amp3'].value\n",
    "        amp4 = model_params['amp4'].value\n",
    "        # amp5 = model_params['amp5'].value\n",
    "        # amp6 = model_params['amp6'].value\n",
    "        # amp7 = model_params['amp7'].value\n",
    "        # amp8 = model_params['amp8'].value\n",
    "        \n",
    "        angphase       = 2 * pi / period * (times-t_secondary)\n",
    "        \n",
    "        phase_curve   = mean + amp1*cos(angphase)   \\\n",
    "                             + amp2*sin(angphase)   \\\n",
    "                             + amp3*cos(2*angphase) \\\n",
    "                             + amp4*sin(2*angphase) # \\\n",
    "                             # + amp5*cos(4*angphase) \\\n",
    "                             # + amp6*sin(4*angphase) \\\n",
    "                             # + amp7*cos(6*angphase) \\\n",
    "                             # + amp8*sin(6*angphase)  \n",
    "    \n",
    "    if phase_method == 'kbs':\n",
    "        mean = model_params['mean'].value\n",
    "        amp1 = model_params['amp1'].value\n",
    "        amp2 = model_params['amp2'].value\n",
    "        \n",
    "        angphase      = 2 * pi / period * (times - t_secondary - amp2)\n",
    "        phase_curve   = mean + amp1*cos(angphase)\n",
    "    \n",
    "    instrumental = intcpt * np.ones(times.size)\n",
    "    if slope  != 0.0:\n",
    "        instrumental += slope * (times-times.mean())\n",
    "    if crvtur != 0.0:\n",
    "        instrumental += crvtur * (times-times.mean())**2.\n",
    "    \n",
    "    batman_model  = partial(batman_lmfit_model, period=period, tCenter=tCenter, inc=inc, aprs=aprs  , \n",
    "                                       times=times, edepth=edepth, tdepth=tdepth  , ecc=ecc, omega=omega)\n",
    "    \n",
    "    transit       = batman_model(transittype=\"primary\"  , ldtype='quadratic', u1=u1, u2=u2)\n",
    "    eclipse       = batman_model(transittype=\"secondary\", ldtype='uniform') - edepth\n",
    "    \n",
    "    # eclipse_1     = (eclipse == 1 - edepth)#*()\n",
    "    # eclipse_2     = (eclipse == 1 - edepth)#*()\n",
    "    # in_eclipse    = eclipse != 1# - edepth\n",
    "    #if in_eclipse.any():\n",
    "    #    # midpoint = int(np.where(in_eclipse)[0].mean())\n",
    "    #    phase_curve[eclipse_1] = 1.0#phase_curve[in_eclipse]#np.mean([phase_curve[in_eclipse]])#eclipse[in_eclipse]#\n",
    "    \n",
    "    phase_curve[eclipse == 1.0 - edepth] = 1.0#np.mean(phase_curve[eclipse == 1.0 - edepth])\n",
    "    clean_model = transit * eclipse * phase_curve\n",
    "    \n",
    "    if baseline == 'PLD':\n",
    "        instrumental += pixel_level_decorrelation_instrument_profile(times, input_features, \n",
    "                        *feature_coeffs)#, intcpt=intcpt, slope=slope, crvtur=crvtur)\n",
    "    \n",
    "    if baseline == 'KRData':\n",
    "        phots_now, gk, nbr_ind = input_features\n",
    "        instrumental *= np.sum((phots_now/clean_model)[nbr_ind]*gk,axis=1)\n",
    "    \n",
    "    if model_only:\n",
    "        return clean_model\n",
    "    else:\n",
    "        return clean_model * instrumental"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def residuals_func(model_params, data, times, input_features, in_eclipse, simple=True, \n",
    "                   spider_params=None, phase_method='fourier', baseline='PLD'):\n",
    "    \"\"\"Returns residuals of current model estimate with input `data`\n",
    "    \n",
    "    Args:\n",
    "        model_params    : set of model parmaeters  -- named to avoid global / local confusion\n",
    "        data            : data to be fit\n",
    "        times           : temporal array to use with model\n",
    "        input_features  : basis vectors for linear fitting to noise array\n",
    "        in_eclipse      : boolean array to identify the eclipse location from initial guess\n",
    "        method          : string to identify if residuals should be taken over \n",
    "                            `fourier` or `spiderman` phase curve models\n",
    "    \n",
    "    Returns:\n",
    "        Returns the array of residuals, that is the data - model\n",
    "    \"\"\"\n",
    "    # Return residuals between this model and the `data` allocated above\n",
    "    model_now = model_selector(model_params, times, input_features, in_eclipse, \n",
    "                               spider_params, phase_method, baseline, simple=simple)\n",
    "    return model_now - data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def load_aor_from_phot_pipeline(planet_dir_name, channel, AORNow, PLD_order=1):\n",
    "    loadfiledir         = environ['HOME']+'/Research/Planets/PhaseCurves/'+planet_dir_name+'/saveFiles/' + channel \n",
    "    loadFileNameHeader  = planet_dir_name+'_'+ AORNow +'_Median'\n",
    "    loadFileType        = '.pickle.save'\n",
    "    \n",
    "    print()\n",
    "    print('Loading from ' + loadfiledir + loadFileNameHeader + loadFileType)\n",
    "    print()\n",
    "    \n",
    "    instance = wanderer(fitsFileDir=loadfitsdir, filetype=filetype, telescope='Spitzer', \n",
    "                                                yguess=yguess, xguess=xguess, method='median', nCores=cpu_count())\n",
    "    \n",
    "    instance.load_data_from_save_files(savefiledir=loadfiledir, \\\n",
    "                                                      saveFileNameHeader=loadFileNameHeader,\\\n",
    "                                                      saveFileType='.pickle.save')\n",
    "    \n",
    "    if PLD_order > 1: instance.extract_PLD_components(order=PLD_order)\n",
    "    \n",
    "    timeCubeLocal       = instance.timeCube\n",
    "    photsLocal          = instance.flux_TSO_df.values\n",
    "    PLDFeatureLocal     = instance.PLD_components.T\n",
    "    \n",
    "    try:\n",
    "        inliers_Phots_local = instance.inliers_Phots.values()\n",
    "    except:\n",
    "        inliers_Phots_local = np.ones(photsLocal.shape)\n",
    "    \n",
    "    try:\n",
    "        inliers_PLD_local = instance.inliers_PLD.values()\n",
    "    except:\n",
    "        inliers_PLD_local = np.ones(PLDFeatureLocal.shape)\n",
    "    \n",
    "    return timeCubeLocal, photsLocal, PLDFeatureLocal, inliers_Phots_local, inliers_PLD_local, instance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Load Planet Params**"
   ]
  },
  {
   "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": [
    "**Starting Position**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "planet_params_name = 'Planet-Name b'# 'HAT-P-23 b'# \n",
    "planet_dir_name    = 'planetName'# 'hat23'# \n",
    "channel = 'ch1/'\n",
    "\n",
    "if   'ch1' in channel:\n",
    "    spitzer_waves_low  = 3.0\n",
    "    spitzer_waves_high = 4.0\n",
    "elif 'ch2' in channel:\n",
    "    spitzer_waves_low  = 4.0\n",
    "    spitzer_waves_high = 5.0\n",
    "else:\n",
    "    raise Exception(\"`channel` must be either `'ch1'` or `'ch2'`\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.randint?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "planet_params = PlanetParams(planet_params_name)\n",
    "\n",
    "iPeriod   = planet_params.per.value\n",
    "iTCenter  = planet_params.tt.value-2400000.5\n",
    "# iBImpact  = planet_params.b.value\n",
    "iApRs     = planet_params.ar.value\n",
    "iInc      = planet_params.i.value\n",
    "# iRsAp     = 1.0/planet_params.ar.value\n",
    "iEdepth   = 3000/ppm # blind guess\n",
    "iTdepth   = planet_params.depth.value\n",
    "iEcc      = planet_params.ecc.value\n",
    "iOmega    = planet_params.om.value*pi/180\n",
    "\n",
    "stellar_radius = planet_params.rstar.value\n",
    "stellar_temp   = planet_params.teff.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pixel Level Decorrelation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "planetDirectory = '/Research/Planets/PhaseCurves/'\n",
    "\n",
    "dataSub = 'bcd/'\n",
    "\n",
    "dataDir     = environ['HOME'] + planetDirectory + planet_dir_name + '/data/raw/' + channel + '/big/'\n",
    "print(dataDir + '*')\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",
    "iAOR        = 1\n",
    "AORNow      = AORs[iAOR]\n",
    "loadfitsdir = dataDir + AORNow + '/' + channel + dataSub\n",
    "print(loadfitsdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeCubeStack       = {}\n",
    "photsStack          = {}\n",
    "PLDFeatureStack     = {}\n",
    "inliers_PhotsStack  = {}\n",
    "inliers_PLDStack    = {}\n",
    "instanceStack       = {}\n",
    "\n",
    "for AORNow in AORs:\n",
    "    timesNow, photsNow, PLDFeaturesNow, inliers_PhotsNow, inliers_PLDNow, instanceNow = \\\n",
    "        load_aor_from_phot_pipeline(planet_dir_name, channel, AORNow, PLD_order=2)\n",
    "                                \n",
    "    timeCubeStack[AORNow]       = timesNow\n",
    "    photsStack[AORNow]          = photsNow\n",
    "    PLDFeatureStack[AORNow]     = PLDFeaturesNow\n",
    "    inliers_PhotsStack0[AORNow] = inliers_PhotsNow\n",
    "    inliers_PLDStack[AORNow]    = inliers_PLDNow\n",
    "    instanceStack[AORNow]       = instanceNow"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for aor_now in instanceStack.keys():\n",
    "    print('Working on AOR {}'.format(aorNow))\n",
    "    instanceStack[aor_now].mp_lmfit_gaussian_centering()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for aor_now in instanceStack.keys():\n",
    "    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/PhaseCurves/'+planet_dir_name+'/SaveFiles/' + channel \n",
    "    saveFileNameHeader  = planet_dir_name+'_'+ aor_now +'_Median'\n",
    "    saveFileType        = '.pickle.save'\n",
    "    \n",
    "    if not path.exists(environ['HOME']+'/Research/Planets/PhaseCurves/'+planet_dir_name+'/SaveFiles/'):\n",
    "        mkdir(environ['HOME']+'/Research/Planets/PhaseCurves/'+planet_dir_name+'/SaveFiles/')\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",
    "    instanceStack[aor_now].save_data_to_save_files(savefiledir=savefiledir, \\\n",
    "                                                    saveFileNameHeader=saveFileNameHeader, \\\n",
    "                                                    saveFileType=saveFileType)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gaussian_Fit_AnnularMask_rad_1.0_0.0\n",
    "\n",
    "staticRad   = '2.0'\n",
    "varRad      = '0.0'\n",
    "phot_select = np.where(instance_0.flux_TSO_df.keys() == 'Gaussian_Fit_AnnularMask_rad_'+staticRad+'_'+varRad)[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    aor_0, aor_1 = timeCubeStack.keys();\n",
    "except:\n",
    "    aor_0, aor_1, aor_2 = timeCubeStack.keys();\n",
    "    print('NEED TO ADJUST FOR 3 AORs')\n",
    "\n",
    "if timeCubeStack[aor_0].min() > timeCubeStack[aor_1].min():\n",
    "    print('Reversing AOR Labels to be Time Ordered')\n",
    "    aor_0, aor_1 = aor_1, aor_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inliersMaster = {}\n",
    "inliersMaster[aor_0]  = array(list(inliers_PhotsStack0[aor_0])).all(axis=0) # Need to Switch `axis=0` for Planet-Name\n",
    "inliersMaster[aor_1]  = array(list(inliers_PhotsStack0[aor_1])).all(axis=0) # Need to Switch `axis=0` for Planet-Name\n",
    "\n",
    "inliersMaster[aor_0]  = inliersMaster[aor_0] * inliers_PLDStack[aor_0].all(axis=1)\n",
    "inliersMaster[aor_1]  = inliersMaster[aor_1] * inliers_PLDStack[aor_1].all(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "instance_0  = instanceStack[aor_0]\n",
    "instance_1  = instanceStack[aor_1]\n",
    "\n",
    "instance_0.__name__ = aor_0\n",
    "instance_1.__name__ = aor_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nSig=6\n",
    "for aorNow in inliersMaster.keys():\n",
    "    if inliersMaster[aorNow].all():\n",
    "        print('Working on AOR {}'.format(aorNow))\n",
    "        cy_now, cx_now        = instanceStack[aorNow].centering_GaussianFit.T\n",
    "        phots_now             = photsStack[aorNow][:,phot_select]\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",
    "        inliersMaster[aorNow] = (phots_clipped == phots_now)*(cy_clipped==cy_now)*(cx_clipped==cx_now)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "times_0 = timeCubeStack[aor_0]#[inliersMaster[aor_0]]\n",
    "times_1 = timeCubeStack[aor_1]#[inliersMaster[aor_1]]\n",
    "\n",
    "phots_0 = photsStack[aor_0][:, phot_select]#[inliersMaster[aor_0], phot_select]\n",
    "phots_1 = photsStack[aor_1][:, phot_select]#[inliersMaster[aor_1], phot_select]\n",
    "\n",
    "PLDfeatures_0 = PLDFeatureStack[aor_0]#[inliersMaster[aor_0]]\n",
    "PLDfeatures_1 = PLDFeatureStack[aor_1]#[inliersMaster[aor_1]]\n",
    "\n",
    "phots_0_err = np.sqrt(abs(phots_0))\n",
    "phots_1_err = np.sqrt(abs(phots_1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_centers_stack = hstack([instance_0.centering_GaussianFit.T[y][inliersMaster[aor_0]], \\\n",
    "                          instance_1.centering_GaussianFit.T[y][inliersMaster[aor_1]]])\n",
    "x_centers_stack = hstack([instance_0.centering_GaussianFit.T[x][inliersMaster[aor_0]], \\\n",
    "                          instance_1.centering_GaussianFit.T[x][inliersMaster[aor_1]]])\n",
    "\n",
    "# planetName_range_ch2 = [[15.65,15.95],[15.85,16.15]]\n",
    "# planetName_range_ch1 = [[15.25,16.0],[15.95,16.2]]\n",
    "wasp77_range_ch1 = [[14.9,15.25],[14.95,15.3]]\n",
    "plt.hist2d(x_centers_stack, y_centers_stack, bins=50, cmap=plt.cm.jet);#, range=wasp77_range_ch1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.plot(*(transpose([instance_0.centering_GaussianFit.T[x],instance_0.centering_GaussianFit.T[y]])).T,'.',ms=1,alpha=0.25);\n",
    "plt.plot(*(transpose([instance_0.centering_GaussianFit.T[x][inliersMaster[aor_0]],\\\n",
    "                      instance_0.centering_GaussianFit.T[y][inliersMaster[aor_0]]])).T,'.',ms=1,alpha=0.25);\n",
    "plt.plot(*(transpose([instance_1.centering_GaussianFit.T[x][inliersMaster[aor_1]],\\\n",
    "                      instance_1.centering_GaussianFit.T[y][inliersMaster[aor_1]]])).T,'.',ms=1,alpha=0.25);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Setup Initial Parameters**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_eclipse_params = TransitParams()\n",
    "\n",
    "init_eclipse_params.per         = iPeriod\n",
    "init_eclipse_params.t0          = iTCenter#+0.0025\n",
    "init_eclipse_params.inc         = iInc\n",
    "init_eclipse_params.a           = iApRs\n",
    "init_eclipse_params.fp          = iEdepth\n",
    "init_eclipse_params.rp          = sqrt(iTdepth)\n",
    "init_eclipse_params.ecc         = iEcc\n",
    "init_eclipse_params.w           = iOmega\n",
    "\n",
    "init_eclipse_params.delta_phase = deltaphase_eclipse(init_eclipse_params.ecc, init_eclipse_params.w) \\\n",
    "                                    if init_eclipse_params.ecc is not 0.0 else 0.5\n",
    "init_eclipse_params.t_secondary = init_eclipse_params.t0 + init_eclipse_params.per*init_eclipse_params.delta_phase\n",
    "\n",
    "# init_eclipse_params.t_secondary = init_eclipse_params.t0 + iPeriod * 0.5\n",
    "init_eclipse_params.limb_dark   = 'uniform'\n",
    "init_eclipse_params.u           = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_1st_eclipse_model          = TransitModel(init_eclipse_params, times_0, transittype=\"secondary\")\n",
    "init_1st_eclipse_lightcurve     = init_1st_eclipse_model.light_curve(init_eclipse_params)\n",
    "\n",
    "init_1st_transit_model          = TransitModel(init_eclipse_params, times_0, transittype=\"primary\")\n",
    "init_1st_transit_lightcurve     = init_1st_transit_model.light_curve(init_eclipse_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "in_1st_eclipse = init_1st_eclipse_lightcurve == 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(in_1st_eclipse == (init_1st_eclipse_lightcurve == init_1st_transit_lightcurve)).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skip = 1000\n",
    "plot(photsStack[aor_0][:,phot_select]);\n",
    "plot(np.arange(phots_0.shape[0])[skip:], phots_0[skip:],'.',ms=2);\n",
    "plot(np.arange(phots_0.shape[0])[:skip], phots_0[:skip],'.',ms=2);\n",
    "\n",
    "med_phots_0 = median(phots_0[inliersMaster[aor_0]])\n",
    "std_phots_0 = std(phots_0[inliersMaster[aor_0]])\n",
    "\n",
    "plot(np.arange(init_1st_eclipse_lightcurve.size), init_1st_transit_lightcurve*med_phots_0);\n",
    "plot(np.arange(init_1st_eclipse_lightcurve.size), init_1st_eclipse_lightcurve*med_phots_0);\n",
    "plot(np.arange(init_1st_eclipse_lightcurve.size)[in_1st_eclipse], \\\n",
    "     init_1st_eclipse_lightcurve[in_1st_eclipse]*med_phots_0);\n",
    "\n",
    "ylim(med_phots_0 - 5*std_phots_0, med_phots_0 + 5*std_phots_0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Choose to Skip First ~30 Minutes or Not**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SKIP_30_MIN = False\n",
    "if SKIP_30_MIN:\n",
    "    nskip   = 1000 # with 2s integrations, 1000 frames ~ 33 minutes\n",
    "    inliersMaster[aor_0][:nskip] = False"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "pca = PCA()\n",
    "ica = FastICA()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "pca_pld_0         = pca.fit_transform(PLDfeatures_0)\n",
    "ica_pld_0         = ica.fit_transform(PLDfeatures_0)\n",
    "\n",
    "pca_scl_pld_0     = pca.fit_transform(scale(PLDfeatures_0))\n",
    "ica_scl_pld_0     = ica.fit_transform(scale(PLDfeatures_0))\n",
    "\n",
    "ica_pca_pld_0     = ica.fit_transform(pca_pld_0)\n",
    "pca_ica_pld_0     = pca.fit_transform(ica_pld_0)\n",
    "\n",
    "ica_pca_scl_pld_0 = ica.fit_transform(pca_scl_pld_0)\n",
    "pca_ica_scl_pld_0 = pca.fit_transform(ica_scl_pld_0)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "input_feature_set_0   = PLDfeatures_0\n",
    "# input_feature_set_0   = scale(PLDfeatures_0)\n",
    "# input_feature_set_0   = pca_pld_0\n",
    "# input_feature_set_0   = ica_pld_0\n",
    "# input_feature_set_0   = pca_scl_pld_0\n",
    "# input_feature_set_0   = ica_scl_pld_0\n",
    "# input_feature_set_0   = ica_pca_pld_0\n",
    "# input_feature_set_0   = pca_ica_pld_0\n",
    "# input_feature_set_0   = ica_pca_scl_pld_0\n",
    "# input_feature_set_0   = pca_ica_scl_pld_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ypos_0, xpos_0  = clipOutlier2D(transpose([instance_0.centering_GaussianFit.T[y][inliersMaster[aor_0]], \\\n",
    "                                           instance_0.centering_GaussianFit.T[x][inliersMaster[aor_0]]])).T\n",
    "                                          # instance_0.centering_GaussianFit[inliersMaster[aor_0]].T\n",
    "ypos_1, xpos_1  = clipOutlier2D(transpose([instance_1.centering_GaussianFit.T[y][inliersMaster[aor_1]], \\\n",
    "                                           instance_1.centering_GaussianFit.T[x][inliersMaster[aor_1]]])).T\n",
    "                                          # instance_1.centering_GaussianFit[inliersMaster[aor_1]].T\n",
    "\n",
    "npix_0          = sqrt(instance_0.effective_widths[inliersMaster[aor_0]])\n",
    "npix_1          = sqrt(instance_1.effective_widths[inliersMaster[aor_1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(figsize=(30,10));\n",
    "ax1 = fig.add_subplot(1,3,1);\n",
    "ax2 = fig.add_subplot(1,3,2);\n",
    "ax3 = fig.add_subplot(1,3,3);\n",
    "ax1.plot(xpos_0, ypos_0,'.',ms=1,alpha=0.25);\n",
    "ax2.plot(xpos_0, npix_0,'.',ms=1,alpha=0.25);\n",
    "ax3.plot(ypos_0, npix_0,'.',ms=1,alpha=0.25);\n",
    "ax1.plot(xpos_1, ypos_1,'.',ms=1,alpha=0.25);\n",
    "ax2.plot(xpos_1, npix_1,'.',ms=1,alpha=0.25);\n",
    "ax3.plot(ypos_1, npix_1,'.',ms=1,alpha=0.25);"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "start = time()\n",
    "gw_0, nbr_0 = mp_find_nbr_qhull(xpos_0, ypos_0, npix_0, \n",
    "                                sm_num = 100, a = 1.0, b = 0.7, c = 1.0, \n",
    "                                print_space = 10000., nCores=cpu_count())\n",
    "print(time() - start)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "start = time()\n",
    "gw_1, nbr_1 = mp_find_nbr_qhull(xpos_1, ypos_1, npix_1, \n",
    "                                sm_num = 100, a = 1.0, b = 0.7, c = 1.0, \n",
    "                                print_space = 10000., nCores=cpu_count())\n",
    "print(time() - start)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "gw_0, nbr_0 = find_nbr_qhull(xpos_0, ypos_0, npix_0, sm_num = 100, a = 1.0, b = 0.7, c = 1.0, print_space = 10000.)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "gw_1, nbr_1 = find_nbr_qhull(xpos_1, ypos_1, npix_1, sm_num = 100, a = 1.0, b = 0.7, c = 1.0, print_space = 10000.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xpos_c  = np.hstack([xpos_0, xpos_1])\n",
    "ypos_c  = np.hstack([ypos_0, ypos_1])\n",
    "npix_c  = np.hstack([npix_0, npix_1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = time()\n",
    "gw_c, nbr_c = mp_find_nbr_qhull(xpos_c, ypos_c, npix_c, sm_num = 100, a = 1.0, b = 0.7, c = 1.0, print_space = 10000.)\n",
    "print(time() - start)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "gw_c, nbr_c = find_nbr_qhull(xpos_c, ypos_c, npix_c, sm_num = 100, a = 1.0, b = 0.7, c = 1.0, print_space = 10000.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_cycle = cycler(rcParams['axes.prop_cycle']).by_key()['color']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "med_eclipse_1 = median(phots_0[in_1st_eclipse == 1.0])"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "coeffs = pywt.wavedec(phots_0 / med_eclipse_1, 'db4', level=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i_u1, i_u2 = 0.1, 0.1\n",
    "xi_0, T_night_0, delta_T_0 = 1e-5, 750.0, 500.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**SPIDERMAN Fitting Procedure**\n",
    "**Fit Spiderman Model with PLD to 1st AOR**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spider_params_hoststar                = sp_ModelParams(brightness_model='zhang')\n",
    "spider_params_hoststar.n_layers       = 20\n",
    "spider_params_hoststar.stellar_radius = stellar_radius\n",
    "spider_params_hoststar.T_s            = stellar_temp\n",
    "spider_params_hoststar.l1             = spitzer_waves_low  # Spitzer IRAC-1 or IRAC-2\n",
    "spider_params_hoststar.l2             = spitzer_waves_high # Spitzer IRAC-1 or IRAC-2\n",
    "\n",
    "stellar_a_au  = iApRs * stellar_radius * R_sun.value / au.value\n",
    "\n",
    "spider_params_hoststar.t0      = iTCenter       # Central time of PRIMARY transit [days]\n",
    "spider_params_hoststar.per     = iPeriod        # Period [days]\n",
    "spider_params_hoststar.a_abs   = stellar_a_au   # The absolute value of the semi-major axis [AU]\n",
    "spider_params_hoststar.inc     = iInc          # orbital inclination (in degrees)\n",
    "spider_params_hoststar.ecc     = iEcc          # Eccentricity\n",
    "spider_params_hoststar.w       = iOmega        # Argument of periastron\n",
    "spider_params_hoststar.rp      = sqrt(iTdepth) # planet radius (in units of stellar radii)\n",
    "spider_params_hoststar.a       = iApRs         # Semi-major axis scaled by stellar radius\n",
    "spider_params_hoststar.p_u1    = i_u1          # Planetary limb darkening parameter\n",
    "spider_params_hoststar.p_u2    = i_u2          # Planetary limb darkening parameter\n",
    "\n",
    "spider_params_hoststar.eclipse = False\n",
    "spider_params_hoststar.xi      = 0.0           # Ratio of radiative to advective timescale             \n",
    "spider_params_hoststar.T_n     = 1000.0        # Temperature of nightside\n",
    "spider_params_hoststar.delta_T = 500.0         # Day-night temperature contrast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Global Fitting Procedure**"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "initialParams_KRData_half_Fourier_PhaseCurve_ch2 = initialParams_KRData_half_Fourier_PhaseCurve\n",
    "del initialParams_KRData_half_Fourier_PhaseCurve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i_u1, i_u2, iEdepth = 0.1, 0.0, 2000/ppm\n",
    "xi_0, T_night_0, delta_T_0 = 1e-5, 750.0, 500.0\n",
    "\n",
    "initialParams_KRData_half_Fourier_PhaseCurve = Parameters()\n",
    "\n",
    "# c1a, c1o, c = 7.8915e-04, 2.9006e-02, 1.0000e+00\n",
    "\n",
    "initialParams_KRData_half_Fourier_PhaseCurve.add_many(\n",
    "    # Planetary Parameters\n",
    "    ('period' , iPeriod  , False),\n",
    "    ('tCenter', iTCenter+0.00565 , True , iTCenter-0.01, iTCenter+0.02),\n",
    "    ('inc'    , iInc     , False, 0.0 ,  90.),\n",
    "    ('aprs'   , iApRs    , False, 0.0 , 100.),\n",
    "    ('edepth' , iEdepth  , True , 0.0 , 1.0 ),\n",
    "    ('tdepth' , iTdepth  , True , 0.0 , 1.0 ),\n",
    "    ('ecc'    , 0.0     , False , 0.0 , 1.0 ),\n",
    "    ('omega'  , iOmega   , False, 0.0 , 360.),\n",
    "    ('u1'     , i_u1     , False, 0.0 , 1.0 ),\n",
    "    ('u2'     , i_u2     , False, 0.0 , 1.0 ),\n",
    "    ('amp1'   , 1e-3     , True ),\n",
    "    ('amp2'   , 1e-3     , True ), # for non-KBS only\n",
    "    # ('amp2'   , 1e-4     , True , 0.0, 0.25*iPeriod ), # for KBS only\n",
    "    ('amp3'   , 1e-4    , True ),\n",
    "    ('amp4'   , 1e-4    , True ),\n",
    "    ('amp5'   , 0.0      , False),\n",
    "    ('amp6'   , 0.0      , False),\n",
    "    ('amp7'   , 0.0      , False),\n",
    "    ('amp8'   , 0.0      , False),\n",
    "    ('mean'   , 1.0      , False),\n",
    "    ('xi'     , xi_0     , False, 0.0, inf),\n",
    "    ('delta_T', delta_T_0, False, 0.0, inf),\n",
    "    ('T_night', T_night_0, False, 0.0, inf),\n",
    "    # Out of transit linear baselines\n",
    "    ('intcpt' , 1.0      , True ),\n",
    "    ('slope'  , 0.0      , True ),\n",
    "    ('crvtur' , 0.0      , False)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Test KR-Data on Whole System**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phots_c = np.hstack([phots_0[inliersMaster[aor_0]], phots_1[inliersMaster[aor_1]]])\n",
    "times_c = np.hstack([times_0[inliersMaster[aor_0]], times_1[inliersMaster[aor_1]]])\n",
    "# phots_c = np.hstack([phots_0, phots_1])\n",
    "# times_c = np.hstack([times_0, times_1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phase_method   = 'fourier'\n",
    "baseline       = 'KRData'\n",
    "partial_residuals = partial(residuals_func, \n",
    "                            data           = phots_c / median(phots_c),\n",
    "                            times          = times_c, \n",
    "                            input_features = [phots_c / median(phots_c), gw_c, nbr_c],#input_feature_set_c, \n",
    "                            in_eclipse     = in_1st_eclipse,\n",
    "                            spider_params  = spider_params_hoststar,\n",
    "                            phase_method   = phase_method,\n",
    "                            baseline       = baseline)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mini = Minimizer(partial_residuals, initialParams_KRData_half_Fourier_PhaseCurve)\n",
    "#fitResults_KRData_half_Fourier_PhaseCurve.params)#\n",
    "start = time()\n",
    "fitResults_KRData_half_Fourier_PhaseCurve = mini.leastsq()\n",
    "print(\"Full phase curve fitting operation took {} seconds\".format(time()-start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "report_errors(fitResults_KRData_half_Fourier_PhaseCurve.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_i_care = ['edepth', 'ecc', 'omega', 'amp1', 'amp2', 'amp3', 'amp4', 'amp5', 'amp6', 'amp7', 'amp8']\n",
    "ppm_or_not = [ppm, 1.0, 1.0, ppm, ppm, ppm,ppm, ppm, ppm, ppm, ppm]\n",
    "unit_list  = ['ppm', ' ', 'deg', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm']\n",
    "\n",
    "for varname in var_i_care:\n",
    "    if varname in fitResults_KRData_half_Fourier_PhaseCurve.var_names:\n",
    "        print('{}:\\t{:.0f}\\tppm'.format(varname, fitResults_KRData_half_Fourier_PhaseCurve.params[varname]*ppm))\n",
    "\n",
    "print('mean-1:\\t{:.0f}\\tppm'.format((fitResults_KRData_half_Fourier_PhaseCurve.params['mean']-1)*ppm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plot_model_over_reduced_lc_half_PLDsq_universal(times_c, \n",
    "                                                phots_c / median(phots_c), \n",
    "                                                sqrt(phots_c) / median(phots_c), \n",
    "                                                # input_feature_set_c, \n",
    "                                                [phots_c / median(phots_c), gw_c, nbr_c], \n",
    "                                                in_1st_eclipse, \n",
    "                                                # initialParams_KRData_half_Fourier_PhaseCurve, \n",
    "                                                fitResults_KRData_half_Fourier_PhaseCurve.params,\n",
    "                                                spider_params=spider_params_hoststar,\n",
    "                                                figsize=None, nbins=200, plotRawData=False,\n",
    "                                                phase_method=phase_method, baseline=baseline,\n",
    "                                                ax = None, returnAx = True)\n",
    "\n",
    "ax.set_title('{} - {} - {:.0f} ppm'.format(planet_params_name, channel[:-1].upper(), \\\n",
    "                                  fitResults_KRData_half_Fourier_PhaseCurve.params['edepth']*ppm));\n",
    "\n",
    "fig   = gcf()\n",
    "# baseline     = 'KRData'\n",
    "# phase_method = 'Fourier4'\n",
    "fig_save_name= planet_params_name.replace(' ','_') + '_full_phase_curve_' \\\n",
    "                                                   + channel.replace('/', '') \\\n",
    "                                                   + '_observations_' \\\n",
    "                                                   + phase_method+'_' \\\n",
    "                                                   + baseline + '.png'\n",
    "\n",
    "axhline(1.0, ls='--')\n",
    "# print('Saving tigure to '  + 'figure_results/' + fig_save_name)\n",
    "\n",
    "# fig = gcf()\n",
    "# fig.savefig('figure_results/' + fig_save_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plot_model_over_reduced_lc_half_PLDsq_universal(times_c, \n",
    "                                                phots_c / median(phots_c), \n",
    "                                                sqrt(phots_c) / median(phots_c), \n",
    "                                                # input_feature_set_c, \n",
    "                                                [phots_c / median(phots_c), gw_c, nbr_c], \n",
    "                                                in_1st_eclipse, \n",
    "                                                # initialParams_KRData_half_Fourier_PhaseCurve, \n",
    "                                                fitResults_KRData_half_Fourier_PhaseCurve.params,\n",
    "                                                spider_params=spider_params_hoststar,\n",
    "                                                figsize=None, nbins=200, plotRawData=False,\n",
    "                                                phase_method=phase_method, baseline=baseline,\n",
    "                                                ax = None, returnAx = True)\n",
    "\n",
    "ax.set_title('{} - {} - {:.0f} ppm'.format(planet_params_name, channel[:-1].upper(), \\\n",
    "                                  fitResults_KRData_half_Fourier_PhaseCurve.params['edepth']*ppm));\n",
    "\n",
    "ax.set_ylim(0.999,1.003)\n",
    "fig   = gcf()\n",
    "# baseline     = 'KRData'\n",
    "# phase_method = 'Fourier4'\n",
    "fig_save_name= planet_params_name.replace(' ','_') + '_top_full_phase_curve_' \\\n",
    "                                                   + channel.replace('/', '') \\\n",
    "                                                   + '_observations_' \\\n",
    "                                                   + phase_method+'_' \\\n",
    "                                                   + baseline + '.png'\n",
    "axhline(1.0, ls='--')\n",
    "\n",
    "\n",
    "# print('Saving tigure to '  + 'figure_results/' + fig_save_name)\n",
    "\n",
    "# fig = gcf()\n",
    "# fig.savefig('figure_results/' + fig_save_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "joblib.dump(fitResults_KRData_half_Fourier_PhaseCurve, 'lmfit_save_files/' + \\\n",
    "            planet_params_name.replace(' ','_') + '_full_phase_curve_' \\\n",
    "                                                + channel.replace('/', '') \\\n",
    "                                                + '_observations_' \\\n",
    "                                                + phase_method+'_' \\\n",
    "                                                + baseline + '.pickle.save')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls lmfit_save_files/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**END TEST**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create lnfunc and starting distribution.\n",
    "lnfunc, guess = create_all(fitResults_KRData_half_Fourier_PhaseCurve)\n",
    "nwalkers, ndim = 100, len(guess)\n",
    "p0 = emcee.utils.sample_ball(guess, 0.1*np.array(guess), nwalkers)\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, lnfunc)\n",
    "n_steps = 1\n",
    "\n",
    "start = time()\n",
    "sampler.run_mcmc(p0, n_steps)\n",
    "print(\"Full phase curve fitting operation took {} seconds\".format(time()-start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params      = fitResults_KRData_half_Fourier_PhaseCurve_ch1.params\n",
    "lmfit_vals  = [i.value for i in params.values() if i.vary] + [0.0]\n",
    "lmfit_names = [i.name  for i in params.values() if i.vary] + ['sigma']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(len(guess), 1, sharex=True, figsize=(8, 2*len(guess)))\n",
    "for (i, name, rv) in zip(range(len(guess)), lmfit_names, lmfit_vals):\n",
    "    # axes[i].hist(sampler.chain[:, :, i].T.ravel(),  alpha=0.05, bins=sampler.chain[:, :, i].T.size//100);\n",
    "    axes[i].plot(sampler.chain[:, :, i].T, color=\"k\", alpha=0.05);\n",
    "    axes[i].yaxis.set_major_locator(plt.MaxNLocator(5));\n",
    "    axes[i].axhline(rv, color=\"#888888\", lw=2);\n",
    "    axes[i].set_ylabel(\"$%s$\" % name);\n",
    "\n",
    "axes[-1].set_xlabel(\"Steps\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i_use, corner_vals, corner_names = np.transpose([(k, val, name) \\\n",
    "                                                 for k, (val, name) in enumerate(zip(lmfit_vals, lmfit_names)) \\\n",
    "                                                 if 'pld' not in name])\n",
    "i_use = int32(i_use)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(corner_names, len(corner_names))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corner_names=[r'tCenter', r'inc', r'aprs', r'edepth', r'tdepth', r'u1', r'u2',\n",
    "              r'amp1', r'amp2', r'amp3', r'amp4', r'mean', r'intcpt', r'slope', r'sigma']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples.shape, dims, len(i_use), len(corner_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "burnin  = int(n_steps * 0.2)\n",
    "samples = sampler.chain[:, burnin:, i_use].reshape((-1, len(i_use)));\n",
    "# corner.corner(samples, labels=corner_names, truths=corner_vals);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Stairstep Plot\n",
    "# labels = corner_names\n",
    "\n",
    "plt.rc('font',size=8)\n",
    "dims = len(corner_names)\n",
    "# fig,axL = plt.subplots(nrows=dims,ncols=dims,figsize=(15,15))\n",
    "\n",
    "corner_kw = dict(\n",
    "    truths=corner_vals,\n",
    "    labels=corner_names,\n",
    "    levels=[0.68,0.95],\n",
    "    plot_datapoints=False,\n",
    "    smooth=True,\n",
    "    bins=30,\n",
    "    )\n",
    "\n",
    "corner.corner(samples, **corner_kw)\n",
    "# plt.show()\n",
    "# fig.savefig('ch'+str(ch)+'_visit'+str(visit)+'_stairstep_indiv_fit_linreg'+exp_opt_name+'.eps')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# fit the data with lmfit\n",
    "partial_residuals = partial(residuals_func, \n",
    "                            data           = phots_0 / med_eclipse_1,\n",
    "                            times          = times_0, \n",
    "                            input_features = [phots_0 / med_eclipse_1, gw_0, nbr_0],#input_feature_set_0, \n",
    "                            in_eclipse     = in_1st_eclipse,\n",
    "                            phase_method   = 'fourier',\n",
    "                            baseline       = 'KRData')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "mini = Minimizer(partial_residuals, initialParams_1st_half_Fourier_PhaseCurve)\n",
    "start = time()\n",
    "fitResults_1st_half_Fourier_PhaseCurve = mini.leastsq()\n",
    "print(\"Full phase curve fitting operation took {} seconds\".format(time()-start))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "report_errors(fitResults_1st_half_Fourier_PhaseCurve.params)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "var_i_care = ['edepth', 'ecc', 'omega', 'amp1', 'amp2', 'amp3', 'amp4', 'amp5', 'amp6', 'amp7', 'amp8']\n",
    "ppm_or_not = [ppm, 1.0, 1.0, ppm, ppm, ppm,ppm, ppm, ppm, ppm, ppm]\n",
    "unit_list  = ['ppm', ' ', 'deg', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm']\n",
    "\n",
    "for varname in var_i_care:\n",
    "    if varname in fitResults_1st_half_Fourier_PhaseCurve.var_names:\n",
    "        print('{}:\\t{:.0f}\\tppm'.format(varname, fitResults_1st_half_Fourier_PhaseCurve.params[varname]*ppm))\n",
    "\n",
    "print('mean-1:\\t{:.0f}\\tppm'.format((fitResults_1st_half_Fourier_PhaseCurve.params['mean']-1)*ppm))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "plot_model_over_reduced_lc_half_PLDsq_universal(times_0, \n",
    "                                                phots_0 / med_eclipse_1, \n",
    "                                                sqrt(phots_0) / med_eclipse_1, \n",
    "                                                # input_feature_set_0, \n",
    "                                                [phots_0 / med_eclipse_1, gw_0, nbr_0], \n",
    "                                                in_1st_eclipse, \n",
    "                                                fitResults_1st_half_Fourier_PhaseCurve.params, \n",
    "                                                figsize=None, nbins=2000, plotRawData=False,\n",
    "                                                phase_method='fourier', baseline='KRData',\n",
    "                                                ax = None, returnAx = False)\n",
    "# ylim(.975,1.01)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "plot_model_over_reduced_lc_half_PLDsq_universal(times_0, \n",
    "                                                phots_0 / med_eclipse_1, \n",
    "                                                sqrt(phots_0) / med_eclipse_1, \n",
    "                                                # input_feature_set_0, \n",
    "                                                [phots_0 / med_eclipse_1, gw_0, nbr_0], \n",
    "                                                in_1st_eclipse, \n",
    "                                                fitResults_1st_half_Fourier_PhaseCurve.params, \n",
    "                                                figsize=None, nbins=2000, plotRawData=False,\n",
    "                                                phase_method='fourier', baseline='KRData',\n",
    "                                                ax = None, returnAx = False)\n",
    "# ylim(.975,1.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Test PLD Correlations**"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "i_use_pld, corner_vals_pld, corner_names_pld = np.transpose([(k, val, name) \\\n",
    "                                                 for k, (val, name) in enumerate(zip(lmfit_vals, lmfit_names)) \\\n",
    "                                                 if 'pld' in name])\n",
    "i_use_pld = int32(i_use_pld)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "corner_names_pld"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for k in range(len(corner_names_pld)):\n",
    "    corner_names_pld[k] = corner_names_pld[k].replace('_','-')\n",
    "corner_names_pld"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "burnin  = int(n_steps * 0.2)\n",
    "samples = sampler.chain[:, burnin:, i_use_pld].reshape((-1, len(i_use_pld)));\n",
    "# corner.corner(samples, labels=corner_names, truths=corner_vals);\n",
    "\n",
    "#Stairstep Plot\n",
    "# labels = corner_names\n",
    "\n",
    "plt.rc('font',size=8)\n",
    "dims = len(corner_names_pld)\n",
    "# fig,axL = plt.subplots(nrows=dims,ncols=dims,figsize=(15,15))\n",
    "\n",
    "corner_kw = dict(\n",
    "    truths=corner_vals_pld,\n",
    "    labels=corner_names_pld,\n",
    "    levels=[0.68,0.95],\n",
    "    plot_datapoints=False,\n",
    "    smooth=True,\n",
    "    bins=30,\n",
    "    )\n",
    "\n",
    "corner.corner(samples, fig=fig,**corner_kw)\n",
    "# plt.show()\n",
    "# fig.savefig('ch'+str(ch)+'_visit'+str(visit)+'_stairstep_indiv_fit_linreg'+exp_opt_name+'.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2nd AOR Fourier**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_2nd_eclipse_model          = TransitModel(init_eclipse_params, times_1, transittype=\"secondary\")\n",
    "init_2nd_eclipse_lightcurve     = init_2nd_eclipse_model.light_curve(init_eclipse_params)\n",
    "\n",
    "in_2nd_eclipse = init_2nd_eclipse_lightcurve == 1\n",
    "med_eclipse_2  = median(phots_1[in_2nd_eclipse == 1.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_pld_1         = pca.fit_transform(PLDfeatures_1)\n",
    "ica_pld_1         = ica.fit_transform(PLDfeatures_1)\n",
    "\n",
    "pca_scl_pld_1     = pca.fit_transform(scale(PLDfeatures_1))\n",
    "ica_scl_pld_1     = ica.fit_transform(scale(PLDfeatures_1))\n",
    "\n",
    "ica_pca_pld_1     = ica.fit_transform(pca_pld_1)\n",
    "pca_ica_pld_1     = pca.fit_transform(ica_pld_1)\n",
    "\n",
    "ica_pca_scl_pld_1 = ica.fit_transform(pca_scl_pld_1)\n",
    "pca_ica_scl_pld_1 = pca.fit_transform(ica_scl_pld_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# input_feature_set_1   = PLDfeatures_1\n",
    "# input_feature_set_1   = scale(PLDfeatures_1)\n",
    "# input_feature_set_1   = pca_pld_1\n",
    "# input_feature_set_1   = ica_pld_1\n",
    "# input_feature_set_1   = pca_scl_pld_1\n",
    "# input_feature_set_1   = ica_scl_pld_1\n",
    "# input_feature_set_1   = ica_pca_pld_1\n",
    "# input_feature_set_1   = pca_ica_pld_1\n",
    "input_feature_set_1   = ica_pca_scl_pld_1\n",
    "# input_feature_set_1   = pca_ica_scl_pld_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitResults_1st_half_Fourier_PhaseCurve.params.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i_u1, i_u2 = 0.1, 0.1\n",
    "\n",
    "fit_params = fitResults_1st_half_Fourier_PhaseCurve.params\n",
    "\n",
    "initialParams_2nd_half_Fourier_PhaseCurve = Parameters()\n",
    "\n",
    "initialParams_2nd_half_Fourier_PhaseCurve.add_many(\n",
    "    # Planetary Parameters\n",
    "    ('period' , iPeriod              , False),\n",
    "    ('tCenter', fit_params['tCenter'], False , iTCenter-0.01, iTCenter+0.01),\n",
    "    ('inc'    , iInc                 , False, 0.0 ,  90.),\n",
    "    ('aprs'   , iApRs                , False, 0.0 , 100.),\n",
    "    ('edepth' , fit_params['edepth'] , False, 0.0 , 1.0 ),\n",
    "    ('tdepth' , fit_params['tdepth'] , False, 0.0 , 1.0 ),\n",
    "    ('ecc'    , iEcc                 , False, 0.0 , 1.0 ),\n",
    "    ('omega'  , iOmega               , False, 0.0 , 360.),\n",
    "    ('u1'     , fit_params['u1']     , False, 0.0 , 1.0 ),\n",
    "    ('u2'     , fit_params['u2']     , False, 0.0 , 1.0 ),\n",
    "    ('amp1'   , fit_params['amp1']   , True ),\n",
    "    ('amp2'   , fit_params['amp2']   , True , 0.0, inf), # for KBS only\n",
    "    ('amp3'   , fit_params['amp3']   , False),\n",
    "    ('amp4'   , fit_params['amp4']   , False),\n",
    "    ('amp5'   , 0.0                  , False),\n",
    "    ('amp6'   , 0.0                  , False),\n",
    "    ('amp7'   , 0.0                  , False),\n",
    "    ('amp8'   , 0.0                  , False),\n",
    "    ('mean'   , fit_params['mean']   , True ),\n",
    "    \n",
    "    # PLD Coefficients for the 0th AOR - Linear\n",
    "    ('pld1_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld2_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld3_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld4_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld5_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld6_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld7_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld8_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld9_l' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    \n",
    "    # PLD Coefficients for the 0th AOR - Quadratic\n",
    "    ('pld1_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld2_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld3_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld4_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld5_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld6_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld7_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld8_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    ('pld9_q' , 0.0      , True ),# , -10.0, 10.0   ),\n",
    "    \n",
    "    # Out of transit linear baselines\n",
    "    ('intcpt' , fit_params['intcpt'], True ),\n",
    "    ('slope'  , fit_params['slope'] , True ),\n",
    "    ('crvtur' , fit_params['crvtur'], False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit the data with lmfit\n",
    "partial_residuals = partial(residuals_func, data           = phots_1 / med_eclipse_2,\n",
    "                                            times          = times_1, \n",
    "                                            input_features = input_feature_set_1, \n",
    "                                            in_eclipse     = in_2nd_eclipse,\n",
    "                                            method         = 'kbs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mini = Minimizer(partial_residuals, initialParams_2nd_half_Fourier_PhaseCurve)\n",
    "start = time()\n",
    "fitResults_2nd_half_Fourier_PhaseCurve = mini.leastsq()\n",
    "print(\"Full phase curve fitting operation took {} seconds\".format(time()-start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "report_errors(fitResults_2nd_half_Fourier_PhaseCurve.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_i_care = ['edepth', 'ecc', 'omega', 'amp1', 'amp2', 'amp3', 'amp4', 'amp5', 'amp6', 'amp7', 'amp8']\n",
    "ppm_or_not = [ppm, 1.0, 1.0, ppm, ppm, ppm,ppm, ppm, ppm, ppm, ppm]\n",
    "unit_list  = ['ppm', ' ', 'deg', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for varname in var_i_care:\n",
    "    if varname in fitResults_2nd_half_Fourier_PhaseCurve.var_names:\n",
    "        print('{}:\\t{:.0f}\\tppm'.format(varname, fitResults_2nd_half_Fourier_PhaseCurve.params[varname]*ppm))\n",
    "\n",
    "print('mean-1:\\t{:.0f}\\tppm'.format((fitResults_2nd_half_Fourier_PhaseCurve.params['mean']-1)*ppm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_model_over_reduced_lc_half_PLDsq_universal(times_1, phots_1 / med_eclipse_2, \n",
    "                                                sqrt(phots_1) / med_eclipse_2, \n",
    "                                                input_feature_set_1, \n",
    "                                                in_2nd_eclipse, \n",
    "                                                fitResults_2nd_half_Fourier_PhaseCurve.params, \n",
    "                                                figsize=None, nbins=2000, plotRawData=False,\n",
    "                                                method='kbs', ax = None, returnAx = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Fit Full Phase Curve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_params_0  = fitResults_1st_half_Fourier_PhaseCurve.params\n",
    "fit_params_1  = fitResults_2nd_half_Fourier_PhaseCurve.params\n",
    "\n",
    "initialParams_Full_Fourier_PhaseCurve = Parameters()\n",
    "\n",
    "initialParams_Full_Fourier_PhaseCurve.add_many(\n",
    "    # Planetary Parameters\n",
    "    ('period' , iPeriod              , False),\n",
    "    ('tCenter', fit_params['tCenter'], True , iTCenter-0.01, iTCenter+0.01),\n",
    "    ('inc'    , iInc                 , False, 0.0 ,  90.),\n",
    "    ('aprs'   , iApRs                , False, 0.0 , 100.),\n",
    "    ('edepth' , fit_params['edepth'] , True , 0.0 , 1.0 ),\n",
    "    ('tdepth' , fit_params['tdepth'] , True , 0.0 , 1.0 ),\n",
    "    ('ecc'    , iEcc                 , False, 0.0 , 1.0 ),\n",
    "    ('omega'  , iOmega               , False, 0.0 , 360.),\n",
    "    ('u1'     , fit_params['u1']     , True , 0.0 , 1.0 ),\n",
    "    ('u2'     , fit_params['u2']     , True , 0.0 , 1.0 ),\n",
    "    ('amp1'   , fit_params['amp1']   , True ),#,-0.01, 0.01 ),\n",
    "    ('amp2'   , fit_params['amp2']   , True , 0.0, inf),#,-0.01, 0.01 ), # For KBS Only\n",
    "    ('amp3'   , fit_params['amp3']   , False),#,-0.01, 0.01 ),\n",
    "    ('amp4'   , fit_params['amp4']   , False),#,-0.01, 0.01 ),\n",
    "    ('amp5'   , 0.0                  , False),#,-0.01, 0.01 ),\n",
    "    ('amp6'   , 0.0                  , False),#,-0.01, 0.01 ),\n",
    "    ('amp7'   , 0.0                  , False),#,-0.01, 0.01 ),\n",
    "    ('amp8'   , 0.0                  , False),#,-0.01, 0.01 ),\n",
    "    ('mean_0' , fit_params_0['mean'] , True ),#,1-1e-3, 1+1e-3 ),\n",
    "    ('mean_1' , fit_params_1['mean'] , True ),#,1-1e-3, 1+1e-3 ),\n",
    "    \n",
    "    # PLD Coefficients for the 0th AOR - Linear\n",
    "    ('pld1_l_1' , fit_params_0['pld1_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld2_l_1' , fit_params_0['pld2_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld3_l_1' , fit_params_0['pld3_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld4_l_1' , fit_params_0['pld4_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld5_l_1' , fit_params_0['pld5_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld6_l_1' , fit_params_0['pld6_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld7_l_1' , fit_params_0['pld7_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld8_l_1' , fit_params_0['pld8_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld9_l_1' , fit_params_0['pld9_l'], True ),# , -10.0, 10.0   ),\n",
    "    \n",
    "    # PLD Coefficients for the 0th AOR - Quadratic\n",
    "    ('pld1_q_1' , fit_params_0['pld1_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld2_q_1' , fit_params_0['pld2_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld3_q_1' , fit_params_0['pld3_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld4_q_1' , fit_params_0['pld4_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld5_q_1' , fit_params_0['pld5_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld6_q_1' , fit_params_0['pld6_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld7_q_1' , fit_params_0['pld7_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld8_q_1' , fit_params_0['pld8_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld9_q_1' , fit_params_0['pld9_q'], True ),# , -10.0, 10.0   ),\n",
    "    \n",
    "    # PLD Coefficients for the 0th AOR - Linear\n",
    "    ('pld1_l_2' , fit_params_1['pld1_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld2_l_2' , fit_params_1['pld2_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld3_l_2' , fit_params_1['pld3_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld4_l_2' , fit_params_1['pld4_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld5_l_2' , fit_params_1['pld5_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld6_l_2' , fit_params_1['pld6_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld7_l_2' , fit_params_1['pld7_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld8_l_2' , fit_params_1['pld8_l'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld9_l_2' , fit_params_1['pld9_l'], True ),# , -10.0, 10.0   ),\n",
    "    \n",
    "    # PLD Coefficients for the 0th AOR - Quadratic\n",
    "    ('pld1_q_2' , fit_params_1['pld1_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld2_q_2' , fit_params_1['pld2_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld3_q_2' , fit_params_1['pld3_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld4_q_2' , fit_params_1['pld4_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld5_q_2' , fit_params_1['pld5_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld6_q_2' , fit_params_1['pld6_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld7_q_2' , fit_params_1['pld7_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld8_q_2' , fit_params_1['pld8_q'], True ),# , -10.0, 10.0   ),\n",
    "    ('pld9_q_2' , fit_params_1['pld9_q'], True ),# , -10.0, 10.0   ),\n",
    "    \n",
    "    # Out of transit linear baselines\n",
    "    ('intcpt' , fit_params_0['intcpt'], True ),\n",
    "    ('slope'  , fit_params_0['slope'] , True ),\n",
    "    ('crvtur' , fit_params_0['crvtur'], False)\n",
    ")\n",
    "#     # Out of transit linear baselines\n",
    "#     ('intcpt_0' , fit_params_0['intcpt'], True ),\n",
    "#     ('slope_0'  , fit_params_0['slope'] , True ),\n",
    "#     ('crvtur_0' , fit_params_0['crvtur'], True ),\n",
    "    \n",
    "#     ('intcpt_1' , fit_params_1['intcpt'], True ),\n",
    "#     ('slope_1'  , fit_params_1['slope'] , True ),\n",
    "#     ('crvtur_1' , fit_params_1['crvtur'], True ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phots_c = hstack([phots_0 / med_eclipse_1, phots_1 / med_eclipse_2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit the data with lmfit\n",
    "partial_residuals_full = partial(residuals_func_full, data        = phots_c,\n",
    "                                                 times_0          = times_0, \n",
    "                                                 times_1          = times_1,\n",
    "                                                 input_features_0 = input_feature_set_0, \n",
    "                                                 input_features_1 = input_feature_set_1, \n",
    "                                                 in_eclipse_0     = in_1st_eclipse,\n",
    "                                                 in_eclipse_1     = in_2nd_eclipse,\n",
    "                                                 method           = 'kbs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mini = Minimizer(partial_residuals_full, initialParams_Full_Fourier_PhaseCurve)\n",
    "start = time()\n",
    "fitResults_Full_Fourier_PhaseCurve = mini.leastsq()\n",
    "print(\"Full phase curve fitting operation took {} seconds\".format(time()-start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "report_errors(fitResults_Full_Fourier_PhaseCurve.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_i_care = ['edepth', 'ecc', 'omega', 'amp1', 'amp2', 'amp3', 'amp4', 'amp5', 'amp6', 'amp7', 'amp8']\n",
    "ppm_or_not = [ppm, 1.0, 1.0, ppm, ppm, ppm,ppm, ppm, ppm, ppm, ppm]\n",
    "unit_list  = ['ppm', ' ', 'deg', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm', 'ppm']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for varname in var_i_care:\n",
    "    if varname in fitResults_Full_Fourier_PhaseCurve.var_names:\n",
    "        print('{}:\\t{:.0f}\\tppm'.format(varname, fitResults_Full_Fourier_PhaseCurve.params[varname]*ppm))\n",
    "\n",
    "print('mean-1:\\t{:.0f}\\tppm'.format((fitResults_Full_Fourier_PhaseCurve.params['mean_0']-1)*ppm))\n",
    "print('mean-1:\\t{:.0f}\\tppm'.format((fitResults_Full_Fourier_PhaseCurve.params['mean_1']-1)*ppm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot_model_over_reduced_lc_full_PLDsq_universal(times_0, \n",
    "                                                times_1, \n",
    "                                                phots_0 / med_eclipse_1,\n",
    "                                                phots_1 / med_eclipse_2, \n",
    "                                                sqrt(phots_0) / med_eclipse_1, \n",
    "                                                sqrt(phots_1) / med_eclipse_2, \n",
    "                                                input_feature_set_0, \n",
    "                                                input_feature_set_1, \n",
    "                                                in_1st_eclipse, \n",
    "                                                in_2nd_eclipse, \n",
    "                                                fitResults_Full_Fourier_PhaseCurve.params, \n",
    "                                                # initialParams_Full_Fourier_PhaseCurve,\n",
    "                                                refTime = times_0.min(),\n",
    "                                                planetName=planet_params_name,figsize=None, \n",
    "                                                nbins=1000, plotRawData=False, \n",
    "                                                model_color=0, data_color=1, alpha=1.0,\n",
    "                                                method='kbs', ax = None, returnAx = False)\n",
    "\n",
    "# fig_save_name = planet_params_name.replace(' ','_') + '_full_phase_curve_'+channel.replace('/', '')+'_observations.png'\n",
    "\n",
    "# print('Saving tigure to '  + 'figure_results/' + fig_save_name)\n",
    "\n",
    "# fig = gcf()\n",
    "# fig.savefig('figure_results/' + fig_save_name)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}