1 /**
2 * \file
3 * \brief Piecewise function class
4 *
5 * Copyright 2007 Michael Sloan <mgsloan@gmail.com>
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it either under the terms of the GNU Lesser General Public
9 * License version 2.1 as published by the Free Software Foundation
10 * (the "LGPL") or, at your option, under the terms of the Mozilla
11 * Public License Version 1.1 (the "MPL"). If you do not alter this
12 * notice, a recipient may use your version of this file under either
13 * the MPL or the LGPL.
14 *
15 * You should have received a copy of the LGPL along with this library
16 * in the file COPYING-LGPL-2.1; if not, output to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 * You should have received a copy of the MPL along with this library
19 * in the file COPYING-MPL-1.1
20 *
21 * The contents of this file are subject to the Mozilla Public License
22 * Version 1.1 (the "License"); you may not use this file except in
23 * compliance with the License. You may obtain a copy of the License at
24 * http://www.mozilla.org/MPL/
25 *
26 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
27 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
28 * the specific language governing rights and limitations.
29 *
30 */
32 #ifndef SEEN_GEOM_PW_SB_H
33 #define SEEN_GEOM_PW_SB_H
35 #include <2geom/sbasis.h>
36 #include <vector>
37 #include <map>
39 #include <2geom/concepts.h>
40 #include <2geom/isnan.h>
41 #include <boost/concept_check.hpp>
43 namespace Geom {
44 /**
45 * %Piecewise function class.
46 * The Piecewise class manages a sequence of elements of a type as segments and
47 * the ’cuts’ between them. These cuts are time values which separate the pieces.
48 * This function representation allows for more interesting functions, as it provides
49 * a viable output for operations such as inversion, which may require multiple
50 * SBasis to properly invert the original.
51 * As for technical details, while the actual SBasis segments begin on the first
52 * cut and end on the last, the function is defined throughout all inputs by ex-
53 * tending the first and last segments. The exact switching between segments is
54 * arbitrarily such that beginnings (t=0) have preference over endings (t=1). This
55 * only matters if it is discontinuous at the location.
56 * \f[
57 * f(t) \rightarrow \left\{
58 * \begin{array}{cc}
59 * s_1,& t <= c_2 \\
60 * s_2,& c_2 <= t <= c_3\\
61 * \ldots
62 * s_n,& c_n <= t
63 * \end{array}\right.
64 * \f]
65 */
66 template <typename T>
67 class Piecewise {
68 BOOST_CLASS_REQUIRE(T, Geom, FragmentConcept);
70 public:
71 std::vector<double> cuts;
72 std::vector<T> segs;
73 //segs[i] stretches from cuts[i] to cuts[i+1].
75 Piecewise() {}
77 explicit Piecewise(const T &s) {
78 push_cut(0.);
79 push_seg(s);
80 push_cut(1.);
81 }
83 unsigned input_dim(){return 1;}
85 typedef typename T::output_type output_type;
87 explicit Piecewise(const output_type & v) {
88 push_cut(0.);
89 push_seg(T(v));
90 push_cut(1.);
91 }
93 inline T operator[](unsigned i) const { return segs[i]; }
94 inline T &operator[](unsigned i) { return segs[i]; }
95 inline output_type operator()(double t) const { return valueAt(t); }
96 inline output_type valueAt(double t) const {
97 unsigned n = segN(t);
98 return segs[n](segT(t, n));
99 }
100 inline output_type firstValue() const {
101 return valueAt(cuts.front());
102 }
103 inline output_type lastValue() const {
104 return valueAt(cuts.back());
105 }
106 std::vector<output_type> valueAndDerivatives(double t, unsigned n_derivs) const {
107 unsigned n = segN(t);
108 std::vector<output_type> ret, val = segs[n].valueAndDerivatives(segT(t, n), n_derivs);
109 double mult = 1;
110 for(unsigned i = 0; i < val.size(); i++) {
111 ret.push_back(val[i] * mult);
112 mult /= cuts[n+1] - cuts[n];
113 }
114 return ret;
115 }
116 //TODO: maybe it is not a good idea to have this?
117 Piecewise<T> operator()(SBasis f);
118 Piecewise<T> operator()(Piecewise<SBasis>f);
120 inline unsigned size() const { return segs.size(); }
121 inline bool empty() const { return segs.empty(); }
122 inline void clear() {
123 segs.clear();
124 cuts.clear();
125 }
127 /**Convenience/implementation hiding function to add segment/cut pairs.
128 * Asserts that basic size and order invariants are correct
129 */
130 inline void push(const T &s, double to) {
131 assert(cuts.size() - segs.size() == 1);
132 push_seg(s);
133 push_cut(to);
134 }
135 //Convenience/implementation hiding function to add cuts.
136 inline void push_cut(double c) {
137 ASSERT_INVARIANTS(cuts.empty() || c > cuts.back());
138 cuts.push_back(c);
139 }
140 //Convenience/implementation hiding function to add segments.
141 inline void push_seg(const T &s) { segs.push_back(s); }
143 /**Returns the segment index which corresponds to a 'global' piecewise time.
144 * Also takes optional low/high parameters to expedite the search for the segment.
145 */
146 inline unsigned segN(double t, int low = 0, int high = -1) const {
147 high = (high == -1) ? size() : high;
148 if(t < cuts[0]) return 0;
149 if(t >= cuts[size()]) return size() - 1;
150 while(low < high) {
151 int mid = (high + low) / 2; //Lets not plan on having huge (> INT_MAX / 2) cut sequences
152 double mv = cuts[mid];
153 if(mv < t) {
154 if(t < cuts[mid + 1]) return mid; else low = mid + 1;
155 } else if(t < mv) {
156 if(cuts[mid - 1] < t) return mid - 1; else high = mid - 1;
157 } else {
158 return mid;
159 }
160 }
161 return low;
162 }
164 /**Returns the time within a segment, given the 'global' piecewise time.
165 * Also takes an optional index parameter which may be used for efficiency or to find the time on a
166 * segment outside its range. If it is left to its default, -1, it will call segN to find the index.
167 */
168 inline double segT(double t, int i = -1) const {
169 if(i == -1) i = segN(t);
170 assert(i >= 0);
171 return (t - cuts[i]) / (cuts[i+1] - cuts[i]);
172 }
174 inline double mapToDomain(double t, unsigned i) const {
175 return (1-t)*cuts[i] + t*cuts[i+1]; //same as: t * (cuts[i+1] - cuts[i]) + cuts[i]
176 }
178 //Offsets the piecewise domain
179 inline void offsetDomain(double o) {
180 assert(IS_FINITE(o));
181 if(o != 0)
182 for(unsigned i = 0; i <= size(); i++)
183 cuts[i] += o;
184 }
186 //Scales the domain of the function by a value. 0 will result in an empty Piecewise.
187 inline void scaleDomain(double s) {
188 assert(s > 0);
189 if(s == 0) {
190 cuts.clear(); segs.clear();
191 return;
192 }
193 for(unsigned i = 0; i <= size(); i++)
194 cuts[i] *= s;
195 }
197 //Retrieves the domain in interval form
198 inline Interval domain() const { return Interval(cuts.front(), cuts.back()); }
200 //Transforms the domain into another interval
201 inline void setDomain(Interval dom) {
202 if(empty()) return;
203 /* dom can not be empty
204 if(dom.isEmpty()) {
205 cuts.clear(); segs.clear();
206 return;
207 }*/
208 double cf = cuts.front();
209 double o = dom.min() - cf, s = dom.extent() / (cuts.back() - cf);
210 for(unsigned i = 0; i <= size(); i++)
211 cuts[i] = (cuts[i] - cf) * s + o;
212 //fix floating point precision errors.
213 cuts[0] = dom.min();
214 cuts[size()] = dom.max();
215 }
217 //Concatenates this Piecewise function with another, offseting time of the other to match the end.
218 inline void concat(const Piecewise<T> &other) {
219 if(other.empty()) return;
221 if(empty()) {
222 cuts = other.cuts; segs = other.segs;
223 return;
224 }
226 segs.insert(segs.end(), other.segs.begin(), other.segs.end());
227 double t = cuts.back() - other.cuts.front();
228 for(unsigned i = 0; i < other.size(); i++)
229 push_cut(other.cuts[i + 1] + t);
230 }
232 //Like concat, but ensures continuity.
233 inline void continuousConcat(const Piecewise<T> &other) {
234 boost::function_requires<AddableConcept<typename T::output_type> >();
235 if(other.empty()) return;
237 if(empty()) {
238 for(unsigned i = 0; i < other.size(); i++)
239 push_seg(other[i]);
240 cuts = other.cuts;
241 return;
242 }
244 typename T::output_type y = segs.back().at1() - other.segs.front().at0();
245 double t = cuts.back() - other.cuts.front();
246 for(unsigned i = 0; i < other.size(); i++)
247 push(other[i] + y, other.cuts[i + 1] + t);
248 }
250 //returns true if the Piecewise<T> meets some basic invariants.
251 inline bool invariants() const {
252 // segs between cuts
253 if(!(segs.size() + 1 == cuts.size() || (segs.empty() && cuts.empty())))
254 return false;
255 // cuts in order
256 for(unsigned i = 0; i < segs.size(); i++)
257 if(cuts[i] >= cuts[i+1])
258 return false;
259 return true;
260 }
262 };
264 template<typename T>
265 inline typename FragmentConcept<T>::BoundsType bounds_fast(const Piecewise<T> &f) {
266 boost::function_requires<FragmentConcept<T> >();
268 if(f.empty()) return typename FragmentConcept<T>::BoundsType();
269 typename FragmentConcept<T>::BoundsType ret(bounds_fast(f[0]));
270 for(unsigned i = 1; i < f.size(); i++)
271 ret.unionWith(bounds_fast(f[i]));
272 return ret;
273 }
275 template<typename T>
276 inline typename FragmentConcept<T>::BoundsType bounds_exact(const Piecewise<T> &f) {
277 boost::function_requires<FragmentConcept<T> >();
279 if(f.empty()) return typename FragmentConcept<T>::BoundsType();
280 typename FragmentConcept<T>::BoundsType ret(bounds_exact(f[0]));
281 for(unsigned i = 1; i < f.size(); i++)
282 ret.unionWith(bounds_exact(f[i]));
283 return ret;
284 }
286 template<typename T>
287 inline typename FragmentConcept<T>::BoundsType bounds_local(const Piecewise<T> &f, const OptInterval &_m) {
288 boost::function_requires<FragmentConcept<T> >();
290 if(f.empty() || !_m) return typename FragmentConcept<T>::BoundsType();
291 Interval const &m = *_m;
292 if(m.isSingular()) return typename FragmentConcept<T>::BoundsType(f(m.min()));
294 unsigned fi = f.segN(m.min()), ti = f.segN(m.max());
295 double ft = f.segT(m.min(), fi), tt = f.segT(m.max(), ti);
297 if(fi == ti) return bounds_local(f[fi], Interval(ft, tt));
299 typename FragmentConcept<T>::BoundsType ret(bounds_local(f[fi], Interval(ft, 1.)));
300 for(unsigned i = fi + 1; i < ti; i++)
301 ret.unionWith(bounds_exact(f[i]));
302 if(tt != 0.) ret.unionWith(bounds_local(f[ti], Interval(0., tt)));
304 return ret;
305 }
307 //returns a portion of a piece of a Piecewise<T>, given the piece's index and a to/from time.
308 template<typename T>
309 T elem_portion(const Piecewise<T> &a, unsigned i, double from, double to) {
310 assert(i < a.size());
311 double rwidth = 1 / (a.cuts[i+1] - a.cuts[i]);
312 return portion( a[i], (from - a.cuts[i]) * rwidth, (to - a.cuts[i]) * rwidth );
313 }
315 /**Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c);
316 * Further subdivides the Piecewise<T> such that there is a cut at every value in c.
317 * Precondition: c sorted lower to higher.
318 *
319 * //Given Piecewise<T> a and b:
320 * Piecewise<T> ac = a.partition(b.cuts);
321 * Piecewise<T> bc = b.partition(a.cuts);
322 * //ac.cuts should be equivalent to bc.cuts
323 */
324 template<typename T>
325 Piecewise<T> partition(const Piecewise<T> &pw, std::vector<double> const &c) {
326 assert(pw.invariants());
327 if(c.empty()) return Piecewise<T>(pw);
329 Piecewise<T> ret = Piecewise<T>();
330 ret.cuts.reserve(c.size() + pw.cuts.size());
331 ret.segs.reserve(c.size() + pw.cuts.size() - 1);
333 if(pw.empty()) {
334 ret.cuts = c;
335 for(unsigned i = 0; i < c.size() - 1; i++)
336 ret.push_seg(T());
337 return ret;
338 }
340 unsigned si = 0, ci = 0; //Segment index, Cut index
342 //if the cuts have something earlier than the Piecewise<T>, add portions of the first segment
343 while(c[ci] < pw.cuts.front() && ci < c.size()) {
344 bool isLast = (ci == c.size()-1 || c[ci + 1] >= pw.cuts.front());
345 ret.push_cut(c[ci]);
346 ret.push_seg( elem_portion(pw, 0, c[ci], isLast ? pw.cuts.front() : c[ci + 1]) );
347 ci++;
348 }
350 ret.push_cut(pw.cuts[0]);
351 double prev = pw.cuts[0]; //previous cut
352 //Loop which handles cuts within the Piecewise<T> domain
353 //Should have the cuts = segs + 1 invariant
354 while(si < pw.size() && ci <= c.size()) {
355 if(ci == c.size() && prev <= pw.cuts[si]) { //cuts exhausted, straight copy the rest
356 ret.segs.insert(ret.segs.end(), pw.segs.begin() + si, pw.segs.end());
357 ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + si + 1, pw.cuts.end());
358 return ret;
359 }else if(ci == c.size() || c[ci] >= pw.cuts[si + 1]) { //no more cuts within this segment, finalize
360 if(prev > pw.cuts[si]) { //segment already has cuts, so portion is required
361 ret.push_seg(portion(pw[si], pw.segT(prev, si), 1.0));
362 } else { //plain copy is fine
363 ret.push_seg(pw[si]);
364 }
365 ret.push_cut(pw.cuts[si + 1]);
366 prev = pw.cuts[si + 1];
367 si++;
368 } else if(c[ci] == pw.cuts[si]){ //coincident
369 //Already finalized the seg with the code immediately above
370 ci++;
371 } else { //plain old subdivision
372 ret.push(elem_portion(pw, si, prev, c[ci]), c[ci]);
373 prev = c[ci];
374 ci++;
375 }
376 }
378 //input cuts extend further than this Piecewise<T>, extend the last segment.
379 while(ci < c.size()) {
380 if(c[ci] > prev) {
381 ret.push(elem_portion(pw, pw.size() - 1, prev, c[ci]), c[ci]);
382 prev = c[ci];
383 }
384 ci++;
385 }
386 return ret;
387 }
389 /**Piecewise<T> portion(const Piecewise<T> &pw, double from, double to);
390 * Returns a Piecewise<T> with a defined domain of [min(from, to), max(from, to)].
391 */
392 template<typename T>
393 Piecewise<T> portion(const Piecewise<T> &pw, double from, double to) {
394 if(pw.empty() || from == to) return Piecewise<T>();
396 Piecewise<T> ret;
398 double temp = from;
399 from = std::min(from, to);
400 to = std::max(temp, to);
402 unsigned i = pw.segN(from);
403 ret.push_cut(from);
404 if(i == pw.size() - 1 || to <= pw.cuts[i + 1]) { //to/from inhabit the same segment
405 ret.push(elem_portion(pw, i, from, to), to);
406 return ret;
407 }
408 ret.push_seg(portion( pw[i], pw.segT(from, i), 1.0 ));
409 i++;
410 unsigned fi = pw.segN(to, i);
411 if (to == pw.cuts[fi]) fi-=1;
413 ret.segs.insert(ret.segs.end(), pw.segs.begin() + i, pw.segs.begin() + fi); //copy segs
414 ret.cuts.insert(ret.cuts.end(), pw.cuts.begin() + i, pw.cuts.begin() + fi + 1); //and their cuts
416 ret.push_seg( portion(pw[fi], 0.0, pw.segT(to, fi)));
417 if(to != ret.cuts.back()) ret.push_cut(to);
418 ret.invariants();
419 return ret;
420 }
422 template<typename T>
423 Piecewise<T> remove_short_cuts(Piecewise<T> const &f, double tol) {
424 if(f.empty()) return f;
425 Piecewise<T> ret;
426 ret.push_cut(f.cuts[0]);
427 for(unsigned i=0; i<f.size(); i++){
428 if (f.cuts[i+1]-f.cuts[i] >= tol || i==f.size()-1) {
429 ret.push(f[i], f.cuts[i+1]);
430 }
431 }
432 return ret;
433 }
435 template<typename T>
436 Piecewise<T> remove_short_cuts_extending(Piecewise<T> const &f, double tol) {
437 if(f.empty()) return f;
438 Piecewise<T> ret;
439 ret.push_cut(f.cuts[0]);
440 double last = f.cuts[0]; // last cut included
441 for(unsigned i=0; i<f.size(); i++){
442 if (f.cuts[i+1]-f.cuts[i] >= tol) {
443 ret.push(elem_portion(f, i, last, f.cuts[i+1]), f.cuts[i+1]);
444 last = f.cuts[i+1];
445 }
446 }
447 return ret;
448 }
450 template<typename T>
451 std::vector<double> roots(const Piecewise<T> &pw) {
452 std::vector<double> ret;
453 for(unsigned i = 0; i < pw.size(); i++) {
454 std::vector<double> sr = roots(pw[i]);
455 for (unsigned j = 0; j < sr.size(); j++) ret.push_back(sr[j] * (pw.cuts[i + 1] - pw.cuts[i]) + pw.cuts[i]);
457 }
458 return ret;
459 }
461 //IMPL: OffsetableConcept
462 template<typename T>
463 Piecewise<T> operator+(Piecewise<T> const &a, typename T::output_type b) {
464 boost::function_requires<OffsetableConcept<T> >();
465 //TODO:empty
466 Piecewise<T> ret = Piecewise<T>();
467 ret.cuts = a.cuts;
468 for(unsigned i = 0; i < a.size();i++)
469 ret.push_seg(a[i] + b);
470 return ret;
471 }
472 template<typename T>
473 Piecewise<T> operator-(Piecewise<T> const &a, typename T::output_type b) {
474 boost::function_requires<OffsetableConcept<T> >();
475 //TODO: empty
476 Piecewise<T> ret = Piecewise<T>();
477 ret.cuts = a.cuts;
478 for(unsigned i = 0; i < a.size();i++)
479 ret.push_seg(a[i] - b);
480 return ret;
481 }
482 template<typename T>
483 Piecewise<T> operator+=(Piecewise<T>& a, typename T::output_type b) {
484 boost::function_requires<OffsetableConcept<T> >();
486 if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
488 for(unsigned i = 0; i < a.size();i++)
489 a[i] += b;
490 return a;
491 }
492 template<typename T>
493 Piecewise<T> operator-=(Piecewise<T>& a, typename T::output_type b) {
494 boost::function_requires<OffsetableConcept<T> >();
496 if(a.empty()) { a.push_cut(0.); a.push(T(b), 1.); return a; }
498 for(unsigned i = 0;i < a.size();i++)
499 a[i] -= b;
500 return a;
501 }
503 //IMPL: ScalableConcept
504 template<typename T>
505 Piecewise<T> operator-(Piecewise<T> const &a) {
506 boost::function_requires<ScalableConcept<T> >();
508 Piecewise<T> ret;
509 ret.cuts = a.cuts;
510 for(unsigned i = 0; i < a.size();i++)
511 ret.push_seg(- a[i]);
512 return ret;
513 }
514 template<typename T>
515 Piecewise<T> operator*(Piecewise<T> const &a, double b) {
516 boost::function_requires<ScalableConcept<T> >();
518 if(a.empty()) return Piecewise<T>();
520 Piecewise<T> ret;
521 ret.cuts = a.cuts;
522 for(unsigned i = 0; i < a.size();i++)
523 ret.push_seg(a[i] * b);
524 return ret;
525 }
526 template<typename T>
527 Piecewise<T> operator/(Piecewise<T> const &a, double b) {
528 boost::function_requires<ScalableConcept<T> >();
530 //FIXME: b == 0?
531 if(a.empty()) return Piecewise<T>();
533 Piecewise<T> ret;
534 ret.cuts = a.cuts;
535 for(unsigned i = 0; i < a.size();i++)
536 ret.push_seg(a[i] / b);
537 return ret;
538 }
539 template<typename T>
540 Piecewise<T> operator*=(Piecewise<T>& a, double b) {
541 boost::function_requires<ScalableConcept<T> >();
543 if(a.empty()) return Piecewise<T>();
545 for(unsigned i = 0; i < a.size();i++)
546 a[i] *= b;
547 return a;
548 }
549 template<typename T>
550 Piecewise<T> operator/=(Piecewise<T>& a, double b) {
551 boost::function_requires<ScalableConcept<T> >();
553 //FIXME: b == 0?
554 if(a.empty()) return Piecewise<T>();
556 for(unsigned i = 0; i < a.size();i++)
557 a[i] /= b;
558 return a;
559 }
561 //IMPL: AddableConcept
562 template<typename T>
563 Piecewise<T> operator+(Piecewise<T> const &a, Piecewise<T> const &b) {
564 boost::function_requires<AddableConcept<T> >();
566 Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
567 Piecewise<T> ret = Piecewise<T>();
568 assert(pa.size() == pb.size());
569 ret.cuts = pa.cuts;
570 for (unsigned i = 0; i < pa.size(); i++)
571 ret.push_seg(pa[i] + pb[i]);
572 return ret;
573 }
574 template<typename T>
575 Piecewise<T> operator-(Piecewise<T> const &a, Piecewise<T> const &b) {
576 boost::function_requires<AddableConcept<T> >();
578 Piecewise<T> pa = partition(a, b.cuts), pb = partition(b, a.cuts);
579 Piecewise<T> ret = Piecewise<T>();
580 assert(pa.size() == pb.size());
581 ret.cuts = pa.cuts;
582 for (unsigned i = 0; i < pa.size(); i++)
583 ret.push_seg(pa[i] - pb[i]);
584 return ret;
585 }
586 template<typename T>
587 inline Piecewise<T> operator+=(Piecewise<T> &a, Piecewise<T> const &b) {
588 a = a+b;
589 return a;
590 }
591 template<typename T>
592 inline Piecewise<T> operator-=(Piecewise<T> &a, Piecewise<T> const &b) {
593 a = a-b;
594 return a;
595 }
597 template<typename T1,typename T2>
598 Piecewise<T2> operator*(Piecewise<T1> const &a, Piecewise<T2> const &b) {
599 //function_requires<MultiplicableConcept<T1> >();
600 //function_requires<MultiplicableConcept<T2> >();
602 Piecewise<T1> pa = partition(a, b.cuts);
603 Piecewise<T2> pb = partition(b, a.cuts);
604 Piecewise<T2> ret = Piecewise<T2>();
605 assert(pa.size() == pb.size());
606 ret.cuts = pa.cuts;
607 for (unsigned i = 0; i < pa.size(); i++)
608 ret.push_seg(pa[i] * pb[i]);
609 return ret;
610 }
612 template<typename T>
613 inline Piecewise<T> operator*=(Piecewise<T> &a, Piecewise<T> const &b) {
614 a = a * b;
615 return a;
616 }
618 Piecewise<SBasis> divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, unsigned k);
619 //TODO: replace divide(a,b,k) by divide(a,b,tol,k)?
620 //TODO: atm, relative error is ~(tol/a)%. Find a way to make it independant of a.
621 //Nota: the result is 'truncated' where b is smaller than 'zero': ~ a/max(b,zero).
622 Piecewise<SBasis>
623 divide(Piecewise<SBasis> const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
624 Piecewise<SBasis>
625 divide(SBasis const &a, Piecewise<SBasis> const &b, double tol, unsigned k, double zero=1.e-3);
626 Piecewise<SBasis>
627 divide(Piecewise<SBasis> const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
628 Piecewise<SBasis>
629 divide(SBasis const &a, SBasis const &b, double tol, unsigned k, double zero=1.e-3);
631 //Composition: functions called compose_* are pieces of compose that are factored out in pw.cpp.
632 std::map<double,unsigned> compose_pullback(std::vector<double> const &cuts, SBasis const &g);
633 int compose_findSegIdx(std::map<double,unsigned>::iterator const &cut,
634 std::map<double,unsigned>::iterator const &next,
635 std::vector<double> const &levels,
636 SBasis const &g);
638 //TODO: add concept check
639 template<typename T>
640 Piecewise<T> compose(Piecewise<T> const &f, SBasis const &g){
641 Piecewise<T> result;
642 if (f.empty()) return result;
643 if (g.isZero()) return Piecewise<T>(f(0));
644 if (f.size()==1){
645 double t0 = f.cuts[0], width = f.cuts[1] - t0;
646 return (Piecewise<T>) compose(f.segs[0],compose(Linear(-t0 / width, (1-t0) / width), g));
647 }
649 //first check bounds...
650 Interval bs = *bounds_fast(g);
651 if (f.cuts.front() > bs.max() || bs.min() > f.cuts.back()){
652 int idx = (bs.max() < f.cuts[1]) ? 0 : f.cuts.size()-2;
653 double t0 = f.cuts[idx], width = f.cuts[idx+1] - t0;
654 return (Piecewise<T>) compose(f.segs[idx],compose(Linear(-t0 / width, (1-t0) / width), g));
655 }
657 std::vector<double> levels;//we can forget first and last cuts...
658 levels.insert(levels.begin(),f.cuts.begin()+1,f.cuts.end()-1);
659 //TODO: use a std::vector<pairs<double,unsigned> > instead of a map<double,unsigned>.
660 std::map<double,unsigned> cuts_pb = compose_pullback(levels,g);
662 //-- Compose each piece of g with the relevant seg of f.
663 result.cuts.push_back(0.);
664 std::map<double,unsigned>::iterator cut=cuts_pb.begin();
665 std::map<double,unsigned>::iterator next=cut; next++;
666 while(next!=cuts_pb.end()){
667 //assert(std::abs(int((*cut).second-(*next).second))<1);
668 //TODO: find a way to recover from this error? the root finder missed some root;
669 // the levels/variations of f might be too close/fast...
670 int idx = compose_findSegIdx(cut,next,levels,g);
671 double t0=(*cut).first;
672 double t1=(*next).first;
674 SBasis sub_g=compose(g, Linear(t0,t1));
675 sub_g=compose(Linear(-f.cuts[idx]/(f.cuts[idx+1]-f.cuts[idx]),
676 (1-f.cuts[idx])/(f.cuts[idx+1]-f.cuts[idx])),sub_g);
677 result.push(compose(f[idx],sub_g),t1);
678 cut++;
679 next++;
680 }
681 return(result);
682 }
684 //TODO: add concept check for following composition functions
685 template<typename T>
686 Piecewise<T> compose(Piecewise<T> const &f, Piecewise<SBasis> const &g){
687 Piecewise<T> result;
688 for(unsigned i = 0; i < g.segs.size(); i++){
689 Piecewise<T> fgi=compose(f, g.segs[i]);
690 fgi.setDomain(Interval(g.cuts[i], g.cuts[i+1]));
691 result.concat(fgi);
692 }
693 return result;
694 }
696 /*
697 Piecewise<D2<SBasis> > compose(D2<SBasis2d> const &sb2d, Piecewise<D2<SBasis> > const &pwd2sb){
698 Piecewise<D2<SBasis> > result;
699 result.push_cut(0.);
700 for(unsigned i = 0; i < pwd2sb.size(); i++){
701 result.push(compose_each(sb2d,pwd2sb[i]),i+1);
702 }
703 return result;
704 }*/
706 template <typename T>
707 Piecewise<T> Piecewise<T>::operator()(SBasis f){return compose((*this),f);}
708 template <typename T>
709 Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){return compose((*this),f);}
711 template<typename T>
712 Piecewise<T> integral(Piecewise<T> const &a) {
713 Piecewise<T> result;
714 result.segs.resize(a.segs.size());
715 result.cuts = a.cuts;
716 typename T::output_type c = a.segs[0].at0();
717 for(unsigned i = 0; i < a.segs.size(); i++){
718 result.segs[i] = integral(a.segs[i])*(a.cuts[i+1]-a.cuts[i]);
719 result.segs[i]+= c-result.segs[i].at0();
720 c = result.segs[i].at1();
721 }
722 return result;
723 }
725 template<typename T>
726 Piecewise<T> derivative(Piecewise<T> const &a) {
727 Piecewise<T> result;
728 result.segs.resize(a.segs.size());
729 result.cuts = a.cuts;
730 for(unsigned i = 0; i < a.segs.size(); i++){
731 result.segs[i] = derivative(a.segs[i])/(a.cuts[i+1]-a.cuts[i]);
732 }
733 return result;
734 }
736 std::vector<double> roots(Piecewise<SBasis> const &f);
738 std::vector<std::vector<double> >multi_roots(Piecewise<SBasis> const &f, std::vector<double> const &values);
740 template<typename T>
741 Piecewise<T> reverse(Piecewise<T> const &f) {
742 Piecewise<T> ret = Piecewise<T>();
743 ret.cuts.resize(f.cuts.size());
744 ret.segs.resize(f.segs.size());
745 double start = f.cuts[0];
746 double end = f.cuts.back();
747 for (unsigned i = 0; i < f.cuts.size(); i++) {
748 double x = f.cuts[f.cuts.size() - 1 - i];
749 ret.cuts[i] = end - (x - start);
750 }
751 for (unsigned i = 0; i < f.segs.size(); i++)
752 ret.segs[i] = reverse(f[f.segs.size() - i - 1]);
753 return ret;
754 }
757 }
759 #endif //SEEN_GEOM_PW_SB_H
760 /*
761 Local Variables:
762 mode:c++
763 c-file-style:"stroustrup"
764 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
765 indent-tabs-mode:nil
766 fill-column:99
767 End:
768 */
769 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :