KrigingInterpolatorDataBundle.java

/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 * 
 * Please see distribution for license.
 */
package com.opengamma.analytics.math.interpolation.data;

import java.util.Arrays;
import java.util.List;

import org.apache.commons.lang.ObjectUtils;
import org.apache.commons.lang.Validate;

import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.interpolation.DistanceCalculator;
import com.opengamma.analytics.math.linearalgebra.Decomposition;
import com.opengamma.analytics.math.linearalgebra.DecompositionFactory;
import com.opengamma.analytics.math.linearalgebra.DecompositionResult;
import com.opengamma.analytics.math.matrix.DoubleMatrix1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix2D;
import com.opengamma.util.tuple.Pair;

/**
 * 
 */
public class KrigingInterpolatorDataBundle extends InterpolatorNDDataBundle {
  private final Decomposition<?> _decomp = DecompositionFactory.LU_COMMONS;
  private final Function1D<Double, Double> _variogram;
  private final double[] _weights;

  /**
   * @param data The data
   * @param beta The beta
   */
  public KrigingInterpolatorDataBundle(final List<Pair<double[], Double>> data, final double beta) {
    super(data);
    Validate.isTrue(beta >= 1 && beta < 2, "Beta was not in acceptable range (1 <= beta < 2");
    _variogram = calculateVariogram(data, beta);
    _weights = calculateWeights(data, _variogram);
  }

  /**
   * Gets the variogram field.
   * @return the variogram
   */
  public Function1D<Double, Double> getVariogram() {
    return _variogram;
  }

  /**
   * Gets the weights field.
   * @return the weights
   */
  public double[] getWeights() {
    return _weights;
  }

  private Function1D<Double, Double> calculateVariogram(final List<Pair<double[], Double>> data, final double beta) {
    final int n = data.size();
    double sum = 0.0;
    double normSum = 0.0;
    double[] x1, x2;
    double y1, y2;
    Pair<double[], Double> dataPoint;
    double rBeta;
    for (int i = 0; i < n; i++) {
      dataPoint = data.get(i);
      x1 = dataPoint.getFirst();
      y1 = dataPoint.getSecond();
      for (int j = i + 1; j < n; j++) {
        dataPoint = data.get(j);
        x2 = dataPoint.getFirst();
        y2 = dataPoint.getSecond();
        rBeta = Math.pow(DistanceCalculator.getDistance(x1, x2), beta);
        sum += (y1 - y2) * (y1 - y2) * rBeta;
        normSum += rBeta * rBeta;
      }
    }

    final double alpha = sum / 2.0 / normSum;

    return new Function1D<Double, Double>() {

      @Override
      public Double evaluate(final Double x) {
        return alpha * Math.pow(x, beta);
      }
    };
  }

  private double[] calculateWeights(final List<Pair<double[], Double>> data, final Function1D<Double, Double> variogram) {
    final int n = data.size();
    final double[] y = new double[n + 1];
    final double[][] v = new double[n + 1][n + 1];
    Pair<double[], Double> dataPoint;
    double[] x1, x2;

    for (int i = 0; i < n; i++) {
      dataPoint = data.get(i);
      x1 = dataPoint.getFirst();
      y[i] = dataPoint.getSecond();
      for (int j = i + 1; j < n; j++) {
        dataPoint = data.get(j);
        x2 = dataPoint.getFirst();
        final double temp = variogram.evaluate(DistanceCalculator.getDistance(x1, x2));
        v[i][j] = temp;
        v[j][i] = temp;
      }
      v[i][n] = 1;
      v[n][i] = 1;
    }
    v[n][n] = 0;
    y[n] = 0;

    double[] res;
    try {
      res = solve(v, y, _decomp);
    } catch (final IllegalArgumentException e) {
      final Decomposition<?> decomp = DecompositionFactory.SV_COMMONS;
      res = solve(v, y, decomp);
    }

    return res;
  }

  private double[] solve(final double[][] v, final double[] y, final Decomposition<?> decomp) {
    final DecompositionResult decompRes = decomp.evaluate(new DoubleMatrix2D(v));
    final DoubleMatrix1D res = decompRes.solve(new DoubleMatrix1D(y));
    return res.getData();
  }

  @Override
  public int hashCode() {
    final int prime = 31;
    int result = super.hashCode();
    result = prime * result + _decomp.hashCode();
    result = prime * result + _variogram.hashCode();
    result = prime * result + Arrays.hashCode(_weights);
    return result;
  }

  @Override
  public boolean equals(final Object obj) {
    if (this == obj) {
      return true;
    }
    if (!super.equals(obj)) {
      return false;
    }
    if (!(obj instanceof KrigingInterpolatorDataBundle)) {
      return false;
    }
    final KrigingInterpolatorDataBundle other = (KrigingInterpolatorDataBundle) obj;
    if (!ObjectUtils.equals(_decomp, other._decomp)) {
      return false;
    }
    if (!ObjectUtils.equals(_variogram, other._variogram)) {
      return false;
    }
    if (!Arrays.equals(_weights, other._weights)) {
      return false;
    }
    return true;
  }

}