Examples/ptBayeslands.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running a Parallel Tempering Bayeslands example\n",
    "\n",
    "In this notebook, we will and run a Bayeslands model using the parallel tempering method.\n",
    "\n",
    "## Bayes\n",
    "\n",
    "Given some data, $D$, we make inference regarding unknown parameters, denoted by  $\\theta$, via the posterior distribution $p(\\theta|D)$,  \n",
    "\n",
    "given by Bayes rule, $p(\\theta|D)=\\frac{p(D|\\theta)p(\\theta)}{P(D)}$\n",
    "\n",
    "$p(D|\\theta)$ is the likelihood of the data given the parameters\n",
    "\n",
    "$p(\\theta)$ is the prior\n",
    "\n",
    "$p(D)$ is a normalizing constant and equal to $\\int p(D|\\theta)p(\\theta)d\\theta$. \n",
    "\n",
    "Bayeslands uses MCMC to obtain samples of $\\theta$ from the posterior distribution by proposing values from some known distribution $q(\\theta)$. The value of $\\theta$ is  then set equal to these proposed values with a probability which ensures the detailed balance condition is met, otherwise the chain remains in its current position. The transition kernel which moves the Markov chain from one point to another in the parameter space is the product of the proposal distribution $q$ and the acceptance probability, $\\alpha$. Under certain conditions, the draws from this transition kernel converge to draws from the stationary distribution, $p(\\theta| D)$.\n",
    "\n",
    "## Parallel Tempering\n",
    "\n",
    "This method requires running Badlands thousands of times to explore the full parameter space. If a model takes 1 second, this is okay, but if a model takes a few minutes, then this very challenging. We use the concept of ***parallel tempering*** to help overcome this problem. Parallel tempering takes its name from thermodynamics of physical systems where the temperature of a system determines the ability to change. Parallel tempering is also known as replica exchange or the Multi-Markov Chain method. It is suitable for multi-modal distributions  by running multiple MCMC chains at different *temperatures* in parallel.\n",
    "\n",
    "   Metropolis Hastings |   MH-Parallel Tempering\n",
    ":-------------------------:|:-------------------------:\n",
    "![](notebook_images/alg_MH.png)  |  ![](notebook_images/alg_MHPT.png)\n",
    "\n",
    "\n",
    "The parallel tempering algorithm can be visualised as follows:\n",
    "\n",
    "<img src=\"notebook_images/multi-core-ptbayeslands.png\" width=500>\n",
    "\n",
    "This comparison between the MCMC and Parallel Tempered exploration of rainfall in a Badlands model. Each sample is trying a slightly different value of rainfall (combined with other parametrs). In this case the ***true value*** of rainfall was 1.5 m/yr.\n",
    "\n",
    "  Rainfall exploration with MCMC |   Rainfall exploration Parallel Tempering\n",
    ":-------------------------:|:-------------------------:\n",
    "![](notebook_images/rain_mcmc.png)  |  ![](notebook_images/rain_pt.png)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set up the environemnt\n",
    "%matplotlib inline\n",
    "import sys\n",
    "sys.path.append('scripts/')\n",
    "\n",
    "#Import the script that produces a start and end topography\n",
    "from ptBayeslands import *\n",
    "\n",
    "# display plots in SVG format\n",
    "#%config InlineBackend.figure_format = 'svg' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "problemfolder = 'crater/'\n",
    "xmlinput = problemfolder + 'crater.xml'\n",
    "simtime = 100000\n",
    "resolu_factor = 1\n",
    "samples = 100\n",
    "\n",
    "likelihood_sediment = True\n",
    "\n",
    "#Set variables\n",
    "#m = 0.5\n",
    "m_min = 0.\n",
    "m_max = 2.\n",
    "\n",
    "#n = 1.\n",
    "n_min = 0.\n",
    "n_max = 2.\n",
    "\n",
    "#rain_real = 1.5\n",
    "rain_min = 0.\n",
    "rain_max = 3.\n",
    "\n",
    "#erod_real = 5.e-6\n",
    "erod_min = 3.e-6\n",
    "erod_max = 7.e-6\n",
    "\n",
    "minlimits_vec=[rain_min,erod_min,m_min,n_min]\n",
    "maxlimits_vec=[rain_max,erod_max,m_max,n_max]    \n",
    "\n",
    "\n",
    "#you can have different ratio values for different parameters depending on the problem. \n",
    "stepsize_ratio  = 0.02 \n",
    "stepratio_vec = [stepsize_ratio, stepsize_ratio, stepsize_ratio, stepsize_ratio] \n",
    "\n",
    "\n",
    "erodep_coords=np.array([[5,5],[10,10],[20,20],[30,30],[40,40],[50,50],[25,25],[37,30],[44,27],[46,10]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#update with additonal parameters. should have as many rows as parameters\n",
    "true_parameter_vec = np.loadtxt(problemfolder + 'data/true_values.txt')\n",
    "vec_parameters = np.random.uniform(minlimits_vec, maxlimits_vec) #  draw intial values for each of the free parameters\n",
    "num_param = vec_parameters.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = problemfolder + 'data/final_elev.txt'\n",
    "groundtruth_elev = np.loadtxt(datapath)\n",
    "groundtruth_erodep = np.loadtxt(problemfolder + 'data/final_erdp.txt')\n",
    "groundtruth_erodep_pts = np.loadtxt(problemfolder + 'data/final_erdp_pts.txt')\n",
    "\n",
    "fname = \"\"\n",
    "run_nb = 0\n",
    "while os.path.exists(problemfolder +'results_%s' % (run_nb)):\n",
    "    run_nb += 1\n",
    "if not os.path.exists(problemfolder +'results_%s' % (run_nb)):\n",
    "    os.makedirs(problemfolder +'results_%s' % (run_nb))\n",
    "    fname = (problemfolder +'results_%s' % (run_nb))\n",
    "    \n",
    "make_directory((fname + '/posterior/pos_parameters')) \n",
    "make_directory((fname + '/posterior/predicted_topo'))\n",
    "make_directory((fname + '/posterior/pos_likelihood'))\n",
    "make_directory((fname + '/posterior/accept_list'))\n",
    "if likelihood_sediment==True:\n",
    "    make_directory((fname + '/posterior/predicted_erodep'))\n",
    "\n",
    "make_directory((fname + '/pred_plots'))\n",
    "run_nb_str = 'results_' + str(run_nb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#-------------------------------------------------------------------------------------\n",
    "# Number of chains of MCMC required to be run\n",
    "# PT is a multicore implementation must num_chains >= 2\n",
    "# Choose a value less than the numbe of core available (avoid context swtiching)\n",
    "#-------------------------------------------------------------------------------------\n",
    "\n",
    "num_chains = 4 \n",
    "swap_ratio = 0.1 \n",
    "burn_in =0.1 \n",
    "\n",
    "#parameters for Parallel Tempering\n",
    "maxtemp = int(num_chains * 5)/2\n",
    "\n",
    "#how ofen you swap neighbours\n",
    "swap_interval = int(swap_ratio * (samples/num_chains)) \n",
    "\n",
    "#How often output is created\n",
    "num_successive_topo = 4\n",
    "sim_interval = np.arange(0,  simtime+1, simtime/num_successive_topo) \n",
    "\n",
    "print(\"Max temp: \", maxtemp)\n",
    "print(\"Samples: \", samples)\n",
    "print(\"swap ratio: \",swap_ratio)\n",
    "print(\"Number of chains: \", num_chains)\n",
    "print(\"Swap interval: \", swap_interval)\n",
    "print(\"Simulation time interval\", sim_interval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timer_start = time.time()\n",
    "\n",
    "#-------------------------------------------------------------------------------------\n",
    "#Create A a Patratellel Tempring object instance \n",
    "#-------------------------------------------------------------------------------------\n",
    "pt = ParallelTempering(  vec_parameters, num_chains, maxtemp, samples,swap_interval,fname, true_parameter_vec, num_param  ,  groundtruth_elev,  groundtruth_erodep_pts , erodep_coords, simtime, sim_interval, resolu_factor, run_nb_str, xmlinput,likelihood_sediment)\n",
    "\n",
    "#-------------------------------------------------------------------------------------\n",
    "# intialize the MCMC chains\n",
    "#-------------------------------------------------------------------------------------\n",
    "pt.initialize_chains( minlimits_vec, maxlimits_vec, stepratio_vec, likelihood_sediment,burn_in)\n",
    "\n",
    "#-------------------------------------------------------------------------------------\n",
    "#run the chains in a sequence in ascending order\n",
    "#-------------------------------------------------------------------------------------\n",
    "pos_param,likehood_rep, accept_list, combined_erodep, pred_elev = pt.run_chains()\n",
    "\n",
    "#print('sucessfully sampled') \n",
    "timer_end = time.time() \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While we wait, explore the function that we a are running in the [script](scripts/ptBayeslands.py) we call.\n",
    "As in the [Badlands Mountain Example](mountain.ipynb), you should be able to see where we manipulate the [XML file](crater/crater.xml) input variables. Look for the \"model = badlandsModel()\" line to see the familiar commands. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "likelihood = likehood_rep[:,0] # just plot proposed likelihood  \n",
    "likelihood = np.asarray(np.split(likelihood,  num_chains ))\n",
    "\n",
    "\n",
    "plt.plot(likelihood.T)\n",
    "plt.savefig( fname+'/likelihood.png')\n",
    "plt.show()\n",
    "plt.clf()\n",
    "plt.plot(accept_list.T)\n",
    "plt.savefig( fname+'/accept_list.png')\n",
    "plt.show()\n",
    "plt.clf()\n",
    "\n",
    "for i in range(sim_interval.size):\n",
    "        pos_ed  = combined_erodep[i, :, :] \n",
    "\n",
    "        #print(pos_ed) \n",
    "        erodep_mean = pos_ed.mean(axis=0)  \n",
    "        erodep_std = pos_ed.std(axis=0)  \n",
    "        plot_erodeposition(erodep_mean, erodep_std, groundtruth_erodep_pts[i,:], sim_interval[i], fname) \n",
    "        #np.savetxt(fname + '/posterior/predicted_erodep/com_erodep_'+str(sim_interval[i]) +'_.txt', pos_ed)\n",
    "\n",
    "pred_erodep = np.zeros(( groundtruth_erodep_pts.shape[0], groundtruth_erodep_pts.shape[1] )) # just to get the right size\n",
    "\n",
    "for i in range(sim_interval.size): \n",
    "        pos_ed  = combined_erodep[i, :, :] # get final one for comparision\n",
    "        pred_erodep[i,:] = pos_ed.mean(axis=0)   \n",
    "\n",
    "rmse, rmse_sed= mean_sqerror(  pred_erodep, pred_elev,  groundtruth_elev,  groundtruth_erodep_pts)\n",
    "\n",
    "\n",
    "print ('time taken  in minutes = ', (timer_end-timer_start)/60)\n",
    "print ('Folder: ', run_nb_str)\n",
    "np.savetxt(fname+'/time_sqerror.txt',[ (timer_end-timer_start)/60,  rmse_sed, rmse], fmt='%1.2f'  )\n",
    "\n",
    "#print(pos_param)\n",
    "\n",
    "mpl_fig = plt.figure()\n",
    "ax = mpl_fig.add_subplot(111)\n",
    "\n",
    "ax.boxplot(pos_param.T) \n",
    "ax.set_xlabel('Badlands parameters')\n",
    "ax.set_ylabel('Posterior') \n",
    "plt.legend(loc='upper right') \n",
    "plt.title(\"Boxplot of Posterior\")\n",
    "plt.savefig(fname+'/badlands_pos.png')\n",
    "plt.savefig(fname+'/badlands_pos.svg', format='svg', dpi=400)\n",
    "\n",
    "print(\"Number of chains, problem, folder, time, RMSE_sed, RMSE\")\n",
    "print (num_chains, problemfolder, run_nb_str, (timer_end-timer_start)/60, rmse_sed, rmse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Check the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.image as mpimg\n",
    "\n",
    "plt.figure(figsize=(12, 16))\n",
    "img=mpimg.imread('crater/results_0/pos_distri_0_pos_.png')\n",
    "imgplot = plt.imshow(img)\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(12, 16))\n",
    "img=mpimg.imread('crater/results_0/x_ymid_opt.png')\n",
    "imgplot = plt.imshow(img)\n",
    "plt.show()"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}