2 ineq-constrained-qp.cc -- implement Ineq_constrained_qp
4 source file of the GNU LilyPond music typesetter
6 (c) 1996, 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
8 #include "ineq-constrained-qp.hh"
11 #include "choleski.hh"
15 May be this also should go into a library
18 const int MAXITER=100; // qlpsolve.hh
23 assume x (idx) == value, and adjust constraints, lin and quad accordingly
28 Ineq_constrained_qp::eliminate_var (int idx, Real value)
30 Vector row (quad_.row (idx));
33 quad_.delete_row (idx);
35 quad_.delete_column (idx);
41 for (int i=0; i < cons_.size(); i++)
43 consrhs_[i] -= cons_[i](idx) *value;
49 Ineq_constrained_qp::add_inequality_cons (Vector c, double r)
55 Ineq_constrained_qp::Ineq_constrained_qp (int novars):
63 Ineq_constrained_qp::OK() const
65 #if !defined (NDEBUG) && defined (PARANOID)
66 assert (cons_.size() == consrhs_.size ());
67 Matrix Qdif= quad_ - quad_.transposed();
68 assert (Qdif.norm()/quad_.norm () < EPS);
74 Ineq_constrained_qp::eval (Vector v)
76 return v * quad_ * v + lin_ * v + const_term_;
81 min_elt_index (Vector v)
85 for (int i = 0; i < v.dim(); i++)
92 assert (v (i) <= infinity_f);
98 /**the numerical solving. Mordecai Avriel, Nonlinear Programming: analysis and methods (1976)
103 This is a "projected gradient" algorithm. Starting from a point x
104 the next point is found in a direction determined by projecting
105 the gradient onto the active constraints. (well, not really the
106 gradient. The optimal solution obeying the active constraints is
107 tried. This is why H = Q^-1 in initialisation))
112 Ineq_constrained_qp::constraint_solve (Vector start) const
117 Active_constraints act (this);
122 Vector gradient=quad_*x+lin_;
125 // Real fvalue = x*quad_*x/2 + lin*x + const_term;// it's no use.
126 Vector last_gradient (gradient);
129 while (iterations++ < MAXITER && act.degenerate_count_i_ < MAXDEGEN)
131 if (experimental_features_global_b)
134 Vector direction= - act.find_active_optimum (gradient);
136 DOUT << "gradient "<< gradient<< "\ndirection " << direction<< '\n';
138 if (direction.norm() > EPS)
140 DOUT << act.status() << '\n';
142 Real unbounded_alfa = 1.0;
143 Real minalf = unbounded_alfa;
145 Inactive_iter minidx (act);
149 we know the optimum on this "hyperplane". Check if we
150 bump into the edges of the simplex
153 for (Inactive_iter ia (act); ia.ok(); ia++)
155 Real dot = ia.vec() * direction;
156 Real mindot = (experimental_features_global_b)
164 Real numerator = ia.rhs () - ia.vec()*x;
169 warning (_f ("Ineq_constrained_qp::solve (): Constraint off by %f", numerator));
170 act.degenerate_count_i_ ++;
177 Real alfa = numerator / dot;
187 Real optimal_step = minalf;
189 Vector deltax = direction * optimal_step;
191 gradient += optimal_step * (quad_ * deltax);
193 DOUT << "step = " << optimal_step << " (|dx| = " <<
194 to_str (deltax.norm()) << ")\n";
196 if (minalf < unbounded_alfa)
198 /* bumped into an edge. try again, in smaller space. */
199 act.add_constraint (minidx.idx());
200 DOUT << "adding cons "<< minidx.idx () << '\n';
203 /*ASSERT: we are at the optimal solution for this "plane"*/
206 Vector lagrange_mult=act.get_lagrange (gradient);
207 int m= min_elt_index (lagrange_mult);
209 if (m>=0 && lagrange_mult (m) > 0)
211 break; // optimal sol.
215 assert (gradient.norm() < EPS) ;
220 DOUT << "dropping cons " << m << '\n';
221 act.drop_constraint (m);
223 if (iterations >= MAXITER)
224 WARN << _ ("didn't converge!") << '\n';
225 if (act.degenerate_count_i_ >= MAXDEGEN)
226 WARN << _ ("Too much degeneracy. ") << '\n';
227 DOUT << ": found " << x << " in " << iterations <<" iterations\n";
234 Ineq_constrained_qp::solve (Vector start) const
236 /* no hassle if no constraints*/
239 Choleski_decomposition chol (quad_);
240 return - chol.solve (lin_);
244 return constraint_solve (start);
249 Ineq_constrained_qp::dim () const
258 Ineq_constrained_qp::assert_solution (Vector sol) const
262 for (int i=0; sol_b && i < cons_.size(); i++)
264 Real R=cons_[i] * sol- consrhs_[i];
271 Ineq_constrained_qp::print() const
274 DOUT << "Quad " << quad_;
275 DOUT << "lin " << lin_ << '\n'
276 << "const " << const_term_<< '\n';
277 for (int i=0; i < cons_.size(); i++)
279 DOUT << "constraint["<<i<<"]: " << cons_[i] << " >= " << consrhs_[i];