Ticket #3813: 3813-anakha-adaptive-plot-v3.patch

a/sage/plot/plot.py

@@ -3414 +3414 @@
-
-    PLOT OPTIONS:
-    The plot options are
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
-        xmin -- starting x value
-        xmax -- ending x value
-        color -- an rgb-tuple (r,g,b) with each of r,g,b between 0 and 1, or
@@ -3461 +3469 @@
-        Graphics object consisting of 1 graphics primitive
-        sage: len(P)     # number of graphics primitives
-        1
-
-00
+ (random)
+24
-        sage: P          # render
-        
-        sage: P = plot(sin, (0,10), plot_points=10); print P
-        Graphics object consisting of 1 graphics primitive
-        sage: len(P[0])  # random output
-80
+32
-        sage: P          # render
-
-    We plot with randomize=False, which makes the initial sample 
-    points evenly spaced (hence always the same).  Adaptive plotting 
-plot_divi
-plot_divi
+adaptive_recur
+adaptive_recur
-        sage: list(p[0])
-        [(0.0, 1.0), (1.0, 1.0), (2.0, 1.0), (3.0, 1.0)]
-
@@ -3488 +3496 @@
-        sage: plot([sin(n*x) for n in [1..4]], (0, pi))
-
-
-
-
-
+
+
-        sage: plot(sin(1/x), (x, -1, 1))
-
-    With the \code{plot_points} option you can increase the number
-    of sample points, to obtain a more accurate plot. 
-        sage: plot(sin(1/x), (x, -1, 1), plot_points=1000)
-
-    Note that the independent variable may be omitted if there is no
-    ambiguity:
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-    The actual sample points are slightly randomized, so the above
-    plots may look slightly different each time you draw them.
@@ -3543 +3564 @@
-    def _reset(self):
-        o = self.options
-        o['plot_points'] = 200
-plot_division'] = 1000 
-max_bend'] = 0.1       
+adaptive_tolerance'] = 0.01
+adaptive_recursion'] = 5
-        o['rgbcolor'] = (0,0,1)   
-
-    def __repr__(self):
-
+r
-        Returns a string representation of this PlotFactory object.
-
-        TESTS:
@@ -3612 +3633 @@
-        plot_points = int(options['plot_points'])
-        del options['plot_points']
-        x, data = var_and_list_of_values(xrange, plot_points)
-        data = list(data)
-        xmin = data[0]
-        xmax = data[-1]
-
@@ -3638 +3658 @@
-            # randomize is true
-            if i > 0 and i < plot_points-1:
-                if randomize:
-
-
-
+
+
-            elif i == plot_points-1: 
-                xi = xmax  # guarantee that we get the last point.
-    
+
-            try:
-                data[i] = (float(xi), float(f(xi)))
-
+
+
+
+
-                sage.misc.misc.verbose("%s\nUnable to compute f(%s)"%(msg, x),1)
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        data = [data[i] for i in range(len(data)) if i not in exception_indices]
-
+
+
+
+
+
-        # adaptive refinement
-        i, j = 0, 0
-
-
-
-
+
+
+
+
+
-        while i < len(data) - 1:
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
-                i += 1
-
+
+
-        if (len(data) == 0 and exceptions > 0) or exceptions > 10:
-            sage.misc.misc.verbose("WARNING: When plotting, failed to evaluate function at %s points."%exceptions, level=0)
-            sage.misc.misc.verbose("Last error message: '%s'"%msg, level=0)
@@ -4449 +4490 @@
-                g.append(fast_float(f, str(xvar), str(yvar)))
-            
-    return g, xstep, ystep, xrange, yrange
+
+def adaptive_refinement(f, p1, p2, adaptive_tolerance=0.01, adaptive_recursion=10, level=0):
+    r"""
+    The adaptive refinement algorithm for plotting a function f. See the
+    docstring for plot or PlotFactory for a description of the algorithm.
+
+    INPUT:
+        f -- a function of one variable
+        p1, p2 -- two points to refine between
+        adaptive_recursion -- (default: 10) how many levels of recursion to go 
+                              before giving up when doing adaptive refinement.
+                              Setting this to 0 disables adaptive refinement.
+        adaptive_tolerance -- (default: 0.01) how large a difference should be
+                              before the adaptive refinement code considers 
+                              it significant.  See the documentation for
+                              plot() for more information.
+    
+    OUTPUT:
+        list -- a list of points to insert between p1 and p2 to get
+                a better linear approximation between them
+
+    TESTS:
+        sage: from sage.plot.plot import adaptive_refinement
+        sage: adaptive_refinement(sin, (0,0), (pi,0), adaptive_tolerance=0.01, level=10)
+        []
+        sage: adaptive_refinement(sin, (0,0), (pi,0), adaptive_tolerance=0.01)
+        [(0.125000000000000*pi, 0.38268343236508978), (0.187500000000000*pi, 0.55557023301960218), (0.250000000000000*pi, 0.707106781186547...), (0.312500000000000*pi, 0.83146961230254524), (0.375000000000000*pi, 0.92387953251128674), (0.437500000000000*pi, 0.98078528040323043), (0.500000000000000*pi, 1.0), (0.562500000000000*pi, 0.98078528040323043), (0.625000000000000*pi, 0.92387953251128674), (0.687500000000000*pi, 0.83146961230254546), (0.750000000000000*pi, 0.70710678118654757), (0.812500000000000*pi, 0.55557023301960218), (0.875000000000000*pi, 0.38268343236508989)]
+
+        This shows that lowering adaptive_tolerance and raising
+        adaptive_recursion both increase the number of subdivision points:
+
+        sage: x = var('x')
+        sage: f = sin(1/x)
+        sage: n1 = len(adaptive_refinement(f, (0,0), (pi,0), adaptive_tolerance=0.01)); n1
+        79
+        sage: n2 = len(adaptive_refinement(f, (0,0), (pi,0), adaptive_recursion=5, adaptive_tolerance=0.01)); n2
+        15
+        sage: n1 > n2
+        True
+
+        sage: n3 = len(adaptive_refinement(f, (0,0), (pi,0), adaptive_tolerance=0.005)); n3
+        88
+        sage: n1 < n3
+        True
+    """
+    if level >= adaptive_recursion:
+        return []
+    x = (p1[0] + p2[0])/2.0
+    try:
+        y = float(f(x))
+    except (ZeroDivisionError, TypeError, ValueError, OverflowError), msg:
+        sage.misc.misc.verbose("%s\nUnable to compute f(%s)"%(msg, x), 1)
+        # give up for this branch
+        return []
+    # this distance calculation is not perfect.
+    if abs((p1[1] + p2[1])/2.0 - y) > adaptive_tolerance:
+        return adaptive_refinement(f, p1, (x, y), 
+                        adaptive_tolerance=adaptive_tolerance,
+                        adaptive_recursion=adaptive_recursion,
+                        level=level+1) \
+               + [(x, y)] + \
+               adaptive_refinement(f, (x, y), p2,
+                        adaptive_tolerance=adaptive_tolerance,
+                        adaptive_recursion=adaptive_recursion,
+                        level=level+1)
+    else:
+        return []