diff --git a/QP_solver/include/CGAL/QP_solver.h b/QP_solver/include/CGAL/QP_solver.h index 867fdae80b5..9544db32098 100644 --- a/QP_solver/include/CGAL/QP_solver.h +++ b/QP_solver/include/CGAL/QP_solver.h @@ -1178,6 +1178,7 @@ private: bool is_solution_optimal_for_auxiliary_problem(); bool is_solution_feasible(); bool is_solution_optimal(); + bool is_value_correct(); bool is_solution_unbounded(); public: diff --git a/QP_solver/include/CGAL/QP_solver/QP_solver.C b/QP_solver/include/CGAL/QP_solver/QP_solver.C index 1b864f452c4..c56096a5d73 100644 --- a/QP_solver/include/CGAL/QP_solver/QP_solver.C +++ b/QP_solver/include/CGAL/QP_solver/QP_solver.C @@ -89,54 +89,68 @@ typename QP_solver::ET QP_solver:: solution_numerator( ) const { - Index_const_iterator i_it; - Value_const_iterator x_i_it, c_it; ET s, z = et0; - int i; + int i, j; - // foreach i - x_i_it = x_B_O.begin(); + if (check_tag(Is_in_standard_form()) || is_phaseI) { + // standard form or phase I; it suffices to go + // through the basic variables; all D- and c-entries + // are obtained through the appropriate iterators + Index_const_iterator i_it; + Value_const_iterator x_i_it, c_it; + + // foreach i + x_i_it = x_B_O.begin(); c_it = minus_c_B .begin(); - for ( i_it = B_O.begin(); i_it != B_O.end(); ++i_it, ++x_i_it, ++c_it){ + for ( i_it = B_O.begin(); i_it != B_O.end(); ++i_it, ++x_i_it, ++c_it){ i = *i_it; // quadratic part s = et0; if ( is_QP && is_phaseII) { - // foreach j < i - s += std::inner_product(x_B_O.begin(), x_i_it, - D_pairwise_iterator( - B_O.begin(), - D_pairwise_accessor( qp_D, i)), - et0); + // foreach j < i + s += std::inner_product(x_B_O.begin(), x_i_it, + D_pairwise_iterator( + B_O.begin(), + D_pairwise_accessor( qp_D, i)), + et0); - // D_{i,i} x_i - s += ET( qp_D[ i][ i]) * *x_i_it; + // D_{i,i} x_i + s += ET( qp_D[ i][ i]) * *x_i_it; } // linear part s -= d * *c_it; // accumulate z += s * *x_i_it; - } - - // linear part of solution numerator with upper bounding - if (!check_tag(Is_in_standard_form())) { - // only in phaseII we need to account for nonbasic original variables, - // since in phaseI only artificial variables have nonzero objective - // function coefficients - // Since we do not have a minus_c_N vector we have to use - // qp_c - if (is_phaseII) { - s = et0; - for (int i = 0; i < qp_n; ++i) { - if (!is_basic(i)) { - s += ET(qp_c[i]) * nonbasic_original_variable_value(i); - } - } - z += d * d * s; - } + } + } else { + // nonstandard form and phase II, + // take all original variables into account; all D- and c-entries + // are obtained from the input data directly + // order in i_it and j_it matches original variable order + if (is_QP) { + // quadratic part + i=0; + for (Variable_numerator_iterator i_it = variables_numerator_begin(); + i_it < variables_numerator_end(); ++i_it, ++i) { + // do something only if *i_it != 0 + if (*i_it == et0) continue; + s = et0; // contribution of i-th row + j=0; + for (Variable_numerator_iterator j_it = variables_numerator_begin(); + j_it < variables_numerator_end(); ++j_it, ++j) + s += ET(qp_D[i][j]) * *j_it; + z += s * *i_it; + } + } + // linear part + j=0; s = et0; + for (Variable_numerator_iterator j_it = variables_numerator_begin(); + j_it < variables_numerator_end(); ++j_it, ++j) + s += ET(qp_c[j]) * *j_it; + z += d * s; } return z; } diff --git a/QP_solver/include/CGAL/QP_solver/Validity.C b/QP_solver/include/CGAL/QP_solver/Validity.C index 8660ab4c675..e00b2b17d3c 100644 --- a/QP_solver/include/CGAL/QP_solver/Validity.C +++ b/QP_solver/include/CGAL/QP_solver/Validity.C @@ -41,6 +41,7 @@ bool QP_solver::is_valid() { const bool f = this->is_solution_feasible(); const bool o = this->is_solution_optimal(); + const bool v = this->is_value_correct(); CGAL_qpe_debug { vout << std::endl << "----------" << std::endl @@ -49,8 +50,9 @@ bool QP_solver::is_valid() vout << " is in phase II: " << is_phaseII << std::endl; vout << " feasible: " << f << std::endl; vout << " optimal: " << o << std::endl; + vout << " correct value: " << v << std::endl; } - return is_phaseII && f && o; + return is_phaseII && f && o && v; } case INFEASIBLE: { @@ -174,6 +176,36 @@ bool QP_solver::is_solution_feasible_for_auxiliary_problem() return true; } +template < class Rep_ > +bool QP_solver::is_value_correct() +{ + // checks whether solution_numerator() returns the right value + // by computing x^T D x + x^T c from scratch; we use the numerators + // of the variables, so x^T D x carries a factor of d^2, while x^T c + // carries a factor of d (must be compensated for); solution_numerator() + // carries a factor of d^2 + CGAL_qpe_assertion(is_phaseII); + // compute x^T D x + d c = sum_i x_i (D_i x + d c_i) + ET z = et0; // for x^T D x + d c + int i = 0; + for (Variable_numerator_iterator i_it = variables_numerator_begin(); + i_it < variables_numerator_end(); ++i_it, ++i) { + if (*i_it == et0) continue; // no contribution from this variable + ET s = et0; // for D_i x + c_i + int j = 0; + if (is_QP) { + for (Variable_numerator_iterator j_it = variables_numerator_begin(); + j_it < variables_numerator_end(); ++j_it, ++j) + s += ET(qp_D[i][j]) * *j_it; // D_ij x_j + } + // now s = D_i x + s += d * qp_c[i]; // add d * c_i + // now s = D_i x + d c_i + z += s * *i_it; // add x_i (D_i x + d c_i) + } + return z == solution_numerator(); +} + template < class Rep_ > bool QP_solver::is_solution_optimal_for_auxiliary_problem() { diff --git a/QP_solver/test/QP_solver/test_solver.C b/QP_solver/test/QP_solver/test_solver.C index 868ef293449..4fb2b4b7caa 100644 --- a/QP_solver/test/QP_solver/test_solver.C +++ b/QP_solver/test/QP_solver/test_solver.C @@ -426,7 +426,6 @@ bool process(const std::string& filename, if (!check_tag(Is_in_standard_form()) && options.find("Strategy")->second != FE) return true; - // solve: CGAL::QP_pricing_strategy *s = create_strategy(options); CGAL::QP_solver solver(qp.number_of_variables(),