lilypond-0.0.78

diff --git a/lily/qlpsolve.cc b/lily/qlpsolve.cc
index 5c4d91b6e67f61ceb7839e90555ae65ddd4f86a1..31ec4a85b4dd10bc3f8ee76c7200b72f5ec6257d 100644 (file)
--- a/lily/qlpsolve.cc
+++ b/lily/qlpsolve.cc
@@ -9,6 +9,7 @@
   try fixed point arithmetic, to speed up lily.
   */
 
+#include "ineq-constrained-qp.hh"
 #include "qlpsolve.hh"
 #include "const.hh"
 #include "debug.hh"
@@ -142,6 +143,7 @@ Active_constraints::Active_constraints(Ineq_constrained_qp const *op)
        inactive.push(i);
     Choleski_decomposition chol(op->quad);
     H=chol.inverse();
+    OK();
 }
 
 /** Find the optimum which is in the planes generated by the active
@@ -153,127 +155,3 @@ Active_constraints::find_active_optimum(Vector g)
     return H*g;
 }
 
-/* *************************************************************** */
-
-int
-min_elt_index(Vector v)
-{
-    Real m=INFTY; int idx=-1;
-    for (int i = 0; i < v.dim(); i++){
-       if (v(i) < m) {
-           idx = i;
-           m = v(i);
-       }
-       assert(v(i) <= INFTY);
-    }
-    return idx;
-}
-
-
-/**the numerical solving. Mordecai Avriel, Nonlinear Programming: analysis and methods (1976)
-    Prentice Hall.
-
-    Section 13.3
-
-    This is a "projected gradient" algorithm. Starting from a point x
-    the next point is found in a direction determined by projecting
-    the gradient onto the active constraints.  (well, not really the
-    gradient. The optimal solution obeying the active constraints is
-    tried. This is why H = Q^-1 in initialisation) )
-
-
-    */
-Vector
-Ineq_constrained_qp::solve(Vector start) const 
-{    
-    if (!dim())
-       return Vector(0);
-    
-    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.
-    
-    Vector last_gradient(gradient);
-    int iterations=0;
-    
-    while (iterations++ < MAXITER) {
-       Vector direction= - act.find_active_optimum(gradient);
-               
-       mtor << "gradient "<< gradient<< "\ndirection " << direction<<"\n";
-       
-       if (direction.norm() > EPS) {
-           mtor << act.status() << '\n';
-           
-           Real minalf = INFTY;
-
-           Inactive_iter minidx(act);
-
-
-           /*
-    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++) {
-
-               if (ia.vec() * direction >= 0)
-                   continue;
-               Real alfa= - (ia.vec()*x - ia.rhs())/
-                   (ia.vec()*direction);
-               
-               if (minalf > alfa) {
-                   minidx = ia;
-                   minalf = alfa;
-               }
-           }
-           Real unbounded_alfa = 1.0;
-           Real optimal_step = min(minalf, unbounded_alfa);
-
-           Vector deltax=direction * optimal_step;
-           x += deltax;            
-           gradient += optimal_step * (quad * deltax);
-           
-           mtor << "step = " << optimal_step<< " (|dx| = " <<
-               deltax.norm() << ")\n";     
-          
-           if (minalf < unbounded_alfa) {
-               /* bumped into an edge. try again, in smaller space. */
-               act.add(minidx.idx());
-               mtor << "adding cons "<< minidx.idx()<<'\n';
-               continue;
-           }
-           /*ASSERT: we are at optimal solution for this "plane"*/
-    
-    
-       }
-       
-       Vector lagrange_mult=act.get_lagrange(gradient);        
-       int m= min_elt_index(lagrange_mult);
-       
-       if (m>=0 && lagrange_mult(m) > 0) {
-           break;              // optimal sol.
-       } else if (m<0) {
-           assert(gradient.norm() < EPS) ;
-           
-           break;
-       }
-       
-       mtor << "dropping cons " << m<<'\n';
-       act.drop(m);
-    }
-    if (iterations >= MAXITER)
-       WARN<<"didn't converge!\n";
-    
-    mtor <<  ": found " << x<<" in " << iterations <<" iterations\n";
-    assert_solution(x);
-    return x;
-} 
-
-