2 qlpsolve.cc -- implement Active_constraints, Inactive_iter
4 source file of the GNU LilyPond music typesetter
6 (c) 1996, 1997--1998 Han-Wen Nienhuys <hanwen@cs.uu.nl>
9 try fixed point arithmetic, to speed up lily.
12 #include "ineq-constrained-qp.hh"
13 #include "qlpsolve.hh"
15 #include "choleski.hh"
17 const Real TOL=1e-1; // roughly 1/30 mm
20 Active_constraints::status() const
22 String s ("Active|Inactive [");
23 for (int i=0; i< active.size(); i++)
25 s += to_str (active[i]) + " ";
29 for (int i=0; i< inactive.size(); i++)
31 s += to_str (inactive[i]) + " ";
39 Active_constraints::OK()
44 assert (active.size() +inactive.size () == opt->cons_.size ());
45 assert (H.dim() == opt->dim ());
46 assert (active.size() == A.rows ());
49 for (int i=0; i < opt->cons_.size(); i++)
51 for (int i=0; i < active.size(); i++)
56 for (int i=0; i < inactive.size(); i++)
61 for (int i=0; i < allcons.size(); i++)
62 assert (allcons[i] == 1);
67 Active_constraints::get_lagrange (Vector gradient)
73 Active_constraints::add_constraint (int k)
78 Vector a (opt->cons_[cidx]);
82 Vector addrow (Ha.dim());
83 bool degenerate = (abs (aHa) < EPS);
87 warning (String ("Active_constraints::add ():")
88 + _("degenerate constraints"));
89 DOUT << "Ha = " << Ha.str () << '\n';
91 a != 0, so if Ha = O(EPS), then
92 Ha * aH / aHa = O(EPS^2/EPS)
94 if H*a == 0, the constraints are dependent.
96 degenerate_count_i_ ++;
101 inactive.swap (k,inactive.size()-1);
104 H -= Matrix (Ha/aHa , Ha);
108 A -= Matrix (A*a, addrow);
109 A.insert_row (addrow,A.rows());
114 Active_constraints::drop_constraint (int k)
116 int q=active.size()-1;
119 Vector a (A.row (q));
123 inactive.push (active[k]);
130 Real aqa = a*opt->quad_*a;
131 Matrix aaq (a,a/aqa);
133 A -= A*opt->quad_*aaq;
137 degenerate_count_i_ ++;
138 warning (String ("Active_constraints::drop ():")
139 + _("degenerate constraints"));
144 Active_constraints::Active_constraints (Ineq_constrained_qp const *op)
149 for (int i=0; i < op->cons_.size(); i++)
151 Choleski_decomposition chol (op->quad_);
159 degenerate_count_i_ = 0;
162 /** Find the optimum which is in the planes generated by the active
166 Active_constraints::find_active_optimum (Vector g)