502976b3 update CI cloud image

c72d1f55 · Repo Updater · cd556732 · c72d1f55
Commit c72d1f55 authored Nov 14, 2022 by Repo Updater
--- a/notebooks/2a--NumPy.ipynb
+++ b/notebooks/2a--NumPy.ipynb
@@ -1144,10 +1144,11 @@
    "det = linalg.det(A)             # determinant\n",
    "P = linalg.pinv(A)              # pseudoinverse (accepts rectangular matrices)\n",
    "\n",
-    "# Further sparse matrices creation\n",
+    "# Further sparse matrix operations\n",
    "A = sparse.random(N, N, density=0.01, format='csr')\n",
    "B = sparse.random(N, N, density=0.01, format='csr')\n",
-    "C, D = sparse.diags(np.arange(4), 1), sparse.kron(B, A)\n",
+    "C = sparse.diags(np.arange(4), 1)\n",
+    "D = sparse.kron(B, A)               # kronecker product of matrices\n",
    "ev, evec = splinalg.eigs(A, 2)      # 2 largest eigenvalues\n",
    "u, s, vh = sparse.linalg.svds(A, 2) # 2 largest singular values"
   ]


 %% Cell type:markdown id: tags:
 # NumPy and SciPy
 **Python for HPC**
 Max Planck Computing and Data Facility, Garching
 %% Cell type:markdown id: tags:
 ## Outline
 * NumPy
 * SciPy examples
 * Input/output with NumPy arrays
 * Bonus: Image processing examples
 %% Cell type:markdown id: tags:
 ## Performance Problem - Loops
 * Looping in Python (and most interpreted languages) can be **very expensive**
 * Numerical codes perform mostly array computation
 * Example: Discrete Fourier Transform
 $$ X_k = \sum_{n=0}^{N-1}x_n\cdot e^{\frac{-2i\pi}{N}kn} $$
 %% Cell type:code id: tags:
 ``` python
 import cmath
 import random
 def DFT_loop(x):
    N = len(x)
    X = [0 for n in range(N)]
    for k in range (N):
        for n in range(N):
            X[k] += x[n]*cmath.exp(-2j * cmath.pi * k * n / N)
    return X
 X = [random.random() for i in range(512)]
 F = DFT_loop(X)
 ```
 %% Cell type:markdown id: tags:
 ### Solution - Numpy Arrays
 <div class="row">
  <div class="col-md-6" markdown="1">
 * Possible solution: turn loop coefficients into arrays and then use linear algebra operations
 * <span style="color:red">Larger memory footprint</span>, but <span style="color:green">higher performance</span>
 * much faster than loops
  </div>
  <div class="col-md-6" markdown="1">
      <center>
        <img src="fig/numpy_fft.svg?1", style="width: 85%">
      </center>
  </div>
 </div>
 $$ X_k = \sum_{n=0}^{N-1}x_n\cdot e^{\frac{-2i\pi}{N}kn} $$
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 def DFT_arrays(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)
 X = np.random.random(512)
 F = DFT_arrays(X)
 ```
 %% Cell type:markdown id: tags:
 ### True Solution - External Call
 * low-level external libraries will almost always be *much* better!
 * calls highly optimized low level (e.g. C, Fortran) functions with OpenMP, cache-hit, vectorization, etc
 * NumPy and SciPy already interface FFT (and many other) libraries
 * How do their performance compare?
 $$ X_k = \sum_{n=0}^{N-1}x_n\cdot e^{kn(-2i\pi /N)} $$
 %% Cell type:code id: tags:
 ``` python
 X = np.random.random(512)
 %timeit Fl = DFT_loop(X)
 %timeit Fa = DFT_arrays(X)
 %timeit Fe = np.fft.fft(X)
 ```
 %% Cell type:markdown id: tags:
 ## Bread & Butter: NumPy & SciPy
 * NumPy:
  * efficient handling of arrays operations via NumPy Arrays
  * executes loops in compiled C code (*much* faster!)
 * SciPy:
  * more advanced mathematical and numerical routines
  * provides several flavours of (sparse) matrices computation
 * Both make calls to BLAS and LAPACK routines
  * `numpy.linalg` and `scipy.linalg`
 %% Cell type:markdown id: tags:
 ## NumPy arrays
 * Main workhorse of NumPy $\to$ stores all kinds of data
 * plain C arrays with the plain values
 * preserves the intuitive syntax from Python
 * may use CPU performance features such as prefetching, caching, vectorization
 * may use high-performance libraries such as MKL (Math Kernel Library by Intel), Lapack, BLAS under the hood
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 my_array = np.array([[1,2,3,4], [5,6,7,8]])
 print(my_array)
 print(type(my_array))
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays - Constructors
 <br>
 * often require "shape" tuple (the parenthesis are mandatory)
 * accepts multi-dimension
 %% Cell type:code id: tags:
 ``` python
 x = np.zeros((3,4))
 x = np.ones((2,3,4))
 x = np.full((2,2),7.2)
 # random arrays
 x = np.random.random((2,2))      # 2x2 matrix with random reals
 x = np.random.randint(10,size=5) # 5 random integers between 0 and 9
 # use uninitialized memory
 x = np.empty((3,2))
 # start, stop, stride
 x = np.arange(10,55,5)
 print("Discrete Range:" + str(x))
 # first element, last element, total number of elements
 x = np.linspace(0,1,9)
 print("Continuous Range:" + str(x))
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays - Datatypes & Compatibility
 * datatypes map to native C types
 * datatype can be specified via the `dtype=` parameter of  constructor
 * operations automatically <span style="color:red">cast compatible types</span>
 * defaults (in most Linux/x86_64): `np.float64` and `np.int64`
 * https://docs.scipy.org/doc/numpy/user/basics.types.html
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 x = np.array([1,2,3])
 y = np.array([1.,2.,3.])
 print("x and y type: %s and %s" % (x.dtype, y.dtype))
 r = x+y
 print("x+y type: %s" % (r.dtype))
 z = np.array([1+2j,2+3j,3+4j])
 sz = z.astype(np.complex64)
 print("z and sz type: %s and %s" % (z.dtype, sz.dtype))
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays - Functions
 Arrays have a large selection of predefined operations. There are too many to list here!
 %% Cell type:code id: tags:
 ``` python
 x = np.random.randint(10,size=8)
 print("* Unsorted array: ", x)
 x.sort()
 print("* Sorted array: ", x)
 # compute some properties
 y = np.random.random((10,3))
 print("* Max, min, mean, std, sum: ", y.max(), y.min(), y.mean(), y.std(), y.sum())
 # may also specify over which axis
 print("* Max of each column: ", y.max(axis=0))
 print("* Max of each row: ", y.max(axis=1))
 # comparing two arrays
 z = y + 1e-16
 print("* Equal: ", np.array_equal(y,z))
 print("* Close: ", np.allclose(y, z, 1e-8))
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays - Functions Performance
 * prefer to use functions if they are there: <span style="color:green">faster</span> and <span style="color:green">cleaner</span>
 %% Cell type:code id: tags:
 ``` python
 # use numpy functions
 def mean_loop(a):
    result = 0
    for i in range(len(a)):
        result += a[i]
    return result/len(a)
 x = np.arange(2**20)
 %timeit mean_loop(x)
 %timeit np.mean(x)
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays - Attributes
 * rather advanced
 * useful if you are writing interfaces between Python and C/Fortran
 * useful for <span style="color:green">asserting contiguity</span> of data
 %% Cell type:code id: tags:
 ``` python
 print(y.ndim)  # number of dimensions
 print(y.size)  # number of elements
 print(y.dtype) # information about the data type
 print(y.flags) # information about memory layout
 print(y.itemsize) # size of one array element in bytes
 print(y.nbytes) # total size of the array in bytes
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays - Basic Slicing
 * syntax identical to lists
 * does <span style="color:red">not guarantee contiguity</span> of data
 %% Cell type:code id: tags:
 ``` python
 x = np.arange(20)
 # basic slice syntax as for python lists is start:end:step
 # for multidimension [start1:end1:step1, start2:end2:step2, ...]
 print("From 2 to 9, in steps 2:", x[2:9:2])
 print("All except 3 last elements:", x[:-3])
 # revert array
 print("Backwards:", x[::-1])
 # Reshape into multi-D arrays
 y = x.reshape((5,4))
 print("Shape of y:", y.shape)
 print("Slice of y:", y[1:,:2])
 # new dimensions can be added
 z = np.arange(3)
 print("z:")
 print(z[:,np.newaxis])
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays - Broadcasting
 * some dimension conversion is automatic, some are not!
 * careful with <span style="color:red">readability</span> of the code
 %% Cell type:markdown id: tags:
 <center>
    <img src="fig/numpy_broadcast.svg" style="width: 60%">
 </center>
 %% Cell type:markdown id: tags:
 ### NumPy Arrays - Broadcasting
 * careful with <span style="color:red">readability</span> of the code
 %% Cell type:code id: tags:
 ``` python
 X = np.arange(12).reshape(4,3)
 Y = np.arange(3)/10
 Z = np.arange(4)/10
 print("X is a 4x3 matrix:\n",X)
 print("\nX+Y is valid:\n",X+Y)
 # X+Z is invalid! We need to add a new dimension:
 print("\nX+Z[:,None] is valid:\n",X+Z[:,None])
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Array - Advanced Indexing
 * the selection object (as in `x[obj]`) can be:
    * a non-tuple sequence object (i.e., list or range) or ndarray of indices
    * or a tuple with at least one non-tuple squence object or ndarray
    * a boolean array $\to$ must have same dimensions
 %% Cell type:code id: tags:
 ``` python
 x = np.arange(10)*10
 print(x)
 # use integer list as indices
 print("X with integer indices: ", x[[2,5,8]])
 print("X with boolean indices: ", x[[False, False, True, False, False, True, False, False, True, False]])
 # create boolean array
 y = x**2 - 750
 index = (x < y)
 print(index)
 print("X with boolean indices: ", x[index])
 print("X with negated indices: ", x[~index])
 print("X with oneliner: ", x[x<y])
 # Set x to 100 for y below 10 - equivalent to:
 # for i in range(len(x)):
 #     if y[i]<10:
 #         x[i] = 100
 # Indexing can be used to assign values to indexed arrays
 x[x < y] = 100
 print("Oneliner x: ", x)
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays vs Views
 * Basic slicing on arrays usually returns a view on the original array
 * Data is **not** copied $\to$ changing the original array also changes the view
 * Views may be <span style="color:red">non-contiguous data</span>
 * use `copy()` for copying the data over
 * Advanced indexing directly returns a copy of the indexed data
 %% Cell type:code id: tags:
 ``` python
 x = np.arange(15)
 print("Array x: ", x)
 y = x[1:8:2]
 print("View y: ", y)
 print(y.flags)
 x[1] = 20
 print("View y after changing x: ", y)
 y[0] *= 10
 print("Original x after changing view y: ", x)
 c = y.copy()
 print("Copy of y: ", c)
 print(c.flags)
 ```
 %% Cell type:markdown id: tags:
 ### NumPy Arrays vs Views: Performance
 * strided data access may <span style="color:red">hinder performance</span>.
 * temporary array **may** be preferable, but: additional cost due to copy!
 * careful with passing views as argument
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 def sum_test(c,x,y):
    # Computes c*x^2 + y
    if len(c) == len(x) == len(y):
        return c*x*x + y
    else:
        print("Wrong shape")
        return None
 def sum_test_copy(c,x,y):
    copy = x.copy()
    return sum_test(c,copy,y)
 x = np.random.random(1000000)
 y = np.random.random(100000)
 c = np.random.random(100000)
 view = x[::10]
 copy = view.copy()
 %timeit res = sum_test(c,view,y)
 %timeit res = sum_test(c,copy,y)
 %timeit res = sum_test_copy(c,view,y)
 ```
 %% Cell type:markdown id: tags:
 ### Example 1:  Midpoint Computation
 $$y_i = \dfrac{x_{i}+x_{i+1}}{2}, \qquad \forall i<N$$
 %% Cell type:code id: tags:
 ``` python
 def midpoints_list(x):    # worst implementation
    y = []
    for i in range(len(x)-1):
        y.append(0.5*(x[i] + x[i+1]))
    return np.array(y)
 def midpoints_loop(x):    # loop forcing numpy array
    N = len(x)
    y = np.zeros(N-1)
    for i in range(N-1):
        y[i] = 0.5*(x[i] + x[i+1])
    return y
 def midpoints_array(x):   # best: no loops at all
    return 0.5*(x[:-1] + x[1:])
 x = np.linspace(0, 10, 10**6)
 %timeit midpoints_list(x)
 %timeit midpoints_loop(x)
 %timeit midpoints_array(x)
 ```
 %% Cell type:markdown id: tags:
 ### Example 1:  Midpoint Computation Explanation
 $$y_i = \dfrac{x_{i}+x_{i+1}}{2}, \qquad \forall i<N$$
 <br>
 * Loop ```for i in range(N-1):``` to left-hand side ```y[0:N-1]```
 * Translate indices on right-hand side to ranges for array
    * ```x[i]``` to ```x[0:N-1]``` and  ```x[i+1]``` to ```x[1:N]```
 %% Cell type:markdown id: tags:
 <center>
    <img src="fig/numpy-indexing.svg", style="width: 50%">
 </center>
 %% Cell type:markdown id: tags:
 ### Example 2: In-Place Manipulations
 * `x[1:] + x[:-1]` will create a temporary array
 * we can avoid it, but does that pay off?
 %% Cell type:code id: tags:
 ``` python
 def temp_arrays(x):
    return 0.5*(x[1:] + x[:-1])
 def in_place(x):
    y = x[1:]
    y += x[:-1]
    y *= 0.5
    return y
 x = np.random.rand(2**15)
 %timeit temp_arrays(x)
 %timeit in_place(x)
 ```
 %% Cell type:markdown id: tags:
 * performance is problem- and size-dependent
 * $\to$ may be faster, but not necessarily!
 %% Cell type:markdown id: tags:
 ### Example 3: Center of Mass
 * For positions $\vec{x}_i$ ($j$-th coordinate $x_{i,j}$), masses $m_i$:
 $$ \vec{x}_\mathsf{com} = \frac{1}{M} \sum_{i=0}^N m_i \vec{x}_i, \qquad\text{where}\qquad M=\sum_i m_i$$
 %% Cell type:code id: tags:
 ``` python
 def center_of_mass_loop(pos, mass):
    com = np.zeros(3)
    for j in range(3):
        for i in range(pos.shape[0]):
            com[j] += pos[i, j] * mass[i]
    total_mass = 0.0
    for i in range(pos.shape[0]):
        total_mass += mass[i]
    for j in range(3):
        com[j] /= total_mass
    return com
 def center_of_mass_arrays(pos, mass):
    com = np.zeros(3)
    for j in range(3):
        com[j] = np.sum(pos[:, j]*mass[:])
    total_mass = np.sum(mass)
    com /= total_mass
    return com
 def center_of_mass_oneline(pos, mass):
    return np.sum(pos * mass[:, None], axis=0)/np.sum(mass)
 # positions and masses of particles
 N = 10**5
 pos = np.random.rand(N, 3)
 mass = np.random.rand(N)
 %timeit center_of_mass_loop(pos, mass)
 %timeit center_of_mass_arrays(pos, mass)
 %timeit center_of_mass_oneline(pos, mass)
 ```
 %% Cell type:markdown id: tags:
 ## NumPy Summary
 * Takeaway messages:
    1. High level manipulations $\to$ in Python
    2. Intensive computation $\to$ delegate to underlying libraries
    3. No dedicated library $\to$ use NumPy arrays and their functions
 <br>
 * Tradeoff: high development productivity versus performance
 * NumPy code *can be* as fast as C or Fortran code, but often it is a factor of 2-4 slower
 $\to$ Implement hotspots in C/Fortran + Python Interface
 %% Cell type:markdown id: tags:
 # SciPy
 * SciPy is a mathematical toolbox built on NumPy, providing functions for
    * FFT
    * integration
    * interpolation
    * linear algebra
    * sparse matrix operations
    * optimization
    * special functions (e.g. Bessel functions)
 * ... and many more
 * most of them delegate <span style="color:red">heavy work</span> to C/Fortran libraries
 * [Documentation](https://docs.scipy.org/doc/scipy/reference/)
 %% Cell type:markdown id: tags:
 ## SciPy Example 1: FFT & Power Spectrum
 * calculate and plot the power spectrum of a signal
 %% Cell type:code id: tags:
 ``` python
 # power spectrum (code taken from the SciPy tutorial)
 import numpy as np
 from scipy.fftpack import fft
 import matplotlib.pyplot as plt
 N = 600  # number of sample points
 T = 1.0 / 800.0  # sample spacing
 x = np.linspace(0.0, N*T, N)
 y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
 yf = fft(y)
 xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
 plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
 plt.grid()
 ```
 %% Cell type:markdown id: tags:
 ## SciPy Example 2: Integration
 * Numerical integration of
    * a given function object
    * a given sample of function values
    * ODEs (ordinary differential equations)
 %% Cell type:code id: tags:
 ``` python
 import scipy.integrate as integrate
 result = integrate.quad(lambda x: x**2, 0, 2.5) # uses Fortran's QUADPACK
 print(result[0], 2.5**3/3.)
 x = np.linspace(0, 2.0, 10)
 y = np.sin(x)
 result = integrate.simps(y, x) # Simpson rule for given sample
 print(result, 1.-np.cos(2.0))
 # simple derivative: exponential
 def yprime(y, t, alpha):
    return np.exp(alpha*t)*alpha
 ts = np.linspace(0, 2.0, 50)
 y0 = 1.0
 y = integrate.odeint(yprime, y0, ts, args=(2.0,))  # uses Fortran's ODEPACK
 print(y[-1], np.exp(2.0 * ts[-1])) # compare to analytical result
 ```
 %% Cell type:markdown id: tags:
 ## SciPy Example 3: Sparse Matrices
 * different sparse formats for different tasks:
    * CSR: compressed sparse row, for row operations/slicing
    * CSC: compressed sparse column, for column operations/slicing
    * LIL: list-of-lists, efficient for building the matrix, slow for most other operations
    * DIA: diagonal matrix, very useful for FD
    * `speye()`: sparse identity
 %% Cell type:markdown id: tags:
 ### SciPy Example 3: Sparse Matrices
 %% Cell type:code id: tags:
 ``` python
 from scipy import sparse
 import numpy as np
 def fill_matrix(M):
    for i in range(300):
        rows = np.random.randint(0,N,300)
        cols = np.random.randint(0,N,300)
        val = np.random.random()
        M[rows,cols] = val
    return M
 N = 10000
 Mlil = sparse.lil_matrix((N,N))
 Mcsr = sparse.csr_matrix((N,N))
 %time Mcsr = fill_matrix(Mcsr)
 %time Mlil = fill_matrix(Mlil)
 %time Mlil = fill_matrix(Mlil).tocsr()
 ```
 %% Cell type:markdown id: tags:
 ## SciPy - Linear Algebra
 * handles dense and sparse matrices
 * delegates computation to lower level libraries
 * e.g. Anaconda Python comes with Intel MKL
 * SIMD instructions ("vectorization")
 * parallel processing ("threading")
 <br>
 * NumPy vs SciPy `linalg` module:
    * NumPy: uses own implementation (slower) if Fortran Compiler unavailable
    * SciPy: requires Fortran Compiler and has more functions
 %% Cell type:markdown id: tags:
 ## SciPy Example 4: Sparse vs Dense Linear Systems
 * dense: `scipy.linalg`, sparse: `scipy.sparse.linalg`
 * sparse matrix operations can be much faster, depending on the density of matrix
 %% Cell type:code id: tags:
 ``` python
 from scipy import sparse
 from scipy import linalg
 import scipy.sparse.linalg as splinalg
 import numpy as np
 def compare_density(density, n=1000):
    b = np.random.random(n)
    A1 = sparse.random(n, n, density=density, format='csr')
    A1 += sparse.eye(n,n, format='csr')
    print("Density: ", density)
    print("Sparse: ")
    %timeit splinalg.spsolve(A1, b)
    A2 = A1.todense()
    print("Dense: ")
    %timeit linalg.solve(A2, b)
 # Density=0.001: sparse is faster, Density=0.01: dense is faster
 compare_density(0.01)
 compare_density(0.001)
 ```
 %% Cell type:markdown id: tags:
 ### SciPy - Further Linear Algebra Functions
 * Matrix-matrix and matrix-vector products
    * careful with preservation of sparsity!
 * Solve linear systems
 * Compute eigenvalues and singular values
 * Compute pseudoinverse
 %% Cell type:code id: tags:
 ``` python
 N = 100
 A = np.random.rand(N, N)
 x = np.random.rand(N)
 u, s, vh = linalg.svd(A)        # singular value decomposition
 values, vectors = linalg.eig(A) # eigenvalue problem
 det = linalg.det(A)             # determinant
 P = linalg.pinv(A)              # pseudoinverse (accepts rectangular matrices)
-# Further sparse matrices creation
+# Further sparse matrix operations
 A = sparse.random(N, N, density=0.01, format='csr')
 B = sparse.random(N, N, density=0.01, format='csr')
-C, D = sparse.diags(np.arange(4), 1), sparse.kron(B, A)
+C = sparse.diags(np.arange(4), 1)
+D = sparse.kron(B, A)               # kronecker product of matrices
 ev, evec = splinalg.eigs(A, 2)      # 2 largest eigenvalues
 u, s, vh = sparse.linalg.svds(A, 2) # 2 largest singular values
 ```
 %% Cell type:markdown id: tags:
 ## SciPy Example 5: Interpolation
 * Interpolation of 1D and multi-D data (structured grid and unstructured points)
 * Splines and other polynomials
 %% Cell type:code id: tags:
 ``` python
 import scipy.interpolate as interpolate
 x = np.linspace(0, 10.0, 12)
 y = np.sin(x)
 y_int1 = interpolate.interp1d(x, y)
 y_int2 = interpolate.PchipInterpolator(x, y)
 import matplotlib.pyplot as plt
 plt.figure(figsize=(20,6))
 plt.subplot(1,2,1)
 plt.plot(x, y)
 x2 = np.linspace(0, 10.0, 50)
 plt.plot(x2, y_int1(x2))
 plt.subplot(1,2,2)
 plt.plot(x, y)
 im = plt.plot(x2, y_int2(x2))
 ```
 %% Cell type:markdown id: tags:
 ## NumPy/SciPy Summary
 * NumPy: efficient handling of arrays
 * SciPy: more advanced mathematical and numerical routines; uses NumPy under the hood *plus* other C/F libraries
 * Performance tips:
    * Work on full arrays (slicing, NumPy routines...)
    * identify hotspots and optimize them
        * profile
        * code in Cython or C/Fortran + Python interface -- or try Numba
 %% Cell type:markdown id: tags:
 ## Futher reading
 * Harris, C.R., Millman, K.J., van der Walt, S.J. et al. *Array programming with NumPy.* **Nature** 585, 357–362 (2020).  (https://doi.org/10.1038/s41586-020-2649-2)
 * Bressert, E. (2012). SciPy and NumPy (1st edition.). O'Reilly Media, Inc.  (https://ebooks.mpdl.mpg.de/ebooks/Record/EB001944176)
 * There are numerous books on the topic available: https://ebooks.mpdl.mpg.de/ebooks/Search/Results?type=AllFields&lookfor=numpy
 (MPG.eBooks work from any Max Planck IP address.)
 %% Cell type:markdown id: tags:
 # Numpy Input and Output
 ## Reading and writing NumPy arrays
 Possibilities:
 * NumPy's native IO facilities
 * HDF5 using H5PY
 %% Cell type:markdown id: tags:
 ## Native NumPy Input/Output
 ### Saving NumPy arrays to files
 %% Cell type:code id: tags:
 ``` python
 import os
 import numpy as np
 x = np.linspace(0.0,100.0,101)
 # save to text file (for debugging and small files only)
 np.savetxt("x.txt", x)
 print("savetxt: ", os.path.getsize('x.txt'), "bytes")
 # save to binary file, suffix ".npy"
 np.save("x", x)
 print("save: ", os.path.getsize('x.npy'), "bytes")
 y = np.linspace(0.0,100.0,201)
 # save several arrays to binary file, suffix ".npz", preserve variable names via kwargs
 np.savez("xy", x=x, y=y)
 print("savez: ", os.path.getsize('xy.npz'), "bytes")
 np.savez_compressed("xy_compressed", x=x, y=y)
 print("savez_compressed: ", os.path.getsize('xy_compressed.npz'), "bytes")
 ```
 %% Cell type:markdown id: tags:
 ### Reading NumPy arrays from files
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 #x = np.loadtxt("x.txt")
 x = np.load("x.npy")
 print("type of x:", type(x))
 xy = np.load("xy.npz")
 print("type of xy:", type(xy))
 # individual NumPy arrays in xy are accessible via a dict-style interface
 print("Keys of xy: ", tuple(xy.keys()))
 print("x and xy['x'] are the same: ", np.allclose(x, xy["x"]))
 print("y and xy['y'] are the same: ", np.allclose(y, xy["y"]))
 ```
 %% Cell type:markdown id: tags:
 ## Native NumPy Input/Output
 * Useful for
    * Quickly saving results (e.g. from interactive sessions)
    * Caching results from expensive calculations locally
    * Data that is read only in the same way
 * Note
    * Binary portability (endianness) across systems is supported
 * Problem
    * Not easily accessible from other programs
 %% Cell type:markdown id: tags:
 ## HDF5 using H5PY
 %% Cell type:markdown id: tags:
 ### HDF5 in a nutshell
 * HDF5 is a hierarchichal data format specification and library implementation
 * HDF5 contains data plus metadata, i.e. it is a self-describing format
 * standard and custom datatypes
 * filters (e.g. checksums, compression, scale offset)
 * parallel IO supported, platform independent
 * claimed to be future proof
 * recommended and widely used in HPC by simulation codes, analysis and visualization tools (e.g. VisIt, ParaView, IDL, Matlab)
 (cf. https://support.hdfgroup.org/HDF5/doc/H5.intro.html)
 %% Cell type:markdown id: tags:
 ### HDF5 schematic
 * tree structure with groups, datasets, and attributes, some analogy to a file system with directories and files
 ![HDF5 schematic](fig/hdf5_structure4.jpg)
 %% Cell type:markdown id: tags:
 ## H5PY: HDF5 bindings for Python
 * easy to use Pythonic high-level interface
    * dictionary syntax (dataset access by name)
    * NumPy array syntax (indexing, slicing, etc.)
    * metadata access: use `attrs` attribute
 * designed from the beginning to integrate well with NumPy
 * low-level Cython interface to the HDF5 C API is available
 %% Cell type:markdown id: tags:
 ### Write a NumPy array to HDF5
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 import h5py
 M = np.random.rand(128,128)
 with h5py.File("data.hdf5", 'w') as fp:
    # a HDF5 file always has a root group, "/" which we use implicitly here
    dset = fp.create_dataset("random_matrix", data=M)
    # add metadata (attributes) describing the dataset via the 'attrs' proxy object
    dset.attrs["whatis"] = "This is a randomly generated 128x128 matrix."
    # or even easier
    fp["random_matrix2"] = M
    fp["random_matrix2"].attrs["whatis"] = "This is another randomly generated 128x128 matrix."
 ```
 %% Cell type:markdown id: tags:
 $\to$ use `with` to ensure file is closed and written!
 %% Cell type:markdown id: tags:
 ### Access a HDF5 dataset
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 import h5py
 with h5py.File("data.hdf5", 'r') as fp:
    print("Available keys:", tuple(fp.keys()))
    # retrieve the dataset
    dset = fp["random_matrix"]
    print(type(dset), dset.shape, dset.dtype)
    # retrieve the metadata that comes with the dataset
    for key,val in dset.attrs.items():
        print("{} : {}".format(key, val))
    # dset can be used similarly to a NumPy array
    elem = dset[0,0]  # indexing, slicing
    dsum = np.sum(dset)  # not dset.sum()!
    # explicit conversion to NumPy array
    M = dset[:]
    print(type(M), M.shape, M.dtype)
    msum = M.sum()
 print("Sums are the same:", dsum == msum)
 ```
 %% Cell type:markdown id: tags:
 ### Parallel Access (Preview)
 %% Cell type:code id: tags:
 ``` python
 from mpi4py import MPI
 import h5py
 rank = MPI.COMM_WORLD.rank  # The process ID (integer 0-3 for 4-process run)
 f = h5py.File('parallel_test.hdf5', 'w', driver='mpio', comm=MPI.COMM_WORLD)
 dset = f.create_dataset('test', (4,), dtype='i')
 dset[rank] = rank
 f.close()
 ```
 %% Cell type:markdown id: tags:
 ## Advanced usage of HDF5-IO on HPC systems
 * checkpointing/resuming of long-running jobs
    * write current state of simulation to HDF5 file
    * resume simulation by reading back-in that state
    * upcoming optional exercise: diffusion example
 * large-scale parallel IO from multiple MPI processes
    * write few, large files from each process
      limit the number of files per directory (max. ~1000)
    * write to a single file using MPI-enabled `h5py` (MPCDF environment module `h5py-mpi`)
 %% Cell type:markdown id: tags:
 ## Summary on HDF5 IO of NumPy data
 * easy to use with `h5py` (alternative module would be `pytables`)
 * hierarchical structure enables you to store complex data in a single file
 * HDF5 files are future-proof and portable data files
 * data can be shared & used by other programs
 $\to$ our recommendation for NumPy-related I/O!
 Further reading:
 * Collette, A. (2013). Python and HDF5 (1st edition.). O'Reilly Media, Inc.
  (https://ebooks.mpdl.mpg.de/ebooks/Record/EB001941719)
 %% Cell type:markdown id: tags:
 # Application: image processing
 %% Cell type:markdown id: tags:
 ## Example 1: Smooth Images
 ![image filter](fig/image-filter-example.svg)
 %% Cell type:markdown id: tags:
 ## Smoothing with image filters
 * Goal: subtract noise
 * Use filtering approach
 * For each pixel, at $i, j$, compute:
 $$\tilde{f}_{i,j} = \alpha f_{i,j} + \frac{(1 - \alpha)}{4} (f_{i-1,j} + f_{i+1,j} + f_{i,j-1} + f_{i,j+1})$$
 * Mix pixel with average of surrounding pixels $\to$ *filter out high-frequency noise*
 * Can also be rewritten using the discrete Laplace operator:
 $$\tilde{f}_{i,j} = f_{i,j} + \frac{(1 - \alpha)}{4} \Delta f_{i,j}$$
 * Strongest smoothing for $\alpha = 0$, no smoothing for $\alpha = 1$
 %% Cell type:markdown id: tags:
 ### Image: alpha Centauri
 * Image from DSS (for copyright see http://archive.stsci.edu/dss/acknowledging.html)
 %% Cell type:code id: tags:
 ``` python
 import numpy as np
 import matplotlib.pyplot as plt
 # convert to float between 0 and 1
 image = plt.imread("fig/dss_alpha_centauri.gif").astype(np.float64)/255.
 plt.imshow(image, cmap='gray')
 print("Shape: ", image.shape)
 ```
 %% Cell type:markdown id: tags:
 ## Loops in python vs. numpy arrays
 %% Cell type:code id: tags:
 ``` python
 def get_filtered_loop(image, alpha=0.25):
    filtered = np.zeros_like(image)
    N, M = image.shape
    for i in range(1, N-1):
        for j in range(1, M-1):
            filtered[i, j] = 0.25*(1.0-alpha)*(image[i-1, j] + image[i+1, j]
                                               + image[i, j-1] + image[i, j+1]) \
                + alpha*image[i, j]
    return filtered
 %timeit get_filtered_loop(image)
 ```
 %% Cell type:code id: tags:
 ``` python
 def get_filtered_array(image, alpha=0.25):
    filtered = np.zeros_like(image)
    N, M = image.shape
    filtered[1:N-1, 1:M-1] = 0.25*(1.0-alpha)*(image[0:N-2, 1:M-1] + image[2:N, 1:M-1]
                                               + image[1:N-1, 0:M-2] + image[1:N-1, 2:M]) \
        + alpha*image[1:N-1, 1:M-1]
    return filtered
 %timeit get_filtered_array(image)
 ```
 %% Cell type:markdown id: tags:
 $\to$ roughly factor 100-150 faster!
 %% Cell type:markdown id: tags:
 ## Transformation step by step
 * Transform for loops to indices on left hand side (same range)
 * From:
 ```python
 for i in range(1, N-1):
    for j in range(1, M-1):
 ```
 To:
 ```python
 filtered[1:N-1, 1:M-1]
 ```
 * Shift indices on right hand side: `i-1` $\to$ `0:N-2`, `i` $\to$ `1:N-1`, `i+1` $\to$ `2:N`
 * E.g.: `image[i-1, j]` $\to$ `image[0:N-2, 1:M-1]`
 * Similar for `j`
 %% Cell type:code id: tags:
 ``` python
 filtered_loop = get_filtered_loop(image)
 filtered_array = get_filtered_array(image)
 f, axs = plt.subplots(2, 3, figsize=(14, 8))
 axs[0, 0].imshow(image, cmap='gray')
 axs[0, 1].imshow(filtered_array, cmap='gray')
 axs[0, 2].imshow(filtered_array - filtered_loop, cmap='gray')
 axs[1, 0].imshow(image[600:700, 600:700], cmap='gray')
 axs[1, 1].imshow(filtered_array[600:700, 600:700], cmap='gray')
 axs[1, 2].imshow((image-filtered_array)[600:700, 600:700], cmap='gray')
 print("Loop and array method yields the same: ", np.allclose(filtered_loop, filtered_array))
 ```
 %% Cell type:markdown id: tags:
 ## Example 2: Get edges
 * Use the following formula:
 $$\tilde{f}_{i,j} = 8 f_{i,j} - (f_{i-1,j} + f_{i+1,j} + f_{i,j-1} + f_{i,j+1} + f_{i-1,j-1} + f_{i-1,j+1} + f_{i+1,j+1} + f_{i+1,j-1})$$
 * Subtract mean of surrounding cells
 * Also use diagonals
 * Truncate results to a certain percentile
 %% Cell type:code id: tags:
 ``` python
 from scipy.misc import ascent
 im = plt.imshow(ascent(), cmap='gray')
 ```
 %% Cell type:markdown id: tags:
 ### Our implementation
 %% Cell type:code id: tags:
 ``` python
 def get_edges(image, percentiles=(25, 99)):
    result = np.zeros_like(image)
    N, M = image.shape
    result[1:N-1, 1:M-1] = -1.0*(image[0:N-2, 1:M-1] + image[2:N, 1:M-1]
                                 + image[1:N-1, 0:M-2] + image[1:N-1, 2:M]
                                 + image[0:N-2, 0:M-2] + image[0:N-2, 2:M]
                                 + image[2:N, 0:M-2] + image[2:N, 2:M])     + 8.*image[1:N-1, 1:M-1]
    # now truncate results
    fmin, fmax = np.percentile(result, percentiles)
    result[result < fmin] = fmin
    result[result > fmax] = fmax
    return result
 f, axs = plt.subplots(1, 2, figsize=(15, 5))
 axs[0].imshow(ascent(), cmap='gray')
 axs[1].imshow(get_edges(ascent()), cmap='gray');None
 ```
 %% Cell type:markdown id: tags:
 ### Predefined filter
 %% Cell type:code id: tags:
 ``` python
 from PIL import Image, ImageFilter
 img = Image.fromarray(ascent().astype("int32"), mode="I").convert("L")
 filtered_img = img.filter(ImageFilter.FIND_EDGES)
 f, axs = plt.subplots(1, 2, figsize=(15, 5))
 axs[0].imshow(ascent(), cmap='gray')
 axs[1].imshow(np.asarray(filtered_img), cmap='gray');None
 ```
 %% Cell type:markdown id: tags:
 $\to$ results very similar!