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--1999 Han-Wen Nienhuys <hanwen@cs.uu.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 += to_str (active[i]) + " ";
  26     }
  27
  28   s+="| ";
  29   for (int i=0; i< inactive.size(); i++)
  30     {
  31       s += to_str (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_constraint (int k)
  74 {
  75   // add indices
  76   int cidx=inactive[k];
  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   bool degenerate =  (abs (aHa) < EPS);
  84
  85   if (degenerate)
  86     {
  87       warning (String ("Active_constraints::add ():")
  88         + _("degenerate constraints"));
  89       DOUT << "Ha = "  << Ha.str () << '\n';
  90       /*
  91         a != 0, so if Ha = O(EPS), then
  92         Ha * aH / aHa = O(EPS^2/EPS)
  93
  94         if H*a == 0, the constraints are dependent.
  95       */
  96       degenerate_count_i_ ++;
  97     }
  98   if (!degenerate)
  99     {
 100       active.push (cidx);
 101       inactive.swap (k,inactive.size()-1);
 102       inactive.pop();
 103
 104       H -= Matrix (Ha/aHa , Ha);
 105
 106       addrow=Ha;
 107       addrow /= aHa;
 108       A -= Matrix (A*a, addrow);
 109       A.insert_row (addrow,A.rows());
 110     }
 111 }
 112
 113 void
 114 Active_constraints::drop_constraint (int k)
 115 {
 116   int q=active.size()-1;
 117
 118
 119   Vector a (A.row (q));
 120   if (a.norm() > EPS)
 121     {
 122       // drop indices
 123       inactive.push (active[k]);
 124       active.swap (k,q);
 125       A.swap_rows (k,q);
 126       active.pop();
 127       /*
 128
 129        */
 130       Real aqa = a*opt->quad_*a;
 131       Matrix aaq (a,a/aqa);
 132       H += aaq;
 133       A -= A*opt->quad_*aaq;
 134
 135       A.delete_row (q);
 136     }else {
 137       degenerate_count_i_ ++;
 138       warning (String ("Active_constraints::drop ():")
 139         + _("degenerate constraints"));
 140     }
 141 }
 142
 143
 144 Active_constraints::Active_constraints (Ineq_constrained_qp const *op)
 145   : A(0,op->dim()),
 146     H(op->dim()),
 147     opt (op)
 148 {
 149   for (int i=0; i < op->cons_.size(); i++)
 150     inactive.push (i);
 151   Choleski_decomposition chol (op->quad_);
 152
 153   /*
 154     ugh.
 155     */
 156   H=chol.inverse();
 157   OK();
 158
 159   degenerate_count_i_ = 0;
 160 }
 161
 162 /** Find the optimum which is in the planes generated by the active
 163   constraints.
 164   */
 165 Vector
 166 Active_constraints::find_active_optimum (Vector g)
 167 {
 168   return H*g;
 169 }