76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * @brief Root finding for sbasis functions.
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * Nathan Hurst <njh@njhurst.com>
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * Copyright 2006-2007 Authors
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * This library is free software; you can redistribute it and/or
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * modify it either under the terms of the GNU Lesser General Public
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * License version 2.1 as published by the Free Software Foundation
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * (the "LGPL") or, at your option, under the terms of the Mozilla
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * Public License Version 1.1 (the "MPL"). If you do not alter this
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * notice, a recipient may use your version of this file under either
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * the MPL or the LGPL.
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * You should have received a copy of the LGPL along with this library
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * in the file COPYING-LGPL-2.1; if not, write to the Free Software
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * You should have received a copy of the MPL along with this library
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * in the file COPYING-MPL-1.1
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * The contents of this file are subject to the Mozilla Public License
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * Version 1.1 (the "License"); you may not use this file except in
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * compliance with the License. You may obtain a copy of the License at
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * OF ANY KIND, either express or implied. See the LGPL or the MPL for
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński * the specific language governing rights and limitations.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * It is more efficient to find roots of f(t) = c_0, c_1, ... all at once, rather than iterating.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * multi-roots using bernstein method, one approach would be:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen take median and find roots of that
6bc0b25077dcb0cce5dea357de5bab735babe891johanengelen whenever a segment lies entirely on one side of the median,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen find the median of the half and recurse.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen in essence we are implementing quicksort on a continuous function
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * the gsl poly roots finder is faster than bernstein too, but we don't use it for 3 reasons:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen a) it requires convertion to poly, which is numerically unstable
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen b) it requires gsl (which is currently not a dependency, and would bring in a whole slew of unrelated stuff)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen c) it finds all roots, even complex ones. We don't want to accidently treat a nearly real root as a real root
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenFrom memory gsl poly roots was about 10 times faster than bernstein in the case where all the roots
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenare in [0,1] for polys of order 5. I spent some time working out whether eigenvalue root finding
981b809bc6ed10a21e89444d9447e5475801874fjohanengelencould be done directly in sbasis space, but the maths was too hard for me. -- njh
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenjfbarraud: eigenvalue root finding could be done directly in sbasis space ?
981b809bc6ed10a21e89444d9447e5475801874fjohanengelennjh: I don't know, I think it should. You would make a matrix whose characteristic polynomial was
981b809bc6ed10a21e89444d9447e5475801874fjohanengelencorrect, but do it by putting the sbasis terms in the right spots in the matrix. normal eigenvalue
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenroot finding makes a matrix that is a diagonal + a row along the top. This matrix has the property
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenthat its characteristic poly is just the poly whose coefficients are along the top row.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenNow an sbasis function is a linear combination of the poly coeffs. So it seems to me that you
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenshould be able to put the sbasis coeffs directly into a matrix in the right spots so that the
981b809bc6ed10a21e89444d9447e5475801874fjohanengelencharacteristic poly is the sbasis. You'll still have problems b) and c).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenWe might be able to lift an eigenvalue solver and include that directly into 2geom. Eigenvalues
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenalso allow you to find intersections of multiple curves but require solving n*m x n*m matrices.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenusing namespace std;
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** Find the smallest interval that bounds a
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param a sbasis function
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \returns inteval
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed GouldOptInterval bounds_exact(SBasisOf<double> const &a) {
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** Find a small interval that bounds a
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param a sbasis function
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \returns inteval
1f37e9d97c3bb8cf02b2cc80af8dcfc9aba7c7b4johanengelen// I have no idea how this works, some clever bounding argument by jfb.
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed GouldOptInterval bounds_fast(const SBasisOf<double> &sb, int order) {
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed GouldOptInterval bounds_fast(const SBasis &sb, int order) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double a=sb[j][0];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen double v, t = 0;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (v>=0 || t<0 || t>1) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (v<=0 || t<0 || t>1) {
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** Find a small interval that bounds a(t) for t in i to order order
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param sb sbasis function
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param i domain interval
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param order number of terms
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould \return interval
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed GouldOptInterval bounds_local(const SBasisOf<double> &sb, const OptInterval &i, int order) {
6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed GouldOptInterval bounds_local(const SBasis &sb, const OptInterval &i, int order) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double a=sb[j][0];
29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen double t = 0;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen hi = std::max(a*(1-t0)+b*t0+hi*t0*(1-t0),a*(1-t1)+b*t1+hi*t1*(1-t1));
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen//-- multi_roots ------------------------------------
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// goal: solve f(t)=c for several c at once.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen/* algo: -compute f at both ends of the given segment [a,b].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen -compute bounds m<df(t)<M for df on the segment.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen let c and C be the levels below and above f(a):
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen going from f(a) down to c with slope m takes at least time (f(a)-c)/m
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen going from f(a) up to C with slope M takes at least time (C-f(a))/M
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Do the same for b: compute some b' such that there are no roots in (b',b].
6bc0b25077dcb0cce5dea357de5bab735babe891johanengelen -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen making things tricky and unpleasant...
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen//TODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic int upper_level(vector<double> const &levels,double x,double tol=0.){
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return(upper_bound(levels.begin(),levels.end(),x-tol)-levels.begin());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic void multi_roots_internal(SBasis const &f,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (idx<(int)levels.size()&&fabs(levels.at(idx))<=vtol){
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// if (f.size()==1){
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// int idxa=upper_level(levels,fa);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// int idxb=upper_level(levels,fb);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// if (fa==fb){
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// if (fa==levels[idxa]){
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// roots[a]=idxa;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// roots[b]=idxa;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// int idx_min=std::min(idxa,idxb);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// int idx_max=std::max(idxa,idxb);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// if (idx_max==levels.size()) idx_max-=1;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// for(int i=idx_min;i<=idx_max; i++){
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// double t=a+(b-a)*(levels[i]-fa)/(fb-fa);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// if(a<t&&t<b) roots[t]=i;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //TODO: use different tol for t and f ?
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //TODO: unsigned idx ? (remove int casts when fixed)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen int idx=std::min(upper_level(levels,fa,vtol),upper_level(levels,fb,vtol));
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if((fa-c)*(fb-c)<=0||fabs(fa-c)<vtol||fabs(fb-c)<vtol){
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //first times when a level (higher or lower) can be reached from a or b.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen ta_hi=ta_lo=b+1;//default values => no root there.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen tb_hi=tb_lo=a-1;//default values => no root there.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (idxa<(int)levels.size() && fabs(fa-levels.at(idxa))<vtol){//a can be considered a root.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //ta_hi=ta_lo=a;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (idxb<(int)levels.size() && fabs(fb-levels.at(idxb))<vtol){//b can be considered a root.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //tb_hi=tb_lo=b;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //hum, rounding errors frighten me! so I add this +tol...
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen multi_roots_internal(f,df,levels,roots,htol,vtol,t0,f(t0),t1,f(t1));
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (idx<(int)levels.size() && fabs(ft-levels.at(idx))<vtol){//t can be considered a root.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //we do not want to count it twice (from the left and from the right)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen multi_roots_internal(f,df,levels,roots,htol,vtol,t0 ,f(t0) ,t_left,ft_left);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen multi_roots_internal(f,df,levels,roots,htol,vtol,t_right,ft_right,t1 ,f(t1) );
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** Solve f(t)=c for several c at once.
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param f sbasis function
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param levels vector of 'y' values
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param htol, vtol
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param a, b left and right bounds
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \returns a vector of vectors, one for each y giving roots
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed GouldEffectively computes:
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gouldresults = roots(f(y_i)) for all y_i
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould* algo: -compute f at both ends of the given segment [a,b].
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould -compute bounds m<df(t)<M for df on the segment.
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould let c and C be the levels below and above f(a):
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould going from f(a) down to c with slope m takes at least time (f(a)-c)/m
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould going from f(a) up to C with slope M takes at least time (C-f(a))/M
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould Do the same for b: compute some b' such that there are no roots in (b',b].
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould making things tricky and unpleasant...
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed GouldTODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstd::vector<std::vector<double> > multi_roots(SBasis const &f,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
6bc0b25077dcb0cce5dea357de5bab735babe891johanengelen multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic bool compareIntervalMin( Interval I, Interval J ){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic bool compareIntervalMax( Interval I, Interval J ){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen//find the first interval whose max is >= x
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic unsigned upper_level(vector<Interval> const &levels, double x ){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen return( lower_bound( levels.begin(), levels.end(), Interval(x,x), compareIntervalMax) - levels.begin() );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic std::vector<Interval> fuseContiguous(std::vector<Interval> const &sets, double tol=0.){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen if ( result.back().max() + tol >= sets[i].min() ){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen/** level_sets internal method.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen* algorithm: (~adaptation of Newton method versus 'accroissements finis')
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen -compute f at both ends of the given segment [a,b].
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen -compute bounds m<df(t)<M for df on the segment.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen Suppose f(a) is between two 'levels' c and C. Then
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen f wont enter c before a + (f(a)-c.max())/m
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen f wont enter C before a + (C.min()-f(a))/M
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen From this we conclude nothing happens before a'=a+min((f(a)-c.max())/m,(C.min()-f(a))/M).
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen We do the same for b: compute some b' such that nothing happens in (b',b].
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen If f(a) or f(b) belongs to some 'level' C, then use the same argument to find a' or b' such
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen that f remains in C on [a,a'] or [b',b]. In case f is monotonic, we also know f won't enter another
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen level before or after some time, allowing us to restrict the search a little more.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen making things tricky and unpleasant...
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic void level_sets_internal(SBasis const &f,
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen if (idx<levels.size() && levels[idx].contains(0.)){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen Interval bs = *bounds_local(df,Interval(a,b));
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //first times when a level (higher or lower) can be reached from a or b.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen double ta_hi; // f remains below next level for t<ta_hi
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen double ta_lo; // f remains above prev level for t<ta_lo
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen double tb_hi; // f remains below next level for t>tb_hi
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen double tb_lo; // f remains above next level for t>tb_lo
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen ta_hi=ta_lo=b+1;//default values => no root there.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen tb_hi=tb_lo=a-1;//default values => no root there.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //--- if f(a) belongs to a level.-------
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen if ( idxa < levels.size() && levels[idxa].contains( fa ) ){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //find the first time when we may exit this level.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen ta_lo = a + ( levels[idxa].min() - fa)/bs.min();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen ta_hi = a + ( levels[idxa].max() - fa)/bs.max();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //move to that time for the next iteration.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen solsets[idxa].push_back( Interval( a, std::min( ta_lo, ta_hi ) ) );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //--- if f(b) does not belong to a level.-------
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen ta_lo = a + ( levels[idxa-1].max() - fa)/bs.min();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen ta_hi = a + ( levels[idxa].min() - fa)/bs.max();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //--- if f(b) belongs to a level.-------
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen if (idxb<levels.size() && levels.at(idxb).contains(fb)){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //find the first time from b when we may exit this level.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen tb_lo = b + ( levels[idxb].min() - fb ) / bs.max();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen tb_hi = b + ( levels[idxb].max() - fb ) / bs.min();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //move to that time for the next iteration.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen solsets[idxb].push_back( Interval( std::max( tb_lo, tb_hi ), b) );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //--- if f(b) does not belong to a level.-------
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen tb_lo = b + ( levels[idxb-1].max() - fb)/bs.max();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen tb_hi = b + ( levels[idxb].min() - fb)/bs.min();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen tb_hi = b + ( levels[idxb ].min() - fb ) / bs.min();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen tb_lo = b + ( levels[idxb-1].max() - fb ) / bs.max();
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //let [t0,t1] be the next interval where to search.
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //if the interval is smaller than our resolution:
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //pretend f simultaneously meets all the levels between f(t0) and f(t1)...
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //push [t0,t1] into all crossed level. Cheat to avoid overlapping intervals on different levels?
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen for (unsigned idx = idxmin; idx < idxmax; idx++){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen if ( idxmax < levels.size() && f_t0t1.intersects( levels[idxmax] ) ){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen solsets[idxmax].push_back( Interval( t0, t1 ) );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen //To make sure we finally exit the level jump at least by tol:
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen level_sets_internal( f, df, levels, solsets, t0, f(t0), t, ft );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen level_sets_internal( f, df, levels, solsets, t, ft, t1, f(t1) );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstd::vector<std::vector<Interval> > level_sets(SBasis const &f,
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen double a, double b, double tol){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen std::vector<std::vector<Interval> > solsets(levels.size(), std::vector<Interval>());
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen level_sets_internal(f,df,levels,solsets,a,f(a),b,f(b),tol);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen // Fuse overlapping intervals...
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen std::sort( solsets[i].begin(), solsets[i].end(), compareIntervalMin );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen solsets[i] = fuseContiguous( solsets[i], tol );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstd::vector<Interval> level_set (SBasis const &f, double level, double vtol, double a, double b, double tol){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen Interval fat_level( level - vtol, level + vtol );
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstd::vector<Interval> level_set (SBasis const &f, Interval const &level, double a, double b, double tol){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen return level_sets(f,levels, a, b, tol).front() ;
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstd::vector<std::vector<Interval> > level_sets (SBasis const &f, std::vector<double> const &levels, double vtol, double a, double b, double tol){
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen std::vector<Interval> fat_levels( levels.size(), Interval());
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen fat_levels[i] = Interval( levels[i]-vtol, levels[i]+vtol);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen//-------------------------------------
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen//-------------------------------------
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return; // no roots here
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double t = s[0][0] / (s[0][0] - s[0][1]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen subdiv_sbasis(compose(s, Linear(0, 0.5)), roots, left, middle);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen subdiv_sbasis(compose(s, Linear(0.5, 1.)), roots, middle, right);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// It is faster to use the bernstein root finder for small degree polynomials (<100?.
6492fcbc5fbe9a07962f989247731b6435fd72b9johanengelen double d = s[0][0] - s[0][1];
6492fcbc5fbe9a07962f989247731b6435fd72b9johanengelen if(d != 0) {
6492fcbc5fbe9a07962f989247731b6435fd72b9johanengelen double r = s[0][0] / d;
98a704c566ce5a750d76d5fc9675ccc804ac65f5johanengelen if(0 <= r && r <= 1)
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstd::vector<double> roots1(SBasis const & s, Interval const ivl) {
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen double d = s[0][0] - s[0][1];
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen double r = s[0][0] / d;
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould/** Find all t s.t s(t) = 0
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \param a sbasis function
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński \see Bezier::roots
01d27eab5fca2dcb8e883011f8be77ae6b78a11cTed Gould \returns vector of zeros (roots)
6492fcbc5fbe9a07962f989247731b6435fd72b9johanengelen switch(s.size()) {
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstd::vector<double> roots(SBasis const & s, Interval const ivl) {
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen switch(s.size()) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Local Variables:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen c-file-style:"stroustrup"
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen indent-tabs-mode:nil
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen fill-column:99
a4030d5ca449e7e384bc699cd249ee704faaeab0Chris Morgan// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :