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