lily/qlpsolve.cc

   1 /*
   2   qlpsolve.cc -- implement Active_constraints, Inactive_iter
   3
   4   source file of the GNU LilyPond music typesetter
   5
   6   (c) 1996, 1997 Han-Wen Nienhuys <hanwen@stack.nl>
   7
   8   TODO:
   9   try fixed point arithmetic, to speed up lily.
  10   */
  11
  12 #include "ineq-constrained-qp.hh"
  13 #include "qlpsolve.hh"
  14 #include "const.hh"
  15 #include "debug.hh"
  16 #include "choleski.hh"
  17
  18 const Real TOL=1e-2;            // roughly 1/10 mm
  19
  20 String
  21 Active_constraints::status() const
  22 {
  23     String s("Active|Inactive [");
  24     for (int i=0; i< active.size(); i++) {
  25         s += String(active[i]) + " ";
  26     }
  27
  28     s+="| ";
  29     for (int i=0; i< inactive.size(); i++) {
  30         s += String(inactive[i]) + " ";
  31     }
  32     s+="]";
  33
  34     return s;
  35 }
  36
  37 void
  38 Active_constraints::OK()
  39 {
  40 #ifndef NDEBUG
  41     H.OK();
  42     A.OK();
  43     assert(active.size() +inactive.size() == opt->cons.size());
  44     assert(H.dim() == opt->dim());
  45     assert(active.size() == A.rows());
  46     Array<int> allcons;
  47
  48     for (int i=0; i < opt->cons.size(); i++)
  49         allcons.push(0);
  50     for (int i=0; i < active.size(); i++) {
  51         int j = active[i];
  52         allcons[j]++;
  53     }
  54     for (int i=0; i < inactive.size(); i++) {
  55         int j = inactive[i];
  56         allcons[j]++;
  57     }
  58     for (int i=0; i < allcons.size(); i++)
  59         assert(allcons[i] == 1);
  60 #endif
  61 }
  62
  63 Vector
  64 Active_constraints::get_lagrange(Vector gradient)
  65 {
  66     Vector l(A*gradient);
  67
  68     return l;
  69 }
  70
  71 void
  72 Active_constraints::add(int k)
  73 {
  74     // add indices
  75     int cidx=inactive[k];
  76     active.push(cidx);
  77
  78     inactive.swap(k,inactive.size()-1);
  79     inactive.pop();
  80
  81     Vector a( opt->cons[cidx] );
  82     // update of matrices
  83     Vector Ha = H*a;
  84     Real aHa = a*Ha;
  85     Vector addrow(Ha.dim());
  86     if (abs(aHa) > EPS) {
  87         /*
  88           a != 0, so if Ha = O(EPS), then
  89           Ha * aH / aHa = O(EPS^2/EPS)
  90
  91           if H*a == 0, the constraints are dependent.
  92           */
  93         H -= Matrix(Ha/aHa , Ha);
  94
  95
  96         /*
  97           sorry, don't know how to justify this. ..
  98           */
  99         addrow=Ha;
 100         addrow/= aHa;
 101         A -= Matrix(A*a, addrow);
 102         A.insert_row(addrow,A.rows());
 103     }else
 104         WARN << "degenerate constraints";
 105 }
 106
 107 void
 108 Active_constraints::drop(int k)
 109 {
 110     int q=active.size()-1;
 111
 112         // drop indices
 113     inactive.push(active[k]);
 114     active.swap(k,q);
 115     A.swap_rows(k,q);
 116     active.pop();
 117
 118     Vector a(A.row(q));
 119     if (a.norm() > EPS) {
 120         /*
 121
 122          */
 123         Real q = a*opt->quad*a;
 124         H += Matrix(a,a/q);
 125         A -= A*opt->quad*Matrix(a,a/q);
 126     }else
 127         WARN << "degenerate constraints";
 128 #ifndef NDEBUG
 129     Vector rem_row(A.row(q));
 130     assert(rem_row.norm() < EPS);
 131 #endif
 132
 133     A.delete_row(q);
 134 }
 135
 136
 137 Active_constraints::Active_constraints(Ineq_constrained_qp const *op)
 138     :       A(0,op->dim()),
 139             H(op->dim()),
 140             opt(op)
 141 {
 142     for (int i=0; i < op->cons.size(); i++)
 143         inactive.push(i);
 144     Choleski_decomposition chol(op->quad);
 145     H=chol.inverse();
 146     OK();
 147 }
 148
 149 /** Find the optimum which is in the planes generated by the active
 150     constraints.
 151     */
 152 Vector
 153 Active_constraints::find_active_optimum(Vector g)
 154 {
 155     return H*g;
 156 }
 157