updated notebooks

chaospy_tut.ipynb

@@ -0,0 +1,300 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<!-- dom:TITLE: A tutorial on chaospy -->\n",
+    "# A tutorial on chaospy\n",
+    "<!-- dom:AUTHOR: Leif Rune Hellevik at Department of Structural Engineering, NTNU -->\n",
+    "<!-- Author: --> **Leif Rune Hellevik**, Department of Structural Engineering, NTNU\n",
+    "\n",
+    "Date: **Jun 21, 2017**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# ipython magic\n",
+    "%matplotlib notebook\n",
+    "%load_ext autoreload\n",
+    "%autoreload 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "\n",
+    "\n",
+    "# plot configuration\n",
+    "import matplotlib\n",
+    "import matplotlib.pyplot as plt\n",
+    "plt.style.use(\"ggplot\")\n",
+    "# import seaborn as sns # sets another style\n",
+    "matplotlib.rcParams['lines.linewidth'] = 3\n",
+    "\n",
+    "# font = {'family' : 'sans-serif',\n",
+    "#         'weight' : 'normal',\n",
+    "#         'size'   : 18.0}\n",
+    "# matplotlib.rc('font', **font)  # pass in the font dict as kwar"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# The chaospy module\n",
+    "\n",
+    "How to import the chaospy module"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "import chaospy as pc"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For convenience we use a very simple model (which may be replaced by\n",
+    "your deterministic model or PDE-solver), namely an exponential decay\n",
+    "function $y(t)$ with two parameters stored in the numpy array 'x':"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "\\begin{equation}\n",
+    "y(t) = x_0 \\, e^{-x_1 t}\n",
+    "\n",
+    "\\end{equation}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which may be implemented in python as:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "def model_solver(t, x):\n",
+    "    # Simple emulator of a PDE solver \n",
+    "    return x[0] * e**(-t*x[1])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and may be plotted by:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "t=linspace(0, 2, 200)\n",
+    "y = model_solver(t, [3,3]) \n",
+    "plot(t,y)\n",
+    "\n",
+    "# Create propability density functions (pdf) for model parameters"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### How to create distributions for model parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "pdf1 = pc.Uniform(0, 1)\n",
+    "pdf2 = pc.Uniform(0, 1)\n",
+    "jpdf = pc.J(pdf1, pdf2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Generate solutions of samples of model parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "nr_samples=300\n",
+    "X=jpdf.sample(nr_samples)\n",
+    "Y=array([model_solver(t, x) for x in X.T ]) #solve for a given time t=0.5\n",
+    "\n",
+    "mu=mean(Y, 0)\n",
+    "p05, p95 = percentile(Y, [5,95], 0)\n",
+    "fill_between(t, p05, p95, alpha=0.5)\n",
+    "plot(t, mu)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Generate statistics based on the sampled solutions  (Monte Carlo)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "nr_samples=400\n",
+    "X=jpdf.sample(nr_samples)\n",
+    "Y=array([model_solver(0.5, x) for x in X.T ]) #solve for a given time t=0.5\n",
+    "nr_samples_list=arange(nr_samples)+1\n",
+    "converge = cumsum(Y, 0)/nr_samples_list\n",
+    "plot(nr_samples_list, converge)\n",
+    "legstr=[]\n",
+    "legstr.append('random')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Compare sampling schemes in the parameter space"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Various sampling schemes\n",
+    "jpdf = pc.J(pc.Uniform(0,1), pc.Uniform(0,1))\n",
+    "ax1=subplot(121)\n",
+    "ax1.set_title('Random')\n",
+    "X1 = jpdf.sample(nr_samples)\n",
+    "scatter(*X1)\n",
+    "ax2=subplot(122)\n",
+    "X2 = jpdf.sample(nr_samples, \"S\")\n",
+    "ax2.set_title('Structured')\n",
+    "scatter(*X2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Impact of sampling strategy on convergence of Monte Carlo simulations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Effect of sampling on convergence of Monte Carlo simulations\n",
+    "X = jpdf.sample(nr_samples, \"S\")\n",
+    "Y = [model_solver(0.5, x) for x in X.T]\n",
+    "converge_structured = cumsum(Y, 0)/nr_samples_list\n",
+    "\n",
+    "# Compare convergence for random and structured sampling\n",
+    "plot(nr_samples_list, converge)\n",
+    "plot(nr_samples_list, converge_structured, \"c\")\n",
+    "legend(['random','structured'])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Polynomial chaos expansions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Polychaos expansions\n",
+    "poly = pc.orth_chol(1, jpdf)\n",
+    "X = jpdf.sample(10, \"S\")\n",
+    "Y = [model_solver(0.5, x) for x in X.T]\n",
+    "\n",
+    "approx = pc.fit_regression(poly, X, Y, rule=\"T\")\n",
+    "\n",
+    "nr_poly_samples=np.arange(20,500,50)\n",
+    "\n",
+    "order = 3\n",
+    "poly = pc.orth_chol(order+1, jpdf)\n",
+    "\n",
+    "mu_values=[]\n",
+    "\n",
+    "for psample in nr_poly_samples:\n",
+    "    X = jpdf.sample(psample, \"S\")\n",
+    "    Y = [model_solver(0.5, x) for x in X.T]\n",
+    "    approx = pc.fit_regression(poly, X, Y, rule=\"T\")\n",
+    "    mu_values.append(pc.E(approx,jpdf))\n",
+    "\n",
+    "plot(nr_poly_samples,mu_values)"
+   ]
+  }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 0
+}