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