{
 "metadata": {
  "name": "ScientificProgrammingI"
 },
 "nbformat": 3,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "source": [
      "Scientific Programming I<br>\n",
      "Berian James / berian@berkeley.edu"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "source": [
      "1. Introduction to SciPy"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "Import NumPy and SciPy"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "import numpy as np"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "import scipy as sp"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Differences in temperament between NumPy and SciPy; see, e.g., http://www.scipy.org/NegativeSquareRoot"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "exp(pi*np.sqrt(-1)) + 1"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "exp(pi*sp.sqrt(-1)) + 1"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "source": [
      "Getting data in and out of SciPy"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "Import Matlab data into Python"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "import scipy.io as sio"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "struct = sio.loadmat('testbox.mat') # Imports Matlab data structure"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "struct"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "box = struct['box'] # Extracts data object from structure"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "**Exercise: What are the dimensions of the imported array?**"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "source": [
      "Vector and array manipulation"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "Construct sequence of n (linearly-spaced) numbers, from a to b"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "a = 50; b = 100-0.1; n = 57;"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "sp.linspace(a,b,n)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Construct sequence of n (base-10 logarithmically-spaced) numbers between $10^a$ and $10^b$"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "a = -1; b = 1; n = 20;"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "sp.logspace(a,b,n)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Construct coordinate array"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "x,y = np.mgrid[0:5,0:5]"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "x"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "y"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "np.sqrt(x**2 + y**2)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Construct tiled array"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "x = np.linspace(0,10,11);"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "x"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "np.c_[x,x**2]"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "source": [
      "Symbolic computation with SymPy"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "Import SymPy (CAUTION!)"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "from sympy import *"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "1/2 + 1/3"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "Rational(1,2) + Rational(1,3)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "Rational(5,6).evalf(6)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "1./2. + 1./3."
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 4,
     "source": [
      "Calculus with symbolic variables"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "x = Symbol(\"x\")"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "limit( sin(x) / x, x, 0)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "diff(erf(x),x)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "integrate(1./sqrt(pi) * exp(-x**2),x)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "integrate(1./sqrt(pi) * exp(-x**2), (x,-oo,0) ).evalf(3)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 4,
     "source": [
      "Matrix compuation"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "M = Matrix( ([1, x], [x, 1]) )\n",
      "pprint(M)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "%load_ext sympyprinting"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "M.inv()"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "M.cholesky()"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "source": [
      "Interfacing with other languages"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "See files in f2py_example. Things to note:\n"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "*i) differences between subroutines and functions in Fortran are not really important here*"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "*ii) passing arrays is fine, but you need to pass their dimensions as well*"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "*iii) In general, SciPy contains many useful functions that work as quickly as Fortran/C. But if you need to crunch something inside a loop, or nested loops, a considerable speed-up is possible*"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "more example.f90"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!f2py --fcompiler=gfortran -c example.f90 -m example"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "import example"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "a = 14; b = 78;"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "example.arithmetic(a, b)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "c, d = example.arithmetic(a, b)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "A weave example (courtesy of Nat Butler) is below"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "from scipy import weave"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "x = np.arange(1000)/1000.\n",
      "y = np.zeros(1000,dtype=np.double)\n",
      "n=len(x)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "code = \"\"\"\n",
      "    int i;\n",
      "    for (i=0;i<n;i++) {\n",
      "      if (*x<0.5) *y = exp(-(*x)*2);\n",
      "      else *y = exp(-(*x));\n",
      "      x++; y++;\n",
      "    }\n",
      "\"\"\""
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "weave.inline(code,['n','x','y'])"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "print y"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "plot(x,y)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "source": [
      "2. SciPy Packages"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "source": [
      "I. Integration and interpolation"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "source": [
      "The SciPy interpolation routines are available through the scipy.interpolate module. "
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "import scipy.interpolate"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "scipy.interpolate."
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "from scipy.interpolate import interp1d"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Create (x,y) points following the function y = x**2"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "x = linspace(0,5,10)\n",
      "y = x**2"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Interpolate routines return a *function* that accepts input at new abscissa (x) locations"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "yfn = interp1d(x,y, kind='cubic')\n",
      "#yfn = interp1d(x,y, kind='cubic', bounds_error = False, fill_value=0)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Plot the interpolated function along a denser line"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "xs = linspace(0,5,100)\n",
      "ys = yfn(xs)\n",
      "plot(x,y,'k.',xs,ys,'k--')"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "The scipy.intergrate module provides numerical integration routines, i.e., for summing functions across a 1-, 2- or 3-dimensional domain."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "from scipy import integrate"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "The following lambda function returns the value $F(x_0) = \\int_0^{x_0} y(x)\\ dx$ for the function $y(x)$ we defined above with interpolate."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "inty = lambda x0: integrate.quad(yfn,0,x0)\n",
      "#inty = lambda x0: integrate.quad(yfn,0,x0)[0]"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "inty(3)[0]"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Alternatively, we could have used trapezoidal integration over the original points, with no iterpolations. Would this make a difference?"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "scipy.integrate.trapz?"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "**II. Random sampling and FFTs**"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "from numpy.random import randn"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "print randn()"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "y = randn(25)\n",
      "print y"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "from scipy.fftpack import fft"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "fy = fft(y)\n",
      "print fy"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "semilogy(range(len(y)),fy * conj(fy))"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "sp.fftpack?"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "**III. Optimization and linear algebra**"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "SciPy provides numerical optimization methods that can tackle a broad range of fitting problems. The simplest generic optimization routine in SciPy 0.10 is fmin (which uses Simplex optimization). "
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "from scipy.optimize import fmin"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Optimization routines operate on functions, finding the extreme (minimum) value location in parameter space by stepping through the space in a clever fashion, testing the value of the function at each point and using that information to progress to the next point."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "x = linspace(-5,5,100)\n",
      "y = lambda x: x**3 - x**2 - x\n",
      "plot(x,y(x))"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "fmin(y, 5)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "fmin(y,-4) # Understand the concept of local minima"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Linear algebra is handled in SciPy by the scipy.linalg module. The module wraps Fortran routines from LAPACK into Python."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "from scipy.io import loadmat\n",
      "data = loadmat('C.mat')\n",
      "C = data['C']\n",
      "imshow(C)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Matrix inversion"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "from scipy.linalg import inv\n",
      "imshow(inv(C))"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Matrix eigendecomposition"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "from scipy.linalg import eigh\n",
      "E,V = eigh(C)"
     ],
     "language": "python",
     "outputs": []
    }
   ]
  }
 ]
}