Directory: | ./ |
---|---|
File: | include/ndcurves/polynomial.h |
Date: | 2025-03-05 17:18:30 |
Exec | Total | Coverage | |
---|---|---|---|
Lines: | 321 | 336 | 95.5% |
Branches: | 377 | 704 | 53.6% |
Line | Branch | Exec | Source |
---|---|---|---|
1 | /** | ||
2 | * \file polynomial.h | ||
3 | * \brief Definition of a cubic spline. | ||
4 | * \author Steve T. | ||
5 | * \version 0.1 | ||
6 | * \date 06/17/2013 | ||
7 | * | ||
8 | * This file contains definitions for the polynomial struct. | ||
9 | * It allows the creation and evaluation of natural | ||
10 | * smooth splines of arbitrary dimension and order | ||
11 | */ | ||
12 | |||
13 | #ifndef _STRUCT_POLYNOMIAL | ||
14 | #define _STRUCT_POLYNOMIAL | ||
15 | |||
16 | #include <algorithm> | ||
17 | #include <functional> | ||
18 | #include <iostream> | ||
19 | #include <stdexcept> | ||
20 | |||
21 | #include "MathDefs.h" | ||
22 | #include "curve_abc.h" | ||
23 | |||
24 | namespace ndcurves { | ||
25 | /// \class polynomial. | ||
26 | /// \brief Represents a polynomial of an arbitrary order defined on the interval | ||
27 | /// \f$[t_{min}, t_{max}]\f$. It follows the equation :<br> | ||
28 | /// \f$ x(t) = a + b(t - t_{min}) + ... + d(t - t_{min})^N \f$<br> | ||
29 | /// where N is the order and \f$ t \in [t_{min}, t_{max}] \f$. | ||
30 | /// | ||
31 | template <typename Time = double, typename Numeric = Time, bool Safe = false, | ||
32 | typename Point = Eigen::Matrix<Numeric, Eigen::Dynamic, 1>, | ||
33 | typename T_Point = | ||
34 | std::vector<Point, Eigen::aligned_allocator<Point> > > | ||
35 | struct polynomial : public curve_abc<Time, Numeric, Safe, Point> { | ||
36 | typedef Point point_t; | ||
37 | typedef T_Point t_point_t; | ||
38 | typedef Time time_t; | ||
39 | typedef Numeric num_t; | ||
40 | typedef curve_abc<Time, Numeric, Safe, Point> curve_abc_t; | ||
41 | typedef Eigen::MatrixXd coeff_t; | ||
42 | typedef Eigen::Ref<coeff_t> coeff_t_ref; | ||
43 | typedef polynomial<Time, Numeric, Safe, Point, T_Point> polynomial_t; | ||
44 | typedef typename curve_abc_t::curve_ptr_t curve_ptr_t; | ||
45 | |||
46 | /* Constructors - destructors */ | ||
47 | public: | ||
48 | /// \brief Empty constructor. Curve obtained this way can not perform other | ||
49 | /// class functions. | ||
50 | /// | ||
51 |
1/2✓ Branch 2 taken 68 times.
✗ Branch 3 not taken.
|
135 | polynomial() : curve_abc_t(), dim_(0), degree_(0), T_min_(0), T_max_(1.0) {} |
52 | |||
53 | /// \brief Constructor. | ||
54 | /// \param coefficients : a reference to an Eigen matrix where each column is | ||
55 | /// a coefficient, from the zero order coefficient, up to the highest order. | ||
56 | /// Spline order is given by the number of the columns -1. \param min : LOWER | ||
57 | /// bound on interval definition of the curve. \param max : UPPER bound on | ||
58 | /// interval definition of the curve. | ||
59 | 112 | polynomial(const coeff_t& coefficients, const time_t min, const time_t max) | |
60 | : curve_abc_t(), | ||
61 | 224 | dim_(coefficients.rows()), | |
62 |
1/2✓ Branch 1 taken 69 times.
✗ Branch 2 not taken.
|
112 | coefficients_(coefficients), |
63 | 112 | degree_(coefficients.cols() - 1), | |
64 | 112 | T_min_(min), | |
65 | 112 | T_max_(max) { | |
66 |
2/2✓ Branch 1 taken 67 times.
✓ Branch 2 taken 2 times.
|
112 | safe_check(); |
67 | 115 | } | |
68 | |||
69 | /// \brief Constructor | ||
70 | /// \param coefficients : a container containing all coefficients of the | ||
71 | /// spline, starting | ||
72 | /// with the zero order coefficient, up to the highest order. Spline order is | ||
73 | /// given by the size of the coefficients. | ||
74 | /// \param min : LOWER bound on interval definition of the spline. | ||
75 | /// \param max : UPPER bound on interval definition of the spline. | ||
76 | 42 | polynomial(const T_Point& coefficients, const time_t min, const time_t max) | |
77 | : curve_abc_t(), | ||
78 | 42 | dim_(coefficients.begin()->size()), | |
79 |
1/2✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
|
42 | coefficients_(init_coeffs(coefficients.begin(), coefficients.end())), |
80 | 42 | degree_(coefficients_.cols() - 1), | |
81 | 42 | T_min_(min), | |
82 | 42 | T_max_(max) { | |
83 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
42 | safe_check(); |
84 | 42 | } | |
85 | |||
86 | /// \brief Constructor. | ||
87 | /// \param zeroOrderCoefficient : an iterator pointing to the first element of | ||
88 | /// a structure containing the coefficients | ||
89 | /// it corresponds to the zero degree coefficient. | ||
90 | /// \param out : an iterator pointing to the last element of a structure | ||
91 | /// ofcoefficients. \param min : LOWER bound on interval definition of the | ||
92 | /// spline. \param max : UPPER bound on interval definition of the spline. | ||
93 | template <typename In> | ||
94 | 163 | polynomial(In zeroOrderCoefficient, In out, const time_t min, | |
95 | const time_t max) | ||
96 | : curve_abc_t(), | ||
97 | 163 | dim_(zeroOrderCoefficient->size()), | |
98 |
1/2✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
|
163 | coefficients_(init_coeffs(zeroOrderCoefficient, out)), |
99 | 163 | degree_(coefficients_.cols() - 1), | |
100 | 163 | T_min_(min), | |
101 | 163 | T_max_(max) { | |
102 |
1/2✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
|
163 | safe_check(); |
103 | 163 | } | |
104 | |||
105 | /// | ||
106 | /// \brief Constructor from boundary condition with C0 : create a polynomial | ||
107 | /// that connect exactly init and end (order 1) \param init the initial point | ||
108 | /// of the curve \param end the final point of the curve \param min : LOWER | ||
109 | /// bound on interval definition of the spline. \param max : UPPER bound on | ||
110 | /// interval definition of the spline. | ||
111 | /// | ||
112 | 80 | polynomial(const Point& init, const Point& end, const time_t min, | |
113 | const time_t max) | ||
114 |
1/2✓ Branch 3 taken 40 times.
✗ Branch 4 not taken.
|
80 | : dim_(init.size()), degree_(1), T_min_(min), T_max_(max) { |
115 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 36 times.
|
80 | if (T_min_ >= T_max_) |
116 |
1/2✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
8 | throw std::invalid_argument("T_min must be strictly lower than T_max"); |
117 |
2/2✓ Branch 2 taken 1 times.
✓ Branch 3 taken 35 times.
|
72 | if (init.size() != end.size()) |
118 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | throw std::invalid_argument( |
119 | "init and end points must have the same dimensions."); | ||
120 |
1/2✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
|
70 | t_point_t coeffs; |
121 |
1/2✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
|
70 | coeffs.push_back(init); |
122 |
4/8✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 35 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 35 times.
✗ Branch 11 not taken.
|
70 | coeffs.push_back((end - init) / (max - min)); |
123 |
1/2✓ Branch 3 taken 35 times.
✗ Branch 4 not taken.
|
70 | coefficients_ = init_coeffs(coeffs.begin(), coeffs.end()); |
124 |
1/2✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
|
70 | safe_check(); |
125 | 90 | } | |
126 | |||
127 | /// | ||
128 | /// \brief Constructor from boundary condition with C1 : | ||
129 | /// create a polynomial that connect exactly init and end and thier first | ||
130 | /// order derivatives(order 3) \param init the initial point of the curve | ||
131 | /// \param d_init the initial value of the derivative of the curve | ||
132 | /// \param end the final point of the curve | ||
133 | /// \param d_end the final value of the derivative of the curve | ||
134 | /// \param min : LOWER bound on interval definition of the spline. | ||
135 | /// \param max : UPPER bound on interval definition of the spline. | ||
136 | /// | ||
137 | 38 | polynomial(const Point& init, const Point& d_init, const Point& end, | |
138 | const Point& d_end, const time_t min, const time_t max) | ||
139 |
1/2✓ Branch 3 taken 19 times.
✗ Branch 4 not taken.
|
38 | : dim_(init.size()), degree_(3), T_min_(min), T_max_(max) { |
140 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 16 times.
|
38 | if (T_min_ >= T_max_) |
141 |
1/2✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
6 | throw std::invalid_argument("T_min must be strictly lower than T_max"); |
142 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
|
32 | if (init.size() != end.size()) |
143 | ✗ | throw std::invalid_argument( | |
144 | "init and end points must have the same dimensions."); | ||
145 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
|
32 | if (init.size() != d_init.size()) |
146 | ✗ | throw std::invalid_argument( | |
147 | "init and d_init points must have the same dimensions."); | ||
148 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
|
32 | if (init.size() != d_end.size()) |
149 | ✗ | throw std::invalid_argument( | |
150 | "init and d_end points must have the same dimensions."); | ||
151 | /* the coefficients [c0 c1 c2 c3] are found by solving the following system | ||
152 | of equation (found from the boundary conditions) : [1 0 0 0 ] [c0] | ||
153 | [ init ] [1 T T^2 T^3 ] x [c1] = [ end ] [0 1 0 0 ] [c2] [d_init] | ||
154 | [0 1 2T 3T^2] [c3] [d_end ] | ||
155 | */ | ||
156 | 32 | double T = max - min; | |
157 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
32 | Eigen::Matrix<double, 4, 4> m; |
158 |
15/30✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 16 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 16 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 16 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 16 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 16 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 16 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 16 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 16 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 16 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 16 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 16 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 16 times.
✗ Branch 44 not taken.
|
32 | m << 1., 0, 0, 0, 1., T, T * T, T * T * T, 0, 1., 0, 0, 0, 1., 2. * T, |
159 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
32 | 3. * T * T; |
160 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
|
32 | Eigen::Matrix<double, 4, 4> m_inv = m.inverse(); |
161 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
32 | Eigen::Matrix<double, 4, 1> bc; // boundary condition vector |
162 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
32 | coefficients_ = coeff_t::Zero( |
163 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
32 | dim_, degree_ + 1); // init coefficient matrix with the right size |
164 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 16 times.
|
128 | for (size_t i = 0; i < dim_; |
165 | ++i) { // for each dimension, solve the boundary condition problem : | ||
166 |
2/4✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
|
96 | bc[0] = init[i]; |
167 |
2/4✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
|
96 | bc[1] = end[i]; |
168 |
2/4✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
|
96 | bc[2] = d_init[i]; |
169 |
2/4✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
|
96 | bc[3] = d_end[i]; |
170 |
4/8✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 48 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 48 times.
✗ Branch 11 not taken.
|
96 | coefficients_.row(i) = (m_inv * bc).transpose(); |
171 | } | ||
172 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
32 | safe_check(); |
173 | 44 | } | |
174 | |||
175 | /// | ||
176 | /// \brief Constructor from boundary condition with C2 : | ||
177 | /// create a polynomial that connect exactly init and end and thier first and | ||
178 | /// second order derivatives(order 5) \param init the initial point of the | ||
179 | /// curve \param d_init the initial value of the derivative of the curve | ||
180 | /// \param d_init the initial value of the second derivative of the curve | ||
181 | /// \param end the final point of the curve | ||
182 | /// \param d_end the final value of the derivative of the curve | ||
183 | /// \param d_end the final value of the second derivative of the curve | ||
184 | /// \param min : LOWER bound on interval definition of the spline. | ||
185 | /// \param max : UPPER bound on interval definition of the spline. | ||
186 | /// | ||
187 | 7882 | polynomial(const Point& init, const Point& d_init, const Point& dd_init, | |
188 | const Point& end, const Point& d_end, const Point& dd_end, | ||
189 | const time_t min, const time_t max) | ||
190 |
1/2✓ Branch 3 taken 7870 times.
✗ Branch 4 not taken.
|
7882 | : dim_(init.size()), degree_(5), T_min_(min), T_max_(max) { |
191 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 7867 times.
|
7882 | if (T_min_ >= T_max_) |
192 |
1/2✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
5 | throw std::invalid_argument("T_min must be strictly lower than T_max"); |
193 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 7867 times.
|
7877 | if (init.size() != end.size()) |
194 | ✗ | throw std::invalid_argument( | |
195 | "init and end points must have the same dimensions."); | ||
196 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 7867 times.
|
7877 | if (init.size() != d_init.size()) |
197 | ✗ | throw std::invalid_argument( | |
198 | "init and d_init points must have the same dimensions."); | ||
199 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 7867 times.
|
7877 | if (init.size() != d_end.size()) |
200 | ✗ | throw std::invalid_argument( | |
201 | "init and d_end points must have the same dimensions."); | ||
202 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 7867 times.
|
7877 | if (init.size() != dd_init.size()) |
203 | ✗ | throw std::invalid_argument( | |
204 | "init and dd_init points must have the same dimensions."); | ||
205 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 7867 times.
|
7877 | if (init.size() != dd_end.size()) |
206 | ✗ | throw std::invalid_argument( | |
207 | "init and dd_end points must have the same dimensions."); | ||
208 | /* the coefficients [c0 c1 c2 c3 c4 c5] are found by solving the following | ||
209 | system of equation (found from the boundary conditions) : [1 0 0 0 0 | ||
210 | 0 ] [c0] [ init ] [1 T T^2 T^3 T^4 T^5 ] [c1] [ end ] [0 | ||
211 | 1 0 0 0 0 ] [c2] [d_init ] [0 1 2T 3T^2 4T^3 5T^4 ] x | ||
212 | [c3] = [d_end ] [0 0 2 0 0 0 ] [c4] [dd_init] [0 0 2 6T | ||
213 | 12T^2 20T^3] [c5] [dd_end ] | ||
214 | */ | ||
215 | 7877 | double T = max - min; | |
216 |
1/2✓ Branch 1 taken 7867 times.
✗ Branch 2 not taken.
|
7877 | Eigen::Matrix<double, 6, 6> m; |
217 |
13/26✓ Branch 1 taken 7867 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7867 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7867 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 7867 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 7867 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 7867 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 7867 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 7867 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 7867 times.
✗ Branch 26 not taken.
✓ Branch 29 taken 7867 times.
✗ Branch 30 not taken.
✓ Branch 33 taken 7867 times.
✗ Branch 34 not taken.
✓ Branch 37 taken 7867 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 7867 times.
✗ Branch 41 not taken.
|
7877 | m << 1., 0, 0, 0, 0, 0, 1., T, T * T, pow(T, 3), pow(T, 4), pow(T, 5), 0, |
218 |
10/20✓ Branch 1 taken 7867 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7867 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7867 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 7867 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 7867 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 7867 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 7867 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 7867 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 7867 times.
✗ Branch 26 not taken.
✓ Branch 29 taken 7867 times.
✗ Branch 30 not taken.
|
7877 | 1., 0, 0, 0, 0, 0, 1., 2. * T, 3. * T * T, 4. * pow(T, 3), |
219 |
12/24✓ Branch 2 taken 7867 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 7867 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 7867 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 7867 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 7867 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 7867 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 7867 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 7867 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 7867 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 7867 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 7867 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 7867 times.
✗ Branch 36 not taken.
|
7877 | 5. * pow(T, 4), 0, 0, 2, 0, 0, 0, 0, 0, 2, 6. * T, 12. * T * T, |
220 |
1/2✓ Branch 2 taken 7867 times.
✗ Branch 3 not taken.
|
7877 | 20. * pow(T, 3); |
221 |
2/4✓ Branch 1 taken 7867 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7867 times.
✗ Branch 5 not taken.
|
7877 | Eigen::Matrix<double, 6, 6> m_inv = m.inverse(); |
222 |
1/2✓ Branch 1 taken 7867 times.
✗ Branch 2 not taken.
|
7877 | Eigen::Matrix<double, 6, 1> bc; // boundary condition vector |
223 |
1/2✓ Branch 1 taken 7867 times.
✗ Branch 2 not taken.
|
7877 | coefficients_ = coeff_t::Zero( |
224 |
1/2✓ Branch 1 taken 7867 times.
✗ Branch 2 not taken.
|
7877 | dim_, degree_ + 1); // init coefficient matrix with the right size |
225 |
2/2✓ Branch 0 taken 23601 times.
✓ Branch 1 taken 7867 times.
|
31508 | for (size_t i = 0; i < dim_; |
226 | ++i) { // for each dimension, solve the boundary condition problem : | ||
227 |
2/4✓ Branch 1 taken 23601 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23601 times.
✗ Branch 5 not taken.
|
23631 | bc[0] = init[i]; |
228 |
2/4✓ Branch 1 taken 23601 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23601 times.
✗ Branch 5 not taken.
|
23631 | bc[1] = end[i]; |
229 |
2/4✓ Branch 1 taken 23601 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23601 times.
✗ Branch 5 not taken.
|
23631 | bc[2] = d_init[i]; |
230 |
2/4✓ Branch 1 taken 23601 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23601 times.
✗ Branch 5 not taken.
|
23631 | bc[3] = d_end[i]; |
231 |
2/4✓ Branch 1 taken 23601 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23601 times.
✗ Branch 5 not taken.
|
23631 | bc[4] = dd_init[i]; |
232 |
2/4✓ Branch 1 taken 23601 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23601 times.
✗ Branch 5 not taken.
|
23631 | bc[5] = dd_end[i]; |
233 |
4/8✓ Branch 1 taken 23601 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23601 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 23601 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 23601 times.
✗ Branch 11 not taken.
|
23631 | coefficients_.row(i) = (m_inv * bc).transpose(); |
234 | } | ||
235 |
1/2✓ Branch 1 taken 7867 times.
✗ Branch 2 not taken.
|
7877 | safe_check(); |
236 | 7887 | } | |
237 | |||
238 | /// \brief Destructor | ||
239 | 32688 | virtual ~polynomial() {} | |
240 | |||
241 | 16200 | polynomial(const polynomial& other) | |
242 | 16200 | : dim_(other.dim_), | |
243 | 16200 | coefficients_(other.coefficients_), | |
244 | 16200 | degree_(other.degree_), | |
245 | 16200 | T_min_(other.T_min_), | |
246 |
1/2✓ Branch 2 taken 8114 times.
✗ Branch 3 not taken.
|
16200 | T_max_(other.T_max_) {} |
247 | |||
248 | // polynomial& operator=(const polynomial& other); | ||
249 | |||
250 | /** | ||
251 | * @brief MinimumJerk Build a polynomial curve connecting p_init to p_final | ||
252 | * minimizing the time integral of the squared jerk with a zero initial and | ||
253 | * final velocity and acceleration | ||
254 | * @param p_init the initial point | ||
255 | * @param p_final the final point | ||
256 | * @param t_min initial time | ||
257 | * @param t_max final time | ||
258 | * @return the polynomial curve | ||
259 | */ | ||
260 | 33 | static polynomial_t MinimumJerk(const point_t& p_init, const point_t& p_final, | |
261 | const time_t t_min = 0., | ||
262 | const time_t t_max = 1.) { | ||
263 |
3/4✓ Branch 3 taken 28 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 26 times.
✓ Branch 7 taken 2 times.
|
39 | polynomial_t out = |
264 | polynomial(coeff_t::Zero(p_init.size(), 6), t_min, t_max); | ||
265 |
2/2✓ Branch 1 taken 24 times.
✓ Branch 2 taken 2 times.
|
30 | MinimumJerk(out, p_init, p_final, t_min, t_max); |
266 | 27 | return out; | |
267 | 3 | } | |
268 | |||
269 | /** | ||
270 | * @brief MinimumJerk Build a polynomial curve connecting p_init to p_final | ||
271 | * minimizing the time integral of the squared jerk with a zero initial and | ||
272 | * final velocity and acceleration. | ||
273 | * @param out The output polynomial needs to be of the correct size. | ||
274 | * @param p_init the initial point | ||
275 | * @param p_final the final point | ||
276 | * @param t_min initial time | ||
277 | * @param t_max final time | ||
278 | * @return the polynomial curve | ||
279 | */ | ||
280 | 37 | static void MinimumJerk(polynomial_t& out, const point_t& p_init, | |
281 | const point_t& p_final, const time_t t_min = 0., | ||
282 | const time_t t_max = 1.) { | ||
283 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 33 times.
|
37 | if (t_min > t_max) |
284 | ✗ | throw std::invalid_argument( | |
285 | "final time should be superior or equal to initial time."); | ||
286 | 37 | const size_t dim(p_init.size()); | |
287 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 31 times.
|
37 | if (static_cast<size_t>(p_final.size()) != dim) |
288 |
1/2✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
3 | throw std::invalid_argument( |
289 | "Initial and final points must have the same dimension."); | ||
290 | 34 | const double T = t_max - t_min; | |
291 | 34 | const double T2 = T * T; | |
292 | 34 | const double T3 = T2 * T; | |
293 | 34 | const double T4 = T3 * T; | |
294 | 34 | const double T5 = T4 * T; | |
295 | |||
296 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 31 times.
|
34 | assert(out.coefficients_.cols() == 6); |
297 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 31 times.
|
34 | assert(out.coefficients_.rows() == static_cast<Eigen::Index>(dim)); |
298 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 31 times.
|
34 | assert(out.dim_ == dim); |
299 |
1/2✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
|
34 | out.coefficients_.fill(0.0); |
300 |
2/4✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 31 times.
✗ Branch 5 not taken.
|
34 | out.coefficients_.col(0) = p_init; |
301 |
5/10✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 31 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 31 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 31 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 31 times.
✗ Branch 14 not taken.
|
34 | out.coefficients_.col(3) = 10 * (p_final - p_init) / T3; |
302 |
5/10✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 31 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 31 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 31 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 31 times.
✗ Branch 14 not taken.
|
34 | out.coefficients_.col(4) = -15 * (p_final - p_init) / T4; |
303 |
5/10✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 31 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 31 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 31 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 31 times.
✗ Branch 14 not taken.
|
34 | out.coefficients_.col(5) = 6 * (p_final - p_init) / T5; |
304 | 34 | out.degree_ = 5; | |
305 | 34 | out.T_min_ = t_min; | |
306 | 34 | out.T_max_ = t_max; | |
307 |
1/2✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
|
34 | out.safe_check(); |
308 | 34 | } | |
309 | |||
310 | private: | ||
311 | 16214 | void safe_check() { | |
312 | if (Safe) { | ||
313 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8113 times.
|
16200 | if (T_min_ > T_max_) { |
314 |
1/2✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
3 | throw std::invalid_argument("Tmin should be inferior to Tmax"); |
315 | } | ||
316 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 8113 times.
|
16197 | if (coefficients_.cols() != int(degree_ + 1)) { |
317 | ✗ | throw std::runtime_error("Spline order and coefficients do not match"); | |
318 | } | ||
319 | } | ||
320 | 16211 | } | |
321 | |||
322 | /* Constructors - destructors */ | ||
323 | |||
324 | /*Operations*/ | ||
325 | public: | ||
326 | /// \brief Evaluation of the cubic spline at time t using horner's scheme. | ||
327 | /// \param t : time when to evaluate the spline. | ||
328 | /// \return \f$x(t)\f$ point corresponding on spline at time t. | ||
329 | 217612 | virtual point_t operator()(const time_t t) const { | |
330 |
1/2✓ Branch 1 taken 111132 times.
✗ Branch 2 not taken.
|
217612 | check_if_not_empty(); |
331 |
4/4✓ Branch 0 taken 110123 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 110120 times.
|
215604 | if ((t < T_min_ || t > T_max_) && Safe) { |
332 |
1/2✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
|
16 | throw std::invalid_argument( |
333 | "error in polynomial : time t to evaluate should be in range [Tmin, " | ||
334 | "Tmax] of the curve"); | ||
335 | } | ||
336 | 217596 | time_t const dt(t - T_min_); | |
337 |
2/4✓ Branch 1 taken 111124 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 111124 times.
✗ Branch 5 not taken.
|
217596 | point_t h = coefficients_.col(degree_); |
338 |
2/2✓ Branch 0 taken 503288 times.
✓ Branch 1 taken 111124 times.
|
1207679 | for (int i = (int)(degree_ - 1); i >= 0; i--) { |
339 |
4/8✓ Branch 1 taken 503288 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 503288 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 503288 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 503288 times.
✗ Branch 11 not taken.
|
990083 | h = dt * h + coefficients_.col(i); |
340 | } | ||
341 | 435192 | return h; | |
342 | } | ||
343 | |||
344 | /** | ||
345 | * @brief isApprox check if other and *this are approximately equals. | ||
346 | * Only two curves of the same class can be approximately equals, for | ||
347 | * comparison between different type of curves see isEquivalent | ||
348 | * @param other the other curve to check | ||
349 | * @param prec the precision threshold, default | ||
350 | * Eigen::NumTraits<Numeric>::dummy_precision() | ||
351 | * @return true is the two curves are approximately equals | ||
352 | */ | ||
353 | 26 | bool isApprox( | |
354 | const polynomial_t& other, | ||
355 | const Numeric prec = Eigen::NumTraits<Numeric>::dummy_precision()) const { | ||
356 |
1/2✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
|
52 | return ndcurves::isApprox<num_t>(T_min_, other.min()) && |
357 | 26 | ndcurves::isApprox<num_t>(T_max_, other.max()) && | |
358 |
4/6✓ Branch 0 taken 13 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 12 times.
✓ Branch 7 taken 1 times.
|
76 | dim_ == other.dim() && degree_ == other.degree() && |
359 |
2/2✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
|
50 | coefficients_.isApprox(other.coefficients_, prec); |
360 | } | ||
361 | |||
362 | 12 | virtual bool isApprox( | |
363 | const curve_abc_t* other, | ||
364 | const Numeric prec = Eigen::NumTraits<Numeric>::dummy_precision()) const { | ||
365 |
1/2✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
|
12 | const polynomial_t* other_cast = dynamic_cast<const polynomial_t*>(other); |
366 |
1/2✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
|
12 | if (other_cast) |
367 | 12 | return isApprox(*other_cast, prec); | |
368 | else | ||
369 | ✗ | return false; | |
370 | } | ||
371 | |||
372 | 14 | virtual bool operator==(const polynomial_t& other) const { | |
373 | 14 | return isApprox(other); | |
374 | } | ||
375 | |||
376 | 6 | virtual bool operator!=(const polynomial_t& other) const { | |
377 | 6 | return !(*this == other); | |
378 | } | ||
379 | |||
380 | /// \brief Evaluation of the derivative of order N of spline at time t. | ||
381 | /// \param t : the time when to evaluate the spline. | ||
382 | /// \param order : order of derivative. | ||
383 | /// \return \f$\frac{d^Nx(t)}{dt^N}\f$ point corresponding on derivative | ||
384 | /// spline at time t. | ||
385 | 59565 | virtual point_t derivate(const time_t t, const std::size_t order) const { | |
386 |
1/2✓ Branch 1 taken 29787 times.
✗ Branch 2 not taken.
|
59565 | check_if_not_empty(); |
387 |
4/4✓ Branch 0 taken 29785 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 29783 times.
|
59565 | if ((t < T_min_ || t > T_max_) && Safe) { |
388 |
1/2✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
8 | throw std::invalid_argument( |
389 | "error in polynomial : time t to evaluate derivative should be in " | ||
390 | "range [Tmin, Tmax] of the curve"); | ||
391 | } | ||
392 | 59557 | time_t const dt(t - T_min_); | |
393 | 59557 | time_t cdt(1); | |
394 |
2/4✓ Branch 1 taken 29783 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 29783 times.
✗ Branch 5 not taken.
|
59557 | point_t currentPoint_ = point_t::Zero(dim_); |
395 |
2/2✓ Branch 0 taken 144881 times.
✓ Branch 1 taken 29783 times.
|
349283 | for (int i = (int)(order); i < (int)(degree_ + 1); ++i, cdt *= dt) { |
396 |
4/8✓ Branch 2 taken 144881 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 144881 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 144881 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 144881 times.
✗ Branch 12 not taken.
|
289726 | currentPoint_ += cdt * coefficients_.col(i) * fact(i, order); |
397 | } | ||
398 | 119114 | return currentPoint_; | |
399 | } | ||
400 | |||
401 | 34 | polynomial_t compute_derivate(const std::size_t order) const { | |
402 |
1/2✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
|
34 | check_if_not_empty(); |
403 |
2/2✓ Branch 0 taken 7 times.
✓ Branch 1 taken 10 times.
|
34 | if (order == 0) { |
404 |
1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
|
14 | return *this; |
405 | } | ||
406 |
2/4✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
|
20 | coeff_t coeff_derivated = deriv_coeff(coefficients_); |
407 |
1/2✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
|
20 | polynomial_t deriv(coeff_derivated, T_min_, T_max_); |
408 |
1/2✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
|
20 | return deriv.compute_derivate(order - 1); |
409 | } | ||
410 | |||
411 | /// \brief Compute the derived curve at order N. | ||
412 | /// \param order : order of derivative. | ||
413 | /// \return A pointer to \f$\frac{d^Nx(t)}{dt^N}\f$ derivative order N of the | ||
414 | /// curve. | ||
415 | 12 | polynomial_t* compute_derivate_ptr(const std::size_t order) const { | |
416 |
1/2✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
|
12 | return new polynomial_t(compute_derivate(order)); |
417 | } | ||
418 | |||
419 | 155 | Eigen::MatrixXd coeff() const { return coefficients_; } | |
420 | |||
421 | 2 | point_t coeffAtDegree(const std::size_t degree) const { | |
422 | 2 | point_t res; | |
423 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if (degree <= degree_) { |
424 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
2 | res = coefficients_.col(degree); |
425 | } | ||
426 | 2 | return res; | |
427 | } | ||
428 | |||
429 | private: | ||
430 | 289726 | num_t fact(const std::size_t n, const std::size_t order) const { | |
431 | 289726 | num_t res(1); | |
432 |
2/2✓ Branch 0 taken 146834 times.
✓ Branch 1 taken 144881 times.
|
583342 | for (std::size_t i = 0; i < std::size_t(order); ++i) { |
433 | 293616 | res *= (num_t)(n - i); | |
434 | } | ||
435 | 289726 | return res; | |
436 | } | ||
437 | |||
438 | 20 | coeff_t deriv_coeff(coeff_t coeff) const { | |
439 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 8 times.
|
20 | if (coeff.cols() == 1) // only the constant part is left, fill with 0 |
440 |
2/4✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
|
4 | return coeff_t::Zero(coeff.rows(), 1); |
441 |
1/2✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
|
16 | coeff_t coeff_derivated(coeff.rows(), coeff.cols() - 1); |
442 |
2/2✓ Branch 1 taken 13 times.
✓ Branch 2 taken 8 times.
|
42 | for (std::size_t i = 0; i < std::size_t(coeff_derivated.cols()); i++) { |
443 |
4/8✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 13 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 13 times.
✗ Branch 11 not taken.
|
26 | coeff_derivated.col(i) = coeff.col(i + 1) * (num_t)(i + 1); |
444 | } | ||
445 | 16 | return coeff_derivated; | |
446 | } | ||
447 | |||
448 | 277211 | void check_if_not_empty() const { | |
449 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 140936 times.
|
277211 | if (coefficients_.size() == 0) { |
450 | ✗ | throw std::runtime_error( | |
451 | "Error in polynomial : there is no coefficients set / did you use " | ||
452 | "empty constructor ?"); | ||
453 | } | ||
454 | 277211 | } | |
455 | /*Operations*/ | ||
456 | |||
457 | public: | ||
458 | /*Helpers*/ | ||
459 | /// \brief Get dimension of curve. | ||
460 | /// \return dimension of curve. | ||
461 | 24414 | std::size_t virtual dim() const { return dim_; }; | |
462 | /// \brief Get the minimum time for which the curve is defined | ||
463 | /// \return \f$t_{min}\f$ lower bound of time range. | ||
464 | 16471 | num_t virtual min() const { return T_min_; } | |
465 | /// \brief Get the maximum time for which the curve is defined. | ||
466 | /// \return \f$t_{max}\f$ upper bound of time range. | ||
467 | 16679 | num_t virtual max() const { return T_max_; } | |
468 | /// \brief Get the degree of the curve. | ||
469 | /// \return \f$degree\f$, the degree of the curve. | ||
470 | 135 | virtual std::size_t degree() const { return degree_; } | |
471 | /*Helpers*/ | ||
472 | |||
473 | 8 | polynomial_t& operator+=(const polynomial_t& p1) { | |
474 | 8 | assert_operator_compatible(p1); | |
475 |
2/2✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
|
5 | if (p1.degree() > degree()) { |
476 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | polynomial_t::coeff_t res = p1.coeff(); |
477 |
2/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
1 | res.block(0, 0, coefficients_.rows(), coefficients_.cols()) += |
478 | 1 | coefficients_; | |
479 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | coefficients_ = res; |
480 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | degree_ = p1.degree(); |
481 | 1 | } else { | |
482 |
4/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
|
4 | coefficients_.block(0, 0, p1.coeff().rows(), p1.coeff().cols()) += |
483 | p1.coeff(); | ||
484 | } | ||
485 | 5 | return *this; | |
486 | } | ||
487 | |||
488 | 8 | polynomial_t& operator-=(const polynomial_t& p1) { | |
489 | 8 | assert_operator_compatible(p1); | |
490 |
2/2✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
|
5 | if (p1.degree() > degree()) { |
491 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
1 | polynomial_t::coeff_t res = -p1.coeff(); |
492 |
2/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
1 | res.block(0, 0, coefficients_.rows(), coefficients_.cols()) += |
493 | 1 | coefficients_; | |
494 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | coefficients_ = res; |
495 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | degree_ = p1.degree(); |
496 | 1 | } else { | |
497 |
4/8✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
|
4 | coefficients_.block(0, 0, p1.coeff().rows(), p1.coeff().cols()) -= |
498 | p1.coeff(); | ||
499 | } | ||
500 | 5 | return *this; | |
501 | } | ||
502 | |||
503 | 5 | polynomial_t& operator+=(const polynomial_t::point_t& point) { | |
504 |
1/2✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
|
5 | coefficients_.col(0) += point; |
505 | 5 | return *this; | |
506 | } | ||
507 | |||
508 | 3 | polynomial_t& operator-=(const polynomial_t::point_t& point) { | |
509 |
1/2✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
3 | coefficients_.col(0) -= point; |
510 | 3 | return *this; | |
511 | } | ||
512 | |||
513 | 3 | polynomial_t& operator/=(const double d) { | |
514 | 3 | coefficients_ /= d; | |
515 | 3 | return *this; | |
516 | } | ||
517 | |||
518 | 4 | polynomial_t& operator*=(const double d) { | |
519 | 4 | coefficients_ *= d; | |
520 | 4 | return *this; | |
521 | } | ||
522 | |||
523 | /// \brief Compute the cross product of the current polynomial by another | ||
524 | /// polynomial. | ||
525 | /// The cross product p1Xp2 of 2 polynomials p1 and p2 is defined such that | ||
526 | /// forall t, p1Xp2(t) = p1(t) X p2(t), with X designing the cross product. | ||
527 | /// This method of course only makes sense for dimension 3 polynomials. | ||
528 | /// \param pOther other polynomial to compute the cross product with. | ||
529 | /// \return a new polynomial defining the cross product between this and | ||
530 | /// pother | ||
531 | 9 | polynomial_t cross(const polynomial_t& pOther) const { | |
532 |
2/2✓ Branch 1 taken 4 times.
✓ Branch 2 taken 5 times.
|
9 | assert_operator_compatible(pOther); |
533 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
|
4 | if (dim() != 3) |
534 | ✗ | throw std::invalid_argument( | |
535 | "Can't perform cross product on polynomials with dimensions != 3 "); | ||
536 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
4 | std::size_t new_degree = degree() + pOther.degree(); |
537 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
4 | coeff_t nCoeffs = Eigen::MatrixXd::Zero(3, new_degree + 1); |
538 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | Eigen::Vector3d currentVec; |
539 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | Eigen::Vector3d currentVecCrossed; |
540 |
2/2✓ Branch 1 taken 17 times.
✓ Branch 2 taken 4 times.
|
21 | for (long i = 0; i < coefficients_.cols(); ++i) { |
541 |
2/4✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 17 times.
✗ Branch 5 not taken.
|
17 | currentVec = coefficients_.col(i); |
542 |
3/4✓ Branch 1 taken 71 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 54 times.
✓ Branch 6 taken 17 times.
|
71 | for (long j = 0; j < pOther.coeff().cols(); ++j) { |
543 |
3/6✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 54 times.
✗ Branch 8 not taken.
|
54 | currentVecCrossed = pOther.coeff().col(j); |
544 |
3/6✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 54 times.
✗ Branch 8 not taken.
|
54 | nCoeffs.col(i + j) += currentVec.cross(currentVecCrossed); |
545 | } | ||
546 | } | ||
547 | // remove last degrees is they are equal to 0 | ||
548 | 4 | long final_degree = new_degree; | |
549 |
8/10✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 13 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 11 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 11 times.
✓ Branch 11 taken 4 times.
|
15 | while (nCoeffs.col(final_degree).norm() <= ndcurves::MARGIN && |
550 | final_degree > 0) { | ||
551 | 11 | --final_degree; | |
552 | } | ||
553 |
5/10✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
|
8 | return polynomial_t(nCoeffs.leftCols(final_degree + 1), min(), max()); |
554 | 4 | } | |
555 | |||
556 | /// \brief Compute the cross product of the current polynomial p by a point | ||
557 | /// point. | ||
558 | /// The cross product pXpoint of is defined such that | ||
559 | /// forall t, pXpoint(t) = p(t) X point, with X designing the cross product. | ||
560 | /// This method of course only makes sense for dimension 3 polynomials. | ||
561 | /// \param point point to compute the cross product with. | ||
562 | /// \return a new polynomial defining the cross product between this and | ||
563 | /// point | ||
564 | 2 | polynomial_t cross(const polynomial_t::point_t& point) const { | |
565 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
|
2 | if (dim() != 3) |
566 | ✗ | throw std::invalid_argument( | |
567 | "Can't perform cross product on polynomials with dimensions != 3 "); | ||
568 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | coeff_t nCoeffs = coefficients_; |
569 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | Eigen::Vector3d currentVec; |
570 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | Eigen::Vector3d pointVec = point; |
571 |
2/2✓ Branch 1 taken 7 times.
✓ Branch 2 taken 2 times.
|
9 | for (long i = 0; i < coefficients_.cols(); ++i) { |
572 |
2/4✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
|
7 | currentVec = coefficients_.col(i); |
573 |
3/6✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
|
7 | nCoeffs.col(i) = currentVec.cross(pointVec); |
574 | } | ||
575 | // remove last degrees is they are equal to 0 | ||
576 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | long final_degree = degree(); |
577 |
4/10✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
|
2 | while (nCoeffs.col(final_degree).norm() <= ndcurves::MARGIN && |
578 | final_degree > 0) { | ||
579 | ✗ | --final_degree; | |
580 | } | ||
581 |
5/10✓ 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.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
|
4 | return polynomial_t(nCoeffs.leftCols(final_degree + 1), min(), max()); |
582 | 2 | } | |
583 | |||
584 | /*Attributes*/ | ||
585 | std::size_t dim_; // const | ||
586 | coeff_t coefficients_; // const | ||
587 | std::size_t degree_; // const | ||
588 | time_t T_min_, T_max_; // const | ||
589 | /*Attributes*/ | ||
590 | |||
591 | private: | ||
592 | 25 | void assert_operator_compatible(const polynomial_t& other) const { | |
593 | 25 | if ((fabs(min() - other.min()) > ndcurves::MARGIN) || | |
594 |
6/6✓ Branch 0 taken 19 times.
✓ Branch 1 taken 6 times.
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 11 times.
✓ Branch 7 taken 14 times.
|
41 | (fabs(max() - other.max()) > ndcurves::MARGIN) || |
595 |
2/2✓ Branch 2 taken 2 times.
✓ Branch 3 taken 14 times.
|
16 | dim() != other.dim()) { |
596 |
1/2✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
|
11 | throw std::invalid_argument( |
597 | "Can't perform base operation (+ - ) on two polynomials with " | ||
598 | "different time ranges or different dimensions"); | ||
599 | } | ||
600 | 14 | } | |
601 | |||
602 | template <typename In> | ||
603 | 278 | coeff_t init_coeffs(In zeroOrderCoefficient, In highestOrderCoefficient) { | |
604 | 278 | std::size_t size = | |
605 |
1/2✓ Branch 1 taken 139 times.
✗ Branch 2 not taken.
|
278 | std::distance(zeroOrderCoefficient, highestOrderCoefficient); |
606 |
1/2✓ Branch 1 taken 139 times.
✗ Branch 2 not taken.
|
278 | coeff_t res = coeff_t(dim_, size); |
607 | 278 | int i = 0; | |
608 |
2/2✓ Branch 1 taken 480 times.
✓ Branch 2 taken 139 times.
|
1238 | for (In cit = zeroOrderCoefficient; cit != highestOrderCoefficient; |
609 | 960 | ++cit, ++i) { | |
610 |
2/4✓ Branch 2 taken 480 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 480 times.
✗ Branch 6 not taken.
|
960 | res.col(i) = *cit; |
611 | } | ||
612 | 556 | return res; | |
613 | } | ||
614 | |||
615 | public: | ||
616 | // Serialization of the class | ||
617 | friend class boost::serialization::access; | ||
618 | |||
619 | template <class Archive> | ||
620 | 152 | void serialize(Archive& ar, const unsigned int version) { | |
621 | if (version) { | ||
622 | // Do something depending on version ? | ||
623 | } | ||
624 |
1/2✓ Branch 3 taken 76 times.
✗ Branch 4 not taken.
|
152 | ar& BOOST_SERIALIZATION_BASE_OBJECT_NVP(curve_abc_t); |
625 |
1/2✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
|
152 | ar& boost::serialization::make_nvp("dim", dim_); |
626 |
1/2✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
|
152 | ar& boost::serialization::make_nvp("coefficients", coefficients_); |
627 |
1/2✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
|
152 | ar& boost::serialization::make_nvp("dim", dim_); |
628 |
1/2✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
|
152 | ar& boost::serialization::make_nvp("degree", degree_); |
629 |
1/2✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
|
152 | ar& boost::serialization::make_nvp("T_min", T_min_); |
630 |
1/2✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
|
152 | ar& boost::serialization::make_nvp("T_max", T_max_); |
631 | } | ||
632 | |||
633 | }; // class polynomial | ||
634 | |||
635 | template <typename T, typename N, bool S, typename P, typename TP> | ||
636 | 6 | polynomial<T, N, S, P, TP> operator+(const polynomial<T, N, S, P, TP>& p1, | |
637 | const polynomial<T, N, S, P, TP>& p2) { | ||
638 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | polynomial<T, N, S, P, TP> res(p1); |
639 |
3/4✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
9 | return res += p2; |
640 | 6 | } | |
641 | |||
642 | template <typename T, typename N, bool S, typename P, typename TP> | ||
643 | 2 | polynomial<T, N, S, P, TP> operator+( | |
644 | const polynomial<T, N, S, P, TP>& p1, | ||
645 | const typename polynomial<T, N, S, P, TP>::point_t& point) { | ||
646 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | polynomial<T, N, S, P, TP> res(p1); |
647 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | return res += point; |
648 | 2 | } | |
649 | |||
650 | template <typename T, typename N, bool S, typename P, typename TP> | ||
651 | 1 | polynomial<T, N, S, P, TP> operator+( | |
652 | const typename polynomial<T, N, S, P, TP>::point_t& point, | ||
653 | const polynomial<T, N, S, P, TP>& p1) { | ||
654 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | polynomial<T, N, S, P, TP> res(p1); |
655 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | return res += point; |
656 | 1 | } | |
657 | |||
658 | template <typename T, typename N, bool S, typename P, typename TP> | ||
659 | 2 | polynomial<T, N, S, P, TP> operator-( | |
660 | const polynomial<T, N, S, P, TP>& p1, | ||
661 | const typename polynomial<T, N, S, P, TP>::point_t& point) { | ||
662 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | polynomial<T, N, S, P, TP> res(p1); |
663 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | return res -= point; |
664 | 2 | } | |
665 | |||
666 | template <typename T, typename N, bool S, typename P, typename TP> | ||
667 | 1 | polynomial<T, N, S, P, TP> operator-( | |
668 | const typename polynomial<T, N, S, P, TP>::point_t& point, | ||
669 | const polynomial<T, N, S, P, TP>& p1) { | ||
670 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | polynomial<T, N, S, P, TP> res(-p1); |
671 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | return res += point; |
672 | 1 | } | |
673 | |||
674 | template <typename T, typename N, bool S, typename P, typename TP> | ||
675 | 3 | polynomial<T, N, S, P, TP> operator-(const polynomial<T, N, S, P, TP>& p1) { | |
676 |
3/6✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
|
3 | typename polynomial<T, N, S, P, TP>::coeff_t res = -p1.coeff(); |
677 |
3/6✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
|
6 | return polynomial<T, N, S, P, TP>(res, p1.min(), p1.max()); |
678 | 3 | } | |
679 | |||
680 | template <typename T, typename N, bool S, typename P, typename TP> | ||
681 | 6 | polynomial<T, N, S, P, TP> operator-(const polynomial<T, N, S, P, TP>& p1, | |
682 | const polynomial<T, N, S, P, TP>& p2) { | ||
683 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | polynomial<T, N, S, P, TP> res(p1); |
684 |
3/4✓ Branch 1 taken 3 times.
✓ Branch 2 taken 3 times.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
9 | return res -= p2; |
685 | 6 | } | |
686 | |||
687 | template <typename T, typename N, bool S, typename P, typename TP> | ||
688 | 2 | polynomial<T, N, S, P, TP> operator/(const polynomial<T, N, S, P, TP>& p1, | |
689 | const double k) { | ||
690 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | polynomial<T, N, S, P, TP> res(p1); |
691 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | return res /= k; |
692 | 2 | } | |
693 | |||
694 | template <typename T, typename N, bool S, typename P, typename TP> | ||
695 | 1 | polynomial<T, N, S, P, TP> operator*(const polynomial<T, N, S, P, TP>& p1, | |
696 | const double k) { | ||
697 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | polynomial<T, N, S, P, TP> res(p1); |
698 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | return res *= k; |
699 | 1 | } | |
700 | |||
701 | template <typename T, typename N, bool S, typename P, typename TP> | ||
702 | 1 | polynomial<T, N, S, P, TP> operator*(const double k, | |
703 | const polynomial<T, N, S, P, TP>& p1) { | ||
704 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | polynomial<T, N, S, P, TP> res(p1); |
705 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | return res *= k; |
706 | 1 | } | |
707 | |||
708 | } // namespace ndcurves | ||
709 | |||
710 | DEFINE_CLASS_TEMPLATE_VERSION( | ||
711 | SINGLE_ARG(typename Time, typename Numeric, bool Safe, typename Point, | ||
712 | typename T_Point), | ||
713 | SINGLE_ARG(ndcurves::polynomial<Time, Numeric, Safe, Point, T_Point>)) | ||
714 | |||
715 | #endif //_STRUCT_POLYNOMIAL | ||
716 |