com.opengamma.analytics.math.interpolation.RationalFunctionInterpolator1D.java Source code

Introduction

Here is the source code for com.opengamma.analytics.math.interpolation.RationalFunctionInterpolator1D.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;

import org.apache.commons.lang.Validate;

import com.opengamma.analytics.math.MathException;
import com.opengamma.analytics.math.interpolation.data.ArrayInterpolator1DDataBundle;
import com.opengamma.analytics.math.interpolation.data.Interpolator1DDataBundle;

/**
 * 
 */
public class RationalFunctionInterpolator1D extends Interpolator1D {
    private static final long serialVersionUID = 1L;
    private final int _degree;
    private final double _eps;

    public RationalFunctionInterpolator1D(final int degree, double eps) {
        Validate.isTrue(degree > 0, "Need a degree of at least one to perform rational function interpolation");
        _degree = degree;
        _eps = eps;
    }

    @Override
    public Double interpolate(final Interpolator1DDataBundle data, final Double value) {
        Validate.notNull(value, "value");
        Validate.notNull(data, "data bundle");
        final int m = _degree + 1;
        if (data.size() < m) {
            throw new IllegalArgumentException(
                    "Need at least " + (_degree + 1) + " data points to perform this interpolation");
        }
        final double[] xArray = data.getKeys();
        final double[] yArray = data.getValues();
        final int n = data.size() - 1;
        if (data.getLowerBoundIndex(value) == n) {
            return yArray[n];
        }
        double diff = Math.abs(value - xArray[0]);
        if (Math.abs(diff) < _eps) {
            return yArray[0];
        }
        double diff1;
        final double[] c = new double[m];
        final double[] d = new double[m];
        int ns = 0;
        for (int i = 0; i < m; i++) {
            diff1 = Math.abs(value - xArray[i]);
            if (diff < _eps) {
                return yArray[i];
            } else if (diff1 < diff) {
                ns = i;
                diff = diff1;
            }
            c[i] = yArray[i];
            d[i] = yArray[i] + _eps;
        }
        double y = yArray[ns--];
        double w, t, dd;
        for (int i = 1; i < m; i++) {
            for (int j = 0; j < m - i; j++) {
                w = c[j + 1] - d[j];
                diff = xArray[i + j] - value;
                t = (xArray[j] - value) * d[j] / diff;
                dd = t - c[j + 1];
                if (Math.abs(dd) < _eps) {
                    throw new MathException("Interpolating function has a pole at x = " + value);
                }
                dd = w / dd;
                d[j] = c[j + 1] * dd;
                c[j] = t * dd;
            }
            y += 2 * (ns + 1) < m - i ? c[ns + 1] : d[ns--];
        }
        return y;
    }

    @Override
    public Interpolator1DDataBundle getDataBundle(final double[] x, final double[] y) {
        return new ArrayInterpolator1DDataBundle(x, y);
    }

    @Override
    public Interpolator1DDataBundle getDataBundleFromSortedArrays(final double[] x, final double[] y) {
        return new ArrayInterpolator1DDataBundle(x, y, true);
    }

    @Override
    public boolean equals(final Object o) {
        if (o == null) {
            return false;
        }
        if (o == this) {
            return true;
        }
        if (!(o instanceof RationalFunctionInterpolator1D)) {
            return false;
        }
        final RationalFunctionInterpolator1D other = (RationalFunctionInterpolator1D) o;
        return _degree == other._degree;
    }

    @Override
    public int hashCode() {
        return getClass().hashCode() * 17 + _degree;
    }

    @Override
    public double[] getNodeSensitivitiesForValue(Interpolator1DDataBundle data, Double value) {
        return getFiniteDifferenceSensitivities(data, value);
    }
}