{ "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.14" }, "name": "", "signature": "sha256:a6848aab80b8f825790dc9339684a6c244ae2e9717fd31b806ca8760d13bfb38" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 1: Hierarchization in Higher Dimensions #\n", "\n", "In this exercise we will implement the multi-recursive algorithm for hierarchization of a multi-dimensional\n", "regular sparse grid. The structure of the code resembles strongly the one-dimensional case, and so we\n", "again have a class (PagodaFunction) representing our grid points.\n" ] }, { "cell_type": "code", "collapsed": true, "input": [ "#%load_ext sympyprinting\n", "%matplotlib inline\n", "from sympy import init_printing\n", "init_printing()\n", "import time\n", "# import the PagodaFunction class and additional helpers for plotting\n", "from PagodaFunction import PagodaFunction as gp\n", "from Plotting import plotHierarchical2d, evaluateSparseGrid, printSparseGrid, plotCombi2d" ], "language": "python", "metadata": { "scrolled": true }, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "def parabola(xvec):\n", " '''The parabola test function in multiple dimensions'''\n", " result = 1.0\n", " for x in xvec:\n", " result *= 4.0*x*(1.0-x)\n", " return result\n", " \n", "def applyRefinementCriterion(sg, dim, criterion, func):\n", " '''\n", " Applies the criterion to the leaves of the sparse grid.\n", " \n", " Calls the \"apply\" member for each leaf of the grid,\n", " then finishes calling the \"finalize\" member passing the\n", " grid and the approximated function as arguments.\n", " We'll also use it to build our sparse grid.\n", " '''\n", " # iterate over all grid points\n", " for (lKey,iKey) in sg:\n", " \n", " # Get the level vector and index vector of this grid point\n", " lvec, ivec = list(lKey), list(iKey)\n", " \n", " # Determine if this is a Leaf grid point by checking if it has exactly 2*d children.\n", " # Missing any child will make it a Leaf!\n", " isLeaf = False\n", " for d in range(dim):\n", " # original level,index values\n", " l_d = lvec[d]\n", " i_d = ivec[d]\n", " # Left child\n", " lvec[d] = l_d + 1\n", " ivec[d] = i_d*2 - 1\n", " leftChildKey = (tuple(lvec), tuple(ivec))\n", " # Right child\n", " ivec[d] = i_d*2 + 1\n", " rightChildKey = (tuple(lvec), tuple(ivec))\n", " # restore values for next dimension iteration\n", " lvec[d] = l_d \n", " ivec[d] = i_d\n", " \n", " if (leftChild not in sg or rightChild not in sg):\n", " isLeaf = True\n", " break # break d-loop only\n", " \n", " # Only if this is a Leaf, we can refine it\n", " if isLeaf:\n", " criterion.apply(sg[key])\n", " \n", " criterion.finalize(sg, func)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (i) ###\n", "Implement the refinement criterion MinLevelCriterion that adds all points up to a specified\n", "level to a given grid.
\n", "Hint: In your grid traversal, try to avoid multiple visits to the same grid points.\n" ] }, { "cell_type": "code", "collapsed": true, "input": [ "class MinLevelCriterion(object):\n", " '''Criterion adding all missing points up to a certain level.'''\n", " def __init__(self, dim, minLevel):\n", " self.dim = dim\n", " self.minLevel = minLevel # the overall level we want our sparse grid to have\n", " def apply(self, af):\n", " pass # work is done in finalize\n", " def finalize(self, sg, func):\n", " '''Inserts the grid points recursively'''\n", " def insertRecursively(minD, l, lvec, ivec, xvec):\n", " \n", " # minD: (1, 2, ..., minD, minD+1, ..., d) \n", " # All dimensions that is < minD is not allowed to be modified.\n", " # minD and minD+ dimensions are allowed to be modified. \n", " # For the purpose of visiting each grid point exactly once.\n", " \n", " # l: (integer) level of the grid point to be added\n", " # lvec: (list) level vector of the grid point to be added, [l1, l2, ..., ld]\n", " # ivec: (list) index vector of the grid point (list), [i1, i2, ..., id]\n", " # xvec: (list) coordinate vector of the grid point, [x1, x2, ..., xd]\n", " \n", " '''\n", " Local recursive helper function for grid traversal.\n", "\n", " Use lvec and ivec (mutable lists) to identify the current point.\n", " Insert it if it does not exist, descend recursively if not\n", " on level minLevel. minD helps to ensure to not pay any point \n", " multiple visits.\n", " '''\n", " \n", " # Construct Key of the grid point to be added\n", " key = (tuple(lvec), tuple(ivec))\n", " \n", " # If this point does not exist, insert it to sg\n", " if not key in sg: \n", " sg[key] = gp(self.dim, lvec, ivec, func(xvec))\n", "\n", " # If this grid point has not reached the desired overall level (minLevel), \n", " # we will try to add its children on each dimension\n", " if (l < self.minLevel):\n", " #TODO\n", " \n", " '''Issue top level call of the recursive function'''\n", " #TODO" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (ii) ###\n", "Implement the function hierarchize efficiently using a recursive approach.
\n", "Hint: The underlying traversal algorithm can be implemented similar to the one above.\n" ] }, { "cell_type": "code", "collapsed": true, "input": [ "def hierarchize(sg, dim):\n", " '''Recursive hierarchization function.'''\n", " def hierarchize1d(workDim, lvec, ivec, f_l, f_r):\n", " '''Recursive 1-D hierarchization to (lvec,ivec) in dimension workDim.'''\n", " key = (tuple(lvec), tuple(ivec))\n", " if not key in sg:\n", " return False\n", " p = sg[key]\n", " f_m = p.getSurplus() # surplus was first initialized with function value\n", " \n", " # original level, index values on workDim dimension\n", " l_w = lvec[workDim]\n", " i_w = ivec[workDim]\n", " # hierarhize left child\n", " lvec[workDim] = l_w+1 \n", " ivec[workDim] = 2*i_w - 1 \n", " hierarchize1d(workDim, lvec, ivec, f_l, f_m) \n", " # hierarhize right child\n", " ivec[workDim] = 2*i_w + 1 \n", " hierarchize1d(workDim, lvec, ivec, f_m, f_r)\n", " # restore values for next dim iteration\n", " lvec[workDim] = l_w\n", " ivec[workDim] = i_w\n", " p.setSurplus(f_m - 0.5*(f_l+f_r)) # Actual hierarchization step: hierarchize self\n", " return True\n", "\n", " def traverseRecursively(workDim, minD, lvec, ivec):\n", " # workDim: is the dimension on which we do hierarchization\n", " \n", " # minD: is the minimum dimension on which we access children. That is, we access children \n", " # on all dimensions ranging from minD to dim-1 \n", " \n", " # lvec, ivec: together defines the grid point on which we want to apply \n", " # the recursive hierarchization algorithm\n", " \n", " '''\n", " Local recursive helper function for traversal of the \"main axis\" w.r.t. workDim.\n", "\n", " Use the same approach as in MinLevelCriterion.finalize to traverse only the main\n", " axis of the grid, i.e. all points with coordinate x_workDim=0.5\n", " Hint: Don't descend into workDim, as hierarchize1d will do that.\n", " '''\n", " key = (tuple(lvec), tuple(ivec))\n", " # apply 1d hierarchization!\n", " if not hierarchize1d(workDim, lvec, ivec, 0.0, 0.0):\n", " return\n", " # being here => point exists\n", " if verbose:\n", " print(\" \" * (2 * (sum(lvec)-dim + 1 - 1)), \"dimension %d: hierarchize\" % workDim, key)\n", "\n", " #TODO implement\n", " pass \n", " \n", " # apply the 1-D hierarchization scheme for all dimensions, one after another\n", " lvec, ivec = [1,]*dim, [1,] * dim\n", " for d in range(dim):\n", " traverseRecursively(d, 0, lvec, ivec)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "'''Main'''\n", "# Funcion to interpolate using the sparse grid\n", "func = parabola\n", "\n", "# A dictionary(map) that contains all sparse grid points\n", "# Dictionary: {Key1:Value1, Key2:Value2, ...}\n", "# Key = ((l1,l2,...,ld),(i1, i2,...,id))\n", "# Value = PagodaFunction object\n", "sg = {}\n", "\n", "# Dimensionality of the problem\n", "dim = 2\n", "\n", "# Minimum Overall level of the Sparse grid\n", "minLevel = 4\n", "\n", "# Create the sparse grid using MinLevelCriterion\n", "applyRefinementCriterion(sg, dim, MinLevelCriterion(dim, minLevel), func)\n", "\n", "# plot\n", "print(\"Created sparse grid of level \", minLevel, \", with \", len(sg), \" grid points.\")\n", "printSparseGrid(sg)\n" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "tic = time.time()\n", "hierarchize(sg, dim)\n", "toc = time.time()\n", "print(\"Hierarchiz\",dim,\"dimensional sparse grid in\",toc-tic,\"seconds\")\n", "printSparseGrid(sg)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "tic = time.time()\n", "for key in sg:\n", " x = sg[key].computeCoordinate()\n", " print(key, x, evaluateSparseGrid(sg, x), func(x))\n", "toc = time.time()\n", "print(\"Evaluation of sparse grid in\", toc-tic, \"seconds\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "tic = time.time()\n", "plotHierarchical2d(sg, showGrid=True, minLevel=minLevel, azimuth=None, elevation=None)\n", "toc = time.time()\n", "print(\"Plotting the sparse grid took\", toc-tic, \"seconds\")" ], "language": "python", "metadata": { "scrolled": true }, "outputs": [], "prompt_number": null }, { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Optional: One-dimensional Sparse Grids---An Adaptiave Implementation #" ] }, { "cell_type": "code", "collapsed": true, "input": [ "# import classes\n", "from AnsatzFunction import AnsatzFunction as gp1d\n", "from Plotting import plotHierarchical1d, plotSG1d" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to use an adaptiave version of Archimedes' approach to approximate the integral $F(f,a,b) = \\int_a^b f(x) \\; dx$ of a function $f: \\mathbb{R} \\rightarrow \\mathbb{R}$ or to approximate the function $f$ itself.
\n", "For the one-dimensional case we want to formalize this approach and generalize it in the following ways:\n", "\n", "* Let $\\phi(x)$ be the mother of all hat functions with" ] }, { "cell_type": "code", "collapsed": true, "input": [ "gp1d(0,0).getPhi()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "* The data structure used to store the hierarchical coefficients is now called Sparse Grid.\n", "* A sparse grid is defined by a particular set of interpolation points $x_{l,i}$ and associated ansatz functions $\\phi_{l,i}(x)$ with\n", "$$\\phi_{l,i}(x) = \\phi ( \\; 2^l \\cdot \\left(x - i \\cdot \\frac{1}{2^l} \\right)\\; ) = \\phi (2^l \\cdot x - i ), \\ l \\in \\mathbb{N}^{+}, 1 \\leq i \\leq 2^l - 1, i \\mbox{ odd}$$\n", "* Archimedes' approach from the lecture corresponds to a regular sparse grid.\n", "* To improve the quality of approximation for arbitrary functions $f$ we introduce spatial adaptivity.\n", "\n", "Import and use the class AnsatzFunction and look at the comments in the provided code snippets for some more details.
\n", "Use it to construct a one-dimensional sparse grid structure, which stores its grid points (resp. ansatz functions of type AnsatzFunction in a hash map.\n" ] }, { "cell_type": "code", "collapsed": true, "input": [ "# the function we want to approximate\n", "def parabola(x):\n", " return 4*x*(1-x)\n", "\n", "def asymmetric(x):\n", " return 8*(-16*x**4 + 40*x**3 -35*x**2 + 11*x)/9\n", "\n", "func = parabola\n", "\n", "def printSparseGrid(sg):\n", " for af in sg.values():\n", " x = af.computeCoordinate()\n", " print(af.getKey(), \": f(%f)=%f\\tu(%f)=%f\" % (x, af.getFunctionValue(), x, af.getSurplus()))\n", "\n", "def evaluateSparseGrid(sg, x):\n", " '''Evaluates a given sparse grid at given point x'''\n", " return sum(map(lambda ansatz: ansatz(x), sg.values()))\n", "\n", "\n", "def applyRefinementCriterion(sg, criterion, func):\n", " '''\n", " Applies the criterion to the leaves of the sparse grid.\n", " \n", " Calls the \"apply\" member for each leaf of the grid,\n", " then finishes calling the \"finalize\" member passing the\n", " grid and the approximated function as arguments.\n", " We'll also use it to build our sparse grid.\n", " '''\n", " # iterate over leaves only\n", " for key in sg:\n", " l = key[0]\n", " i = key[1]\n", " \n", " leftChild = (l+1, 2*i-1)\n", " rightChild = (l+1, 2*i+1)\n", " \n", " if not (leftChild in sg and rightChild in sg):\n", " criterion.apply(sg[key])\n", " criterion.finalize(sg, func)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start by creating an empty sparse grid." ] }, { "cell_type": "code", "collapsed": true, "input": [ "sg = {}" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "We implement a refinement criterion that rather works as a grid generator and fills the sparse grid with points up to a given level.
\n", "The apply member does not really need to do anything. finalize then inserts all missing points up to level minLevel." ] }, { "cell_type": "code", "collapsed": true, "input": [ "class MinLevelCriterion(object):\n", " '''Criterion adding all missing points up to a certain level.'''\n", " def __init__(self, minLevel):\n", " self.minLevel = minLevel\n", " def apply(self, af):\n", " pass\n", " def finalize(self, sg, func):\n", " for l in range(1, self.minLevel+1):\n", " for i in range(1, 2**l, 2):\n", " if not (l,i) in sg:\n", " x = float(i)/2**l\n", " f_x = func(x)\n", " sg[(l,i)] = gp1d(l, i, f_x, f_x)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "# fill the sparse grid\n", "applyRefinementCriterion(sg, MinLevelCriterion(3), func)\n", "printSparseGrid(sg)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "We implement the hierarchization method.\n", "Try to only use the surpluses stored with the ansatz functions, not the function values.
\n", "Hint: In Python functions can be defined in every scope. You could make use of this feature and define and call to a recursive sub-function inside hierarchize." ] }, { "cell_type": "code", "collapsed": true, "input": [ "def hierarchize(sg):\n", " '''Recursive hierarchization function.'''\n", " # TODO\n", " pass" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "# apply hierarchization\n", "hierarchize(sg)\n", "printSparseGrid(sg)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a more complex refinement criterion looking for and refining surpluses which are larger than a given $\\varepsilon$.
\n", "Hint: This criterion can only be applied to a hierarchized grid. Remember this when ensuring a consistent state during the insertion of new grid points." ] }, { "cell_type": "code", "collapsed": true, "input": [ "class EpsilonCriterion(object):\n", " '''Criterion refining points with large surpluses.'''\n", " def __init__(self, eps):\n", " # TODO\n", " pass\n", " def apply(self, sg):\n", " # TODO\n", " pass\n", " def finalize(self, sg, func):\n", " # TODO \n", " pass" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "# apply the criterion\n", "applyRefinementCriterion(sg, EpsilonCriterion(0.06), func)\n", "printSparseGrid(sg)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, have a look at the sparse grid." ] }, { "cell_type": "code", "collapsed": true, "input": [ "plotHierarchical1d(sg, flat=True)\n", "print(\"Basis functions with surpluses on each level\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [ "plotSG1d(sg, func)\n", "print(\"Function interpolation with this sparse grid\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null }, { "cell_type": "code", "collapsed": true, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": null } ], "metadata": {} } ] }