K2 difference image demo

Difference_image_demo_with_strategic_scene_modeling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# K2 difference imaging for transients\n",
    "Jan 16, 2019  \n",
    "Revised March 2019  \n",
    "At NASA Ames  \n",
    "M. Gully-Santiago and RRH."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load /Users/obsidian/Desktop/defaults.py\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "import lightkurve as lk\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ids = pd.read_csv('GO6077-targets.csv')['EPIC ID'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def heuristic_scene_model(frameno, dist_matrix):\n",
    "    this_dist = dist_matrix[frameno, :]\n",
    "    sort_indices = np.argsort(this_dist)\n",
    "    \n",
    "    ### Tuning parameters\n",
    "    # Keep the top ~100 closest cadences\n",
    "    top_N_dist = 100\n",
    "    # Of those, keep the 11 most-distant-in-time cadences\n",
    "    top_N_time = 11\n",
    "    # Make sure that the kept cadences are at least within a minimal proximity\n",
    "    minimal_proximity = 0.1 #pixels\n",
    "    # Make sure that the difference image is composed of frames that are at at least \n",
    "    #    10 days earlier or later\n",
    "    minimal_time_difference = 10.0 #days\n",
    "\n",
    "    # Sort by both distance and then time\n",
    "    gi = sort_indices[1:top_N_dist]\n",
    "    nearest_distances = this_dist[gi]\n",
    "    \n",
    "    g_time = tpf.time[gi]\n",
    "    time_sort_indices = np.argsort(np.abs(g_time - tpf.time[frameno]))[::-1]\n",
    "    \n",
    "    subset = gi[time_sort_indices[0:top_N_time]]\n",
    "    #Safety net: our scene model is not satisfiable!\n",
    "    if any(( this_dist[subset] > minimal_proximity) | \n",
    "            (np.abs(tpf.time[subset] - tpf.time[frameno]) < minimal_time_difference)):\n",
    "        median_frame = tpf.flux[frameno] *np.NaN\n",
    "        std_frame = tpf.flux[frameno] *np.NaN\n",
    "        return median_frame, std_frame\n",
    "                    \n",
    "    median_frame = np.nanmedian(tpf.flux[subset], axis=(0))\n",
    "    std_frame = np.nanstd(tpf.flux[subset], axis=(0))\n",
    "    return median_frame, std_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def diff_im_event_detector_wide_time(tpf):\n",
    "    dist_matrix = np.sqrt( (tpf.pos_corr1 - tpf.pos_corr1[:, np.newaxis])**2 + \n",
    "                           (tpf.pos_corr2 - tpf.pos_corr2[:, np.newaxis])**2 )\n",
    "    n_cad, _, _ = tpf.shape\n",
    "    diffs = np.zeros(tpf.shape)\n",
    "    sigmas = np.zeros(tpf.shape)\n",
    "    for frameno in tqdm(np.arange(0, n_cad, 1, dtype = np.int)):\n",
    "        median_frame, sigmas[frameno, :, :]  = heuristic_scene_model(frameno, dist_matrix)\n",
    "        diffs[frameno, :, :] = tpf.flux[frameno, :, :] - median_frame\n",
    "        \n",
    "    sigma_datacube = diffs/sigmas\n",
    "    detection_threshold = 3\n",
    "    detection_datacube = sigma_datacube > detection_threshold\n",
    "    return detection_datacube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for i, target in enumerate(ids[0:4]):\n",
    "    print(target)\n",
    "    try:\n",
    "        result = lk.search_targetpixelfile(target)\n",
    "        tpf = result[0].download()\n",
    "        good_mask = tpf.pos_corr1 > -10\n",
    "        tpf = tpf[good_mask]\n",
    "        detection_cube = diff_im_event_detector_wide_time(tpf)\n",
    "        detection_counts = np.nansum(np.abs(detection_cube), axis=(0))\n",
    "        fail_counts = detection_cube!=detection_cube\n",
    "        print(np.sum(fail_counts[:, 0, 0]))\n",
    "\n",
    "        fig, ax1 = plt.subplots(figsize=(4, 4), ncols=1)\n",
    "        l = ax1.imshow(detection_counts, vmin=0, interpolation=None)\n",
    "        fig.colorbar(l)\n",
    "        plt.show()\n",
    "    except:\n",
    "        print(\"{} didn't work!\".format(target))\n",
    "        raise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}