def as_scipy_functional(func, return_gradient=False):
    """Wrap ``op`` as a function operating on linear arrays.

    This is intended to be used with the `scipy solvers
    <https://docs.scipy.org/doc/scipy/reference/optimize.html>`_.

    Parameters
    ----------
    func : `Functional`.
        A functional that should be wrapped
    return_gradient : bool, optional
        ``True`` if the gradient of the functional should also be returned,
        ``False`` otherwise.

    Returns
    -------
    function : ``callable``
        The wrapped functional.
    gradient : ``callable``, optional
        The wrapped gradient. Only returned if ``return_gradient`` is true.

    Examples
    --------
    Wrap functional and solve simple problem
    (here toy problem ``min_x ||x||^2``):

    >>> func = odl.solvers.L2NormSquared(odl.rn(3))
    >>> scipy_func = odl.as_scipy_functional(func)
    >>> from scipy.optimize import minimize
    >>> result = minimize(scipy_func, x0=[0, 1, 0])
    >>> np.allclose(result.x, [0, 0, 0])
    True

    The gradient (jacobian) can also be provided:

    >>> func = odl.solvers.L2NormSquared(odl.rn(3))
    >>> scipy_func, scipy_grad = odl.as_scipy_functional(func, True)
    >>> from scipy.optimize import minimize
    >>> result = minimize(scipy_func, x0=[0, 1, 0], jac=scipy_grad)
    >>> np.allclose(result.x, [0, 0, 0])
    True

    Notes
    -----
    If the data representation of ``op``'s domain is of type
    `NumpyTensorSpace`, this incurs no significant overhead. If the space type
    is ``CudaFn`` or some other nonlocal type, the overhead is significant.
    """
    def func_call(arr):
        return func(np.asarray(arr).reshape(func.domain.shape))

    if return_gradient:
        def func_gradient_call(arr):
            return np.asarray(
                func.gradient(np.asarray(arr).reshape(func.domain.shape)))

        return func_call, func_gradient_call
    else:
        return func_call
Wrap ``op`` as a function operating on linear arrays.

    This is intended to be used with the `scipy solvers
    <https://docs.scipy.org/doc/scipy/reference/optimize.html>`_.

    Parameters
    ----------
    func : `Functional`.
        A functional that should be wrapped
    return_gradient : bool, optional
        ``True`` if the gradient of the functional should also be returned,
        ``False`` otherwise.

    Returns
    -------
    function : ``callable``
        The wrapped functional.
    gradient : ``callable``, optional
        The wrapped gradient. Only returned if ``return_gradient`` is true.

    Examples
    --------
    Wrap functional and solve simple problem
    (here toy problem ``min_x ||x||^2``):

    >>> func = odl.solvers.L2NormSquared(odl.rn(3))
    >>> scipy_func = odl.as_scipy_functional(func)
    >>> from scipy.optimize import minimize
    >>> result = minimize(scipy_func, x0=[0, 1, 0])
    >>> np.allclose(result.x, [0, 0, 0])
    True

    The gradient (jacobian) can also be provided:

    >>> func = odl.solvers.L2NormSquared(odl.rn(3))
    >>> scipy_func, scipy_grad = odl.as_scipy_functional(func, True)
    >>> from scipy.optimize import minimize
    >>> result = minimize(scipy_func, x0=[0, 1, 0], jac=scipy_grad)
    >>> np.allclose(result.x, [0, 0, 0])
    True

    Notes
    -----
    If the data representation of ``op``'s domain is of type
    `NumpyTensorSpace`, this incurs no significant overhead. If the space type
    is ``CudaFn`` or some other nonlocal type, the overhead is significant.