mirror of https://github.com/CGAL/cgal
- added a farkas-Lemma proof of infeasibility
- added test with random small QP's
- fixed some bugs discovered by this:
- output routine for programs is now printing symmetric D
- internal comparison routine is now ignoring bound value if it
was specified as infinite
- fixed sign error in subrotuine of ratio_test
- removed constant term from consideration in phase I
This commit is contained in:
Bernd Gärtner
parent
4910117ccb
commit
9f5e838e6f
|
|
@ -1340,7 +1340,8 @@ private:
|
|||
// validity checks:
|
||||
bool is_solution_feasible_for_auxiliary_problem() const;
|
||||
bool is_solution_optimal_for_auxiliary_problem() const;
|
||||
bool is_solution_feasible() const;
|
||||
bool is_solution_feasible() const;
|
||||
bool is_solution_infeasible() const;
|
||||
bool is_solution_optimal() const;
|
||||
bool is_value_correct() const;
|
||||
bool is_solution_unbounded() const;
|
||||
|
|
|
|||
|
|
@ -195,7 +195,7 @@ init_x_O_v_i()
|
|||
// and so we initialize them to zero (if the bound on the variable
|
||||
// allows it), or to the variable's lower or upper bound:
|
||||
for (int i = 0; i < qp_n; ++i) {
|
||||
CGAL_qpe_assertion(qp_l[i]<=qp_u[i] || !*(qp_fl+i) || !*(qp_fu+i));
|
||||
CGAL_qpe_assertion( !*(qp_fl+i) || !*(qp_fu+i) || qp_l[i]<=qp_u[i]);
|
||||
|
||||
if (*(qp_fl+i)) // finite lower bound?
|
||||
if (*(qp_fu+i)) // finite lower and finite upper bound?
|
||||
|
|
|
|||
|
|
@ -135,7 +135,8 @@ namespace QP_functions_detail {
|
|||
int n = p.n();
|
||||
bool empty_D = true;
|
||||
for (int i=0; i<n; ++i, ++it) {
|
||||
for (int j=0; j<n; ++j)
|
||||
// handle only entries on/below diagonal
|
||||
for (int j=0; j<i+1; ++j)
|
||||
if (!CGAL::is_zero(*(*it + j))) {
|
||||
if (empty_D) {
|
||||
// first time we see a nonzero entry
|
||||
|
|
@ -143,6 +144,9 @@ namespace QP_functions_detail {
|
|||
empty_D = false;
|
||||
}
|
||||
out << " x" << i << " x" << j << " " << *(*it + j) << "\n";
|
||||
// QMATRIX format prescribes symmetric matrix
|
||||
if (i != j)
|
||||
out << " x" << j << " x" << i << " " << *(*it + j) << "\n";
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -198,42 +202,41 @@ namespace QP_functions_detail {
|
|||
<< *r1 << " vs. " << *r2 << std::endl;
|
||||
return_val = false;
|
||||
}
|
||||
// check fl
|
||||
// check fl, l
|
||||
typename Quadratic_program1::FL_iterator fl1 = qp1.fl();
|
||||
typename Quadratic_program2::FL_iterator fl2 = qp2.fl();
|
||||
for (int j=0; j<n; ++j, ++fl1, ++fl2)
|
||||
typename Quadratic_program1::L_iterator l1 = qp1.l();
|
||||
typename Quadratic_program2::L_iterator l2 = qp2.l();
|
||||
for (int j=0; j<n; ++j, ++fl1, ++fl2, ++l1, ++l2) {
|
||||
if (*fl1 != *fl2) {
|
||||
std::cerr << "Equality test fails with fl[" << j << "]: "
|
||||
<< *fl1 << " vs. " << *fl2 << std::endl;
|
||||
return_val = false;
|
||||
}
|
||||
// check l
|
||||
typename Quadratic_program1::L_iterator l1 = qp1.l();
|
||||
typename Quadratic_program2::L_iterator l2 = qp2.l();
|
||||
for (int j=0; j<n; ++j, ++l1, ++l2)
|
||||
if (*l1 != *l2) {
|
||||
if ((*fl1 == true) && (*l1 != *l2)) {
|
||||
std::cerr << "Equality test fails with l[" << j << "]: "
|
||||
<< *l1 << " vs. " << *l2 << std::endl;
|
||||
return_val = false;
|
||||
}
|
||||
// check fu
|
||||
}
|
||||
|
||||
// check fu, u
|
||||
typename Quadratic_program1::FU_iterator fu1 = qp1.fu();
|
||||
typename Quadratic_program2::FU_iterator fu2 = qp2.fu();
|
||||
for (int j=0; j<n; ++j, ++fu1, ++fu2)
|
||||
typename Quadratic_program1::U_iterator u1 = qp1.u();
|
||||
typename Quadratic_program2::U_iterator u2 = qp2.u();
|
||||
for (int j=0; j<n; ++j, ++fu1, ++fu2, ++u1, ++u2) {
|
||||
if (*fu1 != *fu2) {
|
||||
std::cerr << "Equality test fails with fu[" << j << "]: "
|
||||
<< *fu1 << " vs. " << *fu2 << std::endl;
|
||||
return_val = false;
|
||||
}
|
||||
// check u
|
||||
typename Quadratic_program1::U_iterator u1 = qp1.u();
|
||||
typename Quadratic_program2::U_iterator u2 = qp2.u();
|
||||
for (int j=0; j<n; ++j, ++u1, ++u2)
|
||||
if (*u1 != *u2) {
|
||||
if ((*fu1 == true) && (*u1 != *u2)) {
|
||||
std::cerr << "Equality test fails with u[" << j << "]: "
|
||||
<< *u1 << " vs. " << *u2 << std::endl;
|
||||
return_val = false;
|
||||
}
|
||||
}
|
||||
// check d
|
||||
typename Quadratic_program1::D_iterator d1 = qp1.d();
|
||||
typename Quadratic_program2::D_iterator d2 = qp2.d();
|
||||
|
|
@ -299,19 +302,28 @@ namespace QP_functions_detail {
|
|||
typename P::C_iterator c = p.c();
|
||||
out << "COLUMNS\n";
|
||||
for (int j=0; j<n; ++j, ++c, ++a) {
|
||||
if (!CGAL_NTS is_zero (*c))
|
||||
// make sure that variable appears here even if it has only
|
||||
// zero coefficients
|
||||
bool written = false;
|
||||
if (!CGAL_NTS is_zero (*c)) {
|
||||
out << " x" << j << " obj " << *c << "\n";
|
||||
for (int i=0; i<m; ++i) {
|
||||
if (!CGAL_NTS is_zero ((*a)[i]))
|
||||
out << " x" << j << " c" << i << " " << (*a)[i] << "\n";
|
||||
written = true;
|
||||
}
|
||||
for (int i=0; i<m; ++i) {
|
||||
if (!CGAL_NTS is_zero ((*a)[i])) {
|
||||
out << " x" << j << " c" << i << " " << (*a)[i] << "\n";
|
||||
written = true;
|
||||
}
|
||||
}
|
||||
if (!written)
|
||||
out << " x" << j << " obj " << 0 << "\n";
|
||||
}
|
||||
|
||||
// RHS section:
|
||||
typename P::B_iterator b = p.b();
|
||||
out << "RHS\n";
|
||||
if (!CGAL_NTS is_zero (p.c0()))
|
||||
out << " rhs obj -" << p.c0() << "\n";
|
||||
out << " rhs obj " << -p.c0() << "\n";
|
||||
for (int i=0; i<m; ++i, ++b)
|
||||
if (!CGAL_NTS is_zero (*b))
|
||||
out << " rhs c" << i << " " << *b << "\n";
|
||||
|
|
|
|||
|
|
@ -161,8 +161,9 @@ solution_numerator( ) const
|
|||
s += et2 * ET(qp_c[j]) * *j_it;
|
||||
z += d * s;
|
||||
}
|
||||
// finally, add the constant term
|
||||
return z += et2 * ET(qp_c0) * d * d;
|
||||
// finally, add the constant term (phase II only)
|
||||
if (is_phaseII) z += et2 * ET(qp_c0) * d * d;
|
||||
return z;
|
||||
}
|
||||
|
||||
// pivot step
|
||||
|
|
@ -929,7 +930,7 @@ test_mixed_bounds_dir_neg(int k, const ET& x_k, const ET& q_k,
|
|||
}
|
||||
}
|
||||
} else { // check for upper bound
|
||||
if ((q_k < et0) && (k < qp_n) && *(qp_fu+k)) {
|
||||
if ((q_k > et0) && (k < qp_n) && *(qp_fu+k)) {
|
||||
ET diff = (d * ET(qp_u[k])) - x_k;
|
||||
if (diff * q_min < d_min * q_k) {
|
||||
i_min = k;
|
||||
|
|
|
|||
|
|
@ -75,6 +75,7 @@ bool QP_solver<Q, ET, Tags>::is_valid() const
|
|||
{
|
||||
const bool f = this->is_solution_feasible_for_auxiliary_problem();
|
||||
const bool o = this->is_solution_optimal_for_auxiliary_problem();
|
||||
const bool i = this->is_solution_infeasible();
|
||||
const bool aux_positive =
|
||||
( (this->solution_numerator() > et0) && (d > et0) );
|
||||
CGAL_qpe_debug {
|
||||
|
|
@ -85,9 +86,10 @@ bool QP_solver<Q, ET, Tags>::is_valid() const
|
|||
vout << " is in phase I: " << is_phaseI << std::endl;
|
||||
vout << " feasible_aux: " << f << std::endl;
|
||||
vout << " optimal_aux: " << o << std::endl;
|
||||
vout << " infeasible: " << i << std::endl;
|
||||
vout << " obj. val. positive: " << aux_positive << std::endl;
|
||||
}
|
||||
return is_phaseI && f && o && aux_positive;
|
||||
return is_phaseI && f && o && i && aux_positive;
|
||||
}
|
||||
case QP_UNBOUNDED:
|
||||
{
|
||||
|
|
@ -127,16 +129,30 @@ bool QP_solver<Q, ET, Tags>::is_solution_feasible_for_auxiliary_problem() const
|
|||
if (*i_it < qp_n) { // original variable?
|
||||
if ((has_finite_lower_bound(*i_it) && (*v_it < lower_bound(*i_it) * d))||
|
||||
(has_finite_upper_bound(*i_it) && (*v_it > upper_bound(*i_it) * d)))
|
||||
return false;
|
||||
{
|
||||
CGAL_qpe_debug {
|
||||
vout4 << "variable " << *i_it << " out of bounds" << std::endl;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
} else // artificial variable?
|
||||
if (*v_it < et0)
|
||||
return false;
|
||||
{
|
||||
CGAL_qpe_debug {
|
||||
vout4 << "artificial variable negative" << std::endl;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
// check nonegativity of slack variables (the basic ones suffice):
|
||||
for (Value_const_iterator v_it = x_B_S.begin(); v_it != x_B_S.end(); ++v_it)
|
||||
if (*v_it < et0)
|
||||
if (*v_it < et0) {
|
||||
CGAL_qpe_debug {
|
||||
vout4 << "slack variable out of bounds" << std::endl;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
// Check whether the current solution is feasible for the auxiliary
|
||||
// problem. For this, we use the fact that the auxiliary problem
|
||||
|
|
@ -192,11 +208,75 @@ bool QP_solver<Q, ET, Tags>::is_solution_feasible_for_auxiliary_problem() const
|
|||
|
||||
// check equality (C1);
|
||||
for (int i=0; i<qp_m; ++i)
|
||||
if (lhs_col[i] != ET(qp_b[i]) * d)
|
||||
if (lhs_col[i] != ET(qp_b[i]) * d) {
|
||||
CGAL_qpe_debug {
|
||||
vout4 << "row " << i << " infeasible" << std::endl;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template < typename Q, typename ET, typename Tags >
|
||||
bool QP_solver<Q, ET, Tags>::is_solution_infeasible() const
|
||||
{
|
||||
// checks whether we have a proof of infeasibilty according
|
||||
// to Farkas Lemma. Namely, the system is infeasible iff
|
||||
// \tau^T = \lambda^T A (A)
|
||||
// \lambda^T b - \sum{j: \tau_j < 0} \tau_j u_j
|
||||
// - \sum{j: \tau_j > 0} \tau_j l_j < 0 (B)
|
||||
// \lambda_i >= 0 for <=-constraints i, (C)
|
||||
// \lambda_i <= 0 for >=-constraints i, (D)
|
||||
|
||||
// get \lambda:
|
||||
// Note: as the \lambda' we have access to is actual d times the \lambda
|
||||
// we will not compute \tau but \tau * d.
|
||||
Optimality_certificate_numerator_iterator
|
||||
itb = optimality_certificate_numerator_begin();
|
||||
Optimality_certificate_numerator_iterator
|
||||
ite = optimality_certificate_numerator_end();
|
||||
Optimality_certificate_iterator
|
||||
itq = optimality_certificate_begin();
|
||||
Values lambda_prime;
|
||||
int row = 0;
|
||||
for (Optimality_certificate_numerator_iterator
|
||||
it = itb; it != ite; ++it, ++itq, ++row) {
|
||||
// simply check whether numerator/denominator gives correct value
|
||||
CGAL_qpe_assertion (*it * (*itq).denominator() == (*itq).numerator() * d);
|
||||
lambda_prime.push_back(*it);
|
||||
// check sign (C) (D)
|
||||
if (qp_r[row] == CGAL::SMALLER)
|
||||
if (lambda_prime[row] < et0) return false;
|
||||
if (qp_r[row] == CGAL::LARGER)
|
||||
if (lambda_prime[row] > et0) return false;
|
||||
}
|
||||
|
||||
// compute \tau^T = \lambda^T A (A):
|
||||
Values tau(qp_n,et0);
|
||||
for (int col = 0; col < qp_n; ++col)
|
||||
for (int i=0; i<qp_m; ++i)
|
||||
tau[col] += lambda_prime[i] * ET((*(qp_A+col))[i]);
|
||||
|
||||
// check condition (B)
|
||||
ET sum = et0;
|
||||
// add \lambda^T b
|
||||
for (int i=0; i<qp_m; ++i)
|
||||
sum += lambda_prime[i] * ET(qp_b[i]); // lambda carries factor of d
|
||||
// subtract the remaining terms
|
||||
for (int col = 0; col < qp_n; ++col) {
|
||||
if (tau[col] < et0) {
|
||||
CGAL_qpe_assertion(has_finite_upper_bound(col));
|
||||
sum -= tau[col] * upper_bound(col); // tau carries factor of d
|
||||
}
|
||||
if (tau[col] > et0) {
|
||||
CGAL_qpe_assertion(has_finite_lower_bound(col));
|
||||
sum -= tau[col] * lower_bound(col); // tau carries factor of d
|
||||
}
|
||||
}
|
||||
return sum < et0;
|
||||
}
|
||||
|
||||
|
||||
template < typename Q, typename ET, typename Tags >
|
||||
bool QP_solver<Q, ET, Tags>::is_value_correct() const
|
||||
{
|
||||
|
|
@ -443,8 +523,12 @@ bool QP_solver<Q, ET, Tags>::is_solution_feasible() const
|
|||
if (var != *it)
|
||||
return false; // original_variables_numerator_begin() inconsistent
|
||||
if ((has_finite_lower_bound(i) && var < lower_bound(i) * d) ||
|
||||
(has_finite_upper_bound(i) && var > upper_bound(i) * d))
|
||||
(has_finite_upper_bound(i) && var > upper_bound(i) * d)) {
|
||||
CGAL_qpe_debug {
|
||||
vout4 << "variable " << i << " out of bounds" << std::endl;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
// compute A x times d:
|
||||
for (int j=0; j<qp_m; ++j)
|
||||
|
|
@ -456,8 +540,12 @@ bool QP_solver<Q, ET, Tags>::is_solution_feasible() const
|
|||
const ET rhs = ET(qp_b[row])*d;
|
||||
if ((qp_r[row] == CGAL::EQUAL && lhs_col[row] != rhs) ||
|
||||
(qp_r[row] == CGAL::SMALLER && lhs_col[row] > rhs) ||
|
||||
(qp_r[row] == CGAL::LARGER && lhs_col[row] < rhs))
|
||||
(qp_r[row] == CGAL::LARGER && lhs_col[row] < rhs)) {
|
||||
CGAL_qpe_debug {
|
||||
vout4 << "row " << row << " infeasible" << std::endl;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
|
|
|
|||
|
|
@ -191,7 +191,7 @@ void create_shifted_instance(CGAL::QP_from_mps
|
|||
for (int i=0; i<n; ++i) {
|
||||
for (int j=0; j<n; ++j)
|
||||
mvTD[i] += qp.d()[j][i] * v[j];
|
||||
mvTD[i] *= -2;
|
||||
mvTD[i] *= -1;
|
||||
}
|
||||
|
||||
// output:
|
||||
|
|
|
|||
|
|
@ -25,7 +25,7 @@
|
|||
#ifndef CGAL_USE_GMP
|
||||
#include <CGAL/MP_Float.h>
|
||||
#else
|
||||
#include <CGAL/Gmpzf.h>
|
||||
#include <CGAL/Gmpz.h>
|
||||
#endif
|
||||
|
||||
#include <CGAL/QP_solver.h>
|
||||
|
|
@ -49,11 +49,11 @@ int main(const int argNr,const char **args) {
|
|||
const int verbosity = argNr < 2? 1 : std::atoi(args[1]);
|
||||
|
||||
// construct QP instance:
|
||||
typedef double IT;
|
||||
typedef int IT;
|
||||
#ifndef CGAL_USE_GMP
|
||||
typedef CGAL::MP_Float ET;
|
||||
#else
|
||||
typedef CGAL::Gmpzf ET;
|
||||
typedef CGAL::Gmpz ET;
|
||||
#endif
|
||||
typedef CGAL::QP_from_mps<IT, CGAL::Tag_false, CGAL::Tag_true> QP;
|
||||
QP qp(std::cin,true,verbosity);
|
||||
|
|
@ -69,7 +69,7 @@ int main(const int argNr,const char **args) {
|
|||
CGAL::print_quadratic_program (cout, qp);
|
||||
cout << std::endl;
|
||||
}
|
||||
typedef CGAL::QP_solver_impl::QP_tags<CGAL::Tag_false,CGAL::Tag_true> Tags;
|
||||
typedef CGAL::QP_solver_impl::QP_tags<CGAL::Tag_false,CGAL::Tag_false> Tags;
|
||||
CGAL::QP_pricing_strategy<QP, ET, Tags> *pricing_strategy =
|
||||
new CGAL::QP_full_exact_pricing<QP, ET, Tags>;
|
||||
typedef CGAL::QP_solver<QP, ET, Tags> Solver;
|
||||
|
|
|
|||
Loading…
Reference in New Issue