patch::: 1.2.9.jcn2

[lilypond.git] / lily / qlpsolve.cc

diff --git a/lily/qlpsolve.cc b/lily/qlpsolve.cc
index d2ffad85d43115d539fdd57dd7b96905bb360e84..3d64b3779289b4ccd58d559e1caa7375784483e9 100644 (file)
--- a/lily/qlpsolve.cc
+++ b/lily/qlpsolve.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>
 
   TODO:
   try fixed point arithmetic, to speed up lily.
@@ -20,15 +20,15 @@ String
 Active_constraints::status() const
 {
   String s ("Active|Inactive [");
-  for (int i=0; i< active.size(); i++) 
+  for (int i=0; i< active.size(); i++)
     {
-       s += String (active[i]) + " ";
+      s += to_str (active[i]) + " ";
     }
 
   s+="| ";
-  for (int i=0; i< inactive.size(); i++) 
+  for (int i=0; i< inactive.size(); i++)
     {
-       s += String (inactive[i]) + " ";
+      s += to_str (inactive[i]) + " ";
     }
   s+="]";
 
@@ -41,25 +41,25 @@ Active_constraints::OK()
 #ifndef NDEBUG
   H.OK();
   A.OK();
-  assert (active.size() +inactive.size () == opt->cons.size ());
+  assert (active.size() +inactive.size () == opt->cons_.size ());
   assert (H.dim() == opt->dim ());
   assert (active.size() == A.rows ());
   Array<int> allcons;
 
-  for (int i=0; i < opt->cons.size(); i++)
-       allcons.push (0);
-  for (int i=0; i < active.size(); i++) 
+  for (int i=0; i < opt->cons_.size(); i++)
+    allcons.push (0);
+  for (int i=0; i < active.size(); i++)
     {
-       int j = active[i];
-       allcons[j]++;
+      int j = active[i];
+      allcons[j]++;
     }
-  for (int i=0; i < inactive.size(); i++) 
+  for (int i=0; i < inactive.size(); i++)
     {
-       int j = inactive[i];
-       allcons[j]++;
+      int j = inactive[i];
+      allcons[j]++;
     }
   for (int i=0; i < allcons.size(); i++)
-       assert (allcons[i] == 1);
+    assert (allcons[i] == 1);
 #endif
 }
 
@@ -70,96 +70,100 @@ Active_constraints::get_lagrange (Vector gradient)
 }
 
 void
-Active_constraints::add (int k)
+Active_constraints::add_constraint (int k)
 {
   // add indices
   int cidx=inactive[k];
-  active.push (cidx);
 
-  inactive.swap (k,inactive.size()-1);
-  inactive.pop();
-
-  Vector a (opt->cons[cidx]);
+  Vector a (opt->cons_[cidx]);
   // update of matrices
   Vector Ha = H*a;
   Real aHa = a*Ha;
   Vector addrow (Ha.dim());
-  if (abs (aHa) > EPS) 
+  bool degenerate =  (abs (aHa) < EPS);
+
+  if (degenerate)
+    {
+      warning (String ("Active_constraints::add ():")
+       + _("degenerate constraints"));
+      DOUT << "Ha = "  << Ha.str () << '\n';
+      /*
+       a != 0, so if Ha = O(EPS), then
+       Ha * aH / aHa = O(EPS^2/EPS)
+
+       if H*a == 0, the constraints are dependent.
+      */
+      degenerate_count_i_ ++;
+    }
+  if (!degenerate)
     {
-       /*
-         a != 0, so if Ha = O(EPS), then
-         Ha * aH / aHa = O(EPS^2/EPS)
-
-         if H*a == 0, the constraints are dependent.
-         */
-       H -= Matrix (Ha/aHa , Ha);
-  
-
-       /*
-         sorry, don't know how to justify this. ..
-         */
-       addrow=Ha;
-      addrow/= aHa;
-       A -= Matrix (A*a, addrow);
-       A.insert_row (addrow,A.rows());
-  }else
-       WARN << "degenerate constraints";
+      active.push (cidx);
+      inactive.swap (k,inactive.size()-1);
+      inactive.pop();
+
+      H -= Matrix (Ha/aHa , Ha);
+      
+      addrow=Ha;
+      addrow /= aHa;
+      A -= Matrix (A*a, addrow);
+      A.insert_row (addrow,A.rows());
+    }
 }
 
 void
-Active_constraints::drop (int k)
+Active_constraints::drop_constraint (int k)
 {
   int q=active.size()-1;
 
-      // drop indices
-  inactive.push (active[k]);
-  active.swap (k,q);
-  A.swap_rows (k,q);
-  active.pop();
 
   Vector a (A.row (q));
-  if (a.norm() > EPS) 
+  if (a.norm() > EPS)
     {
-       /*
-        
-        */
-      Real q = a*opt->quad*a;
-       Matrix aaq (a,a/q);
-       H += aaq;
-       A -= A*opt->quad*aaq;
-  }else
-       WARN << "degenerate constraints";
-#ifndef NDEBUG
-  Vector rem_row (A.row (q));
-  assert (rem_row.norm() < EPS);
-#endif
-   
-  A.delete_row (q);
+      // drop indices
+      inactive.push (active[k]);
+      active.swap (k,q);
+      A.swap_rows (k,q);
+      active.pop();
+      /*
+
+       */
+      Real aqa = a*opt->quad_*a;
+      Matrix aaq (a,a/aqa);
+      H += aaq;
+      A -= A*opt->quad_*aaq;
+ 
+      A.delete_row (q);
+    }else {
+      degenerate_count_i_ ++;
+      warning (String ("Active_constraints::drop ():")
+       + _("degenerate constraints"));
+    }
 }
 
 
 Active_constraints::Active_constraints (Ineq_constrained_qp const *op)
-  :       A(0,op->dim()),
-           H(op->dim()),
-           opt (op)
+  : A(0,op->dim()),
+    H(op->dim()),
+    opt (op)
 {
-  for (int i=0; i < op->cons.size(); i++)
-       inactive.push (i);
-  Choleski_decomposition chol (op->quad);
+  for (int i=0; i < op->cons_.size(); i++)
+    inactive.push (i);
+  Choleski_decomposition chol (op->quad_);
 
   /*
     ugh.
-   */
+    */
   H=chol.inverse();
   OK();
+
+  degenerate_count_i_ = 0;
 }
 
 /** Find the optimum which is in the planes generated by the active
-  constraints.        
+  constraints.
   */
 Vector
 Active_constraints::find_active_optimum (Vector g)
 {
   return H*g;
 }
-