fixed bug that led to wrong objective function values

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:

QP_solver/include/CGAL/QP_solver/QP_solver.C

@@ -89,54 +89,68 @@ typename QP_solver<Rep_>::ET
 QP_solver<Rep_>::
 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;
 }

QP_solver/include/CGAL/QP_solver/Validity.C

@@ -41,6 +41,7 @@ bool QP_solver<Rep_>::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<Rep_>::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<Rep_>::is_solution_feasible_for_auxiliary_problem()
   return true;
 }
 
+template < class Rep_ >
+bool QP_solver<Rep_>::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<Rep_>::is_solution_optimal_for_auxiliary_problem()
 {

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<Traits> *s = create_strategy<Traits>(options);
   CGAL::QP_solver<Traits> solver(qp.number_of_variables(),