001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math3.analysis.interpolation; 018 019import org.apache.commons.math3.analysis.UnivariateFunction; 020import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; 021import org.apache.commons.math3.exception.DimensionMismatchException; 022import org.apache.commons.math3.exception.NoDataException; 023import org.apache.commons.math3.exception.NonMonotonicSequenceException; 024import org.apache.commons.math3.exception.NumberIsTooSmallException; 025import org.apache.commons.math3.util.MathArrays; 026 027/** 028 * Generates a bicubic interpolating function. Due to numerical accuracy issues this should not 029 * be used. 030 * 031 * @since 2.2 032 * @deprecated as of 3.4 replaced by {@link org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolator} 033 */ 034@Deprecated 035public class BicubicSplineInterpolator 036 implements BivariateGridInterpolator { 037 /** Whether to initialize internal data used to compute the analytical 038 derivatives of the splines. */ 039 private final boolean initializeDerivatives; 040 041 /** 042 * Default constructor. 043 * The argument {@link #BicubicSplineInterpolator(boolean) initializeDerivatives} 044 * is set to {@code false}. 045 */ 046 public BicubicSplineInterpolator() { 047 this(false); 048 } 049 050 /** 051 * Creates an interpolator. 052 * 053 * @param initializeDerivatives Whether to initialize the internal data 054 * needed for calling any of the methods that compute the partial derivatives 055 * of the {@link BicubicSplineInterpolatingFunction function} returned from 056 * the call to {@link #interpolate(double[],double[],double[][]) interpolate}. 057 */ 058 public BicubicSplineInterpolator(boolean initializeDerivatives) { 059 this.initializeDerivatives = initializeDerivatives; 060 } 061 062 /** 063 * {@inheritDoc} 064 */ 065 public BicubicSplineInterpolatingFunction interpolate(final double[] xval, 066 final double[] yval, 067 final double[][] fval) 068 throws NoDataException, DimensionMismatchException, 069 NonMonotonicSequenceException, NumberIsTooSmallException { 070 if (xval.length == 0 || yval.length == 0 || fval.length == 0) { 071 throw new NoDataException(); 072 } 073 if (xval.length != fval.length) { 074 throw new DimensionMismatchException(xval.length, fval.length); 075 } 076 077 MathArrays.checkOrder(xval); 078 MathArrays.checkOrder(yval); 079 080 final int xLen = xval.length; 081 final int yLen = yval.length; 082 083 // Samples (first index is y-coordinate, i.e. subarray variable is x) 084 // 0 <= i < xval.length 085 // 0 <= j < yval.length 086 // fX[j][i] = f(xval[i], yval[j]) 087 final double[][] fX = new double[yLen][xLen]; 088 for (int i = 0; i < xLen; i++) { 089 if (fval[i].length != yLen) { 090 throw new DimensionMismatchException(fval[i].length, yLen); 091 } 092 093 for (int j = 0; j < yLen; j++) { 094 fX[j][i] = fval[i][j]; 095 } 096 } 097 098 final SplineInterpolator spInterpolator = new SplineInterpolator(); 099 100 // For each line y[j] (0 <= j < yLen), construct a 1D spline with 101 // respect to variable x 102 final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; 103 for (int j = 0; j < yLen; j++) { 104 ySplineX[j] = spInterpolator.interpolate(xval, fX[j]); 105 } 106 107 // For each line x[i] (0 <= i < xLen), construct a 1D spline with 108 // respect to variable y generated by array fY_1[i] 109 final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen]; 110 for (int i = 0; i < xLen; i++) { 111 xSplineY[i] = spInterpolator.interpolate(yval, fval[i]); 112 } 113 114 // Partial derivatives with respect to x at the grid knots 115 final double[][] dFdX = new double[xLen][yLen]; 116 for (int j = 0; j < yLen; j++) { 117 final UnivariateFunction f = ySplineX[j].derivative(); 118 for (int i = 0; i < xLen; i++) { 119 dFdX[i][j] = f.value(xval[i]); 120 } 121 } 122 123 // Partial derivatives with respect to y at the grid knots 124 final double[][] dFdY = new double[xLen][yLen]; 125 for (int i = 0; i < xLen; i++) { 126 final UnivariateFunction f = xSplineY[i].derivative(); 127 for (int j = 0; j < yLen; j++) { 128 dFdY[i][j] = f.value(yval[j]); 129 } 130 } 131 132 // Cross partial derivatives 133 final double[][] d2FdXdY = new double[xLen][yLen]; 134 for (int i = 0; i < xLen ; i++) { 135 final int nI = nextIndex(i, xLen); 136 final int pI = previousIndex(i); 137 for (int j = 0; j < yLen; j++) { 138 final int nJ = nextIndex(j, yLen); 139 final int pJ = previousIndex(j); 140 d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - 141 fval[pI][nJ] + fval[pI][pJ]) / 142 ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])); 143 } 144 } 145 146 // Create the interpolating splines 147 return new BicubicSplineInterpolatingFunction(xval, yval, fval, 148 dFdX, dFdY, d2FdXdY, 149 initializeDerivatives); 150 } 151 152 /** 153 * Computes the next index of an array, clipping if necessary. 154 * It is assumed (but not checked) that {@code i >= 0}. 155 * 156 * @param i Index. 157 * @param max Upper limit of the array. 158 * @return the next index. 159 */ 160 private int nextIndex(int i, int max) { 161 final int index = i + 1; 162 return index < max ? index : index - 1; 163 } 164 /** 165 * Computes the previous index of an array, clipping if necessary. 166 * It is assumed (but not checked) that {@code i} is smaller than the size 167 * of the array. 168 * 169 * @param i Index. 170 * @return the previous index. 171 */ 172 private int previousIndex(int i) { 173 final int index = i - 1; 174 return index >= 0 ? index : 0; 175 } 176}