"""Routines for numerical differentiation."""

from __future__ import division

import numpy as np
from numpy.linalg import norm

from scipy.sparse.linalg import LinearOperator
from ..sparse import issparse, csc_matrix, csr_matrix, coo_matrix, find
from ._group_columns import group_dense, group_sparse

EPS = np.finfo(np.float64).eps


def _adjust_scheme_to_bounds(x0, h, num_steps, scheme, lb, ub):
    """Adjust final difference scheme to the presence of bounds.

    Parameters
    ----------
    x0 : ndarray, shape (n,)
        Point at which we wish to estimate derivative.
    h : ndarray, shape (n,)
        Desired finite difference steps.
    num_steps : int
        Number of `h` steps in one direction required to implement finite
        difference scheme. For example, 2 means that we need to evaluate
        f(x0 + 2 * h) or f(x0 - 2 * h)
    scheme : {'1-sided', '2-sided'}
        Whether steps in one or both directions are required. In other
        words '1-sided' applies to forward and backward schemes, '2-sided'
        applies to center schemes.
    lb : ndarray, shape (n,)
        Lower bounds on independent variables.
    ub : ndarray, shape (n,)
        Upper bounds on independent variables.

    Returns
    -------
    h_adjusted : ndarray, shape (n,)
        Adjusted step sizes. Step size decreases only if a sign flip or
        switching to one-sided scheme doesn't allow to take a full step.
    use_one_sided : ndarray of bool, shape (n,)
        Whether to switch to one-sided scheme. Informative only for
        ``scheme='2-sided'``.
    """
    if scheme == '1-sided':
        use_one_sided = np.ones_like(h, dtype=bool)
    elif scheme == '2-sided':
        h = np.abs(h)
        use_one_sided = np.zeros_like(h, dtype=bool)
    else:
        raise ValueError("`scheme` must be '1-sided' or '2-sided'.")

    if np.all((lb == -np.inf) & (ub == np.inf)):
        return h, use_one_sided

    h_total = h * num_steps
    h_adjusted = h.copy()

    lower_dist = x0 - lb
    upper_dist = ub - x0

    if scheme == '1-sided':
        x = x0 + h_total
        violated = (x < lb) | (x > ub)
        fitting = np.abs(h_total) <= np.maximum(lower_dist, upper_dist)
        h_adjusted[violated & fitting] *= -1

        forward = (upper_dist >= lower_dist) & ~fitting
        h_adjusted[forward] = upper_dist[forward] / num_steps
        backward = (upper_dist < lower_dist) & ~fitting
        h_adjusted[backward] = -lower_dist[backward] / num_steps
    elif scheme == '2-sided':
        central = (lower_dist >= h_total) & (upper_dist >= h_total)

        forward = (upper_dist >= lower_dist) & ~central
        h_adjusted[forward] = np.minimum(
            h[forward], 0.5 * upper_dist[forward] / num_steps)
        use_one_sided[forward] = True

        backward = (upper_dist < lower_dist) & ~central
        h_adjusted[backward] = -np.minimum(
            h[backward], 0.5 * lower_dist[backward] / num_steps)
        use_one_sided[backward] = True

        min_dist = np.minimum(upper_dist, lower_dist) / num_steps
        adjusted_central = (~central & (np.abs(h_adjusted) <= min_dist))
        h_adjusted[adjusted_central] = min_dist[adjusted_central]
        use_one_sided[adjusted_central] = False

    return h_adjusted, use_one_sided


relative_step = {"2-point": EPS**0.5,
                 "3-point": EPS**(1/3),
                 "cs": EPS**0.5}


def _compute_absolute_step(rel_step, x0, method):
    if rel_step is None:
        rel_step = relative_step[method]
    sign_x0 = (x0 >= 0).astype(float) * 2 - 1
    return rel_step * sign_x0 * np.maximum(1.0, np.abs(x0))


def _prepare_bounds(bounds, x0):
    lb, ub = [np.asarray(b, dtype=float) for b in bounds]
    if lb.ndim == 0:
        lb = np.resize(lb, x0.shape)

    if ub.ndim == 0:
        ub = np.resize(ub, x0.shape)

    return lb, ub


def group_columns(A, order=0):
    """Group columns of a 2-D matrix for sparse finite differencing [1]_.

    Two columns are in the same group if in each row at least one of them
    has zero. A greedy sequential algorithm is used to construct groups.

    Parameters
    ----------
    A : array_like or sparse matrix, shape (m, n)
        Matrix of which to group columns.
    order : int, iterable of int with shape (n,) or None
        Permutation array which defines the order of columns enumeration.
        If int or None, a random permutation is used with `order` used as
        a random seed. Default is 0, that is use a random permutation but
        guarantee repeatability.

    Returns
    -------
    groups : ndarray of int, shape (n,)
        Contains values from 0 to n_groups-1, where n_groups is the number
        of found groups. Each value ``groups[i]`` is an index of a group to
        which ith column assigned. The procedure was helpful only if
        n_groups is significantly less than n.

    References
    ----------
    .. [1] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
           sparse Jacobian matrices", Journal of the Institute of Mathematics
           and its Applications, 13 (1974), pp. 117-120.
    """
    if issparse(A):
        A = csc_matrix(A)
    else:
        A = np.atleast_2d(A)
        A = (A != 0).astype(np.int32)

    if A.ndim != 2:
        raise ValueError("`A` must be 2-dimensional.")

    m, n = A.shape

    if order is None or np.isscalar(order):
        rng = np.random.RandomState(order)
        order = rng.permutation(n)
    else:
        order = np.asarray(order)
        if order.shape != (n,):
            raise ValueError("`order` has incorrect shape.")

    A = A[:, order]

    if issparse(A):
        groups = group_sparse(m, n, A.indices, A.indptr)
    else:
        groups = group_dense(m, n, A)

    groups[order] = groups.copy()

    return groups


def approx_derivative(fun, x0, method='3-point', rel_step=None, f0=None,
                      bounds=(-np.inf, np.inf), sparsity=None,
                      as_linear_operator=False, args=(), kwargs={}):
    """Compute finite difference approximation of the derivatives of a
    vector-valued function.

    If a function maps from R^n to R^m, its derivatives form m-by-n matrix
    called the Jacobian, where an element (i, j) is a partial derivative of
    f[i] with respect to x[j].

    Parameters
    ----------
    fun : callable
        Function of which to estimate the derivatives. The argument x
        passed to this function is ndarray of shape (n,) (never a scalar
        even if n=1). It must return 1-D array_like of shape (m,) or a scalar.
    x0 : array_like of shape (n,) or float
        Point at which to estimate the derivatives. Float will be converted
        to a 1-D array.
    method : {'3-point', '2-point', 'cs'}, optional
        Finite difference method to use:
            - '2-point' - use the first order accuracy forward or backward
                          difference.
            - '3-point' - use central difference in interior points and the
                          second order accuracy forward or backward difference
                          near the boundary.
            - 'cs' - use a complex-step finite difference scheme. This assumes
                     that the user function is real-valued and can be
                     analytically continued to the complex plane. Otherwise,
                     produces bogus results.
    rel_step : None or array_like, optional
        Relative step size to use. The absolute step size is computed as
        ``h = rel_step * sign(x0) * max(1, abs(x0))``, possibly adjusted to
        fit into the bounds. For ``method='3-point'`` the sign of `h` is
        ignored. If None (default) then step is selected automatically,
        see Notes.
    f0 : None or array_like, optional
        If not None it is assumed to be equal to ``fun(x0)``, in  this case
        the ``fun(x0)`` is not called. Default is None.
    bounds : tuple of array_like, optional
        Lower and upper bounds on independent variables. Defaults to no bounds.
        Each bound must match the size of `x0` or be a scalar, in the latter
        case the bound will be the same for all variables. Use it to limit the
        range of function evaluation. Bounds checking is not implemented
        when `as_linear_operator` is True.
    sparsity : {None, array_like, sparse matrix, 2-tuple}, optional
        Defines a sparsity structure of the Jacobian matrix. If the Jacobian
        matrix is known to have only few non-zero elements in each row, then
        it's possible to estimate its several columns by a single function
        evaluation [3]_. To perform such economic computations two ingredients
        are required:

        * structure : array_like or sparse matrix of shape (m, n). A zero
          element means that a corresponding element of the Jacobian
          identically equals to zero.
        * groups : array_like of shape (n,). A column grouping for a given
          sparsity structure, use `group_columns` to obtain it.

        A single array or a sparse matrix is interpreted as a sparsity
        structure, and groups are computed inside the function. A tuple is
        interpreted as (structure, groups). If None (default), a standard
        dense differencing will be used.

        Note, that sparse differencing makes sense only for large Jacobian
        matrices where each row contains few non-zero elements.
    as_linear_operator : bool, optional
        When True the function returns an `scipy.sparse.linalg.LinearOperator`.
        Otherwise it returns a dense array or a sparse matrix depending on
        `sparsity`. The linear operator provides an efficient way of computing
        ``J.dot(p)`` for any vector ``p`` of shape (n,), but does not allow
        direct access to individual elements of the matrix. By default
        `as_linear_operator` is False.
    args, kwargs : tuple and dict, optional
        Additional arguments passed to `fun`. Both empty by default.
        The calling signature is ``fun(x, *args, **kwargs)``.

    Returns
    -------
    J : {ndarray, sparse matrix, LinearOperator}
        Finite difference approximation of the Jacobian matrix.
        If `as_linear_operator` is True returns a LinearOperator
        with shape (m, n). Otherwise it returns a dense array or sparse
        matrix depending on how `sparsity` is defined. If `sparsity`
        is None then a ndarray with shape (m, n) is returned. If
        `sparsity` is not None returns a csr_matrix with shape (m, n).
        For sparse matrices and linear operators it is always returned as
        a 2-D structure, for ndarrays, if m=1 it is returned
        as a 1-D gradient array with shape (n,).

    See Also
    --------
    check_derivative : Check correctness of a function computing derivatives.

    Notes
    -----
    If `rel_step` is not provided, it assigned to ``EPS**(1/s)``, where EPS is
    machine epsilon for float64 numbers, s=2 for '2-point' method and s=3 for
    '3-point' method. Such relative step approximately minimizes a sum of
    truncation and round-off errors, see [1]_.

    A finite difference scheme for '3-point' method is selected automatically.
    The well-known central difference scheme is used for points sufficiently
    far from the boundary, and 3-point forward or backward scheme is used for
    points near the boundary. Both schemes have the second-order accuracy in
    terms of Taylor expansion. Refer to [2]_ for the formulas of 3-point
    forward and backward difference schemes.

    For dense differencing when m=1 Jacobian is returned with a shape (n,),
    on the other hand when n=1 Jacobian is returned with a shape (m, 1).
    Our motivation is the following: a) It handles a case of gradient
    computation (m=1) in a conventional way. b) It clearly separates these two
    different cases. b) In all cases np.atleast_2d can be called to get 2-D
    Jacobian with correct dimensions.

    References
    ----------
    .. [1] W. H. Press et. al. "Numerical Recipes. The Art of Scientific
           Computing. 3rd edition", sec. 5.7.

    .. [2] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
           sparse Jacobian matrices", Journal of the Institute of Mathematics
           and its Applications, 13 (1974), pp. 117-120.

    .. [3] B. Fornberg, "Generation of Finite Difference Formulas on
           Arbitrarily Spaced Grids", Mathematics of Computation 51, 1988.

    Examples
    --------
    >>> import numpy as np
    >>> from scipy.optimize import approx_derivative
    >>>
    >>> def f(x, c1, c2):
    ...     return np.array([x[0] * np.sin(c1 * x[1]),
    ...                      x[0] * np.cos(c2 * x[1])])
    ...
    >>> x0 = np.array([1.0, 0.5 * np.pi])
    >>> approx_derivative(f, x0, args=(1, 2))
    array([[ 1.,  0.],
           [-1.,  0.]])

    Bounds can be used to limit the region of function evaluation.
    In the example below we compute left and right derivative at point 1.0.

    >>> def g(x):
    ...     return x**2 if x >= 1 else x
    ...
    >>> x0 = 1.0
    >>> approx_derivative(g, x0, bounds=(-np.inf, 1.0))
    array([ 1.])
    >>> approx_derivative(g, x0, bounds=(1.0, np.inf))
    array([ 2.])
    """
    if method not in ['2-point', '3-point', 'cs']:
        raise ValueError("Unknown method '%s'. " % method)

    x0 = np.atleast_1d(x0)
    if x0.ndim > 1:
        raise ValueError("`x0` must have at most 1 dimension.")

    lb, ub = _prepare_bounds(bounds, x0)

    if lb.shape != x0.shape or ub.shape != x0.shape:
        raise ValueError("Inconsistent shapes between bounds and `x0`.")

    if as_linear_operator and not (np.all(np.isinf(lb))
                                   and np.all(np.isinf(ub))):
        raise ValueError("Bounds not supported when "
                         "`as_linear_operator` is True.")

    def fun_wrapped(x):
        f = np.atleast_1d(fun(x, *args, **kwargs))
        if f.ndim > 1:
            raise RuntimeError("`fun` return value has "
                               "more than 1 dimension.")
        return f

    if f0 is None:
        f0 = fun_wrapped(x0)
    else:
        f0 = np.atleast_1d(f0)
        if f0.ndim > 1:
            raise ValueError("`f0` passed has more than 1 dimension.")

    if np.any((x0 < lb) | (x0 > ub)):
        raise ValueError("`x0` violates bound constraints.")

    if as_linear_operator:
        if rel_step is None:
            rel_step = relative_step[method]

        return _linear_operator_difference(fun_wrapped, x0,
                                           f0, rel_step, method)
    else:
        h = _compute_absolute_step(rel_step, x0, method)

        if method == '2-point':
            h, use_one_sided = _adjust_scheme_to_bounds(
                x0, h, 1, '1-sided', lb, ub)
        elif method == '3-point':
            h, use_one_sided = _adjust_scheme_to_bounds(
                x0, h, 1, '2-sided', lb, ub)
        elif method == 'cs':
            use_one_sided = False

        if sparsity is None:
            return _dense_difference(fun_wrapped, x0, f0, h,
                                     use_one_sided, method)
        else:
            if not issparse(sparsity) and len(sparsity) == 2:
                structure, groups = sparsity
            else:
                structure = sparsity
                groups = group_columns(sparsity)

            if issparse(structure):
                structure = csc_matrix(structure)
            else:
                structure = np.atleast_2d(structure)

            groups = np.atleast_1d(groups)
            return _sparse_difference(fun_wrapped, x0, f0, h,
                                      use_one_sided, structure,
                                      groups, method)


def _linear_operator_difference(fun, x0, f0, h, method):
    m = f0.size
    n = x0.size

    if method == '2-point':
        def matvec(p):
            if np.array_equal(p, np.zeros_like(p)):
                return np.zeros(m)
            dx = h / norm(p)
            x = x0 + dx*p
            df = fun(x) - f0
            return df / dx

    elif method == '3-point':
        def matvec(p):
            if np.array_equal(p, np.zeros_like(p)):
                return np.zeros(m)
            dx = 2*h / norm(p)
            x1 = x0 - (dx/2)*p
            x2 = x0 + (dx/2)*p
            f1 = fun(x1)
            f2 = fun(x2)
            df = f2 - f1
            return df / dx

    elif method == 'cs':
        def matvec(p):
            if np.array_equal(p, np.zeros_like(p)):
                return np.zeros(m)
            dx = h / norm(p)
            x = x0 + dx*p*1.j
            f1 = fun(x)
            df = f1.imag
            return df / dx

    else:
        raise RuntimeError("Never be here.")

    return LinearOperator((m, n), matvec)


def _dense_difference(fun, x0, f0, h, use_one_sided, method):
    m = f0.size
    n = x0.size
    J_transposed = np.empty((n, m))
    h_vecs = np.diag(h)

    for i in range(h.size):
        if method == '2-point':
            x = x0 + h_vecs[i]
            dx = x[i] - x0[i]  # Recompute dx as exactly representable number.
            df = fun(x) - f0
        elif method == '3-point' and use_one_sided[i]:
            x1 = x0 + h_vecs[i]
            x2 = x0 + 2 * h_vecs[i]
            dx = x2[i] - x0[i]
            f1 = fun(x1)
            f2 = fun(x2)
            df = -3.0 * f0 + 4 * f1 - f2
        elif method == '3-point' and not use_one_sided[i]:
            x1 = x0 - h_vecs[i]
            x2 = x0 + h_vecs[i]
            dx = x2[i] - x1[i]
            f1 = fun(x1)
            f2 = fun(x2)
            df = f2 - f1
        elif method == 'cs':
            f1 = fun(x0 + h_vecs[i]*1.j)
            df = f1.imag
            dx = h_vecs[i, i]
        else:
            raise RuntimeError("Never be here.")

        J_transposed[i] = df / dx

    if m == 1:
        J_transposed = np.ravel(J_transposed)

    return J_transposed.T


def _sparse_difference(fun, x0, f0, h, use_one_sided,
                       structure, groups, method):
    m = f0.size
    n = x0.size
    row_indices = []
    col_indices = []
    fractions = []

    n_groups = np.max(groups) + 1
    for group in range(n_groups):
        # Perturb variables which are in the same group simultaneously.
        e = np.equal(group, groups)
        h_vec = h * e
        if method == '2-point':
            x = x0 + h_vec
            dx = x - x0
            df = fun(x) - f0
            # The result is  written to columns which correspond to perturbed
            # variables.
            cols, = np.nonzero(e)
            # Find all non-zero elements in selected columns of Jacobian.
            i, j, _ = find(structure[:, cols])
            # Restore column indices in the full array.
            j = cols[j]
        elif method == '3-point':
            # Here we do conceptually the same but separate one-sided
            # and two-sided schemes.
            x1 = x0.copy()
            x2 = x0.copy()

            mask_1 = use_one_sided & e
            x1[mask_1] += h_vec[mask_1]
            x2[mask_1] += 2 * h_vec[mask_1]

            mask_2 = ~use_one_sided & e
            x1[mask_2] -= h_vec[mask_2]
            x2[mask_2] += h_vec[mask_2]

            dx = np.zeros(n)
            dx[mask_1] = x2[mask_1] - x0[mask_1]
            dx[mask_2] = x2[mask_2] - x1[mask_2]

            f1 = fun(x1)
            f2 = fun(x2)

            cols, = np.nonzero(e)
            i, j, _ = find(structure[:, cols])
            j = cols[j]

            mask = use_one_sided[j]
            df = np.empty(m)

            rows = i[mask]
            df[rows] = -3 * f0[rows] + 4 * f1[rows] - f2[rows]

            rows = i[~mask]
            df[rows] = f2[rows] - f1[rows]
        elif method == 'cs':
            f1 = fun(x0 + h_vec*1.j)
            df = f1.imag
            dx = h_vec
            cols, = np.nonzero(e)
            i, j, _ = find(structure[:, cols])
            j = cols[j]
        else:
            raise ValueError("Never be here.")

        # All that's left is to compute the fraction. We store i, j and
        # fractions as separate arrays and later construct coo_matrix.
        row_indices.append(i)
        col_indices.append(j)
        fractions.append(df[i] / dx[j])

    row_indices = np.hstack(row_indices)
    col_indices = np.hstack(col_indices)
    fractions = np.hstack(fractions)
    J = coo_matrix((fractions, (row_indices, col_indices)), shape=(m, n))
    return csr_matrix(J)


def check_derivative(fun, jac, x0, bounds=(-np.inf, np.inf), args=(),
                     kwargs={}):
    """Check correctness of a function computing derivatives (Jacobian or
    gradient) by comparison with a finite difference approximation.

    Parameters
    ----------
    fun : callable
        Function of which to estimate the derivatives. The argument x
        passed to this function is ndarray of shape (n,) (never a scalar
        even if n=1). It must return 1-D array_like of shape (m,) or a scalar.
    jac : callable
        Function which computes Jacobian matrix of `fun`. It must work with
        argument x the same way as `fun`. The return value must be array_like
        or sparse matrix with an appropriate shape.
    x0 : array_like of shape (n,) or float
        Point at which to estimate the derivatives. Float will be converted
        to 1-D array.
    bounds : 2-tuple of array_like, optional
        Lower and upper bounds on independent variables. Defaults to no bounds.
        Each bound must match the size of `x0` or be a scalar, in the latter
        case the bound will be the same for all variables. Use it to limit the
        range of function evaluation.
    args, kwargs : tuple and dict, optional
        Additional arguments passed to `fun` and `jac`. Both empty by default.
        The calling signature is ``fun(x, *args, **kwargs)`` and the same
        for `jac`.

    Returns
    -------
    accuracy : float
        The maximum among all relative errors for elements with absolute values
        higher than 1 and absolute errors for elements with absolute values
        less or equal than 1. If `accuracy` is on the order of 1e-6 or lower,
        then it is likely that your `jac` implementation is correct.

    See Also
    --------
    approx_derivative : Compute finite difference approximation of derivative.

    Examples
    --------
    >>> import numpy as np
    >>> from scipy.optimize import check_derivative
    >>>
    >>>
    >>> def f(x, c1, c2):
    ...     return np.array([x[0] * np.sin(c1 * x[1]),
    ...                      x[0] * np.cos(c2 * x[1])])
    ...
    >>> def jac(x, c1, c2):
    ...     return np.array([
    ...         [np.sin(c1 * x[1]),  c1 * x[0] * np.cos(c1 * x[1])],
    ...         [np.cos(c2 * x[1]), -c2 * x[0] * np.sin(c2 * x[1])]
    ...     ])
    ...
    >>>
    >>> x0 = np.array([1.0, 0.5 * np.pi])
    >>> check_derivative(f, jac, x0, args=(1, 2))
    2.4492935982947064e-16
    """
    J_to_test = jac(x0, *args, **kwargs)
    if issparse(J_to_test):
        J_diff = approx_derivative(fun, x0, bounds=bounds, sparsity=J_to_test,
                                   args=args, kwargs=kwargs)
        J_to_test = csr_matrix(J_to_test)
        abs_err = J_to_test - J_diff
        i, j, abs_err_data = find(abs_err)
        J_diff_data = np.asarray(J_diff[i, j]).ravel()
        return np.max(np.abs(abs_err_data) /
                      np.maximum(1, np.abs(J_diff_data)))
    else:
        J_diff = approx_derivative(fun, x0, bounds=bounds,
                                   args=args, kwargs=kwargs)
        abs_err = np.abs(J_to_test - J_diff)
        return np.max(abs_err / np.maximum(1, np.abs(J_diff)))

[((12, 6, 12, 26), 'numpy.finfo', 'np.finfo', ({(12, 15, 12, 25): 'np.float64'}, {}), '(np.float64)', True, 'import numpy as np\n'), ((54, 7, 54, 47), 'numpy.all', 'np.all', ({(54, 14, 54, 46): '((lb == -np.inf) & (ub == np.inf))'}, {}), '((lb == -np.inf) & (ub == np.inf))', True, 'import numpy as np\n'), ((333, 9, 333, 26), 'numpy.atleast_1d', 'np.atleast_1d', ({(333, 23, 333, 25): 'x0'}, {}), '(x0)', True, 'import numpy as np\n'), ((361, 7, 361, 36), 'numpy.any', 'np.any', ({(361, 14, 361, 35): '((x0 < lb) | (x0 > ub))'}, {}), '((x0 < lb) | (x0 > ub))', True, 'import numpy as np\n'), ((441, 11, 441, 41), 'scipy.sparse.linalg.LinearOperator', 'LinearOperator', ({(441, 26, 441, 32): '(m, n)', (441, 34, 441, 40): 'matvec'}, {}), '((m, n), matvec)', False, 'from scipy.sparse.linalg import LinearOperator\n'), ((447, 19, 447, 35), 'numpy.empty', 'np.empty', ({(447, 28, 447, 34): '(n, m)'}, {}), '((n, m))', True, 'import numpy as np\n'), ((448, 13, 448, 23), 'numpy.diag', 'np.diag', ({(448, 21, 448, 22): 'h'}, {}), '(h)', True, 'import numpy as np\n'), ((557, 18, 557, 40), 'numpy.hstack', 'np.hstack', ({(557, 28, 557, 39): 'row_indices'}, {}), '(row_indices)', True, 'import numpy as np\n'), ((558, 18, 558, 40), 'numpy.hstack', 'np.hstack', ({(558, 28, 558, 39): 'col_indices'}, {}), '(col_indices)', True, 'import numpy as np\n'), ((559, 16, 559, 36), 'numpy.hstack', 'np.hstack', ({(559, 26, 559, 35): 'fractions'}, {}), '(fractions)', True, 'import numpy as np\n'), ((47, 24, 47, 51), 'numpy.ones_like', 'np.ones_like', (), '', True, 'import numpy as np\n'), ((107, 14, 107, 40), 'numpy.asarray', 'np.asarray', (), '', True, 'import numpy as np\n'), ((109, 13, 109, 36), 'numpy.resize', 'np.resize', ({(109, 23, 109, 25): 'lb', (109, 27, 109, 35): 'x0.shape'}, {}), '(lb, x0.shape)', True, 'import numpy as np\n'), ((112, 13, 112, 36), 'numpy.resize', 'np.resize', ({(112, 23, 112, 25): 'ub', (112, 27, 112, 35): 'x0.shape'}, {}), '(ub, x0.shape)', True, 'import numpy as np\n'), ((150, 12, 150, 28), 'numpy.atleast_2d', 'np.atleast_2d', ({(150, 26, 150, 27): 'A'}, {}), '(A)', True, 'import numpy as np\n'), ((158, 24, 158, 42), 'numpy.isscalar', 'np.isscalar', ({(158, 36, 158, 41): 'order'}, {}), '(order)', True, 'import numpy as np\n'), ((159, 14, 159, 42), 'numpy.random.RandomState', 'np.random.RandomState', ({(159, 36, 159, 41): 'order'}, {}), '(order)', True, 'import numpy as np\n'), ((162, 16, 162, 33), 'numpy.asarray', 'np.asarray', ({(162, 27, 162, 32): 'order'}, {}), '(order)', True, 'import numpy as np\n'), ((357, 13, 357, 30), 'numpy.atleast_1d', 'np.atleast_1d', ({(357, 27, 357, 29): 'f0'}, {}), '(f0)', True, 'import numpy as np\n'), ((479, 23, 479, 45), 'numpy.ravel', 'np.ravel', ({(479, 32, 479, 44): 'J_transposed'}, {}), '(J_transposed)', True, 'import numpy as np\n'), ((492, 15, 492, 29), 'numpy.max', 'np.max', ({(492, 22, 492, 28): 'groups'}, {}), '(groups)', True, 'import numpy as np\n'), ((495, 12, 495, 35), 'numpy.equal', 'np.equal', ({(495, 21, 495, 26): 'group', (495, 28, 495, 34): 'groups'}, {}), '(group, groups)', True, 'import numpy as np\n'), ((638, 18, 638, 44), 'numpy.abs', 'np.abs', ({(638, 25, 638, 43): 'J_to_test - J_diff'}, {}), '(J_to_test - J_diff)', True, 'import numpy as np\n'), ((49, 12, 49, 21), 'numpy.abs', 'np.abs', ({(49, 19, 49, 20): 'h'}, {}), '(h)', True, 'import numpy as np\n'), ((50, 24, 50, 52), 'numpy.zeros_like', 'np.zeros_like', (), '', True, 'import numpy as np\n'), ((66, 18, 66, 33), 'numpy.abs', 'np.abs', ({(66, 25, 66, 32): 'h_total'}, {}), '(h_total)', True, 'import numpy as np\n'), ((66, 37, 66, 71), 'numpy.maximum', 'np.maximum', ({(66, 48, 66, 58): 'lower_dist', (66, 60, 66, 70): 'upper_dist'}, {}), '(lower_dist, upper_dist)', True, 'import numpy as np\n'), ((77, 30, 78, 62), 'numpy.minimum', 'np.minimum', ({(78, 12, 78, 22): 'h[forward]', (78, 24, 78, 61): '0.5 * upper_dist[forward] / num_steps'}, {}), '(h[forward], 0.5 * upper_dist[forward] / num_steps)', True, 'import numpy as np\n'), ((103, 48, 103, 58), 'numpy.abs', 'np.abs', ({(103, 55, 103, 57): 'x0'}, {}), '(x0)', True, 'import numpy as np\n'), ((397, 21, 397, 42), 'numpy.atleast_1d', 'np.atleast_1d', ({(397, 35, 397, 41): 'groups'}, {}), '(groups)', True, 'import numpy as np\n'), ((503, 20, 503, 33), 'numpy.nonzero', 'np.nonzero', ({(503, 31, 503, 32): 'e'}, {}), '(e)', True, 'import numpy as np\n'), ((82, 32, 83, 64), 'numpy.minimum', 'np.minimum', ({(83, 12, 83, 23): 'h[backward]', (83, 25, 83, 63): '(0.5 * lower_dist[backward] / num_steps)'}, {}), '(h[backward], 0.5 * lower_dist[backward] / num_steps)', True, 'import numpy as np\n'), ((86, 19, 86, 53), 'numpy.minimum', 'np.minimum', ({(86, 30, 86, 40): 'upper_dist', (86, 42, 86, 52): 'lower_dist'}, {}), '(upper_dist, lower_dist)', True, 'import numpy as np\n'), ((395, 28, 395, 52), 'numpy.atleast_2d', 'np.atleast_2d', ({(395, 42, 395, 51): 'structure'}, {}), '(structure)', True, 'import numpy as np\n'), ((409, 33, 409, 49), 'numpy.zeros_like', 'np.zeros_like', ({(409, 47, 409, 48): 'p'}, {}), '(p)', True, 'import numpy as np\n'), ((410, 23, 410, 34), 'numpy.zeros', 'np.zeros', ({(410, 32, 410, 33): 'm'}, {}), '(m)', True, 'import numpy as np\n'), ((411, 21, 411, 28), 'numpy.linalg.norm', 'norm', ({(411, 26, 411, 27): 'p'}, {}), '(p)', False, 'from numpy.linalg import norm\n'), ((522, 17, 522, 28), 'numpy.zeros', 'np.zeros', ({(522, 26, 522, 27): 'n'}, {}), '(n)', True, 'import numpy as np\n'), ((529, 20, 529, 33), 'numpy.nonzero', 'np.nonzero', ({(529, 31, 529, 32): 'e'}, {}), '(e)', True, 'import numpy as np\n'), ((534, 17, 534, 28), 'numpy.empty', 'np.empty', ({(534, 26, 534, 27): 'm'}, {}), '(m)', True, 'import numpy as np\n'), ((632, 22, 632, 46), 'numpy.asarray', 'np.asarray', ({(632, 33, 632, 45): 'J_diff[i, j]'}, {}), '(J_diff[i, j])', True, 'import numpy as np\n'), ((633, 22, 633, 42), 'numpy.abs', 'np.abs', ({(633, 29, 633, 41): 'abs_err_data'}, {}), '(abs_err_data)', True, 'import numpy as np\n'), ((87, 40, 87, 58), 'numpy.abs', 'np.abs', ({(87, 47, 87, 57): 'h_adjusted'}, {}), '(h_adjusted)', True, 'import numpy as np\n'), ((342, 42, 342, 54), 'numpy.isinf', 'np.isinf', ({(342, 51, 342, 53): 'lb'}, {}), '(lb)', True, 'import numpy as np\n'), ((343, 46, 343, 58), 'numpy.isinf', 'np.isinf', ({(343, 55, 343, 57): 'ub'}, {}), '(ub)', True, 'import numpy as np\n'), ((418, 33, 418, 49), 'numpy.zeros_like', 'np.zeros_like', ({(418, 47, 418, 48): 'p'}, {}), '(p)', True, 'import numpy as np\n'), ((419, 23, 419, 34), 'numpy.zeros', 'np.zeros', ({(419, 32, 419, 33): 'm'}, {}), '(m)', True, 'import numpy as np\n'), ((420, 23, 420, 30), 'numpy.linalg.norm', 'norm', ({(420, 28, 420, 29): 'p'}, {}), '(p)', False, 'from numpy.linalg import norm\n'), ((545, 20, 545, 33), 'numpy.nonzero', 'np.nonzero', ({(545, 31, 545, 32): 'e'}, {}), '(e)', True, 'import numpy as np\n'), ((634, 36, 634, 55), 'numpy.abs', 'np.abs', ({(634, 43, 634, 54): 'J_diff_data'}, {}), '(J_diff_data)', True, 'import numpy as np\n'), ((639, 46, 639, 60), 'numpy.abs', 'np.abs', ({(639, 53, 639, 59): 'J_diff'}, {}), '(J_diff)', True, 'import numpy as np\n'), ((430, 33, 430, 49), 'numpy.zeros_like', 'np.zeros_like', ({(430, 47, 430, 48): 'p'}, {}), '(p)', True, 'import numpy as np\n'), ((431, 23, 431, 34), 'numpy.zeros', 'np.zeros', ({(431, 32, 431, 33): 'm'}, {}), '(m)', True, 'import numpy as np\n'), ((432, 21, 432, 28), 'numpy.linalg.norm', 'norm', ({(432, 26, 432, 27): 'p'}, {}), '(p)', False, 'from numpy.linalg import norm\n')]