Opened 6 years ago

# Wrong result from delsarte_bound_additive_hamming_space with GLPK exact simplex

Reported by: Owned by: mkoeppe major sage-wishlist coding theory lp dimpase N/A #20406

`delsarte_bound_additive_hamming_space` 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.

```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
glp_exact: 54 rows, 20 columns, 795 non-zeros
...
2
```

### comment:1 Changed 6 years ago by mkoeppe

(See #20406 where this was discovered.)

### comment:2 Changed 6 years ago by dimpase

• Dependencies set to #20406

### comment:3 Changed 6 years ago by dimpase

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:=7d, 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.

```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
```

here it should have been done, and return 3 (as 73=343), but it is not!

```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
```

### comment:4 Changed 6 years ago by dimpase

• Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.
• Status changed from new to needs_info

in its infinite wisdom, GLPK returns the result of an exact computation as a float. And then the line

```      if q_base**(m+1) == bd:
```

compares, for equality, the exact value, 343=73, with the float 342.999999999999943, and founds them unequal. Whereas for the correct functioning of the algorithm, they must be equal.

Not sure whether this is an upstream bug, or an upstream feature...

### comment:5 follow-up: ↓ 6 Changed 6 years ago by mkoeppe

How come you never run into trouble with the floating-point based solvers with this kind of code?

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.

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.

### comment:6 in reply to: ↑ 5 ; follow-up: ↓ 7 Changed 6 years ago by dimpase

How come you never run into trouble with the floating-point based solvers with this kind of code?

of course I did, and docs explicitly say that as soon as you specify another solver, you are on your own.

```- 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!
```

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.

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.

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.

this design is called lazyness, in CS, and not only :-)

### comment:7 in reply to: ↑ 6 ; follow-up: ↓ 8 Changed 6 years ago by mkoeppe

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.

You can query the `base_ring` of the MIP and then ask `is_exact`.

### comment:8 in reply to: ↑ 7 ; follow-up: ↓ 9 Changed 6 years ago by dimpase

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.

You can query the `base_ring` of the MIP and then ask `is_exact`.

well, `is_exact` 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. `RealField(2000)`.

Besides, I don't think there are places in Sage that forbid doing things cause they are "risky".

### comment:9 in reply to: ↑ 8 Changed 6 years ago by mkoeppe

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.

You can query the `base_ring` of the MIP and then ask `is_exact`.

well, `is_exact` 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. `RealField(2000)`.

Perhaps you mean a `RealIntervalField` here, because certainly `RealField(2000)` 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.

Besides, I don't think there are places in Sage that forbid doing things cause they are "risky".

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.

### comment:10 Changed 6 years ago by mkoeppe

The relevant code is in GLPK's glpapi07.c; it's using `mpq_get_d` to store the rational results in a double.

### comment:11 Changed 6 years ago by mkoeppe

• Milestone changed from sage-7.2 to sage-wishlist
• Report Upstream changed from Not yet reported upstream; Will do shortly. to N/A

Actually, storing the rational results as doubles is fine, in fact I have a test in `glpk_backend` for that. The error is happening when the big coefficients are converted into doubles, and then reconstructed by GLPK's exact simplex.

So there's no upstream bug to be reported -- except that we still really want GLPK to make the interface to glpssx.h public (#18765).

Guarding your code against floating-point fuzz is a wishlist item -- I'm marking this ticket as such.

### comment:12 Changed 5 years ago by jdemeyer

So isn't this essentially a duplicate of #18765 then?

### comment:13 Changed 5 years ago by mkoeppe

• Description modified (diff)

No, #18765 is about making a "proper" rational GLPK API available.

I've changed the description of the present ticket to say what should be done.

### comment:14 Changed 5 years ago by mkoeppe

• Status changed from needs_info to needs_work
Note: See TracTickets for help on using tickets.