difftools.py

# difftools.py - Several useful utilities for numerical differentiation
# 
# Author: Stefan Fuertinger [stefan.fuertinger@gmx.at]
# Created: June 13 2012
# Last modified: <2017-09-21 12:08:48>

from __future__ import division
import numpy as np
from scipy.sparse import spdiags, eye, kron 

##########################################################################################
def fidop2d(N, drt='xy', fds='c'):
    """
    Compute 2D finite difference operators

    Parameters
    ----------
    N : int
        Positive integer holding the grid dimension (has to be square, i.e.
        `N`-by-`N`!)
    drt : str
        String determining the direction of the discrete derivatives. 
        Can be either 'x', 'y' or 'xy'. 
    fds : str
        String determining the employed finite difference scheme. Can be 
        either 'c' for centered, 'f' for forward or 'b' for backward. 
       
    Returns
    -------
    D : SciPy sparse matrix/matrices
        Depending on the given input one or two sparse matrices 
        corresponding to the wanted discrete derivative operators are 
        returned. 

    Examples
    --------
    The command 

    >>> Dx,Dy = fidop2d(N) 

    returns the sparse `N**2`-by-`N**2` matrices `Dx` and `Dy` such that `Dh` defined by

    >>> import numpy as np
    >>> Dh = np.hstack([Dx,Dy])

    is the discrete gradient operator in lexicographic 
    ordering for a function `F` given on a grid of size `N`-by-`N`.
    The matrix `Dx` is a discrete approximation of the first-order derivative operator 
    in `x`-direction, Analogously, `Dy` discretizes the first order derivative operator 
    in `y`-direction. The spacing between points in each direction is assumed to be one 
    (i.e. the step size `h = 1`). 
    
    The command 

    >>> D = fidop2d(N,drt) 

    returns the sparse `N**2`-by-`N**2` matrix `D` corresponding to the discrete derivative 
    operator in direction `drt` in lexicographic ordering. The string `drt` has to be either 'x' (default)
    , 'y' or 'xy'. Note: only 'xy' will return two matrices, namely `Dx`, `Dy`. 

    The command 

    >>> D = fidop2d(N,drt,fds) 

    returns sparse the `N**2`-by-`N**2` matrix `D` corresponding to the discrete derivative 
    operator in direction `drt` in lexicographic ordering using the finite difference scheme `fds`. 
    The string `fds` can either be 'c' (centered differences, in cells, default), 
    'f' (forward differences, in interfaces) or 'b' (backward differences, in interfaces). 

    To discretize first order differential operators on a 4-by-4 grid in both `x`- and `y`-directions 
    using forward differences use

    >>> Dx,Dy = fidop2d(4,fds='f')

    To get only the discrete first order differential operator in `y`-direction on a grid of size `N` use

    >>> Dy = fidop2d(N,'y')

    Notes
    -----
    For clarity here some comments on the rationale behind the code. 

    *The Laplacian*

    If `Dxf` and `Dyf` are computed using forward differences, i.e.

    >>> Dxf,Dyf = fidop2d(N,'xy','f')

    then the discrete Laplacian `Lh` based on centered differences is
    obtained by 

    >>> Lh = -(Dxf.T*Dxf + Dyf.T*Dyf)

    or analogously for backward differences, 

    >>> Dxb,Dyb = fidop2d(N,'xy','b')

    then

    >>> Lh = Dxb.T*Dxb + Dyb.T*Dyb
    
    *The Divergence*

    Note adjoints: for :math:`p` in 
    :math:`H_0(\\mathrm{div}):=\\{p \\in L^2(\\Omega)` | div :math:`(p) \\in L^2(\\Omega), pn = 0` 
    on the boundary of :math:`\\Omega\\}` we have

    .. math:: 

       \\begin{align}
       \\int_{\\Omega} \\nabla u\\cdot p dx &=  \\int_{\\Omega} u p \\cdot  n dS - \\int_{\\Omega} u \\textrm{ div }(p) dx\\\\
                                     & = 0 -\\int_{\\Omega} u \\textrm{ div }(p) dx
       \\end{align}

    or in short :math:`\\textrm{div}^\\ast = -\\nabla`. Hence for a spatial discretization take 
    `Div_h ~ -Dh.T` as approximation for the divergence. 
    """

    # Sanity checks for `N`
    if not np.isscalar(N) or not plt.is_numlike(N) or not np.isreal(N).all():
        raise TypeError("Input `N` must be a real scalar!")
    if not np.isfinite(N):
        raise TypeError("Input `N` must be finite!")
    if round(N) != N:
        raise TypeError("`N` has to be a positive integer!")
    if N <= 1:
        raise ValueError("`N` has to be greater than 1!")
    
    # Sanity checks for `drt`
    if not isinstance(drt,(str,unicode)):
        raise TypeError("Input `drt` has to be a string!")
    if drt != 'x' and drt != 'y' and drt != 'xy':
        raise ValueError("Input `drt` has to be either 'x', 'y' or 'xy'")
    elif drt == 'x' or drt == 'y':
        myout = 1
    else:
        myout = 2

    # Sanity checks for `fds`
    if not isinstance(fds,(str,unicode)):
        raise TypeError("Input `fds` has to be a string!")
    if fds != 'c' and fds != 'b' and fds != 'f':
        raise ValueError("Input `fds` has to be either 'c' (centered), 'b' (backward) or 'f' (forward)")

    # Initialize vector of ones needed to build matrices
    e = np.ones((1,N));

    # Building blocks for the operators
    if fds=='c':
        # Centered differences
        A = spdiags(np.r_[-e,e],np.array([-1,1]),N,N,format='lil');
        A[0,0] = -2.0; A[0,1] = 2.0; A[N-1,N-2] = -2.0; A[N-1,N-1] = 2.0;
        A = 0.5*A;
        A = A.tocsr()
    elif fds=='f':
        # Forward differences
        A = spdiags(np.r_[-e,e],np.array([0,1]),N,N,format='lil');
        A[N-1,:] = 0.0;
        A = A.tocsr()
    else:
        # Backward differences
        A = spdiags(np.r_[-e,e],np.array([-1,0]),N,N,format='lil');
        A[0,:] = 0.0;
        A = A.tocsr()

    # Second building block needed: the identity
    IN = eye(N,N,format='csr');

    # Compute output matrix/matrices
    if myout==1:
        # Only one output: determine direction
        if drt=='x':
            # x-direction
            return kron(IN,A,format='csr')
        else:
            # y-direction
            return kron(A,IN,format='csr')
    else:
        # Two outputs: order of directions is always d/dx,d/dy]
        return kron(IN,A,format='csr'), kron(A,IN,format='csr')

##########################################################################################
def myff2n(n):
    """
    Give factor settings for a two-level full factorial design

    Parameters
    ----------
    n : int
        Number of factors
       
    Returns
    -------
    dff2 : NumPy 2darray
        An `m`-by-`n` array of factor settings 

    Examples    
    --------
    >>> dFF2 = ff2n(3)
    >>> dFF2
        array([[ 0.,  0.,  0.],
               [ 0.,  0.,  1.],
               [ 0.,  1.,  0.],
               [ 0.,  1.,  1.],
               [ 1.,  0.,  0.],
               [ 1.,  0.,  1.],
               [ 1.,  1.,  0.],
               [ 1.,  1.,  1.]])
    """

    # Check correctness of input
    if not np.isscalar(n) or not plt.is_numlike(n) or not np.isreal(n).all():
        raise TypeError("Input `n` must be a real scalar!")
    if not np.isfinite(n):
        raise TypeError("Input `n` must be finite!")
    if round(n) != n:
        raise TypeError("`n` has to be a positive integer!")
    if n <= 0:
        raise ValueError("`n` has to be greater than 0!")

    # Output array dff2 has 2^n rows 
    rows = 2**n
    ncycles = rows
    dff2 = np.zeros((rows,n))

    # This is adapted from the MATLAB source file ff2n.m
    for k in xrange(0,n):
        settings = np.arange(0,2)
        settings.shape = (1,2)
        ncycles = ncycles/2.0
        nreps = rows/(2*ncycles)
        settings = settings[np.zeros((1,nreps)).astype(int),:]
        settings = settings.flatten(1)
        settings.shape = (settings.size,1)
        settings = settings[:,np.zeros((1,ncycles)).astype(int)]
        dff2[:,n-k-1] = settings.flatten(1)

    return dff2