1 /*
2 * piecewise.h - 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 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 SEEN_GEOM_PW_SB_H
32 #define SEEN_GEOM_PW_SB_H
33
34 #include "sbasis.h"
35 #include <vector>
36 #include <map>
37
38 #include "concepts.h"
39 #include "isnan.h"
40 #include <boost/concept_check.hpp>
41
42 namespace Geom {
43
44 template <typename T>
45 class Piecewise {
46 BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept);
47
48 public:
49 std::vector<double> cuts;
50 std::vector<T> segs;
51 //segs[i] stretches from cuts[i] to cuts[i+1].
52
Piecewise()53 Piecewise() {}
54
Piecewise(const T & s)55 explicit Piecewise(const T &s) {
56 push_cut(0.);
57 push_seg(s);
58 push_cut(1.);
59 }
60
61 typedef typename T::output_type output_type;
62
Piecewise(const output_type & v)63 explicit Piecewise(const output_type & v) {
64 push_cut(0.);
65 push_seg(T(v));
66 push_cut(1.);
67 }
68
69 inline T operator[](unsigned i) const { return segs[i]; }
70 inline T &operator[](unsigned i) { return segs[i]; }
operator()71 inline output_type operator()(double t) const { return valueAt(t); }
valueAt(double t)72 inline output_type valueAt(double t) const {
73 unsigned n = segN(t);
74 return segs[n](segT(t, n));
75 }
76 //TODO: maybe it is not a good idea to have this?
77 Piecewise<T> operator()(SBasis f);
78 Piecewise<T> operator()(Piecewise<SBasis>f);
79
size()80 inline unsigned size() const { return segs.size(); }
empty()81 inline bool empty() const { return segs.empty(); }
82
83 /**Convenience/implementation hiding function to add segment/cut pairs.
84 * Asserts that basic size and order invariants are correct
85 */
push(const T & s,double to)86 inline void push(const T &s, double to) {
87 assert(cuts.size() - segs.size() == 1);
88 push_seg(s);
89 push_cut(to);
90 }
91 //Convenience/implementation hiding function to add cuts.
push_cut(double c)92 inline void push_cut(double c) {
93 assert_invariants(cuts.empty() || c > cuts.back());
94 cuts.push_back(c);
95 }
96 //Convenience/implementation hiding function to add segments.
push_seg(const T & s)97 inline void push_seg(const T &s) { segs.push_back(s); }
98
99 /**Returns the segment index which corresponds to a 'global' piecewise time.
100 * Also takes optional low/high parameters to expedite the search for the segment.
101 */
102 inline unsigned segN(double t, int low = 0, int high = -1) const {
103 high = (high == -1) ? size() : high;
104 if(t < cuts[0]) return 0;
105 if(t >= cuts[size()]) return size() - 1;
106 while(low < high) {
107 int mid = (high + low) / 2; //Lets not plan on having huge (> INT_MAX / 2) cut sequences
108 double mv = cuts[mid];
109 if(mv < t) {
110 if(t < cuts[mid + 1]) return mid; else low = mid + 1;
111 } else if(t < mv) {
112 if(cuts[mid - 1] < t) return mid - 1; else high = mid - 1;
113 } else {
114 return mid;
115 }
116 }
117 return low;
118 }
119
120 /**Returns the time within a segment, given the 'global' piecewise time.
121 * Also takes an optional index parameter which may be used for efficiency or to find the time on a
122 * segment outside its range. If it is left to its default, -1, it will call segN to find the index.
123 */
124 inline double segT(double t, int i = -1) const {
125 if(i == -1) i = segN(t);
126 assert(i >= 0);
127 return (t - cuts[i]) / (cuts[i+1] - cuts[i]);
128 }
129
mapToDomain(double t,unsigned i)130 inline double mapToDomain(double t, unsigned i) const {
131 return (1-t)*cuts[i] + t*cuts[i+1]; //same as: t * (cuts[i+1] - cuts[i]) + cuts[i]
132 }
133
134 //Offsets the piecewise domain
offsetDomain(double o)135 inline void offsetDomain(double o) {
136 assert(is_finite(o));
137 if(o != 0)
138 for(unsigned i = 0; i <= size(); i++)
139 cuts[i] += o;
140 }
141
142 //Scales the domain of the function by a value. 0 will result in an empty Piecewise.
scaleDomain(double s)143 inline void scaleDomain(double s) {
144 assert(s > 0);
145 if(s == 0) {
146 cuts.clear(); segs.clear();
147 return;
148 }
149 for(unsigned i = 0; i <= size(); i++)
150 cuts[i] *= s;
151 }
152
153 //Retrieves the domain in interval form
domain()154 inline Interval domain() const { return Interval(cuts.front(), cuts.back()); }
155
156 //Transforms the domain into another interval
setDomain(Interval dom)157 inline void setDomain(Interval dom) {
158 if(empty()) return;
159 if(dom.isEmpty()) {
160 cuts.clear(); segs.clear();
161 return;
162 }
163 double cf = cuts.front();
164 double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);
165 for(unsigned i = 0; i <= size(); i++)
166 cuts[i] = (cuts[i] - cf) * s + o;
167 }
168
169 //Concatenates this Piecewise function with another, offseting time of the other to match the end.
concat(const Piecewise<T> & other)170 inline void concat(const Piecewise<T> &other) {
171 if(other.empty()) return;
172
173 if(empty()) {
174 cuts = other.cuts; segs = other.segs;
175 return;
176 }
177
178 segs.insert(segs.end(), other.segs.begin(), other.segs.end());
179 double t = cuts.back() - other.cuts.front();
180 for(unsigned i = 0; i < other.size(); i++)
181 push_cut(other.cuts[i + 1] + t);
182 }
183
184 //Like concat, but ensures continuity.
continuousConcat(const Piecewise<T> & other)185 inline void continuousConcat(const Piecewise<T> &other) {
186 boost::function_requires<AddableConcept<typename T::output_type> >();
187 if(other.empty()) return;
188 typename T::output_type y = segs.back().at1() - other.segs.front().at0();
189
190 if(empty()) {
191 for(unsigned i = 0; i < other.size(); i++)
192 push_seg(other[i] + y);
193 cuts = other.cuts;
194 return;
195 }
196
197 double t = cuts.back() - other.cuts.front();
198 for(unsigned i = 0; i < other.size(); i++)
199 push(other[i] + y, other.cuts[i + 1] + t);
200 }
201
202 //returns true if the Piecewise<T> meets some basic invariants.
invariants()203 inline bool invariants() const {
204 // segs between cuts
205 if(!(segs.size() + 1 == cuts.size() || (segs.empty() && cuts.empty())))
206 return false;
207 // cuts in order
208 for(unsigned i = 0; i < segs.size(); i++)
209 if(cuts[i] >= cuts[i+1])
210 return false;
211 return true;
212 }
213
214 };
215
216 template<typename T>
bounds_fast(const Piecewise<T> & f)217 inline typename FragmentConcept<T>::BoundsType bounds_fast(const Piecewise<T> &f) {
218 boost::function_requires<FragmentConcept<T> >();
219
220 if(f.empty()) return typename FragmentConcept<T>::BoundsType();
221 typename FragmentConcept<T>::BoundsType ret(bounds_fast(f[0]));
222 for(unsigned i = 1; i < f.size(); i++)
223 ret.unionWith(bounds_fast(f[i]));
224 return ret;
225 }
226
227 template<typename T>
bounds_exact(const Piecewise<T> & f)228 inline typename FragmentConcept<T>::BoundsType bounds_exact(const Piecewise<T> &f) {
229 boost::function_requires<FragmentConcept<T> >();
230
231 if(f.empty()) return typename FragmentConcept<T>::BoundsType();
232 typename FragmentConcept<T>::BoundsType ret(bounds_exact(f[0]));
233 for(unsigned i = 1; i < f.size(); i++)
234 ret.unionWith(bounds_exact(f[i]));
235 return ret;
236 }
237
238 template<typename T>
bounds_local(const Piecewise<T> & f,const Interval & m)239 inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const Interval &m) {
240 boost::function_requires<FragmentConcept<T> >();
241
242 if(f.empty()) return typename FragmentConcept<T>::BoundsType();
243 if(m.isEmpty()) return typename FragmentConcept<T>::BoundsType(f(m.min()));
244
245 unsigned fi = f.segN(m.min()), ti = f.segN(m.max());
246 double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);
247
248 if(fi == ti) return bounds_local(f[fi], Interval(ft, tt));
249
250 typename FragmentConcept<T>::BoundsType ret(bounds_local(f[fi], Interval(ft, 1.)));
251 for(unsigned i = fi + 1; i < ti; i++)
252 ret.unionWith(bounds_exact(f[i]));
253 if(tt != 0.) ret.unionWith(bounds_local(f[ti], Interval(0., tt)));
254
255 return ret;
256 }
257
258 //returns a portion of a piece of a Piecewise<T>, given the piece's index and a to/from time.
259 template<typename T>
elem_portion(const Piecewise<T> & a,unsigned i,double from,double to)260 T elem_portion(const Piecewise<T> &a, unsigned i, double from, double to) {
261 assert(i < a.size());
262 double rwidth = 1 / (a.cuts[i+1] - a.cuts[i]);
263 return portion( a[i], (from - a.cuts[i]) * rwidth, (to - a.cuts[i]) * rwidth );
264 }
265
266 /**Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c);
267 * Further subdivides the Piecewise<T> such that there is a cut at every value in c.
268 * Precondition: c sorted lower to higher.
269 *
270 * //Given Piecewise<T> a and b:
271 * Piecewise<T> ac = a.partition(b.cuts);
272 * Piecewise<T> bc = b.partition(a.cuts);
273 * //ac.cuts should be equivalent to bc.cuts
274 */
275 template<typename T>
partition(const Piecewise<T> & pw,std::vector<double> const & c)276 Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c) {
277 assert(pw.invariants());
278 if(c.empty()) return Piecewise<T>(pw);
279
280 Piecewise<T> ret = Piecewise<T>();
281 ret.cuts.reserve(c.size() + pw.cuts.size());
282 ret.segs.reserve(c.size() + pw.cuts.size() - 1);
283
284 if(pw.empty()) {
285 ret.cuts = c;
286 for(unsigned i = 0; i < c.size() - 1; i++)
287 ret.push_seg(T());
288 return ret;
289 }
290
291 unsigned si = 0, ci = 0; //Segment index, Cut index
292
293 //if the cuts have something earlier than the Piecewise<T>, add portions of the first segment
294 while(ci < c.size() && c[ci] < pw.cuts.front()) {
295 bool isLast = (ci == c.size()-1 || c[ci + 1] >= pw.cuts.front());
296 ret.push_cut(c[ci]);
297 ret.push_seg( elem_portion(pw, 0, c[ci], isLast ? pw.cuts.front() : c[ci + 1]) );
298 ci++;
299 }
300
301 ret.push_cut(pw.cuts[0]);
302 double prev = pw.cuts[0]; //previous cut
303 //Loop which handles cuts within the Piecewise<T> domain
304 //Should have the cuts = segs + 1 invariant
305 while(si < pw.size() && ci <= c.size()) {
306 if(ci == c.size() && prev <= pw.cuts[si]) { //cuts exhausted, straight copy the rest
307 ret.segs.insert(ret.segs.end(), pw.segs.begin() + si, pw.segs.end());
308 ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + si + 1, pw.cuts.end());
309 return ret;
310 }else if(ci == c.size() || c[ci] >= pw.cuts[si + 1]) { //no more cuts within this segment, finalize
311 if(prev > pw.cuts[si]) { //segment already has cuts, so portion is required
312 ret.push_seg(portion(pw[si], pw.segT(prev, si), 1.0));
313 } else { //plain copy is fine
314 ret.push_seg(pw[si]);
315 }
316 ret.push_cut(pw.cuts[si + 1]);
317 prev = pw.cuts[si + 1];
318 si++;
319 } else if(c[ci] == pw.cuts[si]){ //coincident
320 //Already finalized the seg with the code immediately above
321 ci++;
322 } else { //plain old subdivision
323 ret.push(elem_portion(pw, si, prev, c[ci]), c[ci]);
324 prev = c[ci];
325 ci++;
326 }
327 }
328
329 //input cuts extend further than this Piecewise<T>, extend the last segment.
330 while(ci < c.size()) {
331 if(c[ci] > prev) {
332 ret.push(elem_portion(pw, pw.size() - 1, prev, c[ci]), c[ci]);
333 prev = c[ci];
334 }
335 ci++;
336 }
337 return ret;
338 }
339
340 /**Piecewise<T> portion(const Piecewise<T> &pw, double from, double to);
341 * Returns a Piecewise<T> with a defined domain of [min(from, to), max(from, to)].
342 */
343 template<typename T>
portion(const Piecewise<T> & pw,double from,double to)344 Piecewise<T> portion(const Piecewise<T> &pw, double from, double to) {
345 if(pw.empty() || from == to) return Piecewise<T>();
346
347 Piecewise<T> ret;
348
349 double temp = from;
350 from = std::min(from, to);
351 to = std::max(temp, to);
352
353 unsigned i = pw.segN(from);
354 ret.push_cut(from);
355 if(i == pw.size() - 1 || to < pw.cuts[i + 1]) { //to/from inhabit the same segment
356 ret.push(elem_portion(pw, i, from, to), to);
357 return ret;
358 }
359 ret.push_seg(portion( pw[i], pw.segT(from, i), 1.0 ));
360 i++;
361 unsigned fi = pw.segN(to, i);
362
363 ret.segs.insert(ret.segs.end(), pw.segs.begin() + i, pw.segs.begin() + fi); //copy segs
364 ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + i, pw.cuts.begin() + fi + 1); //and their cuts
365
366 ret.push_seg( portion(pw[fi], 0.0, pw.segT(to, fi)));
367 if(to != ret.cuts.back()) ret.push_cut(to);
368 ret.invariants();
369 return ret;
370 }
371
372 template<typename T>
remove_short_cuts(Piecewise<T> const & f,double tol)373 Piecewise<T> remove_short_cuts(Piecewise<T> const &f, double tol) {
374 if(f.empty()) return f;
375 Piecewise<T> ret;
376 ret.push_cut(f.cuts[0]);
377 for(unsigned i=0; i<f.size(); i++){
378 if (f.cuts[i+1]-f.cuts[i] >= tol || i==f.size()-1) {
379 ret.push(f[i], f.cuts[i+1]);
380 }
381 }
382 return ret;
383 }
384
385 template<typename T>
remove_short_cuts_extending(Piecewise<T> const & f,double tol)386 Piecewise<T> remove_short_cuts_extending(Piecewise<T> const &f, double tol) {
387 if(f.empty()) return f;
388 Piecewise<T> ret;
389 ret.push_cut(f.cuts[0]);
390 double last = f.cuts[0]; // last cut included
391 for(unsigned i=0; i<f.size(); i++){
392 if (f.cuts[i+1]-f.cuts[i] >= tol) {
393 ret.push(elem_portion(f, i, last, f.cuts[i+1]), f.cuts[i+1]);
394 last = f.cuts[i+1];
395 }
396 }
397 return ret;
398 }
399
400 template<typename T>
roots(const Piecewise<T> & pw)401 std::vector<double> roots(const Piecewise<T> &pw) {
402 std::vector<double> ret;
403 for(unsigned i = 0; i < pw.size(); i++) {
404 std::vector<double> sr = roots(pw[i]);
405 for (unsigned j = 0; j < sr.size(); j++) ret.push_back(sr[j] * (pw.cuts[i + 1] - pw.cuts[i]) + pw.cuts[i]);
406
407 }
408 return ret;
409 }
410
411 //IMPL: OffsetableConcept
412 template<typename T>
413 Piecewise<T> operator+(Piecewise<T> const &a, typename T::output_type b) {
414 boost::function_requires<OffsetableConcept<T> >();
415 //TODO:empty
416 Piecewise<T> ret = Piecewise<T>();
417 ret.cuts = a.cuts;
418 for(unsigned i = 0; i < a.size();i++)
419 ret.push_seg(a[i] + b);
420 return ret;
421 }
422 template<typename T>
423 Piecewise<T> operator-(Piecewise<T> const &a, typename T::output_type b) {
424 boost::function_requires<OffsetableConcept<T> >();
425 //TODO: empty
426 Piecewise<T> ret = Piecewise<T>();
427 ret.cuts = a.cuts;
428 for(unsigned i = 0; i < a.size();i++)
429 ret.push_seg(a[i] - b);
430 return ret;
431 }
432 template<typename T>
433 Piecewise<T> operator+=(Piecewise<T>& a, typename T::output_type b) {
434 boost::function_requires<OffsetableConcept<T> >();
435
436 if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
437
438 for(unsigned i = 0; i < a.size();i++)
439 a[i] += b;
440 return a;
441 }
442 template<typename T>
443 Piecewise<T> operator-=(Piecewise<T>& a, typename T::output_type b) {
444 boost::function_requires<OffsetableConcept<T> >();
445
446 if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
447
448 for(unsigned i = 0;i < a.size();i++)
449 a[i] -= b;
450 return a;
451 }
452
453 //IMPL: ScalableConcept
454 template<typename T>
455 Piecewise<T> operator-(Piecewise<T> const &a) {
456 boost::function_requires<ScalableConcept<T> >();
457
458 Piecewise<T> ret;
459 ret.cuts = a.cuts;
460 for(unsigned i = 0; i < a.size();i++)
461 ret.push_seg(- a[i]);
462 return ret;
463 }
464 template<typename T>
465 Piecewise<T> operator*(Piecewise<T> const &a, double b) {
466 boost::function_requires<ScalableConcept<T> >();
467
468 if(a.empty()) return Piecewise<T>();
469
470 Piecewise<T> ret;
471 ret.cuts = a.cuts;
472 for(unsigned i = 0; i < a.size();i++)
473 ret.push_seg(a[i] * b);
474 return ret;
475 }
476 template<typename T>
477 Piecewise<T> operator/(Piecewise<T> const &a, double b) {
478 boost::function_requires<ScalableConcept<T> >();
479
480 //FIXME: b == 0?
481 if(a.empty()) return Piecewise<T>();
482
483 Piecewise<T> ret;
484 ret.cuts = a.cuts;
485 for(unsigned i = 0; i < a.size();i++)
486 ret.push_seg(a[i] / b);
487 return ret;
488 }
489 template<typename T>
490 Piecewise<T> operator*=(Piecewise<T>& a, double b) {
491 boost::function_requires<ScalableConcept<T> >();
492
493 if(a.empty()) return Piecewise<T>();
494
495 for(unsigned i = 0; i < a.size();i++)
496 a[i] *= b;
497 return a;
498 }
499 template<typename T>
500 Piecewise<T> operator/=(Piecewise<T>& a, double b) {
501 boost::function_requires<ScalableConcept<T> >();
502
503 //FIXME: b == 0?
504 if(a.empty()) return Piecewise<T>();
505
506 for(unsigned i = 0; i < a.size();i++)
507 a[i] /= b;
508 return a;
509 }
510
511 //IMPL: AddableConcept
512 template<typename T>
513 Piecewise<T> operator+(Piecewise<T> const &a, Piecewise<T> const &b) {
514 boost::function_requires<AddableConcept<T> >();
515
516 Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
517 Piecewise<T> ret = Piecewise<T>();
518 assert(pa.size() == pb.size());
519 ret.cuts = pa.cuts;
520 for (unsigned i = 0; i < pa.size(); i++)
521 ret.push_seg(pa[i] + pb[i]);
522 return ret;
523 }
524 template<typename T>
525 Piecewise<T> operator-(Piecewise<T> const &a, Piecewise<T> const &b) {
526 boost::function_requires<AddableConcept<T> >();
527
528 Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
529 Piecewise<T> ret = Piecewise<T>();
530 assert(pa.size() == pb.size());
531 ret.cuts = pa.cuts;
532 for (unsigned i = 0; i < pa.size(); i++)
533 ret.push_seg(pa[i] - pb[i]);
534 return ret;
535 }
536 template<typename T>
537 inline Piecewise<T> operator+=(Piecewise<T> &a, Piecewise<T> const &b) {
538 a = a+b;
539 return a;
540 }
541 template<typename T>
542 inline Piecewise<T> operator-=(Piecewise<T> &a, Piecewise<T> const &b) {
543 a = a-b;
544 return a;
545 }
546
547 template<typename T1,typename T2>
548 Piecewise<T2> operator*(Piecewise<T1> const &a, Piecewise<T2> const &b) {
549 //function_requires<MultiplicableConcept<T1> >();
550 //function_requires<MultiplicableConcept<T2> >();
551
552 Piecewise<T1> pa = partition(a, b.cuts);
553 Piecewise<T2> pb = partition(b, a.cuts);
554 Piecewise<T2> ret = Piecewise<T2>();
555 assert(pa.size() == pb.size());
556 ret.cuts = pa.cuts;
557 for (unsigned i = 0; i < pa.size(); i++)
558 ret.push_seg(pa[i] * pb[i]);
559 return ret;
560 }
561
562 template<typename T>
563 inline Piecewise<T> operator*=(Piecewise<T> &a, Piecewise<T> const &b) {
564 a = a * b;
565 return a;
566 }
567
568 Piecewise<SBasis> divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, unsigned k);
569 //TODO: replace divide(a,b,k) by divide(a,b,tol,k)?
570 //TODO: atm, relative error is ~(tol/a)%. Find a way to make it independant of a.
571 //Nota: the result is 'truncated' where b is smaller than 'zero': ~ a/max(b,zero).
572 Piecewise<SBasis>
573 divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
574 Piecewise<SBasis>
575 divide(SBasis const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
576 Piecewise<SBasis>
577 divide(Piecewise<SBasis> const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
578 Piecewise<SBasis>
579 divide(SBasis const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
580
581 //Composition: functions called compose_* are pieces of compose that are factored out in pw.cpp.
582 std::map<double,unsigned> compose_pullback(std::vector<double> const &cuts, SBasis const &g);
583 int compose_findSegIdx(std::map<double,unsigned>::iterator const &cut,
584 std::map<double,unsigned>::iterator const &next,
585 std::vector<double> const &levels,
586 SBasis const &g);
587
588 //TODO: add concept check
589 template<typename T>
compose(Piecewise<T> const & f,SBasis const & g)590 Piecewise<T> compose(Piecewise<T> const &f, SBasis const &g){
591 Piecewise<T> result;
592 if (f.empty()) return result;
593 if (g.isZero()) return Piecewise<T>(f(0));
594 if (f.size()==1){
595 double t0 = f.cuts[0], width = f.cuts[1] - t0;
596 return (Piecewise<T>) compose(f.segs[0],compose(Linear(-t0 / width, (1-t0) / width), g));
597 }
598
599 //first check bounds...
600 Interval bs = bounds_fast(g);
601 if (f.cuts.front() > bs.max() || bs.min() > f.cuts.back()){
602 int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;
603 double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;
604 return (Piecewise<T>) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g));
605 }
606
607 std::vector<double> levels;//we can forget first and last cuts...
608 levels.insert(levels.begin(),f.cuts.begin()+1,f.cuts.end()-1);
609 //TODO: use a std::vector<pairs<double,unsigned> > instead of a map<double,unsigned>.
610 std::map<double,unsigned> cuts_pb = compose_pullback(levels,g);
611
612 //-- Compose each piece of g with the relevant seg of f.
613 result.cuts.push_back(0.);
614 std::map<double,unsigned>::iterator cut=cuts_pb.begin();
615 std::map<double,unsigned>::iterator next=cut; next++;
616 while(next!=cuts_pb.end()){
617 //assert(std::abs(int((*cut).second-(*next).second))<1);
618 //TODO: find a way to recover from this error? the root finder missed some root;
619 // the levels/variations of f might be too close/fast...
620 int idx = compose_findSegIdx(cut,next,levels,g);
621 double t0=(*cut).first;
622 double t1=(*next).first;
623
624 SBasis sub_g=compose(g, Linear(t0,t1));
625 sub_g=compose(Linear(-f.cuts[idx]/(f.cuts[idx+1]-f.cuts[idx]),
626 (1-f.cuts[idx])/(f.cuts[idx+1]-f.cuts[idx])),sub_g);
627 result.push(compose(f[idx],sub_g),t1);
628 cut++;
629 next++;
630 }
631 return(result);
632 }
633
634 //TODO: add concept check for following composition functions
635 template<typename T>
compose(Piecewise<T> const & f,Piecewise<SBasis> const & g)636 Piecewise<T> compose(Piecewise<T> const &f, Piecewise<SBasis> const &g){
637 Piecewise<T> result;
638 for(unsigned i = 0; i < g.segs.size(); i++){
639 Piecewise<T> fgi=compose(f, g.segs[i]);
640 fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1]));
641 result.concat(fgi);
642 }
643 return result;
644 }
645
646 template <typename T>
operator()647 Piecewise<T> Piecewise<T>::operator()(SBasis f){return compose((*this),f);}
648 template <typename T>
operator()649 Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){return compose((*this),f);}
650
651 template<typename T>
integral(Piecewise<T> const & a)652 Piecewise<T> integral(Piecewise<T> const &a) {
653 Piecewise<T> result;
654 result.segs.resize(a.segs.size());
655 result.cuts = a.cuts;
656 typename T::output_type c = a.segs[0].at0();
657 for(unsigned i = 0; i < a.segs.size(); i++){
658 result.segs[i] = integral(a.segs[i])*(a.cuts[i+1]-a.cuts[i]);
659 result.segs[i]+= c-result.segs[i].at0();
660 c = result.segs[i].at1();
661 }
662 return result;
663 }
664
665 template<typename T>
derivative(Piecewise<T> const & a)666 Piecewise<T> derivative(Piecewise<T> const &a) {
667 Piecewise<T> result;
668 result.segs.resize(a.segs.size());
669 result.cuts = a.cuts;
670 for(unsigned i = 0; i < a.segs.size(); i++){
671 result.segs[i] = derivative(a.segs[i])/(a.cuts[i+1]-a.cuts[i]);
672 }
673 return result;
674 }
675
676 std::vector<double> roots(Piecewise<SBasis> const &f);
677
678 }
679
680 #endif //SEEN_GEOM_PW_SB_H
681 /*
682 Local Variables:
683 mode:c++
684 c-file-style:"stroustrup"
685 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
686 indent-tabs-mode:nil
687 fill-column:99
688 End:
689 */
690 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
691