Python scipy.linalg.solve_banded() Examples

def _radau(alpha, beta, xr):
    """From <http://www.scientificpython.net/pyblog/radau-quadrature>:
    Compute the Radau nodes and weights with the preassigned node xr.

    Based on the section 7 of the paper

        Some modified matrix eigenvalue problems,
        Gene Golub,
        SIAM Review Vol 15, No. 2, April 1973, pp.318--334.
    """
    from scipy.linalg import solve_banded

    n = len(alpha) - 1
    f = numpy.zeros(n)
    f[-1] = beta[-1]
    A = numpy.vstack((numpy.sqrt(beta), alpha - xr))
    J = numpy.vstack((A[:, 0:-1], A[0, 1:]))
    delta = solve_banded((1, 1), J, f)
    alphar = alpha.copy()
    alphar[-1] = xr + delta[-1]
    x, w = scheme_from_rc(alphar, beta, mode="numpy")
    return x, w 

def _get_natural_f(knots):
    """Returns mapping of natural cubic spline values to 2nd derivatives.

    .. note:: See 'Generalized Additive Models', Simon N. Wood, 2006, pp 145-146

    :param knots: The 1-d array knots used for cubic spline parametrization,
     must be sorted in ascending order.
    :return: A 2-d array mapping natural cubic spline values at
     knots to second derivatives.

    :raise ImportError: if scipy is not found, required for
     ``linalg.solve_banded()``
    """
    try:
        from scipy import linalg
    except ImportError: # pragma: no cover
        raise ImportError("Cubic spline functionality requires scipy.")

    h = knots[1:] - knots[:-1]
    diag = (h[:-1] + h[1:]) / 3.
    ul_diag = h[1:-1] / 6.
    banded_b = np.array([np.r_[0., ul_diag], diag, np.r_[ul_diag, 0.]])
    d = np.zeros((knots.size - 2, knots.size))
    for i in range(knots.size - 2):
        d[i, i] = 1. / h[i]
        d[i, i + 2] = 1. / h[i + 1]
        d[i, i + 1] = - d[i, i] - d[i, i + 2]

    fm = linalg.solve_banded((1, 1), banded_b, d)

    return np.vstack([np.zeros(knots.size), fm, np.zeros(knots.size)])


# Cyclic Cubic Regression Splines 

def __interpolateSpline__(self):
        """Calculate natural spline coefficients that fit the data exactly."""
        # Get x and y values for all supplied points
        x = self.source[:, 0]
        y = self.source[:, 1]
        mdim = len(x)
        h = [x[i + 1] - x[i] for i in range(0, mdim - 1)]
        # Initialize the matrix
        Ab = np.zeros((3, mdim))
        # Construct the Ab banded matrix and B vector
        Ab[1, 0] = 1  # A[0, 0] = 1
        B = [0]
        for i in range(1, mdim - 1):
            Ab[2, i - 1] = h[i - 1]  # A[i, i - 1] = h[i - 1]
            Ab[1, i] = 2 * (h[i] + h[i - 1])  # A[i, i] = 2*(h[i] + h[i - 1])
            Ab[0, i + 1] = h[i]  # A[i, i + 1] = h[i]
            B.append(3 * ((y[i + 1] - y[i]) / (h[i]) - (y[i] - y[i - 1]) / (h[i - 1])))
        Ab[1, mdim - 1] = 1  # A[-1, -1] = 1
        B.append(0)
        # Solve the system for c coefficients
        c = linalg.solve_banded((1, 1), Ab, B, True, True)
        # Calculate other coefficients
        b = [
            ((y[i + 1] - y[i]) / h[i] - h[i] * (2 * c[i] + c[i + 1]) / 3)
            for i in range(0, mdim - 1)
        ]
        d = [(c[i + 1] - c[i]) / (3 * h[i]) for i in range(0, mdim - 1)]
        # Store coefficients
        self.__splineCoefficients__ = np.array([y[0:-1], b, c[0:-1], d]) 

def _lobatto(alpha, beta, xl1, xl2):
    """Compute the Lobatto nodes and weights with the preassigned node xl1, xl2.
    Based on the section 7 of the paper

        Some modified matrix eigenvalue problems,
        Gene Golub,
        SIAM Review Vol 15, No. 2, April 1973, pp.318--334,

    and

        http://www.scientificpython.net/pyblog/radau-quadrature
    """
    from scipy.linalg import solve_banded, solve

    n = len(alpha) - 1
    en = numpy.zeros(n)
    en[-1] = 1
    A1 = numpy.vstack((numpy.sqrt(beta), alpha - xl1))
    J1 = numpy.vstack((A1[:, 0:-1], A1[0, 1:]))
    A2 = numpy.vstack((numpy.sqrt(beta), alpha - xl2))
    J2 = numpy.vstack((A2[:, 0:-1], A2[0, 1:]))
    g1 = solve_banded((1, 1), J1, en)
    g2 = solve_banded((1, 1), J2, en)
    C = numpy.array(((1, -g1[-1]), (1, -g2[-1])))
    xl = numpy.array((xl1, xl2))
    ab = solve(C, xl)

    alphal = alpha
    alphal[-1] = ab[0]
    betal = beta
    betal[-1] = ab[1]
    x, w = scheme_from_rc(alphal, betal, mode="numpy")
    return x, w 

def SolveBanded(A, D, bw=3):
    # Find the diagonals
    ud = np.insert(np.diag(A,1), 0, 0) # upper diagonal
    d = np.diag(A) # main diagonal
    ld = np.insert(np.diag(A,-1), len(d)-1, 0) # lower diagonal
    # simplified matrix
    ab = np.matrix([
        ud,
        d,
        ld,
    ])
    #ds9(ab)
    return solve_banded((1, 1), ab, D ) 

def _get_natural_f(knots):
    """Returns mapping of natural cubic spline values to 2nd derivatives.

    .. note:: See 'Generalized Additive Models', Simon N. Wood, 2006, pp 145-146

    :param knots: The 1-d array knots used for cubic spline parametrization,
     must be sorted in ascending order.
    :return: A 2-d array mapping natural cubic spline values at
     knots to second derivatives.

    :raise ImportError: if scipy is not found, required for
     ``linalg.solve_banded()``
    """
    try:
        from scipy import linalg
    except ImportError: # pragma: no cover
        raise ImportError("Cubic spline functionality requires scipy.")

    h = knots[1:] - knots[:-1]
    diag = (h[:-1] + h[1:]) / 3.
    ul_diag = h[1:-1] / 6.
    banded_b = np.array([np.r_[0., ul_diag], diag, np.r_[ul_diag, 0.]])
    d = np.zeros((knots.size - 2, knots.size))
    for i in range(knots.size - 2):
        d[i, i] = 1. / h[i]
        d[i, i + 2] = 1. / h[i + 1]
        d[i, i + 1] = - d[i, i] - d[i, i + 2]

    fm = linalg.solve_banded((1, 1), banded_b, d)

    return np.vstack([np.zeros(knots.size), fm, np.zeros(knots.size)])


# Cyclic Cubic Regression Splines 

def __init__(self, x, y):
        x = np.asarray(x, dtype=float)
        y = np.asarray(y, dtype=float)

        n = x.shape[0]
        dx = np.diff(x)
        dy = np.diff(y, axis=0)
        dxr = dx.reshape([dx.shape[0]] + [1] * (y.ndim - 1))

        c = np.empty((3, n - 1) + y.shape[1:])
        if n > 2:
            A = np.ones((2, n))
            b = np.empty((n,) + y.shape[1:])
            b[0] = 0
            b[1:] = 2 * dy / dxr
            s = solve_banded((1, 0), A, b, overwrite_ab=True, overwrite_b=True,
                             check_finite=False)
            c[0] = np.diff(s, axis=0) / (2 * dxr)
            c[1] = s[:-1]
            c[2] = y[:-1]
        else:
            c[0] = 0
            c[1] = dy / dxr
            c[2] = y[:-1]

        super(_QuadraticSpline, self).__init__(c, x) 

def _lobatto(coefficients, preassigned):
    """
    Compute the Lobatto nodes and weights with the preassigned value pair.
    Based on the section 7 of the paper

        Some modified matrix eigenvalue problems,
        Gene Golub,
        SIAM Review Vol 15, No. 2, April 1973, pp.318--334,

    and

        http://www.scientificpython.net/pyblog/radau-quadrature

    Args:
        coefficients (numpy.ndarray):
            Three terms recurrence coefficients.
        preassigned (Tuple[float, float]):
            Values that are assume to be fixed.
    """
    alpha = numpy.array(coefficients[0])
    beta = numpy.array(coefficients[1])
    vec_en = numpy.zeros(len(alpha)-1)
    vec_en[-1] = 1
    mat_a1 = numpy.vstack((numpy.sqrt(beta), alpha-preassigned[0]))
    mat_j1 = numpy.vstack((mat_a1[:, 0:-1], mat_a1[0, 1:]))
    mat_a2 = numpy.vstack((numpy.sqrt(beta), alpha - preassigned[1]))
    mat_j2 = numpy.vstack((mat_a2[:, 0:-1], mat_a2[0, 1:]))
    mat_g1 = solve_banded((1, 1), mat_j1, vec_en)
    mat_g2 = solve_banded((1, 1), mat_j2, vec_en)
    mat_c = numpy.array(((1, -mat_g1[-1]), (1, -mat_g2[-1])))
    alpha[-1], beta[-1] = solve(mat_c, preassigned)

    return numpy.array([alpha, beta])