Ticket #3813: sage-adaptive-plot.patch

File sage-adaptive-plot.patch, 6.7 kB (added by anakha, 5 months ago)

Better adaptive rendering

  • a/sage/plot/plot.py

    old new  
    34163416    The plot options are 
    34173417        plot_points -- the number of points to initially plot before 
    34183418                       doing adaptive refinement 
    3419         plot_division -- the maximum number of subdivisions to 
    3420                          introduce in adaptive refinement. 
    3421         max_bend      -- parameter that affects adaptive refinement 
     3419        adaptive_recursion -- how many levels of recursion to go before 
     3420                              giving up when doing adaptive refinement 
     3421                              Setting this to 0 disables adaptive refinement 
     3422        adaptive_tolerance -- how large a difference should be before the 
     3423                              adaptive refinement code considers it significant. 
     3424                              This depends on the interval you use by default. 
    34223425        xmin -- starting x value 
    34233426        xmax -- ending x value 
    34243427        color -- an rgb-tuple (r,g,b) with each of r,g,b between 0 and 1, or 
    …… 
    34613464        Graphics object consisting of 1 graphics primitive 
    34623465        sage: len(P)     # number of graphics primitives 
    34633466        1 
    3464         sage: len(P[0])  # how many points were computed 
    3465         200 
     3467        sage: len(P[0])  # how many points were computed (random) 
     3468        224 
    34663469        sage: P          # render 
    34673470         
    34683471        sage: P = plot(sin, (0,10), plot_points=10); print P 
    34693472        Graphics object consisting of 1 graphics primitive 
    34703473        sage: len(P[0])  # random output 
    3471         80 
     3474        32 
    34723475        sage: P          # render 
    34733476 
    34743477    We plot with randomize=False, which makes the initial sample  
    34753478    points evenly spaced (hence always the same).  Adaptive plotting  
    3476     might insert other points, however, unless plot_division=0.  
    3477         sage: p=plot(lambda x: 1, (x,0,3), plot_points=4, randomize=False, plot_division=0) 
     3479    might insert other points, however, unless adaptive_recursion=0.  
     3480        sage: p=plot(lambda x: 1, (x,0,3), plot_points=4, randomize=False, adaptive_recursion=0) 
    34783481        sage: list(p[0]) 
    34793482        [(0.0, 1.0), (1.0, 1.0), (2.0, 1.0), (3.0, 1.0)] 
    34803483 
    …… 
    34933496    and plots extra points near the origin. 
    34943497        sage: plot(sin(1/x), (x, -1, 1)) 
    34953498 
    3496     With the \code{plot_points} option you can increase the number 
    3497     of sample points, to obtain a more accurate plot.  
    3498         sage: plot(sin(1/x), (x, -1, 1), plot_points=1000) 
    3499  
    35003499    Note that the independent variable may be omitted if there is no 
    35013500    ambiguity: 
    3502         sage: plot(sin(1/x), (-1, 1), plot_points=1000) 
     3501        sage: plot(sin(1/x), (-1, 1)) 
     3502 
     3503    The algorithm used to insert extra points is actually pretty simple. On 
     3504    the picture drawn by the lines below: 
     3505        sage: p = plot(x^2, (0, 1)) + line([(0,0), (1,1)], rgbcolor='green') 
     3506        sage: p += line([(0.5, 0.5), (0.5, 0.5^2)], rgbcolor='purple') 
     3507        sage: p += point(((0, 0), (0.5, 0.5), (0.5, 0.5^2), (1, 1)), rgbcolor='red', pointsize=20) 
     3508        sage: p += text('A', (-0.05, 0.1), rgbcolor='red') 
     3509        sage: p += text('B', (1.01, 1.1), rgbcolor='red') 
     3510        sage: p += text('C', (0.48, 0.57), rgbcolor='red') 
     3511        sage: p += text('D', (0.53, 0.18), rgbcolor='red') 
     3512        sage: p 
     3513     
     3514    You have the function (in blue) its approximation (in green) passing by the 
     3515    points A and B.  The algorithm find the middle point in x of A and B, which 
     3516    is C.  It then computes the distance in y between C and D  and adds D  
     3517    to the curve points if it exceeds some treshold. 
     3518     
     3519    If D was added, it recursivly applies the same algorithm between A and D, 
     3520    and D and B.  It does this up to 10 levels deep by default. 
     3521 
    35033522 
    35043523    The actual sample points are slightly randomized, so the above 
    35053524    plots may look slightly different each time you draw them. 
    …… 
    35433562    def _reset(self): 
    35443563        o = self.options 
    35453564        o['plot_points'] = 200 
    3546         o['plot_division'] = 1000  
    3547         o['max_bend'] = 0.1        
     3565        o['adaptive_recursion'] = 10       
    35483566        o['rgbcolor'] = (0,0,1)    
    35493567 
    35503568    def __repr__(self): 
    3551         """ 
     3569        r""" 
    35523570        Returns a string representation of this PlotFactory object. 
    35533571 
    35543572        TESTS: 
    …… 
    36543672             
    36553673        # adaptive refinement 
    36563674        i, j = 0, 0 
    3657         max_bend = float(options['max_bend']) 
    3658         del options['max_bend'] 
    3659         plot_division = int(options['plot_division']) 
    3660         del options['plot_division'] 
     3675        if 'adaptive_tolerance' in options: 
     3676            adaptive_tolerance = float(options['adaptive_tolerance']) 
     3677            del options['adaptive_tolerance'] 
     3678        else: 
     3679            adaptive_tolerance = delta * 0.01 
     3680        adaptive_recursion = int(options['adaptive_recursion']) 
     3681        del options['adaptive_recursion'] 
     3682         
     3683        def refine(f, p1, p2, level=0): 
     3684            if level >= adaptive_recursion: 
     3685                return [] 
     3686            x = float((p1[0] + p2[0])/2) 
     3687            try: 
     3688                y = float(f(x)) 
     3689            except (ZeroDivisionError, TypeError, ValueError), msg: 
     3690                sage.misc.misc.verbose("%s\nUnable to compute f(%s)"%(msg, x), 1) 
     3691                # give up for this branch 
     3692                return [] 
     3693            # this distance calculation is not perfect. 
     3694            if abs((p1[1] + p2[1])/2 - y) > adaptive_tolerance: 
     3695                return refine(f, p1, (x, y), level+1) + [(x, y)] +\ 
     3696                        refine(f, (x, y), p2, level+1) 
     3697            else: 
     3698                return [] 
     3699         
    36613700        while i < len(data) - 1: 
    3662             if abs(data[i+1][1] - data[i][1]) > max_bend: 
    3663                 x = float((data[i+1][0] + data[i][0])/2) 
    3664                 try: 
    3665                     y = float(f(x)) 
    3666                     data.insert(i+1, (x, y)) 
    3667                 except (ZeroDivisionError, TypeError, ValueError), msg: 
    3668                     sage.misc.misc.verbose("%s\nUnable to compute f(%s)"%(msg, x),1) 
    3669                     exceptions += 1 
    3670                 j += 1 
    3671                 if j > plot_division: 
    3672                     break 
    3673             else: 
     3701            for p in refine(f, data[i], data[i+1]): 
     3702                data.insert(i+1, p) 
    36743703                i += 1 
    3675  
     3704            i += 1 
     3705         
    36763706        if (len(data) == 0 and exceptions > 0) or exceptions > 10: 
    36773707            sage.misc.misc.verbose("WARNING: When plotting, failed to evaluate function at %s points."%exceptions, level=0) 
    36783708            sage.misc.misc.verbose("Last error message: '%s'"%msg, level=0)