001// ***** This file is automatically generated from LinearTernaryCore.java.jpp
002
003package daikon.inv.ternary.threeScalar;
004
005import daikon.*;
006import daikon.inv.*;
007import daikon.inv.binary.twoScalar.*;
008import daikon.inv.unary.scalar.*;
009import java.io.Serializable;
010import java.util.*;
011import java.util.logging.Level;
012import java.util.logging.Logger;
013import org.checkerframework.checker.lock.qual.GuardSatisfied;
014import org.checkerframework.checker.nullness.qual.Nullable;
015import org.checkerframework.dataflow.qual.Pure;
016import org.checkerframework.dataflow.qual.SideEffectFree;
017import org.plumelib.util.ArraysPlume;
018import org.plumelib.util.MathPlume;
019import typequals.prototype.qual.NonPrototype;
020import typequals.prototype.qual.Prototype;
021
022/**
023 * The LinearTernaryCore class is acts as the backend for the invariant (ax + by + cz + d = 0) by
024 * processing samples and computing coefficients. The resulting coefficients a, b, c, and d are
025 * mutually relatively prime, and the coefficient a is always p.
026 */
027
028// Originally, LinearTernaryCore code was not dealing with degenerate
029// linear ternary invariants correct, namely when the plane
030// (ax+by+cz+d=0) is parallel to one of the axes, so one of the
031// coefficients was 0.  In this case, the invariant can be described
032// by a binary invariant (eliminate the 0 coefficient variable).
033// LinearTernaryCore's inability to deal with this degenerate form was
034// (and still is) masked by the suppression of a linear ternary
035// invariant via a binary invariant and a constant (the code threw an
036// exception when seeing this case and falsifies the invariant).  As
037// more samples are processed, the ternary invariant may "evolve" from
038// a degenerate form to a full ternary invariant.  The
039// linearintervention method deals with the degenerate form and the
040// planarintervention focuses on the full ternary invariant.
041
042@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
043    "UnnecessaryParentheses"  // generated code, parens are sometimes necessary
044    })
045public final class LinearTernaryCore implements Serializable, Cloneable {
046  // We are Serializable, so we specify a version to allow changes to
047  // method signatures without breaking serialization.  If you add or
048  // remove fields, you should change this number to the current date.
049  static final long serialVersionUID = 20030822L;
050
051  /** Debug tracer. */
052  static final Logger debug = Logger.getLogger("daikon.inv.ternary.threeScalar.LinearTernaryCore");
053
054  // ax + by + cz + d = 0; first argument is x, second is y, third is z
055  // Although the coefficients are guaranteed to be integers, we use doubles to
056  // store the values because regression tests have shown that a long is
057  // not sufficient
058  public double a = 0, b = 0, c = 0, d = 0;
059  public double min_a, max_a, min_b, max_b, min_c, max_c, min_d, max_d;
060  public double separation = 0;
061
062  // coefficients for the invariant in 3D space
063  public double[] coefficients = new double[4];
064
065  // The values of the flags determine indicate whether a line exists
066  // with the current samples (NO_LINE) and which variables are
067  // constant (constant variables indicate that the 'plane' is
068  // parallel to that axis).  For example, XY_CONSTANT means that the
069  // so far, the x and y values have all been constant and only z has
070  // been varying so so the plane is actually a line.
071  // NO_CONSTANT_VARIABLES indicates that the values have not been
072  // constant but do for a plane between variables.
073  public enum Flag {NO_LINE, XY_CONSTANT, XZ_CONSTANT, YZ_CONSTANT,
074                    X_CONSTANT, Y_CONSTANT, Z_CONSTANT, NO_CONSTANT_VARIABLES};
075
076  // The flag indicates whether a line exists in the 3D space instead
077  // of a plane If the flag value is not NO_LINE, then a line does
078  // currently exist among the values seen and the flag specifies
079  // which variables are constant to make the invariant a line rather
080  // than a plane.  The value NO_LINE indicates either that not enough
081  // values have been seen to make any decision or that a plane
082  // exists.
083  public Flag line_flag = Flag.NO_LINE;
084
085  // Inner class to store points in for more convenient access.
086  public static class Point implements Serializable, Cloneable {
087    static final long serialVersionUID = 20050923L;
088
089    public long x;
090    public long y;
091    public long z;
092    public Point() {
093      this(0, 0, 0);
094    }
095    public Point(long x, long y, long z) {
096      this.x = x; this.y = y; this.z = z;
097    }
098    @Pure
099    public boolean equals(long x, long y, long z) {
100      return (this.x == x) && (this.y == y) && (this.z == z);
101    }
102    @SideEffectFree
103    @Override
104    protected Point clone(@GuardSatisfied Point this) throws CloneNotSupportedException {
105      return (Point) super.clone();
106    }
107    @SideEffectFree
108    @Override
109    public String toString(@GuardSatisfied Point this) {
110      return "(" + x + ", " + y + " ," + z + ")";
111    }
112  }
113
114  // Points that define the plane.  The first three are normally the
115  // current definition.  The next is a new point for consideration.
116  // However, any slot may be set to null at any time, and points may be
117  // duplicated at any time.
118  // (The invariant appears to be a bit stronger than that, since there are
119  // places that dereference the array unconditionally, assuming that its
120  // elements are non-null.  This should be clarified/refactored.)
121  public @Nullable Point[] def_points = new Point[MINTRIPLES];
122
123  public Invariant wrapper;
124
125  // The number of distinct values (not samples) seen.
126  public int values_seen = 0;
127
128  static final int MINTRIPLES = 5;
129
130  public LinearTernaryCore(Invariant wrapper) {
131    this.wrapper = wrapper;
132  }
133
134  @SideEffectFree
135  @Override
136  public LinearTernaryCore clone(@GuardSatisfied LinearTernaryCore this) {
137    try {
138      LinearTernaryCore result = (LinearTernaryCore) super.clone();
139      result.def_points = new Point[MINTRIPLES];
140      for (int i = 0; i < MINTRIPLES; i++) {
141        Point p = def_points[i];
142        if (p != null) {
143          result.def_points[i] = p.clone();
144        }
145      }
146      return result;
147    } catch (CloneNotSupportedException e) {
148      throw new Error(); // can't happen
149    }
150  }
151
152  /**
153   * Reorganize our already-seen state as if the variables had shifted order underneath us
154   * (rearrangement given by the permutation).
155   */
156  public void permute(int[] permutation) {
157    assert permutation.length == 3;
158    assert ArraysPlume.fnIsPermutation(permutation);
159    // Fix a, b, c
160    double[] clever = new double[] { a, b, c };
161    double[] pclever = new double[3];
162    pclever[permutation[0]] = clever[0];
163    pclever[permutation[1]] = clever[1];
164    pclever[permutation[2]] = clever[2];
165    if (isActive()) {
166      // Seen enough values, so permuting a, b, c is useful
167      if ( a == 0 || b == 0 || c == 0) {
168        throw new Error("Active LTs never have 0 coefficients");
169        // // Can't handle this form once rotated.  But if we've seen
170        // // enough values, yet a, b or c is zero, then this is covered by
171        // // a LinearBinary or a constant, so let me just destroy myself.
172        // values_seen = Integer.MAX_VALUE;
173        // a = b = c = d = 0;
174        // if (debug.isLoggable(Level.FINE)) {
175        //   debug.fine ("  Ternary invariant destroyed because a, b, or c = 0");
176        // }
177        // return;
178      } else {
179        // double d = -1.0 / pclever[2];
180        a = pclever[0];
181        b = pclever[1];
182        c = pclever[2];
183
184        // we still need to guarantee a positive x coefficient
185        if (a < 0) {
186          a = -a;
187          b = -b;
188          c = -c;
189          d = -d;
190        }
191      }
192    }
193
194    // Fix caches
195    {
196      long[] temp = new long[3];
197      for (int i = 0; i < MINTRIPLES; i++) {
198        Point p = def_points[i];
199        if (p == null) {
200          continue;
201        }
202        wrapper.log("orig def_points[%s] = %s", i, p);
203        temp[permutation[0]] = p.x;
204        temp[permutation[1]] = p.y;
205        temp[permutation[2]] = p.z;
206        p.x = temp[0];
207        p.y = temp[1];
208        p.z = temp[2];
209        wrapper.log("permuted def_points[%s] = %s", i, p);
210      }
211    }
212    // Assert that caches sync with a,b,c?
213    // This would be a good sanity check, but it would be nontrivial
214    // because we don't keep track of when a, b, and c are really
215    // valid for all the of the points we've cached. If we were
216    // falsified and then resurrected, we might have samples that
217    // didn't fit even before the permutation. -smcc
218  }
219
220  /**
221   * Returns whether or not the invariant is currently active. We become active after MINTRIPLES
222   * values have been seen and a line is not calculated. In addition, the coefficients (a, b, c)
223   * must all be non-zero, else the invariant would be described by a LinearBinary or constant
224   * invariant. Before that, a, b, and c are uninitialized.
225   */
226  // this is the check used throughout the file as a shortcut because
227  // the invariant must have seen enough samples and not be a line for all the
228  // other calculations and comparisons done to be sensible, meaning that
229  // a full linear ternary invariant exists, not a degenerate one
230  @Pure
231  public boolean isActive(@GuardSatisfied LinearTernaryCore this) {
232    return line_flag == Flag.NO_LINE && values_seen >= MINTRIPLES && a != 0 && b != 0 && c != 0;
233  }
234
235  /**
236   * Sets up the invariant from a LinearBinary invariant and a constant value for the third
237   * variable. Points are taken from the LinearBinary cache and its min_x/y and max_x/y points and
238   * combined with the constant value.
239   *
240   * @return InvariantStatus.NO_CHANGE if the invariant is valid, InvariantStatus.FALSIFIED if one
241   * of the points invalidated the LinearTernary invariant
242   */
243  public InvariantStatus setup(LinearBinary lb, VarInfo con_var,
244                                long con_val) {
245
246    if (Debug.logOn()) {
247      wrapper.log("setup from lb %s con var %s con_val %s",
248                   lb, con_var, con_val);
249    }
250
251    int con_index = con_var.varinfo_index;
252    int lb_v1_index = lb.ppt.var_infos[0].varinfo_index;
253    int lb_v2_index = lb.ppt.var_infos[1].varinfo_index;
254
255    long[] con_vals = new long [lb.core.values_seen + 2];
256    long[] lb1_vals = new long [lb.core.values_seen + 2];
257    long[] lb2_vals = new long [lb.core.values_seen + 2];
258
259    int mi = lb.core.values_seen;
260    for (int i = 0; i < mi; i++) {
261      con_vals[i] = con_val;
262      lb1_vals[i] = lb.core.x_cache[i];
263      lb2_vals[i] = lb.core.y_cache[i];
264    }
265    if (lb.isActive()) {
266      con_vals[mi] = con_val;
267      lb1_vals[mi] = lb.core.min_x;
268      lb2_vals[mi] = lb.core.min_y;
269      con_vals[mi + 1] = con_val;
270      lb1_vals[mi + 1] = lb.core.max_x;
271      lb2_vals[mi + 1] = lb.core.max_y;
272      mi += 2;
273    }
274
275    InvariantStatus sts = InvariantStatus.NO_CHANGE;
276
277    for (int i = 0; i < mi; i++ ) {
278      if (con_index < lb_v1_index) {
279        sts = add_modified(con_vals[i], lb1_vals[i], lb2_vals[i], 1);
280
281      } else if (con_index < lb_v2_index) {
282        sts = add_modified(lb1_vals[i], con_vals[i], lb2_vals[i], 1);
283
284      } else {
285        sts = add_modified(lb1_vals[i], lb2_vals[i], con_vals[i], 1);
286
287      }
288      if (sts != InvariantStatus.NO_CHANGE) {
289        break;
290      }
291    }
292
293    if (sts != InvariantStatus.NO_CHANGE) {
294      System.out.println("lb.core.values_seen=" + lb.core.values_seen);
295      for (int i = 0; i < mi; i++ ) {
296        System.out.printf("LTCore: vals %s %s %s%n", con_vals[i], lb1_vals[i],
297                           lb2_vals[i]);
298      }
299      System.out.println("in inv " + wrapper.format() + " " + wrapper.ppt);
300      assert sts == InvariantStatus.NO_CHANGE;
301    }
302
303    return sts;
304  }
305
306  /**
307   * Sets up the invariant from a OneOf and two constants. Points are taken from the OneOf cache
308   * and the constant values.
309   *
310   * @return InvariantStatus.NO_CHANGE if the invariant is valid, or InvariantStatus.FALSIFIED if
311   * one of the points invalidated the LinearTernary invariant
312   */
313
314  public InvariantStatus setup(OneOfScalar oo, VarInfo v1, long con1,
315                                VarInfo v2, long con2) {
316
317    int oo_index = oo.ppt.var_infos[0].varinfo_index;
318    int con1_index = v1.varinfo_index;
319    int con2_index = v2.varinfo_index;
320
321    InvariantStatus sts = InvariantStatus.NO_CHANGE;
322    for (int i = 0; i < oo.num_elts(); i++ ) {
323      if (oo_index < con1_index) {
324        sts = add_modified(((Long)oo.elt(i)).longValue(), con1, con2, 1);
325      } else if (oo_index < con2_index) {
326        sts = add_modified(con1, ((Long)oo.elt(i)).longValue(), con2, 1);
327      } else {
328        sts = add_modified(con1, con2, ((Long)oo.elt(i)).longValue(), 1);
329      }
330      if (sts != InvariantStatus.NO_CHANGE) {
331        break;
332      }
333    }
334
335    if (Debug.logOn()) {
336      wrapper.log("setup from OneOf %s v1=%s v2=%s status = %s",
337                   oo, con1, con2, sts);
338    }
339
340    return sts;
341  }
342
343  /**
344   * Looks for points that define a plane (ax + by + cz + d = 0). Collects MINTRIPLE points before
345   * attempting to define a line through the points. Once the equation for the line is found, each
346   * subsequent point is compared to it. If the point does not match, recalculates the maximally
347   * separated points and attempts to fit the points to the new line. If the points do not fit the
348   * line, attempts to define the plane (to hopefully get at least some spread between the points,
349   * so that small errors don't get magnified). Once the equation for the plane is found, each
350   * subsequent point is compared to it. If the point does not match the point is examined to see if
351   * it would is maximally separated when compared to the points originally used to define the
352   * plane. If it is, it is used to recalcalulate the coefficients (a, b, c). If those
353   * coefficients are relatively close to the original coefficients (within the ratio defined by
354   * Global.fuzzy) then the new coefficients are used.
355   *
356   * @see org.plumelib.util.FuzzyFloat
357   */
358  public InvariantStatus add_modified(long x, long y, long z, int count) {
359
360    if (Debug.logDetail()) {
361      wrapper.log("Adding point, x=%s y=%s z=%s to invariant", x, y, z);
362    }
363
364    if (values_seen < MINTRIPLES) {
365      // We delay computation of a and b until we have seen several triples
366      // so that we can compute a and b based on a far-separated triple.  If
367      // the points in a triple are nearby, then roundoff errors in the
368      // computation of the slope can be non-negligible.
369
370      // skip points we've already seen
371      for (int i = 0; i < values_seen; i++) {
372        if (def_points[i].equals(x, y, z)) {
373          return InvariantStatus.NO_CHANGE;
374        }
375      }
376
377      def_points[values_seen] = new Point(x, y, z);
378      wrapper.log("Add: (%s, %s, %s)", x, y, z);
379      values_seen++;
380      wrapper.log("Values seen: %s", values_seen);
381
382      // check to see if we've seen enough values to create the equation
383      if (values_seen == MINTRIPLES) {
384
385        // Try to fit the points to a line first.
386        // The points may form a degenerate linear ternary invariants
387        // especially if one of the variable is constant
388        linearIntervention(def_points);
389
390        // if line can not be formed, try to form a plane
391        if (line_flag == Flag.NO_LINE) {
392          InvariantStatus stat = planarIntervention(def_points);
393          return stat;
394        } else {
395
396          // points fit in a line
397          return InvariantStatus.NO_CHANGE;
398        }
399      } else {
400
401        // still haven't seen enough values
402        return InvariantStatus.NO_CHANGE;
403      }
404    } else {
405
406      // at this point either a line or a plane must have been formed
407
408      // if line already broken, fit to plane equation
409      if (line_flag == Flag.NO_LINE) {
410        // If the new value doesn't fit the equation
411        if (!Global.fuzzy.eq(-d, a * x + b * y + c * z)) {
412
413          // Try to find small changes that will fit better.
414          if (!try_new_equation(x, y, z)) {
415            if (Invariant.logOn() || debug.isLoggable(Level.FINE)) {
416              wrapper.log(debug, "destroying  (" + wrapper.format()
417                           + ") where x=" + x + " y=" + y + " z=" + z
418                           + " a=" + a + " b=" + b + " c=" + c
419                           + " values_seen=" + values_seen);
420            }
421            // jiggling the values did not work
422            return InvariantStatus.FALSIFIED;
423          }
424        } else {
425
426          // point fits the current plane equation
427          return InvariantStatus.NO_CHANGE;
428        }
429      } else {
430
431        // try to fit to the current line
432        if (!try_points_linear((double)(x), (double)(y), (double)(z))) {
433          // put the crucial point in, so that it can be tested against the new line
434          def_points[MINTRIPLES - 1] = new Point(x, y, z);
435
436          // try to fit the points onto a line again (not sure if this
437          // is really necessary, maybe the fuzzy equals really does
438          // make a difference)
439          linearIntervention(def_points);
440        }
441
442        // new point has broken the line, try to fit to a plane
443        if (line_flag == Flag.NO_LINE) {
444          InvariantStatus stat = planarIntervention(def_points);
445          return stat;
446        } else {
447
448          // fuzzy equals does make a difference
449          return InvariantStatus.NO_CHANGE;
450        }
451
452      }
453
454    }
455    // should never get here because values_seen is either < or >= MINTRIPLES
456    // return InvariantStatus.NO_CHANGE;
457    throw new Error("at end of add_modified");
458  }
459
460  /**
461   * Attempts to fit the points in def_point to a plane and calculates the coefficients (a, b, c, d)
462   * if possible.
463   *
464   * @param def_point array of points. Must have at least 3 elements.
465   *
466   * @return the status of the invariant
467   *     (whether the plane fitting was successful)
468   */
469  private InvariantStatus planarIntervention(@Nullable Point[] def_point) {
470
471    // make a copy of the array so that the points don't get clobbered
472    // because we need to check all the points in def_points against the
473    // coefficients calculated
474
475    @Nullable Point[] dummy = new @Nullable Point[def_point.length];
476    for (int i = 0; i < def_point.length; i++) {
477      dummy[i] = def_point[i];
478    }
479
480    // Find the three points with the maximum separation
481    maxsep_triples(dummy);
482    // Now, first 3 values of dummy are the maximum-separated 3 points.
483
484    // Null points only show up in the cross-checker, not the regression
485    // tests, as of 5/25/2010.
486    // Sometimes the points are still null.
487    if (dummy[0] == null) {
488      return InvariantStatus.FALSIFIED;
489    }
490
491    // Calculate the coefficients of the equation (a, b, c, d).
492    // Requires that the first 3 elements of dummy must be non-null.
493    double[] coef = calc_tri_linear(dummy);
494
495    a = coef[0];
496    b = coef[1];
497    c = coef[2];
498    d = coef[3];
499
500    wrapper.log("Calculated tri linear coeff: (%s), b (%s), c(%s), and d(%s)",
501                 a, b, c, d);
502    // If one of these coefficients is zero (except for d), this
503    // should be a LinearBinary, not a LinearTernary, term.  (It
504    // might not show up as LinearBinary because there might not
505    // have been enough samples; but a random varying third
506    // variable can create enough samples.)  Also, throw out the
507    // invariant if any of the coefficients would be not-a-number.
508    if ( a == 0|| b == 0 || c == 0
509        || Double.isNaN(a) || Double.isNaN(b) || Double.isNaN(c)) {
510
511      wrapper.log("problematic coefficient: invariant falsified");
512      return InvariantStatus.FALSIFIED;
513    }
514
515    // Check all values against a, b, and c.
516    if (!wrapper.is_false()) {
517      for (int i = 0; i < values_seen; i++) {
518
519        // Global.fuzzy.eq works by comparing ratios between numbers,
520        // so comparing the sum of everything to 0 won't work.
521        //  wrapper.log("point: " + def_point[i].x + " "
522        //  + def_point[i].y + " " + def_point[i].z);
523        //  wrapper.log("compare: " + -c * def_point[i].z);
524        //  wrapper.log( "to: " + (a*def_point[i].x
525        //  + b*def_point[i].y + d));
526
527        wrapper.log("Trying at index %s: 0 != %s*%s+%s*%s+%s*%s+%s",
528                     i, a, def_point[i].x,
529                        b, def_point[i].y,
530                        c, def_point[i].z,
531                        d);
532
533        if (!Global.fuzzy.eq(-d, a * def_point[i].x + b * def_point[i].y + c * def_point[i].z)) {
534
535          if (Debug.logOn() || debug.isLoggable(Level.FINE)) {
536            wrapper.log(debug, "Destroying at index " + i + ": 0 != "
537                         + a + "*" + def_point[i].x
538                         + "+" + b + "*" + def_point[i].y + "+" + c
539                         + "*" + def_point[i].z
540                         + "+" + d);
541          }
542
543          return InvariantStatus.FALSIFIED;
544        }
545      }
546      if (Debug.logOn()) {
547        wrapper.log("equation = %s*x %s*y %s*z = %s", a, b, c, -d);
548      }
549
550      // Discard the points not used to define the coefficients.
551      for (int ii = 3; ii < MINTRIPLES; ii++) {
552        def_point[ii] = null;
553      }
554    }
555
556    // if it passes all the checks, then no change to the Invariant
557    return InvariantStatus.NO_CHANGE;
558  }
559
560  /**
561   * Attempts to fit the points in def_point to a line and calculates the corresponding line
562   * classification (line_flag) and coefficients (coefficients[]) if possible. If not (points in
563   * def_point form a plane) resets line_flag to Flag.NO_LINE and coefficients[] elements to -999.
564   *
565   * @param def_point array of points. Must have at least 2 elements.
566   *
567   */
568  private void linearIntervention(@Nullable Point[] def_point) {
569
570    // find the max separation points
571    Point[] maxsep_doubles = maxsep_doubles(def_point);
572
573    // maxsep_doubles[0] could be null when enough points have a NaN
574    if (maxsep_doubles == null) {
575      // line_flag == Flag.NO_LINE so the invariant
576      // will try to fit a plane to the points and fail
577      return;
578    } else {
579      // find the equation (coefficients now stored in coefficients[])
580      calc_bi_linear(maxsep_doubles[0], maxsep_doubles[1]);
581
582      // This seems to imply that all elements of def_points are non-null.
583      // Is that true?
584      // check the other points against the coefficients
585      for (int i = 0; i < def_point.length; i++) {
586        // try to fit the current point to the current line
587        boolean ok = try_points_linear((double)(def_point[i].x), (double)(def_point[i].y), (double)(def_point[i].z));
588
589        if (!ok) {
590          line_flag = Flag.NO_LINE;
591          break;
592        }
593      }
594    }
595
596  }
597
598  /**
599   * Attempts to fit the new point (x,y,z) onto the current line (indicated by line_flag and
600   * coefficients[]). Method should only be called if the current points define a line (line_flag is
601   * not Flag.NO_LINE)
602   *
603   * @param x x component of the point
604   * @param x y component of the point
605   * @param x z component of the point
606   */
607  private boolean try_points_linear(double x, double y, double z) {
608
609    switch(line_flag) {
610    case XY_CONSTANT:
611      return Global.fuzzy.eq(coefficients[0], x) && Global.fuzzy.eq(coefficients[1], y);
612
613    case XZ_CONSTANT:
614      return Global.fuzzy.eq(coefficients[0], x) && Global.fuzzy.eq(coefficients[1], z);
615
616    case YZ_CONSTANT:
617      return Global.fuzzy.eq(coefficients[0], y) && Global.fuzzy.eq(coefficients[1], z);
618
619    case X_CONSTANT:
620      return Global.fuzzy.eq(coefficients[0], x) && Global.fuzzy.eq(z, coefficients[1] * y + coefficients[2]);
621
622    case Y_CONSTANT:
623      return Global.fuzzy.eq(coefficients[0], y) && Global.fuzzy.eq(z, coefficients[1] * x + coefficients[2]);
624
625    case Z_CONSTANT:
626      return Global.fuzzy.eq(coefficients[0], z) && Global.fuzzy.eq(y, coefficients[1] * x + coefficients[2]);
627
628    case NO_CONSTANT_VARIABLES:
629      return Global.fuzzy.eq(y, coefficients[0] * x + coefficients[1]) && Global.fuzzy.eq(z, coefficients[2] * x + coefficients[3]);
630
631    // error if the code gets here because that means, we are trying to fit the points to a linear
632    // when the points already form a plane
633    default:
634      throw new RuntimeException("try_points_linear called when points already form a plane (line_flag = NO_LINE)");
635    }
636  }
637
638  /**
639   * Returns the two points that have the maximum separation in pa.
640   *
641   * @param pa array of points. Must have at least 2 elements. Can be any length and can contain
642   *      nulls (which will be ignored).
643   * @return a 2-element array containing the most-separated elements. Returns null if no pair of
644   *      elements is NaN-free.
645   */
646  private Point @Nullable [] maxsep_doubles(@Nullable Point[] pa) {
647
648    Point p1 = null;
649    Point p2 = null;
650
651    // max_separation is the separation metric for the most separated
652    // pair of points.  We use the sum of the separations as the
653    // metric.
654    double max_separation = Double.MIN_VALUE;
655
656    for (int i = 0; i < pa.length - 1; i++) {
657      for (int j = i + 1; j < pa.length; j++) {
658        double separation = separation(pa[i], pa[j]);
659        if (separation > max_separation) {
660          max_separation = separation;
661          p1 = pa[i];
662          p2 = pa[j];
663        }
664      }
665    }
666
667    // This can happen if enough values are NaN.
668    // if (p1 == null) {
669    //   throw new IllegalArgumentException(Arrays.toString(pa));
670    // }
671
672    if (Debug.logDetail()) {
673      wrapper.log("maxsep_doubles = %s %s", p1, p2);
674    }
675    if (p1 == null) {
676      return null;
677    }
678    return new Point[] { p1, p2 };
679  }
680
681  /**
682   * Calculates the coefficients for the line in 3D
683   * and sets the coefficients[] by side effect.
684   *
685   * @param p0 first point on the line
686   * @param p1 second point on the line
687   */
688  //
689  //  the line can fall under 3 categories
690  //
691  //  1. parallel to an axis which means that 2 of the following is true
692  //  x = A
693  //  y = B
694  //  z = C
695  //
696  //  and the two that are true form a line in the 3rd dimension (the one that's
697  //  not constant)
698  //
699  //  2. parallel to the plane of an axis, which means that one of the following
700  //  is true
701  //  x = A
702  //  y = B
703  //  z = C
704  //
705  //  and the line is just a line in the non-constacoefficientble plane projected
706  //  onto the third dimension
707  //  and the equation for that line is:
708  //  x = A [suppose that x is the constant one]
709  //  y = Bz + C
710  //
711  //
712  //  3. none of the above and the equation of the line is just:
713  //  y = Ax + B
714  //  z = Cx + D
715  //
716  //  Note that A, B, C, and D are coefficients that are stored in coefficients[]
717
718  private void calc_bi_linear(Point p0, Point p1) {
719
720    if ((p0.x == p1.x) && (p0.y == p1.y) && (p0.z != p1.z)) {
721      // x and y have constant values
722      // line is defined by x = A, y = B
723      line_flag = Flag.XY_CONSTANT;
724      coefficients[0] = p0.x;
725      coefficients[1] = p0.y;
726      coefficients[2] = 0;
727      coefficients[3] = 0;
728      return;
729    }
730    if ((p0.x == p1.x) && (p0.z == p1.z) && (p0.y != p1.y)) {
731      // x and z have constant values
732      // line is defined by x = A, z = C
733      line_flag = Flag.XZ_CONSTANT;
734      coefficients[0] = p0.x;
735      coefficients[1] = p0.z;
736      coefficients[2] = 0;
737      coefficients[3] = 0;
738      return;
739    }
740    if ((p0.y == p1.y) && (p0.z == p1.z) && (p0.x != p1.x)) {
741      // y and z have constant values
742      // line is defined by y = B, z = C
743      line_flag = Flag.YZ_CONSTANT;
744      coefficients[0] = p0.y;
745      coefficients[1] = p0.z;
746      coefficients[2] = 0;
747      coefficients[3] = 0;
748      return;
749    }
750    if ((p0.x == p1.x) && (p0.y != p1.y) && (p0.z != p1.z)) {
751      // x has a constant value
752      // plane (degenerate) is defined by x = A, z = By + C
753      line_flag = Flag.X_CONSTANT;
754      coefficients[0] = p0.x;
755      coefficients[1] = (p1.z - p0.z)/(double)(p1.y - p0.y);
756      coefficients[2] = (p0.z * p1.y - p0.y * p1.z)/(double)(p1.y - p0.y);
757      coefficients[3] = 0;
758      return;
759    }
760    if ((p0.y == p1.y) && (p0.x != p1.x) && (p0.z != p1.z)) {
761      // y has a constant value
762      // plane (degenerate) is defined by y = B, z = Ax + C
763      line_flag = Flag.Y_CONSTANT;
764      coefficients[0] = p0.y;
765      coefficients[1] = (p1.z - p0.z)/(double)(p1.x - p0.x);
766      coefficients[2] = (p0.z * p1.x - p0.x * p1.z)/(double)(p1.x - p0.x);
767      coefficients[3] = 0;
768      return;
769    }
770    if ((p0.z == p1.z) && (p0.x != p1.x) && (p0.y != p1.y)) {
771      // z has a constant value
772      // plane (degenerate) is defined by z = C, y = Ax + B
773      line_flag = Flag.Z_CONSTANT;
774      coefficients[0] = p0.z;
775      coefficients[1] = (p1.y - p0.y)/(double)(p1.x - p0.x);
776      coefficients[2] = (p0.y * p1.x - p0.x * p1.y)/(double)(p1.x - p0.x);
777      coefficients[3] = 0;
778      return;
779    }
780    if ((p0.x != p1.x) && (p0.y != p1.y) && (p0.z != p1.z)) {
781      // no variables have a constant value
782      // plane (degenerate) is defined by y = Ax + B, z = Cx + D
783      line_flag = Flag.NO_CONSTANT_VARIABLES;
784      coefficients[0] = (p1.y - p0.y)/(double)(p1.x - p0.x);
785      coefficients[1] = (p0.y * p1.x - p0.x * p1.y)/(double)(p1.x - p0.x);
786      coefficients[2] = (p1.z - p0.z)/(double)(p1.x - p0.x);
787      coefficients[3] = (p0.z * p1.x - p0.x * p1.z)/(double)(p1.x - p0.x);
788    }
789  }
790
791  /**
792   *  Calculates new coefficients that for the new point. Uses the
793   *  new coefficients if they are relatively close to to the previous
794   *  ones. Kills off the invariant if they are not.
795   *
796   *  @return true if the new equation worked, false otherwise
797   */
798  public boolean try_new_equation(long x, long y, long z) {
799
800    // Calculate max separation using this point and the existing 3 points.
801    def_points[3] = new Point(x, y, z);
802    double sep = maxsep_triples(def_points);
803
804    // If this point increased the separation, recalculate a, b, and c.
805    if (sep > separation) {
806      separation = sep;
807      double[] coef;
808      try {
809        // Requires that the first 3 elements of def_points must be non-null.
810        coef = calc_tri_linear(def_points);
811        if (Debug.logDetail()) {
812          wrapper.log("Calc new plane with points %s %s %s %s",
813                       def_points[0], def_points[1], def_points[2], def_points[3]);
814        }
815      } catch (Exception e) {
816        return false;
817      }
818
819      // if the a, b, or c is a new min/max remember it.
820      if (coef[0] < min_a) min_a = coef[0];
821      if (coef[0] > max_a) max_a = coef[0];
822      if (coef[1] < min_b) min_b = coef[1];
823      if (coef[1] > max_b) max_b = coef[1];
824      if (coef[2] < min_c) min_c = coef[2];
825      if (coef[2] > max_c) max_c = coef[2];
826      if (coef[3] < min_d) min_d = coef[3];
827      if (coef[3] > max_d) {
828        max_d = coef[3];
829      }
830
831      // Pick a new a, b, and c as the average of their endpoints
832      a = (min_a + max_a) / 2;
833      b = (min_b + max_b) / 2;
834      c = (min_c + max_c) / 2;
835      d = (min_d + max_d) / 2;
836      if (Invariant.logOn() || debug.isLoggable(Level.FINE)) {
837        wrapper.log(debug, wrapper.ppt.name() + ": Trying new a (" + a
838                     + "), b (" + b + "), c (" + c + "), and d (" + d + ")");
839      }
840      wrapper.log("min max: Trying new a (%s), b (%s), c (%s), and d (%s)",
841                   a, b, c, d);
842      // if the new coefficients are 'equal' to their min and max and
843      // this point fits, then this new equation is good enough both
844      // for existing points and the new point.
845      if (Global.fuzzy.eq(a, min_a) && Global.fuzzy.eq(a, max_a)
846          && Global.fuzzy.eq(b, min_b) && Global.fuzzy.eq(b, max_b)
847          && Global.fuzzy.eq(c, min_c) && Global.fuzzy.eq(c, max_c)
848          && Global.fuzzy.eq(d, min_d) && Global.fuzzy.eq(d, max_d)
849          && Global.fuzzy.eq(-d, a * x + b * y + c * z)) {
850        if (debug.isLoggable(Level.FINE)) {
851          debug.fine(wrapper.ppt.name() + ": New a (" + a + ") and b ("
852                      + b + ") and c (" + c + ")");
853        }
854        return true;
855      } else {
856        return false;
857      }
858    } else { // this point doesn't increase the separation
859
860      return false;
861    }
862  }
863
864  /**
865   * Calculates the separation between p1 and p2.
866   *
867   * @param p1 first point
868   * @param p2 second point
869   *
870   * @return the distance between p1 and p2, or 0 if either point is null, or NaN if either point
871   * contains a NaN
872   */
873  double separation(Point p1, Point p2) {
874
875      // These variable Types are double so the values won't wrap around.
876
877    // make sure both points are specified
878    if ((p1 == null) || (p2 == null)) {
879      return 0;
880    }
881
882    double xsep = (p1.x - p2.x);
883    double ysep = (p1.y - p2.y);
884    double zsep = (p1.z - p2.z);
885    return Math.sqrt(xsep * xsep + ysep * ysep + zsep * zsep);
886  }
887
888  /**
889   * Calculates the three points that have the maximum separation in pa and places them as the first
890   * three elements of pa.
891   *
892   * @param pa array of points. Must have at least 3 elements. Can be any length and can contain
893   *     nulls (which will be ignored). Is side-effected so that the first three elements contain
894   *     the points with the maximum total separation; this may introduce duplicates into the array.
895   *     The first three elements are null if no triple of elements is NaN-free.
896   * @return the maximum separation found
897   */
898  double maxsep_triples(@Nullable Point[] pa) {
899
900    // cache values for the (square of the) distance between each pair of
901    // points, to avoid duplicating work
902    double[][] separations = new double[pa.length][pa.length];
903    for (int i = 0; i < pa.length - 1; i++) {
904      for (int j = i + 1; j < pa.length; j++) {
905        separations[i][j] = separation(pa[i], pa[j]);
906      }
907    }
908
909    Point p1 = null;
910    Point p2 = null;
911    Point p3 = null;
912
913    // max_separation is the separation metric for the most separated
914    // triple of points.  We use the sum of the separations as the
915    // metric.  (The metric "min of the three separations" does not work
916    // because it doesn't choose a unique set of points, making the
917    // result dependent on the order in which the points are seen; more
918    // seriously, it may choose a collinear set of points even when a
919    // non-collinear set exists.)
920    double max_separation = Double.MIN_VALUE;
921
922    for (int i = 0; i < pa.length - 2; i++) {
923      for (int j = i + 1; j < pa.length - 1; j++) {
924        double sep_i_j = separations[i][j];
925        for (int k = j + 1; k < pa.length; k++) {
926
927          // using the sum of separations is valid, as long as
928          // separation is the actual distance between points and not
929          // the square of the distance
930          double separation = sep_i_j + separations[i][k] + separations[j][k];
931          if (separation > max_separation) {
932            max_separation = separation;
933            p1 = pa[i];
934            p2 = pa[j];
935            p3 = pa[k];
936          }
937        }
938      }
939    }
940
941    pa[0] = p1;
942    pa[1] = p2;
943    pa[2] = p3;
944
945    if (Debug.logDetail()) {
946      wrapper.log("maxsep_triples = %s %s %s", pa[0], pa[1], pa[2]);
947    }
948    return max_separation;
949  }
950
951  // Given ((x0,y0,z0),(x1,y1,z1), (x2,y2,z2), calculate a, b, c and d
952  // such that ax + by + cz + d = 0.
953  // a, b, c and d are mutually prime, integers and a is positive.
954  //
955  // A visual explanation of the math can be found in:
956  // http://www.efm.leeds.ac.uk/CIVE/CIVE2599/cive2599-summary-overheads-full.pdf
957  //
958  // The standard form of a plane can be calculated using a vector
959  // normal "n" <n1, n2, n3> to the plane and a vector from the origin
960  // to the plane "r" e.g. <x0, y0, z0>.
961  //
962  // The vector normal to the plane can be calculated by doing the cross product
963  // of two line segments on the plane (<point2 - point0> and <point1 - point0>)
964  //
965  // Given the normal vector and a point on the plane, the standard form
966  // of the plane is:
967  // n1*x + n2*y + n3*z = p where p is the dot product of r and n
968  //
969  //
970
971  /**
972   * Calculates the coefficients (a, b, c) and constant (d) for the equation ax + by + cz + d = 0
973   * using the first three points in points.
974   *
975   * @param points array of points to use to calculate the coefficents. Only the first three points
976   *     (all of which must be non-null) are used.
977   *
978   * @return a four element array where a is the first element, b the second, c the third, and d is
979   * the fourth. All elements are mutually prime, integers and a is positive.
980   */
981  // TODO: should just pass in the first three elements rather than passing
982  // in the point array.
983  public double[] calc_tri_linear(@Nullable Point[] points) { // TODO: remove annotation after refactoring method
984
985    Point p0 = points[0];
986    Point p1 = points[1];
987    Point p2 = points[2];
988
989    // the following code is taken from the source page of
990    // http://jblanco_60.tripod.com/plane.html (the math behind the
991    // implementation is detailed above)
992
993    long px = p0.x;
994    long py = p0.y;
995    long pz = p0.z;
996    long qx = p1.x;
997    long qy = p1.y;
998    long qz = p1.z;
999    long rx = p2.x;
1000    long ry = p2.y;
1001    long rz = p2.z;
1002
1003    // define the i, j, k components of the two line segments on the plane
1004    long a1 = (qx - px);
1005    long a2 = (qy - py);
1006    long a3 = (qz - pz);
1007
1008    long b1 = (rx - px);
1009    long b2 = (ry - py);
1010    long b3 = (rz - pz);
1011
1012    // calculate the i, j, k components of the cross product
1013    long i_var = (a2 * b3) - (b2 * a3);
1014    long j_var = -((a1 * b3) - (b1 * a3));
1015    long k_var = (a1 * b2) - (b1 * a2);
1016
1017    // Calculate the value of the constant.  Note that given the format
1018    // of the point-normal form, we multiply by the negative of P(x,y,z)
1019
1020    long Stand_const = ((i_var*(-px)) + (j_var*(-py)) + (k_var*(-pz)));
1021
1022    // Put the coefficients in lowest terms by dividing by the gcd.
1023    // If the gcd is 1, no harm done -- divide by 1.
1024    // Type is double to ensure floating-point division.
1025    double myGCD = MathPlume.gcd(MathPlume.gcd(i_var, j_var) , MathPlume.gcd(k_var, Stand_const));
1026    double standard_i_var = i_var / myGCD;
1027    double standard_j_var = j_var / myGCD;
1028    double standard_k_var = k_var / myGCD;
1029    double standard_Stand_const = Stand_const / myGCD;
1030
1031    if (wrapper != null) {
1032      wrapper.log("Preliminary: %s %s GCD: %s", i_var, j_var, MathPlume.gcd(i_var, j_var));
1033      wrapper.log("Preliminary: %s %s GCD: %s", k_var, Stand_const, MathPlume.gcd(k_var, Stand_const));
1034      wrapper.log("GCD: %s", myGCD);
1035    }
1036
1037    // Ensure postive x coefficient by negating if "x" term is negative.
1038    if (i_var < 0) {
1039      standard_i_var = - standard_i_var;
1040      standard_j_var = - standard_j_var;
1041      standard_k_var = - standard_k_var;
1042      standard_Stand_const = - standard_Stand_const;
1043    }
1044
1045    double[] coef = new double[] {standard_i_var,
1046                                  standard_j_var,
1047                                  standard_k_var, standard_Stand_const };
1048
1049    return coef;
1050  }
1051
1052  public boolean enoughSamples(@GuardSatisfied LinearTernaryCore this) {
1053    return isActive();
1054  }
1055
1056  // If the values form a line,
1057  // we have no confidence that it's a plane.
1058  // The line is less interesting because the invariant
1059  // is already defined by a linear binary invariant.
1060  public double computeConfidence() {
1061    if (isActive()) {
1062      return Invariant.conf_is_ge(values_seen, MINTRIPLES);
1063    } else {
1064      return 0;
1065    }
1066
1067  }
1068
1069  public String repr(@GuardSatisfied LinearTernaryCore this) {
1070    return "LinearTernaryCore" + wrapper.varNames() + ": a=" + a
1071      + ",b=" + b
1072      + ",c=" + c
1073      + ",d=" + d
1074      + ",values_seen=" + values_seen;
1075  }
1076
1077  public String point_repr(Point p) {
1078    if (p == null) {
1079      return "null";
1080    } else {
1081      return "<" + p.x + "," + p.y + "," + p.z + ">";
1082    }
1083  }
1084
1085  public String cache_repr() {
1086    StringBuilder result = new StringBuilder();
1087    for (int i = 0; i < MINTRIPLES; i++) {
1088      if (i != 0) result.append("; ");
1089      result.append(point_repr(def_points[i]));
1090    }
1091    return result.toString();
1092  }
1093
1094  // In this class for convenience (avoid prefixing "LinearBinaryCore").
1095  static String formatTerm(double coeff, @Nullable String varname, boolean first) {
1096    return LinearBinaryCore.formatTerm(coeff, varname, first);
1097  }
1098
1099  @SideEffectFree
1100  public String format_using(OutputFormat format,
1101                             String vix, String viy, String viz,
1102                             double a, double b, double c, double d) {
1103
1104    if ((a == 0) && (b == 0) && (c == 0) && (d == 0)) {
1105      return wrapper.format_too_few_samples(format, null);
1106    }
1107
1108    if (format == OutputFormat.SIMPLIFY) {
1109      return format_simplify(vix, viy, viz, a, b, c, d);
1110    }
1111
1112    if ((format.isJavaFamily())) {
1113
1114        return formatTerm(a, vix, true)
1115          + formatTerm(b, viy, (a == 0))
1116          + formatTerm(c, viz, (a == 0) && (b == 0))
1117          + formatTerm(d, null, (a == 0) && (b == 0) && (c == 0))
1118          + " == 0";
1119
1120      }
1121
1122    if ((format == OutputFormat.DAIKON)
1123        || (format == OutputFormat.ESCJAVA)
1124        || (format == OutputFormat.CSHARPCONTRACT)) {
1125        String eq = " == ";
1126        return formatTerm(a, vix, true)
1127          + formatTerm(b, viy, (a == 0))
1128          + formatTerm(c, viz, (a == 0) && (b == 0))
1129          + formatTerm(d, null, (a == 0) && (b == 0) && (c == 0))
1130          + eq + "0";
1131      }
1132
1133    throw new Error("unrecognized output format " + format);
1134    // return null;
1135  }
1136
1137  public static String format_simplify(String vix, String viy, String viz,
1138                                       double da, double db, double dc, double dd) {
1139    int ia = (int) da;
1140    int ib = (int) db;
1141    int ic = (int) dc;
1142    int id = (int) dd;
1143
1144    String str_a, str_b, str_c, str_d;
1145    if (ia != da || ib != db || ic != dc || id != dd) {
1146      // floating point
1147
1148      // Disabled for the moment, since mixing integers and floating
1149      // literals seems to give Simplify indigestion. For instance:
1150      // (BG_PUSH (<= w 3))
1151      // (BG_PUSH (EQ 0 (* w 2.0d0)))
1152      // (<= w 1)
1153
1154      // str_a = Invariant.simplify_format_double(da);
1155      // str_b = Invariant.simplify_format_double(db);
1156      // str_c = Invariant.simplify_format_double(dc);
1157      // str_d = Invariant.simplify_format_double(dd);
1158      return "(AND)"; // really s/b format_unimplemented
1159    } else {
1160      // integer
1161      str_a = Invariant.simplify_format_long(ia);
1162      str_b = Invariant.simplify_format_long(ib);
1163      str_c = Invariant.simplify_format_long(ic);
1164      str_d = Invariant.simplify_format_long(id);
1165    }
1166
1167    String str_z = viz;
1168    String str_x = vix;
1169    String str_y = viy;
1170    String str_ax = (da == 1.0) ? str_x : "(* " + str_a + " " + str_x + ")";
1171    String str_by = (db == 1.0) ? str_y : "(* " + str_b + " " + str_y + ")";
1172    String str_cz = (dc == 1.0) ? str_z : "(* " + str_c + " " + str_z + ")";
1173    String str_axPbyPcz = "(+ " + str_ax + " " + str_by + " " + str_cz + ")";
1174    String str_axPbyPczPd = (dd == 0.0) ? str_axPbyPcz :
1175      "(+ " + str_axPbyPcz + " " + str_d + ")";
1176    return "(EQ 0 " + str_axPbyPczPd + ")";
1177  }
1178
1179  @SideEffectFree
1180  public String format_using(OutputFormat format,
1181                             String vix, String viy, String viz
1182                             ) {
1183    String result = format_using(format, vix, viy, viz, a, b, c, d);
1184    if (result != null) {
1185      return result;
1186    }
1187
1188    return wrapper.format_unimplemented(format);
1189  }
1190
1191  // // Format as "x = cy+d" instead of as "y = ax+b".
1192  // public String format_reversed(String x, String y) {
1193  //   assert a == 1 || a == -1;
1194  //   return format(y, x, a, -b/a);
1195  // }
1196
1197  @Pure
1198  public boolean isSameFormula(LinearTernaryCore other) {
1199
1200    // If both have yet to see enough values
1201    if (!isActive() && !other.isActive()) {
1202
1203      // Same formula if all of the points match
1204      if (values_seen != other.values_seen) {
1205        return false;
1206      }
1207      // This and elsewhere seem to assume that values_seen < MINTRIPLES;
1208      // is that necessarily true?
1209      for (int ii = 0; ii < values_seen; ii++) {
1210        // used to use !=, but that test seems wrong
1211        if (!def_points[ii].equals(other.def_points[ii])) {
1212          return false;
1213        }
1214      }
1215      return true;
1216    } else {
1217      return ((values_seen >= MINTRIPLES)
1218              && (other.values_seen >= MINTRIPLES)
1219              && (a == other.a)
1220              && (b == other.b)
1221              && (c == other.c)
1222              && (d == other.d));
1223    }
1224  }
1225
1226  @Pure
1227  public boolean isExclusiveFormula(LinearTernaryCore other) {
1228    if (!isActive()
1229        || !other.isActive()) {
1230      return false;
1231    }
1232
1233    return ((a == other.a)
1234            && (b != other.b)
1235            && (c != other.c)
1236            && (d != other.d));
1237  }
1238
1239  /**
1240   * In general, we can't merge formulas, but we can merge invariants with too few samples to have
1241   * formed a plane with invariants with enough samples. And those will appear to have different
1242   * formulas.
1243   */
1244  public boolean mergeFormulasOk() {
1245    return true;
1246  }
1247
1248  /**
1249   * Merge the linear ternary cores in cores to form a new core. Any core in the list that has seen
1250   * enough points to define a plane, must define the same plane. Any cores that have not yet seen
1251   * enough points, will have each of their points applied to that invariant. The merged core is
1252   * returned. Null is returned if the cores don't describe the same plane
1253   *
1254   * @param cores list of LinearTernary cores to merge. They should all be permuted to match the
1255   *     variable order in ppt.
1256   */
1257  public @Nullable LinearTernaryCore merge(List<LinearTernaryCore> cores, Invariant wrapper) {
1258
1259    // Look for any active planes.  All must define the same plane
1260    LinearTernaryCore first = null;
1261    for (LinearTernaryCore c : cores) {
1262      if (!c.isActive()) {
1263        continue;
1264      }
1265      if (first == null) {
1266        first = c.clone();
1267      } else {
1268        if (!Global.fuzzy.eq(first.a, c.a)
1269            || !Global.fuzzy.eq(first.b, c.b)
1270            || !Global.fuzzy.eq(first.c, c.c)
1271            || !Global.fuzzy.eq(first.d, c.d))
1272          return null;
1273      }
1274    }
1275
1276    // If no active planes were found, created an empty core
1277    if (first == null) {
1278      first = new LinearTernaryCore(wrapper);
1279    } else {
1280      first.wrapper = wrapper;
1281    }
1282
1283    // Merge in any points from non-active cores
1284    // Note that while null points can't show up in a core that is not active
1285    // because it hasn't seen enough values, they can show up if the core
1286    // is not active because it is a line.
1287    for (LinearTernaryCore c : cores) {
1288      if (c.isActive()) {
1289        continue;
1290      }
1291      for (int j = 0; j < c.values_seen; j++) {
1292        Point cp = c.def_points[j];
1293        if (cp == null) {
1294          continue;
1295        }
1296        if (Debug.logDetail()) {
1297          wrapper.log("Adding point %s from %s", cp, c.wrapper.ppt);
1298        }
1299        InvariantStatus stat = first.add_modified(cp.x, cp.y, cp.z, 1);
1300        if (stat == InvariantStatus.FALSIFIED) {
1301          return null;
1302        }
1303        // if (wrapper.is_false())
1304        //  return null;
1305      }
1306    }
1307
1308    return first;
1309  }
1310}