src/main/java/org/apache/datasketches/hll/CubicInterpolation.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one
 * or more contributor license agreements.  See the NOTICE file
 * distributed with this work for additional information
 * regarding copyright ownership.  The ASF 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.datasketches.hll;

import org.apache.datasketches.SketchesArgumentException;

/**
 * @author Lee Rhodes
 * @author Kevin Lang
 */
final class CubicInterpolation {

  /**
   * Cubic interpolation using interpolation X and Y tables.
   *
   * @param xArr xArr
   * @param yArr yArr
   * @param x x
   * @return cubic interpolation
   */
  //Used by AbstractCoupons
  //In C: again-two-registers cubic_interpolate_using_table L1377
  static double usingXAndYTables(final double[] xArr, final double[] yArr,
      final double x) {
    assert (xArr.length >= 4) && (xArr.length == yArr.length);
    if ((x < xArr[0]) || (x > xArr[xArr.length - 1])) {
      throw new SketchesArgumentException("X value out of range: " + x);
    }
    if (x == xArr[xArr.length - 1]) {
      return yArr[yArr.length - 1]; // corner case
    }
    final int offset = findStraddle(xArr, x); //uses recursion
    assert (offset >= 0) && (offset <= (xArr.length - 2));
    if (offset == 0) {
      return interpolateUsingXAndYTables(xArr, yArr, offset, x); // corner case
    }
    if (offset == (xArr.length - 2)) {
      return interpolateUsingXAndYTables(xArr, yArr, offset - 2, x); // corner case
    }
    return interpolateUsingXAndYTables(xArr, yArr, offset - 1, x);
  }

  // In C: again-two-registers cubic_interpolate_aux L1368
  private static double interpolateUsingXAndYTables(final double[] xArr, final double[] yArr,
      final int offset, final double x) {
    return cubicInterpolate(xArr[offset], yArr[offset], xArr[offset + 1], yArr[offset + 1],
        xArr[offset + 2], yArr[offset + 2], xArr[offset + 3], yArr[offset + 3], x);
  }

  //Interpolate using X table and Y stride

  /**
   * Cubic interpolation using interpolation X table and Y stride.
   *
   * @param xArr The x array
   * @param yStride The y stride
   * @param x The value x
   * @return cubic interpolation
   */
  //In C: again-two-registers cubic_interpolate_with_x_arr_and_y_stride L1411
  // Used by HllEstimators
  static double usingXArrAndYStride(
      final double[] xArr, final double yStride, final double x) {
    final int xArrLen = xArr.length;
    final int xArrLenM1 = xArrLen - 1;

    final int offset;
    assert ((xArrLen >= 4) && (x >= xArr[0]) && (x <= xArr[xArrLenM1]));

    if (x ==  xArr[xArrLenM1]) { /* corner case */
      return (yStride * (xArrLenM1));
    }

    offset = findStraddle(xArr, x); //uses recursion
    final int xArrLenM2 = xArrLen - 2;
    assert ((offset >= 0) && (offset <= (xArrLenM2)));

    if (offset == 0) { /* corner case */
      return (interpolateUsingXArrAndYStride(xArr, yStride, (offset - 0), x));
    }
    else if (offset == xArrLenM2) { /* corner case */
      return (interpolateUsingXArrAndYStride(xArr, yStride, (offset - 2), x));
    }
    /* main case */
    return (interpolateUsingXArrAndYStride(xArr, yStride, (offset - 1), x));
  }

  //In C: again-two-registers cubic_interpolate_with_x_arr_and_y_stride_aux L1402
  private static double interpolateUsingXArrAndYStride(final double[] xArr, final double yStride,
      final int offset, final double x) {
    return cubicInterpolate(
        xArr[offset + 0], yStride * (offset + 0),
        xArr[offset + 1], yStride * (offset + 1),
        xArr[offset + 2], yStride * (offset + 2),
        xArr[offset + 3], yStride * (offset + 3),
        x);
  }

  //Cubic Interpolation used by both methods

  // Interpolate using the cubic curve that passes through the four given points, using the
  // Lagrange interpolation formula.
  // In C: again-two-registers cubic_interpolate_aux_aux L1346
  private static double cubicInterpolate(final double x0, final double y0, final double x1,
      final double y1, final double x2, final double y2, final double x3, final double y3,
      final double x) {
    final double l0Numer = (x - x1) * (x - x2) * (x - x3);
    final double l1Numer = (x - x0) * (x - x2) * (x - x3);
    final double l2Numer = (x - x0) * (x - x1) * (x - x3);
    final double l3Numer = (x - x0) * (x - x1) * (x - x2);

    final double l0Denom = (x0 - x1) * (x0 - x2) * (x0 - x3);
    final double l1Denom = (x1 - x0) * (x1 - x2) * (x1 - x3);
    final double l2Denom = (x2 - x0) * (x2 - x1) * (x2 - x3);
    final double l3Denom = (x3 - x0) * (x3 - x1) * (x3 - x2);

    final double term0 = (y0 * l0Numer) / l0Denom;
    final double term1 = (y1 * l1Numer) / l1Denom;
    final double term2 = (y2 * l2Numer) / l2Denom;
    final double term3 = (y3 * l3Numer) / l3Denom;

    return term0 + term1 + term2 + term3;
  }

  //In C: again-two-registers.c find_straddle L1335
  private static int findStraddle(final double[] xArr, final double x) {
    assert ((xArr.length >= 2) && (x >= xArr[0]) && (x <= xArr[xArr.length - 1]));
    return (recursiveFindStraddle(xArr, 0, xArr.length - 1, x));
  }

  //In C: again-two-registers.c find_straddle_aux L1322
  private static int recursiveFindStraddle(final double[] xArr, final int left, final int right,
      final double x) {
    final int middle;
    assert (left < right);
    assert ((xArr[left] <= x) && (x < xArr[right])); /* the invariant */
    if ((left + 1) == right) {
      return (left);
    }
    middle = left + ((right - left) / 2);
    if (xArr[middle] <= x) {
      return (recursiveFindStraddle(xArr, middle, right, x));
    } else {
      return (recursiveFindStraddle(xArr, left, middle, x));
    }
  }

}