partial: 0.0.39-1.jcn

[lilypond.git] / src / qlpsolve.cc

diff --git a/src/qlpsolve.cc b/src/qlpsolve.cc
index 88b703bbbbad7965b56c9dfb6d7eee5e90d703f7..5738be2f3b0223c9236af8e481e72685f2f0bc17 100644 (file)
--- a/src/qlpsolve.cc
+++ b/src/qlpsolve.cc
@@ -9,12 +9,12 @@ String
 Active_constraints::status() const
 {
     String s("Active|Inactive [");
-    for (int i=0; i< active.sz(); i++) {
+    for (int i=0; i< active.size(); i++) {
        s += String(active[i]) + " ";
     }
 
     s+="| ";
-    for (int i=0; i< inactive.sz(); i++) {
+    for (int i=0; i< inactive.size(); i++) {
        s += String(inactive[i]) + " ";
     }
     s+="]";
@@ -23,26 +23,29 @@ Active_constraints::status() const
 }
 
 void
-Active_constraints::OK() {
+Active_constraints::OK()
+{
+    #ifndef NDEBUG
     H.OK();
     A.OK();
-    assert(active.sz() +inactive.sz() == opt->cons.sz());
+    assert(active.size() +inactive.size() == opt->cons.size());
     assert(H.dim() == opt->dim());
-    assert(active.sz() == A.rows());
-    svec<int> allcons;
+    assert(active.size() == A.rows());
+    Array<int> allcons;
 
-    for (int i=0; i < opt->cons.sz(); i++)
-       allcons.add(0);
-    for (int i=0; i < active.sz(); 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]++;
     }
-    for (int i=0; i < inactive.sz(); i++) {
+    for (int i=0; i < inactive.size(); i++) {
        int j = inactive[i];
        allcons[j]++;
     }
-    for (int i=0; i < allcons.sz(); i++)
+    for (int i=0; i < allcons.size(); i++)
        assert(allcons[i] == 1);
+#endif
 }
 
 Vector
@@ -58,9 +61,9 @@ Active_constraints::add(int k)
 {
     // add indices
     int cidx=inactive[k];
-    active.add(cidx);
+    active.push(cidx);
 
-    inactive.swap(k,inactive.sz()-1);
+    inactive.swap(k,inactive.size()-1);
     inactive.pop();
 
     Vector a( opt->cons[cidx] );
@@ -92,10 +95,10 @@ Active_constraints::add(int k)
 void
 Active_constraints::drop(int k)
 {
-    int q=active.sz()-1;
+    int q=active.size()-1;
 
         // drop indices
-    inactive.add(active[k]);
+    inactive.push(active[k]);
     active.swap(k,q);
     A.swap_rows(k,q);
     active.pop();
@@ -124,13 +127,13 @@ Active_constraints::Active_constraints(Ineq_constrained_qp const *op)
            H(op->dim()),
            opt(op)
 {
-    for (int i=0; i < op->cons.sz(); i++)
-       inactive.add(i);
+    for (int i=0; i < op->cons.size(); i++)
+       inactive.push(i);
     Choleski_decomposition chol(op->quad);
     H=chol.inverse();
 }
 
-/* Find the optimum which is in the planes generated by the active
+/** Find the optimum which is in the planes generated by the active
     constraints.        
     */
 Vector
@@ -139,7 +142,7 @@ Active_constraints::find_active_optimum(Vector g)
     return H*g;
 }
 
-/****************************************************************/
+/* *************************************************************** */
 
 int
 min_elt_index(Vector v)
@@ -155,7 +158,20 @@ min_elt_index(Vector v)
     return idx;
 }
 
-///the numerical solving
+
+/**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 
 {    
@@ -167,8 +183,9 @@ Ineq_constrained_qp::solve(Vector start) const
     
     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;
     
@@ -245,17 +262,4 @@ Ineq_constrained_qp::solve(Vector start) const
     return x;
 } 
 
-/** 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) )
-
-
-    */