def _array_newton(func, x0, fprime, args, tol, maxiter, fprime2,
                  converged=False):
    """
    A vectorized version of Newton, Halley, and secant methods for arrays. Do
    not use this method directly. This method is called from :func:`newton`
    when ``np.isscalar(x0)`` is true. For docstring, see :func:`newton`.
    """
    try:
        p = np.asarray(x0, dtype=float)
    except TypeError:  # can't convert complex to float
        p = np.asarray(x0)
    failures = np.ones_like(p, dtype=bool)  # at start, nothing converged
    nz_der = np.copy(failures)
    if fprime is not None:
        # Newton-Raphson method
        for iteration in range(maxiter):
            # first evaluate fval
            fval = np.asarray(func(p, *args))
            # If all fval are 0, all roots have been found, then terminate
            if not fval.any():
                failures = fval.astype(bool)
                break
            fder = np.asarray(fprime(p, *args))
            nz_der = (fder != 0)
            # stop iterating if all derivatives are zero
            if not nz_der.any():
                break
            # Newton step
            dp = fval[nz_der] / fder[nz_der]
            if fprime2 is not None:
                fder2 = np.asarray(fprime2(p, *args))
                dp = dp / (1.0 - 0.5 * dp * fder2[nz_der] / fder[nz_der])
            # only update nonzero derivatives
            p[nz_der] -= dp
            failures[nz_der] = np.abs(dp) >= tol  # items not yet converged
            # stop iterating if there aren't any failures, not incl zero der
            if not failures[nz_der].any():
                break
    else:
        # Secant method
        dx = np.finfo(float).eps**0.33
        p1 = p * (1 + dx) + np.where(p >= 0, dx, -dx)
        q0 = np.asarray(func(p, *args))
        q1 = np.asarray(func(p1, *args))
        active = np.ones_like(p, dtype=bool)
        for iteration in range(maxiter):
            nz_der = (q1 != q0)
            # stop iterating if all derivatives are zero
            if not nz_der.any():
                p = (p1 + p) / 2.0
                break
            # Secant Step
            dp = (q1 * (p1 - p))[nz_der] / (q1 - q0)[nz_der]
            # only update nonzero derivatives
            p[nz_der] = p1[nz_der] - dp
            active_zero_der = ~nz_der & active
            p[active_zero_der] = (p1 + p)[active_zero_der] / 2.0
            active &= nz_der  # don't assign zero derivatives again
            failures[nz_der] = np.abs(dp) >= tol  # not yet converged
            # stop iterating if there aren't any failures, not incl zero der
            if not failures[nz_der].any():
                break
            p1, p = p, p1
            q0 = q1
            q1 = np.asarray(func(p1, *args))
    zero_der = ~nz_der & failures  # don't include converged with zero-ders
    if zero_der.any():
        # secant warnings
        if fprime is None:
            nonzero_dp = (p1 != p)
            # non-zero dp, but infinite newton step
            zero_der_nz_dp = (zero_der & nonzero_dp)
            if zero_der_nz_dp.any():
                rms = np.sqrt(
                    sum((p1[zero_der_nz_dp] - p[zero_der_nz_dp]) ** 2)
                )
                warnings.warn('RMS of {:g} reached'.format(rms),
                              RuntimeWarning)
        # newton or halley warnings
        else:
            all_or_some = 'all' if zero_der.all() else 'some'
            msg = '{:s} derivatives were zero'.format(all_or_some)
            warnings.warn(msg, RuntimeWarning)
    elif failures.any():
        all_or_some = 'all' if failures.all() else 'some'
        msg = '{0:s} failed to converge after {1:d} iterations'.format(
            all_or_some, maxiter
        )
        if failures.all():
            raise RuntimeError(msg)
        warnings.warn(msg, RuntimeWarning)
    if converged:
        result = namedtuple('result', ('root', 'converged', 'zero_der'))
        p = result(p, ~failures, zero_der)
    return p
def _array_newton ( func , x0 , fprime , args , tol , maxiter , fprime2 , converged = False ) : try : p = np . asarray ( x0 , dtype = float ) except TypeError : p = np . asarray ( x0 ) failures = np . ones_like ( p , dtype = bool ) nz_der = np . copy ( failures ) if fprime is not None : for iteration in range ( maxiter ) : fval = np . asarray ( func ( p , * args ) ) if not fval . any ( ) : failures = fval . astype ( bool ) break fder = np . asarray ( fprime ( p , * args ) ) nz_der = ( fder != 0 ) if not nz_der . any ( ) : break dp = fval [ nz_der ] / fder [ nz_der ] if fprime2 is not None : fder2 = np . asarray ( fprime2 ( p , * args ) ) dp = dp / ( 1.0 - 0.5 * dp * fder2 [ nz_der ] / fder [ nz_der ] ) p [ nz_der ] -= dp failures [ nz_der ] = np . abs ( dp ) >= tol if not failures [ nz_der ] . any ( ) : break else : dx = np . finfo ( float ) . eps ** 0.33 p1 = p * ( 1 + dx ) + np . where ( p >= 0 , dx , - dx ) q0 = np . asarray ( func ( p , * args ) ) q1 = np . asarray ( func ( p1 , * args ) ) active = np . ones_like ( p , dtype = bool ) for iteration in range ( maxiter ) : nz_der = ( q1 != q0 ) if not nz_der . any ( ) : p = ( p1 + p ) / 2.0 break dp = ( q1 * ( p1 - p ) ) [ nz_der ] / ( q1 - q0 ) [ nz_der ] p [ nz_der ] = p1 [ nz_der ] - dp active_zero_der = ~ nz_der & active p [ active_zero_der ] = ( p1 + p ) [ active_zero_der ] / 2.0 active &= nz_der failures [ nz_der ] = np . abs ( dp ) >= tol if not failures [ nz_der ] . any ( ) : break p1 , p = p , p1 q0 = q1 q1 = np . asarray ( func ( p1 , * args ) ) zero_der = ~ nz_der & failures if zero_der . any ( ) : if fprime is None : nonzero_dp = ( p1 != p ) zero_der_nz_dp = ( zero_der & nonzero_dp ) if zero_der_nz_dp . any ( ) : rms = np . sqrt ( sum ( ( p1 [ zero_der_nz_dp ] - p [ zero_der_nz_dp ] ) ** 2 ) ) warnings . warn ( 'RMS of {:g} reached' . format ( rms ) , RuntimeWarning ) else : all_or_some = 'all' if zero_der . all ( ) else 'some' msg = '{:s} derivatives were zero' . format ( all_or_some ) warnings . warn ( msg , RuntimeWarning ) elif failures . any ( ) : all_or_some = 'all' if failures . all ( ) else 'some' msg = '{0:s} failed to converge after {1:d} iterations' . format ( all_or_some , maxiter ) if failures . all ( ) : raise RuntimeError ( msg ) warnings . warn ( msg , RuntimeWarning ) if converged : result = namedtuple ( 'result' , ( 'root' , 'converged' , 'zero_der' ) ) p = result ( p , ~ failures , zero_der ) return p
A vectorized version of Newton, Halley, and secant methods for arrays. Do
    not use this method directly. This method is called from :func:`newton`
    when ``np.isscalar(x0)`` is true. For docstring, see :func:`newton`.