Fix use of scipy rtol parameter

[infinity0]

Otherwise it doesn't work with Debian's version of scipy

Component: interfaces

Author: Ximin Luo

Branch/Commit: e9ebab2

Reviewer: Jeroen Demeyer

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

[infinity0]

Branch: u/infinity0/fix_use_of_scipy_rtol_parameter

[infinity0]

Author: Ximin Luo

[infinity0]

Commit: 50a0c64

[infinity0]

New commits:

50a0c64

Fix Sage's use of scipy rtol parameter

[jdemeyer]

comment:4

This seems a bit over-engineered to me. What exactly is the problem?

[jdemeyer]

comment:5

I mean: why not just change the default value of rtol directly?

[infinity0]

comment:6

Oh, that is because sage.numerical.optimize is not imported until we actually run that function. But the default value of a function has to be available when the function is defined, i.e. before it is run.

An alternative is to define default_rtol in some other module and have both sage/numerical/optimize.py and src/sage/symbolic/expression.pyx import this module at the top, unconditionally. Do you have a suggestion on which module this could be?

[jdemeyer]

comment:7

Replying to @infinity0:

An alternative is to define default_rtol

Why is that needed?

Why not just change the default value and nothing else?

[infinity0]

comment:8

To change the default value in both locations as you're suggesting, I would have to write 4*finfo(float).eps in two modules.

The definition comes from the documentation of scipy yet is a calculation on a numpy architecture-dependent property. (The docs say 8.88..e-16, giving away the architecture of the build machine.) I was unsure if this might change in the future, so I only wanted to define it once, and avoid coupling two Sage modules to what looks like a numpy implementation detail but is actually one of scipy's.

It's unfortunate that scipy don't expose the value themselves. Actually, I will send a patch to them for that, when I next get some time. I just noticed again that it's available as from scipy.optimize.zeros import _rtol, i.e. a "private" variable, but this could be properly exposed.

In the meantime, shall I do it your way and duplicate the definitions?

[jdemeyer]

comment:9

Are you sure that this value is architecture-dependent? I think it's just the constant 2.0 ** -50.

[infinity0]

comment:10

I'm not confident how architecture-dependent this is; numpy themselves perform quite a long complex calculation in numpy.core.getlimits and numpy.core.machar, and also Sage currently thinks it's a different value. (The code says 4.5e-16).

I think it would be safer to use the symbolic value as indicated in the scipy documentation, than try to guess that it might be architecture-dependent. (It's true that Sage only works on x86 and x64 right now, but there's no need to add more obstacles to portability.)

[jdemeyer]

comment:11

Replying to @infinity0:

also Sage currently thinks it's a different value. (The code says 4.5e-16).

Sage doesn't "think" it's a different value, Sage just hardcodes a different number.

So why does this default value for rtol even need to be equal to whatever scipy uses?

[infinity0]

comment:12

It needs to be greater than the documented minimum, otherwise we get errors.

From the documentation:

rtol [..] cannot be smaller than its default value of 4*np.finfo(float).eps.

[jdemeyer]

comment:13

Replying to @infinity0:

I'm not confident how architecture-dependent this is

According to Wikipedia, a C double (which is how Python implements float) is required by C99 to be an IEEE-754 double precision number. So that would make it not architecture dependent.

[jdemeyer]

comment:14

In other words: even though it might theoretically be architecture-dependent, in practice it's not.

[jdemeyer]

comment:15

Replying to @infinity0:

It's true that Sage only works on x86 and x64 right now

For the record: we do consider that a bug: #22280. In the past, Sage has worked on SPARC, Itanium, various kinds of PowerPCs, ARM (some of these with a few issues).

[infinity0]

comment:16

OK, I see

typedef struct {
    PyObject_HEAD
    double ob_fval;
} PyFloatObject;

in python's Include/floatobject.h so I guess this is fine. So you want me to rewrite the patch to hard-code 2.0 ** -50 instead? I'll add some comments too.

[jdemeyer]

comment:17

Replying to @infinity0:

So you want me to rewrite the patch to hard-code 2.0 ** -50 instead? I'll add some comments too.

Sounds good.

[sagetrac-git]

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

e9ebab2

Fix Sage's use of scipy rtol parameter

[sagetrac-git]

Changed commit from 50a0c64 to e9ebab2

[infinity0]

comment:20

Oh, I dug a bit deeper and found that scipy changed the meaning of this flag between 0.17.1 and 0.18.1. Sage's current usage is also OK for the older version but not the newer version:

scipy/scipy@cbd91b9#diff-9167bb8d33bd7d9eb5e07e0ebe18e05e

I'm unsure if this patch will cause any extra doctest failures for Sage upstream and scipy 0.17.1, but I guess we'll see that when a build is run.

(On Debian, we are also patching some doctests to expect slightly different float results, but I can't remember which of these are caused by scipy. The differences are minute, of course.)

[infinity0]

comment:21

I guess it may be best to put off this change until you guys upgrade to scipy >= 0.18.1. At least, the extra doc comment I added in that commit is wrong for scipy 0.17.1.

[jdemeyer]

Reviewer: Jeroen Demeyer

[vbraun]

Changed branch from u/infinity0/fix_use_of_scipy_rtol_parameter to e9ebab2