Adaptive Orientation

Orientation.java

/*
 * Licensed to Elasticsearch under one or more contributor
 * license agreements. See the NOTICE file distributed with
 * this work for additional information regarding copyright
 * ownership. Elasticsearch licenses this file to you under
 * the Apache License, Version 2.0 (the "License"); you may
 * not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *    http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing,
 * software distributed under the License is distributed on an
 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
 * KIND, either express or implied.  See the License for the
 * specific language governing permissions and limitations
 * under the License.
 */
package org.apache.lucene.geo;

public class Orientation {

  private static double check;
  private static double lastcheck;
  private static boolean everyOther;
  private static double epsilon;
  private static double splitter;

  private static final double half = 0.5d;
  private static final double resulterrbound;
  private static final double ccwerrboundA;
  private static final double ccwerrboundB;
  private static final double ccwerrboundC;

  static {
    everyOther = true;
    epsilon = 1.0;
    splitter = 1.0;
    check = 1.0;
    /* Repeatedly divide `epsilon' by two until it is too small to add to    */
    /*   one without causing roundoff.  (Also check if the sum is equal to   */
    /*   the previous sum, for machines that round up instead of using exact */
    /*   rounding.  Not that this library will work on such machines anyway. */
    do {
      lastcheck = check;
      epsilon *= half;
      if (everyOther == true) {
        splitter *= 2.0;
      }
      everyOther = !everyOther;
      check = 1.0 + epsilon;
    } while ((check != 1.0) && (check != lastcheck));
    splitter += 1.0;

    /* Error bounds for orientation and incircle tests. */
    resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
    ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
    ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
    ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
  }

  public static double orient2d(double ax, double ay, double bx, double by, double cx, double cy) {
    double detleft, detright, det;
    double detsum, errbound;

    detleft = (ax - cx) * (by - cy);
    detright = (ay - cy) * (bx - cx);
    det = detleft - detright;

    if (detleft > 0.0) {
      if (detright <= 0.0) {
        return det;
      } else {
        detsum = detleft + detright;
      }
    } else if (detleft < 0.0) {
      if (detright >= 0.0) {
        return det;
      } else {
        detsum = -detleft - detright;
      }
    } else {
      return det;
    }

    errbound = ccwerrboundA * detsum;
    if ((det >= errbound) || (-det >= errbound)) {
      return det;
    }

    return orient2dadapt(ax, ay, bx, by, cx, cy, detsum);
  }

  private static double orient2dadapt(double ax, double ay, double bx, double by, double cx, double cy, double detsum)
  {
    double acx = ax - cx;
    double acy = ay - cy;
    double bcx = bx - cx;
    double bcy = by - cy;

    double acxtail, acytail, bcxtail, bcytail;
    double detleft, detright;
    double detlefttail, detrighttail;
    double det, errbound;
    double[] B;
    double[] C1 = new double[8];
    double[] C2 = new double[12];
    double[] D = new double[16];
    int C1length, C2length, Dlength;
    double[] u;

    detleft = acx * bcy;
    detlefttail = twoProductTail(acx, bcy, detleft);
    detright = acy * bcx;
    detrighttail = twoProductTail(acy, bcx, detright);

    B = twoTwoDiff(detleft, detlefttail, detright, detrighttail);

    det = estimate(B);
    errbound = ccwerrboundB * detsum;
    if ((det >= errbound) || (-det >= errbound)) {
      return det;
    }

    acxtail = twoDiffTail(ax, cx, acx);
    bcxtail = twoDiffTail(bx, cx, bcx);
    acytail = twoDiffTail(ay, cy, acy);
    bcytail = twoDiffTail(by, cy, bcy);

    if ((acxtail == 0.0) && (acytail == 0.0)
        && (bcxtail == 0.0) && (bcytail == 0.0)) {
      return det;
    }

    errbound = ccwerrboundC * detsum + resulterrbound * Math.abs(det);
    det += (acx * bcytail + bcy * acxtail)
        - (acy * bcxtail + bcx * acytail);
    if ((det >= errbound) || (-det >= errbound)) {
      return det;
    }

    double[] s = twoProduct(acxtail, bcy);
    double[] t = twoProduct(acytail, bcx);
    u = twoTwoDiff(s[0], s[1], t[0], t[1]);
    C1length = fastExpansionSumZeroelim(4, B, 4, u, C1);

    s = twoProduct(acx, bcytail);
    t = twoProduct(acy, bcxtail);
    u = twoTwoDiff(s[0], s[1], t[0], t[1]);
    C2length = fastExpansionSumZeroelim(C1length, C1, 4, u, C2);

    s = twoProduct(acxtail, bcytail);
    t = twoProduct(acytail, bcxtail);
    u = twoTwoDiff(s[0], s[1], t[0], t[1]);
    Dlength = fastExpansionSumZeroelim(C2length, C2, 4, u, D);

    return(D[Dlength - 1]);
  }

  private static int fastExpansionSumZeroelim(int elen, double[] e, int flen, double[] f, double[] h)  /* h cannot be e or f. */
  {
    double Q;
    double Qnew;
    double hh;
    int eindex, findex, hindex;
    double enow, fnow;

    enow = e[0];
    fnow = f[0];
    eindex = findex = 0;
    if ((fnow > enow) == (fnow > -enow)) {
      Q = enow;
      enow = e[++eindex];
    } else {
      Q = fnow;
      fnow = f[++findex];
    }
    hindex = 0;
    if ((eindex < elen) && (findex < flen)) {
      if ((fnow > enow) == (fnow > -enow)) {
        double[] temp = fastTwoSum(enow, Q);
        Qnew = temp[0];
        hh = temp[1];
        enow = e[++eindex];
      } else {
        double[] temp = fastTwoSum(fnow, Q);
        Qnew = temp[0];
        hh = temp[1];
        fnow = f[++findex];
      }
      Q = Qnew;
      if (hh != 0.0) {
        h[hindex++] = hh;
      }
      while ((eindex < elen) && (findex < flen)) {
        if ((fnow > enow) == (fnow > -enow)) {
          double[] temp;
          temp = twoSum(Q, enow);
          Qnew = temp[0];
          hh = temp[1];
          enow = e[++eindex];
        } else {
          double[] temp;
          temp = twoSum(Q, fnow);
          Qnew = temp[0];
          hh = temp[1];
          fnow = f[++findex];
        }
        Q = Qnew;
        if (hh != 0.0) {
          h[hindex++] = hh;
        }
      }
    }
    while (eindex < elen) {
      double[] temp;
      temp = twoSum(Q, enow);
      Qnew = temp[0];
      hh = temp[1];
      enow = e[++eindex];
      Q = Qnew;
      if (hh != 0.0) {
        h[hindex++] = hh;
      }
    }
    while (findex < flen) {
      double[] temp = twoSum(Q, fnow);
      Qnew = temp[0];
      hh = temp[1];
      fnow = f[++findex];
      Q = Qnew;
      if (hh != 0.0) {
        h[hindex++] = hh;
      }
    }
    if ((Q != 0.0) || (hindex == 0)) {
      h[hindex++] = Q;
    }
    return hindex;
  }


  private static double[] twoProduct(double a, double b) {
    double[] result = new double[2];
    result[0] = a * b;
    result[1] = twoProductTail(a, b, result[0]);
    return result;
  }

  private static double estimate(double[] e) {
    double q = e[0];
    for (int eindex = 1; eindex < e.length; ++eindex) {
      q += e[eindex];
    }
    return q;
  }

  private static double twoProductTail(double a, double b, double x) {
    double[] aSplit = split(a);
    double[] bSplit = split(b);
    double err1 = x - (aSplit[0] * bSplit[0]);
    double err2 = err1 - (aSplit[1] * bSplit[0]);
    double err3 = err2 - (aSplit[0] * bSplit[1]);
    return (aSplit[1] * bSplit[1]) - err3;
  }

  /** splits into double[0] = aHi, double[1] = aLo */
  private static double[] split(double a) {
    double c = splitter * a;
    double abig = (c - a);
    double[] result = new double[2];
    result[0] = c - abig;
    result[1] = a - result[0];
    return result;
  }

  private static double twoDiffTail(double a, double b, double x) {
    double bvirt = a - x;
    double avirt = x + bvirt;
    double bround = bvirt - b;
    double around = a - avirt;
    return around + bround;
  }

  private static double[] twoDiff(double a, double b) {
    double[] result = new double[2];
    result[0] = a - b;
    result[1] = twoDiffTail(a, b, result[0]);
    return result;
  }

  private static double twoSumTail(double a, double b, double x) {
    double bvirt = x - a;
    double avirt = x - bvirt;
    double bround = b - bvirt;
    double around = a - avirt;
    return around + bround;
  }

  private static double[] twoSum(double a, double b) {
    double[] result = new double[2];
    result[0] = a + b;
    result[1] = twoSumTail(a, b, result[0]);
    return result;
  }

  private static double fastTwoSumTail(double a, double b, double x) {
    double bvirt = x - a;
    return b - bvirt;
  }

  private static double[] fastTwoSum(double a, double b) {
    double[] result = new double[2];
    result[0] = a + b;
    result[1] = fastTwoSumTail(a, b, result[0]);
    return result;
  }

  private static double[] twoOneDiff(double a1, double a0, double b) {
    double[] twoDiff = twoDiff(a0, b);
    double[] twoSum = twoDiff(a1, twoDiff[0]);
    return new double[] {twoDiff[1], twoSum[0], twoSum[1]};
  }

  private static double[] twoTwoDiff(double a1, double a0, double b1, double b0) {
    double[] twoOneDiffA = twoOneDiff(a1, a0, b0);
    double[] twoOneDiffB = twoOneDiff(twoOneDiffA[0], twoOneDiffA[1], b1);
    return new double[] {twoOneDiffB[0], twoOneDiffB[1], twoOneDiffB[2], twoOneDiffA[2]};
  }
}

It seems there is a small bug, I get in some situation an AAOB exception. I think in method fastExpansionSumZeroelim the conditions should be like:

while ((eindex < elen - 1) && (findex < flen - 1)) {
....
}

I have run some performance test for points:

||Approach||Index time (sec)||Force merge time (sec)||Index size (GB)||Reader heap (MB)||
          ||Dev||Base||Diff ||Dev  ||Base  ||diff   ||Dev||Base||Diff||Dev||Base||Diff ||
|points|231.3s|224.7s| 3%|0.0s|0.0s| 0%|0.51|0.51| 0%|1.72|1.72| 0%|
()
||Approach||Shape||M hits/sec      ||QPS            ||Hit count      ||
                 ||Dev||Base ||Diff||Dev||Base||Diff||Dev||Base||Diff||
|points|poly 10|84.62|87.94|-4%|53.51|55.61|-4%|355809475|355809475| 0%|
|points|distance|79.76|78.90| 1%|46.86|46.36| 1%|382961957|382961957| 0%|
|points|polyRussia|17.93|19.14|-6%|5.11|5.45|-6%|3508846|3508846| 0%|
|points|polyMedium|9.68|9.64| 0%|118.65|118.10| 0%|2693559|2693559| 0%|
|points|nearest 10|0.02|0.02|-1%|1555.49|1564.59|-1%|60844404|60844404| 0%|
|points|box|79.82|79.79| 0%|81.22|81.19| 0%|221118844|221118844| 0%|
|points|sort|42.58|41.64| 2%|43.33|42.37| 2%|221118844|221118844| 0%|

and for shapes:

|Approach||Index time (sec)||Force merge time (sec)||Index size (GB)||Reader heap (MB)||
          ||Dev||Base||Diff ||Dev  ||Base  ||diff   ||Dev||Base||Diff||Dev||Base||Diff ||
|shapes|1337.8s|1294.1s| 3%|0.0s|0.0s| 0%|2.15|2.15| 0%|1.82|1.82| 0%|
()
||Approach||Shape||M hits/sec      ||QPS            ||Hit count      ||
                 ||Dev||Base ||Diff||Dev||Base||Diff||Dev||Base||Diff||
|shapes|poly 10|22.12|23.03|-4%|13.99|14.56|-4%|355809475|355809475| 0%|
|shapes|polyRussia|2.99|2.84| 5%|0.85|0.81| 5%|3508846|3508846| 0%|
|shapes|polyMedium|1.42|1.37| 4%|17.43|16.81| 4%|2693559|2693559| 0%|
|shapes|box|37.23|36.06| 3%|37.88|36.69| 3%|221118844|221118844| 0%|

Maybe we can assume that the performance hit is around ~5%, not too bad.