GCC Code Coverage Report

Directory: ./
File: include/2geom/piecewise.h
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 217 370 58.6%
Functions: 48 89 53.9%
Branches: 163 519 31.4%
Line Branch Exec Source
1 /** @file
2 * @brief Piecewise function class
3 *//*
4 * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it either under the terms of the GNU Lesser General Public
8 * License version 2.1 as published by the Free Software Foundation
9 * (the "LGPL") or, at your option, under the terms of the Mozilla
10 * Public License Version 1.1 (the "MPL"). If you do not alter this
11 * notice, a recipient may use your version of this file under either
12 * the MPL or the LGPL.
13 *
14 * You should have received a copy of the LGPL along with this library
15 * in the file COPYING-LGPL-2.1; if not, output to the Free Software
16 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 * You should have received a copy of the MPL along with this library
18 * in the file COPYING-MPL-1.1
19 *
20 * The contents of this file are subject to the Mozilla Public License
21 * Version 1.1 (the "License"); you may not use this file except in
22 * compliance with the License. You may obtain a copy of the License at
23 * http://www.mozilla.org/MPL/
24 *
25 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
26 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
27 * the specific language governing rights and limitations.
28 *
29 */
30
31 #ifndef LIB2GEOM_SEEN_PIECEWISE_H
32 #define LIB2GEOM_SEEN_PIECEWISE_H
33
34 #include <vector>
35 #include <map>
36 #include <utility>
37 #include <boost/concept_check.hpp>
38 #include <2geom/concepts.h>
39 #include <2geom/math-utils.h>
40 #include <2geom/sbasis.h>
41
42
43 namespace Geom {
44 /**
45 * @brief Function defined as discrete pieces.
46 *
47 * The Piecewise class manages a sequence of elements of a type as segments and
48 * the ’cuts’ between them. These cuts are time values which separate the pieces.
49 * This function representation allows for more interesting functions, as it provides
50 * a viable output for operations such as inversion, which may require multiple
51 * SBasis to properly invert the original.
52 *
53 * As for technical details, while the actual SBasis segments begin on the first
54 * cut and end on the last, the function is defined throughout all inputs by ex-
55 * tending the first and last segments. The exact switching between segments is
56 * arbitrarily such that beginnings (t=0) have preference over endings (t=1). This
57 * only matters if it is discontinuous at the location.
58 * \f[
59 * f(t) \rightarrow \left\{
60 * \begin{array}{cc}
61 * s_1,& t <= c_2 \\
62 * s_2,& c_2 <= t <= c_3\\
63 * \ldots \\
64 * s_n,& c_n <= t
65 * \end{array}\right.
66 * \f]
67 *
68 * @ingroup Fragments
69 */
70 template <typename T>
71 class Piecewise {
72 BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept);
73
74 public:
75 std::vector<double> cuts;
76 std::vector<T> segs;
77 //segs[i] stretches from cuts[i] to cuts[i+1].
78
79 296433 Piecewise() {}
80
81 30900 explicit Piecewise(const T &s) {
82
1/2
✓ Branch 1 taken 30900 times.
✗ Branch 2 not taken.
 30900 push_cut(0.);
83
1/2
✓ Branch 1 taken 30900 times.
✗ Branch 2 not taken.
 30900 push_seg(s);
84
1/2
✓ Branch 1 taken 30900 times.
✗ Branch 2 not taken.
 30900 push_cut(1.);
85 30900 }
86
87 unsigned input_dim(){return 1;}
88
89 typedef typename T::output_type output_type;
90
91 ✗ explicit Piecewise(const output_type & v) {
92 ✗ push_cut(0.);
93 ✗ push_seg(T(v));
94 ✗ push_cut(1.);
95 ✗ }
96
97 118800 inline void reserve(unsigned i) { segs.reserve(i); cuts.reserve(i + 1); }
98
99 528840 inline T const& operator[](unsigned i) const { return segs[i]; }
100 325800 inline T& operator[](unsigned i) { return segs[i]; }
101 35100 inline output_type operator()(double t) const { return valueAt(t); }
102 35100 inline output_type valueAt(double t) const {
103 35100 unsigned n = segN(t);
104 35100 return segs[n](segT(t, n));
105 }
106 ✗ inline output_type firstValue() const {
107 ✗ return valueAt(cuts.front());
108 }
109 ✗ inline output_type lastValue() const {
110 ✗ return valueAt(cuts.back());
111 }
112
113 /**
114 * The size of the returned vector equals n_derivs+1.
115 */
116 std::vector<output_type> valueAndDerivatives(double t, unsigned n_derivs) const {
117 unsigned n = segN(t);
118 std::vector<output_type> ret, val = segs[n].valueAndDerivatives(segT(t, n), n_derivs);
119 double mult = 1;
120 for(unsigned i = 0; i < val.size(); i++) {
121 ret.push_back(val[i] * mult);
122 mult /= cuts[n+1] - cuts[n];
123 }
124 return ret;
125 }
126
127 //TODO: maybe it is not a good idea to have this?
128 Piecewise<T> operator()(SBasis f);
129 Piecewise<T> operator()(Piecewise<SBasis>f);
130
131 2274108 inline unsigned size() const { return segs.size(); }
132 502800 inline bool empty() const { return segs.empty(); }
133 inline void clear() {
134 segs.clear();
135 cuts.clear();
136 }
137
138 /**Convenience/implementation hiding function to add segment/cut pairs.
139 * Asserts that basic size and order invariants are correct
140 */
141 35100 inline void push(const T &s, double to) {
142
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 35100 times.
 35100 assert(cuts.size() - segs.size() == 1);
143 35100 push_seg(s);
144 35100 push_cut(to);
145 35100 }
146 145098 inline void push(T &&s, double to) {
147
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 72582 times.
 145098 assert(cuts.size() - segs.size() == 1);
148 145098 push_seg(s);
149 145098 push_cut(to);
150 145098 }
151 //Convenience/implementation hiding function to add cuts.
152 869931 inline void push_cut(double c) {
153
3/6
✓ Branch 1 taken 325182 times.
✓ Branch 2 taken 109827 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 325182 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 869931 ASSERT_INVARIANTS(cuts.empty() || c > cuts.back());
154 869931 cuts.push_back(c);
155 869931 }
156 //Convenience/implementation hiding function to add segments.
157 578298 inline void push_seg(const T &s) { segs.push_back(s); }
158 115200 inline void push_seg(T &&s) { segs.emplace_back(s); }
159
160 /**Returns the segment index which corresponds to a 'global' piecewise time.
161 * Also takes optional low/high parameters to expedite the search for the segment.
162 */
163 70200 inline unsigned segN(double t, int low = 0, int high = -1) const {
164
1/2
✓ Branch 0 taken 35100 times.
✗ Branch 1 not taken.
 70200 high = (high == -1) ? size() : high;
165
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 35100 times.
 70200 if(t < cuts[0]) return 0;
166
2/2
✓ Branch 2 taken 19200 times.
✓ Branch 3 taken 15900 times.
 70200 if(t >= cuts[size()]) return size() - 1;
167
2/2
✓ Branch 0 taken 18900 times.
✓ Branch 1 taken 2700 times.
 43200 while(low < high) {
168 37800 int mid = (high + low) / 2; //Lets not plan on having huge (> INT_MAX / 2) cut sequences
169 37800 double mv = cuts[mid];
170
2/2
✓ Branch 0 taken 2700 times.
✓ Branch 1 taken 16200 times.
 37800 if(mv < t) {
171
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2700 times.
 5400 if(t < cuts[mid + 1]) return mid; else low = mid + 1;
172
2/2
✓ Branch 0 taken 3000 times.
✓ Branch 1 taken 13200 times.
 32400 } else if(t < mv) {
173
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3000 times.
 6000 if(cuts[mid - 1] < t) return mid - 1; else high = mid - 1;
174 } else {
175 26400 return mid;
176 }
177 }
178 5400 return low;
179 }
180
181 /**Returns the time within a segment, given the 'global' piecewise time.
182 * Also takes an optional index parameter which may be used for efficiency or to find the time on a
183 * segment outside its range. If it is left to its default, -1, it will call segN to find the index.
184 */
185 82200 inline double segT(double t, int i = -1) const {
186
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 41100 times.
 82200 if(i == -1) i = segN(t);
187
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 41100 times.
 82200 assert(i >= 0);
188 82200 return (t - cuts[i]) / (cuts[i+1] - cuts[i]);
189 }
190
191 ✗ inline double mapToDomain(double t, unsigned i) const {
192 ✗ return (1-t)*cuts[i] + t*cuts[i+1]; //same as: t * (cuts[i+1] - cuts[i]) + cuts[i]
193 }
194
195 //Offsets the piecewise domain
196 ✗ inline void offsetDomain(double o) {
197 ✗ assert(std::isfinite(o));
198 ✗ if(o != 0)
199 ✗ for(unsigned i = 0; i <= size(); i++)
200 ✗ cuts[i] += o;
201 ✗ }
202
203 //Scales the domain of the function by a value. 0 will result in an empty Piecewise.
204 ✗ inline void scaleDomain(double s) {
205 ✗ assert(s > 0);
206 ✗ if(s == 0) {
207 ✗ cuts.clear(); segs.clear();
208 ✗ return;
209 }
210 ✗ for(unsigned i = 0; i <= size(); i++)
211 ✗ cuts[i] *= s;
212 }
213
214 //Retrieves the domain in interval form
215
1/2
✓ Branch 4 taken 9600 times.
✗ Branch 5 not taken.
 9600 inline Interval domain() const { return Interval(cuts.front(), cuts.back()); }
216
217 //Transforms the domain into another interval
218 47400 inline void setDomain(Interval dom) {
219
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 47400 times.
 47400 if(empty()) return;
220 /* dom can not be empty
221 if(dom.empty()) {
222 cuts.clear(); segs.clear();
223 return;
224 }*/
225 47400 double cf = cuts.front();
226 47400 double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);
227
2/2
✓ Branch 1 taken 116100 times.
✓ Branch 2 taken 47400 times.
 163500 for(unsigned i = 0; i <= size(); i++)
228 116100 cuts[i] = (cuts[i] - cf) * s + o;
229 //fix floating point precision errors.
230 47400 cuts[0] = dom.min();
231 47400 cuts[size()] = dom.max();
232 }
233
234 //Concatenates this Piecewise function with another, offsetting time of the other to match the end.
235 45300 inline void concat(const Piecewise<T> &other) {
236
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 45300 times.
 45300 if(other.empty()) return;
237
238
2/2
✓ Branch 1 taken 10500 times.
✓ Branch 2 taken 34800 times.
 45300 if(empty()) {
239 10500 cuts = other.cuts; segs = other.segs;
240 10500 return;
241 }
242
243
1/2
✓ Branch 7 taken 34800 times.
✗ Branch 8 not taken.
 34800 segs.insert(segs.end(), other.segs.begin(), other.segs.end());
244 34800 double t = cuts.back() - other.cuts.front();
245 34800 cuts.reserve(cuts.size() + other.size());
246
2/2
✓ Branch 1 taken 54000 times.
✓ Branch 2 taken 34800 times.
 88800 for(unsigned i = 0; i < other.size(); i++)
247 54000 push_cut(other.cuts[i + 1] + t);
248 }
249
250 //Like concat, but ensures continuity.
251 ✗ inline void continuousConcat(const Piecewise<T> &other) {
252 ✗ boost::function_requires<AddableConcept<typename T::output_type> >();
253 ✗ if(other.empty()) return;
254
255 ✗ if(empty()) { segs = other.segs; cuts = other.cuts; return; }
256
257 ✗ typename T::output_type y = segs.back().at1() - other.segs.front().at0();
258 ✗ double t = cuts.back() - other.cuts.front();
259 ✗ reserve(size() + other.size());
260 ✗ for(unsigned i = 0; i < other.size(); i++)
261 ✗ push(other[i] + y, other.cuts[i + 1] + t);
262 }
263
264 //returns true if the Piecewise<T> meets some basic invariants.
265 136800 inline bool invariants() const {
266 // segs between cuts
267
2/8
✗ Branch 2 not taken.
✓ Branch 3 taken 68400 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 68400 times.
 136800 if(!(segs.size() + 1 == cuts.size() || (segs.empty() && cuts.empty())))
268 ✗ return false;
269 // cuts in order
270
2/2
✓ Branch 1 taken 142200 times.
✓ Branch 2 taken 68400 times.
 421200 for(unsigned i = 0; i < segs.size(); i++)
271
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 142200 times.
 284400 if(cuts[i] >= cuts[i+1])
272 ✗ return false;
273 136800 return true;
274 }
275
276 };
277
278 /**
279 * ...
280 * \return ...
281 * \relates Piecewise
282 */
283 template<typename T>
284 ✗ inline typename FragmentConcept<T>::BoundsType bounds_fast(const Piecewise<T> &f) {
285 ✗ boost::function_requires<FragmentConcept<T> >();
286
287 ✗ if(f.empty()) return typename FragmentConcept<T>::BoundsType();
288 ✗ typename FragmentConcept<T>::BoundsType ret(bounds_fast(f[0]));
289 ✗ for(unsigned i = 1; i < f.size(); i++)
290 ✗ ret.unionWith(bounds_fast(f[i]));
291 ✗ return ret;
292 }
293
294 /**
295 * ...
296 * \return ...
297 * \relates Piecewise
298 */
299 template<typename T>
300 300 inline typename FragmentConcept<T>::BoundsType bounds_exact(const Piecewise<T> &f) {
301 300 boost::function_requires<FragmentConcept<T> >();
302
303
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 300 times.
 300 if(f.empty()) return typename FragmentConcept<T>::BoundsType();
304
1/2
✓ Branch 3 taken 300 times.
✗ Branch 4 not taken.
 300 typename FragmentConcept<T>::BoundsType ret(bounds_exact(f[0]));
305
2/2
✓ Branch 1 taken 11700 times.
✓ Branch 2 taken 300 times.
 12000 for(unsigned i = 1; i < f.size(); i++)
306
2/4
✓ Branch 3 taken 11700 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 11700 times.
✗ Branch 7 not taken.
 11700 ret.unionWith(bounds_exact(f[i]));
307 300 return ret;
308 }
309
310 /**
311 * ...
312 * \return ...
313 * \relates Piecewise
314 */
315 template<typename T>
316 ✗ inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const OptInterval &_m) {
317 ✗ boost::function_requires<FragmentConcept<T> >();
318
319 ✗ if(f.empty() || !_m) return typename FragmentConcept<T>::BoundsType();
320 ✗ Interval const &m = *_m;
321 ✗ if(m.isSingular()) return typename FragmentConcept<T>::BoundsType(f(m.min()));
322
323 ✗ unsigned fi = f.segN(m.min()), ti = f.segN(m.max());
324 ✗ double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);
325
326 ✗ if(fi == ti) return bounds_local(f[fi], Interval(ft, tt));
327
328 ✗ typename FragmentConcept<T>::BoundsType ret(bounds_local(f[fi], Interval(ft, 1.)));
329 ✗ for(unsigned i = fi + 1; i < ti; i++)
330 ✗ ret.unionWith(bounds_exact(f[i]));
331 ✗ if(tt != 0.) ret.unionWith(bounds_local(f[ti], Interval(0., tt)));
332
333 ✗ return ret;
334 }
335
336 /**
337 * Returns a portion of a piece of a Piecewise<T>, given the piece's index and a to/from time.
338 * \relates Piecewise
339 */
340 template<typename T>
341 18000 T elem_portion(const Piecewise<T> &a, unsigned i, double from, double to) {
342
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9000 times.
 18000 assert(i < a.size());
343 18000 double rwidth = 1 / (a.cuts[i+1] - a.cuts[i]);
344 18000 return portion( a[i], (from - a.cuts[i]) * rwidth, (to - a.cuts[i]) * rwidth );
345 }
346
347 /**Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c);
348 * Further subdivides the Piecewise<T> such that there is a cut at every value in c.
349 * Precondition: c sorted lower to higher.
350 *
351 * //Given Piecewise<T> a and b:
352 * Piecewise<T> ac = a.partition(b.cuts);
353 * Piecewise<T> bc = b.partition(a.cuts);
354 * //ac.cuts should be equivalent to bc.cuts
355 *
356 * \relates Piecewise
357 */
358 template<typename T>
359 136800 Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c) {
360
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 68400 times.
 136800 assert(pw.invariants());
361
3/4
✓ Branch 1 taken 9600 times.
✓ Branch 2 taken 58800 times.
✓ Branch 4 taken 9600 times.
✗ Branch 5 not taken.
 136800 if(c.empty()) return Piecewise<T>(pw);
362
363 117600 Piecewise<T> ret = Piecewise<T>();
364
1/2
✓ Branch 3 taken 58800 times.
✗ Branch 4 not taken.
 117600 ret.reserve(c.size() + pw.cuts.size() - 1);
365
366
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 58800 times.
 117600 if(pw.empty()) {
367 ✗ ret.cuts = c;
368 ✗ for(unsigned i = 0; i < c.size() - 1; i++)
369 ✗ ret.push_seg(T());
370 ✗ return ret;
371 }
372
373 117600 unsigned si = 0, ci = 0; //Segment index, Cut index
374
375 //if the cuts have something earlier than the Piecewise<T>, add portions of the first segment
376
3/6
✓ Branch 1 taken 58800 times.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 58800 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 58800 times.
 117600 while(ci < c.size() && c[ci] < pw.cuts.front()) {
377 ✗ bool isLast = (ci == c.size()-1 || c[ci + 1] >= pw.cuts.front());
378 ✗ ret.push_cut(c[ci]);
379 ✗ ret.push_seg( elem_portion(pw, 0, c[ci], isLast ? pw.cuts.front() : c[ci + 1]) );
380 ✗ ci++;
381 }
382
383
1/2
✓ Branch 2 taken 58800 times.
✗ Branch 3 not taken.
 117600 ret.push_cut(pw.cuts[0]);
384 117600 double prev = pw.cuts[0]; //previous cut
385 //Loop which handles cuts within the Piecewise<T> domain
386 //Should have the cuts = segs + 1 invariant
387
5/6
✓ Branch 1 taken 265200 times.
✓ Branch 2 taken 58800 times.
✓ Branch 4 taken 265200 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 265200 times.
✓ Branch 7 taken 58800 times.
 648000 while(si < pw.size() && ci <= c.size()) {
388
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 265200 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 265200 times.
 530400 if(ci == c.size() && prev <= pw.cuts[si]) { //cuts exhausted, straight copy the rest
389 ✗ ret.segs.insert(ret.segs.end(), pw.segs.begin() + si, pw.segs.end());
390 ✗ ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + si + 1, pw.cuts.end());
391 ✗ return ret;
392
5/6
✓ Branch 1 taken 265200 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 132600 times.
✓ Branch 6 taken 132600 times.
✓ Branch 7 taken 132600 times.
✓ Branch 8 taken 132600 times.
 530400 }else if(ci == c.size() || c[ci] >= pw.cuts[si + 1]) { //no more cuts within this segment, finalize
393
2/2
✓ Branch 1 taken 6000 times.
✓ Branch 2 taken 126600 times.
 265200 if(prev > pw.cuts[si]) { //segment already has cuts, so portion is required
394
3/9
✓ Branch 2 taken 6000 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 6000 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6000 times.
✗ Branch 10 not taken.
 12000 ret.push_seg(portion(pw[si], pw.segT(prev, si), 1.0));
395 } else { //plain copy is fine
396
1/2
✓ Branch 2 taken 126600 times.
✗ Branch 3 not taken.
 253200 ret.push_seg(pw[si]);
397 }
398
1/2
✓ Branch 2 taken 132600 times.
✗ Branch 3 not taken.
 265200 ret.push_cut(pw.cuts[si + 1]);
399 265200 prev = pw.cuts[si + 1];
400 265200 si++;
401
2/2
✓ Branch 2 taken 123600 times.
✓ Branch 3 taken 9000 times.
 265200 } else if(c[ci] == pw.cuts[si]){ //coincident
402 //Already finalized the seg with the code immediately above
403 247200 ci++;
404 } else { //plain old subdivision
405
2/4
✓ Branch 4 taken 9000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9000 times.
✗ Branch 8 not taken.
 18000 ret.push(elem_portion(pw, si, prev, c[ci]), c[ci]);
406 18000 prev = c[ci];
407 18000 ci++;
408 }
409 }
410
411 //input cuts extend further than this Piecewise<T>, extend the last segment.
412
2/2
✓ Branch 1 taken 58800 times.
✓ Branch 2 taken 58800 times.
 235200 while(ci < c.size()) {
413
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 58800 times.
 117600 if(c[ci] > prev) {
414 ✗ ret.push(elem_portion(pw, pw.size() - 1, prev, c[ci]), c[ci]);
415 ✗ prev = c[ci];
416 }
417 117600 ci++;
418 }
419 117600 return ret;
420 117600 }
421
422 /**
423 * Returns a Piecewise<T> with a defined domain of [min(from, to), max(from, to)].
424 * \relates Piecewise
425 */
426 template<typename T>
427 ✗ Piecewise<T> portion(const Piecewise<T> &pw, double from, double to) {
428 ✗ if(pw.empty() || from == to) return Piecewise<T>();
429
430 ✗ Piecewise<T> ret;
431
432 ✗ double temp = from;
433 ✗ from = std::min(from, to);
434 ✗ to = std::max(temp, to);
435
436 ✗ unsigned i = pw.segN(from);
437 ✗ ret.push_cut(from);
438 ✗ if(i == pw.size() - 1 || to <= pw.cuts[i + 1]) { //to/from inhabit the same segment
439 ✗ ret.push(elem_portion(pw, i, from, to), to);
440 ✗ return ret;
441 }
442 ✗ ret.push_seg(portion( pw[i], pw.segT(from, i), 1.0 ));
443 ✗ i++;
444 ✗ unsigned fi = pw.segN(to, i);
445 ✗ ret.reserve(fi - i + 1);
446 ✗ if (to == pw.cuts[fi]) fi-=1;
447
448 ✗ ret.segs.insert(ret.segs.end(), pw.segs.begin() + i, pw.segs.begin() + fi); //copy segs
449 ✗ ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + i, pw.cuts.begin() + fi + 1); //and their cuts
450
451 ✗ ret.push_seg( portion(pw[fi], 0.0, pw.segT(to, fi)));
452 ✗ if(to != ret.cuts.back()) ret.push_cut(to);
453 ✗ ret.invariants();
454 ✗ return ret;
455 ✗ }
456
457 //TODO: seems like these should be mutating
458 /**
459 * ...
460 * \return ...
461 * \relates Piecewise
462 */
463 template<typename T>
464 600 Piecewise<T> remove_short_cuts(Piecewise<T> const &f, double tol) {
465
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 600 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 600 if(f.empty()) return f;
466 600 Piecewise<T> ret;
467
1/2
✓ Branch 2 taken 600 times.
✗ Branch 3 not taken.
 600 ret.reserve(f.size());
468
1/2
✓ Branch 2 taken 600 times.
✗ Branch 3 not taken.
 600 ret.push_cut(f.cuts[0]);
469
2/2
✓ Branch 1 taken 23400 times.
✓ Branch 2 taken 600 times.
 24000 for(unsigned i=0; i<f.size(); i++){
470
2/6
✗ Branch 2 not taken.
✓ Branch 3 taken 23400 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 23400 times.
✗ Branch 8 not taken.
 23400 if (f.cuts[i+1]-f.cuts[i] >= tol || i==f.size()-1) {
471
1/2
✓ Branch 3 taken 23400 times.
✗ Branch 4 not taken.
 23400 ret.push(f[i], f.cuts[i+1]);
472 }
473 }
474 600 return ret;
475 600 }
476
477 //TODO: seems like these should be mutating
478 /**
479 * ...
480 * \return ...
481 * \relates Piecewise
482 */
483 template<typename T>
484 Piecewise<T> remove_short_cuts_extending(Piecewise<T> const &f, double tol) {
485 if(f.empty()) return f;
486 Piecewise<T> ret;
487 ret.reserve(f.size());
488 ret.push_cut(f.cuts[0]);
489 double last = f.cuts[0]; // last cut included
490 for(unsigned i=0; i<f.size(); i++){
491 if (f.cuts[i+1]-f.cuts[i] >= tol) {
492 ret.push(elem_portion(f, i, last, f.cuts[i+1]), f.cuts[i+1]);
493 last = f.cuts[i+1];
494 }
495 }
496 return ret;
497 }
498
499 /**
500 * ...
501 * \return ...
502 * \relates Piecewise
503 */
504 template<typename T>
505 std::vector<double> roots(const Piecewise<T> &pw) {
506 std::vector<double> ret;
507 for(unsigned i = 0; i < pw.size(); i++) {
508 std::vector<double> sr = roots(pw[i]);
509 for (double & j : sr) ret.push_back(j * (pw.cuts[i + 1] - pw.cuts[i]) + pw.cuts[i]);
510
511 }
512 return ret;
513 }
514
515 //IMPL: OffsetableConcept
516 /**
517 * ...
518 * \return \f$ a + b = \f$
519 * \relates Piecewise
520 */
521 template<typename T>
522 ✗ Piecewise<T> operator+(Piecewise<T> const &a, typename T::output_type b) {
523 ✗ boost::function_requires<OffsetableConcept<T> >();
524 //TODO:empty
525 ✗ Piecewise<T> ret;
526 ✗ ret.segs.reserve(a.size());
527 ✗ ret.cuts = a.cuts;
528 ✗ for(unsigned i = 0; i < a.size();i++)
529 ✗ ret.push_seg(a[i] + b);
530 ✗ return ret;
531 ✗ }
532 template<typename T>
533 Piecewise<T> operator-(Piecewise<T> const &a, typename T::output_type b) {
534 boost::function_requires<OffsetableConcept<T> >();
535 //TODO: empty
536 Piecewise<T> ret;
537 ret.segs.reserve(a.size());
538 ret.cuts = a.cuts;
539 for(unsigned i = 0; i < a.size();i++)
540 ret.push_seg(a[i] - b);
541 return ret;
542 }
543 template<typename T>
544 ✗ Piecewise<T>& operator+=(Piecewise<T>& a, typename T::output_type b) {
545 ✗ boost::function_requires<OffsetableConcept<T> >();
546
547 ✗ if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
548
549 ✗ for(unsigned i = 0; i < a.size();i++)
550 ✗ a[i] += b;
551 ✗ return a;
552 }
553 template<typename T>
554 10200 Piecewise<T>& operator-=(Piecewise<T>& a, typename T::output_type b) {
555 10200 boost::function_requires<OffsetableConcept<T> >();
556
557
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10200 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 10200 if(a.empty()) { a.push_cut(0.); a.push(T(-b), 1.); return a; }
558
559
2/2
✓ Branch 1 taken 35700 times.
✓ Branch 2 taken 10200 times.
 45900 for(unsigned i = 0;i < a.size();i++)
560 35700 a[i] -= b;
561 10200 return a;
562 }
563
564 //IMPL: ScalableConcept
565 /**
566 * ...
567 * \return \f$ -a = \f$
568 * \relates Piecewise
569 */
570 template<typename T>
571 ✗ Piecewise<T> operator-(Piecewise<T> const &a) {
572 ✗ boost::function_requires<ScalableConcept<T> >();
573
574 ✗ Piecewise<T> ret;
575 ✗ ret.segs.reserve(a.size());
576 ✗ ret.cuts = a.cuts;
577 ✗ for(unsigned i = 0; i < a.size();i++)
578 ✗ ret.push_seg(- a[i]);
579 ✗ return ret;
580 ✗ }
581 /**
582 * ...
583 * \return \f$ a * b = \f$
584 * \relates Piecewise
585 */
586 template<typename T>
587 ✗ Piecewise<T> operator*(Piecewise<T> const &a, double b) {
588 ✗ boost::function_requires<ScalableConcept<T> >();
589
590 ✗ if(a.empty()) return Piecewise<T>();
591
592 ✗ Piecewise<T> ret;
593 ✗ ret.segs.reserve(a.size());
594 ✗ ret.cuts = a.cuts;
595 ✗ for(unsigned i = 0; i < a.size();i++)
596 ✗ ret.push_seg(a[i] * b);
597 ✗ return ret;
598 ✗ }
599 /**
600 * ...
601 * \return \f$ a * b = \f$
602 * \relates Piecewise
603 */
604 template<typename T>
605 Piecewise<T> operator*(Piecewise<T> const &a, T b) {
606 boost::function_requires<ScalableConcept<T> >();
607
608 if(a.empty()) return Piecewise<T>();
609
610 Piecewise<T> ret;
611 ret.segs.reserve(a.size());
612 ret.cuts = a.cuts;
613 for(unsigned i = 0; i < a.size();i++)
614 ret.push_seg(a[i] * b);
615 return ret;
616 }
617 /**
618 * ...
619 * \return \f$ a / b = \f$
620 * \relates Piecewise
621 */
622 template<typename T>
623 Piecewise<T> operator/(Piecewise<T> const &a, double b) {
624 boost::function_requires<ScalableConcept<T> >();
625
626 //FIXME: b == 0?
627 if(a.empty()) return Piecewise<T>();
628
629 Piecewise<T> ret;
630 ret.segs.reserve(a.size());
631 ret.cuts = a.cuts;
632 for(unsigned i = 0; i < a.size();i++)
633 ret.push_seg(a[i] / b);
634 return ret;
635 }
636 template<typename T>
637 600 Piecewise<T>& operator*=(Piecewise<T>& a, double b) {
638 600 boost::function_requires<ScalableConcept<T> >();
639
640
2/2
✓ Branch 1 taken 24000 times.
✓ Branch 2 taken 600 times.
 24600 for(unsigned i = 0; i < a.size();i++)
641
0/2
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 24000 a[i] *= b;
642 600 return a;
643 }
644 template<typename T>
645 ✗ Piecewise<T>& operator/=(Piecewise<T>& a, double b) {
646 ✗ boost::function_requires<ScalableConcept<T> >();
647
648 //FIXME: b == 0?
649
650 ✗ for(unsigned i = 0; i < a.size();i++)
651 ✗ a[i] /= b;
652 ✗ return a;
653 }
654
655 //IMPL: AddableConcept
656 /**
657 * ...
658 * \return \f$ a + b = \f$
659 * \relates Piecewise
660 */
661 template<typename T>
662 300 Piecewise<T> operator+(Piecewise<T> const &a, Piecewise<T> const &b) {
663 300 boost::function_requires<AddableConcept<T> >();
664
665
2/4
✓ Branch 2 taken 300 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 300 times.
✗ Branch 7 not taken.
 300 Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
666 300 Piecewise<T> ret;
667
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 300 times.
 300 assert(pa.size() == pb.size());
668
1/2
✓ Branch 2 taken 300 times.
✗ Branch 3 not taken.
 300 ret.segs.reserve(pa.size());
669
1/2
✓ Branch 1 taken 300 times.
✗ Branch 2 not taken.
 300 ret.cuts = pa.cuts;
670
2/2
✓ Branch 1 taken 21000 times.
✓ Branch 2 taken 300 times.
 21300 for (unsigned i = 0; i < pa.size(); i++)
671
2/4
✓ Branch 4 taken 21000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 21000 times.
✗ Branch 8 not taken.
 21000 ret.push_seg(pa[i] + pb[i]);
672 600 return ret;
673 300 }
674 /**
675 * ...
676 * \return \f$ a - b = \f$
677 * \relates Piecewise
678 */
679 template<typename T>
680 9600 Piecewise<T> operator-(Piecewise<T> const &a, Piecewise<T> const &b) {
681 9600 boost::function_requires<AddableConcept<T> >();
682
683
2/4
✓ Branch 2 taken 9600 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 9600 times.
✗ Branch 7 not taken.
 9600 Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
684 9600 Piecewise<T> ret = Piecewise<T>();
685
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 9600 times.
 9600 assert(pa.size() == pb.size());
686
1/2
✓ Branch 2 taken 9600 times.
✗ Branch 3 not taken.
 9600 ret.segs.reserve(pa.size());
687
1/2
✓ Branch 1 taken 9600 times.
✗ Branch 2 not taken.
 9600 ret.cuts = pa.cuts;
688
2/2
✓ Branch 1 taken 9600 times.
✓ Branch 2 taken 9600 times.
 19200 for (unsigned i = 0; i < pa.size(); i++)
689
2/4
✓ Branch 4 taken 9600 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9600 times.
✗ Branch 8 not taken.
 9600 ret.push_seg(pa[i] - pb[i]);
690 19200 return ret;
691 9600 }
692 template<typename T>
693 inline Piecewise<T>& operator+=(Piecewise<T> &a, Piecewise<T> const &b) {
694 a = a+b;
695 return a;
696 }
697 template<typename T>
698 inline Piecewise<T>& operator-=(Piecewise<T> &a, Piecewise<T> const &b) {
699 a = a-b;
700 return a;
701 }
702
703 /**
704 * ...
705 * \return \f$ a \cdot b = \f$
706 * \relates Piecewise
707 */
708 template<typename T1,typename T2>
709 300 Piecewise<T2> operator*(Piecewise<T1> const &a, Piecewise<T2> const &b) {
710 //function_requires<MultiplicableConcept<T1> >();
711 //function_requires<MultiplicableConcept<T2> >();
712
713
1/2
✓ Branch 2 taken 300 times.
✗ Branch 3 not taken.
 300 Piecewise<T1> pa = partition(a, b.cuts);
714
1/2
✓ Branch 2 taken 300 times.
✗ Branch 3 not taken.
 300 Piecewise<T2> pb = partition(b, a.cuts);
715 300 Piecewise<T2> ret = Piecewise<T2>();
716
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 300 times.
 300 assert(pa.size() == pb.size());
717
1/2
✓ Branch 2 taken 300 times.
✗ Branch 3 not taken.
 300 ret.segs.reserve(pa.size());
718
1/2
✓ Branch 1 taken 300 times.
✗ Branch 2 not taken.
 300 ret.cuts = pa.cuts;
719
2/2
✓ Branch 1 taken 21000 times.
✓ Branch 2 taken 300 times.
 21300 for (unsigned i = 0; i < pa.size(); i++)
720
2/4
✓ Branch 4 taken 21000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 21000 times.
✗ Branch 8 not taken.
 21000 ret.push_seg(pa[i] * pb[i]);
721 600 return ret;
722 300 }
723
724 /**
725 * ...
726 * \return \f$ a \cdot b \f$
727 * \relates Piecewise
728 */
729 template<typename T>
730 ✗ inline Piecewise<T>& operator*=(Piecewise<T> &a, Piecewise<T> const &b) {
731 ✗ a = a * b;
732 ✗ return a;
733 }
734
735 Piecewise<SBasis> divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, unsigned k);
736 //TODO: replace divide(a,b,k) by divide(a,b,tol,k)?
737 //TODO: atm, relative error is ~(tol/a)%. Find a way to make it independent of a.
738 //Nota: the result is 'truncated' where b is smaller than 'zero': ~ a/max(b,zero).
739 Piecewise<SBasis>
740 divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
741 Piecewise<SBasis>
742 divide(SBasis const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
743 Piecewise<SBasis>
744 divide(Piecewise<SBasis> const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
745 Piecewise<SBasis>
746 divide(SBasis const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
747
748 //Composition: functions called compose_* are pieces of compose that are factored out in pw.cpp.
749 std::map<double,unsigned> compose_pullback(std::vector<double> const &cuts, SBasis const &g);
750 int compose_findSegIdx(std::map<double,unsigned>::iterator const &cut,
751 std::map<double,unsigned>::iterator const &next,
752 std::vector<double> const &levels,
753 SBasis const &g);
754
755 /**
756 * ...
757 * \return ...
758 * \relates Piecewise
759 */
760 template<typename T>
761 24000 Piecewise<T> compose(Piecewise<T> const &f, SBasis const &g){
762 /// \todo add concept check
763 24000 Piecewise<T> result;
764
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24000 times.
 24000 if (f.empty()) return result;
765
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 24000 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 24000 if (g.isZero()) return Piecewise<T>(f(0));
766
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24000 times.
 24000 if (f.size()==1){
767 ✗ double t0 = f.cuts[0], width = f.cuts[1] - t0;
768 ✗ return (Piecewise<T>) compose(f.segs[0],compose(Linear(-t0 / width, (1-t0) / width), g));
769 }
770
771 //first check bounds...
772
1/2
✓ Branch 3 taken 24000 times.
✗ Branch 4 not taken.
 24000 Interval bs = *bounds_fast(g);
773
3/6
✓ Branch 2 taken 24000 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 24000 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 24000 times.
 24000 if (f.cuts.front() > bs.max() || bs.min() > f.cuts.back()){
774 ✗ int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;
775 ✗ double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;
776 ✗ return (Piecewise<T>) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g));
777 }
778
779 24000 std::vector<double> levels;//we can forget first and last cuts...
780
1/2
✓ Branch 11 taken 24000 times.
✗ Branch 12 not taken.
 24000 levels.insert(levels.begin(),f.cuts.begin()+1,f.cuts.end()-1);
781 //TODO: use a std::vector<pairs<double,unsigned> > instead of a map<double,unsigned>.
782
1/2
✓ Branch 2 taken 24000 times.
✗ Branch 3 not taken.
 24000 std::map<double,unsigned> cuts_pb = compose_pullback(levels,g);
783
784 //-- Compose each piece of g with the relevant seg of f.
785
1/2
✓ Branch 2 taken 24000 times.
✗ Branch 3 not taken.
 24000 result.cuts.push_back(0.);
786 24000 std::map<double,unsigned>::iterator cut=cuts_pb.begin();
787 24000 std::map<double,unsigned>::iterator next=cut; next++;
788
2/2
✓ Branch 4 taken 42000 times.
✓ Branch 5 taken 24000 times.
 66000 while(next!=cuts_pb.end()){
789 //assert(std::abs(int((*cut).second-(*next).second))<1);
790 //TODO: find a way to recover from this error? the root finder missed some root;
791 // the levels/variations of f might be too close/fast...
792
1/2
✓ Branch 1 taken 42000 times.
✗ Branch 2 not taken.
 42000 int idx = compose_findSegIdx(cut,next,levels,g);
793 42000 double t0=(*cut).first;
794 42000 double t1=(*next).first;
795
796
1/2
✓ Branch 1 taken 42000 times.
✗ Branch 2 not taken.
 42000 if (!are_near(t0,t1,EPSILON*EPSILON)) { // prevent adding cuts that are extremely close together and that may cause trouble with rounding e.g. when reversing the path
797
2/4
✓ Branch 5 taken 42000 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 42000 times.
✗ Branch 9 not taken.
 42000 SBasis sub_g=compose(g, Linear(t0,t1));
798
3/6
✓ Branch 7 taken 42000 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 42000 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 42000 times.
✗ Branch 14 not taken.
 84000 sub_g=compose(Linear(-f.cuts[idx]/(f.cuts[idx+1]-f.cuts[idx]),
799 42000 (1-f.cuts[idx])/(f.cuts[idx+1]-f.cuts[idx])),sub_g);
800
2/4
✓ Branch 3 taken 42000 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 42000 times.
✗ Branch 7 not taken.
 42000 result.push(compose(f[idx],sub_g),t1);
801 42000 }
802
803 42000 cut++;
804 42000 next++;
805 }
806 24000 return(result);
807 24000 }
808
809 /**
810 * ...
811 * \return ...
812 * \relates Piecewise
813 */
814 template<typename T>
815 600 Piecewise<T> compose(Piecewise<T> const &f, Piecewise<SBasis> const &g){
816 /// \todo add concept check
817 600 Piecewise<T> result;
818
2/2
✓ Branch 2 taken 24000 times.
✓ Branch 3 taken 600 times.
 48600 for(unsigned i = 0; i < g.segs.size(); i++){
819
1/2
✓ Branch 3 taken 24000 times.
✗ Branch 4 not taken.
 24000 Piecewise<T> fgi=compose(f, g.segs[i]);
820
1/4
✓ Branch 4 taken 24000 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 24000 fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1]));
821
1/2
✓ Branch 1 taken 24000 times.
✗ Branch 2 not taken.
 24000 result.concat(fgi);
822 }
823 600 return result;
824 ✗ }
825
826 /*
827 Piecewise<D2<SBasis> > compose(D2<SBasis2d> const &sb2d, Piecewise<D2<SBasis> > const &pwd2sb){
828 /// \todo add concept check
829 Piecewise<D2<SBasis> > result;
830 result.push_cut(0.);
831 for(unsigned i = 0; i < pwd2sb.size(); i++){
832 result.push(compose_each(sb2d,pwd2sb[i]),i+1);
833 }
834 return result;
835 }*/
836
837 /** Compose an SBasis with the inverse of another.
838 * WARNING: It's up to the user to check that the second SBasis is indeed
839 * invertible (i.e. strictly increasing or decreasing).
840 * \return \f$ f \cdot g^{-1} \f$
841 * \relates Piecewise
842 */
843 Piecewise<SBasis> pw_compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double zero);
844
845
846
847 template <typename T>
848 Piecewise<T> Piecewise<T>::operator()(SBasis f){return compose((*this),f);}
849 template <typename T>
850 Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){return compose((*this),f);}
851
852 /**
853 * ...
854 * \return ...
855 * \relates Piecewise
856 */
857 template<typename T>
858 9600 Piecewise<T> integral(Piecewise<T> const &a) {
859 9600 Piecewise<T> result;
860
1/2
✓ Branch 2 taken 9600 times.
✗ Branch 3 not taken.
 9600 result.segs.resize(a.segs.size());
861
1/2
✓ Branch 1 taken 9600 times.
✗ Branch 2 not taken.
 9600 result.cuts = a.cuts;
862 9600 typename T::output_type c = a.segs[0].at0();
863
2/2
✓ Branch 1 taken 11700 times.
✓ Branch 2 taken 9600 times.
 21300 for(unsigned i = 0; i < a.segs.size(); i++){
864
3/6
✓ Branch 6 taken 11700 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 11700 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 11700 times.
✗ Branch 14 not taken.
 11700 result.segs[i] = integral(a.segs[i])*(a.cuts[i+1]-a.cuts[i]);
865
2/4
✓ Branch 2 taken 11700 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 11700 times.
✗ Branch 7 not taken.
 11700 result.segs[i]+= c-result.segs[i].at0();
866
1/2
✓ Branch 2 taken 11700 times.
✗ Branch 3 not taken.
 11700 c = result.segs[i].at1();
867 }
868 9600 return result;
869 ✗ }
870
871 /**
872 * ...
873 * \return ...
874 * \relates Piecewise
875 */
876 template<typename T>
877 9900 Piecewise<T> derivative(Piecewise<T> const &a) {
878 9900 Piecewise<T> result;
879
1/2
✓ Branch 2 taken 9900 times.
✗ Branch 3 not taken.
 9900 result.segs.resize(a.segs.size());
880
1/2
✓ Branch 1 taken 9900 times.
✗ Branch 2 not taken.
 9900 result.cuts = a.cuts;
881
2/2
✓ Branch 1 taken 21300 times.
✓ Branch 2 taken 9900 times.
 31200 for(unsigned i = 0; i < a.segs.size(); i++){
882
3/6
✓ Branch 6 taken 21300 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 21300 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 21300 times.
✗ Branch 14 not taken.
 21300 result.segs[i] = derivative(a.segs[i])/(a.cuts[i+1]-a.cuts[i]);
883 }
884 9900 return result;
885 ✗ }
886
887 std::vector<double> roots(Piecewise<SBasis> const &f);
888
889 std::vector<std::vector<double> >multi_roots(Piecewise<SBasis> const &f, std::vector<double> const &values);
890
891 //TODO: implement level_sets directly for pwsb instead of sb (and derive it fo sb).
892 //It should be faster than the reverse as the algorithm may jump over full cut intervals.
893 std::vector<Interval> level_set(Piecewise<SBasis> const &f, Interval const &level, double tol=1e-5);
894 std::vector<Interval> level_set(Piecewise<SBasis> const &f, double v, double vtol, double tol=1e-5);
895 //std::vector<Interval> level_sets(Piecewise<SBasis> const &f, std::vector<Interval> const &levels, double tol=1e-5);
896 //std::vector<Interval> level_sets(Piecewise<SBasis> const &f, std::vector<double> &v, double vtol, double tol=1e-5);
897
898
899 /**
900 * ...
901 * \return ...
902 * \relates Piecewise
903 */
904 template<typename T>
905 ✗ Piecewise<T> reverse(Piecewise<T> const &f) {
906 ✗ Piecewise<T> ret = Piecewise<T>();
907 ✗ ret.reserve(f.size());
908 ✗ double start = f.cuts[0];
909 ✗ double end = f.cuts.back();
910 ✗ for (unsigned i = 0; i < f.cuts.size(); i++) {
911 ✗ double x = f.cuts[f.cuts.size() - 1 - i];
912 ✗ ret.push_cut(end - (x - start));
913 }
914 ✗ for (unsigned i = 0; i < f.segs.size(); i++)
915 ✗ ret.push_seg(reverse(f[f.segs.size() - i - 1]));
916 ✗ return ret;
917 ✗ }
918
919 /**
920 * Interpolates between a and b.
921 * \return a if t = 0, b if t = 1, or an interpolation between a and b for t in [0,1]
922 * \relates Piecewise
923 */
924 template<typename T>
925 Piecewise<T> lerp(double t, Piecewise<T> const &a, Piecewise<T> b) {
926 // Make sure both paths have the same number of segments and cuts at the same locations
927 b.setDomain(a.domain());
928 Piecewise<T> pA = partition(a, b.cuts);
929 Piecewise<T> pB = partition(b, a.cuts);
930
931 return (pA*(1-t) + pB*t);
932 }
933
934 }
935 #endif //LIB2GEOM_SEEN_PIECEWISE_H
936 /*
937 Local Variables:
938 mode:c++
939 c-file-style:"stroustrup"
940 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
941 indent-tabs-mode:nil
942 fill-column:99
943 End:
944 */
945 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
946