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