From 3b7daf43bf5af4aaeb58801ebd4b064fb719ba63 Mon Sep 17 00:00:00 2001
From: fred <fred>
Date: Sun, 24 Mar 2002 19:49:50 +0000
Subject: [PATCH] lilypond-0.0.78

---
 lily/qlpsolve.cc | 126 +----------------------------------------------
 1 file changed, 2 insertions(+), 124 deletions(-)

diff --git a/lily/qlpsolve.cc b/lily/qlpsolve.cc
index 5c4d91b6e6..31ec4a85b4 100644
--- 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;
-} 
-
-    
-- 
2.39.5