release: 0.1.24

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

diff --git a/lily/ineq-constrained-qp.cc b/lily/ineq-constrained-qp.cc
index 4455a195b4590b23ee1128352e8e7bb8d41b8313..f14697dbb6792d59dcedf6db7c67053eb5b398b5 100644 (file)
--- a/lily/ineq-constrained-qp.cc
+++ b/lily/ineq-constrained-qp.cc
@@ -35,7 +35,7 @@ Ineq_constrained_qp::eliminate_var (int idx, Real value)
   row.del (idx);
   lin +=row ;
 
-  for (int i=0; i < cons.size(); i++) 
+  for (int i=0; i < cons.size(); i++)
     {
       consrhs[i] -= cons[i](idx) *value;
       cons[i].del (idx);
@@ -63,9 +63,9 @@ Ineq_constrained_qp::OK() const
   assert (cons.size() == consrhs.size ());
   Matrix Qdif= quad - quad.transposed();
   assert (Qdif.norm()/quad.norm () < EPS);
-#endif    
+#endif
 }
-   
+
 
 Real
 Ineq_constrained_qp::eval (Vector v)
@@ -77,11 +77,11 @@ Ineq_constrained_qp::eval (Vector v)
 int
 min_elt_index (Vector v)
 {
-  Real m=infinity_f; 
+  Real m=infinity_f;
   int idx=-1;
   for (int i = 0; i < v.dim(); i++)
     {
-      if (v (i) < m) 
+      if (v (i) < m)
        {
          idx = i;
          m = v (i);
@@ -106,19 +106,19 @@ min_elt_index (Vector v)
 
   */
 Vector
-Ineq_constrained_qp::constraint_solve (Vector start) const 
-{    
+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();    
+  act.OK();
+
 
-  
   Vector x (start);
   Vector gradient=quad*x+lin;
   //    Real fvalue = x*quad*x/2 + lin*x + const_term;
@@ -126,17 +126,17 @@ Ineq_constrained_qp::constraint_solve (Vector start) const
 
   Vector last_gradient (gradient);
   int iterations=0;
-  
-  while (iterations++ < MAXITER) 
+
+  while (iterations++ < MAXITER)
     {
       Vector direction= - act.find_active_optimum (gradient);
-       
+
       DOUT << "gradient "<< gradient<< "\ndirection " << direction<<"\n";
-       
-      if (direction.norm() > EPS) 
+
+      if (direction.norm() > EPS)
        {
          DOUT << act.status() << '\n';
-           
+
          Real minalf = infinity_f;
 
          Inactive_iter minidx (act);
@@ -146,16 +146,16 @@ Ineq_constrained_qp::constraint_solve (Vector start) const
            we know the optimum on this "hyperplane". Check if we
            bump into the edges of the simplex
            */
-  
-         for (Inactive_iter ia (act); ia.ok(); ia++) 
+
+         for (Inactive_iter ia (act); ia.ok(); ia++)
            {
 
              if (ia.vec() * direction >= 0)
                continue;
              Real alfa= - (ia.vec()*x - ia.rhs ())/
                (ia.vec()*direction);
-               
-             if (minalf > alfa) 
+
+             if (minalf > alfa)
                {
                  minidx = ia;
                  minalf = alfa;
@@ -165,13 +165,13 @@ Ineq_constrained_qp::constraint_solve (Vector start) const
          Real optimal_step = min (minalf, unbounded_alfa);
 
          Vector deltax=direction * optimal_step;
-         x += deltax;      
+         x += deltax;
          gradient += optimal_step * (quad * deltax);
-           
+
          DOUT << "step = " << optimal_step<< " (|dx| = " <<
-           deltax.norm() << ")\n";         
-          
-         if (minalf < unbounded_alfa) 
+           deltax.norm() << ")\n";
+
+         if (minalf < unbounded_alfa)
            {
              /* bumped into an edge. try again, in smaller space. */
              act.add (minidx.idx());
@@ -179,46 +179,46 @@ Ineq_constrained_qp::constraint_solve (Vector start) const
              continue;
            }
          /*ASSERT: we are at optimal solution for this "plane"*/
-  
-  
+
+
        }
-       
-      Vector lagrange_mult=act.get_lagrange (gradient);        
+
+      Vector lagrange_mult=act.get_lagrange (gradient);
       int m= min_elt_index (lagrange_mult);
-       
-      if (m>=0 && lagrange_mult (m) > 0) 
+
+      if (m>=0 && lagrange_mult (m) > 0)
        {
          break;                // optimal sol.
        }
-      else if (m<0) 
+      else if (m<0)
        {
          assert (gradient.norm() < EPS) ;
-           
+
          break;
        }
-       
+
       DOUT << "dropping cons " << m<<'\n';
       act.drop (m);
     }
   if (iterations >= MAXITER)
-    WARN<<"didn't converge!\n";
-  
+    WARN<<_("didn't converge!\n");
+
   DOUT <<  ": found " << x<<" in " << iterations <<" iterations\n";
   assert_solution (x);
   return x;
-} 
+}
+
 
-  
 Vector
 Ineq_constrained_qp::solve (Vector start) const
-{ 
+{
   /* no hassle if no constraints*/
-  if (! cons.size()) 
+  if (! cons.size())
     {
       Choleski_decomposition chol (quad);
       return - chol.solve (lin);
     }
-  else 
+  else
     {
       return constraint_solve (start);
     }