Rev 18123 | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
12065 | fschmid | 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 |
||
18123 | mrdocs | 16 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
12065 | fschmid | 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 | |||
83 | /**Convenience/implementation hiding function to add segment/cut pairs. |
||
84 | * Asserts that basic size and order invariants are correct |
||
85 | */ |
||
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 | |||
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 | |||
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 | |||
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> |
||
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 |
||
22493 | jghali | 294 | while(ci < c.size() && c[ci] < pw.cuts.front()) { |
12065 | fschmid | 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> |
||
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 : |