db_permute_kic965.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from toolkit import slice_spectrum, concatenate_spectra, bands_TiO, match_spectra, SimpleSpectrum\n",
    "import astropy.units as u\n",
    "from scipy.optimize import fmin_l_bfgs_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_spliced_spectrum(observed_spectrum, model_flux, other_model=None):\n",
    "    n_chunks = len(observed_spectrum.wavelength_splits)\n",
    "    fig, ax = plt.subplots(n_chunks, 1, figsize=(8, 10))\n",
    "\n",
    "    for i, inds in enumerate(observed_spectrum.wavelength_splits):\n",
    "        min_ind, max_ind = inds\n",
    "        \n",
    "        ax[i].errorbar(observed_spectrum.wavelength[min_ind:max_ind].value, \n",
    "                       observed_spectrum.flux[min_ind:max_ind], \n",
    "                       0.025*np.ones(max_ind-min_ind))\n",
    "        ax[i].plot(observed_spectrum.wavelength[min_ind:max_ind], \n",
    "                   model_flux[min_ind:max_ind])\n",
    "        \n",
    "        if other_model is not None:\n",
    "            ax[i].plot(observed_spectrum.wavelength[min_ind:max_ind], \n",
    "                       other_model[min_ind:max_ind], alpha=0.4)\n",
    "        \n",
    "        ax[i].set_xlim([observed_spectrum.wavelength[min_ind].value,\n",
    "                        observed_spectrum.wavelength[max_ind-1].value])\n",
    "        ax[i].set_ylim([0.9*observed_spectrum.flux[min_ind:max_ind].min(), \n",
    "                        1.1])\n",
    "\n",
    "    return fig, ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "\n",
    "archive = h5py.File('/Users/bmmorris/git/aesop/notebooks/spectra.hdf5', 'r+')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['2017-09-05T06:04:52.570',\n",
       " '2017-09-05T06:09:54.789',\n",
       " '2017-09-11T01:56:48.809',\n",
       " '2017-09-11T02:02:20.940']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(archive['51Peg'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets = list(archive)\n",
    "spectrum1 = archive['KIC9652680']['2016-07-05T09:58:18.760']\n",
    "wavelength1 = spectrum1['wavelength'][:]\n",
    "flux1 = spectrum1['flux'][:]\n",
    "target = SimpleSpectrum(wavelength1, flux1, dispersion_unit=u.Angstrom)\n",
    "\n",
    "spec_band = []\n",
    "\n",
    "first_n_bands = 5\n",
    "width = 3\n",
    "\n",
    "for band in bands_TiO[:first_n_bands]:\n",
    "    target_slice = slice_spectrum(target, band.min-width*u.Angstrom, band.max+width*u.Angstrom)\n",
    "    target_slice.flux /= target_slice.flux.max()\n",
    "    spec_band.append(target_slice)\n",
    "\n",
    "target_slices = concatenate_spectra(spec_band)\n",
    "\n",
    "spectrum2 = archive['51Peg']['2017-09-11T01:56:48.809']#archive[target][time]\n",
    "wavelength2 = spectrum2['wavelength'][:]\n",
    "flux2 = spectrum2['flux'][:]\n",
    "source1 = SimpleSpectrum(wavelength2, flux2, dispersion_unit=u.Angstrom)\n",
    "\n",
    "spec_band = []\n",
    "for band, inds in zip(bands_TiO[:first_n_bands], target_slices.wavelength_splits):\n",
    "    target_slice = slice_spectrum(source1, band.min-width*u.Angstrom, band.max+width*u.Angstrom, \n",
    "                                  force_length=abs(np.diff(inds))[0])\n",
    "    target_slice.flux /= target_slice.flux.max()\n",
    "    spec_band.append(target_slice)\n",
    "\n",
    "source1_slices = concatenate_spectra(spec_band)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi2(p, target, comp1, comp2):\n",
    "    model, residuals = instr_model(target, comp1, comp2, *p)\n",
    "    return residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/bmmorris/anaconda/lib/python3.5/site-packages/scipy/signal/windows.py:1217: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  w = np.exp(-n ** 2 / sig2)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GJ4099 0.273652152867 0.026156648259\n",
      "GJ4099 0.376290139184 1.02404502467e-13\n",
      "GJ702B 0.275963237122 0.0855269638279\n",
      "GJ9781A 0.272833269119 0.130178276445\n",
      "HATP11 0.278672932811 0.00019095726489\n",
      "HATP11 0.278661499135 0.000208422637521\n",
      "HATP11 0.278651446581 1.38773112835e-06\n",
      "HATP11 0.278651400327 1.32475635597e-06\n",
      "HATP11 0.278651523315 9.56440103778e-06\n",
      "HATP11 0.278665132305 7.42036937231e-07\n",
      "HATP11 0.334050056326 0.088977723921\n",
      "HATP11 0.27865139539 6.23519836344e-07\n",
      "HATP11 0.278665073365 1.20967661284e-06\n",
      "HATP11 0.336609394055 0.000434148550282\n",
      "HATP11 0.278651504294 3.33003944561e-06\n",
      "HATP11 0.278712917853 0.00129382965263\n",
      "HATP11 0.278651394502 4.23187343605e-07\n",
      "HATP11 0.278665174804 1.41563484884e-05\n",
      "HATP11 0.382699259367 0.0407034047833\n",
      "HATP11 0.278677645683 0.000157068489392\n",
      "HATP11 0.278665730009 6.53814327875e-06\n",
      "HATP11 0.278665113991 1.42001232849e-07\n",
      "HATP11 0.278651224195 9.35762296884e-14\n",
      "HATP11 0.278651408928 1.09953441812e-09\n",
      "HD10697 0.273263535426 0.304643786395\n",
      "HD110833 0.276365962593 0.0937953421964\n",
      "HD110833 0.276623653082 0.0875812158434\n",
      "HD113827 0.317032709745 0.251295121635\n",
      "HD122120 0.27720579306 0.0641572351455\n",
      "HD127506 0.270261518276 0.181252276823\n",
      "HD129333 0.251884190596 0.388692716968\n",
      "HD134319 0.294059109745 0.4668635179\n",
      "HD14039 0.336626898416 9.03635278077e-08\n",
      "HD14039 0.278664994894 2.78269146368e-06\n",
      "HD148467 0.331238781729 0.128322947583\n",
      "HD149957 0.273826347326 0.0841574380502\n",
      "HD151288 0.271315112875 0.0924917219057\n",
      "HD175742 0.270745888362 0.270770537727\n",
      "HD178126 0.278576406904 0.00754195365724\n",
      "HD182488 0.277493334177 0.0841371411927\n",
      "HD182488 0.276492910526 0.11038800228\n",
      "HD182488 0.27813102339 0.085505560132\n",
      "HD182488 0.278196327919 0.0616718717743\n",
      "HD200560 0.380731342144 0.134877810227\n",
      "HD209290 0.27062158484 0.0546387854629\n",
      "HD210277 0.276877545441 0.101658633243\n",
      "HD210277 0.318628157386 0.344742966278\n",
      "HD220182 0.266616843857 0.265328423008\n",
      "HD221639 0.278681101799 0.000121196687314\n",
      "HD221639 0.278665115057 5.22871187424e-06\n",
      "HD222107 0.276656135289 0.122964053781\n",
      "HD222107 0.276187613575 0.134387209807\n",
      "HD266611 0.269979787074 0.0924448106367\n",
      "HD266611 0.269745947307 0.0934112227727\n",
      "HD266611 0.270013407562 0.0924577807907\n",
      "HD34411 0.321184609674 0.535135699557\n",
      "HD34411 0.25555108269 0.310685880564\n",
      "HD38230 0.27840629456 0.0491781547124\n",
      "HD38230 0.277507841408 0.102763611252\n",
      "HD39587 0.246809249948 0.420588836135\n",
      "HD39587 0.307364409532 0.707302613143\n",
      "HD41593 0.268518567394 0.263966997881\n",
      "HD41593 0.268822476463 0.24554081072\n",
      "HD42250 0.275204723384 0.153003946934\n",
      "HD42250 0.27392269516 0.170550438356\n",
      "HD42250 0.271332412531 0.197890790841\n",
      "HD42250 0.301521531662 0.151471631554\n",
      "HD45088 0.30938700688 0.383629278256\n",
      "HD45088 0.265386197782 0.256573638271\n",
      "HD47752 0.272036099022 0.170582588842\n",
      "HD47752 0.272238763588 0.16754930878\n",
      "HD47752 0.273155468999 0.14893888522\n",
      "HD47752 0.273455040101 0.148692845664\n",
      "HD50692 0.256494234302 0.419001797342\n",
      "HD5857 0.278665117308 1.04360965725e-06\n",
      "HD62613 0.260419775541 0.451113252617\n",
      "HD6497 0.278651395614 8.17915426107e-08\n",
      "HD6497 0.278651405755 8.16160075569e-08\n",
      "HD67767 0.278221729478 0.0568233622928\n",
      "HD67767 0.27822518037 0.0629518845416\n",
      "HD68017 0.262587531863 0.27356600043\n",
      "HD68017 0.324070900165 0.592184400306\n",
      "HD73667 0.304872672664 0.516374749818\n",
      "HD79555 0.264811976519 0.214952827034\n",
      "HD79555 0.26484975139 0.218801978891\n",
      "HD79555 0.3136299944 0.296329502925\n",
      "HD82106 0.26871501862 0.182705475421\n",
      "HD82106 0.269453289136 0.17537762978\n",
      "HD82106 0.269384151076 0.17446670771\n",
      "HD86728 0.263340326258 0.229371126\n",
      "HD86728 0.265254064534 0.244675421256\n",
      "HD87884 0.266772947054 0.234099750461\n",
      "HD87884 0.308250236546 0.368672227782\n",
      "HD87884 0.265643048281 0.243575439927\n",
      "HD88230 0.272650793589 0.077341944803\n",
      "HD89744 0.250203466225 0.516621318119\n",
      "HD98230 0.2534825712 0.34164134165\n",
      "HD98230 0.250254675589 0.358736673515\n",
      "HD98230 0.252541650838 0.354538822904\n",
      "HR8832 0.275623731308 0.0351148908694\n",
      "HR8832 0.278651399319 2.13361280512e-11\n",
      "sigmaDraconis 0.265929473982 0.265189117268\n",
      "sigmaDraconis 0.311050458324 0.374534381465\n"
     ]
    }
   ],
   "source": [
    "from toolkit import instr_model\n",
    "from astropy.utils.console import ProgressBar\n",
    "bfgs_options_fast = dict(epsilon=1e-3, approx_grad=True,\n",
    "                         m=10, maxls=20)\n",
    "bfgs_options_precise = dict(epsilon=1e-3, approx_grad=True,\n",
    "                            m=30, maxls=50)\n",
    "\n",
    "target_list = []\n",
    "time_list = []\n",
    "residuals_list = []\n",
    "area_list = []\n",
    "counter = 0\n",
    "for target in archive: \n",
    "    for time in archive[target]:\n",
    "        counter +=1 \n",
    "\n",
    "for target in archive: \n",
    "    for time in archive[target]:\n",
    "        if target != 'KIC9652680' and target != '51Peg':\n",
    "\n",
    "            spectrum3 = archive[target][time]\n",
    "            wavelength3 = spectrum3['wavelength'][:]\n",
    "            flux3 = spectrum3['flux'][:]\n",
    "            source2 = SimpleSpectrum(wavelength3, flux3, dispersion_unit=u.Angstrom)\n",
    "\n",
    "            spec_band = []\n",
    "            for band, inds in zip(bands_TiO[:first_n_bands], target_slices.wavelength_splits):\n",
    "                target_slice = slice_spectrum(source2, band.min-width*u.Angstrom, band.max+width*u.Angstrom, \n",
    "                                              force_length=abs(np.diff(inds))[0])\n",
    "                target_slice.flux /= target_slice.flux.max()\n",
    "                spec_band.append(target_slice)\n",
    "\n",
    "            source2_slices = concatenate_spectra(spec_band)\n",
    "\n",
    "            bounds = [[-30, 0], [0, 15]] + first_n_bands*[[-2, 2]]\n",
    "            initp = [-0.5, 1] + first_n_bands*[0.0]\n",
    "\n",
    "            result = fmin_l_bfgs_b(chi2, initp, bounds=bounds, \n",
    "                                   args=(target_slices, source1_slices, source2_slices),\n",
    "                                   **bfgs_options_precise)\n",
    "\n",
    "            model, resid = instr_model(target_slices, source1_slices, source2_slices, *result[0])\n",
    "\n",
    "            target_list.append(target)\n",
    "            time_list.append(time)\n",
    "            residuals_list.append(resid)\n",
    "            area_list.append(np.exp(result[0][0]))\n",
    "            print(target, resid, np.exp(result[0][0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_percentile = np.percentile(residuals_list, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HD39587 2016-12-02T09:35:04.980 0.246809249948 0.420588836135\n",
      "HD89744 2017-11-06T09:45:40.470 0.250203466225 0.516621318119\n",
      "HD98230 2017-03-17T06:05:03.709 0.250254675589 0.358736673515\n"
     ]
    }
   ],
   "source": [
    "for target, time, resid, area in zip(target_list, time_list, residuals_list, area_list):\n",
    "    if first_percentile > resid: \n",
    "        print(target, time, resid, area)"
   ]
  },
  {
   "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}