sbasis-roots.cpp revision 8001ba81cb851b38d86650a2fef5817facffb763
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh/** root finding for sbasis functions.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * Copyright 2006 N Hurst
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * Copyright 2007 JF Barraud
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * It is more efficient to find roots of f(t) = c_0, c_1, ... all at once, rather than iterating.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * multi-roots using bernstein method, one approach would be:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh take median and find roots of that
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh whenever a segment lies entirely on one side of the median,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh find the median of the half and recurse.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh in essence we are implementing quicksort on a continuous function
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * the gsl poly roots finder is faster than bernstein too, but we don't use it for 3 reasons:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh a) it requires convertion to poly, which is numerically unstable
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh b) it requires gsl (which is currently not a dependency, and would bring in a whole slew of unrelated stuff)
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh c) it finds all roots, even complex ones. We don't want to accidently treat a nearly real root as a real root
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshFrom memory gsl poly roots was about 10 times faster than bernstein in the case where all the roots
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshare in [0,1] for polys of order 5. I spent some time working out whether eigenvalue root finding
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshcould be done directly in sbasis space, but the maths was too hard for me. -- njh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshjfbarraud: eigenvalue root finding could be done directly in sbasis space ?
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshnjh: I don't know, I think it should. You would make a matrix whose characteristic polynomial was
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshcorrect, but do it by putting the sbasis terms in the right spots in the matrix. normal eigenvalue
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshroot finding makes a matrix that is a diagonal + a row along the top. This matrix has the property
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshthat its characteristic poly is just the poly whose coefficients are along the top row.
3711b3e25395437ee0a09dbbb2a76d999c4ef322mikloshNow an sbasis function is a linear combination of the poly coeffs. So it seems to me that you
3711b3e25395437ee0a09dbbb2a76d999c4ef322mikloshshould be able to put the sbasis coeffs directly into a matrix in the right spots so that the
fba63a357654d8b3e84c60007e40aa698cd45d19mikloshcharacteristic poly is the sbasis. You'll still have problems b) and c).
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshWe might be able to lift an eigenvalue solver and include that directly into 2geom. Eigenvalues
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshalso allow you to find intersections of multiple curves but require solving n*m x n*m matrices.
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshusing namespace std;
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double a=sb[j][0];
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double v, t = 0;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (v>=0 || t<0 || t>1) {
dc4f69a188c203f2fdc65f22d0d57904a8c52dd7miklosh if (v<=0 || t<0 || t>1) {
68664e00e2372534b4df2fdc5f54f836bafece18mikloshInterval bounds_local(const SBasis &sb, const Interval &i, int order) {
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh double a=sb[j][0];
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double t = 0;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
1667116521643e2475184b048e0abb77a2aa9735miklosh hi = std::max(a*(1-t0)+b*t0+hi*t0*(1-t0),a*(1-t1)+b*t1+hi*t1*(1-t1));
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh//-- multi_roots ------------------------------------
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh// goal: solve f(t)=c for several c at once.
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh/* algo: -compute f at both ends of the given segment [a,b].
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh -compute bounds m<df(t)<M for df on the segment.
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh let c and C be the levels below and above f(a):
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh going from f(a) down to c with slope m takes at least time (f(a)-c)/m
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh going from f(a) up to C with slope M takes at least time (C-f(a))/M
68664e00e2372534b4df2fdc5f54f836bafece18miklosh From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
68664e00e2372534b4df2fdc5f54f836bafece18miklosh Do the same for b: compute some b' such that there are no roots in (b',b].
68664e00e2372534b4df2fdc5f54f836bafece18miklosh -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
68664e00e2372534b4df2fdc5f54f836bafece18miklosh unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
68664e00e2372534b4df2fdc5f54f836bafece18miklosh making things tricky and unpleasant...
68664e00e2372534b4df2fdc5f54f836bafece18miklosh//TODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
68664e00e2372534b4df2fdc5f54f836bafece18mikloshstatic int upper_level(vector<double> const &levels,double x,double tol=0.){
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh return(upper_bound(levels.begin(),levels.end(),x-tol)-levels.begin());
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (f.size()==0){
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (idx<(int)levels.size()&&fabs(levels.at(idx))<=vtol){
1667116521643e2475184b048e0abb77a2aa9735miklosh////usefull?
1667116521643e2475184b048e0abb77a2aa9735miklosh// if (f.size()==1){
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh// int idxa=upper_level(levels,fa);
1667116521643e2475184b048e0abb77a2aa9735miklosh// int idxb=upper_level(levels,fb);
1667116521643e2475184b048e0abb77a2aa9735miklosh// if (fa==fb){
1667116521643e2475184b048e0abb77a2aa9735miklosh// if (fa==levels[idxa]){
1667116521643e2475184b048e0abb77a2aa9735miklosh// roots[a]=idxa;
1667116521643e2475184b048e0abb77a2aa9735miklosh// roots[b]=idxa;
1667116521643e2475184b048e0abb77a2aa9735miklosh// int idx_min=std::min(idxa,idxb);
1667116521643e2475184b048e0abb77a2aa9735miklosh// int idx_max=std::max(idxa,idxb);
1667116521643e2475184b048e0abb77a2aa9735miklosh// if (idx_max==levels.size()) idx_max-=1;
1667116521643e2475184b048e0abb77a2aa9735miklosh// for(int i=idx_min;i<=idx_max; i++){
1667116521643e2475184b048e0abb77a2aa9735miklosh// double t=a+(b-a)*(levels[i]-fa)/(fb-fa);
1667116521643e2475184b048e0abb77a2aa9735miklosh// if(a<t&&t<b) roots[t]=i;
1667116521643e2475184b048e0abb77a2aa9735miklosh if ((b-a)<htol){
1667116521643e2475184b048e0abb77a2aa9735miklosh //TODO: use different tol for t and f ?
1667116521643e2475184b048e0abb77a2aa9735miklosh //TODO: unsigned idx ? (remove int casts when fixed)
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh int idx=std::min(upper_level(levels,fa,vtol),upper_level(levels,fb,vtol));
1cda9431ef400135f5e1bd899a94b921bdad0eafmiklosh if((fa-c)*(fb-c)<=0||fabs(fa-c)<vtol||fabs(fb-c)<vtol){
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh //first times when a level (higher or lower) can be reached from a or b.
1667116521643e2475184b048e0abb77a2aa9735miklosh if (idxa<(int)levels.size() && fabs(fa-levels.at(idxa))<vtol){//a can be considered a root.
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh //ta_hi=ta_lo=a;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (idxb<(int)levels.size() && fabs(fb-levels.at(idxb))<vtol){//b can be considered a root.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //tb_hi=tb_lo=b;
d27f5758e12c3107ee69e66702043931e0756f6bmiklosh //hum, rounding errors frighten me! so I add this +tol...
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh multi_roots_internal(f,df,levels,roots,htol,vtol,t0,f(t0),t1,f(t1));
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (idx<(int)levels.size() && fabs(ft-levels.at(idx))<vtol){//t can be considered a root.
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh //we do not want to count it twice (from the left and from the right)
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh multi_roots_internal(f,df,levels,roots,htol,vtol,t0 ,f(t0) ,t_left,ft_left);
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh multi_roots_internal(f,df,levels,roots,htol,vtol,t_right,ft_right,t1 ,f(t1) );
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshstd::vector<std::vector<double> > multi_roots(SBasis const &f,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh//-------------------------------------
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double n = p.size();
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh while(a > tol) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "xk = " << xk << std::endl;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh b = xk*b + p[j];
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh err *= 1e-7; // magic epsilon for convergence, should be computed from tol
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double px = b;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //if(std::norm(px) < tol*tol)
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh // return xk;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double G = d / px;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double H = G*G - f / px;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "G = " << G << "H = " << H;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //assert(radicand.real() > 0);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "radicand = " << radicand << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "a = " << a << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh a = n / (a + G);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "a = " << a << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "xk = " << xk << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return; // no roots here
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double t = s[0][0] / (s[0][0] - s[0][1]);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh subdiv_sbasis(compose(s, Linear(0, 0.5)), roots, left, middle);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh subdiv_sbasis(compose(s, Linear(0.5, 1.)), roots, middle, right);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh// It is faster to use the bernstein root finder for small degree polynomials (<100?.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double d = s[0][0] - s[0][1];
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if(d != 0) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double r = s[0][0] / d;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh switch(s.size()) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh Local Variables:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh c-file-style:"stroustrup"
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh indent-tabs-mode:nil
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh fill-column:99
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :