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
116 Active_constraints
act (this);
121 Vector gradient
=quad_
*x
+lin_
;
122 // Real fvalue = x*quad_*x/2 + lin*x + const_term;
125 Vector
last_gradient (gradient
);
128 while (iterations
++ < MAXITER
)
131 if (experimental_features_global_b
)
135 Vector direction
= - act
.find_active_optimum (gradient
);
137 DOUT
<< "gradient "<< gradient
<< "\ndirection " << direction
<<"\n";
139 if (direction
.norm() > EPS
)
141 DOUT
<< act
.status() << '\n';
143 Real unbounded_alfa
= 1.0;
144 Real minalf
= unbounded_alfa
;
146 Inactive_iter
minidx (act
);
150 we know the optimum on this "hyperplane". Check if we
151 bump into the edges of the simplex
154 for (Inactive_iter
ia (act
); ia
.ok(); ia
++)
156 Real dot
= ia
.vec() * direction
;
161 Real numerator
= ia
.rhs () - ia
.vec()*x
;
165 warning (String ("Ineq_constrained_qp::solve (): Constraint off by ") + numerator
);
171 Real alfa
= numerator
/ dot
;
181 Real optimal_step
= minalf
;
183 Vector deltax
=direction
* optimal_step
;
185 gradient
+= optimal_step
* (quad_
* deltax
);
187 DOUT
<< "step = " << optimal_step
<< " (|dx| = " <<
188 deltax
.norm() << ")\n";
190 if (minalf
< unbounded_alfa
)
192 /* bumped into an edge. try again, in smaller space. */
193 act
.add (minidx
.idx());
194 DOUT
<< "adding cons "<< minidx
.idx()<<'\n';
197 /*ASSERT: we are at optimal solution for this "plane"*/
202 Vector lagrange_mult
=act
.get_lagrange (gradient
);
203 int m
= min_elt_index (lagrange_mult
);
205 if (m
>=0 && lagrange_mult (m
) > 0)
207 break; // optimal sol.
211 assert (gradient
.norm() < EPS
) ;
216 DOUT
<< "dropping cons " << m
<<'\n';
219 if (iterations
>= MAXITER
)
220 WARN
<<_("didn't converge!\n");
222 DOUT
<< ": found " << x
<<" in " << iterations
<<" iterations\n";
229 Ineq_constrained_qp::solve (Vector start
) const
231 /* no hassle if no constraints*/
234 Choleski_decomposition
chol (quad_
);
235 return - chol
.solve (lin_
);
239 return constraint_solve (start
);
244 Ineq_constrained_qp::dim () const