HPF CALS data: ratio of twilight flats for fiber throughput calibration

03-CALS_can_we_flat_field_with_twilight_flats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5c558293",
   "metadata": {},
   "source": [
    "# Can we flat field with twilight flats\n",
    "\n",
    "Monday, August 16, 2021\n",
    "\n",
    "We downloaded many GB of CALS data from TACC today."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5548b8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from muler.hpf import HPFSpectrumList, HPFSpectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e8f899c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import pandas as pd\n",
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc9c4354",
   "metadata": {},
   "outputs": [],
   "source": [
    "fns = sorted(glob.glob('../../data/HPF/goldilocks/CALS/Goldilocks_20*.spectra.fits'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "947f03c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(fns)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1c99ae5",
   "metadata": {},
   "source": [
    "We have about 2000 calibration files.  Let's inspect a random one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fd44847",
   "metadata": {},
   "outputs": [],
   "source": [
    "fn = fns[1234]\n",
    "fn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fc10ed8",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec = HPFSpectrum(file=fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8f6ff0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "info_dict = {}\n",
    "for fn in tqdm(fns):\n",
    "    spec = HPFSpectrum(file=fn)\n",
    "    info_dict[fn] = spec.meta['header']['OBJECT']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64e001a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "objects = pd.Series(info_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42068696",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_objects = objects.to_frame(name='object').reset_index().rename(columns={'index':'fn'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ceb0878",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.max_rows = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b4fe31a",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_objects.object.value_counts().head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f67bf2d4",
   "metadata": {},
   "source": [
    "### Construct the flat field spectrum\n",
    "We want the spectra that say `Quartz FCU`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7f7a7a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format='retina'\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21b0782c",
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_fns = df_objects[df_objects.object == 'Quartz FCU'].fn.values\n",
    "spec_list = HPFSpectrumList.read(flat_fns[10])\n",
    "spec = spec_list[10]\n",
    "beta = spec / spec.sky\n",
    "\n",
    "ax = spec.plot()\n",
    "spec.sky.plot(ax=ax)\n",
    "plt.ylim(0, 5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e00b86b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = spec.sky.divide(spec, compare_wcs='first_found')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ed5440f",
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_list = []\n",
    "flux_list = []\n",
    "for fn in tqdm(flat_fns):\n",
    "    spec_list = HPFSpectrumList.read(fn)\n",
    "    spec = spec_list[10]\n",
    "    beta = spec / spec.sky\n",
    "    wave_list.append(beta.wavelength.value)\n",
    "    flux_list.append(beta.flux.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd404282",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import binned_statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dda6398",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8571e89d",
   "metadata": {},
   "outputs": [],
   "source": [
    "waves = np.hstack(wave_list)\n",
    "fluxes = np.hstack(flux_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fcab796",
   "metadata": {},
   "outputs": [],
   "source": [
    "func25 = lambda x: np.nanpercentile(x, 97.5)\n",
    "func75 = lambda x: np.nanpercentile(x, 2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "815c2a1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_bins = binned_statistic(waves, waves, statistic='mean', bins=2048).statistic\n",
    "mean_beta = binned_statistic(waves, fluxes, statistic='mean', bins=2048).statistic\n",
    "beta_25 = binned_statistic(waves, fluxes, statistic=func25, bins=2048).statistic\n",
    "beta_75 = binned_statistic(waves, fluxes, statistic=func75, bins=2048).statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af4840a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['axes.facecolor'] = 'white'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e884f010",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(flat_fns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec5b4938",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(wave_bins, mean_beta)\n",
    "plt.fill_between(wave_bins, beta_25, beta_75, alpha=0.5, label='$95\\%$ region')\n",
    "plt.axhline(1.0, linestyle='dashed', color='k')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('$\\lambda \\;(\\AA)$', fontsize=12)\n",
    "plt.ylabel(r'$\\beta \\equiv f_{\\star}/f_{\\mathrm{sky}}$', fontsize=12);\n",
    "plt.savefig('../../figures/HPF_flat_field_demo.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2321dcea",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_list.normalize().plot(ylo=0, yhi=10);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee754dc5",
   "metadata": {},
   "source": [
    "### Construct the twilight flat field spectrum\n",
    "We want the spectra that say `twilight`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74998ae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils.manipulation import FluxConservingResampler, LinearInterpolatedResampler, SplineInterpolatedResampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9face2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_fns = df_objects[df_objects.object == 'twilight'].fn.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "675adad2",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_list = HPFSpectrumList.read(flat_fns[25])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "558069e7",
   "metadata": {},
   "source": [
    "We can't just divide the spectra because their slightly different wavelength values make a difference.  We need to resample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42654db4",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec = spec_list[25].normalize()\n",
    "\n",
    "resampler = LinearInterpolatedResampler()\n",
    "new_sky = resampler(spec.sky, spec.spectral_axis) \n",
    "\n",
    "beta = spec / new_sky"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb36a1e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = spec.plot()\n",
    "ax.plot(new_sky.wavelength.value, new_sky.flux.value)\n",
    "plt.ylim(0, 1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74eb80f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#spec.meta['header']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1ccd4d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(beta.wavelength.value, beta.flux.value, label='2021 Aug 8 Twilight Flat')\n",
    "plt.ylim(0.88, 1.02)\n",
    "plt.axhline(1.0, linestyle='dashed', color='k')\n",
    "plt.xlabel('$\\lambda \\;(\\AA)$', fontsize=12)\n",
    "plt.ylabel(r'$\\beta \\equiv f_{\\star}/f_{\\mathrm{sky}}$', fontsize=12);\n",
    "plt.legend()\n",
    "plt.savefig('../../figures/HPF_twilight_2021Aug8.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec82e20e",
   "metadata": {},
   "outputs": [],
   "source": [
    "resampler = LinearInterpolatedResampler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09448002",
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_list = []\n",
    "flux_list = []\n",
    "this_wave = []\n",
    "this_flux = []\n",
    "date_list = []\n",
    "for fn in tqdm(flat_fns):\n",
    "    spec_list = HPFSpectrumList.read(fn)\n",
    "    for j in range(0, 28):\n",
    "        spec = spec_list[j]\n",
    "        new_sky = resampler(spec.sky, spec.spectral_axis) \n",
    "        beta = spec / new_sky\n",
    "        this_wave.append(beta.wavelength.value)\n",
    "        this_flux.append(beta.flux.value)\n",
    "\n",
    "    wave_list.append(np.hstack(this_wave))\n",
    "    flux_list.append(np.hstack(this_flux))\n",
    "    date_list.append(spec.meta['header']['DATE'][0:10])\n",
    "    this_wave = []\n",
    "    this_flux = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e228de7",
   "metadata": {},
   "outputs": [],
   "source": [
    "waves = np.vstack(wave_list)\n",
    "wave_bins = waves.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc6aee8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "one_micron_indices = np.abs(wave_bins - 10_500) < 10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d29ff4d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes = np.vstack(flux_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de92289c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8fd8e88",
   "metadata": {},
   "outputs": [],
   "source": [
    "by_order_fluxes = fluxes.reshape((26, 28, 2048))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41685acf",
   "metadata": {},
   "outputs": [],
   "source": [
    "by_order_wavelengths = waves.reshape((26, 28, 2048))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74cbf27e",
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_groups = by_order_wavelengths.mean(axis=2)\n",
    "flux_groups = np.nanstd(by_order_fluxes, axis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8f9cbdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "bi = flux_groups[:, -1] > 0.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ae1e3d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes = fluxes[~bi]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7033c70",
   "metadata": {},
   "outputs": [],
   "source": [
    "dates = np.array(date_list)[~bi]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e91f61e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "pal = sns.color_palette(\"rocket\", n_colors=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7adcd050",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "for j in range(25):\n",
    "    plt.plot(wave_bins, fluxes[j], alpha=0.1, color=pal[j], label=None)\n",
    "for j in range(25):\n",
    "    plt.plot(wave_bins, fluxes[j]*0.0, alpha=1, color=pal[j], label=dates[j])\n",
    "plt.ylim(0.8, 1.2)\n",
    "plt.axhline(1.0, linestyle='dashed', color='k')\n",
    "plt.axhline(0.93, linestyle='dotted', color='k', label='93%')\n",
    "plt.legend(ncol=7);\n",
    "#plt.xlim(10440, 10_950);\n",
    "plt.xlabel('$\\lambda \\;(\\AA)$', fontsize=12)\n",
    "plt.ylabel(r'$\\beta \\equiv$ (target/sky)', fontsize=12)\n",
    "plt.title('Ratio of HPF Twilight Flat Spectra from target and sky fibers')\n",
    "plt.savefig('../../figures/HPF_twilight_ratio_overview.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c85020d4",
   "metadata": {},
   "source": [
    "Nice!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfb36fc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "for j in range(25):\n",
    "    plt.plot(wave_bins, fluxes[j], alpha=0.1, color=pal[j], label=None)\n",
    "for j in range(25):\n",
    "    plt.plot(wave_bins, fluxes[j]*0.0, alpha=1, color=pal[j], label=dates[j])\n",
    "plt.ylim(0.91, 0.97)\n",
    "plt.axhline(1.0, linestyle='dashed', color='k')\n",
    "plt.axhline(0.93, linestyle='dotted', color='k', label='93%')\n",
    "plt.legend(ncol=7);\n",
    "plt.xlim(10440, 10_950);\n",
    "plt.xlabel('$\\lambda \\;(\\AA)$', fontsize=12)\n",
    "plt.ylabel(r'$\\beta \\equiv$ (target/sky)', fontsize=12)\n",
    "plt.title('Ratio of HPF Twilight Flat Spectra from target and sky fibers')\n",
    "plt.savefig('../../figures/HPF_twilight_ratio_over_time.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f59a9013",
   "metadata": {},
   "outputs": [],
   "source": [
    "#! open ../../figures/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c418968",
   "metadata": {},
   "outputs": [],
   "source": [
    "#by_order_fluxes = by_order_fluxes / np.median(by_order_fluxes, axis=2)[:,:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62d06a0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f28e94b",
   "metadata": {},
   "source": [
    "## Construct the binned statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08cda7f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "waves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fd953b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "waves = waves[~bi]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e89046a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_edges = wave_bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "066ca502",
   "metadata": {},
   "outputs": [],
   "source": [
    "#wave_bins_alt = binned_statistic(waves, waves, statistic='mean', bins=2048).statistic\n",
    "mean_beta = binned_statistic(waves.reshape(-1), fluxes.reshape(-1), statistic=np.nanmedian, bins=bin_edges).statistic\n",
    "beta_25 = binned_statistic(waves.reshape(-1), fluxes.reshape(-1), statistic=func25, bins=bin_edges).statistic\n",
    "beta_75 = binned_statistic(waves.reshape(-1), fluxes.reshape(-1), statistic=func75, bins=bin_edges).statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae4d8e1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.plot(wave_bins[1:], mean_beta)\n",
    "plt.fill_between(wave_bins[1:], beta_25, beta_75, alpha=0.5, label='$95\\%$ region')\n",
    "plt.axhline(1.0, linestyle='dashed', color='k')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('$\\lambda \\;(\\AA)$', fontsize=12)\n",
    "plt.ylabel(r'$\\beta \\equiv f_{\\star}/f_{\\mathrm{sky}}$', fontsize=12);\n",
    "plt.ylim(0.9, 1.05)\n",
    "#plt.savefig('../../figures/HPF_twilight_flat_demo.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7cc7130",
   "metadata": {},
   "source": [
    "That's nice and all, but let's smooth this out..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aafeeee4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import celerite2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0ce980f",
   "metadata": {},
   "outputs": [],
   "source": [
    "wave1d = waves.reshape(-1)\n",
    "flux1d = fluxes.reshape(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36b72185",
   "metadata": {},
   "outputs": [],
   "source": [
    "from celerite2 import terms\n",
    "\n",
    "# Quasi-periodic term\n",
    "term1 = terms.SHOTerm(sigma=0.001, rho=1000.0, tau=1000.0)\n",
    "\n",
    "# Non-periodic component\n",
    "term2 = terms.SHOTerm(sigma=0.01, rho=1000, Q=0.25)\n",
    "kernel = term1 + term2\n",
    "\n",
    "# Setup the GP\n",
    "gp = celerite2.GaussianProcess(kernel, mean=0.94)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "487ae91a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted_inds = np.argsort(wave_bins[0:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28b74cd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "w1d = wave_bins[0:-1][sorted_inds]\n",
    "f1d = mean_beta[sorted_inds]\n",
    "s1d = (beta_75-beta_25)[sorted_inds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f43aa6b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "gi = (f1d == f1d) & (s1d == s1d) & (f1d < 0.975) & (f1d > 0.92)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "586cd7b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "w1d = w1d[gi]\n",
    "f1d = f1d[gi]\n",
    "s1d = s1d[gi]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e38ea9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "s1d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa35acf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "gp.compute(w1d, yerr=0.0001)\n",
    "\n",
    "print(\"Initial log likelihood: {0}\".format(gp.log_likelihood(f1d)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95837ee2",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.plot(w1d, gp.sample(include_mean=True))\n",
    "#plt.xlim(10_000, 10_500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4eff9b",
   "metadata": {},
   "source": [
    "# Iterative Sigma clipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b492b479",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_sorted_inds = np.argsort(wave1d)\n",
    "\n",
    "wave1d_sorted, flux1d_sorted = wave1d[all_sorted_inds], flux1d[all_sorted_inds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70708a6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "flux_smooth = gp.predict(f1d, wave1d_sorted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ef6c10a",
   "metadata": {},
   "outputs": [],
   "source": [
    "residual = flux1d_sorted - flux_smooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a5cffaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_bins_resid = binned_statistic(wave1d_sorted, wave1d_sorted, statistic='mean', bins=bin_edges).statistic\n",
    "resid_std = binned_statistic(wave1d_sorted, residual, statistic=np.nanstd, bins=bin_edges).statistic\n",
    "resid_median = binned_statistic(wave1d_sorted, residual, statistic=np.nanmedian, bins=bin_edges).statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fc7a861",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.plot(wave1d_sorted, flux1d_sorted, '.',alpha=0.01, ms=1)\n",
    "plt.plot(wave1d_sorted, flux_smooth, '.',alpha=0.01, ms=1)\n",
    "#plt.fill_between(wave_bins[1:], beta_25, beta_75, alpha=0.5, label='$95\\%$ region')\n",
    "plt.axhline(1.0, linestyle='dashed', color='k')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('$\\lambda \\;(\\AA)$', fontsize=12)\n",
    "plt.ylabel(r'$\\beta \\equiv f_{\\star}/f_{\\mathrm{sky}}$', fontsize=12);\n",
    "#plt.ylim(0, 0.05)\n",
    "plt.ylim(0.9, 1.0)\n",
    "#plt.savefig('../../figures/HPF_twilight_flat_demo.png', dpi=300, bbox_inches='tight')\n",
    "#plt.xlim(11_000, 11_500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f4721a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.plot(wave1d_sorted, flux1d_sorted, '.',alpha=0.03, ms=2)\n",
    "plt.plot(wave1d_sorted, flux_smooth, '.',alpha=0.03, ms=2)\n",
    "#plt.fill_between(wave_bins[1:], beta_25, beta_75, alpha=0.5, label='$95\\%$ region')\n",
    "plt.axhline(1.0, linestyle='dashed', color='k')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('$\\lambda \\;(\\AA)$', fontsize=12)\n",
    "plt.ylabel(r'$\\beta \\equiv f_{\\star}/f_{\\mathrm{sky}}$', fontsize=12);\n",
    "#plt.ylim(0, 0.05)\n",
    "plt.ylim(0.9, 1.0)\n",
    "#plt.savefig('../../figures/HPF_twilight_flat_demo.png', dpi=300, bbox_inches='tight')\n",
    "plt.xlim(11_000, 12_000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4804b2d",
   "metadata": {},
   "source": [
    "Looks good!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd25a9f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "noisy_inds = (np.abs(residual) > 0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0347481d",
   "metadata": {},
   "outputs": [],
   "source": [
    "good_inds = ~noisy_inds & (flux1d_sorted==flux1d_sorted) & (residual == residual)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9680a78",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.plot(wave_bins_resid, resid_median)\n",
    "#plt.fill_between(wave_bins[1:], beta_25, beta_75, alpha=0.5, label='$95\\%$ region')\n",
    "plt.axhline(1.0, linestyle='dashed', color='k')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('$\\lambda \\;(\\AA)$', fontsize=12)\n",
    "plt.ylabel(r'$\\beta \\equiv f_{\\star}/f_{\\mathrm{sky}}$', fontsize=12);\n",
    "#plt.ylim(0, 0.05)\n",
    "plt.ylim(-0.02, 0.02)\n",
    "#plt.savefig('../../figures/HPF_twilight_flat_demo.png', dpi=300, bbox_inches='tight')\n",
    "#plt.xlim(11_000, 11_500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1343cb5",
   "metadata": {},
   "source": [
    "refine_again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a23bfdc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = wave1d_sorted[good_inds]\n",
    "y = flux1d_sorted[good_inds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f22c7947",
   "metadata": {},
   "outputs": [],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ec4b7af",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "gp.compute(x, yerr=0.0001)\n",
    "\n",
    "print(\"Initial log likelihood: {0}\".format(gp.log_likelihood(y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5eefffa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "gp.sample().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a39c9e2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "\n",
    "plt.plot(x, y, '.', alpha=0.01, ms=1)\n",
    "plt.plot(x, gp.predict(y, x), '.', alpha=0.01, ms=1)\n",
    "plt.plot(x, gp.sample(), '.', alpha=0.01, ms=1)\n",
    "plt.ylim(0.9, 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7374a84e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "\n",
    "plt.plot(x, y, '.', alpha=0.02, ms=2)\n",
    "plt.plot(x, gp.predict(y, x), '.', alpha=0.02, ms=2)\n",
    "plt.ylim(0.9, 1.0)\n",
    "plt.xlim(10_500, 10_950)\n",
    "plt.xlabel('$\\lambda \\; (\\AA)$')\n",
    "plt.ylabel(r'$\\beta$')\n",
    "plt.savefig('../../figures/HPF_twilight_ratio_smoothed.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c7ccdb3",
   "metadata": {},
   "source": [
    "Perfect!  Let's use that!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76014df4",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.nanmedian(resid_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fd48063",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.plot(wave_bins[1:], resid_std)\n",
    "plt.ylim(0, 0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3b371e8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1604193",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_beta_estimator = gp.predict(y, wave_bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a7f634f",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_beta_estimator.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b046d06f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af54ce06",
   "metadata": {},
   "outputs": [],
   "source": [
    "std_dev_out = np.hstack([resid_std, 0.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d0c61d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "nan_vals = (std_dev_out !=std_dev_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a6debfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "std_dev_out[nan_vals] = 0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9935659",
   "metadata": {},
   "outputs": [],
   "source": [
    "std_dev_out[std_dev_out>0.8] = 0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d45bd7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_out = {'wave_Ang':wave_bins, 'beta_estimator':final_beta_estimator, 'uncertainty':std_dev_out}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7cde997",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_out = pd.DataFrame(dict_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5729500",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ace1d20",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_out.to_csv('../../../muler/src/muler/templates/HPF_sci_to_sky_ratio_beta.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e01ec9d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}