Opened 3 years ago

Closed 3 years ago

#22439 closed defect (fixed)

Fix use of scipy rtol parameter

Reported by: infinity0 Owned by:
Priority: major Milestone: sage-7.6
Component: interfaces Keywords:
Cc: Merged in:
Authors: Ximin Luo Reviewers: Jeroen Demeyer
Report Upstream: N/A Work issues:
Branch: e9ebab2 (Commits) Commit: e9ebab2e0fe8202c35a2f6070bf1ebd876136e05
Dependencies: Stopgaps:

Description (last modified by infinity0)

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

Change History (23)

comment:1 Changed 3 years ago by infinity0

  • Branch set to u/infinity0/fix_use_of_scipy_rtol_parameter

comment:2 Changed 3 years ago by infinity0

  • Authors set to Ximin Luo
  • Commit set to 50a0c646c900eaf24daca8e1aec54bd4eb373145
  • Component changed from PLEASE CHANGE to interfaces
  • Description modified (diff)
  • Status changed from new to needs_review

New commits:

50a0c64Fix Sage's use of scipy rtol parameter

comment:3 Changed 3 years ago by infinity0

  • Type changed from PLEASE CHANGE to defect

comment:4 Changed 3 years ago by jdemeyer

  • Status changed from needs_review to needs_info

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

comment:5 Changed 3 years ago by jdemeyer

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

comment:6 follow-up: Changed 3 years ago by infinity0

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?

comment:7 in reply to: ↑ 6 Changed 3 years ago by jdemeyer

Replying to infinity0:

An alternative is to define default_rtol

Why is that needed?

Why not just change the default value and nothing else?

comment:8 Changed 3 years ago by infinity0

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?

comment:9 Changed 3 years ago by jdemeyer

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

comment:10 follow-ups: Changed 3 years ago by infinity0

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.)

comment:11 in reply to: ↑ 10 Changed 3 years ago by jdemeyer

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?

comment:12 Changed 3 years ago by infinity0

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.

comment:13 in reply to: ↑ 10 Changed 3 years ago by jdemeyer

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.

comment:14 Changed 3 years ago by jdemeyer

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

comment:15 in reply to: ↑ 10 Changed 3 years ago by jdemeyer

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).

comment:16 follow-up: Changed 3 years ago by infinity0

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.

comment:17 in reply to: ↑ 16 Changed 3 years ago by jdemeyer

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.

comment:18 Changed 3 years ago by git

  • Commit changed from 50a0c646c900eaf24daca8e1aec54bd4eb373145 to e9ebab2e0fe8202c35a2f6070bf1ebd876136e05

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

e9ebab2Fix Sage's use of scipy rtol parameter

comment:19 Changed 3 years ago by infinity0

  • Status changed from needs_info to needs_review

comment:20 Changed 3 years ago by infinity0

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:

https://github.com/scipy/scipy/commit/cbd91b93c752eb2e947c2b657fe0471a12759b83#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.)

comment:21 Changed 3 years ago by infinity0

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.

comment:22 Changed 3 years ago by jdemeyer

  • Reviewers set to Jeroen Demeyer
  • Status changed from needs_review to positive_review

comment:23 Changed 3 years ago by vbraun

  • Branch changed from u/infinity0/fix_use_of_scipy_rtol_parameter to e9ebab2e0fe8202c35a2f6070bf1ebd876136e05
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.