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
);