com.opengamma.analytics.math.interpolation.data.KrigingInterpolatorDataBundle.java Source code

Introduction

Here is the source code for com.opengamma.analytics.math.interpolation.data.KrigingInterpolatorDataBundle.java

Source

/**
 * 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;
    }

}