2 ineq-constrained-qp.cc -- implement Ineq_constrained_qp
4 source file of the GNU LilyPond music typesetter
6 (c) 1996,1997 Han-Wen Nienhuys <hanwen@stack.nl>
8 #include "ineq-constrained-qp.hh"
13 MAy be this also should go into a library
16 const int MAXITER=100; // qlpsolve.hh
19 assume x(idx) == value, and adjust constraints, lin and quad accordingly
24 Ineq_constrained_qp::eliminate_var(int idx, Real value)
26 Vector row(quad.row(idx));
31 quad.delete_column(idx);
37 for (int i=0; i < cons.size(); i++) {
38 consrhs[i] -= cons[i](idx) *value;
44 Ineq_constrained_qp::add_inequality_cons(Vector c, double r)
50 Ineq_constrained_qp::Ineq_constrained_qp(int novars):
58 Ineq_constrained_qp::OK() const
60 #if !defined(NDEBUG) && defined(PARANOID)
61 assert(cons.size() == consrhs.size());
62 Matrix Qdif= quad - quad.transposed();
63 assert(Qdif.norm()/quad.norm() < EPS);
69 Ineq_constrained_qp::eval (Vector v)
71 return v * quad * v + lin * v + const_term;
76 min_elt_index(Vector v)
80 for (int i = 0; i < v.dim(); i++){
85 assert(v(i) <= INFTY_f);
91 /**the numerical solving. Mordecai Avriel, Nonlinear Programming: analysis and methods (1976)
96 This is a "projected gradient" algorithm. Starting from a point x
97 the next point is found in a direction determined by projecting
98 the gradient onto the active constraints. (well, not really the
99 gradient. The optimal solution obeying the active constraints is
100 tried. This is why H = Q^-1 in initialisation) )
105 Ineq_constrained_qp::solve(Vector start) const
110 // quad.try_set_band();
112 Active_constraints act(this);
119 Vector gradient=quad*x+lin;
120 // Real fvalue = x*quad*x/2 + lin*x + const_term;
123 Vector last_gradient(gradient);
126 while (iterations++ < MAXITER) {
127 Vector direction= - act.find_active_optimum(gradient);
129 mtor << "gradient "<< gradient<< "\ndirection " << direction<<"\n";
131 if (direction.norm() > EPS) {
132 mtor << act.status() << '\n';
134 Real minalf = INFTY_f;
136 Inactive_iter minidx(act);
140 we know the optimum on this "hyperplane". Check if we
141 bump into the edges of the simplex
144 for (Inactive_iter ia(act); ia.ok(); ia++) {
146 if (ia.vec() * direction >= 0)
148 Real alfa= - (ia.vec()*x - ia.rhs())/
149 (ia.vec()*direction);
156 Real unbounded_alfa = 1.0;
157 Real optimal_step = min(minalf, unbounded_alfa);
159 Vector deltax=direction * optimal_step;
161 gradient += optimal_step * (quad * deltax);
163 mtor << "step = " << optimal_step<< " (|dx| = " <<
164 deltax.norm() << ")\n";
166 if (minalf < unbounded_alfa) {
167 /* bumped into an edge. try again, in smaller space. */
168 act.add(minidx.idx());
169 mtor << "adding cons "<< minidx.idx()<<'\n';
172 /*ASSERT: we are at optimal solution for this "plane"*/
177 Vector lagrange_mult=act.get_lagrange(gradient);
178 int m= min_elt_index(lagrange_mult);
180 if (m>=0 && lagrange_mult(m) > 0) {
181 break; // optimal sol.
183 assert(gradient.norm() < EPS) ;
188 mtor << "dropping cons " << m<<'\n';
191 if (iterations >= MAXITER)
192 WARN<<"didn't converge!\n";
194 mtor << ": found " << x<<" in " << iterations <<" iterations\n";