6_Error_Propagation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Propogate Errors\n",
    "\n",
    "This notebook takes you through the steps of how to propogate errors for through the neural network model\n",
    "\n",
    "* required packages: `numpy h5py keras`\n",
    "* data files: \n",
    "    - starnet_cnn.h5\n",
    "    - mean_and_std.npy\n",
    "    - test_data.h5\n",
    "    - apStar_combined_main.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from keras.models import load_model\n",
    "import h5py\n",
    "import tensorflow as tf\n",
    "import time\n",
    "import keras.backend as K\n",
    "import subprocess\n",
    "\n",
    "\n",
    "datadir= \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define path variables for your keras model, denormalization data, and test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model_path = datadir + 'starnet_cnn.h5'\n",
    "denormalization_path = datadir + 'mean_and_std.npy'\n",
    "test_data_path = datadir + 'test_data.h5'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Define functions to:**\n",
    "\n",
    "1. compute the jacobian matrix\n",
    "2. compute the covariance\n",
    "3. compute the variance\n",
    "\n",
    "Note: these functions can be combined into one, but they are separated here to allow users to extract intermediate results for analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calc_jacobian(model,spectra,denormalize=None):\n",
    "        \n",
    "    if denormalize==None:\n",
    "        y_list = tf.unstack(model.output)\n",
    "    else:\n",
    "        y_list = tf.unstack(denormalize(model.output[0]))\n",
    "\n",
    "    J = [tf.gradients(y, model.input) for y in y_list]\n",
    "\n",
    "\n",
    "    jacobian_func = [K.function([model.input, K.learning_phase()], j_) for j_ in J]\n",
    "\n",
    "    jacobian = np.array([jf([spectra,False]) for jf in jacobian_func])[:,0,0,:]\n",
    "    '''\n",
    "    for i in range(len(spectra)):\n",
    "        jacobian = np.array([jf([spectra,False]) for jf in jacobian_func])[:,:,0,:,0]\n",
    "        np.save('temp/temp_jacobian_'+str(i)+'.npy',jacobian)\n",
    "        if i%int(0.1*len(spectra))==0:\n",
    "            print('Jacobians completed: '+str(i))\n",
    "    \n",
    "    for i in range(len(spectra)):\n",
    "        if i==0:\n",
    "            jacobian = np.load('temp/temp_jacobian_'+str(i)+'.npy')\n",
    "        else:\n",
    "            jacobian = np.concatenate((jacobian,np.load('temp/temp_jacobian_'+str(i)+'.npy')))\n",
    "        subprocess.check_output(['rm','temp/temp_jacobian_'+str(i)+'.npy'])\n",
    "    '''\n",
    "    return jacobian\n",
    "\n",
    "def calc_covariance(model,spectra,err_spectra,denormalize=None):\n",
    "    jac_matrix = calc_jacobian(model,spectra,denormalize)\n",
    "    err_spectra[err_spectra > 6] = 0\n",
    "    jac_matrix = np.nan_to_num(jac_matrix)\n",
    "    covariance = np.einsum('ij,jl->il',(jac_matrix*(err_spectra**2)),jac_matrix.T)\n",
    "    return covariance\n",
    "\n",
    "def calc_variance(model,spectra,err_spectra,denormalize=None):\n",
    "    covariance = calc_covariance(model,spectra,err_spectra,denormalize)\n",
    "    return np.diagonal(covariance, offset=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Create a denormalization function **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mean_and_std = np.load(denormalization_path)\n",
    "mean_labels = mean_and_std[0]\n",
    "std_labels = mean_and_std[1]\n",
    "num_labels = mean_and_std.shape[1]\n",
    "\n",
    "def denormalize(lb_norm):\n",
    "    return ((lb_norm*std_labels)+mean_labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Load the StarNet model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = load_model(model_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Load Test Data **\n",
    "\n",
    "The error propagation technique takes some time, so for the purpose of example, we will only use the first 100 spectra in the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test set contains 300 stars\n"
     ]
    }
   ],
   "source": [
    "num_test = 300\n",
    "\n",
    "f = h5py.File(test_data_path, 'r')\n",
    "test_spectra = f['spectrum']\n",
    "test_err_spectra = f['error_spectrum']\n",
    "test_ap_ids = f['Ap_ID'][0:num_test]\n",
    "test_labels = np.column_stack((f['TEFF'][0:num_test],f['LOGG'][0:num_test],f['FE_H'][0:num_test]))\n",
    "\n",
    "print('Test set contains '  + str(len(test_ap_ids))+' stars')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Compute predictions and errors for the test set **\n",
    "\n",
    "**Steps:**\n",
    "1. compute predictions\n",
    "\n",
    "    \\begin{equation}\n",
    "    h_(\\textbf{x},\\textbf{W}) =  h_{1}(\\textbf{x},\\textbf{W}),...,h_{j}(\\textbf{x},\\textbf{W}))\n",
    "    \\end{equation} \n",
    "\n",
    "        j = 3\n",
    "\n",
    "2. compute jacobian matrix\n",
    "\n",
    "    \\begin{equation}\n",
    "    Jac = \\frac{\\partial h_{j}(\\textbf{x},\\textbf{W})}{\\partial \\textbf{x}} =  (\\frac{\\partial h_{j}(\\textbf{x},\\textbf{W})}{\\partial x_{1}},...,\\frac{\\partial h_{j}(\\textbf{x},\\textbf{W})}{\\partial x_{n}})\n",
    "    \\end{equation} \n",
    "\n",
    "        j = 1,...,3\n",
    "\n",
    "        n = 7214\n",
    "\n",
    "3. compute covariance matrix\n",
    "\n",
    "    \\begin{equation}\n",
    "    Cov = Jac \\times \\Delta \\textbf{x}^2 \\times Jac^T\n",
    "    \\end{equation}\n",
    "    \n",
    "\n",
    "4. obtain propagated variance due to error spectrum from the diagonal of the covariance matrix\n",
    "\n",
    "    \\begin{equation}\n",
    "    \\sigma_{\\mathrm{prop}}^2 \\approx diag(Cov)\n",
    "    \\end{equation}\n",
    "    \n",
    "\n",
    "5. determine which region of the label-space the labels are within to obtain the intrinsic scatter in the corresponding bin. These values have been predetermined from training StarNet on synthetic data and applying it to a test set of synthetic data\n",
    "\n",
    "    \\begin{equation}\n",
    "    \\sigma_{\\mathrm{int}}\n",
    "    \\end{equation}\n",
    "    \n",
    "6. combine propagated error with the intrinsic scatter term\n",
    "\n",
    "    \\begin{equation}\n",
    "    \\Delta h_{j} = \\sqrt{\\sigma_{\\mathrm{prop}}^2  + \\sigma_{\\mathrm{int}}^2}\n",
    "    \\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Making predictions and computing propagated variance for 300 spectra\n",
      "\n",
      "1 completed.\n",
      "0.43507194519 seconds elapsed.\n",
      "\n",
      "31 completed.\n",
      "17.3262429237 seconds elapsed.\n",
      "\n",
      "61 completed.\n",
      "47.8530139923 seconds elapsed.\n",
      "\n",
      "91 completed.\n",
      "93.3576600552 seconds elapsed.\n",
      "\n",
      "121 completed.\n",
      "153.251788855 seconds elapsed.\n",
      "\n",
      "151 completed.\n",
      "230.940063953 seconds elapsed.\n",
      "\n",
      "181 completed.\n",
      "327.530734062 seconds elapsed.\n",
      "\n",
      "211 completed.\n",
      "441.528791904 seconds elapsed.\n",
      "\n",
      "241 completed.\n",
      "572.318181038 seconds elapsed.\n",
      "\n",
      "271 completed.\n",
      "721.773078918 seconds elapsed.\n",
      "\n",
      "All 300 completed.\n",
      "884.853347063 seconds elapsed.\n"
     ]
    }
   ],
   "source": [
    "variance = np.zeros((len(test_labels),3))\n",
    "predictions = np.zeros(test_labels.shape)\n",
    "print('Making predictions and computing propagated variance for '+str(len(test_labels))+' spectra')\n",
    "time_start = time.time()\n",
    "for i in range(len(test_labels)):\n",
    "    spectrum = test_spectra[i:i+1]\n",
    "    err_spectrum = test_err_spectra[i:i+1]\n",
    "    variance[i] = calc_variance(model,spectrum,err_spectrum,denormalize)\n",
    "    predictions[i] = denormalize(model.predict(spectrum))\n",
    "    if i%int(0.1*len(test_labels))==0:\n",
    "        print('\\n'+str(i+1)+' completed.\\n'+str(time.time()-time_start)+' seconds elapsed.')\n",
    "print('\\nAll '+str(i+1)+' completed.\\n'+str(time.time()-time_start)+' seconds elapsed.')\n",
    "f.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Create intrinsic scatter arrays (predetermined) **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "scatter_terms = np.array([[  2.85209088e+01,   2.30193645e+01,   2.10676180e+01,\n",
    "          1.91357425e+01,   1.72090644e+01,   1.58693655e+01,\n",
    "          1.52684102e+01,   1.42387830e+01,   1.64239293e+01,\n",
    "          2.18981017e+01],\n",
    "       [  3.86073715e-02,   3.04916170e-02,   2.44161726e-02,\n",
    "          2.25093310e-02,   2.35929675e-02,   2.36922221e-02,\n",
    "          2.58764773e-02,   2.80946934e-02,   3.34534390e-02,\n",
    "          3.56641714e-02],\n",
    "       [  3.90793092e-02,   2.43149947e-02,   2.25292707e-02,\n",
    "          1.81974298e-02,   1.58638867e-02,   1.46142515e-02,\n",
    "          1.36038125e-02,   1.25392930e-02,   1.24740228e-02,\n",
    "          1.53680421e-02]])\n",
    "scatter_ranges = np.array([[  3.50000000e+03,   3.95000000e+03,   4.40000000e+03,\n",
    "          4.85000000e+03,   5.30000000e+03,   5.75000000e+03,\n",
    "          6.20000000e+03,   6.65000000e+03,   7.10000000e+03,\n",
    "          7.55000000e+03,   8.00000000e+03],\n",
    "       [  0.00000000e+00,   5.00000000e-01,   1.00000000e+00,\n",
    "          1.50000000e+00,   2.00000000e+00,   2.50000000e+00,\n",
    "          3.00000000e+00,   3.50000000e+00,   4.00000000e+00,\n",
    "          4.50000000e+00,   5.00000000e+00],\n",
    "       [ -2.50000000e+00,  -2.20000000e+00,  -1.90000000e+00,\n",
    "         -1.60000000e+00,  -1.30000000e+00,  -1.00000000e+00,\n",
    "         -7.00000000e-01,  -4.00000000e-01,  -1.00000000e-01,\n",
    "          2.00000000e-01,   5.00000000e-01]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** assign each spectrum an intrinsic scatter term depending on which region of the parameter-space the prediction lies **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "scatter_errs = np.empty(test_labels.shape)\n",
    "\n",
    "for i in range(scatter_terms.shape[0]):\n",
    "    for j in range(scatter_terms.shape[1]):\n",
    "        current_min = scatter_ranges[i,j]\n",
    "        current_max = scatter_ranges[i,j+1]\n",
    "        current_scatter = scatter_terms[i,j]\n",
    "        index = np.where((test_labels[:,i]>current_min)&(test_labels[:,i]<current_max))[0]\n",
    "        scatter_errs[index,i]=current_scatter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** combine the propagated error (or the square root of the variance) and intrinsic error in quadrature **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "total_errors = np.sqrt(variance+np.square(scatter_errs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean total statistical errors: \n",
      "\n",
      "Teff  :  33.947 K\n",
      "log(g):  0.076 dex\n",
      "[Fe/H]:  0.024 dex\n"
     ]
    }
   ],
   "source": [
    "# label names\n",
    "label_names = ['Teff  ','log(g)','[Fe/H]']\n",
    "units = ['K','dex','dex']\n",
    "\n",
    "mean_err_total = np.mean(total_errors, axis=0)\n",
    "print('Mean total statistical errors: \\n')\n",
    "for i, err in enumerate(mean_err_total):\n",
    "      print(label_names[i]+':  '+\"{0:.3f}\".format(err)+' '+units[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}