{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Faster Computations with Numba" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Some notes mostly for myself, but could be useful to you" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Altough Python is fast compared to other high-level languages, it still is not as fast as C, C++ or Fortran. Luckily, two open source projects [Numba](http://numba.pydata.org/) and [Cython](http://cython.org/) can be used to speed-up computations. [Numba](http://numba.pydata.org/) is sponsored by the producer of [Anaconda](https://store.continuum.io/cshop/anaconda/), [Continuum Analytics](http://continuum.io/). Both projects allow you to convert your code to interpreted language so that it runs faster. Here I will use Numba, given its ease and Just-In-Time nature, although I still have to figure out how to save and use the compiled functions (newer versions of [Numba](http://numba.pydata.org/) seem to have introduced [Ahead-of-Time](https://numba.pydata.org/numba-doc/dev/user/pycc.html) compilation using pycc and static libraries it seems.). Cython on the other hand needs you to understand much more of C and Ctypes. But once you figure out a way to run your code and compile it, you have a library ready to use and import. See the code for the repository for my paper [Optimal consumption under uncertainty, liquidity constraints, and bounded rationality](http://ozak.github.io/BoundedConsumption/) where I used Cython for the dynsysf.pyx function." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Numba" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Numba is quite easy to use. Start by importing it\n", "\n", " import numba as nb\n", "\n", "or some of its functions\n", "\n", " from numba import jit, autojit, njit\n", " \n", "Then define your function and notate it using the Numba commands jit, njit, autojit or their decorators @jit, @njit." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "from numba import jit, njit, autojit, jitclass\n", "import numba as nb\n", "import math\n", "import warnings\n", "\n", "with warnings.catch_warnings():\n", " warnings.simplefilter('ignore', nb.errors.NumbaDeprecationWarning)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Examples:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.37 µs ± 231 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n", "6.19 µs ± 304 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "# A simple function that computes the maximum distance between two vectors\n", "@nb.jit('f8(f8[:],f8[:])')\n", "def errnb(V0,V1):\n", " maxi=0.0\n", " for i in range(V0.shape[0]):\n", " m = abs(V1[i]-V0[i])\n", " if m>=maxi:\n", " maxi = m\n", " return maxi\n", "\n", "x = np.random.random((1000))\n", "y = np.random.random((1000))\n", "\n", "%timeit errnb(x,y)\n", "%timeit np.max(np.abs(x-y))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "2.0" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Let's use Numba to compute a Cobb-Douglas function\n", "@jit\n", "def CobbDouglas(K, L, A, alpha, beta):\n", " return A * K**alpha * L**beta\n", "\n", "K, L, A, alpha, beta = 2, 2, 1, .3, .7\n", "\n", "CobbDouglas(K,L,A,alpha,beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's time it and compare with Numpy or simple Python" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "147 ns ± 2.49 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n" ] } ], "source": [ "%timeit A * K**alpha * L**beta" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "225 ns ± 12.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglas(K,L,A,alpha,beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Well not fast enough...why?" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CobbDouglas (int64, int64, int64, float64, float64)\n", "--------------------------------------------------------------------------------\n", "# File: \n", "# --- LINE 2 --- \n", "# label 0\n", "\n", "@jit\n", "\n", "# --- LINE 3 --- \n", "\n", "def CobbDouglas(K, L, A, alpha, beta):\n", "\n", " # --- LINE 4 --- \n", " # K = arg(0, name=K) :: int64\n", " # L = arg(1, name=L) :: int64\n", " # A = arg(2, name=A) :: int64\n", " # alpha = arg(3, name=alpha) :: float64\n", " # beta = arg(4, name=beta) :: float64\n", " # $0.4 = K ** alpha :: float64\n", " # del alpha\n", " # del K\n", " # $0.5 = A * $0.4 :: float64\n", " # del A\n", " # del $0.4\n", " # $0.8 = L ** beta :: float64\n", " # del beta\n", " # del L\n", " # $0.9 = $0.5 * $0.8 :: float64\n", " # del $0.8\n", " # del $0.5\n", " # $0.10 = cast(value=$0.9) :: float64\n", " # del $0.9\n", " # return $0.10\n", "\n", " return A * K**alpha * L**beta\n", "\n", "\n", "================================================================================\n" ] } ], "source": [ "CobbDouglas.inspect_types()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "OK...it is correctly compiled and it is using Numba types and not PyObjects. So perhaps we cannot gain much in this simple computation...but that is ok. Let's try something a little more complex, computing a CobbDouglas for vectors K and L" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Python function\n", "def CobbDouglas(K, L, A, alpha, beta):\n", " return A * K**alpha * L**beta\n", "\n", "CobbDouglas_nb = jit(CobbDouglas)\n", "CobbDouglas_nb2 = njit(CobbDouglas)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "K = np.random.random((100,1))\n", "L = np.random.random((100,1))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.52 µs ± 575 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglas(K,L,A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.7 µs ± 59.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglas_nb(K,L,A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.94 µs ± 621 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglas_nb2(K,L,A,alpha,beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Ooops what went wrong? Up to V.0.12.2 Numba could not create arrays in ``njit`` mode, i.e. without using Python Objects. But things seem to have improved a lot since then. But even for the old days, we do not need to despair, we can create a fast CobbDouglas function by iterating over our vectors. Here I first create the `CobbDouglasNB` function which takes floating point numbers and tells Numba not to use PyObjects, and then I create the vectorized version `CobbDouglasVecNB`, which iterates over the elements of K and L (here we assume both are vectors) and returns our output Y. I know it is strange that we have to give the output as part of the function, but it allows us to speed up the computations." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "@nb.njit('f8(f8, f8, f8, f8, f8)')\n", "def CobbDouglasNB(K, L, A, alpha, beta):\n", " return A * K**alpha * L**beta\n", "\n", "@nb.jit('f8[:](f8[:], f8[:], f8[:], f8, f8, f8)')\n", "def CobbDouglasVecNB(Y, K, L, A, alpha, beta):\n", " for i in range(K.shape[0]):\n", " Y[i]=CobbDouglasNB(K[i],L[i],A,alpha,beta)\n", " return Y" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "0.48966064683582006" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "K.mean()" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "9.17 µs ± 1.26 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglas(K[:,0],L[:,0],A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6.46 µs ± 157 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "Y=np.zeros_like(K)\n", "%timeit CobbDouglasVecNB(Y[:,0],K[:,0],L[:,0],A,alpha,beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "## Almost twice as fast as Numpy!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's generalize the function to arrays of 2 dimensions" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "@nb.jit('f8[:,:](f8[:,:], f8[:,:], f8[:,:], f8, f8, f8)')\n", "def CobbDouglasVecNB(Y, K, L, A, alpha, beta):\n", " for i in range(K.shape[0]):\n", " for j in range(K.shape[1]):\n", " Y[i,j]=CobbDouglasNB(K[i,j], L[i,j], A, alpha, beta)\n", " return Y" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "48.9 ms ± 3.69 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "K=np.random.random((1000,1000))\n", "L=np.random.random((1000,1000))\n", "%timeit CobbDouglas(K,L,A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "44 ms ± 573 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "Y=np.zeros_like(K)\n", "%timeit CobbDouglasVecNB(Y,K,L,A,alpha,beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "### 20% faster than Numpy...so we get the idea." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "While this is great to speed up, we see that there are some issues one has to think about in order to get those gains.\n", "\n", "Now let's use the vectorize option, which creates Numpy Ufunctions, that is functions that operate on vectors (or at leats that s what I gather)." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "@nb.vectorize('f8(f8, f8, f8, f8, f8,)')\n", "def CobbDouglasNB2(K, L, A, alpha, beta):\n", " return A * K**alpha * L**beta" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "45.4 ms ± 741 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "%timeit CobbDouglasNB2(K,L,A,alpha,beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Let's compare these functions on different sizes of matrices $K$ and $L$" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "K = np.random.random((10000, 10000))\n", "L = np.random.random((10000, 10000))\n", "Y = np.zeros_like(K)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "478 µs ± 18.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglas(K[1:100,1:100],L[1:100,1:100],A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "450 µs ± 9.08 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglasVecNB(Y[1:100,1:100],K[1:100,1:100],L[1:100,1:100],A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "428 µs ± 40.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglasNB2(K[1:100,1:100],L[1:100,1:100],A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "40.4 ms ± 2.11 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "%timeit CobbDouglas(K[1:1000,1:1000],L[1:1000,1:1000],A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "37.1 ms ± 2.33 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "%timeit CobbDouglasVecNB(Y[1:1000,1:1000],K[1:1000,1:1000],L[1:1000,1:1000],A,alpha,beta)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "39 ms ± 2.81 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "%timeit CobbDouglasNB2(K[1:1000,1:1000],L[1:1000,1:1000],A,alpha,beta)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Now let's try the CRRA utility function" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "@nb.jit('f8(f8, f8)')\n", "def U_numba(c, sigma):\n", " '''This function returns the value of utility when the CRRA\n", " coefficient is sigma. I.e. \n", " u(c,sigma)=(c**(1-sigma)-1)/(1-sigma) if sigma!=1 \n", " and \n", " u(c,sigma)=ln(c) if sigma==1\n", " Usage: u(c,sigma)\n", " '''\n", " if sigma!=1:\n", " u = (c**(1-sigma)-1)/(1-sigma)\n", " else:\n", " u = math.log(c)\n", " return u" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "@nb.jit('f8[:](f8[:], f8[:], f8)')\n", "def Uvec(u,c,sigma):\n", " if sigma!=1:\n", " for i in range(c.shape[0]):\n", " u[i] = U_numba(c[i], sigma)\n", " else:\n", " u = np.log(c)\n", " return u" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def Unp(c,sigma):\n", " if sigma!=1:\n", " u=(c**(1-sigma)-1)/(1-sigma)\n", " else:\n", " u=np.log(c)\n", " return u" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "@nb.vectorize('f8(f8,f8)')\n", "def Unb(c,sigma):\n", " if sigma!=1:\n", " u = (c**(1-sigma)-1)/(1-sigma)\n", " else:\n", " u = math.log(c)\n", " return u" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "@nb.jit('f8[:](f8[:], f8)')\n", "def U(c, sigma):\n", " if sigma!=1:\n", " u = Unb(c,sigma)#(c**(1-sigma)-1)/(1-sigma)\n", " else:\n", " u = np.log(c)\n", " return u" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# Grid of values for state variable over which function will be approximated\n", "gridmin, gridmax, gridsize = 0.1, 5, 300\n", "grid = np.linspace(gridmin, gridmax**1e-1, gridsize)**10" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.49 µs ± 113 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit U(grid,1)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.45 µs ± 63.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit Unp(grid,1)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.71 µs ± 152 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit Unb(grid,1)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.42 µs ± 181 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "u = np.zeros_like(grid)\n", "%timeit Uvec(u, grid, 1)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.93 µs ± 327 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit U(grid,3)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6.33 µs ± 58.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit Unb(grid,3)" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "8.68 µs ± 43.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit Unp(grid,3)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.73 µs ± 162 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "u = np.zeros_like(grid)\n", "%timeit Uvec(u,grid,3)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Optimizing the code for Dynamic Programming" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Optimal Growth" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Parameters Optimal Growth\n", "alpha = .3\n", "beta = .9\n", "sigma = 1\n", "delta = 1\n", "A = 1" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "tol=0.000001\n", "(300,)\n" ] } ], "source": [ "# Grid of values for state variable over which function will be approximated\n", "gridmin, gridmax, gridsize = 0.1, 5, 300\n", "grid = np.linspace(gridmin, gridmax**1e-1, gridsize)**10\n", "\n", "# Parameters for the optimization procedures\n", "count=0\n", "maxiter=1000\n", "tol=1e-6\n", "print('tol=%f' % tol)\n", "print(grid.shape)" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "import numba as nb\n", "from scipy.optimize import fminbound\n", "from scipy import interp\n", "\n", "# Auxiliary functions \n", "# Maximize function V on interval [a,b]\n", "def maximum(V, a, b, args=[]):\n", " return float(V(fminbound(lambda x: -V(x), a, b,args=args)))\n", "\n", "# Return Maximizer of function V on interval [a,b]\n", "def maximizer(V, a, b, args=[]):\n", " return float(fminbound(lambda x: -V(x), a, b, args=args))\n", "\n", "# Interpolation functions Class\n", "class LinInterp:\n", " \"Provides linear interpolation in one dimension.\"\n", "\n", " def __init__(self, X, Y):\n", " \"\"\"Parameters: X and Y are sequences or arrays\n", " containing the (x,y) interpolation points.\n", " \"\"\"\n", " self.X, self.Y = X, Y\n", "\n", " def __call__(self, z):\n", " \"\"\"Parameters: z is a number, sequence or array.\n", " This method makes an instance f of LinInterp callable,\n", " so f(z) returns the interpolation value(s) at z.\n", " \"\"\"\n", " if isinstance(z, int) or isinstance(z, float):\n", " return interp ([z], self.X, self.Y)[0]\n", " else:\n", " return interp(z, self.X, self.Y)" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "@nb.jit('f8[:](f8[:], f8)')\n", "def U(c,sigma):\n", " if sigma!=1:\n", " u = Unb(c,sigma)\n", " else:\n", " u = np.log(c)\n", " return u\n", "\n", "@nb.vectorize('f8(f8,f8)')\n", "def Unb(c,sigma):\n", " if sigma!=1:\n", " u = (c**(1-sigma)-1)/(1-sigma)\n", " else:\n", " u = math.log(c)\n", " return u" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "@nb.vectorize('f8(f8, f8, f8, f8, f8,)')\n", "def CobbDouglasNB2(K, L, A, alpha, beta):\n", " '''CobbDouglasNB2(K, L, A, alpha, beta)'''\n", " return A * K**alpha * L**beta\n", "\n", "@nb.vectorize('f8(f8, f8, f8)')\n", "def F_nb(K, alpha, A):\n", " '''\n", " Cobb-Douglas production function\n", " F(K)=A* K^alpha\n", " '''\n", " return A * K**alpha\n", "\n", "def F_np(K, alpha, A):\n", " '''\n", " Cobb-Douglas production function\n", " F(K)=A* K^alpha\n", " '''\n", " return A * K**alpha" ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "9.68 µs ± 382 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit CobbDouglasNB2(grid,1,A,alpha,0)" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.28 µs ± 357 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit F_nb(grid,alpha,A)" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.6 µs ± 538 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit F_np(grid,alpha,A)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "V0 = LinInterp(grid,U(grid,sigma))" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def bellman(x,w):\n", " \"\"\"The approximate Bellman operator.\n", " Parameters: w is a LinInterp object (i.e., a \n", " callable object which acts pointwise on arrays).\n", " Returns: An instance of LinInterp that represents the optimal operator.\n", " w is a function defined on the state space.\n", " \"\"\"\n", " vals = []\n", " for k in x:\n", " kmax=F_nb(k,alpha,A)\n", " h = lambda kp: U_numba(kmax + (1-delta) * k - kp,sigma) + beta * w(kp)\n", " vals.append(maximum(h, 0, kmax))\n", " return LinInterp(grid, vals)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def policy(x,w):\n", " \"\"\"\n", " For each function w, policy(w) returns the function that maximizes the \n", " RHS of the Bellman operator.\n", " Replace w for the Value function to get optimal policy.\n", " The approximate optimal policy operator w-greedy (See Stachurski (2009)). \n", " Parameters: w is a LinInterp object (i.e., a \n", " callable object which acts pointwise on arrays).\n", " Returns: An instance of LinInterp that captures the optimal policy.\n", " \"\"\"\n", " vals = []\n", " for k in x:\n", " kmax=F_nb(k,alpha,A)\n", " h = lambda kp: U_numba(kmax + (1-delta) * k - kp, sigma) + beta * w(kp)\n", " vals.append(maximizer(h, 0, kmax))\n", " return LinInterp(grid, vals)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def solve():\n", " count=0\n", " V0=LinInterp(grid,U(grid,sigma))\n", " while count:4: NumbaWarning: \n", "Compilation is falling back to object mode WITH looplifting enabled because Function \"bellmannb\" failed type inference due to: Untyped global name 'LinInterp': cannot determine Numba type of \n", "\n", "File \"\", line 12:\n", "def bellmannb(x,w0):\n", " \n", " \"\"\"\n", " w = LinInterp(x,w0)\n", " ^\n", "\n", " @nb.jit('f8[:](f8[:], f8[:])')\n", ":4: NumbaWarning: \n", "Compilation is falling back to object mode WITHOUT looplifting enabled because Function \"bellmannb\" failed type inference due to: Untyped global name 'LinInterp': cannot determine Numba type of \n", "\n", "File \"\", line 12:\n", "def bellmannb(x,w0):\n", " \n", " \"\"\"\n", " w = LinInterp(x,w0)\n", " ^\n", "\n", " @nb.jit('f8[:](f8[:], f8[:])')\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:742: NumbaWarning: Function \"bellmannb\" was compiled in object mode without forceobj=True, but has lifted loops.\n", "\n", "File \"\", line 5:\n", "@nb.jit('f8[:](f8[:], f8[:])')\n", "def bellmannb(x,w0):\n", "^\n", "\n", " self.func_ir.loc))\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:751: NumbaDeprecationWarning: \n", "Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.\n", "\n", "For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit\n", "\n", "File \"\", line 5:\n", "@nb.jit('f8[:](f8[:], f8[:])')\n", "def bellmannb(x,w0):\n", "^\n", "\n", " warnings.warn(errors.NumbaDeprecationWarning(msg, self.func_ir.loc))\n", ":19: NumbaWarning: \n", "Compilation is falling back to object mode WITH looplifting enabled because Function \"policynb\" failed type inference due to: Untyped global name 'LinInterp': cannot determine Numba type of \n", "\n", "File \"\", line 27:\n", "def policynb(x,w0):\n", " \n", " \"\"\"\n", " w = LinInterp(x,w0)\n", " ^\n", "\n", " @nb.jit('f8[:](f8[:],f8[:])')\n", ":19: NumbaWarning: \n", "Compilation is falling back to object mode WITHOUT looplifting enabled because Function \"policynb\" failed type inference due to: Untyped global name 'LinInterp': cannot determine Numba type of \n", "\n", "File \"\", line 27:\n", "def policynb(x,w0):\n", " \n", " \"\"\"\n", " w = LinInterp(x,w0)\n", " ^\n", "\n", " @nb.jit('f8[:](f8[:],f8[:])')\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:742: NumbaWarning: Function \"policynb\" was compiled in object mode without forceobj=True, but has lifted loops.\n", "\n", "File \"\", line 20:\n", "@nb.jit('f8[:](f8[:],f8[:])')\n", "def policynb(x,w0):\n", "^\n", "\n", " self.func_ir.loc))\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:751: NumbaDeprecationWarning: \n", "Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.\n", "\n", "For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit\n", "\n", "File \"\", line 20:\n", "@nb.jit('f8[:](f8[:],f8[:])')\n", "def policynb(x,w0):\n", "^\n", "\n", " warnings.warn(errors.NumbaDeprecationWarning(msg, self.func_ir.loc))\n", ":43: NumbaWarning: \n", "Compilation is falling back to object mode WITH looplifting enabled because Function \"solvenb\" failed type inference due to: Invalid use of type(CPUDispatcher()) with parameters (readonly array(float64, 1d, C), Literal[int](1))\n", "Known signatures:\n", " * (array(float64, 1d, A), float64) -> array(float64, 1d, A)\n", " * parameterized\n", "[1] During: resolving callee type: type(CPUDispatcher())\n", "[2] During: typing of call at (46)\n", "\n", "\n", "File \"\", line 46:\n", "def solvenb():\n", " \n", " count=0.0\n", " V0=U(grid,sigma)\n", " ^\n", "\n", " @nb.jit('f8[:]()')\n", ":43: NumbaWarning: \n", "Compilation is falling back to object mode WITHOUT looplifting enabled because Function \"solvenb\" failed type inference due to: cannot determine Numba type of \n", "\n", "File \"\", line 46:\n", "def solvenb():\n", " \n", " count=0.0\n", " V0=U(grid,sigma)\n", " ^\n", "\n", " @nb.jit('f8[:]()')\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:742: NumbaWarning: Function \"solvenb\" was compiled in object mode without forceobj=True, but has lifted loops.\n", "\n", "File \"\", line 44:\n", "@nb.jit('f8[:]()')\n", "def solvenb():\n", "^\n", "\n", " self.func_ir.loc))\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:751: NumbaDeprecationWarning: \n", "Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.\n", "\n", "For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit\n", "\n", "File \"\", line 44:\n", "@nb.jit('f8[:]()')\n", "def solvenb():\n", "^\n", "\n", " warnings.warn(errors.NumbaDeprecationWarning(msg, self.func_ir.loc))\n" ] } ], "source": [ "def Utrf(kp,kmax,k,sigma,w):\n", " return -U_numba(kmax + (1-delta) * k - kp,sigma)-beta*w(kp)\n", "\n", "@nb.jit('f8[:](f8[:], f8[:])')\n", "def bellmannb(x,w0):\n", " \"\"\"The approximate Bellman operator.\n", " Parameters: w is a LinInterp object (i.e., a \n", " callable object which acts pointwise on arrays).\n", " Returns: An instance of LinInterp that represents the optimal operator.\n", " w is a function defined on the state space.\n", " \"\"\"\n", " w = LinInterp(x,w0)\n", " vals = np.array([])\n", " for k in x:\n", " kmax = F_nb(k,alpha,A)\n", " vals = np.append(vals, -Utrf(fminbound(Utrf, 0, kmax, args=[kmax,k,sigma,w]),kmax,k,sigma,w))\n", " return vals\n", "\n", "@nb.jit('f8[:](f8[:],f8[:])')\n", "def policynb(x,w0):\n", " \"\"\"The approximate Bellman operator.\n", " Parameters: w is a LinInterp object (i.e., a \n", " callable object which acts pointwise on arrays).\n", " Returns: An instance of LinInterp that represents the optimal operator.\n", " w is a function defined on the state space.\n", " \"\"\"\n", " w = LinInterp(x,w0)\n", " vals = np.array([])\n", " for k in x:\n", " kmax = F_nb(k,alpha,A)\n", " vals = np.append(vals,fminbound(Utrf, 0, kmax, args=[kmax,k,sigma,w]))\n", " return vals\n", "\n", "@nb.jit('f8(f8[:],f8[:])')\n", "def errnb(V0,V1):\n", " maxi=0.0\n", " for i in range(V0.shape[0]):\n", " m=abs(V1[i]-V0[i])\n", " if m>=maxi:\n", " maxi=m\n", " return maxi\n", "\n", "@nb.jit('f8[:]()')\n", "def solvenb():\n", " count=0.0\n", " V0=U(grid,sigma)\n", " while count:4: NumbaWarning: \n", "Compilation is falling back to object mode WITHOUT looplifting enabled because Function \"bellmannb\" failed type inference due to: Untyped global name 'Utrf': cannot determine Numba type of \n", "\n", "File \"\", line 16:\n", "def bellmannb(x,w0):\n", " \n", " kmax = F_nb(k,alpha,A)\n", " vals = np.append(vals, -Utrf(fminbound(Utrf, 0, kmax, args=[kmax,k,sigma,w]),kmax,k,sigma,w))\n", " ^\n", "\n", " @nb.jit('f8[:](f8[:], f8[:])')\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:742: NumbaWarning: Function \"bellmannb\" was compiled in object mode without forceobj=True.\n", "\n", "File \"\", line 14:\n", "def bellmannb(x,w0):\n", " \n", " vals = np.array([])\n", " for k in x:\n", " ^\n", "\n", " self.func_ir.loc))\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:751: NumbaDeprecationWarning: \n", "Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.\n", "\n", "For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit\n", "\n", "File \"\", line 14:\n", "def bellmannb(x,w0):\n", " \n", " vals = np.array([])\n", " for k in x:\n", " ^\n", "\n", " warnings.warn(errors.NumbaDeprecationWarning(msg, self.func_ir.loc))\n", ":19: NumbaWarning: \n", "Compilation is falling back to object mode WITHOUT looplifting enabled because Function \"policynb\" failed type inference due to: Untyped global name 'fminbound': cannot determine Numba type of \n", "\n", "File \"\", line 31:\n", "def policynb(x,w0):\n", " \n", " kmax = F_nb(k,alpha,A)\n", " vals = np.append(vals,fminbound(Utrf, 0, kmax, args=[kmax,k,sigma,w]))\n", " ^\n", "\n", " @nb.jit('f8[:](f8[:],f8[:])')\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:742: NumbaWarning: Function \"policynb\" was compiled in object mode without forceobj=True.\n", "\n", "File \"\", line 29:\n", "def policynb(x,w0):\n", " \n", " vals = np.array([])\n", " for k in x:\n", " ^\n", "\n", " self.func_ir.loc))\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:751: NumbaDeprecationWarning: \n", "Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.\n", "\n", "For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit\n", "\n", "File \"\", line 29:\n", "def policynb(x,w0):\n", " \n", " vals = np.array([])\n", " for k in x:\n", " ^\n", "\n", " warnings.warn(errors.NumbaDeprecationWarning(msg, self.func_ir.loc))\n" ] } ], "source": [ "U0 = U(grid, sigma)\n", "\n", "U0 = bellmannb(grid, U0)\n", "\n", "c0 = policynb(grid, U0)" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "130 ms ± 2.09 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n", "127 ms ± 1.45 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":43: NumbaWarning: \n", "Compilation is falling back to object mode WITHOUT looplifting enabled because Function \"solvenb\" failed type inference due to: Invalid use of type(CPUDispatcher()) with parameters (readonly array(float64, 1d, C), array(float64, 1d, C))\n", " * parameterized\n", "[1] During: resolving callee type: type(CPUDispatcher())\n", "[2] During: typing of call at (48)\n", "\n", "\n", "File \"\", line 48:\n", "def solvenb():\n", " \n", " while count\", line 46:\n", "def solvenb():\n", " \n", " count=0.0\n", " V0=U(grid,sigma)\n", " ^\n", "\n", " self.func_ir.loc))\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:751: NumbaDeprecationWarning: \n", "Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.\n", "\n", "For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit\n", "\n", "File \"\", line 46:\n", "def solvenb():\n", " \n", " count=0.0\n", " V0=U(grid,sigma)\n", " ^\n", "\n", " warnings.warn(errors.NumbaDeprecationWarning(msg, self.func_ir.loc))\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "138.0\n", "138.0\n", "138.0\n", "138.0\n", "138.0\n", "138.0\n", "138.0\n", "138.0\n", "16.6 s ± 688 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit bellmannb(grid, U0)\n", "%timeit policynb(grid, U0)\n", "%timeit solvenb()" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAat0lEQVR4nO3de5Bc5Xnn8e/T9+65ayQhoZEQF2FzibHRICC+YQLGcWyDWRyjreyytne1cdkue12uZDHYa2/hcnCStZNK/rAWs1suZyM76xDjJQEjxyR2KgLEzSBksAAJjQah0WikufW9n/1jRmIGtRhJPWd65p3fp6pruvuced/niNJPL+95zznm7oiISJhizS5ARESio5AXEQmYQl5EJGAKeRGRgCnkRUQClmh2AVMtXbrU165d2+wyREQWlMcee+yguy+rty3ykDez9wF/BsSBu9z9j06079q1a9m+fXvUJYmIBMXM9pxoW6TTNWYWB/4S+G3gQmCjmV0YZZ8iIvKaqOfkNwC73P1Fdy8BW4DrI+5TREQmRR3yq4C9Uz73TX4nIiJzIOo5eavz3bT7KJjZJmATwJo1a47buVwu09fXR6FQiKTAuZDJZOjp6SGZTDa7FBFZZKIO+T5g9ZTPPUD/1B3cfTOwGaC3t/e4G+n09fXR1tbG2rVrMav3b8b85u4MDg7S19fH2Wef3exyRGSRiXq65lFgnZmdbWYp4Gbg3lNpoFAo0N3dvSADHsDM6O7uXtD/JyIiC1ekI3l3r5jZp4EHmFhCebe77zjVdhZqwB+10OsXkYUr8nXy7v73wN9H3Y+IyFzzWo1yuUSpmKdcLFAuFSgXx6lMvq+UC1RLRWrlAtVSgWq5iFcK1MpFvHL0VYBKidTKC1j//k/Meo3z6orX+eiqq67i1ltv5brrrjv23be+9S2ef/55Lr/8cu644w4Abr/9dm655ZZmlSmy6NSqVUrFPMVCnnJhnGJhnEpxnHIpT6WYp1Icp1rOUy0VqJUK1Mp5vFzAKwW8XIBKAasUsWoJqxaxWolYtUSsViJeKxOvlYjXSiS8RMLLx15JyqS8RIoKaSuTAlKzcDyPHb4aFPJzb+PGjWzZsmVayG/ZsoU777yTj33sY2zfvh0zY/369XzoQx+iq6uridWKzD2v1SgWxinmxyjkRynlRykV8pQLoxNBWxqfDNo8tXLhBEFbxCoFYtUi8VqRWLVIolYkcSxkS6RqxYmApUTay6StTAbInGbdVTeKpChbgjJJSpaiYkkqJKnEUlQtSSWephhrpxpLUYul8HiKWjyFx1J4Io3F03gihSXSEE9jiTSxZIZY8rWf8WSaRDJDPJ0lkcqQSGVIHv2ZzpLOZEmlMqyPx2fzP8sxCvkZ3HTTTdx+++0Ui0XS6TS7d++mv7+f/v5+rr32WpYsWQLAtddey/3338/GjRubXLHIhEq5RCE/RmF8lFJ+jHJhjFJh4melODHqrRXHqJXz1ErjeGl8InTLeaySxyoF4tU88WqBeLVAslYkWSuQ9CKpWpE0RdJeIkOJjDkZoOMUazwatEVLUyJJ2VITr1iaiqUoxXPk411UY2lq8YmXx9N4IguJNCTSWDKDJbOToTrxSqRzxFNZEuksyXSWRDpHMp0llc6SzraQzuRIJFPkoviDn2cWVMh/9cc7eLZ/eFbbvPDMdv7bBy864fbu7m42bNjA/fffz/XXX8+WLVv46Ec/yr59+1i9+rXVoT09Pezbt29Wa5Pw1apV8uMj5MdGKI4PUxwfpZQfoZwfo1IYoVoYpVocw0ujEyFcGidWHiVWyROvjBOvjE8EbzVP2vOkawUyFMh6gZRVaQVaT7GmvKcoWJoSKUqWphTLULY0lXiaQrJjInATWWqJLB5PQzKHJ7MTQZvKEUvliKezxFMtxNM5kpncjEG7GMK2WRZUyDfL0SmboyF/991389Of/vS4/bSKJnyVcomxkSOMjxyiMDZMcewIpbEjVPLDVAoj1AojeGEYSqPESqPEy6MTQVwdJ1nNk6oVSHuejBfIeJGcFWkBWk6y/5InyFuaPFmKsQwly1CKZxlPLWEknqWSaMETWWrJFkhmsVQOm/wZT+eIp3Ik0jkSmRaS6RypbAupbCvpTAuZXCvpTI5sLEY2yj9EmVMLKuTfaMQdpRtuuIHPf/7zPP744+TzeS699FKee+45HnrooWP79PX1cdVVVzWlPplZqVhg5PBBxkeGyA8fojh6iNLoYar5w9QKw3hxFCuOECuPEi+PkaiMkayMka6Nk66Nk/U8LT5Oxsp0MPO0RMkTjFmOvGUpWJZSPEcp3sJ4ehnVeJZqMocncniqBUvmsHQrsXQL8XQLiUwbiWzrRPjm2kjn2knn2si1tJFKpUmdRP8iRy2okG+W1tZWrrrqKj7+8Y8fm3O/7rrr+OIXv8jQ0BAAP/nJT/j617/ezDKDVhgfnQjp4UPkRw5RHB2iPDZEdfwItfxhKAwTKx0hURohWR4hXRkhUxslVxuj1cfIWoluoPsN+hj3NOOToVyI5SjGc4ykljOUaKGaasWTrXiqFcu0EUu3Ec+2kci2k851kGrpINPSQa6tk1xbJ6l0hhSg0/DSbAr5k7Rx40ZuvPFGtmzZAsCSJUv40pe+xGWXXQbAl7/85WMnYeXESsUCw0MHGB06wPiRgxSHD1IZHaQ6Nojnh4gXhkgUD5MuHyFXGaalNky7j5CZXElxwnY9zqi1MGat5OOtFOOtjGWW82qynVq6A0+3Ect2Es91kmjpItXSSbati2zbErJtXbS0tpNbJCfiZHEx9+NuF9M0vb29/vqHhuzcuZMLLrigSRXNnlCOYyqv1RgbPcKRg/2MDL5Cfmg/5eEDVEcPYGMHSRYGSZWGyFaGaalOhHWLnfj2DiWPM2xtjMbaGI93UEy2U051Ust04dkuYrku4rkOUrkuUm1dZNu6aWnrorWzm3Qmh8X0NEtZnMzsMXfvrbdNI3mZxms1hg8PMvTqHkYG9lIYeoXqyAF8bIBEfpBUcZBceYi2ymE6/TCtVq67emPEsxyJdTIW72AstZSh1LlU00fDegmJ1m5Sbd1kO5bR0rmMtq7ltLR2sDQWY+mcH7VIuBTyi8j46BEG97/M8IG95A/1UTncDyOvkBx/lWxhgI7KAN21Q3RMnlycquQJhqyDkXgn48klHG45h93ZpdCylHjbclIdZ5DrWkFb90o6l66kLdtCW1OOUkSmWhAh7+4LenniXEyJea3G4IF9HOp/gZH9L1Ae3IMd2UtmbB/txf101wZoZ/y4Oee8pxiMdTOcXMorrRezN7cc2leS6DyTXPdq2rpX0rGsh7b2Ls6IxTgj8iMRkdk070M+k8kwODi4YG83fPR+8pnM6V58/ZqxkcPsf+lZDvftpHzgBWLDL5MZ76eztJ/l1QMstfK0qY5hWhiIL+dwZhUHcr3U2laS6DiTzJIe2pf10LVyLW3tXfRoLlskWPM+5Ht6eujr62NgYKDZpZy2o0+GOhmF/Bj7d+9kaO+vKL76PLFDL9AytodlpX0s5xDnTtl3kA4GE2dwMHcu/S3vwrrOIr10Le0rzmFpz3m0d3bTHs0hicgCMe9DPplMBvlEpUq5xL4Xn2XwxSco9j9D+tBzLBt/gTNrr7DWnLWT+w3SwYHkKvZ0Xs4LXeeQXr6OztUXsmLtm+lu7XjDdd8iIvM+5ENQyI/x8s5HGdr1CPbKkywZ/hWrKy9zlpU5i4mbNPXHVjKQO4++Je8nufx8Onou4IyzL6K7s1tBLiKnTSE/y7xWo3/3c/T/8qf4y9voPrKDNZU9nG9VAIZoY2/mfJ5YdgXxlRfTtfYSetZdwupc67SH4YqIzAaFfINq1SovPfsoB3f8jETfNlaP/pJVHGIVEyc+d2fezPYz3knmrPWsePMVrFi9ji6d6BSROaKQPw2Dr/bx0sM/hl1bOXf4Ec5lmHOBV+nm5ba38dLqK1h+8Xs4602X8paIHgQgInIyFPInad+LO9n7L3/Nkj3/wHnlX9NtziHaeaH9cvyc97Dqrdewcs06ztAoXUTmkchC3sz+GPggUAJeAD7m7oej6i8KB/fvZdfWu1j60o85r/oCq4BfJ9bx8Nr/zNK3/g7nvuXt9GqkLiLzWJQj+QeBW929YmZ3ArcCfxhhf7OiVq3y9D/9kNr2/83FY9u4wqo8nzifbed9jjXv+LesW/sm1jW7SBGRkxRZyLv7T6Z83AbcFFVfs6GQH+OX932bFc/exSW1fQzSwWMrb2blVf+R8998abPLExE5LXM1J/9x4Pv1NpjZJmATwJo1a+aonNeUS0Ue/7s/55xn/5INDLErfi7b13+Dt7z3Fq5IN34rAhGRZmoo5M1sK7Cizqbb3P1Hk/vcBlSAv6rXhrtvBjbDxP3kG6nnVD31j1vo/vlXuNxfYWfyQl59959z0W9+QPclF5FgNBTy7n7NG203s1uADwC/5fPo6SQH9+9lz/c+zfrRh9gTW82T7/g2l7zndxXuIhKcKFfXvI+JE63vdvfxqPo5Vc/84l5WbP0MF/sY/3r2J1m/8SucpWkZEQlUlHPyfwGkgQcnbxG8zd1/P8L+ZrTt/9zBhuf+hL3xVYzc9H+58sLLmlmOiEjkolxdc15UbZ8qr9XYdtdnubL/uzze+k7e9Pvfo6Wts9lliYhEblFc8brtf36GK1/5Hg9330DvJ79DPLEoDltEJPyQf/j7dx4L+A2f+l86uSoii0rQiffrJ3/O2569k6eyl9P7ye8o4EVk0Qk29Qr5MTI/2sSQdXLWJ76rKRoRWZSCDfknfvA1Vns/B67+UzqX1rteS0QkfEGG/NDAK1zy4l080fIOfuNdH252OSIiTRNkyP/q3j8lZ0WWfOC/N7sUEZGmCi7ky6Uib977fZ7MXclZF6xvdjkiIk0VXMjv+Pk9dDEM6/9Ds0sREWm64EK+8tTfMEQbF71Tc/EiIkGFfK1a5dzhh9nV8XaSqXSzyxERabqgQv6lZx+lixE4593NLkVEZF4IKuQHnn4QgDXrr2tyJSIi80NQIR9/9WkOsIQzes5tdikiIvNCUCHfNbqL/Zlzml2GiMi8EUzIV8olVldeZrzzTc0uRURk3ggm5A/se4m0lYktW9fsUkRE5o1gQv7Iq7sByC49q7mFiIjMI8GE/NjAHgDalyvkRUSOijzkzewLZuZmtjTKfipDfQAsOfPsKLsREVlQIg15M1sNXAu8HGU/AIwdpOBJ2jqWRN6ViMhCEfVI/pvAHwAecT/EikcYsdaouxERWVAiC3kz+xCwz92fmmG/TWa23cy2DwwMnHZ/ifIw4zGFvIjIVA09+NTMtgL1nq13G/BF4L0zteHum4HNAL29vac94k+WR8gr5EVEpmko5N39mnrfm9lvAGcDT5kZQA/wuJltcPf9jfR5IpnqKOOJziiaFhFZsBoK+RNx96eB5Uc/m9luoNfdD0bRH0C2OspwdnVUzYuILEjBrJNPeolqPNPsMkRE5pVIRvKv5+5ro+4jSRmPp6LuRkRkQQlmJJ+ggseSzS5DRGReCSbkU14GjeRFRKYJJuSTVDRdIyLyOkGEfKVcIm4OCT28W0RkqiBCvlwqTrzRSF5EZJogQr5ULABgCYW8iMhUQYR8uTQZ8hrJi4hME0TIe6068SY2J8v+RUQWjCBCXkRE6lPIi4gETCEvIhKwIELePfIHT4mILEhBhPxRk/euFxGRSUGFvIiITKeQFxEJmEJeRCRgQYS8e63ZJYiIzEtBhPwxOvEqIjJNWCEvIiLTRBryZvYZM3vOzHaY2Tei7EtERI4X2R29zOw9wPXAW9y9aGbLo+pLRETqi3Ik/0ngj9y9CODuB6LqSFe8iojUF2XInw+808weNrN/MrPL6u1kZpvMbLuZbR8YGGiwS514FRGZqqHpGjPbCqyos+m2yba7gCuAy4AfmNk5/rpht7tvBjYD9Pb2akguIjKLGgp5d7/mRNvM7JPA306G+iNmVgOWAo0O10VE5CRFOV3zd8DVAGZ2PpACDkbYn4iIvE6Uz8u7G7jbzJ4BSsAtr5+qmTU68SoiUldkIe/uJeD3omq/Ht1qWERkOl3xKiISMIW8iEjAggh53YVSRKS+IEL+KNecvIjINEGFvIiITKeQFxEJmEJeRCRgQYS8roUSEakviJB/jU68iohMFVjIi4jIVAp5EZGAKeRFRAIWRsjrzKuISF1hhPwk3YVSRGS6oEJeRESmU8iLiARMIS8iErAgQj6qpwqKiCx0QYT8a3TiVURkqshC3szeambbzOxJM9tuZhui6ktEROqLciT/DeCr7v5W4MuTn0VEZA5FGfIOtE++7wD6I+xLRETqSETY9ueAB8zsT5j4x+Q36+1kZpuATQBr1qw5za70jFcRkXoaCnkz2wqsqLPpNuC3gP/i7j80s98FvgNc8/od3X0zsBmgt7e3sWUyuuJVRGSahkLe3Y8L7aPM7LvAZyc//g1wVyN9iYjIqYtyTr4fePfk+6uBX0fYl4iI1BHlnPx/Av7MzBJAgcl5dxERmTuRhby7/wJYH1X70/qq6YpXEZF6grriVbcaFhGZLqiQFxGR6RTyIiIBCyLkHc3Ji4jUE0TIv0Zz8iIiUwUW8iIiMpVCXkQkYAp5EZGAhRHyrrtQiojUE0bIH6WLoUREpgkr5EVEZBqFvIhIwBTyIiIBCyLk3XXFq4hIPUGE/DE68SoiMk1YIS8iItMo5EVEAqaQFxEJWBAhr/OuIiL1NRTyZvYRM9thZjUz633dtlvNbJeZPWdm1zVW5kkXNCfdiIgsFI0+yPsZ4Ebg21O/NLMLgZuBi4Azga1mdr67VxvsT0RETkFDI3l33+nuz9XZdD2wxd2L7v4SsAvY0EhfIiJy6qKak18F7J3yuW/yu+OY2SYz225m2wcGBiIqR0RkcZpxusbMtgIr6my6zd1/dKJfq/Nd3dOj7r4Z2AzQ29t7eqdQdathEZG6Zgx5d7/mNNrtA1ZP+dwD9J9GO6fE9IxXEZFpopquuRe42czSZnY2sA54JKK+RETkBBpdQvlhM+sDrgTuM7MHANx9B/AD4FngfuBTWlkjIjL3GlpC6e73APecYNvXgK810v4pFDIn3YiILDRBXPF6jC6GEhGZJqyQFxGRaRTyIiIBU8iLiAQsiJD3+tdZiYgsekGE/GsCOxwRkQYpFUVEAqaQFxEJmEJeRCRgQYS864pXEZG6ggj5o3TBq4jIdEGFvIiITKeQFxEJmEJeRCRgYYS8Hv8nIlJXGCF/lM68iohME1bIi4jINAp5EZGAKeRFRALW6IO8P2JmO8ysZma9U76/1sweM7OnJ39e3XipJ6YLXkVE6mvoQd7AM8CNwLdf9/1B4IPu3m9mFwMPAKsa7GtGphOvIiLTNBTy7r4Tjg9Xd39iyscdQMbM0u5ebKQ/ERE5NXMxJ/9vgCcU8CIic2/GkbyZbQVW1Nl0m7v/aIbfvQi4E3jvG+yzCdgEsGbNmpnKqUt3oRQRqW/GkHf3a06nYTPrAe4B/r27v/AG7W8GNgP09vY2lNauxUIiItNEkopm1gncB9zq7v8SRR8iIjKzRpdQftjM+oArgfvM7IHJTZ8GzgO+ZGZPTr6WN1iriIicokZX19zDxJTM67+/A7ijkbZFRKRxYUxi6y6UIiJ1hRHyk3QxlIjIdEGFvIiITKeQFxEJmEJeRCRgQYS8rngVEakviJAXEZH6FPIiIgFTyIuIBEwhLyISsEBCXideRUTqCSTkJ5gFdTgiIg1TKoqIBEwhLyISMIW8iEjAggh5XfEqIlJfECF/jG41LCIyTVghLyIi0yjkRUQCppAXEQlYQyFvZh8xsx1mVjOz3jrb15jZqJl9oZF+ZqQTryIidTU6kn8GuBH45xNs/ybwDw32cfJ04lVEZJpEI7/s7juh/gO0zewG4EVgrJE+RETk9EUyJ29mLcAfAl89iX03mdl2M9s+MDAQRTkiIovWjCFvZlvN7Jk6r+vf4Ne+CnzT3Udnat/dN7t7r7v3Llu27FRqn9rGaf2eiEjoZpyucfdrTqPdy4GbzOwbQCdQM7OCu//FabR10jQlLyIyXUNz8ifi7u88+t7MvgKMRh3wIiJyvEaXUH7YzPqAK4H7zOyB2SlLRERmQ6Ora+4B7plhn6800oeIiJy+IK54zbV383jru2np7ml2KSIi80okc/Jzree8i+n5wr3NLkNEZN4JYiQvIiL1KeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYDafbtNrZgPAngaaWAocnKVyFoLFdrygY14sdMyn5ix3r3uv9nkV8o0ys+3uftyzZkO12I4XdMyLhY559mi6RkQkYAp5EZGAhRbym5tdwBxbbMcLOubFQsc8S4KakxcRkelCG8mLiMgUCnkRkYAFEfJm9j4ze87MdpnZf212PVEzs7vN7ICZPdPsWuaKma02s5+Z2U4z22Fmn212TVEzs4yZPWJmT00e81ebXdNcMLO4mT1hZv+v2bXMFTPbbWZPm9mTZrZ9Vtte6HPyZhYHngeuBfqAR4GN7v5sUwuLkJm9CxgFvuvuFze7nrlgZiuBle7+uJm1AY8BNwT+39mAFncfNbMk8Avgs+6+rcmlRcrMPg/0Au3u/oFm1zMXzGw30Ovus34BWAgj+Q3ALnd/0d1LwBbg+ibXFCl3/2fgULPrmEvu/oq7Pz75fgTYCaxqblXR8gmjkx+Tk6+FPSqbgZn1AL8D3NXsWkIRQsivAvZO+dxH4H/5FzszWwu8DXi4uZVEb3Lq4kngAPCgu4d+zN8C/gCoNbuQOebAT8zsMTPbNJsNhxDyVue7oEc7i5mZtQI/BD7n7sPNridq7l5197cCPcAGMwt2es7MPgAccPfHml1LE7zd3S8Ffhv41OSU7KwIIeT7gNVTPvcA/U2qRSI0OS/9Q+Cv3P1vm13PXHL3w8BDwPuaXEqU3g58aHJ+egtwtZl9r7klzQ1375/8eQC4h4lp6FkRQsg/Cqwzs7PNLAXcDNzb5Jpklk2ehPwOsNPd/0ez65kLZrbMzDon32eBa4BfNbeq6Lj7re7e4+5rmfh7/I/u/ntNLityZtYyuZgAM2sB3gvM2sq5BR/y7l4BPg08wMTJuB+4+47mVhUtM/tr4F+BN5lZn5l9otk1zYG3A/+OidHdk5Ov9ze7qIitBH5mZr9kYjDzoLsvmmWFi8gZwC/M7CngEeA+d79/thpf8EsoRUTkxBb8SF5ERE5MIS8iEjCFvIhIwBTyIiIBU8iLiARMIS8iEjCFvIhIwP4/G1tRxEUvzZ8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(grid, V0(grid), label='V0')\n", "plt.plot(grid, U0)\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfGklEQVR4nO3deXxU9b3/8dcnIRAIqyTKEiAIWAw7DJDU21a90mpbBWt/Fqx7K8UWq/f3aGvr7XLt9VZbf23tQkup2vvrdaHaorUWl9ZqbStLAgKyC5ElLBLCmoSs87l/JNAAQYYwyZkz834+HjzMnHPmzHseyttvvvOdc8zdERGR8EsLOoCIiMSHCl1EJEmo0EVEkoQKXUQkSajQRUSSRIegXjg7O9vz8vKCenkRkVBatmzZXnfPaWlfYIWel5dHcXFxUC8vIhJKZrb1VPs05SIikiRU6CIiSUKFLiKSJAKbQ29JXV0dpaWlVFdXBx2lVTIzM8nNzSUjIyPoKCKSghKq0EtLS+nWrRt5eXmYWdBxzoi7U15eTmlpKYMHDw46joikoISacqmurqZ3796hK3MAM6N3796h/e1CRMIvoQodCGWZHxXm7CISfgk15SIikmjqamvYvOJ19q9/HautjMs5s4b9C6M+9Im4nKs5FXoLdu/ezV133UVRURGdOnUiLy+Phx56iEWLFnHfffcB8PWvf52bbrop4KQiEm8ejbJt4wp2vfkimdtfZ2jlCobbEQCiHp/fwpdE60CF3vbcnauvvpqbbrqJ+fPnA7BixQp27drFvffeS3FxMWbGhAkTuOqqq+jVq1fAiUXkbO0v28WmN56BktcYdLCIQexjEFBqfViT/WEyhl7C+ROvoGd2n7i8XmFcznIyFfoJXn31VTIyMpg1a9axbWPHjuXJJ59kypQpnHPOOQBMmTKFF198kRkzZgQVVUTOUn1dLcVPf5eRG+Yw0Y6wn+6UdJvAlrwPkTv+CnIHDyc36JBnIGEL/d4/rGHtzkNxPWd+v+5868oR73nM6tWrmTBhwknbd+zYwYABA449zs3NZceOHXHNJyLtZ/2Sl+n40pcpiG5hVecInT/yLYaMvogJ6elBR2u1hC30RNPSvVe1qkUkfMrfLaXkyS8x8cAL7CabNwt/zNgpN2BpCbfo74wlbKGfbiTdVkaMGMFvf/vbk7bn5uby2muvHXtcWlrKxRdf3H7BROSsNNTXU7zgB1y49iHGejWL+t3ImE/fR5+uPYKOFjcJW+hBufTSS7nnnnv45S9/yW233QZAUVER/fr14+WXX2b//v0AvPzyy9x///1BRhVJGQf3lRFtqG/18999ZzXpL3+NyfVvs6bjGLp+4iEKh4+PY8LEoEI/gZnxzDPPcNddd/HAAw+QmZl5bNniN77xDSZOnAjAN7/5zWMfkIpI29i9fRPvPvF5xhxZclbn6QWU0YviyINM+Ohnk2J6pSUq9Bb069ePp5566qTtw4YN49Zbbw0gkUhqiTY0UPS77zNizQ/oTpRFA24lreu5rT6fZWRy4WU3EemR3IMwFbqIJJRtG1dQ8fQXmFy3mrcyx9F7+lwKBw8POlYoqNBFJCHU1dZQ/OS9jC+ZR0/ryNIx9zFx6heSdnqkLSRcobt7aJcDtrS0UUROb9PKv2PP3UFhQwnLu32QgdfPYVKfgUHHCp2EKvTMzEzKy8tDeQndo9dDz8zMDDqKSMJrqK/n7Tdf48DKhfTe/TpD6jaxz3qwvPCnjP/IDUHHC62YCt3MLgd+BKQDD7v7Ayfsvxj4PfBO06YF7v7tMw2Tm5tLaWkpZWVlZ/rUhHD0jkUicrK9u7dRsuj3pJe8wtDDSxlOJQ1uvN3xQpbkfY78qV9i/Dk5QccMtdMWupmlA3OAKUApUGRmz7n72hMO/Zu7f/xswmRkZOhuPyIJZu/ubbyz+DnSNv+Z3MOrSKfhjM9hONkcJBvYS0829vwgacMuY2jBlQzvfV78Q6eoWEbok4BN7l4CYGbzganAiYUuIkmgvq6Wjcv+wsG3XiBn998Y2rCZbBrXcW/tPoGGjKxWnXdj91xyxn6M80cWkB3i66UkslgKvT+wvdnjUmByC8cVmtlKYCfwJXdfc+IBZjYTmAkwcKA+8BBJNMXPz2NI8bfJ5zD1nsbGjvksGvgFzh33cc4fWUCOVpwktFgKvaVPJ09czrEcGOTuFWb2UeBZYNhJT3KfB8wDiEQiWhIikiCqKg6y+uFZTDqwkPUdLuSdCZ9jSMGV5PfKDjqanIFYCr0UGNDscS6No/Bj3P1Qs58XmtnPzCzb3ffGJ6aItJWS1UtIX3ArkYYdLMq9hYk3f48OGR2DjiWtEEuhFwHDzGwwsAOYDlzX/AAz6wO86+5uZpNovPl0ebzDikj8eDTK0qcfZOzaBzlsWay97P9T+IGpQceSs3DaQnf3ejObDbxE47LFR919jZnNato/F/gkcLuZ1QNHgOmub9mIJKyD+8rY/MgtTK78G6s6T6T/Lf/NyPO05DbsLKjejUQiXlxcHMhri4TNW39dwLmvfolefiAu50snSpQ0lg27g0kzvkmaVp2Ehpktc/dIS/sS6puiInK8aEMDS359D5O3/IJt6QMo6XNWX/X4JzOyJ36SgrEfiM/5JCGo0EUS1MF9Zbzzy+spPLKY4h6XkT/zUfKS6O46En8qdJEEtHnVG2Q+czMjontZcuFXmXTt3brqoJyWCl0kgdTV1vDmH37G6FX/xSHrxuaPP8XkiZcFHUtCQoUuErCG+nrWL3mJiuW/4YLyvzCJw6zpNJrzbn2C4X0GnP4EIk1U6CJtbH/ZLrav/gd+whesvb6O6rdfZciePzGC/VR5J9Z1v4i0Udcw6pJr9eUeOWMqdJE2Un2kkjef+g6jSh5htB1p8Zgaz2Bt18lsG/EJLvzgJ5mgDz3lLKjQReLMo1GWv/Ao/YoeoJAy3sx6PxkXfZ6MTidfpbDv0DGM69k7gJSSjFToInG0ofgv+Iv3MKF+HZvTB7P64u8zTl+nl3aiQheJg4P7ytjw2L8xad8f2EtPlo66lwlTZ5PeQX/FpP3ovzaRs+DRKMsWPszg4vsY74dZ3GcGIz99P5O69wo6mqQgFbpIK+0oWUf5U7OJVBezscMF7L/ySQrGXBR0LElhKnSRJvv27GDT/LvJOfBWTMf3bdhBDzqwePjdTPw/X9H0igRO/wVKyos2NFD87E+44K0HGedHWNNlIp52+r8aezLHknf1NynIHdIOKUVOT4UuKW3rumVULvgik+pWsy5jBJ0/8WPGXtjilUlFEp4KXVJSbU01yx77dyZs+xVV1pmlo+4lMu0OXRdcQk2FLiln86o3sGdvpzC6heIeUzj/+h8x6dz+QccSOWsqdEkZdbU1FD/2DSJbH+agdWPFRT8nMuW60z9RJCRU6BI6u7dvYtfaf1BXsY9o1X78yH7SaitO+7zs/SsobCihuMdlDL1xDmOz+7RDWpH2o0KXUNm08h+ct+AaxjW72FWtp1NlnXHsPZ9baV15s/DHRD5yU1vHFAmECl1CY/vbKznnmelUWhY7PvYYPc8bRLdeOXTJ6k7PGO7m0wvQfe0lmanQJRR2b99ExuPXAFD36QUMHzYm4EQiiUc3KZSEt2/PDmp/dRVZXsG+q+czQGUu0iKN0CVh7ShZw/Ylv6fP20/Qp2EPJVc8Rr6ulSJySip0CdSBvbtZ94cfYtUHjm2z+iP021fEAN9Jf2C79WPjJb9gdMHlwQUVCQEVugTCo1GKfj+HYSu/y2SvoIrMY/uilsaWzAvZMegGcidexYChI9GtkkVOT4Uu7W7rumVUPnMnk2rfYn1GPgemPcTgEZOPO2Z0QNlEwkyFLu1qxSvzGf76bGqsI0tH/QeRaV/U9VNE4iSmVS5mdrmZbTCzTWb21fc4bqKZNZjZJ+MXUZJF0bNzGPn67WzvMIi6WYuZdM2/qcxF4ui0hW5m6cAc4AogH5hhZvmnOO67wEvxDinht/jxe5m44h7WZ46m751/JrvPwKAjiSSdWKZcJgGb3L0EwMzmA1OBtSccdwfwO2BiXBNKKJTt3MK2Fa9St+dt0g+UkHlkN4YD0KGhhoL6dSzv+kFGzP4NnTK7BJxWJDnFUuj9ge3NHpcCx32CZWb9gauBS3mPQjezmcBMgIEDNUJLFjtK1tH11//KBCoBKKMX+zqcR9Qap1Pc0ljU7yYm3foD3aZNpA3F8rerpSse+QmPHwLudvcGs1NfIMnd5wHzACKRyInnkBA6UnmYmsem0xVYf8XTDMifRE63nuQEHUwkBcVS6KVw3DLgXGDnCcdEgPlNZZ4NfNTM6t392biklITk0ShrfnEz4xu2svrihxk9+cNBRxJJabEUehEwzMwGAzuA6cBxdwVw98FHfzaz/waeV5knvyXzv0PBoT+zaPDtFF6ihU0iQTttobt7vZnNpnH1SjrwqLuvMbNZTfvntnFGaWdHKg+zv+zEX8IaHXp3K/vX/5Uuu5cSqSrmzayLmHzDf7VzQhFpSUyfULn7QmDhCdtaLHJ3v/nsY0lQNq38B+c8M51+HGpxf7+mf25Ny2VZzjTyb/i+1pKLJAgtOZBj1he/Qr/nr6eKLJaOuhfSTi7qjKxeDBp7CYPO7c+gADKKyKmp0AWANW8sJO+lW9if1pP0m59j0qD3BR1JRM6QbnAhLFv4K4a+dCN703PodNtL9FWZi4SSRugpzKNRFj/2LQpLfsy6jvn0mbmAXjl9g44lIq2kQk8Bb6/4GxUv30+a1x23vVN9BYV1a1nW7VJGfP4xMjtnBZRQROJBhZ7k3i3dTK9nryebKHs7nHfS/kWDZjH5xu9opYpIElChJ7EjlYc59Ktr6es1lE//I8MunHDSMcMCyCUibUMfiiap+rpa1s69gSH1m9n0gR8yqIUyF5HkohF6Etq7exvvPvppJtSuYtGQL1J42YygI4lIO1ChJ4H6ulqWPfsTopVlEI0ydNtTnO+VFI37DoXTvhB0PBFpJyr0kPNolGVzb2Ny+T+vhbY1LZeKa37DxBNuvCwiyU2FHnJLnriXgvJnWdT3eiK3/ACAgR0ysDR9PCKSalToIVb8x19SsOkhlnf9EJM/+2MtPRRJcRrGhVTRs3MYt/TLrM0YSf7nn1CZi4hG6Inu8MF9rH3pEby++p8bD+2iYPfjrM4cy/l3PEdml67BBRSRhKFCT2C1NdVsmzONybUrT9r3Zpf3c+Hsp1XmInKMCj1BeTTKyp/fzMTalSwdcx/vu/ifa8nNjHE9eweYTkQSkQo9QS3+9T0UHniBRQNnUnj1HUHHEZEQ0IeiCaj4+XkUbvk5xd2nUHDzd4OOIyIhoRF6wLZuWMGuV36GResbN3gD4/Y+z9pOoxj1+V9rPbmIxEyFHqDtb68k68mpjPcKqizz2PbNHS+g/+d+R6fMLgGmE5GwUaEHZM+Od8h4/BOkEWXXda8w6H1jj+3rGWAuEQkv/T4fAI9G2f0/t9HdD7P/mqePK3MRkdZSoQeg6PdzGF1dxKrhdzFkVEHQcUQkSWjKpZ1UHNrP23Ovo2/VRsb6AdZ1HMGka+8OOpaIJBEVejs4fHAfO376MUbVrmdFj39la8duDLrqHl1/RUTiSoXeDjY8chtjajewqvCHRC6/Oeg4IpKkYppDN7PLzWyDmW0ys6+2sH+qma0ysxVmVmxm/xL/qOG0fumfiBz6M8UDbmK8ylxE2tBpR+hmlg7MAaYApUCRmT3n7mubHfYK8Jy7u5mNBp4ChrdF4ERXX1fLqleeIFrbeHXEnivnsYdzGP2pbwWcTESSXSxTLpOATe5eAmBm84GpwLFCd/eKZsdnAR7PkGGy7HffZ/L6B449jrqxfOL3iHTT6nIRaVuxFHp/YHuzx6XASTerNLOrgfuBc4GPtXQiM5sJzAQYOHDgmWZNeHW1NQxc/wjrM/LJ+tQ8ADK7dCPSLy/YYCKSEmKZQ7cWtp00Anf3Z9x9ODAN+M+WTuTu89w94u6RnJycM0saAiteeIS+lFFTeCcDho5iwNBR5KjMRaSdxFLopcCAZo9zgZ2nOtjdXweGmFn2WWYLlfVFf2bom/fzTloeoy++Nug4IpKCYplyKQKGmdlgYAcwHbiu+QFmNhTY3PSh6HigI1Ae77CJpK62hq3fLWRow2ag8RPgUutDhxm6QqKIBOO0he7u9WY2G3gJSAcedfc1Zjaraf9c4BrgRjOrA44An3L3pP5gdPnvf8rkhs0s6T2NaJds6JDJ8I/OpldO36CjiUiKsqB6NxKJeHFxcSCvfbaOVB7m8IOjKc/ow/Cv/UMjchFpN2a2zN0jLe1TE7XCql99kXPZh112r8pcRBKGvvp/BooW/Iiea/+HyfVvs/i8GRRM/nDQkUREjtHwMkYbl7/GuJX/QUa0hiW9pzH25u8HHUlE5Dgaoccg2tBAxvN3stfOodcdr5LXK6VWZIpISGiEHoNVrz3N4OgWSsd/iR4qcxFJUBqhv4eqioOsfuR28vYvYjc5jLn81qAjiYickgr9FKINDWz+yTQmVL/Jxo75HIncTp+OnYKOJSJySir0U1hf9CdG1Sxn8fCvUDDj34OOIyJyWppDP4XDSx+nyjsx+srZQUcREYmJRujNHCx/l5JHbqZX9XbGNuzmrR4XE+naI+hYIiIxUaE3U/qLaxhRs47VXd9Pedowel9+d9CRRERipkJvsvOd9YyofYtFQ75I4Y0tXs5dRCShaQ4dqDi0n21/exyAAe//VMBpRERaJ+VH6MXPzyNS/GUKgC1pA8gbOjLoSCIirZLyhd5h7QL2cA4lF3yGnJGXBh1HRKTVUrrQq6sqeF/lMlblXEnBdV8POo6IyFlJ2Tn0LeuKWf+Ta+hstWSNuzroOCIiZy1lR+hVC+5kbN1qFp97LQUXXRl0HBGRs5aSI/TqI5UMrV3Pkt5TKfj8L4OOIyISFylZ6JuXv0ZHq6fj+3THIRFJHik15VJVcZANP5vOuKo3OOydOT/ykaAjiYjETUqN0N9a+AvGVb3B0l4fY9+MhfQ4JyfoSCIicZNSI/S0HcWU04OJdzyGpaXU/8tEJAWkRKs11Nez6NEvc/7BxWzvkq8yF5GklBLNtvaNP1C4bR7pNFA/fFrQcURE2kRKTLlUrVhApWeSefcGIl26Bh1HRKRNJP0IfdfWDUze9xzru7+fTJW5iCSxmArdzC43sw1mtsnMvtrC/k+b2aqmP2+Y2Zj4R22dsidvByB9zLUBJxERaVunLXQzSwfmAFcA+cAMM8s/4bB3gA+5+2jgP4F58Q7aGh6NMrB6AyszJzLmUl3nXESSWywj9EnAJncvcfdaYD4wtfkB7v6Gu+9vergYyI1vzNZZ+Zff0JMKqodcrpUtIpL0Ymm5/sD2Zo9Lm7adymeAF1raYWYzzazYzIrLyspiT9kKWzesYOzfZ3HYOzP0g9Pb9LVERBJBLIVuLWzzFg80u4TGQm/x7sruPs/dI+4eyclp229p7nztEQD2fPJZep+XEL8wiIi0qViWLZYCA5o9zgV2nniQmY0GHgaucPfy+MRrnUMHyhmy63k2driAC0YVBBlFRKTdxDJCLwKGmdlgM+sITAeea36AmQ0EFgA3uPvG+Mc8MyW/uI5z2cf+XqODjiIi0m5OO0J393ozmw28BKQDj7r7GjOb1bR/LvBNoDfwMzMDqHf3SNvFfo+80SiDj7zFO2mDyL/+wSAiiIgEIqZvirr7QmDhCdvmNvv5s8Bn4xutdVa++hRjqWTDiK8wuMc5QccREWk3SbeWL+uNxlH5kH/5ZMBJRETaV9IV+rkNu1ie9QGtbBGRlJNUhb7q1d/Sg0pq+04IOoqISLtLqkLv9I/vAdC/QNMtIpJ6kqbQPRqlT30pS3pPY8DQUUHHERFpd0lT6Eue+DY9qMT6jQ06iohIIJKm0HM3z2c32Yz48C1BRxERCURSFHpDfT25vostfT5CVreeQccREQlEUhT6yleebPyhe99gg4iIBCgpCr2mdAUAo668I+AkIiLBSYpCH7H9Cbam5Wq6RURSWugL/eC+MrpTRXnnwUFHEREJVOgLfeNrTwCQPuH6gJOIiAQr9IXeecMCAPoPLww4iYhIsEJf6Fl1B1iVOYHsfoOCjiIiEqhQF3rZzi3kNmynqscFQUcREQlcqAv93XfeIsMa6DryiqCjiIgELtSFfnj9XwHo0qtPwElERIIX6kJPO1IOwKDhuv65iEhoC92jUSJlz7A5fTDpHWK6NaqISFILbaHXVFeRbk5lRu+go4iIJITQFnr57u0AHBn84YCTiIgkhvAW+tY1AHTNHRlwEhGRxBDaQq/eVwpAr/5DA04iIpIYQlvo0bKN1HgGOf3ygo4iIpIQQlvoWfvXsb3DQDI6dgo6iohIQghtoefUbOdA1vlBxxARSRihLHSPRunpB6jvcm7QUUREEkZMhW5ml5vZBjPbZGZfbWH/cDNbZGY1Zval+Mc8XmXFQTKtDrKy2/qlRERC47SFbmbpwBzgCiAfmGFm+Scctg/4IvD/4p6wBQf37gQgvZtG6CIiR8UyQp8EbHL3EnevBeYDU5sf4O573L0IqGuDjCc5XL4LgE49dFEuEZGjYin0/sD2Zo9Lm7adMTObaWbFZlZcVlbWmlMAcGT/bgC6nKNCFxE5KpZCtxa2eWtezN3nuXvE3SM5OTmtOQUAdYf2ANCtd99Wn0NEJNnEUuilwIBmj3OBnW0TJzYNFY2F3jNbhS4iclQshV4EDDOzwWbWEZgOPNe2sd6bVe7lEF3olNklyBgiIgnltBcSd/d6M5sNvASkA4+6+xozm9W0f66Z9QGKge5A1MzuAvLd/VBbhM6oLueg9aR7W5xcRCSkYrozhLsvBBaesG1us5930zgV0y7OrVhPRYee7fVyIiKhEMpvikYtnYxoTdAxREQSSigLPcNr2N91SNAxREQSSigLvZPXEO3QOegYIiIJJaSFXounZwYdQ0QkoYSy0DOpwTM0QhcRaS50hV5XW0MHi0IHjdBFRJoLXaEfqaoAwDrqS0UiIs2FrtBrqysBME25iIgcJ3SFXlPVWOhpGqGLiBwndIVeV3O00DVCFxFpLnyFXl0FQLpG6CIixwldoTfU1wKQlp4RcBIRkcQSukInGgXA0tMDDiIiklhCV+hRjzb9FLroIiJtKnyt6I13v7O0lu6MJyKSukJX6N405YKp0EVEmgtfodM0h26aQxcRaS50hX7sQ1GN0EVEjhO6Qm+aQsfSQhddRKRNha8VvQEAQyN0EZHmQlfofnSIrhG6iMhxQteK7lrlIiLSktAV+tFJ9LQ0rXIREWkudIV+bISuOXQRkeOErtD/+U3R8EUXEWlLoWtFjzatctEcuojIcUJX6EeZhTa6iEibiKkVzexyM9tgZpvM7Kst7Dcz+3HT/lVmNj7+UZs0rUNHhS4icpzTtqI1XjRlDnAFkA/MMLP8Ew67AhjW9Gcm8PM45zzGj61y0ZSLiEhzsQxzJwGb3L3E3WuB+cDUE46ZCvzaGy0GeppZ3zhnBcCjx7773xanFxEJrVhasT+wvdnj0qZtZ3oMZjbTzIrNrLisrOxMswKQlTOA5V0/ROduvVr1fBGRZNUhhmNamtvwVhyDu88D5gFEIpGT9sdi+MTLYOJlrXmqiEhSi2WEXgoMaPY4F9jZimNERKQNxVLoRcAwMxtsZh2B6cBzJxzzHHBj02qXAuCgu++Kc1YREXkPp51ycfd6M5sNvASkA4+6+xozm9W0fy6wEPgosAmoAm5pu8giItKSWObQcfeFNJZ2821zm/3swBfiG01ERM6E1v6JiCQJFbqISJJQoYuIJAkVuohIkrBj9+hs7xc2KwO2tvLp2cDeOMYJA73n1KD3nBrO5j0PcveclnYEVuhnw8yK3T0SdI72pPecGvSeU0NbvWdNuYiIJAkVuohIkghroc8LOkAA9J5Tg95zamiT9xzKOXQRETlZWEfoIiJyAhW6iEiSCF2hn+6G1cnGzB41sz1mtjroLO3FzAaY2atmts7M1pjZnUFnamtmlmlmS81sZdN7vjfoTO3BzNLN7E0zez7oLO3BzLaY2VtmtsLMiuN+/jDNoTfdsHojMIXGm2oUATPcfW2gwdqQmX0QqKDxnq0jg87THpruR9vX3ZebWTdgGTAtyf89G5Dl7hVmlgH8Hbiz6R69ScvM/i8QAbq7+8eDztPWzGwLEHH3NvkiVdhG6LHcsDqpuPvrwL6gc7Qnd9/l7subfj4MrKOFe9Qmk6YbrFc0Pcxo+hOe0VYrmFku8DHg4aCzJIuwFXpMN6OW5GFmecA4YEmwSdpe0/TDCmAP8Cd3T/b3/BDwFSAadJB25MDLZrbMzGbG++RhK/SYbkYtycHMugK/A+5y90NB52lr7t7g7mNpvCfvJDNL2ik2M/s4sMfdlwWdpZ1d5O7jgSuALzRNqcZN2ApdN6NOEU3zyL8DHnf3BUHnaU/ufgB4Dbg84Cht6SLgqqY55fnApWb2WLCR2p6772z65x7gGRqnkeMmbIUeyw2rJeSaPiB8BFjn7j8IOk97MLMcM+vZ9HNn4DJgfbCp2o67f83dc909j8a/x39x9+sDjtWmzCyr6UN+zCwL+DAQ19VroSp0d68Hjt6weh3wlLuvCTZV2zKzJ4FFwPvMrNTMPhN0pnZwEXADjaO2FU1/Php0qDbWF3jVzFbROHD5k7unxFK+FHIe8HczWwksBf7o7i/G8wVCtWxRREROLVQjdBEROTUVuohIklChi4gkCRW6iEiSUKGLiCQJFbqISJJQoYuIJIn/BaXuUVAhThPjAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(grid, C0(grid), label='C0')\n", "plt.plot(grid, c0)\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":10: NumbaWarning: \n", "Compilation is falling back to object mode WITH looplifting enabled because Function \"solvenb\" failed type inference due to: Invalid use of type(CPUDispatcher()) with parameters (readonly array(float64, 1d, C), Literal[int](1))\n", "Known signatures:\n", " * (array(float64, 1d, A), float64) -> array(float64, 1d, A)\n", " * parameterized\n", "[1] During: resolving callee type: type(CPUDispatcher())\n", "[2] During: typing of call at (13)\n", "\n", "\n", "File \"\", line 13:\n", "def solvenb():\n", " \n", " count=0.0\n", " V0 = U(grid,sigma)\n", " ^\n", "\n", " @nb.jit('f8[:]()')\n", ":10: NumbaWarning: \n", "Compilation is falling back to object mode WITHOUT looplifting enabled because Function \"solvenb\" failed type inference due to: cannot determine Numba type of \n", "\n", "File \"\", line 13:\n", "def solvenb():\n", " \n", " count=0.0\n", " V0 = U(grid,sigma)\n", " ^\n", "\n", " @nb.jit('f8[:]()')\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:742: NumbaWarning: Function \"solvenb\" was compiled in object mode without forceobj=True, but has lifted loops.\n", "\n", "File \"\", line 11:\n", "@nb.jit('f8[:]()')\n", "def solvenb():\n", "^\n", "\n", " self.func_ir.loc))\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:751: NumbaDeprecationWarning: \n", "Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.\n", "\n", "For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit\n", "\n", "File \"\", line 11:\n", "@nb.jit('f8[:]()')\n", "def solvenb():\n", "^\n", "\n", " warnings.warn(errors.NumbaDeprecationWarning(msg, self.func_ir.loc))\n", ":10: NumbaWarning: \n", "Compilation is falling back to object mode WITHOUT looplifting enabled because Function \"solvenb\" failed type inference due to: Invalid use of type(CPUDispatcher()) with parameters (readonly array(float64, 1d, C), array(float64, 1d, C))\n", " * parameterized\n", "[1] During: resolving callee type: type(CPUDispatcher())\n", "[2] During: typing of call at (15)\n", "\n", "\n", "File \"\", line 15:\n", "def solvenb():\n", " \n", " while count\", line 13:\n", "def solvenb():\n", " \n", " count=0.0\n", " V0 = U(grid,sigma)\n", " ^\n", "\n", " self.func_ir.loc))\n", "/Users/ozak/anaconda3/envs/GeoPython36env/lib/python3.6/site-packages/numba/compiler.py:751: NumbaDeprecationWarning: \n", "Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.\n", "\n", "For more information visit http://numba.pydata.org/numba-doc/latest/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit\n", "\n", "File \"\", line 13:\n", "def solvenb():\n", " \n", " count=0.0\n", " V0 = U(grid,sigma)\n", " ^\n", "\n", " warnings.warn(errors.NumbaDeprecationWarning(msg, self.func_ir.loc))\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "138.0\n", "138.0\n", "138.0\n", "138.0\n", "138.0\n", "138.0\n", "138.0\n", "138.0\n", "16.2 s ± 72.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "@nb.jit('f8(f8[:], f8[:])')\n", "def errnb(V0, V1):\n", " maxi=0.0\n", " for i in range(V0.shape[0]):\n", " m=abs(V1[i]-V0[i])\n", " if m>=maxi:\n", " maxi=m\n", " return maxi\n", "\n", "@nb.jit('f8[:]()')\n", "def solvenb():\n", " count=0.0\n", " V0 = U(grid,sigma)\n", " while count