main_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variational inference in an SDE tumor growth model using a recurrent neural network\n",
    "\n",
    "The system of SDEs is given by:\n",
    "\\begin{align}\n",
    "\\log Y_{ij} &=  \\log V_{ij} + \\epsilon_{ij}, \\hspace{20pt} j = 1,2,\\dots,m_i; \\hspace{5pt} i=1, 2, \\dots, n_{subjects}\\\\\n",
    "V_{ij} &= V^s_{ij} + V^k_{ij} \\\\\n",
    "dV^s_{i,t} &= (\\beta_i - \\alpha_{k_i})V^s_{i,t} dt + \\gamma V^s_{i,t} dB^1_{i,t} \\\\\n",
    "dV^k_{i,t} &= (\\alpha_{k_i} V^s_{i,t} - \\delta_i  V^k_{i,t})dt + \\tau V^k_{i,t} dB^2_{i,t}\n",
    "\\end{align}\n",
    "where $m_i$ is the number of observations for subject $i$, $k_i \\in \\{1,2,\\dots, n_{groups}\\}$ maps the $ith$ patient to their treatment group index.\n",
    "\n",
    "The noise is assumed to be distributed according to:\n",
    "- $\\log \\epsilon_{ij} \\sim \\mathcal{N}(0, \\sigma^2),$\n",
    "\n",
    "which may equivalently be written as\n",
    "- $\\tilde{\\epsilon}_{ij}^{-1} \\sim \\tilde{\\epsilon}_{ij} \\sim \\mathcal{LN}\\left(e^{\\sigma^2/2}, e^{\\sigma^2}(e^{\\sigma^2}-1)\\right)$.\n",
    "\n",
    "\n",
    "The assumed prior distributions:\n",
    "- $p(\\alpha_k) = \\mathcal{LN}(\\alpha_k |\\mu_{\\alpha_k}, \\sigma^2_{\\alpha_k}) \\hspace{20pt} k=1\\dots n_{groups}$\n",
    "\n",
    "where $\\mu_{\\alpha_k}$ and $\\sigma_{\\alpha_k}$ are treatment group dependent hyperparameters, also\n",
    "- $p(\\beta_i) = \\mathcal{LN}(\\beta_i | 1, 3), \\hspace{20pt} i=1\\dots n_{subjects},$\n",
    "- $p(\\delta_i) = \\mathcal{LN}(\\delta_i | 1, 3),$\n",
    "- $p(\\gamma) = \\mathcal{LN}(\\gamma | 1, 3),$\n",
    "- $p(\\tau) = \\mathcal{LN}(\\tau | 1, 3).$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import data\n",
    "import features\n",
    "from importlib import reload\n",
    "import distributions\n",
    "\n",
    "tfd = tf.contrib.distributions\n",
    "tfb = tfd.bijectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = reload(data)\n",
    "filename = 'dataset1.csv'\n",
    "n_groups = 3\n",
    "dt = 0.01\n",
    "P = 50\n",
    "obs_times, obs, y__0  = data.load_data(filename, n_groups)\n",
    "n = obs.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = reload(features)\n",
    "next_times, next_obs, prev_obs = features.extract_features(obs, obs_times, dt)\n",
    "next_times = tf.tile(tf.expand_dims(tf.expand_dims(next_times, axis = 0), axis = 0), [P, n, 1])\n",
    "next_obs = tf.tile(tf.expand_dims(next_obs, axis = 0), [P, 1, 1])\n",
    "prev_obs = tf.tile(tf.expand_dims(prev_obs, axis = 0), [P, 1, 1])\n",
    "\n",
    "features_dict = {'next_times': next_times, 'next_obs': next_obs, 'first_obs': tf.expand_dims(prev_obs[...,0], axis=2)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To - do: set x0 prior\n",
    "# Place priors on the parameters\n",
    "priors = {'gamma_mean': 1., 'gamma_std': 3., 'tau_mean': 1., 'tau_std': 3., 'alphas_mean': [0.1, 1., 3.],\n",
    "            'alphas_std': [0.5, 2., 3.], 'betas_mean': 1., 'betas_std': 3., 'deltas_mean': 1., 'deltas_std': 3.}\n",
    "# Network settings\n",
    "network_params = {'n_inputs': 6, 'n_outputs': 5 , 'n_hidden_layers': 4, 'hidden_layer_width': 100}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise variational distributions with the same parameters as the priors\n",
    "gamma_vdistr, gamma_vmean, gamma_vstd = distributions.init_param(priors['gamma_mean'], priors['gamma_std'])\n",
    "tau_vdistr, tau_vmean, tau_vstd = distributions.init_param(priors['tau_mean'], priors['tau_std'])\n",
    "alphas_vdistr, alphas_vmean, alphas_vstd = distributions.init_param(priors['alphas_mean'], priors['alphas_std'])\n",
    "betas_vdistr, betas_vmean, betas_vstd = distributions.init_param(priors['betas_mean'], priors['betas_std'])\n",
    "deltas_vdistr, deltas_vmean, deltas_vstd = distributions.init_param(priors['deltas_mean'], priors['deltas_std'])\n",
    "x0_vdistr, x0_vmean, x0_vstd = distributions.init_param(priors['x0_mean'], priors['x0_std'])\n",
    "# x0 = tf.contrib.distributions.MultivariateNormalDiag(loc = obs[:,0], scale_diag = sigma * np.ones(n))\n",
    "# To - do: set x0 prior\n",
    "params = {#'x0' : x0,\n",
    "          'gamma_vdistr': gamma_vdistr, 'gamma_vmean': gamma_vmean, 'gamma_vstd': gamma_vstd,\n",
    "          'alphas_vdistr': alphas_vdistr, 'alphas_vmean': alphas_vmean, 'alphas_vstd': alphas_vstd,\n",
    "          'betas_vdistr': betas_vdistr, 'betas_vmean': betas_vmean, 'betas_vstd': betas_vstd,\n",
    "          'deltas_vdistr': deltas_vdistr, 'deltas_vmean': deltas_vmean, 'deltas_vstd': deltas_vstd}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Not functional yet\n",
    "# sde_model = Model(obs, features_dict, priors, sigma, dt, T, P, params, network_params)"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}