GCC Code Coverage Report

Directory: ./
File: include/ndcurves/quadratic_variable.h
Date: 2025-03-05 17:18:30
Exec Total Coverage
Lines: 64 72 88.9%
Branches: 57 122 46.7%
Line Branch Exec Source
1 /**
2 * \file quadratic_variable.h
3 * \brief storage for variable points of the form p_i = x' A_i x + B_i x + c_i
4 * \author Steve T.
5 * \version 0.1
6 * \date 07/02/2019
7 */
8
9 #ifndef _CLASS_QUADRATIC_VARIABLE
10 #define _CLASS_QUADRATIC_VARIABLE
11
12 #include <math.h>
13
14 #include <Eigen/Core>
15 #include <stdexcept>
16 #include <vector>
17
18 #include "MathDefs.h"
19 #include "ndcurves/curve_abc.h"
20 #include "ndcurves/linear_variable.h"
21
22 namespace ndcurves {
23
24 template <typename Numeric = double>
25 struct quadratic_variable {
26 typedef Eigen::Matrix<Numeric, Eigen::Dynamic, Eigen::Dynamic> matrix_x_t;
27 typedef Eigen::Matrix<Numeric, Eigen::Dynamic, 1> point_t;
28 typedef quadratic_variable<Numeric> quadratic_variable_t;
29
30
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 quadratic_variable() {
31 2 c_ = 0.;
32
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 b_ = point_t::Zero(1);
33
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 A_ = matrix_x_t::Zero(1, 1);
34 2 zero = true;
35 2 }
36
37 66 quadratic_variable(const matrix_x_t& A, const point_t& b, const Numeric c = 0)
38
1/2
✓ Branch 2 taken 66 times.
✗ Branch 3 not taken.
 66 : c_(c), b_(b), A_(A), zero(false) {
39
3/6
✓ Branch 2 taken 66 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 66 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 66 times.
 66 if (A.cols() != b.rows() || A.cols() != A.rows()) {
40 ✗ throw std::invalid_argument("The dimensions of A and b are incorrect.");
41 }
42 66 }
43
44 2 quadratic_variable(const point_t& b, const Numeric c = 0)
45 2 : c_(c),
46 2 b_(b),
47
2/4
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 2 A_(matrix_x_t::Zero((int)(b.rows()), (int)(b.rows()))),
48 2 zero(false) {}
49
50 static quadratic_variable_t Zero(size_t /*dim*/ = 0) {
51 return quadratic_variable_t();
52 }
53
54 // linear evaluation
55 Numeric operator()(const Eigen::Ref<const point_t>& val) const {
56 if (isZero()) {
57 throw std::runtime_error("Not initialized! (isZero)");
58 }
59 return val.transpose() * A() * val + b().transpose() * val + c();
60 }
61
62 32 quadratic_variable& operator+=(const quadratic_variable& w1) {
63
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 32 times.
 32 if (w1.isZero()) return *this;
64
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 32 times.
 32 if (isZero()) {
65 ✗ this->A_ = w1.A_;
66 ✗ this->b_ = w1.b_;
67 ✗ this->c_ = w1.c_;
68 ✗ zero = false;
69 } else {
70 32 this->A_ += w1.A_;
71 32 this->b_ += w1.b_;
72 32 this->c_ += w1.c_;
73 }
74 32 return *this;
75 }
76 quadratic_variable& operator-=(const quadratic_variable& w1) {
77 if (w1.isZero()) return *this;
78 if (isZero()) {
79 this->A_ = -w1.A_;
80 this->b_ = -w1.b_;
81 this->c_ = -w1.c_;
82 zero = false;
83 } else {
84 this->A_ -= w1.A_;
85 this->b_ -= w1.b_;
86 this->c_ -= w1.c_;
87 }
88 return *this;
89 }
90
91 2 quadratic_variable& operator/=(const double d) {
92 // handling zero case
93
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 if (!isZero()) {
94 2 this->A_ /= d;
95 2 this->b_ /= d;
96 2 this->c_ /= d;
97 }
98 2 return *this;
99 }
100 32 quadratic_variable& operator*=(const double d) {
101 // handling zero case
102
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 if (!isZero()) {
103 32 this->A_ *= d;
104 32 this->b_ *= d;
105 32 this->c_ *= d;
106 }
107 32 return *this;
108 }
109
110 34 const matrix_x_t& A() const {
111
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 34 if (isZero()) {
112 ✗ throw std::runtime_error("Not initialized! (isZero)");
113 }
114 34 return A_;
115 }
116 34 const point_t& b() const {
117
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 34 if (isZero()) {
118 ✗ throw std::runtime_error("Not initialized! (isZero)");
119 }
120 34 return b_;
121 }
122 34 const Numeric c() const {
123
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 34 if (isZero()) {
124 ✗ throw std::runtime_error("Not initialized! (isZero)");
125 }
126 34 return c_;
127 }
128 200 bool isZero() const { return zero; }
129 std::size_t size() const {
130 return zero ? 0 : std::max(A_.cols(), (std::max(b_.cols(), c_.size())));
131 }
132
133 private:
134 Numeric c_;
135 point_t b_;
136 matrix_x_t A_;
137 bool zero;
138 };
139
140 /// \brief Transforms a vector into a diagonal matrix
141 template <typename N>
142 64 Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> to_diagonal(
143 const Eigen::Ref<const Eigen::Matrix<N, Eigen::Dynamic, 1> > vec) {
144 typedef typename Eigen::Matrix<N, Eigen::Dynamic, Eigen::Dynamic> matrix_t;
145
2/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
 128 return vec.asDiagonal();
146 matrix_t res(matrix_t::Zero(vec.rows(), vec.rows()));
147 for (int i = 0; i < vec.rows(); ++i) res(i, i) = vec(i);
148 return res;
149 }
150
151 // only works with diagonal linear variables
152 template <typename N>
153 32 inline quadratic_variable<N> operator*(const linear_variable<N>& w1,
154 const linear_variable<N>& w2) {
155 typedef quadratic_variable<N> quad_var_t;
156 typedef linear_variable<N> lin_var_t;
157 typedef typename quad_var_t::matrix_x_t matrix_x_t;
158 typedef typename quad_var_t::point_t point_t;
159 typedef typename lin_var_t::vector_x_t point_dim_t;
160
2/4
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 32 times.
✗ Branch 7 not taken.
 32 point_dim_t ones = point_dim_t::Ones(w1.c().size());
161
6/12
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 32 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 32 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 32 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 32 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 32 times.
✗ Branch 19 not taken.
 32 point_t b1 = w1.B().transpose() * ones, b2 = w2.B().transpose() * ones;
162
2/4
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
 32 matrix_x_t B1 = to_diagonal<N>(b1);
163
2/4
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
 32 matrix_x_t B2 = to_diagonal<N>(b2); // b1.array().square()
164 // omitting all transposes since A matrices are diagonals
165
3/6
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 32 times.
✗ Branch 8 not taken.
 32 matrix_x_t A = B1.transpose() * B2;
166
6/12
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 32 times.
✗ Branch 7 not taken.
✓ Branch 11 taken 32 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 32 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 32 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 32 times.
✗ Branch 21 not taken.
 32 point_t b = w1.c().transpose() * w2.B() + w2.c().transpose() * w1.B();
167
3/6
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 32 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 32 times.
✗ Branch 10 not taken.
 32 N c = w1.c().transpose() * w2.c();
168
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 return quad_var_t(A, b, c);
169 32 }
170
171 template <typename N>
172 inline quadratic_variable<N> operator+(const quadratic_variable<N>& w1,
173 const quadratic_variable<N>& w2) {
174 quadratic_variable<N> res(w1.A(), w1.b(), w1.c());
175 return res += w2;
176 }
177
178 template <typename N>
179 quadratic_variable<N> operator-(const quadratic_variable<N>& w1,
180 const quadratic_variable<N>& w2) {
181 quadratic_variable<N> res(w1.A(), w1.b(), w1.c());
182 return res -= w2;
183 }
184
185 template <typename N>
186 quadratic_variable<N> operator*(const double k,
187 const quadratic_variable<N>& w) {
188 quadratic_variable<N> res(w.A(), w.b(), w.c());
189 return res *= k;
190 }
191
192 template <typename N>
193 32 quadratic_variable<N> operator*(const quadratic_variable<N>& w,
194 const double k) {
195
4/8
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 32 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 32 times.
✗ Branch 11 not taken.
 32 quadratic_variable<N> res(w.A(), w.b(), w.c());
196
2/4
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
 64 return res *= k;
197 32 }
198
199 template <typename N>
200 2 quadratic_variable<N> operator/(const quadratic_variable<N>& w,
201 const double k) {
202
4/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
 2 quadratic_variable<N> res(w.A(), w.b(), w.c());
203
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 return res /= k;
204 2 }
205
206 } // namespace ndcurves
207 #endif //_CLASS_QUADRATIC_VARIABLE
208