Sage: Ticket #20447: Wrong result from delsarte_bound_additive_hamming_space with GLPK exact simplex
https://trac.sagemath.org/ticket/20447
<p>
<code>delsarte_bound_additive_hamming_space</code> should be guarded, using epsilons, against floating point fuzz that will appear with numerical solvers and even the GLPK exact solver because of the nature of its interface.
</p>
<pre class="wiki">sage: delsarte_bound_additive_hamming_space(19,15,7,isinteger=False)
3
sage: from sage.numerical.backends.generic_backend import get_solver
sage: def glpk_exact_solver():
b = get_solver(solver="GLPK")
b.solver_parameter("simplex_or_intopt", "exact_simplex_only")
return b
sage: delsarte_bound_additive_hamming_space(19,15,7,solver=glpk_exact_solver,isinteger=False)
glp_exact: 54 rows, 20 columns, 795 non-zeros
...
2
</pre>en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/20447
Trac 1.1.6mkoeppeThu, 14 Apr 2016 22:53:11 GMT
https://trac.sagemath.org/ticket/20447#comment:1
https://trac.sagemath.org/ticket/20447#comment:1
<p>
(See <a class="closed ticket" href="https://trac.sagemath.org/ticket/20406" title="enhancement: get_solver should allow passing a function (a solver factory) as the ... (closed: fixed)">#20406</a> where this was discovered.)
</p>
TicketdimpaseFri, 15 Apr 2016 10:30:34 GMTdependencies set
https://trac.sagemath.org/ticket/20447#comment:2
https://trac.sagemath.org/ticket/20447#comment:2
<ul>
<li><strong>dependencies</strong>
set to <em>#20406</em>
</li>
</ul>
TicketdimpaseFri, 15 Apr 2016 11:06:31 GMT
https://trac.sagemath.org/ticket/20447#comment:3
https://trac.sagemath.org/ticket/20447#comment:3
<p>
the function performs the following, in a loop: maximises the sum of variables, gets the value of this maximum, say M, then rounds M down to the closest m:=7<sup>d</sup>, adds a constraint that the sum of the variables is at most m, and repeats. In the case of rational variables, this process obviously must stop after one iteration (it need not stop if all the variables are integers). When it stops, it returns the latest d.
</p>
<pre class="wiki">glp_exact: 54 rows, 20 columns, 795 non-zeros
GNU MP bignum library is being used
0: infsum = 1 (14)
1: infsum = 0 (14)
* 1: objval = 1 (14)
* 19: objval = 1671.30573248408 (0)
OPTIMAL SOLUTION FOUND
glp_exact: 54 rows, 20 columns, 795 non-zeros
GNU MP bignum library is being used
0: infsum = 1 (14)
1: infsum = 0 (14)
* 1: objval = 1 (14)
* 17: objval = 343 (0)
OPTIMAL SOLUTION FOUND
</pre><p>
here it should have been done, and return 3 (as 7<sup>3</sup>=343), but it is not!
</p>
<pre class="wiki">glp_exact: 54 rows, 20 columns, 795 non-zeros
GNU MP bignum library is being used
0: infsum = 1 (14)
1: infsum = 0 (14)
* 1: objval = 1 (14)
* 16: objval = 49 (0)
OPTIMAL SOLUTION FOUND
2
</pre>
TicketdimpaseFri, 15 Apr 2016 11:36:53 GMTstatus, upstream changed
https://trac.sagemath.org/ticket/20447#comment:4
https://trac.sagemath.org/ticket/20447#comment:4
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_info</em>
</li>
<li><strong>upstream</strong>
changed from <em>N/A</em> to <em>Not yet reported upstream; Will do shortly.</em>
</li>
</ul>
<p>
in its infinite wisdom, GLPK returns the result of an exact computation as a float. And then the line
</p>
<pre class="wiki"> if q_base**(m+1) == bd:
</pre><p>
compares, for equality, the exact value, 343=7<sup>3</sup>, with the float 342.999999999999943, and founds them unequal. Whereas for the correct functioning of the algorithm, they must be equal.
</p>
<p>
Not sure whether this is an upstream bug, or an upstream feature...
</p>
TicketmkoeppeFri, 15 Apr 2016 20:08:56 GMT
https://trac.sagemath.org/ticket/20447#comment:5
https://trac.sagemath.org/ticket/20447#comment:5
<p>
How come you never run into trouble with the floating-point based solvers with this kind of code?
</p>
<p>
It is certainly ironic that the exact solver has more floating-point fuzz than the floating-point solvers. But your code cannot be correct if solver is a floating-point solver.
</p>
<p>
And there is an apparent upstream issue with GLPK. I wouldn't call it a bug, but given the design decision to return the exact values as double-floats, it should at least be improved to return all integers that have an exact representation in double-float as such.
</p>
TicketdimpaseFri, 15 Apr 2016 20:19:26 GMT
https://trac.sagemath.org/ticket/20447#comment:6
https://trac.sagemath.org/ticket/20447#comment:6
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20447#comment:5" title="Comment 5">mkoeppe</a>:
</p>
<blockquote class="citation">
<p>
How come you never run into trouble with the floating-point based solvers with this kind of code?
</p>
</blockquote>
<p>
of course I did, and docs explicitly say that as soon as you specify another solver, you are on your own.
</p>
<pre class="wiki">- solver – the LP/ILP solver to be used. Defaults to PPL. It is arbitrary precision, thus there will be no rounding errors. With other solvers (see MixedIntegerLinearProgram for the list), you are on your own!
</pre><blockquote class="citation">
<p>
It is certainly ironic that the exact solver has more floating-point fuzz than the floating-point solvers. But your code cannot be correct if solver is a floating-point solver.
</p>
</blockquote>
<p>
I don't think there is an API to check whether a solver is exact, and so I never bothered to check this in the code.
</p>
<blockquote class="citation">
<p>
And there is an apparent upstream issue with GLPK. I wouldn't call it a bug, but given the design decision to return the exact values as double-floats, it should at least be improved to return all integers that have an exact representation in double-float as such.
</p>
</blockquote>
<p>
this design is called lazyness, in CS, and not only :-)
</p>
TicketmkoeppeFri, 15 Apr 2016 20:22:04 GMT
https://trac.sagemath.org/ticket/20447#comment:7
https://trac.sagemath.org/ticket/20447#comment:7
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20447#comment:6" title="Comment 6">dimpase</a>:
</p>
<blockquote class="citation">
<p>
I don't think there is an API to check whether a solver is exact, and so I never bothered to check this in the code.
</p>
</blockquote>
<p>
You can query the <code>base_ring</code> of the MIP and then ask <code>is_exact</code>.
</p>
TicketdimpaseFri, 15 Apr 2016 20:39:59 GMT
https://trac.sagemath.org/ticket/20447#comment:8
https://trac.sagemath.org/ticket/20447#comment:8
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20447#comment:7" title="Comment 7">mkoeppe</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20447#comment:6" title="Comment 6">dimpase</a>:
</p>
<blockquote class="citation">
<p>
I don't think there is an API to check whether a solver is exact, and so I never bothered to check this in the code.
</p>
</blockquote>
<p>
You can query the <code>base_ring</code> of the MIP and then ask <code>is_exact</code>.
</p>
</blockquote>
<p>
well, <code>is_exact</code> is a bit too much to ask, one merely needs extra precision.
(one can figure out how much, in fact).
So you can use a backend that allows you to set the base ring to e.g. <code>RealField(2000)</code>.
</p>
<p>
Besides, I don't think there are places in Sage that forbid doing things cause they are "risky".
</p>
TicketmkoeppeFri, 15 Apr 2016 23:05:39 GMT
https://trac.sagemath.org/ticket/20447#comment:9
https://trac.sagemath.org/ticket/20447#comment:9
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20447#comment:8" title="Comment 8">dimpase</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20447#comment:7" title="Comment 7">mkoeppe</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20447#comment:6" title="Comment 6">dimpase</a>:
</p>
<blockquote class="citation">
<p>
I don't think there is an API to check whether a solver is exact, and so I never bothered to check this in the code.
</p>
</blockquote>
<p>
You can query the <code>base_ring</code> of the MIP and then ask <code>is_exact</code>.
</p>
</blockquote>
<p>
well, <code>is_exact</code> is a bit too much to ask, one merely needs extra precision.
(one can figure out how much, in fact).
So you can use a backend that allows you to set the base ring to e.g. <code>RealField(2000)</code>.
</p>
</blockquote>
<p>
Perhaps you mean a <code>RealIntervalField</code> here, because certainly <code>RealField(2000)</code> does not guarantee that the result of some unspecified numerical algorithm such as the implementation of the simplex method in the solver leads to results with 2000 correct bits, just like double-floats don't guarantee 53 correct bits.
</p>
<blockquote class="citation">
<p>
Besides, I don't think there are places in Sage that forbid doing things cause they are "risky".
</p>
</blockquote>
<p>
I think there more places should forbid things like that, or at least display a warning. For example, the polyhedral code in Sage is written in a way that it assumes exact arithmetic -- and if fed with floating point numbers, leads to mysterious errors. One would usually assume from algorithms that accept floating point numbers that they have some (however naive) accommodation for floating point fuzz, in the form of some epsilons.
</p>
TicketmkoeppeSat, 16 Apr 2016 07:36:14 GMT
https://trac.sagemath.org/ticket/20447#comment:10
https://trac.sagemath.org/ticket/20447#comment:10
<p>
The relevant code is in GLPK's glpapi07.c; it's using <code>mpq_get_d</code> to store the rational results in a double.
</p>
TicketmkoeppeSat, 16 Apr 2016 16:08:27 GMTupstream, milestone changed
https://trac.sagemath.org/ticket/20447#comment:11
https://trac.sagemath.org/ticket/20447#comment:11
<ul>
<li><strong>upstream</strong>
changed from <em>Not yet reported upstream; Will do shortly.</em> to <em>N/A</em>
</li>
<li><strong>milestone</strong>
changed from <em>sage-7.2</em> to <em>sage-wishlist</em>
</li>
</ul>
<p>
Actually, storing the rational results as doubles is fine, in fact I have a test in <code>glpk_backend</code> for that.
The error is happening when the big coefficients are converted into doubles, and then reconstructed by GLPK's exact simplex.
</p>
<p>
So there's no upstream bug to be reported -- except that we still really want GLPK to make the interface to glpssx.h public (<a class="new ticket" href="https://trac.sagemath.org/ticket/18765" title="enhancement: Add Cython wrappers for GLPK's interface glpssx.h (exact rational simplex) (new)">#18765</a>).
</p>
<p>
Guarding your code against floating-point fuzz is a wishlist item -- I'm marking this ticket as such.
</p>
TicketjdemeyerThu, 31 Aug 2017 15:18:57 GMT
https://trac.sagemath.org/ticket/20447#comment:12
https://trac.sagemath.org/ticket/20447#comment:12
<p>
So isn't this essentially a duplicate of <a class="new ticket" href="https://trac.sagemath.org/ticket/18765" title="enhancement: Add Cython wrappers for GLPK's interface glpssx.h (exact rational simplex) (new)">#18765</a> then?
</p>
TicketmkoeppeFri, 01 Sep 2017 05:44:39 GMTdescription changed
https://trac.sagemath.org/ticket/20447#comment:13
https://trac.sagemath.org/ticket/20447#comment:13
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/20447?action=diff&version=13">diff</a>)
</li>
</ul>
<p>
No, <a class="new ticket" href="https://trac.sagemath.org/ticket/18765" title="enhancement: Add Cython wrappers for GLPK's interface glpssx.h (exact rational simplex) (new)">#18765</a> is about making a "proper" rational GLPK API available.
</p>
<p>
I've changed the description of the present ticket to say what should be done.
</p>
TicketmkoeppeFri, 01 Sep 2017 05:47:21 GMTstatus changed
https://trac.sagemath.org/ticket/20447#comment:14
https://trac.sagemath.org/ticket/20447#comment:14
<ul>
<li><strong>status</strong>
changed from <em>needs_info</em> to <em>needs_work</em>
</li>
</ul>
Ticket