Possible bug in find_maximum_on_interval

[fidelbc]

Hello all,

I was trying to answer a question on ask sage and I might have ran
into a bug.

http://ask.sagemath.org/question/1308/get_minmax_data-too-generosous

The following statements work with no problem:

sage: f(x) = sin(x)
sage: find_minimum_on_interval(f,0,1)
(6.5078162602101728e-09, 6.5078162602101728e-09)
sage: find_maximum_on_interval(f,0,1)
(0.841470969301, 0.99999997130024143)
sage: find_minimum_on_interval(f,0,2)
(8.740260678362169e-09, 8.740260678362169e-09)

the error message arises when executing

sage: find_maximum_on_interval(f,0,2)
[some output ommited]
RuntimeError: ECL says: THROW: The catch MACSYMA-QUIT is undefined.

Is this a bug?

Thanks,
Fidel

[Nils Bruin]

Yes, that's a bug. Something goes wrong in what gets fed to the
function:

def f(z):
v = sin(x)(x=z)
print "evaluating sin(%s)=%s"%(z,v)
return vsage: find_maximum_on_interval(f,0,2)
evaluating sin(0.7639320225)=0.691766183638
evaluating sin(1.2360679775)=0.944499585162
evaluating sin(1.527864045)=0.999078551136
evaluating sin(1.59062900282 == 1.58711797464)=sin(1.59062900282 ==
1.58711797464)

As you see, for some reason the evaluation point becomes an equality
relation. From that point it is of course just a question *which*
piece of code throws an error.

These routines use scipy.optimize.fminbound. This doesn't act too
nicely with the SymbolicRing. In the source, we find the following
lines:

...
si = numpy.sign(xm-xf) + ((xm-xf)==0)
...
si = numpy.sign(rat) + (rat == 0)
x = xf + si*max([abs(rat), tol1])
...

I think this code is highly suspicious: They're comparing what should
be floats to 0??? and they're adding 0 or 1 to something based on
whether something is exactly 0? I'd say the scipy guys should
reconsider their implementation, but it may be they can justify their
steps.

Anyway, the problem comes from the fact that the symbolic ring is
quite conservative in deciding whether floats are equal to zero and
has the habit of returning something that is not readily 0 or 1:

sage: 1.1 == 0
False
sage: SR(1.1) == 0
1.10000000000000 == 0

So at the least, the sage interface prepares input to scipy that
doesn't behave as scipy expects. It may also be that the scipy
implementation is suspect.

The nasty error comes from the fact that the syntactically invalid
expression sin(1.59062900282 == 1.58711797464) doesn't get rejected
upon construction.

[Nils Bruin]

> ...
>                     si = numpy.sign(xm-xf) + ((xm-xf)==0)
> ...
>         si = numpy.sign(rat) + (rat == 0)
>         x = xf + si*max([abs(rat), tol1])
> ...

Ah shoot. That code is quite understandable. What they need is:

....
si = numpy.sign(xm-xf)
if si == 0:
si = 1
...
si = numpy.sign(rat)
if si == 0:
si = 1

x = xf + si*max([abs(rat), tol1])
...

Their assumption that numpy.sign(rat) == 0 iff (rat == 0) evaluates to
1 and that (rat == 0) evaluates to 0 otherwise is not unreasonable in
their world. The above code is of course a little more robust, but I
doubt scipy is willing to change it just because sage's SR acts funny.
So the fix: sage should interface a little more carefully with scipy
and take care no SR objects leak into it.

[Eviatar]

Is there a ticket for this? I've also encountered this before.

[Dan Drake]

On Sun, 08 Apr 2012 at 05:41PM -0700, Eviatar wrote:
> Is there a ticket for this? I've also encountered this before.

Let me add a "me too". I ran into this just a few days ago. Looks like a
ticket for it is: http://trac.sagemath.org/sage_trac/ticket/12032

Dan

--
--- Dan Drake
----- http://mathsci.kaist.ac.kr/~drake
-------

[Keshav Kini]

Would replacing that code with your code be detrimental to scipy? Why
not submit a pull request? I submitted a "please be more robust in this
weird corner case so that Sage doesn't break you with its custom
objects" type pull request to simplejson [1], for example, and it was
kindly accepted by the maintainer. Maybe the same would happen with
scipy.

In any case, one could also argue that implicitly casting bools to
integers is not very Pythonic, if one were feeling a bit impish :P Then
again maybe there are some performance-related reasons, invisible to me,
for the way the code is.

-Keshav

[1]: http://github.com/simplejson/simplejson/pull/29/files