Doc tests - why not write the exact of high precision result?

7 views
Skip to first unread message

Dr. David Kirkby

unread,
Feb 20, 2010, 6:29:59 PM2/20/10
to sage-...@googlegroups.com
The thread 'Sage 4.3.3.alpha1 released' which covers loads of stuff. The
particular bit I'm referring to is this test failure failure reported by
Robert Marik.

----------------------------------------------------
----------------------------------------------------
----------------------------------------------------
sage -t "devel/sage/sage/misc/functional.py"
**********************************************************************
File "/opt/sage-4.3.3.alpha1/devel/sage/sage/misc/functional.py", line
705:
sage: h.n()
Expected:
0.33944794097891573
Got:
0.33944794097891567
**********************************************************************
1 items had failures:
1 of 14 in __main__.example_28

Robert Marik
----------------------------------------------------
----------------------------------------------------
----------------------------------------------------

The test reads like this:

##############################################
Numerical approximation::

sage: h = integral(sin(x)/x^2, (x, 1, pi/2)); h
integrate(sin(x)/x^2, x, 1, 1/2*pi)
sage: h.n()
0.33944794097891573
#################################################


The test certainly does *not* test Sage's ability to correctly calculate that
integral. It tests the documentation gives the same answers each, but that is
not a lot of use of those answers are incorrect. Nobody appears to have bothered
to evaluate if they are correct or not.

It might verify that Sage gives the same results each time, but without knowing
what the exact result is, this is not a very good test.

I think the doc test would be a *lot* more useful if there were comments in the
test like

# In Mathematica

# In[11]:= N[Integrate[ Sin[x]/x^2,{x,1,Pi/2}],50]
# Out[11]= 0.33944794097891567969192717186521861799447698826918

# This can also be computed in Wolfram Alpha, which uses Mathematica
# as a back end.
#
http://www.wolframalpha.com/input/?i=N[Integrate[+Sin[x]%2Fx^2%2C{x%2C1%2CPi%2F2}]%2C50]

# With mpmath

# sage: import mpmath
# sage: mpmath.mp.dps = 50
# sage: mpmath.quad(lambda x: mpmath.sin(x)/x**2, [1, mpmath.pi/2])
# mpf('0.33944794097891567969192717186521861799447698826917531')

# With Maple:
# some Result computed by Maple to arbitrary precision if
# someone knows how to do it.


Had this been done, then nobody in their right mind would have put the expected
result as 0.33944794097891573, since they would have known those last two digits
were inaccurate.

I'm aware this might not be possible in all cases. If that is so, then some
explanation of why not, and why the result seems 'reasonable' would be
appropriate. I'm not a mathematician, but I'm well aware that for some
calculations you can put an upper and lower limit on a result given theories X
and Y.

I for one would have increased confidence in Sage if I could see that people had
gone to the trouble of making independent checks of the results in the
documentation.

Dave

Message has been deleted

Alex Ghitza

unread,
Feb 20, 2010, 7:12:36 PM2/20/10
to Dr. David Kirkby, sage-...@googlegroups.com
On Sat, 20 Feb 2010 23:29:59 +0000, "Dr. David Kirkby" <david....@onetel.net> wrote:
> I think the doc test would be a *lot* more useful if there were comments in the
> test like
>
> # In Mathematica
>
> # In[11]:= N[Integrate[ Sin[x]/x^2,{x,1,Pi/2}],50]
> # Out[11]= 0.33944794097891567969192717186521861799447698826918
>
> # This can also be computed in Wolfram Alpha, which uses Mathematica
> # as a back end.
> #
> http://www.wolframalpha.com/input/?i=N[Integrate[+Sin[x]%2Fx^2%2C{x%2C1%2CPi%2F2}]%2C50]
>

We really need to be careful about this. Wolfram|Alpha has very
specific terms of use, which I think could be interpreted as forbidding
precisely the use that you suggest. I quote from one of the relevant
sections (there might be others):


Ways You May Use Our Free Service and Its Results

The free Wolfram|Alpha service is available for ad hoc, personal,
non-commercial use only. For such use you are welcome to download
results, print copies, store downloaded content on your computer, and
reference this information in your documents. You may use it to get
information for your own use for any purpose, including occasional
purposes related to your job, as long as you aren't specifically being
paid to use Wolfram|Alpha.

Systematic professional or commercial use of the website, or use for
which you are being specifically paid, requires a commercial license. If
you are interested in a version of Wolfram|Alpha that allows such use,
please contact us.

It is permitted to use and post individual, incidental results or small
groups of results from Wolfram|Alpha on non-commercial websites and
blogs, provided those results are properly credited to Wolfram|Alpha
(see Attribution and Licensing below). It is not permitted to extract
multiple pieces of data and reassemble them into anything with the
characteristics of a data service, database, or source of large and
systematic collections of data. You are not allowed to use Wolfram|Alpha
to create something that is likely or intended to be reused as a data
source for further processing, or that in some other way serves as a
replacement or alternative to using Wolfram|Alpha itself. This applies
whether what you create is in electronic or print form.


One could argue that we are creating something that is intended to
replace or serve as an alternative to Wolfram|Alpha. In which case we
better not use their results in Sage (not even in doctests) or we're in
trouble.

(Note that similar terms of use might apply to Mathematica itself or
other commercial math software packages; I don't have any available at
this very moment to check.)


We can however use mpmath, octave, or whatever other suitably-licensed
software for independent verifications of numerical results.


Best,
Alex

--
Alex Ghitza -- http://aghitza.org/
Lecturer in Mathematics -- The University of Melbourne -- Australia

Dr. David Kirkby

unread,
Feb 20, 2010, 8:34:33 PM2/20/10
to sage-...@googlegroups.com

It would be a bit of a stretch of the imagination to suggest Sage is a
replacement or alternative to using Wolfram Alpha. I can't see anyone wishing to
download 260 MB of source code, or a 500 MB binary to replace or serve as an
alternative to Wolfram Alpha.

It should also be noted that Sage existed before Wolfram Alpha did and Sage's
mission has not been revised in any way since Wolfram Alpha was released.

But I can see you may have a point. (As always, a legal opinion is needed on
this). It is certainly good you bought the point up though.

It would be good if William emailed Wolfram Research and asked about this, as
it's not clear.

It could be a useful to give increased confidence.

> (Note that similar terms of use might apply to Mathematica itself or
> other commercial math software packages; I don't have any available at
> this very moment to check.)

I can't find a similar term.

http://www.wolfram.com/terms/MathematicaLicenseAgreement.pdf

I suspect a lot of people developing Sage have access to Mathematica or other
similar software, and could verify results. Or if not, ask on sage-devel for
someone else to verify a result using whatever (hardware, software, brain power
...) to verify the results.

> We can however use mpmath, octave, or whatever other suitably-licensed
> software for independent verifications of numerical results.

IMHO, the usefulness of a test is dramatically reduced if the expected result is
not independently verified.

> Best,
> Alex
>

Dave

Mike Hansen

unread,
Feb 20, 2010, 8:42:42 PM2/20/10
to sage-...@googlegroups.com
On Sat, Feb 20, 2010 at 3:29 PM, Dr. David Kirkby
<david....@onetel.net> wrote:
> The test reads like this:
>
> ##############################################
>    Numerical approximation::
>
>        sage: h = integral(sin(x)/x^2, (x, 1, pi/2)); h
>        integrate(sin(x)/x^2, x, 1, 1/2*pi)
>        sage: h.n()
>        0.33944794097891573
> #################################################
>
>
> The test certainly does *not* test Sage's ability to correctly calculate
> that integral. It tests the documentation gives the same answers each, but
> that is not a lot of use of those answers are incorrect. Nobody appears to
> have bothered to evaluate if they are correct or not.

This isn't supposed to test Sage's ability to correctly calculate that
integral. It's supposed to test the code which calls the numerical
integration routine when you call .n() on an unevaluated integral.

Tests for correctness should be done at a lower level.

--Mike

Robert Bradshaw

unread,
Feb 20, 2010, 9:00:10 PM2/20/10
to sage-...@googlegroups.com

Part of refereeing a patch is making sure that the tests in the
documentation are correct--if people are not doing this then they are
not doing a proper job refereeing what goes in. This means more than
just trying out the relevant command and seeing that it gives the
predicted answer. (Note that here the test is really that .n() works
on an integral, not that the numerical routines used are stable.)

> It might verify that Sage gives the same results each time, but
> without knowing what the exact result is, this is not a very good
> test.

The exact answer can't actually be written down. I wouldn't say that
this answer is incorrect, just not as precise as you might want. Now
you might want every printed digit to be correct (which I think is a
valid concern, and I think n() should have this property) but it's not
a distinction that is currently made in Sage, and it's not so easy to
make in general. For example, every digit in

sage: N = 1000
sage: sin(RR(N))
0.826879540532003

is correct. However, for sufficiently large N, the results

sage: sum(sin(RR(n)) for n in range(0, N))
-0.0129099064588385
sage: sum(sin(RealField(100)(n)) for n in range(0, N))
-0.012909906458836372437031510981

are not correct in every digit. Any inexact numerical computation
should be taken with a grain of salt.

> I think the doc test would be a *lot* more useful if there were
> comments in the test like
>
> # In Mathematica
>
> # In[11]:= N[Integrate[ Sin[x]/x^2,{x,1,Pi/2}],50]
> # Out[11]= 0.33944794097891567969192717186521861799447698826918
>
> # This can also be computed in Wolfram Alpha, which uses Mathematica
> # as a back end.
> # http://www.wolframalpha.com/input/?i=N[Integrate[+Sin[x]%2Fx^2%2C{x
> %2C1%2CPi%2F2}]%2C50]
>
> # With mpmath
>
> # sage: import mpmath
> # sage: mpmath.mp.dps = 50
> # sage: mpmath.quad(lambda x: mpmath.sin(x)/x**2, [1, mpmath.pi/2])
> # mpf('0.33944794097891567969192717186521861799447698826917531')
>
> # With Maple:
> # some Result computed by Maple to arbitrary precision if
> # someone knows how to do it.
>
>
> Had this been done, then nobody in their right mind would have put
> the expected result as 0.33944794097891573, since they would have
> known those last two digits were inaccurate.

The answer "0.339447940978915..." succinctly lets the user know there
are some precision issues here.

> I'm aware this might not be possible in all cases. If that is so,
> then some explanation of why not, and why the result seems
> 'reasonable' would be appropriate. I'm not a mathematician, but I'm
> well aware that for some calculations you can put an upper and lower
> limit on a result given theories X and Y.
>
> I for one would have increased confidence in Sage if I could see
> that people had gone to the trouble of making independent checks of
> the results in the documentation.

I agree here--the best doctests are those that have internal
consistency checks. This could involve trying out the numerical
routines against known closed-form solutions, verifying properties
such as linearity and additivity (especially if those are not obvious
consequences of the algorithms used), computing the same result to
higher precision, etc. In number theory, for example, there are often
strong consistency checks that are consequences of deep theorems that
can be demonstrated and verified (and that provides a more educational
docstring as well).

- Robert

Dr. David Kirkby

unread,
Feb 20, 2010, 9:24:55 PM2/20/10
to sage-...@googlegroups.com

It's not clear to me how one can test whether .n() is working on an unevaluated
integral correctly, unless you have checked that the result returned by .n() on
that unevaluated integral is correct.

Dave

Reply all
Reply to author
Forward
0 new messages