example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import desispec\n",
    "import desitarget\n",
    "import desimodel.io\n",
    "import pickle\n",
    "import fitsio\n",
    "from astropy.table import Table\n",
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qso_templates import qsoproc\n",
    "from qso_templates import outlier_detection\n",
    "from qso_templates import run_hmf\n",
    "from qso_templates import empca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rest_offset = 2.6\n",
    "nbins = 13500\n",
    "min_per_lam=50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pre-process spectra ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "fluxes, ivars, qsocat = qsoproc.efficient_pipeline(100, nkeep=2, rlb=rest_offset, nbins=nbins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxes.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(qsocat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "with open(\"data/final_data2/fluxes2.pkl\",\"wb\") as file:\n",
    "    pickle.dump(fluxes,file)\n",
    "    \n",
    "with open(\"data/final_data2/ivars2.pkl\", \"wb\") as file:\n",
    "    pickle.dump(ivars, file)\n",
    "    \n",
    "qsocat.write('data/final_data2/qsocat_selected2.fits', format='fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(fluxes, aspect='auto', vmin=0, vmax=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quality cuts ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('data/final_data2/fluxes2.pkl',\"rb\") as file:\n",
    "    fluxes = pickle.load(file)\n",
    "\n",
    "with open('data/final_data2/ivars2.pkl',\"rb\") as file:\n",
    "    ivars = pickle.load(file)\n",
    "\n",
    "qsocat = Table(fits.getdata('data/final_data2/qsocat_selected2.fits'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Removes wavelengths with < min_spectra nonzero observations\n",
    "fluxes, ivars = outlier_detection.enough_obs(fluxes, ivars, min_spectra=min_per_lam)\n",
    "print(fluxes.shape)\n",
    "imshow(ivars, aspect='auto', vmin=0, vmax=2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "old_dist = outlier_detection.calc_chisq_dist(fluxes, ivars)\n",
    "outlier_detection.plot_chisq_dist(old_dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_chisq_dist(old_dist, hist=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_f, new_i, new_dist, qsocat = outlier_detection.outlier_detection(fluxes, ivars, qsocat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outlier_detection.plot_chisq_dist(new_dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outlier_detection.plot_chisq_dist(new_dist, hist=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(new_f.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/final_data2/fluxes2_good.pkl\",\"wb\") as file:\n",
    "    pickle.dump(new_f,file)\n",
    "with open(\"data/final_data2/ivars2_good.pkl\",\"wb\") as file:\n",
    "    pickle.dump(new_i,file)\n",
    "qsocat.write('data/final_data2/qsocat_selected2_good.fits', format='fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HMF-factorization ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/final_data2/fluxes2_good.pkl\",\"rb\") as file:\n",
    "    fluxes = pickle.load(file)\n",
    "    \n",
    "with open(\"data/final_data2/ivars2_good.pkl\",\"rb\") as file:\n",
    "    ivars = pickle.load(file)\n",
    "\n",
    "fluxes.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f2, i2 = run_hmf.enough_obs(fluxes, ivars, min_per_lam)\n",
    "f2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/final_data2/final_fluxes.pkl\",\"wb\") as file:\n",
    "    pickle.dump(f2, file)\n",
    "    \n",
    "with open(\"data/final_data2/final_ivars.pkl\",\"wb\") as file:\n",
    "    pickle.dump(i2, file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow(i2, aspect='auto', vmin=0, vmax=2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get and store loglams\n",
    "bin_range = np.arange(f2.shape[1])\n",
    "rest_loglam_diffs = bin_range * 0.0001\n",
    "rest_loglams = np.array(rest_loglam_diffs + rest_offset)\n",
    "rest_loglams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rest_loglams.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelengths = np.power(10, rest_loglams)\n",
    "wavelengths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/final_data2/final_rest_loglams.pkl\",\"wb\") as file:\n",
    "    pickle.dump(rest_loglams, file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "V, C, chistats = run_hmf.hmf_weighted(fluxes, ivars, 5, num_iter=5,)\n",
    "\n",
    "with open(\"data/final_data2/model_V.pkl\",\"wb\") as file:\n",
    "    pickle.dump(V,file)\n",
    "    \n",
    "with open(\"data/final_data2/model_C.pkl\", \"wb\") as file:\n",
    "    pickle.dump(C, file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/final_data2/model_V.pkl\",\"rb\") as file:\n",
    "    V = pickle.load(file)\n",
    "    \n",
    "with open(\"data/final_data2/model_C.pkl\", \"rb\") as file:\n",
    "    C = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_hmf = np.dot(V, C).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_hmf.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/final_data2/model_hmf.pkl\",\"wb\") as file:\n",
    "    pickle.dump(model_hmf,file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EMPCA "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract model from file\n",
    "with open(\"data/final_data2/model_hmf.pkl\",\"rb\") as file:\n",
    "    model_hmf = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create uniform weights to go along with it\n",
    "weights = np.ones(model_hmf.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(np.arange(model_hmf.shape[1]), np.average(model_hmf, axis=0, weights=weights));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_model = empca(data=model_hmf, weights=weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/final_data2/final_model.pkl\",\"wb\") as file:\n",
    "    pickle.dump(final_model, file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/final_data2/rest_loglams.pkl\",\"rb\") as file:\n",
    "    rest_loglams = pickle.load(file)\n",
    "    \n",
    "wavelengths = 10**rest_loglams"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**These are the plots of all the eigenvectors:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(wavelengths, final_model.eigvec[0]);\n",
    "xlim(4500, 5500);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(wavelengths, final_model.eigvec[1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(wavelengths, final_model.eigvec[2]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(wavelengths, final_model.eigvec[3]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(wavelengths, final_model.eigvec[4]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0\n",
    "plot(model_hmf.data[i], alpha=0.6);\n",
    "plot(final_model.model[i])"
   ]
  }
 ],
 "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.7.4"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}