release: 1.1.29

[lilypond.git] / lily / ineq-constrained-qp.cc

diff --git a/lily/ineq-constrained-qp.cc b/lily/ineq-constrained-qp.cc
index f1cfebfdd2089a239ebeb34e8d7658234561b8b8..2863799a5984600b020ef3da8cada051de0689a6 100644 (file)
--- a/lily/ineq-constrained-qp.cc
+++ b/lily/ineq-constrained-qp.cc
@@ -3,7 +3,7 @@
 
   source file of the GNU LilyPond music typesetter
 
-  (c) 1996,1997 Han-Wen Nienhuys <hanwen@stack.nl>
+  (c) 1996, 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
 */
 #include "ineq-constrained-qp.hh"
 #include "qlpsolve.hh"
@@ -12,11 +12,12 @@
 #include "main.hh"
 
 /*
-  MAy be this also should go into a library
+  May be this also should go into a library
  */
 
 const int MAXITER=100;         // qlpsolve.hh
 
+const int MAXDEGEN=5;
 
 /*
   assume x (idx) == value, and adjust constraints, lin and quad accordingly
@@ -113,32 +114,26 @@ Ineq_constrained_qp::constraint_solve (Vector start) const
   if (!dim())
     return Vector (0);
 
-  // experimental
-  if (quad_.dim() > 10)
-    quad_.try_set_band();
-
   Active_constraints act (this);
   act.OK();
 
 
   Vector x (start);
   Vector gradient=quad_*x+lin_;
-  //    Real fvalue = x*quad_*x/2 + lin*x + const_term;
-  // it's no use.
 
+  
+  //    Real fvalue = x*quad_*x/2 + lin*x + const_term;// it's no use.
   Vector last_gradient (gradient);
   int iterations=0;
 
-  while (iterations++ < MAXITER)
+  while (iterations++ < MAXITER && act.degenerate_count_i_ < MAXDEGEN)
     {
-      //#ifdef PARANOID
       if (experimental_features_global_b)
        assert_solution (x);
-      //#endif
       
       Vector direction= - act.find_active_optimum (gradient);
 
-      DOUT << "gradient "<< gradient<< "\ndirection " << direction<<"\n";
+      DOUT << "gradient "<< gradient<< "\ndirection " << direction<< '\n';
 
       if (direction.norm() > EPS)
        {
@@ -158,15 +153,22 @@ Ineq_constrained_qp::constraint_solve (Vector start) const
          for (Inactive_iter ia (act); ia.ok(); ia++)
            {
              Real dot = ia.vec() * direction;
-             if (dot >= 0)
+             Real mindot =  (experimental_features_global_b)
+               ? -EPS
+               : 0;
+             
+             if (dot >= mindot)
                continue;
-
+             
              
              Real numerator = ia.rhs () - ia.vec()*x;
              if (numerator >= 0)
                {
                  if (numerator > EPS)
-                   warning (String ("Ineq_constrained_qp::solve (): Constraint off by ") + numerator);
+                   {
+                     warning (_f ("Ineq_constrained_qp::solve (): Constraint off by %f", numerator));
+                     act.degenerate_count_i_ ++;
+                   }
                  minalf = -numerator;
                  minidx = ia;
                  break;
@@ -184,23 +186,21 @@ Ineq_constrained_qp::constraint_solve (Vector start) const
 
          Real optimal_step = minalf;
 
-         Vector deltax=direction * optimal_step;
+         Vector deltax = direction * optimal_step;
          x += deltax;
          gradient += optimal_step * (quad_ * deltax);
 
-         DOUT << "step = " << optimal_step<< " (|dx| = " <<
-           deltax.norm() << ")\n";
+         DOUT << "step = " << optimal_step << " (|dx| = " <<
+           to_str (deltax.norm()) << ")\n";
 
          if (minalf < unbounded_alfa)
            {
              /* bumped into an edge. try again, in smaller space. */
-             act.add (minidx.idx());
-             DOUT << "adding cons "<< minidx.idx()<<'\n';
+             act.add_constraint (minidx.idx());
+             DOUT << "adding cons "<< minidx.idx () << '\n';
              continue;
            }
-         /*ASSERT: we are at optimal solution for this "plane"*/
-
-
+         /*ASSERT: we are at the optimal solution for this "plane"*/
        }
 
       Vector lagrange_mult=act.get_lagrange (gradient);
@@ -217,13 +217,14 @@ Ineq_constrained_qp::constraint_solve (Vector start) const
          break;
        }
 
-      DOUT << "dropping cons " << m<<'\n';
-      act.drop (m);
+      DOUT << "dropping cons " << m << '\n';
+      act.drop_constraint (m);
     }
   if (iterations >= MAXITER)
-    WARN<<_("didn't converge!\n");
-
-  DOUT <<  ": found " << x<<" in " << iterations <<" iterations\n";
+    WARN << _ ("didn't converge!") << '\n';
+  if (act.degenerate_count_i_ >= MAXDEGEN)
+    WARN << _ ("Too much degeneracy. ") << '\n';
+  DOUT <<  ": found " << x << " in " << iterations <<" iterations\n";
   assert_solution (x);
   return x;
 }
@@ -250,3 +251,33 @@ Ineq_constrained_qp::dim () const
   return lin_.dim();
 }
 
+
+
+
+void
+Ineq_constrained_qp::assert_solution (Vector sol) const
+{
+  bool sol_b=true;
+
+  for (int i=0; sol_b && i < cons_.size(); i++) 
+    {
+      Real R=cons_[i] * sol- consrhs_[i];
+      if (R> -EPS)
+       sol_b = false;
+    }
+}
+
+void
+Ineq_constrained_qp::print() const
+{
+#ifndef NPRINT
+  DOUT << "Quad " << quad_;
+  DOUT << "lin " << lin_ << '\n'
+       << "const " << const_term_<< '\n';
+  for (int i=0; i < cons_.size(); i++) 
+    {
+      DOUT << "constraint["<<i<<"]: " << cons_[i] << " >= " << consrhs_[i];
+      DOUT << '\n';
+    }
+#endif
+}