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