// ***** This file is automatically generated from LinearTernaryCore.java.jpp

package daikon.inv.ternary.threeScalar;

import daikon.*;
import daikon.inv.*;
import daikon.inv.binary.twoScalar.*;
import daikon.inv.unary.scalar.*;
import java.io.Serializable;
import java.util.*;
import java.util.logging.Level;
import java.util.logging.Logger;
import org.checkerframework.checker.lock.qual.GuardSatisfied;
import org.checkerframework.checker.nullness.qual.Nullable;
import org.checkerframework.dataflow.qual.Pure;
import org.checkerframework.dataflow.qual.SideEffectFree;
import org.plumelib.util.ArraysPlume;
import org.plumelib.util.MathPlume;
import typequals.prototype.qual.NonPrototype;
import typequals.prototype.qual.Prototype;

/**
 * The LinearTernaryCore class is acts as the backend for the invariant (ax + by + cz + d = 0) by
 * processing samples and computing coefficients. The resulting coefficients a, b, c, and d are
 * mutually relatively prime, and the coefficient a is always p.
 */

// Originally, LinearTernaryCore code was not dealing with degenerate
// linear ternary invariants correct, namely when the plane
// (ax+by+cz+d=0) is parallel to one of the axes, so one of the
// coefficients was 0.  In this case, the invariant can be described
// by a binary invariant (eliminate the 0 coefficient variable).
// LinearTernaryCore's inability to deal with this degenerate form was
// (and still is) masked by the suppression of a linear ternary
// invariant via a binary invariant and a constant (the code threw an
// exception when seeing this case and falsifies the invariant).  As
// more samples are processed, the ternary invariant may "evolve" from
// a degenerate form to a full ternary invariant.  The
// linearintervention method deals with the degenerate form and the
// planarintervention focuses on the full ternary invariant.

@SuppressWarnings({"nullness", // ***** TEMPORARY:  def_points is confused (is packed initially with non-nulls at beginning, but is not all non-nulls; and later may not be packed); needs review/refactoring/rewriting
    "UnnecessaryParentheses"  // generated code, parens are sometimes necessary
    })
public final class LinearTernaryCoreFloat implements Serializable, Cloneable {
  static final long serialVersionUID = 20030822L;

  /** Debug tracer. */
  static final Logger debug = Logger.getLogger("daikon.inv.ternary.threeScalar.LinearTernaryCoreFloat");

  // ax + by + cz + d = 0; first argument is x, second is y, third is z
  // Although the coefficients are guaranteed to be integers, we use doubles to
  // store the values because regression tests have shown that a long is
  // not sufficient
  public double a = 0, b = 0, c = 0, d = 0;
  public double min_a, max_a, min_b, max_b, min_c, max_c, min_d, max_d;
  public double separation = 0;

  // coefficients for the invariant in 3D space
  public double[] coefficients = new double[4];

  // The values of the flags determine indicate whether a line exists
  // with the current samples (NO_LINE) and which variables are
  // constant (constant variables indicate that the 'plane' is
  // parallel to that axis).  For example, XY_CONSTANT means that the
  // so far, the x and y values have all been constant and only z has
  // been varying so so the plane is actually a line.
  // NO_CONSTANT_VARIABLES indicates that the values have not been
  // constant but do for a plane between variables.
  public enum Flag {NO_LINE, XY_CONSTANT, XZ_CONSTANT, YZ_CONSTANT,
                    X_CONSTANT, Y_CONSTANT, Z_CONSTANT, NO_CONSTANT_VARIABLES};

  // The flag indicates whether a line exists in the 3D space instead
  // of a plane If the flag value is not NO_LINE, then a line does
  // currently exist among the values seen and the flag specifies
  // which variables are constant to make the invariant a line rather
  // than a plane.  The value NO_LINE indicates either that not enough
  // values have been seen to make any decision or that a plane
  // exists.
  public Flag line_flag = Flag.NO_LINE;

  // Inner class to store points in for more convenient access.
  public static class Point implements Serializable, Cloneable {
    static final long serialVersionUID = 20050923L;

    public double x;
    public double y;
    public double z;
    public Point() {
      this(0, 0, 0);
    }
    public Point(double x, double y, double z) {
      this.x = x; this.y = y; this.z = z;
    }
    @Pure
    public boolean equals(double x, double y, double z) {
      return (this.x == x) && (this.y == y) && (this.z == z);
    }
    @SideEffectFree
    @Override
    protected Point clone(@GuardSatisfied Point this) throws CloneNotSupportedException {
      return (Point) super.clone();
    }
    @SideEffectFree
    @Override
    public String toString(@GuardSatisfied Point this) {
      return "(" + x + ", " + y + " ," + z + ")";
    }
  }

  // Points that define the plane.  The first three are normally the
  // current definition.  The next is a new point for consideration.
  // However, any slot may be set to null at any time, and points may be
  // duplicated at any time.
  // (The invariant appears to be a bit stronger than that, since there are
  // places that dereference the array unconditionally, assuming that its
  // elements are non-null.  This should be clarified/refactored.)
  public @Nullable Point[] def_points = new Point[MINTRIPLES];

  public Invariant wrapper;

  // The number of distinct values (not samples) seen.
  public int values_seen = 0;

  static final int MINTRIPLES = 5;

  public LinearTernaryCoreFloat(Invariant wrapper) {
    this.wrapper = wrapper;
  }

  @SideEffectFree
  @Override
  public LinearTernaryCoreFloat clone(@GuardSatisfied LinearTernaryCoreFloat this) {
    try {
      LinearTernaryCoreFloat result = (LinearTernaryCoreFloat) super.clone();
      result.def_points = new Point[MINTRIPLES];
      for (int i = 0; i < MINTRIPLES; i++) {
        Point p = def_points[i];
        if (p != null) {
          result.def_points[i] = p.clone();
        }
      }
      return result;
    } catch (CloneNotSupportedException e) {
      throw new Error(); // can't happen
    }
  }

  /**
   * Reorganize our already-seen state as if the variables had shifted order underneath us
   * (rearrangement given by the permutation).
   */
  public void permute(int[] permutation) {
    assert permutation.length == 3;
    assert ArraysPlume.fnIsPermutation(permutation);
    // Fix a, b, c
    double[] clever = new double[] { a, b, c };
    double[] pclever = new double[3];
    pclever[permutation[0]] = clever[0];
    pclever[permutation[1]] = clever[1];
    pclever[permutation[2]] = clever[2];
    if (isActive()) {
      // Seen enough values, so permuting a, b, c is useful
      if ( a == 0 || b == 0 || c == 0) {
        throw new Error("Active LTs never have 0 coefficients");
        // // Can't handle this form once rotated.  But if we've seen
        // // enough values, yet a, b or c is zero, then this is covered by
        // // a LinearBinary or a constant, so let me just destroy myself.
        // values_seen = Integer.MAX_VALUE;
        // a = b = c = d = 0;
        // if (debug.isLoggable(Level.FINE)) {
        //   debug.fine ("  Ternary invariant destroyed because a, b, or c = 0");
        // }
        // return;
      } else {
        // double d = -1.0 / pclever[2];
        a = pclever[0];
        b = pclever[1];
        c = pclever[2];

        // we still need to guarantee a positive x coefficient
        if (a < 0) {
          a = -a;
          b = -b;
          c = -c;
          d = -d;
        }
      }
    }

    // Fix caches
    {
      double[] temp = new double[3];
      for (int i = 0; i < MINTRIPLES; i++) {
        Point p = def_points[i];
        if (p == null) {
          continue;
        }
        wrapper.log("orig def_points[%s] = %s", i, p);
        temp[permutation[0]] = p.x;
        temp[permutation[1]] = p.y;
        temp[permutation[2]] = p.z;
        p.x = temp[0];
        p.y = temp[1];
        p.z = temp[2];
        wrapper.log("permuted def_points[%s] = %s", i, p);
      }
    }
    // Assert that caches sync with a,b,c?
    // This would be a good sanity check, but it would be nontrivial
    // because we don't keep track of when a, b, and c are really
    // valid for all the of the points we've cached. If we were
    // falsified and then resurrected, we might have samples that
    // didn't fit even before the permutation. -smcc
  }

  /**
   * Returns whether or not the invariant is currently active. We become active after MINTRIPLES
   * values have been seen and a line is not calculated. In addition, the coefficients (a, b, c)
   * must all be non-zero, else the invariant would be described by a LinearBinary or constant
   * invariant. Before that, a, b, and c are uninitialized.
   */
  // this is the check used throughout the file as a shortcut because
  // the invariant must have seen enough samples and not be a line for all the
  // other calculations and comparisons done to be sensible, meaning that
  // a full linear ternary invariant exists, not a degenerate one
  @Pure
  public boolean isActive(@GuardSatisfied LinearTernaryCoreFloat this) {
    return line_flag == Flag.NO_LINE && values_seen >= MINTRIPLES && a != 0 && b != 0 && c != 0;
  }

  /**
   * Sets up the invariant from a LinearBinary invariant and a constant value for the third
   * variable. Points are taken from the LinearBinary cache and its min_x/y and max_x/y points and
   * combined with the constant value.
   *
   * @return InvariantStatus.NO_CHANGE if the invariant is valid, InvariantStatus.FALSIFIED if one
   * of the points invalidated the LinearTernary invariant
   */
  public InvariantStatus setup(LinearBinaryFloat lb, VarInfo con_var,
                                double con_val) {

    if (Debug.logOn()) {
      wrapper.log("setup from lb %s con var %s con_val %s",
                   lb, con_var, con_val);
    }

    int con_index = con_var.varinfo_index;
    int lb_v1_index = lb.ppt.var_infos[0].varinfo_index;
    int lb_v2_index = lb.ppt.var_infos[1].varinfo_index;

    double[] con_vals = new double [lb.core.values_seen + 2];
    double[] lb1_vals = new double [lb.core.values_seen + 2];
    double[] lb2_vals = new double [lb.core.values_seen + 2];

    int mi = lb.core.values_seen;
    for (int i = 0; i < mi; i++) {
      con_vals[i] = con_val;
      lb1_vals[i] = lb.core.x_cache[i];
      lb2_vals[i] = lb.core.y_cache[i];
    }
    if (lb.isActive()) {
      con_vals[mi] = con_val;
      lb1_vals[mi] = lb.core.min_x;
      lb2_vals[mi] = lb.core.min_y;
      con_vals[mi + 1] = con_val;
      lb1_vals[mi + 1] = lb.core.max_x;
      lb2_vals[mi + 1] = lb.core.max_y;
      mi += 2;
    }

    InvariantStatus sts = InvariantStatus.NO_CHANGE;

    for (int i = 0; i < mi; i++ ) {
      if (con_index < lb_v1_index) {
        sts = add_modified(con_vals[i], lb1_vals[i], lb2_vals[i], 1);

      } else if (con_index < lb_v2_index) {
        sts = add_modified(lb1_vals[i], con_vals[i], lb2_vals[i], 1);

      } else {
        sts = add_modified(lb1_vals[i], lb2_vals[i], con_vals[i], 1);

      }
      if (sts != InvariantStatus.NO_CHANGE) {
        break;
      }
    }

    if (sts != InvariantStatus.NO_CHANGE) {
      System.out.println("lb.core.values_seen=" + lb.core.values_seen);
      for (int i = 0; i < mi; i++ ) {
        System.out.printf("LTCore: vals %s %s %s%n", con_vals[i], lb1_vals[i],
                           lb2_vals[i]);
      }
      System.out.println("in inv " + wrapper.format() + " " + wrapper.ppt);
      assert sts == InvariantStatus.NO_CHANGE;
    }

    return sts;
  }

  /**
   * Sets up the invariant from a OneOf and two constants. Points are taken from the OneOf cache
   * and the constant values.
   *
   * @return InvariantStatus.NO_CHANGE if the invariant is valid, or InvariantStatus.FALSIFIED if
   * one of the points invalidated the LinearTernary invariant
   */

  public InvariantStatus setup(OneOfFloat oo, VarInfo v1, double con1,
                                VarInfo v2, double con2) {

    int oo_index = oo.ppt.var_infos[0].varinfo_index;
    int con1_index = v1.varinfo_index;
    int con2_index = v2.varinfo_index;

    InvariantStatus sts = InvariantStatus.NO_CHANGE;
    for (int i = 0; i < oo.num_elts(); i++ ) {
      if (oo_index < con1_index) {
        sts = add_modified(((Double)oo.elt(i)).doubleValue(), con1, con2, 1);
      } else if (oo_index < con2_index) {
        sts = add_modified(con1, ((Double)oo.elt(i)).doubleValue(), con2, 1);
      } else {
        sts = add_modified(con1, con2, ((Double)oo.elt(i)).doubleValue(), 1);
      }
      if (sts != InvariantStatus.NO_CHANGE) {
        break;
      }
    }

    if (Debug.logOn()) {
      wrapper.log("setup from OneOf %s v1=%s v2=%s status = %s",
                   oo, con1, con2, sts);
    }

    return sts;
  }

  /**
   * Looks for points that define a plane (ax + by + cz + d = 0). Collects MINTRIPLE points before
   * attempting to define a line through the points. Once the equation for the line is found, each
   * subsequent point is compared to it. If the point does not match, recalculates the maximally
   * separated points and attempts to fit the points to the new line. If the points do not fit the
   * line, attempts to define the plane (to hopefully get at least some spread between the points,
   * so that small errors don't get magnified). Once the equation for the plane is found, each
   * subsequent point is compared to it. If the point does not match the point is examined to see if
   * it would is maximally separated when compared to the points originally used to define the
   * plane. If it is, it is used to recalcalulate the coefficients (a, b, c). If those
   * coefficients are relatively close to the original coefficients (within the ratio defined by
   * Global.fuzzy) then the new coefficients are used.
   *
   * @see org.plumelib.util.FuzzyFloat
   */
  public InvariantStatus add_modified(double x, double y, double z, int count) {

    if (Debug.logDetail()) {
      wrapper.log("Adding point, x=%s y=%s z=%s to invariant", x, y, z);
    }

    if (values_seen < MINTRIPLES) {
      // We delay computation of a and b until we have seen several triples
      // so that we can compute a and b based on a far-separated triple.  If
      // the points in a triple are nearby, then roundoff errors in the
      // computation of the slope can be non-negligible.

      // skip points we've already seen
      for (int i = 0; i < values_seen; i++) {
        if (def_points[i].equals(x, y, z)) {
          return InvariantStatus.NO_CHANGE;
        }
      }

      def_points[values_seen] = new Point(x, y, z);
      wrapper.log("Add: (%s, %s, %s)", x, y, z);
      values_seen++;
      wrapper.log("Values seen: %s", values_seen);

      // check to see if we've seen enough values to create the equation
      if (values_seen == MINTRIPLES) {

        // Try to fit the points to a line first.
        // The points may form a degenerate linear ternary invariants
        // especially if one of the variable is constant
        linearIntervention(def_points);

        // if line can not be formed, try to form a plane
        if (line_flag == Flag.NO_LINE) {
          InvariantStatus stat = planarIntervention(def_points);
          return stat;
        } else {

          // points fit in a line
          return InvariantStatus.NO_CHANGE;
        }
      } else {

        // still haven't seen enough values
        return InvariantStatus.NO_CHANGE;
      }
    } else {

      // at this point either a line or a plane must have been formed

      // if line already broken, fit to plane equation
      if (line_flag == Flag.NO_LINE) {
        // If the new value doesn't fit the equation
        if (!Global.fuzzy.eq(-d, a * x + b * y + c * z)) {

          // Try to find small changes that will fit better.
          if (!try_new_equation(x, y, z)) {
            if (Invariant.logOn() || debug.isLoggable(Level.FINE)) {
              wrapper.log(debug, "destroying  (" + wrapper.format()
                           + ") where x=" + x + " y=" + y + " z=" + z
                           + " a=" + a + " b=" + b + " c=" + c
                           + " values_seen=" + values_seen);
            }
            // jiggling the values did not work
            return InvariantStatus.FALSIFIED;
          }
        } else {

          // point fits the current plane equation
          return InvariantStatus.NO_CHANGE;
        }
      } else {

        // try to fit to the current line
        if (!try_points_linear((x), (y), (z))) {
          // put the crucial point in, so that it can be tested against the new line
          def_points[MINTRIPLES - 1] = new Point(x, y, z);

          // try to fit the points onto a line again (not sure if this
          // is really necessary, maybe the fuzzy equals really does
          // make a difference)
          linearIntervention(def_points);
        }

        // new point has broken the line, try to fit to a plane
        if (line_flag == Flag.NO_LINE) {
          InvariantStatus stat = planarIntervention(def_points);
          return stat;
        } else {

          // fuzzy equals does make a difference
          return InvariantStatus.NO_CHANGE;
        }

      }

    }
    // should never get here because values_seen is either < or >= MINTRIPLES
    // return InvariantStatus.NO_CHANGE;
    throw new Error("at end of add_modified");
  }

  /**
   * Attempts to fit the points in def_point to a plane and calculates the coefficients (a, b, c, d)
   * if possible.
   *
   * @param def_point array of points. Must have at least 3 elements.
   *
   * @return the status of the invariant
   *     (whether the plane fitting was successful)
   */
  private InvariantStatus planarIntervention(@Nullable Point[] def_point) {

    // make a copy of the array so that the points don't get clobbered
    // because we need to check all the points in def_points against the
    // coefficients calculated

    @Nullable Point[] dummy = new @Nullable Point[def_point.length];
    for (int i = 0; i < def_point.length; i++) {
      dummy[i] = def_point[i];
    }

    // Find the three points with the maximum separation
    maxsep_triples(dummy);
    // Now, first 3 values of dummy are the maximum-separated 3 points.

    // Null points only show up in the cross-checker, not the regression
    // tests, as of 5/25/2010.
    // Sometimes the points are still null.
    if (dummy[0] == null) {
      return InvariantStatus.FALSIFIED;
    }

    // Calculate the coefficients of the equation (a, b, c, d).
    // Requires that the first 3 elements of dummy must be non-null.
    double[] coef = calc_tri_linear(dummy);

    a = coef[0];
    b = coef[1];
    c = coef[2];
    d = coef[3];

    wrapper.log("Calculated tri linear coeff: (%s), b (%s), c(%s), and d(%s)",
                 a, b, c, d);
    // If one of these coefficients is zero (except for d), this
    // should be a LinearBinary, not a LinearTernary, term.  (It
    // might not show up as LinearBinary because there might not
    // have been enough samples; but a random varying third
    // variable can create enough samples.)  Also, throw out the
    // invariant if any of the coefficients would be not-a-number.
    if ( a == 0|| b == 0 || c == 0
        || Double.isNaN(a) || Double.isNaN(b) || Double.isNaN(c)) {

      wrapper.log("problematic coefficient: invariant falsified");
      return InvariantStatus.FALSIFIED;
    }

    // Check all values against a, b, and c.
    if (!wrapper.is_false()) {
      for (int i = 0; i < values_seen; i++) {

        // Global.fuzzy.eq works by comparing ratios between numbers,
        // so comparing the sum of everything to 0 won't work.
        //  wrapper.log("point: " + def_point[i].x + " "
        //  + def_point[i].y + " " + def_point[i].z);
        //  wrapper.log("compare: " + -c * def_point[i].z);
        //  wrapper.log( "to: " + (a*def_point[i].x
        //  + b*def_point[i].y + d));

        wrapper.log("Trying at index %s: 0 != %s*%s+%s*%s+%s*%s+%s",
                     i, a, def_point[i].x,
                        b, def_point[i].y,
                        c, def_point[i].z,
                        d);

        if (!Global.fuzzy.eq(-d, a * def_point[i].x + b * def_point[i].y + c * def_point[i].z)) {

          if (Debug.logOn() || debug.isLoggable(Level.FINE)) {
            wrapper.log(debug, "Destroying at index " + i + ": 0 != "
                         + a + "*" + def_point[i].x
                         + "+" + b + "*" + def_point[i].y + "+" + c
                         + "*" + def_point[i].z
                         + "+" + d);
          }

          return InvariantStatus.FALSIFIED;
        }
      }
      if (Debug.logOn()) {
        wrapper.log("equation = %s*x %s*y %s*z = %s", a, b, c, -d);
      }

      // Discard the points not used to define the coefficients.
      for (int ii = 3; ii < MINTRIPLES; ii++) {
        def_point[ii] = null;
      }
    }

    // if it passes all the checks, then no change to the Invariant
    return InvariantStatus.NO_CHANGE;
  }

  /**
   * Attempts to fit the points in def_point to a line and calculates the corresponding line
   * classification (line_flag) and coefficients (coefficients[]) if possible. If not (points in
   * def_point form a plane) resets line_flag to Flag.NO_LINE and coefficients[] elements to -999.
   *
   * @param def_point array of points. Must have at least 2 elements.
   *
   */
  private void linearIntervention(@Nullable Point[] def_point) {

    // find the max separation points
    Point[] maxsep_doubles = maxsep_doubles(def_point);

    // maxsep_doubles[0] could be null when enough points have a NaN
    if (maxsep_doubles == null) {
      // line_flag == Flag.NO_LINE so the invariant
      // will try to fit a plane to the points and fail
      return;
    } else {
      // find the equation (coefficients now stored in coefficients[])
      calc_bi_linear(maxsep_doubles[0], maxsep_doubles[1]);

      // This seems to imply that all elements of def_points are non-null.
      // Is that true?
      // check the other points against the coefficients
      for (int i = 0; i < def_point.length; i++) {
        // try to fit the current point to the current line
        boolean ok = try_points_linear((def_point[i].x), (def_point[i].y), (def_point[i].z));

        if (!ok) {
          line_flag = Flag.NO_LINE;
          break;
        }
      }
    }

  }

  /**
   * Attempts to fit the new point (x,y,z) onto the current line (indicated by line_flag and
   * coefficients[]). Method should only be called if the current points define a line (line_flag is
   * not Flag.NO_LINE)
   *
   * @param x x component of the point
   * @param x y component of the point
   * @param x z component of the point
   */
  private boolean try_points_linear(double x, double y, double z) {

    switch(line_flag) {
    case XY_CONSTANT:
      return Global.fuzzy.eq(coefficients[0], x) && Global.fuzzy.eq(coefficients[1], y);

    case XZ_CONSTANT:
      return Global.fuzzy.eq(coefficients[0], x) && Global.fuzzy.eq(coefficients[1], z);

    case YZ_CONSTANT:
      return Global.fuzzy.eq(coefficients[0], y) && Global.fuzzy.eq(coefficients[1], z);

    case X_CONSTANT:
      return Global.fuzzy.eq(coefficients[0], x) && Global.fuzzy.eq(z, coefficients[1] * y + coefficients[2]);

    case Y_CONSTANT:
      return Global.fuzzy.eq(coefficients[0], y) && Global.fuzzy.eq(z, coefficients[1] * x + coefficients[2]);

    case Z_CONSTANT:
      return Global.fuzzy.eq(coefficients[0], z) && Global.fuzzy.eq(y, coefficients[1] * x + coefficients[2]);

    case NO_CONSTANT_VARIABLES:
      return Global.fuzzy.eq(y, coefficients[0] * x + coefficients[1]) && Global.fuzzy.eq(z, coefficients[2] * x + coefficients[3]);

    // error if the code gets here because that means, we are trying to fit the points to a linear
    // when the points already form a plane
    default:
      throw new RuntimeException("try_points_linear called when points already form a plane (line_flag = NO_LINE)");
    }
  }

  /**
   * Returns the two points that have the maximum separation in pa.
   *
   * @param pa array of points. Must have at least 2 elements. Can be any length and can contain
   *      nulls (which will be ignored).
   * @return a 2-element array containing the most-separated elements. Returns null if no pair of
   *      elements is NaN-free.
   */
  private Point @Nullable [] maxsep_doubles(@Nullable Point[] pa) {

    Point p1 = null;
    Point p2 = null;

    // max_separation is the separation metric for the most separated
    // pair of points.  We use the sum of the separations as the
    // metric.
    double max_separation = Double.MIN_VALUE;

    for (int i = 0; i < pa.length - 1; i++) {
      for (int j = i + 1; j < pa.length; j++) {
        double separation = separation(pa[i], pa[j]);
        if (separation > max_separation) {
          max_separation = separation;
          p1 = pa[i];
          p2 = pa[j];
        }
      }
    }

    // This can happen if enough values are NaN.
    // if (p1 == null) {
    //   throw new IllegalArgumentException(Arrays.toString(pa));
    // }

    if (Debug.logDetail()) {
      wrapper.log("maxsep_doubles = %s %s", p1, p2);
    }
    if (p1 == null) {
      return null;
    }
    return new Point[] { p1, p2 };
  }

  /**
   * Calculates the coefficients for the line in 3D
   * and sets the coefficients[] by side effect.
   *
   * @param p0 first point on the line
   * @param p1 second point on the line
   */
  //
  //  the line can fall under 3 categories
  //
  //  1. parallel to an axis which means that 2 of the following is true
  //  x = A
  //  y = B
  //  z = C
  //
  //  and the two that are true form a line in the 3rd dimension (the one that's
  //  not constant)
  //
  //  2. parallel to the plane of an axis, which means that one of the following
  //  is true
  //  x = A
  //  y = B
  //  z = C
  //
  //  and the line is just a line in the non-constacoefficientble plane projected
  //  onto the third dimension
  //  and the equation for that line is:
  //  x = A [suppose that x is the constant one]
  //  y = Bz + C
  //
  //
  //  3. none of the above and the equation of the line is just:
  //  y = Ax + B
  //  z = Cx + D
  //
  //  Note that A, B, C, and D are coefficients that are stored in coefficients[]

  private void calc_bi_linear(Point p0, Point p1) {

    if ((p0.x == p1.x) && (p0.y == p1.y) && (p0.z != p1.z)) {
      // x and y have constant values
      // line is defined by x = A, y = B
      line_flag = Flag.XY_CONSTANT;
      coefficients[0] = p0.x;
      coefficients[1] = p0.y;
      coefficients[2] = 0;
      coefficients[3] = 0;
      return;
    }
    if ((p0.x == p1.x) && (p0.z == p1.z) && (p0.y != p1.y)) {
      // x and z have constant values
      // line is defined by x = A, z = C
      line_flag = Flag.XZ_CONSTANT;
      coefficients[0] = p0.x;
      coefficients[1] = p0.z;
      coefficients[2] = 0;
      coefficients[3] = 0;
      return;
    }
    if ((p0.y == p1.y) && (p0.z == p1.z) && (p0.x != p1.x)) {
      // y and z have constant values
      // line is defined by y = B, z = C
      line_flag = Flag.YZ_CONSTANT;
      coefficients[0] = p0.y;
      coefficients[1] = p0.z;
      coefficients[2] = 0;
      coefficients[3] = 0;
      return;
    }
    if ((p0.x == p1.x) && (p0.y != p1.y) && (p0.z != p1.z)) {
      // x has a constant value
      // plane (degenerate) is defined by x = A, z = By + C
      line_flag = Flag.X_CONSTANT;
      coefficients[0] = p0.x;
      coefficients[1] = (p1.z - p0.z)/(p1.y - p0.y);
      coefficients[2] = (p0.z * p1.y - p0.y * p1.z)/(p1.y - p0.y);
      coefficients[3] = 0;
      return;
    }
    if ((p0.y == p1.y) && (p0.x != p1.x) && (p0.z != p1.z)) {
      // y has a constant value
      // plane (degenerate) is defined by y = B, z = Ax + C
      line_flag = Flag.Y_CONSTANT;
      coefficients[0] = p0.y;
      coefficients[1] = (p1.z - p0.z)/(p1.x - p0.x);
      coefficients[2] = (p0.z * p1.x - p0.x * p1.z)/(p1.x - p0.x);
      coefficients[3] = 0;
      return;
    }
    if ((p0.z == p1.z) && (p0.x != p1.x) && (p0.y != p1.y)) {
      // z has a constant value
      // plane (degenerate) is defined by z = C, y = Ax + B
      line_flag = Flag.Z_CONSTANT;
      coefficients[0] = p0.z;
      coefficients[1] = (p1.y - p0.y)/(p1.x - p0.x);
      coefficients[2] = (p0.y * p1.x - p0.x * p1.y)/(p1.x - p0.x);
      coefficients[3] = 0;
      return;
    }
    if ((p0.x != p1.x) && (p0.y != p1.y) && (p0.z != p1.z)) {
      // no variables have a constant value
      // plane (degenerate) is defined by y = Ax + B, z = Cx + D
      line_flag = Flag.NO_CONSTANT_VARIABLES;
      coefficients[0] = (p1.y - p0.y)/(p1.x - p0.x);
      coefficients[1] = (p0.y * p1.x - p0.x * p1.y)/(p1.x - p0.x);
      coefficients[2] = (p1.z - p0.z)/(p1.x - p0.x);
      coefficients[3] = (p0.z * p1.x - p0.x * p1.z)/(p1.x - p0.x);
    }
  }

  /**
   *  Calculates new coefficients that for the new point. Uses the
   *  new coefficients if they are relatively close to to the previous
   *  ones. Kills off the invariant if they are not.
   *
   *  @return true if the new equation worked, false otherwise
   */
  public boolean try_new_equation(double x, double y, double z) {

    // Calculate max separation using this point and the existing 3 points.
    def_points[3] = new Point(x, y, z);
    double sep = maxsep_triples(def_points);

    // If this point increased the separation, recalculate a, b, and c.
    if (sep > separation) {
      separation = sep;
      double[] coef;
      try {
        // Requires that the first 3 elements of def_points must be non-null.
        coef = calc_tri_linear(def_points);
        if (Debug.logDetail()) {
          wrapper.log("Calc new plane with points %s %s %s %s",
                       def_points[0], def_points[1], def_points[2], def_points[3]);
        }
      } catch (Exception e) {
        return false;
      }

      // if the a, b, or c is a new min/max remember it.
      if (coef[0] < min_a) min_a = coef[0];
      if (coef[0] > max_a) max_a = coef[0];
      if (coef[1] < min_b) min_b = coef[1];
      if (coef[1] > max_b) max_b = coef[1];
      if (coef[2] < min_c) min_c = coef[2];
      if (coef[2] > max_c) max_c = coef[2];
      if (coef[3] < min_d) min_d = coef[3];
      if (coef[3] > max_d) max_d = coef[3];

      // Pick a new a, b, and c as the average of their endpoints
      a = (min_a + max_a) / 2;
      b = (min_b + max_b) / 2;
      c = (min_c + max_c) / 2;
      d = (min_d + max_d) / 2;
      if (Invariant.logOn() || debug.isLoggable(Level.FINE)) {
        wrapper.log(debug, wrapper.ppt.name() + ": Trying new a (" + a
                     + "), b (" + b + "), c (" + c + "), and d (" + d + ")");
      }
      wrapper.log("min max: Trying new a (%s), b (%s), c (%s), and d (%s)",
                   a, b, c, d);
      // if the new coefficients are 'equal' to their min and max and
      // this point fits, then this new equation is good enough both
      // for existing points and the new point.
      if (Global.fuzzy.eq(a, min_a) && Global.fuzzy.eq(a, max_a)
          && Global.fuzzy.eq(b, min_b) && Global.fuzzy.eq(b, max_b)
          && Global.fuzzy.eq(c, min_c) && Global.fuzzy.eq(c, max_c)
          && Global.fuzzy.eq(d, min_d) && Global.fuzzy.eq(d, max_d)
          && Global.fuzzy.eq(-d, a * x + b * y + c * z)) {
        if (debug.isLoggable(Level.FINE)) {
          debug.fine(wrapper.ppt.name() + ": New a (" + a + ") and b ("
                      + b + ") and c (" + c + ")");
        }
        return true;
      } else {
        return false;
      }
    } else { // this point doesn't increase the separation

      return false;
    }
  }

  /**
   * Calculates the separation between p1 and p2.
   *
   * @param p1 first point
   * @param p2 second point
   *
   * @return the distance between p1 and p2, or 0 if either point is null, or NaN if either point
   * contains a NaN
   */
  double separation(Point p1, Point p2) {

      // Potential problem:  This may return 0, for two reasons.  First, x1-x2
      // might be zero due to underflow, even if x1!=x2.  Second, squaring a
      // small x1-x2 might result in 0, even if x1-x2!=0.

    // make sure both points are specified
    if ((p1 == null) || (p2 == null)) {
      return 0;
    }

    double xsep = (p1.x - p2.x);
    double ysep = (p1.y - p2.y);
    double zsep = (p1.z - p2.z);
    return Math.sqrt(xsep * xsep + ysep * ysep + zsep * zsep);
  }

  /**
   * Calculates the three points that have the maximum separation in pa and places them as the first
   * three elements of pa.
   *
   * @param pa array of points. Must have at least 3 elements. Can be any length and can contain
   *     nulls (which will be ignored). Is side-effected so that the first three elements contain
   *     the points with the maximum total separation; this may introduce duplicates into the array.
   *     The first three elements are null if no triple of elements is NaN-free.
   * @return the maximum separation found
   */
  double maxsep_triples(@Nullable Point[] pa) {

    // cache values for the (square of the) distance between each pair of
    // points, to avoid duplicating work
    double[][] separations = new double[pa.length][pa.length];
    for (int i = 0; i < pa.length - 1; i++) {
      for (int j = i + 1; j < pa.length; j++) {
        separations[i][j] = separation(pa[i], pa[j]);
      }
    }

    Point p1 = null;
    Point p2 = null;
    Point p3 = null;

    // max_separation is the separation metric for the most separated
    // triple of points.  We use the sum of the separations as the
    // metric.  (The metric "min of the three separations" does not work
    // because it doesn't choose a unique set of points, making the
    // result dependent on the order in which the points are seen; more
    // seriously, it may choose a collinear set of points even when a
    // non-collinear set exists.)
    double max_separation = Double.MIN_VALUE;

    for (int i = 0; i < pa.length - 2; i++) {
      for (int j = i + 1; j < pa.length - 1; j++) {
        double sep_i_j = separations[i][j];
        for (int k = j + 1; k < pa.length; k++) {

          // using the sum of separations is valid, as long as
          // separation is the actual distance between points and not
          // the square of the distance
          double separation = sep_i_j + separations[i][k] + separations[j][k];
          if (separation > max_separation) {
            max_separation = separation;
            p1 = pa[i];
            p2 = pa[j];
            p3 = pa[k];
          }
        }
      }
    }

    pa[0] = p1;
    pa[1] = p2;
    pa[2] = p3;

    if (Debug.logDetail()) {
      wrapper.log("maxsep_triples = %s %s %s", pa[0], pa[1], pa[2]);
    }
    return max_separation;
  }

  // Given ((x0,y0,z0),(x1,y1,z1), (x2,y2,z2), calculate a, b, c and d
  // such that ax + by + cz + d = 0.
  // a, b, c and d are mutually prime, integers and a is positive.
  //
  // A visual explanation of the math can be found in:
  // http://www.efm.leeds.ac.uk/CIVE/CIVE2599/cive2599-summary-overheads-full.pdf
  //
  // The standard form of a plane can be calculated using a vector
  // normal "n" <n1, n2, n3> to the plane and a vector from the origin
  // to the plane "r" e.g. <x0, y0, z0>.
  //
  // The vector normal to the plane can be calculated by doing the cross product
  // of two line segments on the plane (<point2 - point0> and <point1 - point0>)
  //
  // Given the normal vector and a point on the plane, the standard form
  // of the plane is:
  // n1*x + n2*y + n3*z = p where p is the dot product of r and n
  //
  //

  /**
   * Calculates the coefficients (a, b, c) and constant (d) for the equation ax + by + cz + d = 0
   * using the first three points in points.
   *
   * @param points array of points to use to calculate the coefficents. Only the first three points
   *     (all of which must be non-null) are used.
   *
   * @return a four element array where a is the first element, b the second, c the third, and d is
   * the fourth. All elements are mutually prime, integers and a is positive.
   */
  // TODO: should just pass in the first three elements rather than passing
  // in the point array.
  public double[] calc_tri_linear(@Nullable Point[] points) { // TODO: remove annotation after refactoring method

    Point p0 = points[0];
    Point p1 = points[1];
    Point p2 = points[2];

    // the following code is taken from the source page of
    // http://jblanco_60.tripod.com/plane.html (the math behind the
    // implementation is detailed above)

    double px = p0.x;
    double py = p0.y;
    double pz = p0.z;
    double qx = p1.x;
    double qy = p1.y;
    double qz = p1.z;
    double rx = p2.x;
    double ry = p2.y;
    double rz = p2.z;

    // define the i, j, k components of the two line segments on the plane
    double a1 = (qx - px);
    double a2 = (qy - py);
    double a3 = (qz - pz);

    double b1 = (rx - px);
    double b2 = (ry - py);
    double b3 = (rz - pz);

    // calculate the i, j, k components of the cross product
    double i_var = (a2 * b3) - (b2 * a3);
    double j_var = -((a1 * b3) - (b1 * a3));
    double k_var = (a1 * b2) - (b1 * a2);

    // Calculate the value of the constant.  Note that given the format
    // of the point-normal form, we multiply by the negative of P(x,y,z)

    double Stand_const = ((i_var*(-px)) + (j_var*(-py)) + (k_var*(-pz)));

    // Put the coefficients in lowest terms by dividing by the gcd.
    // If the gcd is 1, no harm done -- divide by 1.
    // Type is double to ensure floating-point division.
    double myGCD = MathPlume.gcd(MathPlume.gcd(i_var, j_var) , MathPlume.gcd(k_var, Stand_const));
    double standard_i_var = i_var / myGCD;
    double standard_j_var = j_var / myGCD;
    double standard_k_var = k_var / myGCD;
    double standard_Stand_const = Stand_const / myGCD;

    if (wrapper != null) {
      wrapper.log("Preliminary: %s %s GCD: %s", i_var, j_var, MathPlume.gcd(i_var, j_var));
      wrapper.log("Preliminary: %s %s GCD: %s", k_var, Stand_const, MathPlume.gcd(k_var, Stand_const));
      wrapper.log("GCD: %s", myGCD);
    }

    // Ensure postive x coefficient by negating if "x" term is negative.
    if (i_var < 0) {
      standard_i_var = - standard_i_var;
      standard_j_var = - standard_j_var;
      standard_k_var = - standard_k_var;
      standard_Stand_const = - standard_Stand_const;
    }

    double[] coef = new double[] {standard_i_var,
                                  standard_j_var,
                                  standard_k_var, standard_Stand_const };

    return coef;
  }

  public boolean enoughSamples(@GuardSatisfied LinearTernaryCoreFloat this) {
    return isActive();
  }

  // If the values form a line,
  // we have no confidence that it's a plane.
  // The line is less interesting because the invariant
  // is already defined by a linear binary invariant.
  public double computeConfidence() {
    if (isActive()) {
      return Invariant.conf_is_ge(values_seen, MINTRIPLES);
    } else {
      return 0;
    }

  }

  public String repr(@GuardSatisfied LinearTernaryCoreFloat this) {
    return "LinearTernaryCoreFloat" + wrapper.varNames() + ": a=" + a
      + ",b=" + b
      + ",c=" + c
      + ",d=" + d
      + ",values_seen=" + values_seen;
  }

  public String point_repr(Point p) {
    if (p == null) {
      return "null";
    } else {
      return "<" + p.x + "," + p.y + "," + p.z + ">";
    }
  }

  public String cache_repr() {
    StringJoiner result = new StringJoiner("; ");
    for (int i = 0; i < MINTRIPLES; i++) {
      result.add(point_repr(def_points[i]));
    }
    return result.toString();
  }

  // In this class for convenience (avoid prefixing "LinearBinaryCore").
  static String formatTerm(double coeff, @Nullable String varname, boolean first) {
    return LinearBinaryCore.formatTerm(coeff, varname, first);
  }

  @SideEffectFree
  public String format_using(OutputFormat format,
                             String vix, String viy, String viz,
                             double a, double b, double c, double d) {

    if ((a == 0) && (b == 0) && (c == 0) && (d == 0)) {
      return wrapper.format_too_few_samples(format, null);
    }

    if (format == OutputFormat.SIMPLIFY) {
      return format_simplify(vix, viy, viz, a, b, c, d);
    }

    if ((format.isJavaFamily())) {

        return
          // Add a 1 to each side of the equation because FuzzyFloat
          // doesn't treat 0 as 'fuzzy' when compared against another
          // number.
          "daikon.Quant.fuzzy.eq("
          + formatTerm(a, vix, true)
          + formatTerm(b, viy, (a == 0))
          + formatTerm(c, viz, (a == 0) && (b == 0))
          + formatTerm(d, null, (a == 0) && (b == 0) && (c == 0))
          + " + 1, 1)";

      }

    if ((format == OutputFormat.DAIKON)
        || (format == OutputFormat.ESCJAVA)
        || (format == OutputFormat.CSHARPCONTRACT)) {
        String eq = " == ";
        return formatTerm(a, vix, true)
          + formatTerm(b, viy, (a == 0))
          + formatTerm(c, viz, (a == 0) && (b == 0))
          + formatTerm(d, null, (a == 0) && (b == 0) && (c == 0))
          + eq + "0";
      }

    throw new Error("unrecognized output format " + format);
    // return null;
  }

  public static String format_simplify(String vix, String viy, String viz,
                                       double da, double db, double dc, double dd) {
    int ia = (int) da;
    int ib = (int) db;
    int ic = (int) dc;
    int id = (int) dd;

    String str_a, str_b, str_c, str_d;
    if (ia != da || ib != db || ic != dc || id != dd) {
      // floating point

      // Disabled for the moment, since mixing integers and floating
      // literals seems to give Simplify indigestion. For instance:
      // (BG_PUSH (<= w 3))
      // (BG_PUSH (EQ 0 (* w 2.0d0)))
      // (<= w 1)

      // str_a = Invariant.simplify_format_double(da);
      // str_b = Invariant.simplify_format_double(db);
      // str_c = Invariant.simplify_format_double(dc);
      // str_d = Invariant.simplify_format_double(dd);
      return "(AND)"; // really s/b format_unimplemented
    } else {
      // integer
      str_a = Invariant.simplify_format_long(ia);
      str_b = Invariant.simplify_format_long(ib);
      str_c = Invariant.simplify_format_long(ic);
      str_d = Invariant.simplify_format_long(id);
    }

    String str_z = viz;
    String str_x = vix;
    String str_y = viy;
    String str_ax = (da == 1.0) ? str_x : "(* " + str_a + " " + str_x + ")";
    String str_by = (db == 1.0) ? str_y : "(* " + str_b + " " + str_y + ")";
    String str_cz = (dc == 1.0) ? str_z : "(* " + str_c + " " + str_z + ")";
    String str_axPbyPcz = "(+ " + str_ax + " " + str_by + " " + str_cz + ")";
    String str_axPbyPczPd = (dd == 0.0) ? str_axPbyPcz :
      "(+ " + str_axPbyPcz + " " + str_d + ")";
    return "(EQ 0 " + str_axPbyPczPd + ")";
  }

  @SideEffectFree
  public String format_using(OutputFormat format,
                             String vix, String viy, String viz
                             ) {
    String result = format_using(format, vix, viy, viz, a, b, c, d);
    if (result != null) {
      return result;
    }

    return wrapper.format_unimplemented(format);
  }

  // // Format as "x = cy+d" instead of as "y = ax+b".
  // public String format_reversed(String x, String y) {
  //   assert a == 1 || a == -1;
  //   return format(y, x, a, -b/a);
  // }

  @Pure
  public boolean isSameFormula(LinearTernaryCoreFloat other) {

    // If both have yet to see enough values
    if (!isActive() && !other.isActive()) {

      // Same formula if all of the points match
      if (values_seen != other.values_seen) {
        return false;
      }
      // This and elsewhere seem to assume that values_seen < MINTRIPLES;
      // is that necessarily true?
      for (int ii = 0; ii < values_seen; ii++) {
        // used to use !=, but that test seems wrong
        if (!def_points[ii].equals(other.def_points[ii])) {
          return false;
        }
      }
      return true;
    } else {
      return ((values_seen >= MINTRIPLES)
              && (other.values_seen >= MINTRIPLES)
              && (a == other.a)
              && (b == other.b)
              && (c == other.c)
              && (d == other.d));
    }
  }

  @Pure
  public boolean isExclusiveFormula(LinearTernaryCoreFloat other) {
    if (!isActive()
        || !other.isActive()) {
      return false;
    }

    return ((a == other.a)
            && (b != other.b)
            && (c != other.c)
            && (d != other.d));
  }

  /**
   * In general, we can't merge formulas, but we can merge invariants with too few samples to have
   * formed a plane with invariants with enough samples. And those will appear to have different
   * formulas.
   */
  public boolean mergeFormulasOk() {
    return true;
  }

  /**
   * Merge the linear ternary cores in cores to form a new core. Any core in the list that has seen
   * enough points to define a plane, must define the same plane. Any cores that have not yet seen
   * enough points, will have each of their points applied to that invariant. The merged core is
   * returned. Null is returned if the cores don't describe the same plane
   *
   * @param cores list of LinearTernary cores to merge. They should all be permuted to match the
   *     variable order in ppt.
   */
  public @Nullable LinearTernaryCoreFloat merge(List<LinearTernaryCoreFloat> cores, Invariant wrapper) {

    // Look for any active planes.  All must define the same plane
    LinearTernaryCoreFloat first = null;
    for (LinearTernaryCoreFloat c : cores) {
      if (!c.isActive()) {
        continue;
      }
      if (first == null) {
        first = c.clone();
      } else {
        if (!Global.fuzzy.eq(first.a, c.a)
            || !Global.fuzzy.eq(first.b, c.b)
            || !Global.fuzzy.eq(first.c, c.c)
            || !Global.fuzzy.eq(first.d, c.d))
          return null;
      }
    }

    // If no active planes were found, created an empty core
    if (first == null) {
      first = new LinearTernaryCoreFloat(wrapper);
    } else {
      first.wrapper = wrapper;
    }

    // Merge in any points from non-active cores
    // Note that while null points can't show up in a core that is not active
    // because it hasn't seen enough values, they can show up if the core
    // is not active because it is a line.
    for (LinearTernaryCoreFloat c : cores) {
      if (c.isActive()) {
        continue;
      }
      for (int j = 0; j < c.values_seen; j++) {
        Point cp = c.def_points[j];
        if (cp == null) {
          continue;
        }
        if (Debug.logDetail()) {
          wrapper.log("Adding point %s from %s", cp, c.wrapper.ppt);
        }
        InvariantStatus stat = first.add_modified(cp.x, cp.y, cp.z, 1);
        if (stat == InvariantStatus.FALSIFIED) {
          return null;
        }
        // if (wrapper.is_false())
        //  return null;
      }
    }

    return first;
  }
}