{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1 (using Chaospy (3.0.5))\n",
    "\n",
    "$ x \\in [-2, 2] $\n",
    "\n",
    "$\n",
    "y = x - 0.5x ^ 3 + 0.1x ^ 4\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# true model\n",
    "def model(x):\n",
    "    return x - 0.5 * x ** 3 + 0.1 * x ** 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# set the domain of interest\n",
    "pdom = np . empty ((2 , 1) )\n",
    "pdom[0 , 0] = -2.0\n",
    "pdom[1 , 0] = 2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the UQ interface\n",
    "from surrogate.chaospy_model import UQChaospy as UQ\n",
    "\n",
    "# setup the UQ analysis ( proj : projection method )\n",
    "uq = UQ(pdom, order=4, method='proj')\n",
    "\n",
    "# create the training data\n",
    "xtrain = uq.generate_quadrature()\n",
    "ytrain = model(xtrain)\n",
    "\n",
    "# train the surrogate model\n",
    "uq.fit(xtrain, ytrain )\n",
    "\n",
    "# predict the model response on validation points\n",
    "xval = np.random.uniform(pdom [0 , 0], pdom[1 , 0], (5 , 1))\n",
    "yval = uq.predict(xval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(12, 6), dpi=150)\n",
    "\n",
    "# evaluate true model\n",
    "x = np.linspace(pdom[0, 0], pdom[1, 0], 100)\n",
    "plt.plot(x, model(x), label='true model')\n",
    "\n",
    "# plot training points\n",
    "plt.plot(xtrain, ytrain, 'ro', label='training points')\n",
    "\n",
    "# plot validation points\n",
    "plt.plot(xval, yval, 'bo', label='predicted points')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f(x)')\n",
    "plt.grid(linestyle='dashed')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Ishigami function (using Chaospy (3.0.5))\n",
    "\n",
    "$\n",
    "m(x_1, x_2, x_3) = \\sin(x_1) + a \\cdot\\sin^2(x_2) + b \\cdot x_3^4\\sin(x_1)\n",
    "$\n",
    "\n",
    "with $x_i\\in [-\\pi, \\pi]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ishigami(x1, x2, x3, a, b):\n",
    "    import numpy as np\n",
    "    return np.sin(x1) + a * np.sin(x2) ** 2 + b * x3 ** 4 * np.sin(x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "l = -np.pi\n",
    "u =  np.pi\n",
    "\n",
    "order = 9\n",
    "\n",
    "a = 7\n",
    "b = 0.05\n",
    "\n",
    "nsamples = 440\n",
    "\n",
    "np.random.seed(42)\n",
    "x1 = np.random.uniform(l, u, nsamples)\n",
    "\n",
    "np.random.seed(47)\n",
    "x2 = np.random.uniform(l, u, nsamples)\n",
    "\n",
    "np.random.seed(53)\n",
    "x3 = np.random.uniform(l, u, nsamples)\n",
    "\n",
    "xtrain = np.vstack((x1, x2, x3)).T\n",
    "\n",
    "y_true = ishigami(xtrain[:, 0], xtrain[:, 1], xtrain[:, 2], a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analytical solution\n",
    "D1  = b * np.pi ** 4 / 5.0 + b ** 2 * np.pi ** 8 / 50.0 + 0.5\n",
    "D2  = a ** 2 / 8.0\n",
    "D13 = (1.0 / 18 - 1.0 / 50.0) * b ** 2 * np.pi ** 8\n",
    "D   = D1 + D2 + D13\n",
    "\n",
    "print ( D1 / D )\n",
    "print ( D2 / D )\n",
    "print ( 0 )\n",
    "print ( (D1 + 0 + D13) / D )\n",
    "print ( (D2 + 0 + 0) / D )\n",
    "print ( (0  + 0 + D13) / D )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from chaospy_model import UQChaospy\n",
    "import numpy as np\n",
    "\n",
    "pdom = np.empty((2, 3))\n",
    "\n",
    "for i in range(3):\n",
    "    pdom[0, i] = l\n",
    "    pdom[1, i] = u\n",
    "    \n",
    "ch = UQChaospy(pdom, order)\n",
    "\n",
    "ch.fit(xtrain, y_true)\n",
    "\n",
    "ytrain = ch.predict(xtrain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sens_m = ch.main_sensitivity()\n",
    "sens_t = ch.total_sensitivity()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ( sens_m )\n",
    "print ( sens_t )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bootstrap.bootstrap import bootstrap_sobol\n",
    "\n",
    "s_m, s_t = bootstrap_sobol(xtrain, y_true, ch, n_boot=20, seed=42)\n",
    "\n",
    "print ( s_m )\n",
    "print ( s_t )"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}