RNN.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from armored.models import *\n",
    "from armored.preprocessing import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from scipy.stats import linregress\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "params = {'legend.fontsize': 18,\n",
    "          'figure.figsize': (8, 6),\n",
    "         'axes.labelsize': 24,\n",
    "         'axes.titlesize':24,\n",
    "         'axes.linewidth':5,\n",
    "         'xtick.labelsize':20,\n",
    "         'ytick.labelsize':20}\n",
    "plt.rcParams.update(params)\n",
    "plt.style.use('seaborn-colorblind')\n",
    "plt.rcParams['pdf.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import default_rng\n",
    "rng = default_rng(seed = 123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define simulation parameters and import full dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of trials \n",
    "n_trials = 30\n",
    "\n",
    "# number of dtl cycles \n",
    "n_dtl  = 5\n",
    "\n",
    "# define number of initial samples to train on\n",
    "n_init = 5\n",
    "\n",
    "# number of samples for next experiment \n",
    "n_test = 5\n",
    "\n",
    "# number of species in model\n",
    "n_s = 5\n",
    "\n",
    "# number of resources\n",
    "n_r = 7\n",
    "\n",
    "# define all system variables \n",
    "species = ['s'+str(i+1) for i in range(n_s)]\n",
    "outputs = ['product']\n",
    "sys_var = species + outputs\n",
    "\n",
    "# define parameters in the objective function\n",
    "obj_params = ['volume']\n",
    "\n",
    "# define subset of controls\n",
    "controls = ['rf'+str(i+1) for i in range(n_r)] + ['feed']\n",
    "system_variables = species + outputs + controls\n",
    "\n",
    "# define an objective function (product of endpoint volume * endpoint product concentration)\n",
    "# where product concentration is the last column of the predicted output\n",
    "objective = lambda pred, vol: pred[-1, -1]*vol[-1] \n",
    "\n",
    "# import data \n",
    "main_df = pd.read_csv(\"Data/reactor_ubiome.csv\")\n",
    "all_exp_names = main_df.Experiments.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# determine random sets of initial experiments\n",
    "initial_exps = [rng.choice(np.unique(all_exp_names), n_init, replace=False) for _ in range(n_trials)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to compute sum of squares error \n",
    "def sse(a, b):\n",
    "    return np.sum((a-b)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loop over each trial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# init dataframe that stores DTL information\n",
    "dtl_df = pd.DataFrame()\n",
    "dtl_df_R = pd.DataFrame()\n",
    "dtl_df_sse = pd.DataFrame()\n",
    "elapsed_time = []\n",
    "\n",
    "for trial in range(n_trials):\n",
    "\n",
    "    # format data \n",
    "    main_data, main_obj_params, unique_exp_names, N_total = format_data(main_df, species, outputs, controls, obj_params=obj_params)\n",
    "\n",
    "    # keep track of objective \n",
    "    objective_found = []\n",
    "\n",
    "    # choose random set of training samples\n",
    "    train_df = main_df.iloc[np.in1d(all_exp_names, initial_exps[trial])].copy()\n",
    "    train_data, train_obj_params, new_experiments, N = format_data(train_df, species, outputs, controls, obj_params=obj_params)\n",
    "\n",
    "    # remove training samples from main dataset so that they're not selected more than once\n",
    "    train_inds = np.in1d(unique_exp_names, new_experiments)\n",
    "    main_data = main_data[~train_inds]\n",
    "    main_obj_params  = main_obj_params[~train_inds]\n",
    "    unique_exp_names = unique_exp_names[~train_inds]\n",
    "\n",
    "    # compute objectives\n",
    "    target_found = []\n",
    "    for sample, train_obj_param in zip(train_data, train_obj_params):\n",
    "        target_found.append(objective(sample[:, :len(sys_var)], train_obj_param))\n",
    "    target_found = np.array(target_found)\n",
    "    objective_found.append(np.max(target_found))\n",
    "    objective_rval = []\n",
    "    objective_sse  = []\n",
    "\n",
    "    # Search over full factorial and update model\n",
    "    for dtl in range(n_dtl):\n",
    "        print(f\"Running trial {trial+1}, cycle {dtl+1}\")\n",
    "\n",
    "        # scale train and design space data\n",
    "        scaler = ZeroMaxScaler().fit(train_data)\n",
    "        train_data = scaler.transform(train_data)\n",
    "        main_data_scaled = scaler.transform(main_data)\n",
    "\n",
    "        # fit model \n",
    "        brnn = RNN(n_species=n_s, n_metabolites=1, n_controls=len(controls), n_hidden=16, N=N)\n",
    "        brnn.fit(train_data)\n",
    "        \n",
    "        # assess prediction performance of end-point product\n",
    "        pred, stdv, cov = brnn.predict(main_data_scaled)\n",
    "        pred = scaler.inverse_transform(pred)\n",
    "        stdv = scaler.inverse_transform(stdv)\n",
    "        rvalue = linregress(np.array(main_data[:, -1, 5]), pred[:, -1, -1]).rvalue\n",
    "        sse_value = sse(np.array(main_data[:, -1, 5]), pred[:, -1, -1])\n",
    "        plt.scatter(np.array(main_data[:, -1, 5]), pred[:, -1, -1], label=\"R = {:.3f}\\nSSE = {:.3f}\".format(rvalue, sse_value))\n",
    "        plt.legend()\n",
    "        plt.show()\n",
    "        objective_rval.append(rvalue)\n",
    "        objective_sse.append(sse_value)\n",
    "        \n",
    "        # randomly search over design space\n",
    "        t0 = time.time()\n",
    "        new_experiment_inds = brnn.fast_search(main_data_scaled, objective, scaler, n_test, P=main_obj_params)\n",
    "        new_experiments = unique_exp_names[new_experiment_inds]\n",
    "        elapsed_time.append(time.time()-t0)\n",
    "\n",
    "        # collect new data \n",
    "        new_df   = main_df.iloc[np.in1d(all_exp_names, new_experiments)].copy()\n",
    "        new_data, new_obj_params, new_experiments, N = format_data(new_df, species, outputs, controls, obj_params=obj_params)\n",
    "\n",
    "        # remove training samples from main dataset\n",
    "        train_inds = np.in1d(unique_exp_names, new_experiments)\n",
    "        main_data = main_data[~train_inds]\n",
    "        main_obj_params  = main_obj_params[~train_inds]\n",
    "        unique_exp_names = unique_exp_names[~train_inds]\n",
    "\n",
    "        # compute objectives\n",
    "        target_found = []\n",
    "        for sample, new_obj_param in zip(new_data, new_obj_params):\n",
    "            target_found.append(objective(sample[:, :len(sys_var)], new_obj_param))\n",
    "        target_found = np.array(target_found)\n",
    "\n",
    "        # store the best objective found (so far)\n",
    "        objective_found.append(np.max([np.max(objective_found), np.max(target_found)]))\n",
    "\n",
    "        # Update dataset\n",
    "        train_df = pd.concat((train_df, new_df))\n",
    "        train_data, train_obj_params, train_experiments, N = format_data(train_df, species, outputs, controls, obj_params=obj_params)\n",
    "        \n",
    "    ### fit model one last time to assess final prediction performance ### \n",
    "    # scale train and design space data\n",
    "    scaler = ZeroMaxScaler().fit(train_data)\n",
    "    train_data = scaler.transform(train_data)\n",
    "    main_data_scaled = scaler.transform(main_data)\n",
    "\n",
    "    # fit model \n",
    "    brnn = RNN(n_species=n_s, n_metabolites=1, n_controls=len(controls), n_hidden=16, N=N)\n",
    "    brnn.fit(train_data)\n",
    "\n",
    "    # assess prediction performance of end-point product\n",
    "    pred, stdv, cov = brnn.predict(main_data_scaled)\n",
    "    pred = scaler.inverse_transform(pred)\n",
    "    stdv = scaler.inverse_transform(stdv)\n",
    "    rvalue = linregress(np.array(main_data[:, -1, 5]), pred[:, -1, -1]).rvalue\n",
    "    sse_value = sse(np.array(main_data[:, -1, 5]), pred[:, -1, -1])\n",
    "    plt.scatter(np.array(main_data[:, -1, 5]), pred[:, -1, -1], label=\"R = {:.3f}\\nSSE = {:.3f}\".format(rvalue, sse_value))\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    objective_rval.append(rvalue)\n",
    "    objective_sse.append(sse_value)\n",
    "        \n",
    "    # save data to dataframe\n",
    "    dtl_df_i = pd.DataFrame()\n",
    "    dtl_df_i['Trial'] = [trial]\n",
    "    for j,obj_found in enumerate(objective_found):\n",
    "        dtl_df_i[f'DTL {j}'] = [obj_found]\n",
    "    dtl_df = pd.concat((dtl_df, dtl_df_i))\n",
    "    \n",
    "    # save data to dataframe\n",
    "    dtl_df_r = pd.DataFrame()\n",
    "    dtl_df_r['Trial'] = [trial]\n",
    "    for j,r_val in enumerate(objective_rval):\n",
    "        dtl_df_r[f'DTL {j}'] = [r_val]\n",
    "    dtl_df_R = pd.concat((dtl_df_R, dtl_df_r))\n",
    "    \n",
    "    # save data to dataframe\n",
    "    dtl_df_e = pd.DataFrame()\n",
    "    dtl_df_e['Trial'] = [trial]\n",
    "    for j,e in enumerate(objective_sse):\n",
    "        dtl_df_e[f'DTL {j}'] = [e]\n",
    "    dtl_df_sse = pd.concat((dtl_df_sse, dtl_df_e))\n",
    "    \n",
    "    dtl_df.to_csv(\"results/RNN_explore_exploit.csv\", index=False)\n",
    "    dtl_df_R.to_csv(\"results/RNN_explore_exploit_rvals.csv\", index=False)\n",
    "    dtl_df_sse.to_csv(\"results/RNN_explore_exploit_sse.csv\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtl_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtl_df_R.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtl_df_sse.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Average   search time {:.3f}\".format(np.mean(elapsed_time)))\n",
    "print(\"Std. dev. search time {:.3f}\".format(np.std(elapsed_time)))"
   ]
  }
 ],
 "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}