content: package com.airbnb.aerosolve.core.function;
  
  import lombok.extern.slf4j.Slf4j;
  import org.apache.commons.lang3.tuple.ImmutablePair;
  import org.apache.commons.lang3.tuple.Pair;
  
  import java.util.Arrays;
  import java.util.List;
  
  @Slf4j
  public class FunctionUtil {
    public static float[] fitPolynomial(float[] data) {
      int numCoeff = 6;
      int iterations = numCoeff * 4;
      float[] initial = new float[numCoeff];
  
      float[] initialStep = new float[numCoeff];
      Arrays.fill(initialStep, 1.0f);
      return optimize(1.0 / 512.0, iterations, initial, initialStep,
          new ImmutablePair<Float, Float>(-10.0f, 10.0f), data);
    }
  
    public static float evaluatePolynomial(float[] coeff, float[] data, boolean overwrite) {
      int len = data.length;
      float err = 0;
      long count = 0;
      for (int i = 0; i < len; i++) {
        float t = (float) i / (len - 1);
        float tinv = 1 - t;
        float diracStart = (i == 0) ? coeff[0] : 0;
        float diracEnd = (i == len - 1) ? coeff[1] : 0;
        double eval = coeff[2] * tinv * tinv * tinv +
            coeff[3] * 3.0 * tinv * tinv * t +
            coeff[4] * 3.0 * tinv * t * t +
            coeff[5] * t * t * t +
            diracStart +
            diracEnd;
        if (data[i] != 0.0) {
          err += Math.abs(eval - data[i]);
          count++;
        }
        if (overwrite) {
          data[i] = (float) eval;
        }
      }
      return err / count;
    }
  
    // CyclicCoordinateDescent
    public static float[] optimize(double tolerance, int iterations,
                                   float[] initial, float[] initialStep,
                                   Pair<Float, Float> bounds, float[] data) {
      float[] best = initial;
      float bestF = evaluatePolynomial(best, data, false);
      int maxDim = initial.length;
      for (int i = 0; i < iterations; ++i) {
        for (int dim = 0; dim < maxDim; ++dim) {
          float step = initialStep[dim];
          while (step > tolerance) {
            float[] left = best.clone();
            left[dim] = Math.max(bounds.getLeft(), best[dim] - step);
            float leftF = evaluatePolynomial(left, data, false);
            float[] right = best.clone();
            right[dim] = Math.min(bounds.getRight(), best[dim] + step);
            float rightF = evaluatePolynomial(right, data, false);
            if (leftF < bestF) {
              best = left;
              bestF = leftF;
            }
            if (rightF < bestF) {
              best = right;
              bestF = rightF;
            }
            step *= 0.5;
          }
        }
      }
      return best;
    }
  
    public static float[] toFloat(List<Double> list) {
      float[] result = new float[list.size()];
      for (int i = 0; i < result.length; i++) {
        result[i] = list.get(i).floatValue();
      }
      return result;
    }
  
    /*
   * @param  tolerance if fitted array's deviation from weights is less than tolerance
   *         use the fitted, otherwise keep original weights.
   * @param  weights the curve you want to smooth
   * @return double errAndCoeff in the weights
     */
    public static double smooth(double tolerance, boolean toleranceIsPercentage, float[] weights) {
      // TODO use apache math's PolynomialCurveFitter
      float[] best = FunctionUtil.fitPolynomial(weights);
      double errAndCoeff = FunctionUtil.evaluatePolynomial(best, weights, false);
      if (toleranceIsPercentage) {
        double absMean = getAbsMean(weights);
        return smoothInternal(errAndCoeff, tolerance * absMean, best, weights) / absMean;
      } else {
        return smoothInternal(errAndCoeff, tolerance, best, weights);
      }
    }
  
    private static double smoothInternal(
        double errAndCoeff, double tolerance, float[] best, float[] weights) {
      if (errAndCoeff < tolerance) {
        FunctionUtil.evaluatePolynomial(best, weights, true);
      }
      return errAndCoeff;
    }
  
    public static double getAbsMean(float[] weights) {
      double sum = 0;
      for (float f : weights) {
        sum += Math.abs(f);
      }
      return  sum / weights.length;
    }
  }
  