36#ifndef VIGRA_QUADPROG_HXX
37#define VIGRA_QUADPROG_HXX
40#include "mathutil.hxx"
42#include "linear_solve.hxx"
43#include "numerictraits.hxx"
44#include "array_vector.hxx"
50template <
class T,
class C1,
class C2,
class C3>
51bool quadprogAddConstraint(MultiArrayView<2, T, C1> & R, MultiArrayView<2, T, C2> & J, MultiArrayView<2, T, C3> & d,
52 int activeConstraintCount,
double& R_norm)
56 linalg::detail::qrGivensStepImpl(0,
subVector(d, activeConstraintCount, n),
57 J.subarray(Shape(activeConstraintCount,0), Shape(n,n)));
58 if (
abs(d(activeConstraintCount,0)) <= NumericTraits<T>::epsilon() * R_norm)
60 R_norm = std::max<T>(R_norm,
abs(d(activeConstraintCount,0)));
62 ++activeConstraintCount;
64 columnVector(R, Shape(0, activeConstraintCount - 1), activeConstraintCount) =
subVector(d, 0, activeConstraintCount);
68template <
class T,
class C1,
class C2,
class C3>
69void quadprogDeleteConstraint(MultiArrayView<2, T, C1> & R, MultiArrayView<2, T, C2> & J, MultiArrayView<2, T, C3> & u,
70 int activeConstraintCount,
int constraintToBeRemoved)
74 int newActiveConstraintCount = activeConstraintCount - 1;
76 if(constraintToBeRemoved == newActiveConstraintCount)
79 std::swap(u(constraintToBeRemoved,0), u(newActiveConstraintCount,0));
81 linalg::detail::qrGivensStepImpl(0, R.subarray(Shape(constraintToBeRemoved, constraintToBeRemoved),
82 Shape(newActiveConstraintCount,newActiveConstraintCount)),
83 J.subarray(Shape(constraintToBeRemoved, 0),
84 Shape(newActiveConstraintCount,newActiveConstraintCount)));
204template <
class T,
class C1,
class C2,
class C3,
class C4,
class C5,
class C6,
class C7>
217 constraintCount = me + mi;
220 "quadraticProgramming(): Matrix shape mismatch between G and g.");
221 vigra_precondition(
rowCount(x) == n,
222 "quadraticProgramming(): Output vector x has illegal shape.");
225 "quadraticProgramming(): Matrix CE has illegal shape.");
228 "quadraticProgramming(): Matrix CI has illegal shape.");
240 T f_value = 0.5 *
dot(g, x);
242 T epsilonZ = NumericTraits<T>::epsilon() *
sq(J.
norm(0)),
243 inf = std::numeric_limits<T>::infinity();
245 Matrix<T> R(n, n), r(constraintCount, 1), u(constraintCount,1);
246 T R_norm = NumericTraits<T>::one();
249 for (
int i=0; i < me; ++i)
260 : (-
dot(np, x) + ce(i,0)) /
dot(z, np);
266 f_value += 0.5 *
sq(step) *
dot(z, np);
268 vigra_precondition(vigra::detail::quadprogAddConstraint(R, J, d, i, R_norm),
269 "quadraticProgramming(): Equality constraints are linearly dependent.");
271 int activeConstraintCount = me;
275 for (
int i = 0; i < mi; ++i)
278 int constraintToBeAdded = 0;
280 for (
int i = activeConstraintCount-me; i < mi; ++i)
282 T s =
dot(
rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
286 constraintToBeAdded = i;
290 int iter = 0, maxIter = 10*mi;
291 while(iter++ < maxIter)
299 Matrix<T> z = transpose(J).subarray(Shape(0, activeConstraintCount), Shape(n,n))*
subVector(d, activeConstraintCount, n);
314 int constraintToBeRemoved = 0;
315 for (
int k = me; k < activeConstraintCount; ++k)
319 if (u(k,0) / r(k,0) < dualStep)
321 dualStep = u(k,0) / r(k,0);
322 constraintToBeRemoved = k;
328 T step = std::min(dualStep, primalStep);
337 if (primalStep == inf)
340 subVector(u, 0, activeConstraintCount) -= step *
subVector(r, 0, activeConstraintCount);
341 vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
342 --activeConstraintCount;
343 std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
350 f_value += 0.5 *
sq(step) *
dot(z, np);
352 subVector(u, 0, activeConstraintCount) -= step *
subVector(r, 0, activeConstraintCount);
353 u(activeConstraintCount,0) = step;
355 if (step == primalStep)
358 vigra::detail::quadprogAddConstraint(R, J, d, activeConstraintCount, R_norm);
359 std::swap(activeSet[constraintToBeAdded], activeSet[activeConstraintCount-me]);
360 ++activeConstraintCount;
365 vigra::detail::quadprogDeleteConstraint(R, J, u, activeConstraintCount, constraintToBeRemoved);
366 --activeConstraintCount;
367 std::swap(activeSet[constraintToBeRemoved-me], activeSet[activeConstraintCount-me]);
372 for (
int i = activeConstraintCount-me; i < mi; ++i)
375 T s =
dot(
rowVector(CI, activeSet[i]), x) - ci(activeSet[i], 0);
379 constraintToBeAdded = i;
Definition array_vector.hxx:514
TinyVector< MultiArrayIndex, N > type
Definition multi_shape.hxx:272
Base class for, and view to, MultiArray.
Definition multi_array.hxx:705
const difference_type & shape() const
Definition multi_array.hxx:1650
Definition matrix.hxx:125
NormTraits< Matrix >::NormType norm() const
Linear algebra functions.
Definition eigensystem.hxx:50
void identityMatrix(MultiArrayView< 2, T, C > &r)
Definition matrix.hxx:845
unsigned int quadraticProgramming(...)
void choleskySolve(MultiArrayView< 2, T, C1 > const &L, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x)
Definition linear_solve.hxx:1168
bool linearSolveUpperTriangular(const MultiArrayView< 2, T, C1 > &r, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition linear_solve.hxx:1068
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:684
NumericTraits< T >::Promote sq(T t)
The square function.
Definition mathutil.hxx:382
bool linearSolveLowerTriangular(const MultiArrayView< 2, T, C1 > &l, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition linear_solve.hxx:1118
MultiArrayView< 2, T, C > rowVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition matrix.hxx:697
bool choleskyDecomposition(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &L)
Definition linear_solve.hxx:959
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:671
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition fftw3.hxx:1002
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition matrix.hxx:727
MultiArrayView< 2, T, C > subVector(MultiArrayView< 2, T, C > const &m, int first, int end)
Definition matrix.hxx:761
PromoteTraits< V1, V2 >::Promote dot(RGBValue< V1, RIDX1, GIDX1, BIDX1 > const &r1, RGBValue< V2, RIDX2, GIDX2, BIDX2 > const &r2)
dot product
Definition rgbvalue.hxx:906