find_local_maximum/minimum() fails with expressions containing complex numbers

[rwst]

sage: find_local_maximum(abs(x+I),-1,1)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-f0919e050ec5> in <module>()
----> 1 find_local_maximum(abs(x+I),-Integer(1),Integer(1))

/home/ralf/sage/src/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport.__call__ (build/cythonized/sage/misc/lazy_import.c:3703)()
    352             True
    353         """
--> 354         return self.get_object()(*args, **kwds)
    355 
    356     def __repr__(self):

/home/ralf/sage/local/lib/python2.7/site-packages/sage/numerical/optimize.pyc in find_local_maximum(f, a, b, tol, maxfun)
    143     """
    144     try:
--> 145         return f.find_local_maximum(a=a, b=b, tol=tol, maxfun=maxfun)
    146     except AttributeError:
    147         pass

/home/ralf/sage/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.find_local_maximum (build/cythonized/sage/symbolic/expression.cpp:66056)()
  11680             (0.561090323458081..., 0.857926501456...)
  11681         """
> 11682         minval, x = (-self).find_local_minimum(a, b, var=var, tol=tol,
  11683                                                      maxfun=maxfun)
  11684         return -minval, x

/home/ralf/sage/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.find_local_minimum (build/cythonized/sage/symbolic/expression.cpp:66379)()
  11739         if var is None:
  11740             var = self.default_variable()
> 11741         return find_local_minimum(self._fast_float_(var),
  11742                                         a=a, b=b, tol=tol, maxfun=maxfun )
  11743 

/home/ralf/sage/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._fast_float_ (build/cythonized/sage/symbolic/expression.cpp:66543)()
  11762         """
  11763         from sage.symbolic.expression_conversions import fast_float
> 11764         return fast_float(self, *vars)
  11765 
  11766     def _fast_callable_(self, etb):

/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in fast_float(ex, *vars)
   1572         1.4142135623730951
   1573     """
-> 1574     return FastFloatConverter(ex, *vars)()
   1575 
   1576 #################

/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
    215             if getattr(self, 'use_fake_div', False) and (operator is _operator.mul or operator is mul_vararg):
    216                 div = self.get_fake_div(ex)
--> 217                 return self.arithmetic(div, div.operator())
    218             return self.arithmetic(ex, operator)
    219         elif operator in relation_operators:

/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in arithmetic(self, ex, operator)
   1513         operands = ex.operands()
   1514         if operator is _operator.neg:
-> 1515             return operator(self(operands[0]))
   1516 
   1517         from sage.rings.all import Rational

/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
    224             return self.tuple(ex)
    225         else:
--> 226             return self.composition(ex, operator)
    227 
    228     def get_fake_div(self, ex):

/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in composition(self, ex, operator)
   1547         """
   1548         f = operator
-> 1549         g = [self(_) for _ in ex.operands()]
   1550         try:
   1551             return f(*g)

/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
    216                 div = self.get_fake_div(ex)
    217                 return self.arithmetic(div, div.operator())
--> 218             return self.arithmetic(ex, operator)
    219         elif operator in relation_operators:
    220             return self.relation(ex, operator)

/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in arithmetic(self, ex, operator)
   1519             from sage.functions.all import sqrt
   1520             return sqrt(self(operands[0]))
-> 1521         fops = map(self, operands)
   1522         if operator == add_vararg:
   1523             operator = _operator.add

/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.pyc in __call__(self, ex)
    206         except TypeError as err:
    207             if 'self must be a numeric expression' not in err.args:
--> 208                 raise err
    209 
    210         operator = ex.operator()

TypeError: unable to coerce to a real number

The rethrown error comes from an attempt to do float(I).

The find_local functions use FastFloatConverter to compute values.

CC: @orlitzky

Component: symbolics

Author: Michael Orlitzky

Branch/Commit: 1fe5795

Reviewer: Matthias Koeppe

Issue created by migration from https://trac.sagemath.org/ticket/24536

[rwst]

comment:2

Apparently "Float" literally means Python (real) "float" not floating-point. Of course then expressions containing I raise errors, even if the outcome is real. The original case was the usage by find_local_maximum---so that never worked, and the restriction was undocumented.

[rwst]

comment:4

The documentation of fast_float states

def fast_float(ex, *vars):
    """
    Returns an object which provides fast floating point evaluation of
    the symbolic expression *ex*.

but

sage: from sage.symbolic.expression_conversions import fast_float
sage: ff = fast_float(abs(x+I))
/home/ralf/sage/local/lib/python2.7/site-packages/sage/symbolic/expression_conversions.py:1574: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...)
See http://trac.sagemath.org/5930 for details.
  return FastFloatConverter(ex, *vars)()

gives the ticket error.

[rwst]

comment:5

Relevant: #5572. Somewhat related: #13559, #16899.

[johanrosenkilde]

comment:6

Note that changing the call from find_local_maximum(abs(x+I),-1,1) to

find_local_maximum(lambda x: abs(x+I),-1,1)

Then everything works fine.

[mkoeppe]

Dependencies: #32234

[orlitzky]

comment:8

The fix for this will be similar to the one in #8450. After #32234, we're using fast_callable() on these expressions with domain=float, which converts all intermediate expression-evaluation results to float. Obviously that causes problems with the x+I in abs(x+I). To work around that, we can probably just run the computation with domain=CDF and then convert the answer to float if its imaginary part is small enough. (But we don't want to return NaN if the end result has a nontrivial imaginary part, like we do when plotting.)

[orlitzky]

comment:9

I'll deal with it once the prereqs start to settle.

[orlitzky]

Branch: u/mjo/ticket/24536

[orlitzky]

comment:13

As promised, I've factored out a superclass of the FastCallablePlotWrapper that won't ignore complex results and will instead raise an error. This is more suitable for functions like find_root() and find_local_minimum() where a complex result can't really be ignored.

Those two functions (and the "maximum" alias) have been updated, but I'm sure there are others where the same strategy can easily be applied.

---

New commits:

925dcf0	Trac #24536: factor out a FastCallableFloatWrapper class from plots.
3b8364b	Trac #24536: support complex expressions in find_{root,local_minimum}.

[orlitzky]

Changed dependencies from #32234 to none

[orlitzky]

Commit: 3b8364b

[orlitzky]

Author: Michael Orlitzky

[orlitzky]

comment:28

Replying to @mkoeppe:

So this should be CDF.coerce_map_from(float).section()?

In spirit, but with a tolerance, so that the the underlying routine doesn't crash when numerical issues cause a 0.0000000000001*I to appear. Additionally, in the plotting class, we choose to return nan upon encountering a complex result so that the corresponding point is "skipped."

[mkoeppe]

comment:29

Yes, so we would have a signaling and a non-signaling version of it

[mkoeppe]

comment:30

Some more random facts:

sage: numpy.real(complex(CDF(1,1)))
1.0
sage: numpy.real(CDF(1,1))
<built-in method real of sage.rings.complex_double.ComplexDoubleElement object at 0x1526bc8c0>

:(

[orlitzky]

comment:31

Replying to @mkoeppe:

NumPy has https://numpy.org/doc/stable/reference/generated/numpy.real_if_close.html

That would basically work here, although I think it might be slower to go from CDF through numpy and back. The main downside to that though is that the error message will be thrown deep in some plotting or optimization routine when it receives a complex number. In the current branch I've chosen to raise an error immediately, in the wrapper, "ValueError: complex fast-callable function result..." that is usually going to better explain the problem.

[mkoeppe]

comment:32

I didn't mean to say that we should use this numpy function. Just collecting existing idioms and interfaces.

I was looking for a name for this operation ("convert to float if close to a real number, NaN otherwise) - and so far haven't found anything in Python, NumPy, C++ reference, ...

[mkoeppe]

comment:33

I guess one would need to figure out what complex infinities should convert to.

[mkoeppe]

comment:34

Also relevant: https://docs.python.org/3/library/cmath.html#cmath.isclose (which accepts both rel_tol and abs_tol). "The IEEE 754 special values of NaN, inf, and -inf will be handled according to IEEE rules. Specifically, [...] inf and -inf are only considered close to themselves." - this forgets to talk about complex infinities.

[mkoeppe]

comment:35

Not sure how well developed the built-in complex type is:

sage: cmath.inf + float(1e-16)*cmath.infj
(nan+infj)

??

[mkoeppe]

comment:36

OK that's https://bugs.python.org/issue25453, open since 2015

[mkoeppe]

comment:37

Here's another one:

sage: cmath.rect(math.inf, 0)
(inf+0j)
sage: cmath.rect(1, math.pi/4)
(0.7071067811865476+0.7071067811865475j)
sage: cmath.rect(math.inf, math.pi/4)
(inf+infj)
sage: cmath.rect(math.inf, pi/2)
(inf+infj)                             # What?
sage: cmath.rect(math.inf, 2*pi)
(inf-infj)                             # Help?!

https://docs.python.org/3/library/cmath.html#cmath.rect

[sagetrac-git]

Changed commit from 40a7012 to 5ef827c

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

5ef827c

Trac #24536: post-refactoring example wording.

[mkoeppe]

comment:39

Replying to @orlitzky:

Replying to @mkoeppe:

Maybe a better name can be found for FastCallablePlotWrapper that expresses better what it does?

So you're not satisfied with solving easy problems all day? =)

It's awkward, and I'm open to suggestions.

Let's postpone this to another ticket. The current branch certainly does what it promises, the implementation looks fine, and tests pass.

[mkoeppe]

Reviewer: Matthias Koeppe

[vbraun]

comment:40

On OSX:

**********************************************************************
File "src/doc/de/thematische_anleitungen/sage_gymnasium.rst", line 534, in doc.de.thematische_anleitungen.sage_gymnasium
Failed example:
    find_root(f, 0.5, 5)
Expected:
    0.6299605249475858
Got:
    0.6299605249475857
**********************************************************************
1 item had failures:
   1 of 294 in doc.de.thematische_anleitungen.sage_gymnasium
    [207 tests, 1 failure, 5.86 s]
**********************************************************************

[sagetrac-git]

Changed commit from 5ef827c to 1fe5795

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

beca559	Trac #24536: factor out a FastCallableFloatWrapper class from plots.
9d5cb78	Trac #24536: support complex expressions in find_{root,local_minimum}.
b9c338e	Trac #24536: add backwards-compatibility case for "complex nan."
b786f52	Trac #24536: post-refactoring example wording.
1fe5795	Trac #24536: add "abs tol" to a finicky doctest.

[orlitzky]

comment:42

Untested (no macOS), but that should take care of it.

[vbraun]

Changed branch from u/mjo/ticket/24536 to 1fe5795