<commit_before>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.coordinates as coords\n",
    "import astropy.units as u\n",
    "from astropy.time import Time\n",
    "from astroML.time_series import \\\n",
    "    lomb_scargle, lomb_scargle_bootstrap\n",
    "\n",
    "from ztf_summerschool import source_lightcurve, barycenter_times\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Hands-On Exercise 3:\n",
    "Period Finding\n",
    "========\n",
    "One of the fundamental tasks of time-domain astronomy is determining if a source is periodic, and if so, measuring the period.  Period measurements are a vital first step for more detailed scientific study, which may include source classification (e.g., RR Lyrae, W Uma), lightcurve modeling (binaries), or luminosity estimation (Cepheids).\n",
    "\n",
    "Binary stars in particular have lightcurves which may show a wide variety of shapes, depending on the nature of the stars and the system geometry.\n",
    "\n",
    "In this workbook we will develop a basic toolset for the generic problem of finding periodic sources.\n",
    "\n",
    "***\n",
    "by <NAME> (2014-2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the relative-photometry corrected light curves we built in Exercise 2. We'll use the utility function `source_lightcurve` to load the columns MJD, magnitude, and magnitude error.  Note that we will use days as our time coordinate throughout the homework."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# point to our previously-saved data\n",
    "reference_catalog = '../data/PTF_Refims_Files/PTF_d022683_f02_c06_u000114210_p12_sexcat.ctlg'\n",
    "outfile = reference_catalog.split('/')[-1].replace('ctlg','shlv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start by loading the data from our favorite star, which has coordinates $\\alpha_\\mathrm{J2000}, \\delta_\\mathrm{J2000} = (312.503802, -0.706603)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ra_fav, dec_fav = (312.503802, -0.706603)\n",
    "mjds, mags, magerrs = source_lightcurve('../data/'+outfile, ra_fav, dec_fav)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 0:  Barycentering\n",
    "\n",
    "Our times are Modified Julian Date *on earth*.  We need to correct them for Earth's motion around the sun (this is called heliocentering or barycentering).  What is the largest this timing error can be if we do not make this correction?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import astropy.constants as const\n",
    "# CALCULATE AN ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have provided a script to barycenter the data--note that it assumes that the data come from the P48.  Use the `bjds` (barycentered modified julian date) variable through the remainder of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bjds = barycenter_times(mjds,ra_fav,dec_fav)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Optional exercise: plot a histogram of the time differences between the barycentered and non-barycentered data.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: Getting started plotting\n",
    "\n",
    "Complete this function for plotting the lightcurve:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define plot function\n",
    "def plot_data( # COMPLETE THIS LINE\n",
    "    plt.errorbar( # COMPLETE THIS LINE\n",
    "        fmt = '_', capsize=0)\n",
    "    plt.xlabel('Date (MJD)')\n",
    "    plt.ylabel('Magnitude')\n",
    "    plt.gca().invert_yaxis()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run plot function\n",
    "plot_data(bjds, mags, magerrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##The Lomb Scargle Periodogram\n",
    "\n",
    "The [Lomb-Scarge Periodogram](http://adsabs.harvard.edu/abs/1982ApJ...263..835S) provides a method for searching for periodicities in time-series data.  It is comparable to the discrete Fourier Transform, but may be applied to irregularly sampled data.  The periodogram gives as output the relative significance of a least-squares sinusoidal fit to the data as a function of frequency.\n",
    "\n",
    "Much of this presentation follows Ch. 10 of [Ivezic et al.](http://press.princeton.edu/titles/10159.html).\n",
    "\n",
    "We use the \"[generalized](http://adsabs.harvard.edu/abs/2009A%26A...496..577Z)\" LS version implemented in `astroML` rather than the \"standard\" version implemented in `scipy`: the generalized version accounts better for cases of poor sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# documentation for the astroML lomb_scargle function\n",
    "help(lomb_scargle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2: Determining the frequency grid\n",
    "\n",
    "One of the challenges of using the LS periodogram is determining the appropriate frequency grid to search.  We have to select the minimum and maximum frequencies as well as the bin size.\n",
    "\n",
    "If we don't include the true frequency in our search range, we can't find the period!\n",
    "\n",
    "If the bins are too coarse, true peaks may be lost.  If the bins are too fine, the periodogram becomes very slow to compute.\n",
    "\n",
    "The first question to ask is what range of frequencies our data is sensitive to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.1\n",
    "\n",
    "What is the smallest angular frequency $\\omega_{\\rm min}$ our data is sensitive to?  (Hint: smallest frequency => largest time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "freq_min =  # COMPLETE\n",
    "print 'The minimum frequency our data is sensitive to is {:.3f} radian/days, corresponding to a period of {:.3f} days'.format(freq_min, 2*np.pi/freq_min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.2\n",
    "\n",
    "Determining the highest frequency we are sensitive to turns out to be complicated.  \n",
    "\n",
    "if $\\Delta t$ is the difference between consecutive observations,\n",
    "$\\pi$/median($\\Delta t$) is a good starting point, although in practice we may be sensitive to frequencies even higher than $2 \\pi$/min($\\Delta t$) depending on the details of the sampling.\n",
    "\n",
    "What is the largest angular frequency $\\omega_{\\rm max}$ our data is sensitive to?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "freq_max =  # COMPLETE\n",
    "print 'The maximum frequency our data is sensitive to is APPROXIMATELY {:.3f} radian/days, corresponding to a period of {:.3f} days'.format(freq_max, 2*np.pi/freq_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.3\n",
    "\n",
    "We need enough bins to resolve the periodogram peaks, which [have](http://www.astroml.org/gatspy/periodic/lomb_scargle.html) frequency width $\\Delta f \\sim 2\\pi/ (t_{\\rm max} - t_{\\rm min}) = \\omega_{\\rm min}$.\n",
    "If we want to have 5 samples of $\\Delta f$, how many bins will be in our periodogram?  Is this computationally feasible?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_bins =  # COMPLETE\n",
    "print n_bins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.4\n",
    "\n",
    "Let's wrap this work up in a convenience function that takes as input a list of observation times and returns a frequency grid with decent defaults."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define frequency function\n",
    "def frequency_grid(times):\n",
    "    freq_min = # COMPLETE\n",
    "    freq_max = # COMPLETE\n",
    "    n_bins =   # COMPLETE\n",
    "    print 'Using {} bins'.format(n_bins)\n",
    "    return np.linspace(freq_min, freq_max, n_bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run frequency function\n",
    "omegas = frequency_grid(bjds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some cases you'll want to generate the frequency grid by hand, either to extend to higher frequencies (shorter periods) than found by default, to avoid generating too many bins, or to get a more precise estimate of the period.  In that case use the following code.  We'll use a large fixed number of bins to smoothly sample the periodogram as we zoom in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# provided alternate frequency function\n",
    "def alt_frequency_grid(Pmin, Pmax, n_bins = 5000):\n",
    "    \"\"\"Generate an angular frequency grid between Pmin and Pmax (assumed to be in days)\"\"\"\n",
    "    freq_max = 2*np.pi / Pmin\n",
    "    freq_min = 2*np.pi / Pmax\n",
    "    return np.linspace(freq_min, freq_max, n_bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3: Computing the Periodogram\n",
    "\n",
    "Calculate the LS periodiogram and plot the power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# calculate and plot LS periodogram\n",
    "P_LS = lomb_scargle( # COMPLETE\n",
    "plt.plot(omegas, P_LS)\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('$P_{LS}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# provided: define function to find best period\n",
    "def LS_peak_to_period(omegas, P_LS):\n",
    "    \"\"\"find the highest peak in the LS periodogram and return the corresponding period.\"\"\"\n",
    "    max_freq = omegas[np.argmax(P_LS)]\n",
    "    return 2*np.pi/max_freq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run function to find best period\n",
    "best_period = LS_peak_to_period(omegas, P_LS)\n",
    "print \"Best period: {} days\".format(best_period)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4: Phase Calculation\n",
    "Complete this function that returns the phase of an observation (in the range 0-1) given its period.  For simplicity set the zero of the phase to be the time of the initial observation.\n",
    "\n",
    "Hint: Consider the python modulus operator, `%`.\n",
    "\n",
    "Add a keyword that allows your function to have an optional user-settable time of zero phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define function to phase lightcurves\n",
    "def phase(time, period, t0 = None):\n",
    "    \"\"\" Given an input array of times and a period, return the corresponding phase.\"\"\"\n",
    "    if t0 is None:\n",
    "        t0 = time[0]\n",
    "    return # COMPLETE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5: Phase Plotting\n",
    "\n",
    "Plot the phased lightcurve at the best-fit period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define function to plot phased lc\n",
    "def plot_phased_lc(mjds, mags, magerrs, period, t0=None):\n",
    "    phases = # COMPLETE\n",
    "    plt.errorbar( #COMPLETE\n",
    "        fmt = '_', capsize=0)\n",
    "    plt.xlabel('Phase')\n",
    "    plt.ylabel('Magnitude')\n",
    "    plt.gca().invert_yaxis()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run function to plot phased lc\n",
    "plot_phased_lc(bjds, mags, magerrs, best_period)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "How does that look?  Do you think you are close to the right period?\n",
    "\n",
    "Try re-running your analysis using the `alt_frequency_grid` command, searching a narrower period range around the best-fit period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "omegas = alt_frequency_grid( # COMPLETE\n",
    "P_LS = lomb_scargle( # COMPLETE\n",
    "plt.plot(omegas, P_LS)\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('$P_{LS}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "best_period = # COMPLETE\n",
    "print \"Best period: {} days\".format(best_period)\n",
    "plot_phased_lc(bjds, mags, magerrs, best_period)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 6: Calculating significance of the period detection\n",
    "\n",
    "Real data may have *aliases*--frequency components that appear because of the sampling of the data, such as once per night.  Bootstrap significance tests, which shuffle the data values around but keep the times the same, can help rule these out.\n",
    "\n",
    "Calculate the chance probability of finding a LS peak higher than the observed value in random data observed at the specified intervals: use `lomb_scargle_bootstrap` and `np.percentile` to find the 95 and 99 percent significance levels and plot them over the LS power.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "D = lomb_scargle_bootstrap( # COMPLETE\n",
    "sig99, sig95 = np.percentile( # COMPLETE\n",
    "plt.plot(omegas, P_LS)\n",
    "plt.plot([omegas[0],omegas[-1]], sig99*np.ones(2),'--')\n",
    "plt.plot([omegas[0],omegas[-1]], sig95*np.ones(2),'--')\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('$P_{LS}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 7: Find periods of other sources\n",
    "\n",
    "Now find the periods of these sources, plot their phased lightcurves, and evaluate the significance of the period you find:\n",
    "\n",
    "* 312.066287628 -0.983790357518 \n",
    "* 311.967177518 -0.886275170839\n",
    "* 312.263445107 -0.342008023626\n",
    "* 312.050877142 -1.0632849268\n",
    "* 312.293550866 -0.783896411315\n",
    "\n",
    "Suggestion: wrap the code you used above in a function that takes ra & dec as input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Challenge] Exercise 8: gatspy\n",
    "\n",
    "Try using the [`gatspy`](http://www.astroml.org/gatspy/) package to search for periods.  It uses a slightly different interface but has several nice features, such as automatic zooming on candidate frequency peaks.  \n",
    "\n",
    "You'll need to read the online documentation or call `help(gatspy.periodic.LombScargleFast())` to learn how to which commands to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import gatspy\n",
    "ls = gatspy.periodic.LombScargleFast()\n",
    "ls.optimizer.period_range = ( # COMPLETE\n",
    "# we have to subtract the t0 time so the model plotting has the correct phase origin\n",
    "ls.fit(bjds-bjds[0],mags,magerrs)\n",
    "gatspy_period = ls. # COMPLETE\n",
    "print gatspy_period\n",
    "plot_phased_lc(bjds, mags, magerrs, gatspy_period)\n",
    "p = np.linspace(0,gatspy_period,100)\n",
    "plt.plot(p/gatspy_period,ls.predict(p,period=gatspy_period))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Challenge] Exercise 9: Alternate Algorithms \n",
    "\n",
    "Lomb-Scargle is equivalent to fitting a sinusoid to the phased data, but many kinds of variable stars do not have phased lightcurves that are well-represented by a sinusoid.  Other algorithms, such as those that attempt to minimize the dispersion within phase bins over a grid of trial phases, may be more successful in such cases.  See [Graham et al (2013)](http://adsabs.harvard.edu/abs/2013MNRAS.434.3423G) for a review.  \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ss = gatspy.periodic.SuperSmoother(fit_period=True)\n",
    "ss.optimizer.period_range = ( #COMPLETE\n",
    "ss.fit( # COMPLETE\n",
    "gatspy_period = ss. # COMPLETE\n",
    "print gatspy_period\n",
    "plot_phased_lc(bjds, mags, magerrs, gatspy_period)\n",
    "p = np.linspace(0,gatspy_period,100)\n",
    "plt.plot(p/gatspy_period,ss.predict(p,period=gatspy_period))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Challenge] Exercise 10: Multi-harmonic fitting \n",
    "\n",
    "Both AstroML and gatspy include code for including multiple Fourier components in the fit, which can better fit lightcurves that don't have a simple sinusoidal shape (like RR Lyrae)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from astroML.time_series import multiterm_periodogram\n",
    "omegas = alt_frequency_grid(.2,1.2)\n",
    "P_mt = multiterm_periodogram( #COMPLETE\n",
    "plt.plot(omegas, P_mt)\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('$P_{mt}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "best_period = # COMPLETE\n",
    "print \"Best period: {} days\".format(best_period)\n",
    "plot_phased_lc(bjds, mags, magerrs, best_period)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ls = gatspy.periodic.LombScargle(Nterms=4)\n",
    "ls.optimizer.period_range = ( # COMPLETE\n",
    "ls.fit( # COMPLETE\n",
    "gatspy_period = ls. # COMPLETE\n",
    "print gatspy_period\n",
    "plot_phased_lc(bjds, mags, magerrs, gatspy_period)\n",
    "p = np.linspace(0,gatspy_period,100)\n",
    "plt.plot(p/gatspy_period,ls.predict(p,period=gatspy_period))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Challenge] Exercise 11: Compute all periods\n",
    "\n",
    "This is a big one: can you compute periods for *all* of the sources in our database with showing evidence of variability?  How will you compute variablity?  How can you tell which sources are likely to have good periods?\n",
    "\n",
    "We're giving you less guidance here than before--see how you can do!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# open the stored data\n",
    "import shelve\n",
    "import astropy\n",
    "shelf = shelve.open('../data/'+outfile)\n",
    "all_mags = shelf['mags']\n",
    "all_mjds = shelf['mjds']\n",
    "all_errs = shelf['magerrs']\n",
    "all_coords = shelf['ref_coords']\n",
    "shelf.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# loop over stars\n",
    "variable_inds = []\n",
    "best_periods = []\n",
    "best_power = []\n",
    "\n",
    "with astropy.utils.console.ProgressBar(all_mags.shape[0],ipython_widget=True) as bar:\n",
    "    for i in range(all_mags.shape[0]):\n",
    "        # make sure there's real data\n",
    "        wgood = (all_mags[i,:].mask == False)\n",
    "        n_obs = np.sum(wgood)\n",
    "        # if we don't have many observations, don't bother computing periods\n",
    "        if n_obs < 40:\n",
    "            continue # the \"continue\" instruction tells python to skip the rest of the loop for this element and continue with the next one\n",
    "            \n",
    "        # COMPLETE: make a cut so you only complete periods on variabile sources\n",
    "        if # source is not variable: \n",
    "            continue\n",
    "        variable_inds.append(i)\n",
    "        bjds = barycenter_times(all_mjds[wgood],all_coords[i].ra.degree,all_coords[i].dec.degree)\n",
    "        # COMPLETE: calculate best period\n",
    "        best_periods.append( # COMPLETE\n",
    "        best_power.append( # COMPLETE: add the LS power here\n",
    "        bar.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# COMPLETE: now find the most promising periods and plot them!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Other effects to consider*\n",
    "\n",
    "\n",
    "\n",
    "Many eclipsing binaries have primary and secondary eclipses, often with comparable depths.  The period found by LS (which fits a single sinusoid) will thus often be only half of the true period.  Plotting the phased lightcurve at double the LS period is often the easiest way to determine the true period. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References and Further Reading\n",
    "\n",
    "<NAME>. 1982, [ApJ 263, 835](http://adsabs.harvard.edu/abs/1982ApJ...263..835S)\n",
    "\n",
    "<NAME>., and <NAME>. 2009, [A&A 496, 577](http://adsabs.harvard.edu/abs/2009A%26A...496..577Z)\n",
    "\n",
    "<NAME>. et al. 2013, [MNRAS 434, 3423](http://adsabs.harvard.edu/abs/2013MNRAS.434.3423G)\n",
    "\n",
    "[Statistics, Data Mining, and Machine Learning in Astronomy](http://press.princeton.edu/titles/10159.html) (Ivezic, Connolly, VanderPlas, & Gray)\n",
    "\n",
    "[gatspy documentation](http://www.astroml.org/gatspy/index.html)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}<commit_msg>remove barycenter calculation, convert print statements to functions<commit_after>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# imports\n",
    "from __future__ import division, print_function, absolute_import, unicode_literals\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.coordinates as coords\n",
    "import astropy.units as u\n",
    "from astropy.time import Time\n",
    "from astroML.time_series import \\\n",
    "    lomb_scargle, lomb_scargle_bootstrap\n",
    "\n",
    "from ztf_summerschool import source_lightcurve, barycenter_times\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hands-On Exercise 3: Period Finding\n",
    "---\n",
    "\n",
    "One of the fundamental tasks of time-domain astronomy is determining if a source is periodic, and if so, measuring the period.  Period measurements are a vital first step for more detailed scientific study, which may include source classification (e.g., RR Lyrae, W Uma), lightcurve modeling (binaries), or luminosity estimation (Cepheids).\n",
    "\n",
    "Binary stars in particular have lightcurves which may show a wide variety of shapes, depending on the nature of the stars and the system geometry.\n",
    "\n",
    "In this workbook we will develop a basic toolset for the generic problem of finding periodic sources.\n",
    "\n",
    "***\n",
    "by <NAME> (2014-2016)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the relative-photometry corrected light curves we built in Exercise 2. We'll use the utility function `source_lightcurve` to load the columns MJD, magnitude, and magnitude error.  Note that we will use days as our time coordinate throughout the homework."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# point to our previously-saved data\n",
    "reference_catalog = '../data/PTF_Refims_Files/PTF_d022683_f02_c06_u000114210_p12_sexcat.ctlg'\n",
    "outfile = reference_catalog.split('/')[-1].replace('ctlg','shlv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start by loading the data from our favorite star, which has coordinates $\\alpha_\\mathrm{J2000}, \\delta_\\mathrm{J2000} = (312.503802, -0.706603)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ra_fav, dec_fav = (312.503802, -0.706603)\n",
    "mjds, mags, magerrs = source_lightcurve('../data/'+outfile, ra_fav, dec_fav)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Barycentering\n",
    "\n",
    "Our times are Modified Julian Date *on earth*.  We need to correct them for Earth's motion around the sun (this is called heliocentering or barycentering).  What is the largest this timing error can be if we do not make this correction?  It's the light travel time over one AU.  We can use `astropy` constants to calculate this easily:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$8.3167464 \\; \\mathrm{min}$"
      ],
      "text/plain": [
       "<Quantity 8.316746397269274 min>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import astropy.constants as const\n",
    "(const.au / const.c).to(u.minute)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have provided a script to barycenter the data--note that it assumes that the data come from the P48.  Use the `bjds` (barycentered modified julian date) variable through the remainder of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading http://maia.usno.navy.mil/ser7/finals2000A.all [Done]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ebellm/Eric/PTF/Meetings/iPTF Summer School 2016/ztf_summerschool_2015/notebooks/BaryTime.py:61: UserWarning: Cannot calculate UT1: using UTC with degraded accuracy\n",
      "  warnings.warn('Cannot calculate UT1: using UTC with degraded accuracy')\n"
     ]
    }
   ],
   "source": [
    "bjds = barycenter_times(mjds,ra_fav,dec_fav)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Optional exercise: plot a histogram of the time differences between the barycentered and non-barycentered data.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: Getting started plotting\n",
    "\n",
    "Complete this function for plotting the lightcurve:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define plot function\n",
    "def plot_data( # COMPLETE THIS LINE\n",
    "    plt.errorbar( # COMPLETE THIS LINE\n",
    "        fmt = '_', capsize=0)\n",
    "    plt.xlabel('Date (MJD)')\n",
    "    plt.ylabel('Magnitude')\n",
    "    plt.gca().invert_yaxis()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run plot function\n",
    "plot_data(bjds, mags, magerrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Lomb Scargle Periodogram\n",
    "\n",
    "The [Lomb-Scarge Periodogram](http://adsabs.harvard.edu/abs/1982ApJ...263..835S) provides a method for searching for periodicities in time-series data.  It is comparable to the discrete Fourier Transform, but may be applied to irregularly sampled data.  The periodogram gives as output the relative significance of a least-squares sinusoidal fit to the data as a function of frequency.\n",
    "\n",
    "Much of this presentation follows Ch. 10 of [Ivezic et al.](http://press.princeton.edu/titles/10159.html).\n",
    "\n",
    "We use the \"[generalized](http://adsabs.harvard.edu/abs/2009A%26A...496..577Z)\" LS version implemented in `astroML` rather than the \"standard\" version implemented in `scipy`: the generalized version accounts better for cases of poor sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# documentation for the astroML lomb_scargle function\n",
    "help(lomb_scargle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2: Determining the frequency grid\n",
    "\n",
    "One of the challenges of using the LS periodogram is determining the appropriate frequency grid to search.  We have to select the minimum and maximum frequencies as well as the bin size.\n",
    "\n",
    "If we don't include the true frequency in our search range, we can't find the period!\n",
    "\n",
    "If the bins are too coarse, true peaks may be lost.  If the bins are too fine, the periodogram becomes very slow to compute.\n",
    "\n",
    "The first question to ask is what range of frequencies our data is sensitive to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.1\n",
    "\n",
    "What is the smallest angular frequency $\\omega_{\\rm min}$ our data is sensitive to?  (Hint: smallest frequency => largest time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "freq_min =  # COMPLETE\n",
    "print('The minimum frequency our data is sensitive to is {:.3f} radian/days, corresponding to a period of {:.3f} days'.format(freq_min, 2*np.pi/freq_min)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.2\n",
    "\n",
    "Determining the highest frequency we are sensitive to turns out to be complicated.  \n",
    "\n",
    "if $\\Delta t$ is the difference between consecutive observations,\n",
    "$\\pi$/median($\\Delta t$) is a good starting point, although in practice we may be sensitive to frequencies even higher than $2 \\pi$/min($\\Delta t$) depending on the details of the sampling.\n",
    "\n",
    "What is the largest angular frequency $\\omega_{\\rm max}$ our data is sensitive to?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "freq_max =  # COMPLETE\n",
    "print('The maximum frequency our data is sensitive to is APPROXIMATELY {:.3f} radian/days, corresponding to a period of {:.3f} days'.format(freq_max, 2*np.pi/freq_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.3\n",
    "\n",
    "We need enough bins to resolve the periodogram peaks, which [have](http://www.astroml.org/gatspy/periodic/lomb_scargle.html) frequency width $\\Delta f \\sim 2\\pi/ (t_{\\rm max} - t_{\\rm min}) = \\omega_{\\rm min}$.\n",
    "If we want to have 5 samples of $\\Delta f$, how many bins will be in our periodogram?  Is this computationally feasible?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_bins =  # COMPLETE\n",
    "print(n_bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.4\n",
    "\n",
    "Let's wrap this work up in a convenience function that takes as input a list of observation times and returns a frequency grid with decent defaults."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define frequency function\n",
    "def frequency_grid(times):\n",
    "    freq_min = # COMPLETE\n",
    "    freq_max = # COMPLETE\n",
    "    n_bins =   # COMPLETE\n",
    "    print('Using {} bins'.format(n_bins))\n",
    "    return np.linspace(freq_min, freq_max, n_bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run frequency function\n",
    "omegas = frequency_grid(bjds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some cases you'll want to generate the frequency grid by hand, either to extend to higher frequencies (shorter periods) than found by default, to avoid generating too many bins, or to get a more precise estimate of the period.  In that case use the following code.  We'll use a large fixed number of bins to smoothly sample the periodogram as we zoom in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# provided alternate frequency function\n",
    "def alt_frequency_grid(Pmin, Pmax, n_bins = 5000):\n",
    "    \"\"\"Generate an angular frequency grid between Pmin and Pmax (assumed to be in days)\"\"\"\n",
    "    freq_max = 2*np.pi / Pmin\n",
    "    freq_min = 2*np.pi / Pmax\n",
    "    return np.linspace(freq_min, freq_max, n_bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3: Computing the Periodogram\n",
    "\n",
    "Calculate the LS periodiogram and plot the power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# calculate and plot LS periodogram\n",
    "P_LS = lomb_scargle( # COMPLETE\n",
    "plt.plot(omegas, P_LS)\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('$P_{LS}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# provided: define function to find best period\n",
    "def LS_peak_to_period(omegas, P_LS):\n",
    "    \"\"\"find the highest peak in the LS periodogram and return the corresponding period.\"\"\"\n",
    "    max_freq = omegas[np.argmax(P_LS)]\n",
    "    return 2*np.pi/max_freq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run function to find best period\n",
    "best_period = LS_peak_to_period(omegas, P_LS)\n",
    "print(\"Best period: {} days\".format(best_period))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4: Phase Calculation\n",
    "Complete this function that returns the phase of an observation (in the range 0-1) given its period.  For simplicity set the zero of the phase to be the time of the initial observation.\n",
    "\n",
    "Hint: Consider the python modulus operator, `%`.\n",
    "\n",
    "Add a keyword that allows your function to have an optional user-settable time of zero phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define function to phase lightcurves\n",
    "def phase(time, period, t0 = None):\n",
    "    \"\"\" Given an input array of times and a period, return the corresponding phase.\"\"\"\n",
    "    if t0 is None:\n",
    "        t0 = time[0]\n",
    "    return # COMPLETE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5: Phase Plotting\n",
    "\n",
    "Plot the phased lightcurve at the best-fit period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define function to plot phased lc\n",
    "def plot_phased_lc(mjds, mags, magerrs, period, t0=None):\n",
    "    phases = # COMPLETE\n",
    "    plt.errorbar( #COMPLETE\n",
    "        fmt = '_', capsize=0)\n",
    "    plt.xlabel('Phase')\n",
    "    plt.ylabel('Magnitude')\n",
    "    plt.gca().invert_yaxis()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run function to plot phased lc\n",
    "plot_phased_lc(bjds, mags, magerrs, best_period)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "How does that look?  Do you think you are close to the right period?\n",
    "\n",
    "Try re-running your analysis using the `alt_frequency_grid` command, searching a narrower period range around the best-fit period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "omegas = alt_frequency_grid( # COMPLETE\n",
    "P_LS = lomb_scargle( # COMPLETE\n",
    "plt.plot(omegas, P_LS)\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('$P_{LS}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "best_period = # COMPLETE\n",
    "print(\"Best period: {} days\".format(best_period))\n",
    "plot_phased_lc(bjds, mags, magerrs, best_period)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 6: Calculating significance of the period detection\n",
    "\n",
    "Real data may have *aliases*--frequency components that appear because of the sampling of the data, such as once per night.  Bootstrap significance tests, which shuffle the data values around but keep the times the same, can help rule these out.\n",
    "\n",
    "Calculate the chance probability of finding a LS peak higher than the observed value in random data observed at the specified intervals: use `lomb_scargle_bootstrap` and `np.percentile` to find the 95 and 99 percent significance levels and plot them over the LS power.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "D = lomb_scargle_bootstrap( # COMPLETE\n",
    "sig99, sig95 = np.percentile( # COMPLETE\n",
    "plt.plot(omegas, P_LS)\n",
    "plt.plot([omegas[0],omegas[-1]], sig99*np.ones(2),'--')\n",
    "plt.plot([omegas[0],omegas[-1]], sig95*np.ones(2),'--')\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('$P_{LS}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 7: Find periods of other sources\n",
    "\n",
    "Now find the periods of these sources, plot their phased lightcurves, and evaluate the significance of the period you find:\n",
    "\n",
    "* 312.066287628 -0.983790357518 \n",
    "* 311.967177518 -0.886275170839\n",
    "* 312.263445107 -0.342008023626\n",
    "* 312.050877142 -1.0632849268\n",
    "* 312.293550866 -0.783896411315\n",
    "\n",
    "Suggestion: wrap the code you used above in a function that takes ra & dec as input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Challenge] Exercise 8: gatspy\n",
    "\n",
    "Try using the [`gatspy`](http://www.astroml.org/gatspy/) package to search for periods.  It uses a slightly different interface but has several nice features, such as automatic zooming on candidate frequency peaks.  \n",
    "\n",
    "You'll need to read the online documentation or call `help(gatspy.periodic.LombScargleFast())` to learn how to which commands to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import gatspy\n",
    "ls = gatspy.periodic.LombScargleFast()\n",
    "ls.optimizer.period_range = ( # COMPLETE\n",
    "# we have to subtract the t0 time so the model plotting has the correct phase origin\n",
    "ls.fit(bjds-bjds[0],mags,magerrs)\n",
    "gatspy_period = ls. # COMPLETE\n",
    "print(gatspy_period)\n",
    "plot_phased_lc(bjds, mags, magerrs, gatspy_period)\n",
    "p = np.linspace(0,gatspy_period,100)\n",
    "plt.plot(p/gatspy_period,ls.predict(p,period=gatspy_period))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Challenge] Exercise 9: Alternate Algorithms \n",
    "\n",
    "Lomb-Scargle is equivalent to fitting a sinusoid to the phased data, but many kinds of variable stars do not have phased lightcurves that are well-represented by a sinusoid.  Other algorithms, such as those that attempt to minimize the dispersion within phase bins over a grid of trial phases, may be more successful in such cases.  See [Graham et al (2013)](http://adsabs.harvard.edu/abs/2013MNRAS.434.3423G) for a review.  \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ss = gatspy.periodic.SuperSmoother(fit_period=True)\n",
    "ss.optimizer.period_range = ( #COMPLETE\n",
    "ss.fit( # COMPLETE\n",
    "gatspy_period = ss. # COMPLETE\n",
    "print(gatspy_period)\n",
    "plot_phased_lc(bjds, mags, magerrs, gatspy_period)\n",
    "p = np.linspace(0,gatspy_period,100)\n",
    "plt.plot(p/gatspy_period,ss.predict(p,period=gatspy_period))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Challenge] Exercise 10: Multi-harmonic fitting \n",
    "\n",
    "Both AstroML and gatspy include code for including multiple Fourier components in the fit, which can better fit lightcurves that don't have a simple sinusoidal shape (like RR Lyrae)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from astroML.time_series import multiterm_periodogram\n",
    "omegas = alt_frequency_grid(.2,1.2)\n",
    "P_mt = multiterm_periodogram( #COMPLETE\n",
    "plt.plot(omegas, P_mt)\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('$P_{mt}$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "best_period = # COMPLETE\n",
    "print(\"Best period: {} days\".format(best_period))\n",
    "plot_phased_lc(bjds, mags, magerrs, best_period)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ls = gatspy.periodic.LombScargle(Nterms=4)\n",
    "ls.optimizer.period_range = ( # COMPLETE\n",
    "ls.fit( # COMPLETE\n",
    "gatspy_period = ls. # COMPLETE\n",
    "print(gatspy_period)\n",
    "plot_phased_lc(bjds, mags, magerrs, gatspy_period)\n",
    "p = np.linspace(0,gatspy_period,100)\n",
    "plt.plot(p/gatspy_period,ls.predict(p,period=gatspy_period))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Challenge] Exercise 11: Compute all periods\n",
    "\n",
    "This is a big one: can you compute periods for *all* of the sources in our database with showing evidence of variability?  How will you compute variablity?  How can you tell which sources are likely to have good periods?\n",
    "\n",
    "We're giving you less guidance here than before--see how you can do!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# open the stored data\n",
    "import shelve\n",
    "import astropy\n",
    "shelf = shelve.open('../data/'+outfile)\n",
    "all_mags = shelf['mags']\n",
    "all_mjds = shelf['mjds']\n",
    "all_errs = shelf['magerrs']\n",
    "all_coords = shelf['ref_coords']\n",
    "shelf.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# loop over stars\n",
    "variable_inds = []\n",
    "best_periods = []\n",
    "best_power = []\n",
    "\n",
    "with astropy.utils.console.ProgressBar(all_mags.shape[0],ipython_widget=True) as bar:\n",
    "    for i in range(all_mags.shape[0]):\n",
    "        # make sure there's real data\n",
    "        wgood = (all_mags[i,:].mask == False)\n",
    "        n_obs = np.sum(wgood)\n",
    "        # if we don't have many observations, don't bother computing periods\n",
    "        if n_obs < 40:\n",
    "            continue # the \"continue\" instruction tells python to skip the rest of the loop for this element and continue with the next one\n",
    "            \n",
    "        # COMPLETE: make a cut so you only complete periods on variabile sources\n",
    "        if # source is not variable: \n",
    "            continue\n",
    "        variable_inds.append(i)\n",
    "        bjds = barycenter_times(all_mjds[wgood],all_coords[i].ra.degree,all_coords[i].dec.degree)\n",
    "        # COMPLETE: calculate best period\n",
    "        best_periods.append( # COMPLETE\n",
    "        best_power.append( # COMPLETE: add the LS power here\n",
    "        bar.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# COMPLETE: now find the most promising periods and plot them!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Other effects to consider*\n",
    "\n",
    "\n",
    "\n",
    "Many eclipsing binaries have primary and secondary eclipses, often with comparable depths.  The period found by LS (which fits a single sinusoid) will thus often be only half of the true period.  Plotting the phased lightcurve at double the LS period is often the easiest way to determine the true period. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References and Further Reading\n",
    "\n",
    "<NAME>. 1982, [ApJ 263, 835](http://adsabs.harvard.edu/abs/1982ApJ...263..835S)\n",
    "\n",
    "<NAME>., and <NAME>. 2009, [A&A 496, 577](http://adsabs.harvard.edu/abs/2009A%26A...496..577Z)\n",
    "\n",
    "<NAME>. et al. 2013, [MNRAS 434, 3423](http://adsabs.harvard.edu/abs/2013MNRAS.434.3423G)\n",
    "\n",
    "[Statistics, Data Mining, and Machine Learning in Astronomy](http://press.princeton.edu/titles/10159.html) (Ivezic, Connolly, VanderPlas, & Gray)\n",
    "\n",
    "[gatspy documentation](http://www.astroml.org/gatspy/index.html)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [summerschool]",
   "language": "python",
   "name": "Python [summerschool]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
{'fraction_non_alphanumeric': 0.1666867107636801, 'fraction_numerical': 0.01723792343154941, 'mean_word_length': 2.864224304856324, 'pattern_counts': {'":': 592, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5977315', 'n_tokens_mistral': 18420, 'n_tokens_neox': 17394, 'n_words': 4637}