Implementing Adaptive function plotting

Question

I'm attempting to implement the Adaptive function plotting algorithm using the psuedocode from these two examples (both examples the same really)

https://www.andr.mu/logs/acquiring-samples-to-plot-a-math-function-adaptive/ http://yacas.readthedocs.io/en/latest/book_of_algorithms/basic.html

The problems I've encountered (from my probably incorrect implementation) are:

(1) Duplicated coordinates are being created. I know this is because of the splitting, where both ends are kept for each spilt, so the end of one interval is the start of the other interval - same x value evaluated twice. But is there a way to setup the algorithm to avoid deplicated coordinates. I could avoid adding the start or end coordinate as a dirty fix (see comment in code below) but then I'd be missing one them for the whole interval.

(2) Some parts of the plot are missing coordinates for what is essentially a symmetrical function, which is strange. The algorithm should work the same way for both sides, yet it doesn't. This starts to happen when depth >= 6, it does an extra split to the right side of the function and nothing for the left side around the origin. The rest of the coordinates away from the origin seem to match up. The right side seems to get more splits than the left side overall.

Problem (2)

My Implementation of algorithm

static List<Double[]> linePoints;

public static void main(String[] args) {
    linePoints =  new ArrayList<>();

    double startX = -50;
    double endX = 50;

    sampling(startX, endX, depth, tolerance);

    /* Print all points to be plotted - x,y coordinates */
    for (Double[] point : linePoints) {
        System.out.println(point[0]+","+point[1]);
    }
}

/* math function */
public static double f(double x){
    return x*Math.sin(x);
}


static int depth = 6; /* 8 */
static double tolerance = 0.005; /* just a guess */

/* Adaptive sampling algorithm */
/* mostly followed along 2st website and used 1st website variable names  */
public static void sampling(double xa, double xc, int depth, double tolerance){
    /* step (1) of 2nd website - determine mid-intervals */
    double xb = (xa+xc)/2;   /* (xc-xa)/2; tried these from 1st website - didn't work out */
    double xab = (xa+xb)/2;  /* (xb-xa)/2; */
    double xbc = (xb+xc)/2;  /* (xc-xb)/2; */

    /* evaluate the above points using math function - store in array */
    double[] points = new double[5];
    points[0] = f(xa); points[1] = f(xab); points[2] = f(xb); points[3] = f(xbc); points[4] = f(xc);

    /* step (2) of 2nd website */
    if (depth <= 0){
        linePoints.add(new Double[]{xa, points[0]});  /* either I comment out this line for dirty fix */
        linePoints.add(new Double[]{xab, points[1]});
        linePoints.add(new Double[]{xb, points[2]});
        linePoints.add(new Double[]{xbc, points[3]});
        linePoints.add(new Double[]{xc, points[4]});  /* or comment out this line */
    } else {
        /* step (3) of 2nd website */
        int counter = 0;

        for (int i = 1; i < points.length-1; i++){
            /* Check if prev, current, next values are infinite or NaN */
            if (    (Double.isInfinite(points[i-1]) || Double.isNaN(points[i-1])) ||
                    (Double.isInfinite(points[i]) || Double.isNaN(points[i])) ||
                    (Double.isInfinite(points[i+1]) || Double.isNaN(points[i+1]))){
                counter++;
                continue;
            }

            /* Determine the fluctuations - if current is < or > both it's left/right neighbours */
            boolean middleLarger = (points[i] > points[i-1]) && (points[i] > points[i+1]);
            boolean middleSmaller = (points[i] < points[i-1]) && (points[i] < points[i+1]);

            if (middleLarger || middleSmaller){
                counter++;
            }
        }

        if (counter <= 2){  /* at most 2 */
            /* Newton-Cotes quadratures - check if smooth enough */
            double f1 = (3d/8d)*points[0]+(19d/24d)*points[1]-(5d/24d)*points[2]+(1d/24d)*points[3];    /* add 'd' to end of number, otherwise get 0 always */
            double f2 = (5d/12d)*points[2]+(2d/3d)*points[3]-(1d/12d)*points[4];

         if (Math.abs(f1-f2) < tolerance * f2){
                linePoints.add(new Double[]{xa, points[0]});
                linePoints.add(new Double[]{xab, points[1]});
                linePoints.add(new Double[]{xb, points[2]});
                linePoints.add(new Double[]{xbc, points[3]});
                linePoints.add(new Double[]{xc, points[4]});
            } else {
                /* not smooth enough - needs more refinement */
                depth--;
                tolerance *= 2;
                sampling(xa, xb, depth, tolerance);
                sampling(xb, xc, depth, tolerance);
            }

        } else {
            /* else (count > 2), that means further splittings are needed to produce more accurate samples */
            depth--;
            tolerance *= 2;
            sampling(xa, xb, depth, tolerance);
            sampling(xb, xc, depth, tolerance);
        }
    }
}

FIX - Modifications to my code

Looking at Gene's example and multiplying the tolerance by 0.5 instead of 2 seemed to fix problem (2)

Genes example is a better and cleaner implementation of this algorithm and handles the duplicated coordinates

Answer 1

I think you've implemented faithfully, but the algorithm is broken. Asymmetric quadrature can easily give asymmetric results. I get the same kinds of weirdness.

But you can eliminate duplicate points by maintaining them in a sorted set and using the invariant that the endpoints have already been inserted when the recursive analyzer runs.

Here's a simplification using modern Java features more fully:

import static java.lang.Double.compare;
import static java.lang.Double.isFinite;
import static java.lang.Math.PI;

import java.util.List;
import java.util.SortedSet;
import java.util.TreeSet;
import java.util.function.DoubleUnaryOperator;
import static java.util.stream.Collectors.toList;
import java.util.stream.DoubleStream;

public class AdaptivePlot {
  private final DoubleUnaryOperator f;
  private final double a;
  private final double c;
  private final SortedSet<Point> plot = new TreeSet<>((s, t) -> compare(s.x, t.x));

  public AdaptivePlot(DoubleUnaryOperator f, double a, double c) {
    this.f = f;
    this.a = a;
    this.c = c;
  }

  public static class Point {
    final double x, y;
    public Point(double x, double y) {
      this.x = x;
      this.y = y;
    }
  }

  public AdaptivePlot computePlot(int depth, double eps) {
    plot.clear();
    Point pa = pointAt(a);
    Point pc = pointAt(c);
    plot.add(pa);
    plot.add(pc);
    computePlot(pa, pc, depth, eps);
    return this;
  }

  public List<Point> getPlot() {
    return plot.stream().collect(toList());
  }

  private Point pointAt(double x) {
    return new Point(x, f.applyAsDouble(x));
  }

  private void computePlot(Point pa, Point pc, int depth, double eps) {
    Point pb = pointAt(0.5 * (pa.x + pc.x));
    Point pa1 = pointAt(0.5 * (pa.x + pb.x));
    Point pb1 = pointAt(0.5 * (pb.x + pc.x));
    plot.add(pb);
    if (depth > 0 && 
        (oscillates(pa.y, pa1.y, pb.y, pb1.y, pc.y) 
            || unsmooth(pa.y, pa1.y, pb.y, pb1.y, pc.y, eps))) {
      computePlot(pa, pb, depth - 1, 2 * eps);
      computePlot(pb, pc, depth - 1, 2 * eps);
    }
    plot.add(pa1);
    plot.add(pb1);
  }

  private static boolean oscillates(
      double ya, double ya1, double yb, double yb1, double yc) {
    return isOscillation(ya, ya1, yb) 
        && isOscillation(ya1, yb, yb1) 
        && isOscillation(yb, yb1, yc);
  }

  private static boolean isOscillation(double ya, double yb, double yc) {
    return !isFinite(ya) || !isFinite(yb) || !isFinite(yc) 
        || (yb > ya && yb > yc) || (yb < ya && yb < yc);
  }

  private static boolean unsmooth(
      double ya, double ya1, double yb, double yb1,double yc, double eps) {
    double y0 = DoubleStream.of(ya, ya1, yb, yb1, yc).min().getAsDouble();
    double [] yg = DoubleStream.of(ya, ya1, yb, yb1, yc).map(y -> y - y0).toArray();
    double q4 = quadrature(yg[0], yg[1], yg[2], yg[3]);
    double q3 = quadrature(yg[2], yg[3], yg[4]);
    return Math.abs(q4 - q3) > eps * q3;
  }

  private static double quadrature(double y0, double y1, double y2, double y3) {
    return 3d/8d * y0 + 19d/24d * y1 - 5d/24d * y2 + 1d/24d * y3;
  }

  private static double quadrature(double y0, double y1, double y2) {
    return 5d/12d * y0 + 2d/3d * y1 - 1d/12d * y2;
  }

  public static void main(String [] args) {
    List<Point> plot = new AdaptivePlot(x -> x * Math.sin(x), -2d * PI, 2d * PI)
        .computePlot(6, 0.005).getPlot();
    for (Point p : plot) {
      System.out.println(p.x + "\t" + p.y);
    }
  }
}

Answer 2

I personally find this method horrible.

You can get nice results with a simple rule: subdivide the interval as long as the distance of the middle point to the chord exceeds a small threshold.

Caution: to avoid scaling effects, you should express all geometric quantities in mm (with appropriate scale factors on both axes).

Subdivide left, middle, right:
    if DistancePointLine (middle, f(middle)), (left, f(left)), (right, f(right)) < Tolerance:
        DrawLine (left, f(left), (right, f(right))
    else
        Subdivide left, (left + middle) / 2, middle
        Subdivide middle, (middle + right) / 2, right

The method can fail in symmetric situations, around inflection points. To cope, you can force one extra subdivision even if the tolerance has been met.