Repository: scipy/scipy
Pull request: FIX: _array_newton should preserve complex inputs
#### Reference issue
Closes gh-11412

#### What does this implement/fix?
Instead of casting the input with `dtype=float`, use `np.promote_types` to promote integer types but leave complex types unchanged. This will also avoid truncating `np.longfloat` input as well.

Problem statement:
Array copying causes unwanted type casting from complex to float in _array_newton
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->
### Problem description
In the first lines of _array_newton where the array p is created by copying the input guess array x0, the handling of complex x0 appears to be broken.  (Lines 377-381)

```   
    try:
        p = np.array(x0, copy=True, dtype=float)
    except TypeError:
        # can't convert complex to float
        p = np.array(x0, copy=True)
```

A complex x0 does not raise a TypeError, but rather you get a ComplexWarning message in the terminal.

```
ComplexWarning: Casting complex values to real discards the imaginary part
```

This means that the p array will be of float type, even when you put in a complex x0. Here is samll example that demonstrates this

```
import numpy as np

t = np.zeros(4,dtype='complex128')
t_1 = 1.0*t
t_2 = np.array(t,copy=True,dtype='float')
print(t.dtype)
print(t_1.dtype)
print(t_2.dtype)
```
The output of this code is
```
demonstrate_type_casting.py:5: ComplexWarning: Casting complex values to real discards the imaginary part
  t_2 = np.array(t,copy=True,dtype='float')
complex128
complex128
float64
```

 The knock-on effect of this is that when the Newton step is performed, you get a TypeError when updating your guesses, since the updates can in general be complex if the function is complex-valued. Here is the error message that I get when I try to run a program that searches for complex roots of a complex-valued function.

#### Code:
```
import numpy as np
import scipy.optimize

def f(x):
    return x + 1+1j
def f_prime(x):
    return 1.0

t = np.ones(4,dtype='complex128')
t *=1j
roots = scipy.optimize.newton(f,t,fprime=f_prime)
print(roots)
```


#### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
/home/edvinolofsson/.local/lib/python3.6/site-packages/scipy/optimize/zeros.py:378: ComplexWarning: Casting complex values to real discards the imaginary part
  p = np.array(x0, copy=True, dtype=float)
Traceback (most recent call last):
  File "demonstrate_type_casting.py", line 11, in <module>
    roots = scipy.optimize.newton(f,t,fprime=f_prime)
  File "/home/edvinolofsson/.local/lib/python3.6/site-packages/scipy/optimize/zeros.py", line 273, in newton
    full_output)
  File "/home/edvinolofsson/.local/lib/python3.6/site-packages/scipy/optimize/zeros.py", line 405, in _array_newton
    p[nz_der] -= dp
numpy.core._exceptions.UFuncTypeError: Cannot cast ufunc 'subtract' output from dtype('complex128') to dtype('float64') with casting rule 'same_kind'

```

### Possible solution
Simply replacing the try except statement with
```
p = np.array(x0, copy=True)
```
appears to fix the issue for me, but I assume that the code was written like this for a reason.


#### Scipy/Numpy/Python version information:
```
Numpy version 1.8.1
Scipy version 1.4.1
Python version 3.6.9
```

Hints from issue discussion:
My assumption is that the dtype=float is to handle the case of starting with an array of integer values, say.  The dtype needed is a type that can store both the provided x0, the desired type of the final result, and the type of any intermediate results.  [To find the zeros of a complex function, the np.array needs to be of complex type, even if all the values are pure real.]
When the code was written perhaps the copy did raise a TypeError but now it only raises ComplexWarning?  If so that would explain why the code was written that way and why it fails now.

Commit message: FIX: _array_newton should handle real x0 and complex f(x0)

Patch to review:
diff --git a/scipy/optimize/tests/test_zeros.py b/scipy/optimize/tests/test_zeros.py
index e3557da28664..6501d9a2e631 100644
--- a/scipy/optimize/tests/test_zeros.py
+++ b/scipy/optimize/tests/test_zeros.py
@@ -341,6 +341,25 @@ def f1_2(x, *a):
         x = zeros.newton(f1, x0, args=args)
         assert_allclose(x, x_expected)
 
+    def test_array_newton_complex(self):
+        def f(x):
+            return x + 1+1j
+
+        def fprime(x):
+            return 1.0
+
+        t = np.full(4, 1j)
+        x = zeros.newton(f, t, fprime=fprime)
+        assert_allclose(f(x), 0.)
+
+        # should work even if x0 is not complex
+        t = np.ones(4)
+        x = zeros.newton(f, t, fprime=fprime)
+        assert_allclose(f(x), 0.)
+
+        x = zeros.newton(f, t)
+        assert_allclose(f(x), 0.)
+
     def test_array_secant_active_zero_der(self):
         """test secant doesn't continue to iterate zero derivatives"""
         x = zeros.newton(lambda x, *a: x*x - a[0], x0=[4.123, 5],
diff --git a/scipy/optimize/zeros.py b/scipy/optimize/zeros.py
index cdf11a516a7f..54834ae43193 100644
--- a/scipy/optimize/zeros.py
+++ b/scipy/optimize/zeros.py
@@ -374,24 +374,24 @@ def _array_newton(func, x0, fprime, args, tol, maxiter, fprime2, full_output):
     """
     # Explicitly copy `x0` as `p` will be modified inplace, but the
     # user's array should not be altered.
-    try:
-        p = np.array(x0, copy=True, dtype=float)
-    except TypeError:
-        # can't convert complex to float
-        p = np.array(x0, copy=True)
+    p = np.array(x0, copy=True)
 
     failures = np.ones_like(p, dtype=bool)
     nz_der = np.ones_like(failures)
     if fprime is not None:
         # Newton-Raphson method
+        fval = np.asarray(func(p, *args))
+        fder = np.asarray(fprime(p, *args))
+        p = np.asarray(p, dtype=np.result_type(p, fval, fder, np.float64))
         for iteration in range(maxiter):
             # first evaluate fval
-            fval = np.asarray(func(p, *args))
+            if iteration != 0:
+                fval = np.asarray(func(p, *args))
+                fder = np.asarray(fprime(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():
@@ -413,6 +413,9 @@ def _array_newton(func, x0, fprime, args, tol, maxiter, fprime2, full_output):
         p1 = p * (1 + dx) + np.where(p >= 0, dx, -dx)
         q0 = np.asarray(func(p, *args))
         q1 = np.asarray(func(p1, *args))
+        dtype = np.result_type(q0, p, np.float64)
+        p = np.asarray(p, dtype=dtype)
+        p1 = np.asarray(p1, dtype=dtype)
         active = np.ones_like(p, dtype=bool)
         for iteration in range(maxiter):
             nz_der = (q1 != q0)


Code review comments:
diff_hunk: @@ -374,24 +374,24 @@ def _array_newton(func, x0, fprime, args, tol, maxiter, fprime2, full_output):
     """
     # Explicitly copy `x0` as `p` will be modified inplace, but the
     # user's array should not be altered.
-    try:
-        p = np.array(x0, copy=True, dtype=float)
-    except TypeError:
-        # can't convert complex to float
-        p = np.array(x0, copy=True)
+    p = np.array(x0, copy=True)
 
     failures = np.ones_like(p, dtype=bool)
     nz_der = np.ones_like(failures)
     if fprime is not None:
         # Newton-Raphson method
+        fval = np.asarray(func(p, *args))
+        fder = np.asarray(fprime(p, *args))
+        p = np.asarray(p, dtype=np.result_type(p, fval, fder, np.float64))
path: scipy/optimize/zeros.py
text: @user1:
I like this way of making the type big enough, but it looks like this only happens on the first iteration.  I.e. If x0, f(x0) and f'(x0) all happen to be real, but a 2nd iteration has complex values, then this won't catch it. 

Does adding it inside the loop instead solve the problem? Perhaps like the following:
```
            ...
            # 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
            # Ensure p has a big enough type
            p = np.asarray(p, dtype=np.result_type(p, dp, np.float64))
            p[nz_der] -= dp
            ...
```

If that works, then iteration 0 wouldn't need to get special cased.


@author:
>I like this way of making the type big enough, but it looks like this only happens on the first iteration. I.e. If x0, f(x0) and f'(x0) all happen to be real, but a 2nd iteration has complex values, then this won't catch it.

Fair enough, unnecessary special cases are generally bad practice.

>Does adding it inside the loop instead solve the problem?

This should be fine. The main issue here is `p[nz_der]` can't do type promotion on `p` so as long as we manually promote `p` before that statement, we should be fine. 

@user1:
LGTM