3_Train-Network_PaperII.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train Network Paper II\n",
    "In this notebook we lay out the neural network used in the second paper in the series (https://arxiv.org/abs/2102.06230). To create the synthetic spectra used in this work, please run the `Create_Spectra_all.py` code after updating the output region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "import tensorflow as tf\n",
    "from keras.backend import clear_session\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense, InputLayer, Flatten, Dropout\n",
    "from keras.layers.convolutional import Conv1D\n",
    "from keras.layers.convolutional import MaxPooling1D\n",
    "from keras.optimizers import Adam\n",
    "from keras.callbacks import EarlyStopping, ReduceLROnPlateau\n",
    "from keras.utils import to_categorical\n",
    "from sklearn import preprocessing\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "from pickle import dump\n",
    "from keras.backend import clear_session\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn import metrics\n",
    "import statistics\n",
    "from scipy import interpolate\n",
    "from scipy.stats import gaussian_kde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#------------------------ INPUTS ----------------------------#\n",
    "# Set Input Parameters\n",
    "sim_names = ['10k-red-broad/BOND']#, '10k-red/PNe_2020', '10k-red/SNR_2008']\n",
    "ref_dir = '/home/carterrhea/Dropbox/CFHT/Analysis-Paper2/SyntheticData/'\n",
    "syn_num = 30000  # Number of Synthetic Data per simulation\n",
    "ensemble_num = 10  # Number of layers in ensemble"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"# -- Read in Reference Spectra -- #\")\n",
    "# IF RESOLUTION IS SAMPLED USE THIS\n",
    "# This gets us the wavenumbers we are going to interpolate on\n",
    "# We use this spectra because it is a nice looking (can easily see the lines -- not relevant here)\n",
    "# spectra at a resolution of 5000 exactly. We want to keep that sampling :)\n",
    "ref_SN3 = fits.open(ref_dir+'Reference-Spectrum-SN3.fits')[1].data\n",
    "channel = []\n",
    "counts = []\n",
    "for chan in ref_SN3:\n",
    "    channel.append(chan[0])\n",
    "    counts.append(chan[1])\n",
    "wavenumbers_syn_SN3 = channel\n",
    "\n",
    "\n",
    "ref_SN2 = fits.open(ref_dir+'Reference-Spectrum-SN2.fits')[1].data\n",
    "channel = []\n",
    "counts = []\n",
    "for chan in ref_SN2:\n",
    "    channel.append(chan[0])\n",
    "    counts.append(chan[1])\n",
    "wavenumbers_syn_SN2 = channel\n",
    "\n",
    "\n",
    "ref_SN1 = fits.open(ref_dir+'Reference-Spectrum-SN1.fits')[1].data\n",
    "channel = []\n",
    "counts = []\n",
    "for chan in ref_SN1:\n",
    "    channel.append(chan[0])\n",
    "    counts.append(chan[1])\n",
    "wavenumbers_syn_SN1 = channel\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "print(\"# -- Read in SN3 -- #\")\n",
    "# Read in Fits SN3\n",
    "Counts_SN3 = []  # List containing the spectra\n",
    "Param_dict_SN3 = {}\n",
    "sim_ct = 0\n",
    "for sim_name in sim_names:\n",
    "    print('We are on simutlation: %s'%sim_name)\n",
    "    output_dir = ref_dir+sim_name+'/'\n",
    "    plot_dir = ref_dir+sim_name+'/'\n",
    "    for spec_ct in range(syn_num):\n",
    "        if spec_ct%1000 == 0:\n",
    "            print('  We are on spectrum number %i'%spec_ct)\n",
    "        spectrum = fits.open(output_dir+'Spectrum_SN3_%i.fits'%spec_ct)\n",
    "        header = spectrum[0].header\n",
    "        spec = spectrum[1].data\n",
    "        channel = []\n",
    "        counts = []\n",
    "        for chan in spec:\n",
    "            channel.append(chan[0])\n",
    "            counts.append(chan[1])\n",
    "        # interpolate\n",
    "        f = interpolate.interp1d(channel, counts, kind='slinear')\n",
    "        # Get fluxes of interest\n",
    "        coun = f(wavenumbers_syn_SN3[2:-2])  # Need to cut the size just a little bit of the interpolating region otherwise we wont be able to interpolate!\n",
    "        Counts_SN3.append(np.array(coun))\n",
    "        Param_dict_SN3[sim_ct*syn_num+spec_ct] = [header['Halpha'], header['NII6548'], header['NII6583'], header['SII6716'], header['SII6731']]\n",
    "    sim_ct += 1\n",
    "print(\"# -- Read in SN2 -- #\")\n",
    "# Read in Fits SN2\n",
    "Counts_SN2 = []  # List containing the spectra\n",
    "Param_dict_SN2 = {}\n",
    "sim_ct = 0\n",
    "for sim_name in sim_names:\n",
    "    print('We are on simutlation: %s'%sim_name)\n",
    "    output_dir = ref_dir+sim_name+'/'\n",
    "    plot_dir = ref_dir+sim_name+'/'\n",
    "    for spec_ct in range(syn_num):\n",
    "        if spec_ct%1000 == 0:\n",
    "            print('We are on spectrum number %i'%spec_ct)\n",
    "        spectrum = fits.open(output_dir+'Spectrum_SN2_%i.fits'%spec_ct)\n",
    "        header = spectrum[0].header\n",
    "        spec = spectrum[1].data\n",
    "        channel = []\n",
    "        counts = []\n",
    "        for chan in spec:\n",
    "            channel.append(chan[0])\n",
    "            counts.append(chan[1])\n",
    "        # interpolate\n",
    "        f = interpolate.interp1d(channel, counts, kind='slinear')\n",
    "        # Get fluxes of interest\n",
    "        coun = f(wavenumbers_syn_SN2[2:-2])\n",
    "        Counts_SN2.append(np.array(coun))\n",
    "        Param_dict_SN2[sim_ct*syn_num+spec_ct] = [header['OIII5007'],header['OIII4959'], header['Hbeta']]\n",
    "    sim_ct += 1\n",
    "print(\"# -- Read in SN1 -- #\")\n",
    "# Read in Fits SN1\n",
    "Counts_SN1 = []  # List containing the spectra\n",
    "Param_dict_SN1 = {}\n",
    "sim_ct = 0\n",
    "for sim_name in sim_names:\n",
    "    print('We are on simutlation: %s'%sim_name)\n",
    "    output_dir = ref_dir+sim_name+'/'\n",
    "    plot_dir = ref_dir+sim_name+'/'\n",
    "    for spec_ct in range(syn_num):\n",
    "        if spec_ct%1000 == 0:\n",
    "            print('We are on spectrum number %i'%spec_ct)\n",
    "        spectrum = fits.open(output_dir+'Spectrum_SN1_%i.fits'%spec_ct)\n",
    "        header = spectrum[0].header\n",
    "        spec = spectrum[1].data\n",
    "        channel = []\n",
    "        counts = []\n",
    "        for chan in spec:\n",
    "            channel.append(chan[0])\n",
    "            counts.append(chan[1])\n",
    "        # interpolate\n",
    "        f = interpolate.interp1d(channel, counts, kind='slinear')\n",
    "        # Get fluxes of interest\n",
    "        coun = f(wavenumbers_syn_SN1[2:-2])\n",
    "        Counts_SN1.append(np.array(coun))\n",
    "        Param_dict_SN1[sim_ct*syn_num+spec_ct] = [header['OII3726'], header['OII3729']]\n",
    "    sim_ct += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine the Filters for each observation\n",
    "Counts = [np.concatenate((sn1,sn2,sn3), axis=0) for sn1,sn2,sn3 in zip(Counts_SN1,Counts_SN2,Counts_SN3)]\n",
    "# Define values of interest\n",
    "s2_s1 = [np.arcsinh(np.round(val[4]/val[3], 2)) for val in list(Param_dict_SN3.values())]  # s2/s1\n",
    "s1s2_ha = [np.arcsinh(np.round((val[3]+val[4])/val[0], 2)) for val in list(Param_dict_SN3.values())]   # (s1+s2)/ha\n",
    "n2_ha = [np.arcsinh(np.round(val[2]/val[0], 2)) for val in list(Param_dict_SN3.values())]  # n2/ha\n",
    "n2_s1s2 = [np.arcsinh(np.round(val[2]/(val[3]+val[4]), 2)) for val in list(Param_dict_SN3.values())]  # n2/(s1+s2)\n",
    "o3_hb = [np.arcsinh(np.round(val[0]/val[2], 2)) for val in list(Param_dict_SN2.values())]  # O3/hb\n",
    "o2_hb = [np.arcsinh(np.round((O2[0]+O2[1])/hb[2], 2)) for O2,hb in list(zip(Param_dict_SN1.values(),Param_dict_SN2.values()))]  # O2/hb\n",
    "o2_o3 = [np.arcsinh(np.round((O2[0]+O2[1])/O3[0], 2)) for O2,O3 in list(zip(Param_dict_SN1.values(),Param_dict_SN2.values()))]  # O2/O3\n",
    "ha_hb = [np.arcsinh(np.round(ha[0]/hb[2], 2)) for ha,hb in list(zip(Param_dict_SN3.values(),Param_dict_SN2.values()))]  # ha/hb\n",
    "n1_n2 = [np.arcsinh(np.round(val[1]/val[2], 2)) for val in list(Param_dict_SN3.values())]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------------------- DEFINE MODEL -------------------------- #\n",
    "\n",
    "activation = 'relu'  # activation function\n",
    "initializer = 'normal'  # model initializer\n",
    "batch_size = 16  # number of data fed into model at once\n",
    "max_epochs = 25  # maximum number of interations\n",
    "early_stopping_min_delta = 0.0001\n",
    "early_stopping_patience = 4\n",
    "reduce_lr_factor = 0.5\n",
    "reuce_lr_epsilon = 0.009\n",
    "reduce_lr_patience = 2\n",
    "reduce_lr_min = 0.00008\n",
    "loss_function = 'Huber'\n",
    "metrics_ = ['accuracy', 'mae', 'mse']\n",
    "\n",
    "\n",
    "\n",
    "clear_session()\n",
    "def create_model():\n",
    "    model = Sequential([\n",
    "        Dense(units=1024, kernel_initializer=initializer, activation=activation,\n",
    "              kernel_regularizer=l2(0.00005)),\n",
    "        Dropout(0.25),\n",
    "        Dense(units=1024, kernel_initializer=initializer, activation=activation,\n",
    "              kernel_regularizer=l2(0.00005)),\n",
    "        Dropout(0.15),\n",
    "        Dense(units=512, kernel_initializer=initializer, activation=activation,\n",
    "              kernel_regularizer=l2(0.00005)),\n",
    "        Dense(9),\n",
    "    ])\n",
    "    return model\n",
    "\n",
    "optimizer = Adam(learning_rate=0.0004)\n",
    "\n",
    "early_stopping = EarlyStopping(monitor='val_loss', min_delta=early_stopping_min_delta,\n",
    "                                       patience=early_stopping_patience, verbose=2, mode='min')\n",
    "\n",
    "reduce_lr = ReduceLROnPlateau(monitor='loss', factor=0.5, epsilon=reuce_lr_epsilon,\n",
    "                                  patience=reduce_lr_patience, min_lr=reduce_lr_min, mode='min', verbose=2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Data sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----------------------------- TRAIN and VALIDATE ------------------------------ #\n",
    "print(\"# -- Train Network -- #\")\n",
    "# Get number of sucessful cloudy runs\n",
    "syn_num_pass = len(Counts)\n",
    "train_div = int(0.7*syn_num_pass)  # Percent of synthetic data to use as training\n",
    "valid_div = int(0.9*syn_num_pass)  # Percent of synthetic data used as training and validation (validation being the difference between this and train_div)\n",
    "test_div = int(1.0*syn_num_pass)\n",
    "# Set training set\n",
    "TraningSet = np.array(Counts[:train_div])\n",
    "Traning_init_labels = np.array((s2_s1[0:train_div], s1s2_ha[0:train_div], n2_ha[0:train_div], n2_s1s2[0:train_div], o3_hb[0:train_div], o2_hb[0:train_div], o2_o3[0:train_div], ha_hb[0:train_div], n1_n2[0:train_div])).T\n",
    "# Set validation set\n",
    "ValidSet = np.array(Counts[train_div:valid_div])\n",
    "Valid_init_labels = np.array((s2_s1[train_div:valid_div], s1s2_ha[train_div:valid_div], n2_ha[train_div:valid_div], n2_s1s2[train_div:valid_div], o3_hb[train_div:valid_div], o2_hb[train_div:valid_div], o2_o3[train_div:valid_div], ha_hb[train_div:valid_div], n1_n2[train_div:valid_div])).T\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model using deep ensembling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for step in range(ensemble_num):\n",
    "    model = create_model()\n",
    "    model.compile(optimizer=optimizer, loss=loss_function, metrics=metrics_)\n",
    "    model.fit(TraningSet, Traning_init_labels, validation_data=(ValidSet, Valid_init_labels),\n",
    "          epochs=max_epochs,batch_size= batch_size, verbose=2, callbacks=[reduce_lr,early_stopping])\n",
    "    model.save_weights('./ensembling_weights/step_%i'%step)\n",
    "\n",
    "\n",
    "model.save('Sitelle-Line-Ratios-Model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test the model -- this is important for later if you want to implement this network in your own code!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"# -- Test Network -- #\")\n",
    "TestSet = np.array(Counts[valid_div:test_div])\n",
    "TestSetLabels = np.array((s2_s1[valid_div:test_div], s1s2_ha[valid_div:test_div], n2_ha[valid_div:test_div], n2_s1s2[valid_div:test_div], o3_hb[valid_div:test_div], o2_hb[valid_div:test_div], o2_o3[valid_div:test_div], ha_hb[valid_div:test_div], n1_n2[valid_div:test_div])).T\n",
    "# Apply Deep ensembling predictions\n",
    "ensembling_predictions = {}  # {ensemble number: list of predictions}\n",
    "for step in range(ensemble_num):\n",
    "    # Restore the weights\n",
    "    model = create_model()\n",
    "    model.compile(optimizer=optimizer, loss=loss_function, metrics=metrics_)\n",
    "    model.load_weights('./ensembling_weights/step_%i'%step)\n",
    "    predictions = model.predict(TestSet)  # If not testing use predictions=model(TestSet) as it is faster\n",
    "    ensembling_predictions[step] = predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = np.zeros((len(TestSet),9))\n",
    "prediction_prob = np.zeros((len(TestSet),9))  # (test_i, [ratios])\n",
    "prediction_std = np.zeros((len(TestSet),9))\n",
    "\n",
    "for test in range(len(TestSet)):  # Get prediction for each test\n",
    "    ens_pred_final = np.zeros((9,ensemble_num))\n",
    "    for ens in range(ensemble_num):  # get prediction of each network\n",
    "        ens_pred = ensembling_predictions[ens]  # Get ensemble number predictions\n",
    "        ens_test_pred = ens_pred[test]  # Get test prediction for given model in ensemble\n",
    "        ens_pred_final[:,ens] = ens_test_pred  # Sum up the predictions\n",
    "    prediction_prob[test] = ens_pred_final.mean(axis=1)\n",
    "    prediction_std[test] = ens_pred_final.std(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error calculations, plotting, and pickle dumps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------------------------ ERROR CALCULATIONS ------------------#\n",
    "print(\"# -- Error Calculations -- #\")\n",
    "errs_s2_s1 = [100*(prediction_prob[i,0] - TestSetLabels[i,0])/TestSetLabels[i,0] for i in range(len(prediction_prob))]\n",
    "errs_s1s2_ha = [100*(prediction_prob[i,1] - TestSetLabels[i,1])/TestSetLabels[i,1] for i in range(len(prediction_prob))]\n",
    "errs_n2_ha = [100*(prediction_prob[i,2] - TestSetLabels[i,2])/TestSetLabels[i,2] for i in range(len(prediction_prob))]\n",
    "errs_n1_s1s2 = [100*(prediction_prob[i,3] - TestSetLabels[i,3])/TestSetLabels[i,3] for i in range(len(prediction_prob))]\n",
    "errs_o3_hb = [100*(prediction_prob[i,4] - TestSetLabels[i,4])/TestSetLabels[i,4] for i in range(len(prediction_prob))]\n",
    "errs_o2_hb = [100*(prediction_prob[i,5] - TestSetLabels[i,5])/TestSetLabels[i,5] for i in range(len(prediction_prob))]\n",
    "errs_o2_o3 = [100*(prediction_prob[i,6] - TestSetLabels[i,6])/TestSetLabels[i,6] for i in range(len(prediction_prob))]\n",
    "errs_ha_hb = [100*(prediction_prob[i,7] - TestSetLabels[i,7])/TestSetLabels[i,7] for i in range(len(prediction_prob))]\n",
    "errs_n1_n2 = [100*(prediction_prob[i,8] - TestSetLabels[i,8])/TestSetLabels[i,8] for i in range(len(prediction_prob))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fig, axs = plt.subplots(2, 4, figsize=(20,8))\n",
    "sns.distplot(errs_s2_s1, hist=False, kde=True,  kde_kws={\"shade\": True},\n",
    "             bins=20, color = 'darkmagenta',\n",
    "             ax=axs[0,0])\n",
    "#axs[0,0].set_xlim(-5,5)\n",
    "axs[0,0].set_title(r'S[II]$\\lambda$6731/S[II]$\\lambda$6716')# -- std:%.2f'%(np.std(errs_s2_s1)))\n",
    "axs[0,0].set_ylabel('Density', fontweight='bold', fontsize=14)\n",
    "\n",
    "\n",
    "sns.distplot(errs_n2_ha, hist=False, kde=True,  kde_kws={\"shade\": True},\n",
    "             bins=20, color = 'darkmagenta',\n",
    "             ax=axs[0,1])\n",
    "#axs[0,1].set_xlim(-50,50)\n",
    "axs[0,1].set_title(r'N[II]$\\lambda$6584/H$\\alpha(6563)\\AA$')#' -- std:%.2f'%(np.std(errs_n2_ha)))\n",
    "\n",
    "sns.distplot(errs_s1s2_ha, hist=False, kde=True,  kde_kws={\"shade\": True},\n",
    "             bins=20, color = 'darkmagenta',\n",
    "             ax=axs[1,0])\n",
    "#axs[1,0].set_xlim(-50,50)\n",
    "axs[1,0].set_title(r'(S[II]$\\lambda6717$+S[II]$\\lambda6731$)/H$\\alpha(6563)\\AA$')# -- std:%.2f'%(np.std(errs_s1s2_ha)))\n",
    "axs[1,0].set_xlabel('Relative Error (%)', fontweight='bold', fontsize=14)\n",
    "axs[1,0].set_ylabel('Density', fontweight='bold', fontsize=14)\n",
    "\n",
    "sns.distplot(errs_n1_s1s2, hist=False, kde=True,  kde_kws={\"shade\": True},\n",
    "             bins=20, color = 'darkmagenta',\n",
    "             ax=axs[1,1])\n",
    "#axs[1,1].set_xlim(-50,50)\n",
    "axs[1,1].set_title(r'N[II]$\\lambda$6584/(S[II]$\\lambda6717$+S[II]$\\lambda6731$)')#' -- std:%.2f'%(np.std(errs_n1_s1s2)))\n",
    "axs[1,1].set_xlabel('Relative Error (%)', fontweight='bold', fontsize=14)\n",
    "\n",
    "sns.distplot(errs_o2_hb, hist=False, kde=True,  kde_kws={\"shade\": True},\n",
    "             bins=20, color = 'darkmagenta',\n",
    "             ax=axs[0,2])\n",
    "#axs[1,1].set_xlim(-50,50)\n",
    "axs[0,2].set_title(r'O[II]$\\lambda$3726/H$\\beta$')#' -- std:%.2f'%(np.std(errs_n1_s1s2)))\n",
    "#axs[0,2].set_xlabel('Relative Error (%)', fontweight='bold', fontsize=14)\n",
    "\n",
    "sns.distplot(errs_o3_hb, hist=False, kde=True,  kde_kws={\"shade\": True},\n",
    "             bins=20, color = 'darkmagenta',\n",
    "             ax=axs[1,2])\n",
    "#axs[1,1].set_xlim(-50,50)\n",
    "axs[1,2].set_title(r'OIII$\\lambda$5007/H$\\beta$')#' -- std:%.2f'%(np.std(errs_n1_s1s2)))\n",
    "axs[1,2].set_xlabel('Relative Error (%)', fontweight='bold', fontsize=14)\n",
    "\n",
    "sns.distplot(errs_o2_o3, hist=False, kde=True,  kde_kws={\"shade\": True},\n",
    "             bins=20, color = 'darkmagenta',\n",
    "             ax=axs[0,3])\n",
    "#axs[1,1].set_xlim(-50,50)\n",
    "axs[0,3].set_title(r'OII$\\lambda$3726/OIII$\\lambda$5007')#' -- std:%.2f'%(np.std(errs_n1_s1s2)))\n",
    "#axs[0,3].set_xlabel('Relative Error (%)', fontweight='bold', fontsize=14)\n",
    "\n",
    "sns.distplot(errs_ha_hb, hist=False, kde=True,  kde_kws={\"shade\": True},\n",
    "             bins=20, color = 'darkmagenta',\n",
    "             ax=axs[1,3])\n",
    "#axs[1,1].set_xlim(-50,50)\n",
    "axs[1,3].set_title(r'H$\\alpha$/H$\\beta$')#' -- std:%.2f'%(np.std(errs_n1_s1s2)))\n",
    "axs[1,3].set_xlabel('Relative Error (%)', fontweight='bold', fontsize=14)\n",
    "\n",
    "\n",
    "\n",
    "plt.savefig(plot_dir+'Ratio-Errors.pdf', type='pdf')\n",
    "\n",
    "error_txt = open(plot_dir+'err_std.txt', 'w+')\n",
    "errs = [errs_s2_s1, errs_n2_ha, errs_s1s2_ha, errs_n1_s1s2, errs_o2_o3, errs_o3_hb, errs_o2_hb, errs_ha_hb, errs_n1_n2]\n",
    "errs_name = ['errs_s2_s1', 'errs_n2_ha', 'errs_s1s2_ha', 'errs_n1_s1s2', 'errs_o2_o3', 'errs_o3_hb', 'errs_o2_hb', 'errs_ha_hb', 'errs_n1_n2']\n",
    "for err,err_name in zip(errs,errs_name):\n",
    "    std_err = np.std(err)\n",
    "    mean_ = np.median(err)\n",
    "    print(\"The standard deviation for %s is %.2f\"%(str(err_name), std_err))\n",
    "    print(\"The mean for %s is %.2f\"%(str(err_name), mean_))\n",
    "    error_txt.write(\"The standard deviation for %s is %.2f\\n\"%(str(err_name), std_err))\n",
    "    error_txt.write(\"The mean for %s is %.2f\\n\"%(str(err_name), mean_))\n",
    "error_txt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pickle.dump(errs_s2_s1, open('errs_s2_s1.pkl', 'wb'))\n",
    "pickle.dump(errs_n2_ha, open('errs_n2_ha.pkl', 'wb'))\n",
    "pickle.dump(errs_s1s2_ha, open('errs_s1s2_ha.pkl', 'wb'))\n",
    "pickle.dump(errs_n1_s1s2, open('errs_n1_s1s2.pkl', 'wb'))\n",
    "pickle.dump(errs_o2_o3, open('errs_o2_o3.pkl', 'wb'))\n",
    "pickle.dump(errs_o3_hb, open('errs_o3_hb.pkl', 'wb'))\n",
    "pickle.dump(errs_o2_hb, open('errs_o2_hb.pkl', 'wb'))\n",
    "pickle.dump(errs_ha_hb, open('errs_ha_hb.pkl', 'wb'))\n",
    "pickle.dump(errs_n1_n2, open('errs_n1_n2.pkl', 'wb'))"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}