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)
54 typedef typename MultiArrayShape<2>::type Shape;
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)
72 typedef typename MultiArrayShape<2>::type Shape;
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>
211 using namespace linalg;
219 vigra_precondition(columnCount(
G) == n && rowCount(
G) == n,
220 "quadraticProgramming(): Matrix shape mismatch between G and g.");
221 vigra_precondition(rowCount(x) == n,
222 "quadraticProgramming(): Output vector x has illegal shape.");
223 vigra_precondition((
me > 0 && columnCount(
CE) == n && rowCount(
CE) ==
me) ||
224 (
me == 0 && columnCount(
CE) == 0),
225 "quadraticProgramming(): Matrix CE has illegal shape.");
226 vigra_precondition((
mi > 0 && columnCount(
CI) == n && rowCount(
CI) ==
mi) ||
227 (
mi == 0 && columnCount(
CI) == 0),
228 "quadraticProgramming(): Matrix CI has illegal shape.");
233 choleskyDecomposition(
G,
L);
235 choleskySolve(
L, -g, x);
237 linearSolveLowerTriangular(
L,
J,
J);
242 T
epsilonZ = NumericTraits<T>::epsilon() * sq(
J.norm(0)),
243 inf = std::numeric_limits<T>::infinity();
246 T
R_norm = NumericTraits<T>::one();
249 for (
int i=0;
i <
me; ++
i)
254 linearSolveUpperTriangular(R.subarray(Shape(0, 0), Shape(
i,
i)),
260 : (-dot(
np, x) +
ce(
i,0)) / dot(z,
np);
264 subVector(u, 0,
i) -=
step * subVector(r, 0,
i);
268 vigra_precondition(vigra::detail::quadprogAddConstraint(R,
J, d,
i,
R_norm),
269 "quadraticProgramming(): Equality constraints are linearly dependent.");
275 for (
int i = 0;
i <
mi; ++
i)
Class for a single RGB value.
Definition rgbvalue.hxx:128
TinyVectorView< VALUETYPE, TO-FROM > subarray() const
Definition tinyvector.hxx:887
Class for fixed size vectors.
Definition tinyvector.hxx:1008
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:684
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
unsigned int quadraticProgramming(...)
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044