sbasis-roots.cpp revision 8001ba81cb851b38d86650a2fef5817facffb763
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen/** root finding for sbasis functions.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Copyright 2006 N Hurst
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Copyright 2007 JF Barraud
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
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen 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;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenInterval 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) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenInterval bounds_local(const SBasis &sb, const Interval &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].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen -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 (f.size()==0){
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) );
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstd::vector<std::vector<double> > multi_roots(SBasis const &f,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen//-------------------------------------
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double n = p.size();
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen while(a > tol) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //std::cout << "xk = " << xk << std::endl;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen b = xk*b + p[j];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen err *= 1e-7; // magic epsilon for convergence, should be computed from tol
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //if(std::norm(px) < tol*tol)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen // return xk;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double G = d / px;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double H = G*G - f / px;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //std::cout << "G = " << G << "H = " << H;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //assert(radicand.real() > 0);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //std::cout << "radicand = " << radicand << std::endl;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //std::cout << "a = " << a << std::endl;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen a = n / (a + G);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //std::cout << "a = " << a << std::endl;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //std::cout << "xk = " << xk << std::endl;
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;
6492fcbc5fbe9a07962f989247731b6435fd72b9johanengelen 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
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :