Subversion Repositories Scribus

Rev

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 :