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"
11 #include "choleski.hh"
15 MAy be this also should go into a library
18 const int MAXITER=100; // qlpsolve.hh
22 assume x (idx) == value, and adjust constraints, lin and quad accordingly
27 Ineq_constrained_qp::eliminate_var (int idx, Real value)
29 Vector row (quad_.row (idx));
32 quad_.delete_row (idx);
34 quad_.delete_column (idx);
40 for (int i=0; i < cons_.size(); i++)
42 consrhs_[i] -= cons_[i](idx) *value;
48 Ineq_constrained_qp::add_inequality_cons (Vector c, double r)
54 Ineq_constrained_qp::Ineq_constrained_qp (int novars):
62 Ineq_constrained_qp::OK() const
64 #if !defined (NDEBUG) && defined (PARANOID)
65 assert (cons_.size() == consrhs_.size ());
66 Matrix Qdif= quad_ - quad_.transposed();
67 assert (Qdif.norm()/quad_.norm () < EPS);
73 Ineq_constrained_qp::eval (Vector v)
75 return v * quad_ * v + lin_ * v + const_term_;
80 min_elt_index (Vector v)
84 for (int i = 0; i < v.dim(); i++)
91 assert (v (i) <= infinity_f);
97 /**the numerical solving. Mordecai Avriel, Nonlinear Programming: analysis and methods (1976)
102 This is a "projected gradient" algorithm. Starting from a point x
103 the next point is found in a direction determined by projecting
104 the gradient onto the active constraints. (well, not really the
105 gradient. The optimal solution obeying the active constraints is
106 tried. This is why H = Q^-1 in initialisation))
111 Ineq_constrained_qp::constraint_solve (Vector start) const
117 if (quad_.dim() > 10)
118 quad_.try_set_band();
120 Active_constraints act (this);
125 Vector gradient=quad_*x+lin_;
126 // Real fvalue = x*quad_*x/2 + lin*x + const_term;
129 Vector last_gradient (gradient);
132 while (iterations++ < MAXITER)
135 if (experimental_features_global_b)
139 Vector direction= - act.find_active_optimum (gradient);
141 DOUT << "gradient "<< gradient<< "\ndirection " << direction<<"\n";
143 if (direction.norm() > EPS)
145 DOUT << act.status() << '\n';
147 Real unbounded_alfa = 1.0;
148 Real minalf = unbounded_alfa;
150 Inactive_iter minidx (act);
154 we know the optimum on this "hyperplane". Check if we
155 bump into the edges of the simplex
158 for (Inactive_iter ia (act); ia.ok(); ia++)
160 Real dot = ia.vec() * direction;
165 Real numerator = ia.rhs () - ia.vec()*x;
169 warning (String ("Ineq_constrained_qp::solve (): Constraint off by ") + numerator);
175 Real alfa = numerator / dot;
185 Real optimal_step = minalf;
187 Vector deltax=direction * optimal_step;
189 gradient += optimal_step * (quad_ * deltax);
191 DOUT << "step = " << optimal_step<< " (|dx| = " <<
192 deltax.norm() << ")\n";
194 if (minalf < unbounded_alfa)
196 /* bumped into an edge. try again, in smaller space. */
197 act.add (minidx.idx());
198 DOUT << "adding cons "<< minidx.idx()<<'\n';
201 /*ASSERT: we are at 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';
223 if (iterations >= MAXITER)
224 WARN<<_("didn't converge!\n");
226 DOUT << ": found " << x<<" in " << iterations <<" iterations\n";
233 Ineq_constrained_qp::solve (Vector start) const
235 /* no hassle if no constraints*/
238 Choleski_decomposition chol (quad_);
239 return - chol.solve (lin_);
243 return constraint_solve (start);
248 Ineq_constrained_qp::dim () const