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"
14 MAy be this also should go into a library
17 const int MAXITER=100; // qlpsolve.hh
20 assume x (idx) == value, and adjust constraints, lin and quad accordingly
25 Ineq_constrained_qp::eliminate_var (int idx, Real value)
27 Vector row (quad.row (idx));
30 quad.delete_row (idx);
32 quad.delete_column (idx);
38 for (int i=0; i < cons.size(); i++)
40 consrhs[i] -= cons[i](idx) *value;
46 Ineq_constrained_qp::add_inequality_cons (Vector c, double r)
52 Ineq_constrained_qp::Ineq_constrained_qp (int novars):
60 Ineq_constrained_qp::OK() const
62 #if !defined (NDEBUG) && defined (PARANOID)
63 assert (cons.size() == consrhs.size ());
64 Matrix Qdif= quad - quad.transposed();
65 assert (Qdif.norm()/quad.norm () < EPS);
71 Ineq_constrained_qp::eval (Vector v)
73 return v * quad * v + lin * v + const_term;
78 min_elt_index (Vector v)
82 for (int i = 0; i < v.dim(); i++)
89 assert (v (i) <= infinity_f);
95 /**the numerical solving. Mordecai Avriel, Nonlinear Programming: analysis and methods (1976)
100 This is a "projected gradient" algorithm. Starting from a point x
101 the next point is found in a direction determined by projecting
102 the gradient onto the active constraints. (well, not really the
103 gradient. The optimal solution obeying the active constraints is
104 tried. This is why H = Q^-1 in initialisation))
109 Ineq_constrained_qp::constraint_solve (Vector start) const
118 Active_constraints act (this);
123 Vector gradient=quad*x+lin;
124 // Real fvalue = x*quad*x/2 + lin*x + const_term;
127 Vector last_gradient (gradient);
130 while (iterations++ < MAXITER)
132 Vector direction= - act.find_active_optimum (gradient);
134 DOUT << "gradient "<< gradient<< "\ndirection " << direction<<"\n";
136 if (direction.norm() > EPS)
138 DOUT << act.status() << '\n';
140 Real minalf = infinity_f;
142 Inactive_iter minidx (act);
146 we know the optimum on this "hyperplane". Check if we
147 bump into the edges of the simplex
150 for (Inactive_iter ia (act); ia.ok(); ia++)
153 if (ia.vec() * direction >= 0)
155 Real alfa= - (ia.vec()*x - ia.rhs ())/
156 (ia.vec()*direction);
164 Real unbounded_alfa = 1.0;
165 Real optimal_step = min (minalf, unbounded_alfa);
167 Vector deltax=direction * optimal_step;
169 gradient += optimal_step * (quad * deltax);
171 DOUT << "step = " << optimal_step<< " (|dx| = " <<
172 deltax.norm() << ")\n";
174 if (minalf < unbounded_alfa)
176 /* bumped into an edge. try again, in smaller space. */
177 act.add (minidx.idx());
178 DOUT << "adding cons "<< minidx.idx()<<'\n';
181 /*ASSERT: we are at optimal solution for this "plane"*/
186 Vector lagrange_mult=act.get_lagrange (gradient);
187 int m= min_elt_index (lagrange_mult);
189 if (m>=0 && lagrange_mult (m) > 0)
191 break; // optimal sol.
195 assert (gradient.norm() < EPS) ;
200 DOUT << "dropping cons " << m<<'\n';
203 if (iterations >= MAXITER)
204 WARN<<"didn't converge!\n";
206 DOUT << ": found " << x<<" in " << iterations <<" iterations\n";
213 Ineq_constrained_qp::solve (Vector start) const
215 /* no hassle if no constraints*/
218 Choleski_decomposition chol (quad);
219 return - chol.solve (lin);
223 return constraint_solve (start);