[PPL-devel] PPL's handling of floating-point coefficients

[Enea Zaffanella]

Sebastian Nowozin wrote:
> 
> Hello,
> 
> thanks for making the PPL library open-source software.  I have been 
> using it for a few days now and it works quite well.
> 
> For solving an optimization problem I use a polytope to 
> outer-approximate a convex set incrementally.  I want to use PPL to find 
> the set of new vertices of this polytope when I add a new halfspace to 
> it (ie. refine the outer-approximation).  The polytope lives in R^n, n 
> being a fixed dimension and I have no prior knowledge of the nature of 
> the coefficients, that is, they are not integer but just some floating 
> point numbers.
> 
> For this, my code looks like (given Ax <= b halfspace representation, 
> what new vertices are created by a new halfspace row'x <= rowB):
> 
>    PPL::Constraint_System cs;
>    for (unsigned int n = 0 ; n < b.size () ; ++n) {
>        PPL::Linear_Expression e;
> 
>        for (unsigned int j = 0 ; j < gmm::mat_ncols (A) ; ++j) {
>            if (j > 0)
>                std::cout << " + ";
> 
>            std::cout << A(n,j) << "*x_" << j;
>            e += mpq_class (A(n,j)) * PPL::Variable (j);
>        }
> 
>        cs.insert (e <= mpq_class (b[n]));
>        std::cout << " <= " << b[n] << std::endl;
>    }
> 
>    /* Add new row'x = rowB constraint
>     */
>    PPL::Linear_Expression e;
>    for (unsigned int j = 0 ; j < gmm::mat_ncols (A) ; ++j) {
>        if (j > 0)
>            std::cout << " + ";
> 
>        std::cout << row[j] << "*x_" << j;
>        e += mpf_class (row[j]) * PPL::Variable (j);
>    }
>    std::cout << " == " << rowB << std::endl;
>    cs.insert (e == mpf_class (rowB));
> 
>    using namespace PPL::IO_Operators;
>    std::cout << "Constraint system:\n" << cs << std::endl;
> 
> However, to my surprise, PPL seems to round the coefficients to integer 
> numbers (as printed by the constraint system).  I did not find any hint 
> in the documentation, except the mention of rational polytopes.  Is PPL 
> able to handle arbitrary floating point coefficients or do I always have 
> to normalize them in some way?  If so, what is the easiest way to do so?

Hello.

The coefficients of constraints and generators are integers
(mpz_class in the default configuration).
The domain of convex polyhedra implements rational polyhedra, but only 
allows for integer coefficients as far as the representation and 
input/output are concerned.

Any constraint having rational coefficients can be normalized so as to 
have integer coefficients (as for generators, the same can be done by 
exploiting the integer divisor, which by default is 1; but you seem to 
be most interested to the constraint representation).

The translation from a rational-coefficient constraint (i.e., your 
representation) to an equivalent integer-coefficient constraint (i.e., 
the PPL representation) is typically done by computing the lcm of all 
the denominators and normalizing the numerators accordingly.
To this end, you can use helper functions computing the lcm (or gcd) of 
two coefficients:

//! Assigns to \p x the least common multiple of \p y and \p z.
/*! \relates GMP_Integer */
void lcm_assign(GMP_Integer& x,
		const GMP_Integer& y, const GMP_Integer& z);


//! Assigns to \p x the greatest common divisor of \p y and \p z.
/*! \relates GMP_Integer */
void gcd_assign(GMP_Integer& x,
		const GMP_Integer& y, const GMP_Integer& z);

[NOTE: GMP_Integer is just an alias for mpz_class.]

I hope this is enough to bring you to the right track.
I am also forwarding this reply to the PPL-devel mailing list, since 
other users/developers may want to provide more details.

Cheers,
Enea Zaffanella.

> 
> Thanks,
> Sebastian Nowozin
> 
>