Scribus
Open source desktop publishing at your fingertips
piecewise.h
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 
53  Piecewise() {}
54 
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 
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]; }
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  //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 
80  inline unsigned size() const { return segs.size(); }
81  inline bool empty() const { return segs.empty(); }
82 
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.
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.
97  inline void push_seg(const T &s) { segs.push_back(s); }
98 
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 
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 
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
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.
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
154  inline Interval domain() const { return Interval(cuts.front(), cuts.back()); }
155 
156  //Transforms the domain into another interval
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.
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.
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.
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>
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>
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>
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>
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 
275 template<typename T>
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(c[ci] < pw.cuts.front() && ci < c.size()) {
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 
343 template<typename T>
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>
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>
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>
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>
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>
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>
647 Piecewise<T> Piecewise<T>::operator()(SBasis f){return compose((*this),f);}
648 template <typename T>
649 Piecewise<T> Piecewise<T>::operator()(Piecewise<SBasis>f){return compose((*this),f);}
650 
651 template<typename T>
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>
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 :
Definition: angle.h:38
Definition: concepts.h:59
void push(const T &s, double to)
Definition: piecewise.h:86
Piecewise< T > partition(const Piecewise< T > &pw, std::vector< double > const &c)
Definition: piecewise.h:276
Definition: piecewise.h:45
unsigned segN(double t, int low=0, int high=-1) const
Definition: piecewise.h:102
Definition: sbasis.h:48
double segT(double t, int i=-1) const
Definition: piecewise.h:124