MATH-1173: Tricubic interpolation New class "TricubicInterpolatingFunction" to replace "TricubicSplineInterpolatingFunction".
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/fd7a4362 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/fd7a4362 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/fd7a4362 Branch: refs/heads/master Commit: fd7a4362caf2cd69368a38f2f71ffe19706e341c Parents: 17f52a2 Author: Gilles <er...@apache.org> Authored: Wed Dec 3 12:25:09 2014 +0100 Committer: Gilles <er...@apache.org> Committed: Wed Dec 3 12:25:09 2014 +0100 ---------------------------------------------------------------------- .../TricubicInterpolatingFunction.java | 510 +++++++++++++++ .../TricubicInterpolatingFunctionTest.java | 631 +++++++++++++++++++ 2 files changed, 1141 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/fd7a4362/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java new file mode 100644 index 0000000..b527996 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java @@ -0,0 +1,510 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.interpolation; + +import org.apache.commons.math3.analysis.TrivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NoDataException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.util.MathArrays; + +/** + * Function that implements the + * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation";> + * tricubic spline interpolation</a>, as proposed in + * <quote> + * Tricubic interpolation in three dimensions<br/> + * F. Lekien and J. Marsden<br/> + * <em>Int. J. Numer. Meth. Eng</em> 2005; <b>63</b>:455-471 + * </quote> + * + * @since 2.2 + * @deprecated in 3.4. + */ +@Deprecated +public class TricubicInterpolatingFunction + implements TrivariateFunction { + /** + * Matrix to compute the spline coefficients from the function values + * and function derivatives values + */ + private static final double[][] AINV = { + { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 }, + {-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,0,0,0,0,0,0,0,0,-2,-2,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 }, + { 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 }, + { -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 }, + { 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 }, + { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 }, + { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 }, + { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 }, + { 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 }, + { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 }, + { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 }, + { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 }, + { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 }, + { 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 }, + { 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 }, + { -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 }, + { 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 } + }; + + /** Samples x-coordinates */ + private final double[] xval; + /** Samples y-coordinates */ + private final double[] yval; + /** Samples z-coordinates */ + private final double[] zval; + /** Set of cubic splines pacthing the whole data grid */ + private final TricubicFunction[][][] splines; + + /** + * @param x Sample values of the x-coordinate, in increasing order. + * @param y Sample values of the y-coordinate, in increasing order. + * @param z Sample values of the y-coordinate, in increasing order. + * @param f Values of the function on every grid point. + * @param dFdX Values of the partial derivative of function with respect to x on every grid point. + * @param dFdY Values of the partial derivative of function with respect to y on every grid point. + * @param dFdZ Values of the partial derivative of function with respect to z on every grid point. + * @param d2FdXdY Values of the cross partial derivative of function on every grid point. + * @param d2FdXdZ Values of the cross partial derivative of function on every grid point. + * @param d2FdYdZ Values of the cross partial derivative of function on every grid point. + * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point. + * @throws NoDataException if any of the arrays has zero length. + * @throws DimensionMismatchException if the various arrays do not contain the expected number of elements. + * @throws NonMonotonicSequenceException if {@code x}, {@code y} or {@code z} are not strictly increasing. + */ + public TricubicInterpolatingFunction(double[] x, + double[] y, + double[] z, + double[][][] f, + double[][][] dFdX, + double[][][] dFdY, + double[][][] dFdZ, + double[][][] d2FdXdY, + double[][][] d2FdXdZ, + double[][][] d2FdYdZ, + double[][][] d3FdXdYdZ) + throws NoDataException, + DimensionMismatchException, + NonMonotonicSequenceException { + final int xLen = x.length; + final int yLen = y.length; + final int zLen = z.length; + + if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) { + throw new NoDataException(); + } + if (xLen != f.length) { + throw new DimensionMismatchException(xLen, f.length); + } + if (xLen != dFdX.length) { + throw new DimensionMismatchException(xLen, dFdX.length); + } + if (xLen != dFdY.length) { + throw new DimensionMismatchException(xLen, dFdY.length); + } + if (xLen != dFdZ.length) { + throw new DimensionMismatchException(xLen, dFdZ.length); + } + if (xLen != d2FdXdY.length) { + throw new DimensionMismatchException(xLen, d2FdXdY.length); + } + if (xLen != d2FdXdZ.length) { + throw new DimensionMismatchException(xLen, d2FdXdZ.length); + } + if (xLen != d2FdYdZ.length) { + throw new DimensionMismatchException(xLen, d2FdYdZ.length); + } + if (xLen != d3FdXdYdZ.length) { + throw new DimensionMismatchException(xLen, d3FdXdYdZ.length); + } + + MathArrays.checkOrder(x); + MathArrays.checkOrder(y); + MathArrays.checkOrder(z); + + xval = x.clone(); + yval = y.clone(); + zval = z.clone(); + + final int lastI = xLen - 1; + final int lastJ = yLen - 1; + final int lastK = zLen - 1; + splines = new TricubicFunction[lastI][lastJ][lastK]; + + for (int i = 0; i < lastI; i++) { + if (f[i].length != yLen) { + throw new DimensionMismatchException(f[i].length, yLen); + } + if (dFdX[i].length != yLen) { + throw new DimensionMismatchException(dFdX[i].length, yLen); + } + if (dFdY[i].length != yLen) { + throw new DimensionMismatchException(dFdY[i].length, yLen); + } + if (dFdZ[i].length != yLen) { + throw new DimensionMismatchException(dFdZ[i].length, yLen); + } + if (d2FdXdY[i].length != yLen) { + throw new DimensionMismatchException(d2FdXdY[i].length, yLen); + } + if (d2FdXdZ[i].length != yLen) { + throw new DimensionMismatchException(d2FdXdZ[i].length, yLen); + } + if (d2FdYdZ[i].length != yLen) { + throw new DimensionMismatchException(d2FdYdZ[i].length, yLen); + } + if (d3FdXdYdZ[i].length != yLen) { + throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen); + } + + final int ip1 = i + 1; + final double xR = xval[ip1] - xval[i]; + for (int j = 0; j < lastJ; j++) { + if (f[i][j].length != zLen) { + throw new DimensionMismatchException(f[i][j].length, zLen); + } + if (dFdX[i][j].length != zLen) { + throw new DimensionMismatchException(dFdX[i][j].length, zLen); + } + if (dFdY[i][j].length != zLen) { + throw new DimensionMismatchException(dFdY[i][j].length, zLen); + } + if (dFdZ[i][j].length != zLen) { + throw new DimensionMismatchException(dFdZ[i][j].length, zLen); + } + if (d2FdXdY[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen); + } + if (d2FdXdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen); + } + if (d2FdYdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen); + } + if (d3FdXdYdZ[i][j].length != zLen) { + throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen); + } + + final int jp1 = j + 1; + final double yR = yval[jp1] - yval[j]; + final double xRyR = xR * yR; + for (int k = 0; k < lastK; k++) { + final int kp1 = k + 1; + final double zR = zval[kp1] - zval[k]; + final double xRzR = xR * zR; + final double yRzR = yR * zR; + final double xRyRzR = xR * yRzR; + + final double[] beta = new double[] { + f[i][j][k], f[ip1][j][k], + f[i][jp1][k], f[ip1][jp1][k], + f[i][j][kp1], f[ip1][j][kp1], + f[i][jp1][kp1], f[ip1][jp1][kp1], + + dFdX[i][j][k] * xR, dFdX[ip1][j][k] * xR, + dFdX[i][jp1][k] * xR, dFdX[ip1][jp1][k] * xR, + dFdX[i][j][kp1] * xR, dFdX[ip1][j][kp1] * xR, + dFdX[i][jp1][kp1] * xR, dFdX[ip1][jp1][kp1] * xR, + + dFdY[i][j][k] * yR, dFdY[ip1][j][k] * yR, + dFdY[i][jp1][k] * yR, dFdY[ip1][jp1][k] * yR, + dFdY[i][j][kp1] * yR, dFdY[ip1][j][kp1] * yR, + dFdY[i][jp1][kp1] * yR, dFdY[ip1][jp1][kp1] * yR, + + dFdZ[i][j][k] * zR, dFdZ[ip1][j][k] * zR, + dFdZ[i][jp1][k] * zR, dFdZ[ip1][jp1][k] * zR, + dFdZ[i][j][kp1] * zR, dFdZ[ip1][j][kp1] * zR, + dFdZ[i][jp1][kp1] * zR, dFdZ[ip1][jp1][kp1] * zR, + + d2FdXdY[i][j][k] * xRyR, d2FdXdY[ip1][j][k] * xRyR, + d2FdXdY[i][jp1][k] * xRyR, d2FdXdY[ip1][jp1][k] * xRyR, + d2FdXdY[i][j][kp1] * xRyR, d2FdXdY[ip1][j][kp1] * xRyR, + d2FdXdY[i][jp1][kp1] * xRyR, d2FdXdY[ip1][jp1][kp1] * xRyR, + + d2FdXdZ[i][j][k] * xRzR, d2FdXdZ[ip1][j][k] * xRzR, + d2FdXdZ[i][jp1][k] * xRzR, d2FdXdZ[ip1][jp1][k] * xRzR, + d2FdXdZ[i][j][kp1] * xRzR, d2FdXdZ[ip1][j][kp1] * xRzR, + d2FdXdZ[i][jp1][kp1] * xRzR, d2FdXdZ[ip1][jp1][kp1] * xRzR, + + d2FdYdZ[i][j][k] * yRzR, d2FdYdZ[ip1][j][k] * yRzR, + d2FdYdZ[i][jp1][k] * yRzR, d2FdYdZ[ip1][jp1][k] * yRzR, + d2FdYdZ[i][j][kp1] * yRzR, d2FdYdZ[ip1][j][kp1] * yRzR, + d2FdYdZ[i][jp1][kp1] * yRzR, d2FdYdZ[ip1][jp1][kp1] * yRzR, + + d3FdXdYdZ[i][j][k] * xRyRzR, d3FdXdYdZ[ip1][j][k] * xRyRzR, + d3FdXdYdZ[i][jp1][k] * xRyRzR, d3FdXdYdZ[ip1][jp1][k] * xRyRzR, + d3FdXdYdZ[i][j][kp1] * xRyRzR, d3FdXdYdZ[ip1][j][kp1] * xRyRzR, + d3FdXdYdZ[i][jp1][kp1] * xRyRzR, d3FdXdYdZ[ip1][jp1][kp1] * xRyRzR, + }; + + splines[i][j][k] = new TricubicFunction(computeCoefficients(beta)); + } + } + } + } + + /** + * {@inheritDoc} + * + * @throws OutOfRangeException if any of the variables is outside its interpolation range. + */ + public double value(double x, double y, double z) + throws OutOfRangeException { + final int i = searchIndex(x, xval); + if (i == -1) { + throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); + } + final int j = searchIndex(y, yval); + if (j == -1) { + throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); + } + final int k = searchIndex(z, zval); + if (k == -1) { + throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]); + } + + final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); + final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); + final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]); + + return splines[i][j][k].value(xN, yN, zN); + } + + /** + * Indicates whether a point is within the interpolation range. + * + * @param x First coordinate. + * @param y Second coordinate. + * @param z Third coordinate. + * @return {@code true} if (x, y, z) is a valid point. + */ + public boolean isValidPoint(double x, double y, double z) { + if (x < xval[0] || + x > xval[xval.length - 1] || + y < yval[0] || + y > yval[yval.length - 1] || + z < zval[0] || + z > zval[zval.length - 1]) { + return false; + } else { + return true; + } + } + + /** + * @param c Coordinate. + * @param val Coordinate samples. + * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1} + * if {@code c} is out of the range defined by the end values of {@code val}. + */ + private int searchIndex(double c, double[] val) { + if (c < val[0]) { + return -1; + } + + final int max = val.length; + for (int i = 1; i < max; i++) { + if (c <= val[i]) { + return i - 1; + } + } + + return -1; + } + + /** + * Compute the spline coefficients from the list of function values and + * function partial derivatives values at the four corners of a grid + * element. They must be specified in the following order: + * <ul> + * <li>f(0,0,0)</li> + * <li>f(1,0,0)</li> + * <li>f(0,1,0)</li> + * <li>f(1,1,0)</li> + * <li>f(0,0,1)</li> + * <li>f(1,0,1)</li> + * <li>f(0,1,1)</li> + * <li>f(1,1,1)</li> + * + * <li>f<sub>x</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>x</sub>(1,1,1)</li> + * + * <li>f<sub>y</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>y</sub>(1,1,1)</li> + * + * <li>f<sub>z</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>z</sub>(1,1,1)</li> + * + * <li>f<sub>xy</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xy</sub>(1,1,1)</li> + * + * <li>f<sub>xz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xz</sub>(1,1,1)</li> + * + * <li>f<sub>yz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>yz</sub>(1,1,1)</li> + * + * <li>f<sub>xyz</sub>(0,0,0)</li> + * <li>... <em>(same order as above)</em></li> + * <li>f<sub>xyz</sub>(1,1,1)</li> + * </ul> + * where the subscripts indicate the partial derivative with respect to + * the corresponding variable(s). + * + * @param beta List of function values and function partial derivatives values. + * @return the spline coefficients. + */ + private double[] computeCoefficients(double[] beta) { + final int sz = 64; + final double[] a = new double[sz]; + + for (int i = 0; i < sz; i++) { + double result = 0; + final double[] row = AINV[i]; + for (int j = 0; j < sz; j++) { + result += row[j] * beta[j]; + } + a[i] = result; + } + + return a; + } +} + +/** + * 3D-spline function. + * + */ +class TricubicFunction + implements TrivariateFunction { + /** Number of points. */ + private static final short N = 4; + /** Coefficients */ + private final double[][][] a = new double[N][N][N]; + + /** + * @param aV List of spline coefficients. + */ + public TricubicFunction(double[] aV) { + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + a[i][j][k] = aV[i + N * (j + N * k)]; + } + } + } + } + + /** + * @param x x-coordinate of the interpolation point. + * @param y y-coordinate of the interpolation point. + * @param z z-coordinate of the interpolation point. + * @return the interpolated value. + * @throws OutOfRangeException if {@code x}, {@code y} or + * {@code z} are not in the interval {@code [0, 1]}. + */ + public double value(double x, double y, double z) + throws OutOfRangeException { + if (x < 0 || x > 1) { + throw new OutOfRangeException(x, 0, 1); + } + if (y < 0 || y > 1) { + throw new OutOfRangeException(y, 0, 1); + } + if (z < 0 || z > 1) { + throw new OutOfRangeException(z, 0, 1); + } + + final double x2 = x * x; + final double x3 = x2 * x; + final double[] pX = { 1, x, x2, x3 }; + + final double y2 = y * y; + final double y3 = y2 * y; + final double[] pY = { 1, y, y2, y3 }; + + final double z2 = z * z; + final double z3 = z2 * z; + final double[] pZ = { 1, z, z2, z3 }; + + double result = 0; + for (int i = 0; i < N; i++) { + for (int j = 0; j < N; j++) { + for (int k = 0; k < N; k++) { + result += a[i][j][k] * pX[i] * pY[j] * pZ[k]; + } + } + } + + return result; + } +} http://git-wip-us.apache.org/repos/asf/commons-math/blob/fd7a4362/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunctionTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunctionTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunctionTest.java new file mode 100644 index 0000000..9760555 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunctionTest.java @@ -0,0 +1,631 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.analysis.interpolation; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.analysis.TrivariateFunction; +import org.apache.commons.math3.distribution.UniformRealDistribution; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well19937c; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; +import org.junit.Assert; +import org.junit.Test; + +/** + * Test case for the bicubic function. + * + */ +public final class TricubicInterpolatingFunctionTest { + /** + * Test preconditions. + */ + @Test + public void testPreconditions() { + double[] xval = new double[] {3, 4, 5, 6.5}; + double[] yval = new double[] {-4, -3, -1, 2.5}; + double[] zval = new double[] {-12, -8, -5.5, -3, 0, 2.5}; + double[][][] fval = new double[xval.length][yval.length][zval.length]; + + @SuppressWarnings("unused") + TrivariateFunction tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, fval, fval, fval); + + double[] wxval = new double[] {3, 2, 5, 6.5}; + try { + tcf = new TricubicInterpolatingFunction(wxval, yval, zval, + fval, fval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected + } + double[] wyval = new double[] {-4, -1, -1, 2.5}; + try { + tcf = new TricubicInterpolatingFunction(xval, wyval, zval, + fval, fval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected + } + double[] wzval = new double[] {-12, -8, -9, -3, 0, 2.5}; + try { + tcf = new TricubicInterpolatingFunction(xval, yval, wzval, + fval, fval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (MathIllegalArgumentException e) { + // Expected + } + double[][][] wfval = new double[xval.length - 1][yval.length - 1][zval.length]; + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + wfval, fval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, wfval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, wfval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, wfval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + wfval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, wfval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, fval, wfval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, fval, fval, wfval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + wfval = new double[xval.length][yval.length - 1][zval.length]; + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + wfval, fval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, wfval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, wfval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, wfval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + wfval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, wfval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, fval, wfval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, fval, fval, wfval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + wfval = new double[xval.length][yval.length][zval.length - 1]; + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + wfval, fval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, wfval, fval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, wfval, fval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, wfval, + fval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + wfval, fval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, wfval, fval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, fval, wfval, fval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + try { + tcf = new TricubicInterpolatingFunction(xval, yval, zval, + fval, fval, fval, fval, + fval, fval, fval, wfval); + Assert.fail("an exception should have been thrown"); + } catch (DimensionMismatchException e) { + // Expected + } + } + + /** + * @param minimumX Lower bound of interpolation range along the x-coordinate. + * @param maximumX Higher bound of interpolation range along the x-coordinate. + * @param minimumY Lower bound of interpolation range along the y-coordinate. + * @param maximumY Higher bound of interpolation range along the y-coordinate. + * @param minimumZ Lower bound of interpolation range along the z-coordinate. + * @param maximumZ Higher bound of interpolation range along the z-coordinate. + * @param numberOfElements Number of data points (along each dimension). + * @param numberOfSamples Number of test points. + * @param f Function to test. + * @param dfdx Partial derivative w.r.t. x of the function to test. + * @param dfdy Partial derivative w.r.t. y of the function to test. + * @param dfdz Partial derivative w.r.t. z of the function to test. + * @param d2fdxdy Second partial cross-derivative w.r.t x and y of the function to test. + * @param d2fdxdz Second partial cross-derivative w.r.t x and z of the function to test. + * @param d2fdydz Second partial cross-derivative w.r.t y and z of the function to test. + * @param d3fdxdy Third partial cross-derivative w.r.t x, y and z of the function to test. + * @param meanRelativeTolerance Allowed average error (mean error on all interpolated values). + * @param maxRelativeTolerance Allowed error on each interpolated value. + * @param onDataMaxAbsoluteTolerance Allowed error on a data point. + */ + private void testInterpolation(double minimumX, + double maximumX, + double minimumY, + double maximumY, + double minimumZ, + double maximumZ, + int numberOfElements, + int numberOfSamples, + TrivariateFunction f, + TrivariateFunction dfdx, + TrivariateFunction dfdy, + TrivariateFunction dfdz, + TrivariateFunction d2fdxdy, + TrivariateFunction d2fdxdz, + TrivariateFunction d2fdydz, + TrivariateFunction d3fdxdydz, + double meanRelativeTolerance, + double maxRelativeTolerance, + double onDataMaxAbsoluteTolerance, + boolean print) { + double expected; + double actual; + double currentX; + double currentY; + double currentZ; + final double deltaX = (maximumX - minimumX) / numberOfElements; + final double deltaY = (maximumY - minimumY) / numberOfElements; + final double deltaZ = (maximumZ - minimumZ) / numberOfElements; + final double[] xValues = new double[numberOfElements]; + final double[] yValues = new double[numberOfElements]; + final double[] zValues = new double[numberOfElements]; + final double[][][] fValues = new double[numberOfElements][numberOfElements][numberOfElements]; + final double[][][] dfdxValues = new double[numberOfElements][numberOfElements][numberOfElements]; + final double[][][] dfdyValues = new double[numberOfElements][numberOfElements][numberOfElements]; + final double[][][] dfdzValues = new double[numberOfElements][numberOfElements][numberOfElements]; + final double[][][] d2fdxdyValues = new double[numberOfElements][numberOfElements][numberOfElements]; + final double[][][] d2fdxdzValues = new double[numberOfElements][numberOfElements][numberOfElements]; + final double[][][] d2fdydzValues = new double[numberOfElements][numberOfElements][numberOfElements]; + final double[][][] d3fdxdydzValues = new double[numberOfElements][numberOfElements][numberOfElements]; + + for (int i = 0; i < numberOfElements; i++) { + xValues[i] = minimumX + deltaX * i; + final double x = xValues[i]; + for (int j = 0; j < numberOfElements; j++) { + yValues[j] = minimumY + deltaY * j; + final double y = yValues[j]; + for (int k = 0; k < numberOfElements; k++) { + zValues[k] = minimumZ + deltaZ * k; + final double z = zValues[k]; + fValues[i][j][k] = f.value(x, y, z); + dfdxValues[i][j][k] = dfdx.value(x, y, z); + dfdyValues[i][j][k] = dfdy.value(x, y, z); + dfdzValues[i][j][k] = dfdz.value(x, y, z); + d2fdxdyValues[i][j][k] = d2fdxdy.value(x, y, z); + d2fdxdzValues[i][j][k] = d2fdxdz.value(x, y, z); + d2fdydzValues[i][j][k] = d2fdydz.value(x, y, z); + d3fdxdydzValues[i][j][k] = d3fdxdydz.value(x, y, z); + } + } + } + + final TrivariateFunction interpolation + = new TricubicInterpolatingFunction(xValues, + yValues, + zValues, + fValues, + dfdxValues, + dfdyValues, + dfdzValues, + d2fdxdyValues, + d2fdxdzValues, + d2fdydzValues, + d3fdxdydzValues); + + for (int i = 0; i < numberOfElements; i++) { + currentX = xValues[i]; + for (int j = 0; j < numberOfElements; j++) { + currentY = yValues[j]; + for (int k = 0; k < numberOfElements; k++) { + currentZ = zValues[k]; + expected = f.value(currentX, currentY, currentZ); + actual = interpolation.value(currentX, currentY, currentZ); + Assert.assertTrue("On data point: " + expected + " != " + actual, + Precision.equals(expected, actual, onDataMaxAbsoluteTolerance)); + } + } + } + + final RandomGenerator rng = new Well19937c(1234567L); + final UniformRealDistribution distX = new UniformRealDistribution(rng, xValues[0], xValues[xValues.length - 1]); + final UniformRealDistribution distY = new UniformRealDistribution(rng, yValues[0], yValues[yValues.length - 1]); + final UniformRealDistribution distZ = new UniformRealDistribution(rng, zValues[0], zValues[zValues.length - 1]); + + double sumError = 0; + for (int i = 0; i < numberOfSamples; i++) { + currentX = distX.sample(); + currentY = distY.sample(); + currentZ = distZ.sample(); + expected = f.value(currentX, currentY, currentZ); + + actual = interpolation.value(currentX, currentY, currentZ); + final double relativeError = FastMath.abs(actual - expected) / FastMath.max(FastMath.abs(actual), FastMath.abs(expected)); + sumError += relativeError; + + if (print) { + System.out.println(currentX + " " + currentY + " " + currentZ + + " " + actual + + " " + expected + + " " + (expected - actual)); + } + + Assert.assertEquals(0, relativeError, maxRelativeTolerance); + } + + final double meanError = sumError / numberOfSamples; + Assert.assertEquals(0, meanError, meanRelativeTolerance); + } + + /** + * Test for a plane. + * <p> + * f(x, y, z) = 2 x - 3 y - 4 z + 5 + * </p> + */ + @Test + public void testPlane() { + final TrivariateFunction f = new TrivariateFunction() { + public double value(double x, double y, double z) { + return 2 * x - 3 * y - 4 * z + 5; + } + }; + + final TrivariateFunction dfdx = new TrivariateFunction() { + public double value(double x, double y, double z) { + return 2; + } + }; + + final TrivariateFunction dfdy = new TrivariateFunction() { + public double value(double x, double y, double z) { + return -3; + } + }; + + final TrivariateFunction dfdz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return -4; + } + }; + + final TrivariateFunction zero = new TrivariateFunction() { + public double value(double x, double y, double z) { + return 0; + } + }; + + testInterpolation(-10, 3, + 4.5, 6, + -150, -117, + 7, + 1000, + f, + dfdx, + dfdy, + dfdz, + zero, + zero, + zero, + zero, + 1e-12, + 1e-11, + 1e-10, + false); + } + + /** + * Test for a quadric. + * <p> + * f(x, y, z) = 2 x<sup>2</sup> - 3 y<sup>2</sup> - 4 z<sup>2</sup> + 5 x y + 6 x z - 2 y z + 3 + * </p> + */ + @Test + public void testQuadric() { + final TrivariateFunction f = new TrivariateFunction() { + public double value(double x, double y, double z) { + return 2 * x * x - 3 * y * y - 4 * z * z + 5 * x * y + 6 * x * z - 2 * y * z + 3; + } + }; + + final TrivariateFunction dfdx = new TrivariateFunction() { + public double value(double x, double y, double z) { + return 4 * x + 5 * y + 6 * z; + } + }; + + final TrivariateFunction dfdy = new TrivariateFunction() { + public double value(double x, double y, double z) { + return -6 * y + 5 * x - 2 * z; + } + }; + + final TrivariateFunction dfdz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return -8 * z + 6 * x - 2 * y; + } + }; + + final TrivariateFunction d2fdxdy = new TrivariateFunction() { + public double value(double x, double y, double z) { + return 5; + } + }; + + final TrivariateFunction d2fdxdz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return 6; + } + }; + + final TrivariateFunction d2fdydz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return -2; + } + }; + + final TrivariateFunction d3fdxdydz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return 0; + } + }; + + testInterpolation(-10, 3, + 4.5, 6, + -150, -117, + 7, + 5000, + f, + dfdx, + dfdy, + dfdz, + d2fdxdy, + d2fdxdz, + d2fdydz, + d3fdxdydz, + 1e-12, + 1e-11, + 1e-8, + false); + } + + /** + * Wave. + * <p> + * f(x, y, z) = a cos (ω z - k<sub>x</sub> x - k<sub>y</sub> y) + * </p> + * with a = 5, ω = 0.3, k<sub>x</sub> = 0.8, k<sub>y</sub> = 1. + */ + @Test + public void testWave() { + final double a = 5; + final double omega = 0.3; + final double kx = 0.8; + final double ky = 1; + + final TrivariateFunction arg = new TrivariateFunction() { + public double value(double x, double y, double z) { + return omega * z - kx * x - ky * y; + } + }; + + final TrivariateFunction f = new TrivariateFunction() { + public double value(double x, double y, double z) { + return a * FastMath.cos(arg.value(x, y, z)); + } + }; + + final TrivariateFunction dfdx = new TrivariateFunction() { + public double value(double x, double y, double z) { + return kx * a * FastMath.sin(arg.value(x, y, z)); + } + }; + + final TrivariateFunction dfdy = new TrivariateFunction() { + public double value(double x, double y, double z) { + return ky * a * FastMath.sin(arg.value(x, y, z)); + } + }; + + final TrivariateFunction dfdz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return -omega * a * FastMath.sin(arg.value(x, y, z)); + } + }; + + final TrivariateFunction d2fdxdy = new TrivariateFunction() { + public double value(double x, double y, double z) { + return -ky * kx * a * FastMath.cos(arg.value(x, y, z)); + } + }; + + final TrivariateFunction d2fdxdz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return omega * kx * a * FastMath.cos(arg.value(x, y, z)); + } + }; + + final TrivariateFunction d2fdydz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return omega * ky * a * FastMath.cos(arg.value(x, y, z)); + } + }; + + final TrivariateFunction d3fdxdydz = new TrivariateFunction() { + public double value(double x, double y, double z) { + return omega * ky * kx * a * FastMath.sin(arg.value(x, y, z)); + } + }; + + testInterpolation(-10, 3, + 4.5, 6, + -150, -117, + 30, + 5000, + f, + dfdx, + dfdy, + dfdz, + d2fdxdy, + d2fdxdz, + d2fdydz, + d3fdxdydz, + 1e-3, + 1e-2, + 1e-12, + false); + } +}
