row.del (idx);
lin +=row ;
- for (int i=0; i < cons.size(); i++)
+ for (int i=0; i < cons.size(); i++)
{
consrhs[i] -= cons[i](idx) *value;
cons[i].del (idx);
assert (cons.size() == consrhs.size ());
Matrix Qdif= quad - quad.transposed();
assert (Qdif.norm()/quad.norm () < EPS);
-#endif
+#endif
}
-
+
Real
Ineq_constrained_qp::eval (Vector v)
int
min_elt_index (Vector v)
{
- Real m=infinity_f;
+ Real m=infinity_f;
int idx=-1;
for (int i = 0; i < v.dim(); i++)
{
- if (v (i) < m)
+ if (v (i) < m)
{
idx = i;
m = v (i);
*/
Vector
-Ineq_constrained_qp::constraint_solve (Vector start) const
-{
+Ineq_constrained_qp::constraint_solve (Vector start) const
+{
if (!dim())
return Vector (0);
-
+
// experimental
if (quad.dim() > 10)
quad.try_set_band();
-
+
Active_constraints act (this);
- act.OK();
+ act.OK();
+
-
Vector x (start);
Vector gradient=quad*x+lin;
// Real fvalue = x*quad*x/2 + lin*x + const_term;
Vector last_gradient (gradient);
int iterations=0;
-
- while (iterations++ < MAXITER)
+
+ while (iterations++ < MAXITER)
{
Vector direction= - act.find_active_optimum (gradient);
-
+
DOUT << "gradient "<< gradient<< "\ndirection " << direction<<"\n";
-
- if (direction.norm() > EPS)
+
+ if (direction.norm() > EPS)
{
DOUT << act.status() << '\n';
-
+
Real minalf = infinity_f;
Inactive_iter minidx (act);
we know the optimum on this "hyperplane". Check if we
bump into the edges of the simplex
*/
-
- for (Inactive_iter ia (act); ia.ok(); ia++)
+
+ for (Inactive_iter ia (act); ia.ok(); ia++)
{
if (ia.vec() * direction >= 0)
continue;
Real alfa= - (ia.vec()*x - ia.rhs ())/
(ia.vec()*direction);
-
- if (minalf > alfa)
+
+ if (minalf > alfa)
{
minidx = ia;
minalf = alfa;
Real optimal_step = min (minalf, unbounded_alfa);
Vector deltax=direction * optimal_step;
- x += deltax;
+ x += deltax;
gradient += optimal_step * (quad * deltax);
-
+
DOUT << "step = " << optimal_step<< " (|dx| = " <<
- deltax.norm() << ")\n";
-
- if (minalf < unbounded_alfa)
+ deltax.norm() << ")\n";
+
+ if (minalf < unbounded_alfa)
{
/* bumped into an edge. try again, in smaller space. */
act.add (minidx.idx());
continue;
}
/*ASSERT: we are at optimal solution for this "plane"*/
-
-
+
+
}
-
- Vector lagrange_mult=act.get_lagrange (gradient);
+
+ Vector lagrange_mult=act.get_lagrange (gradient);
int m= min_elt_index (lagrange_mult);
-
- if (m>=0 && lagrange_mult (m) > 0)
+
+ if (m>=0 && lagrange_mult (m) > 0)
{
break; // optimal sol.
}
- else if (m<0)
+ else if (m<0)
{
assert (gradient.norm() < EPS) ;
-
+
break;
}
-
+
DOUT << "dropping cons " << m<<'\n';
act.drop (m);
}
if (iterations >= MAXITER)
- WARN<<"didn't converge!\n";
-
+ WARN<<_("didn't converge!\n");
+
DOUT << ": found " << x<<" in " << iterations <<" iterations\n";
assert_solution (x);
return x;
-}
+}
+
-
Vector
Ineq_constrained_qp::solve (Vector start) const
-{
+{
/* no hassle if no constraints*/
- if (! cons.size())
+ if (! cons.size())
{
Choleski_decomposition chol (quad);
return - chol.solve (lin);
}
- else
+ else
{
return constraint_solve (start);
}