From e4a1d4679443e55d787d5ba5ce90396b1565cdff Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Sat, 14 Sep 2024 14:12:50 +0200 Subject: [PATCH] Added FieldBivariateGridInterpolator, FieldBilinearInterpolator and FieldBilinearInterpolatingFunction. Fixes #354 --- hipparchus-core/src/changes/changes.xml | 4 +- .../FieldBilinearInterpolatingFunction.java | 136 +++++++++++ .../FieldBilinearInterpolator.java | 48 ++++ .../FieldBivariateGridInterpolator.java | 48 ++++ .../analysis/interpolation/FieldGridAxis.java | 207 +++++++++++++++++ .../analysis/interpolation/GridAxis.java | 4 +- ...ieldBilinearInterpolatingFunctionTest.java | 215 ++++++++++++++++++ 7 files changed, 658 insertions(+), 4 deletions(-) create mode 100644 hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolatingFunction.java create mode 100644 hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolator.java create mode 100644 hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBivariateGridInterpolator.java create mode 100644 hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldGridAxis.java create mode 100644 hipparchus-core/src/test/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolatingFunctionTest.java diff --git a/hipparchus-core/src/changes/changes.xml b/hipparchus-core/src/changes/changes.xml index 29f84f585..d1f276e3d 100644 --- a/hipparchus-core/src/changes/changes.xml +++ b/hipparchus-core/src/changes/changes.xml @@ -50,8 +50,8 @@ If the output is not quite correct, check for invisible trailing spaces! - - Fixed {Field}Complex.atan for real values. + + Added FieldBivariateGridInterpolator, FieldBilinearInterpolator and FieldBilinearInterpolatingFunction. Added getAddendum() to CalculusFieldElement interface. diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolatingFunction.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolatingFunction.java new file mode 100644 index 000000000..d0326ed86 --- /dev/null +++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolatingFunction.java @@ -0,0 +1,136 @@ +/* + * Licensed to the Hipparchus project under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The Hipparchus project 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 + * + * https://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.hipparchus.analysis.interpolation; + +import org.hipparchus.CalculusFieldElement; +import org.hipparchus.Field; +import org.hipparchus.analysis.CalculusFieldBivariateFunction; +import org.hipparchus.exception.MathIllegalArgumentException; +import org.hipparchus.util.MathArrays; + +/** + * Interpolate grid data using bi-linear interpolation. + *

+ * This interpolator is thread-safe. + *

+ * @param Type of the field elements. + * @since 4.0 + */ +public class FieldBilinearInterpolatingFunction> + implements CalculusFieldBivariateFunction { + + /** Grid along the x axis. */ + private final FieldGridAxis xGrid; + + /** Grid along the y axis. */ + private final FieldGridAxis yGrid; + + /** Grid size along the y axis. */ + private final int ySize; + + /** Values of the interpolation points on all the grid knots (in a flatten array). */ + private final T[] fVal; + + /** Simple constructor. + * @param xVal All the x-coordinates of the interpolation points, sorted + * in increasing order. + * @param yVal All the y-coordinates of the interpolation points, sorted + * in increasing order. + * @param fVal The values of the interpolation points on all the grid knots: + * {@code fVal[i][j] = f(xVal[i], yVal[j])}. + * @exception MathIllegalArgumentException if grid size is smaller than 2 + * or if the grid is not sorted in strict increasing order + */ + public FieldBilinearInterpolatingFunction(final T[] xVal, final T[] yVal, final T[][] fVal) + throws MathIllegalArgumentException { + final Field field = fVal[0][0].getField(); + this.xGrid = new FieldGridAxis<>(xVal, 2); + this.yGrid = new FieldGridAxis<>(yVal, 2); + this.ySize = yVal.length; + this.fVal = MathArrays.buildArray(field, xVal.length * ySize); + int k = 0; + for (int i = 0; i < xVal.length; ++i) { + final T[] fi = fVal[i]; + for (int j = 0; j < ySize; ++j) { + this.fVal[k++] = fi[j]; + } + } + } + + /** Get the lowest grid x coordinate. + * @return lowest grid x coordinate + */ + public T getXInf() { + return xGrid.node(0); + } + + /** Get the highest grid x coordinate. + * @return highest grid x coordinate + */ + public T getXSup() { + return xGrid.node(xGrid.size() - 1); + } + + /** Get the lowest grid y coordinate. + * @return lowest grid y coordinate + */ + public T getYInf() { + return yGrid.node(0); + } + + /** Get the highest grid y coordinate. + * @return highest grid y coordinate + */ + public T getYSup() { + return yGrid.node(yGrid.size() - 1); + } + + /** {@inheritDoc} */ + @Override + public T value(T x, T y) { + + // get the interpolation nodes + final int i = xGrid.interpolationIndex(x); + final int j = yGrid.interpolationIndex(y); + final T x0 = xGrid.node(i); + final T x1 = xGrid.node(i + 1); + final T y0 = yGrid.node(j); + final T y1 = yGrid.node(j + 1); + + // get the function values at interpolation nodes + final int k0 = i * ySize + j; + final int k1 = k0 + ySize; + final T z00 = fVal[k0]; + final T z01 = fVal[k0 + 1]; + final T z10 = fVal[k1]; + final T z11 = fVal[k1 + 1]; + + // interpolate + final T dx0 = x.subtract(x0); + final T mdx1 = x.subtract(x1); + final T dx10 = x1.subtract(x0); + final T dy0 = y.subtract(y0); + final T mdy1 = y.subtract(y1); + final T dy10 = y1.subtract(y0); + + return dy0.multiply(z11).subtract(mdy1.multiply(z10)).multiply(dx0). + subtract(dy0.multiply(z01).subtract(mdy1.multiply(z00)).multiply(mdx1)). + divide(dx10.multiply(dy10)); + + } + +} diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolator.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolator.java new file mode 100644 index 000000000..c89e38dce --- /dev/null +++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolator.java @@ -0,0 +1,48 @@ +/* + * Licensed to the Hipparchus project under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The Hipparchus project 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 + * + * https://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.hipparchus.analysis.interpolation; + +import org.hipparchus.CalculusFieldElement; +import org.hipparchus.exception.MathIllegalArgumentException; + +/** + * Interpolate grid data using bi-linear interpolation. + * @param type of the field elements + * @since 4.0 + */ +public class FieldBilinearInterpolator> + implements FieldBivariateGridInterpolator { + + /** Empty constructor. + *

+ * This constructor is not strictly necessary, but it prevents spurious + * javadoc warnings with JDK 18 and later. + *

+ * @since 3.0 + */ + public FieldBilinearInterpolator() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy + // nothing to do + } + + /** {@inheritDoc} */ + @Override + public FieldBilinearInterpolatingFunction interpolate(final T[] xval, final T[] yval, final T[][] fval) + throws MathIllegalArgumentException { + return new FieldBilinearInterpolatingFunction<>(xval, yval, fval); + } + +} diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBivariateGridInterpolator.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBivariateGridInterpolator.java new file mode 100644 index 000000000..33873909d --- /dev/null +++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldBivariateGridInterpolator.java @@ -0,0 +1,48 @@ +/* + * Licensed to the Hipparchus project under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The Hipparchus project 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 + * + * https://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.hipparchus.analysis.interpolation; + +import org.hipparchus.CalculusFieldElement; +import org.hipparchus.analysis.CalculusFieldBivariateFunction; +import org.hipparchus.exception.MathIllegalArgumentException; + +/** + * Interface representing a bivariate field interpolating function where the + * sample points must be specified on a regular grid. + * @param type of the field elements + * @since 4.0 + */ +public interface FieldBivariateGridInterpolator> { + /** + * Compute an interpolating function for the dataset. + * + * @param xval All the x-coordinates of the interpolation points, sorted + * in increasing order. + * @param yval All the y-coordinates of the interpolation points, sorted + * in increasing order. + * @param fval The values of the interpolation points on all the grid knots: + * {@code fval[i][j] = f(xval[i], yval[j])}. + * @return a function which interpolates the dataset. + * @throws MathIllegalArgumentException if any of the arrays has zero length. + * @throws MathIllegalArgumentException if the array lengths are inconsistent. + * @throws MathIllegalArgumentException if the array is not sorted. + * @throws MathIllegalArgumentException if the number of points is too small for + * the order of the interpolation + */ + CalculusFieldBivariateFunction interpolate(T[] xval, T[] yval, T[][] fval) + throws MathIllegalArgumentException; +} diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldGridAxis.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldGridAxis.java new file mode 100644 index 000000000..be624cb22 --- /dev/null +++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/FieldGridAxis.java @@ -0,0 +1,207 @@ +/* + * Licensed to the Hipparchus project under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The Hipparchus project 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 + * + * https://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.hipparchus.analysis.interpolation; + +import org.hipparchus.CalculusFieldElement; +import org.hipparchus.exception.LocalizedCoreFormats; +import org.hipparchus.exception.MathIllegalArgumentException; +import org.hipparchus.util.FastMath; +import org.hipparchus.util.MathArrays; + +import java.util.concurrent.atomic.AtomicInteger; + +/** + * Helper for finding interpolation nodes along one axis of grid data. + *

+ * This class is intended to be used for interpolating inside grids. + * It works on any sorted data without duplication and size at least + * {@code n} where {@code n} is the number of points required for + * interpolation (i.e. 2 for linear interpolation, 3 for quadratic...) + *

+ *

+ * The method uses linear interpolation to select the nodes indices. + * It should be O(1) for sufficiently regular data, therefore much faster + * than bisection. It also features caching, which improves speed when + * interpolating several points in row in the close locations, i.e. when + * successive calls have a high probability to return the same interpolation + * nodes. This occurs for example when scanning with small steps a loose + * grid. The method also works on non-regular grids, but may be slower in + * this case. + *

+ *

+ * This class is thread-safe. + *

+ * @param Type of the field elements. + * @since 4.0 + */ +public class FieldGridAxis> { + + /** All the coordinates of the interpolation points, sorted in increasing order. */ + private final T[] grid; + + /** Number of points required for interpolation. */ + private final int n; + + /** Cached value of last x index. */ + private final AtomicInteger cache; + + /** Simple constructor. + * @param grid coordinates of the interpolation points, sorted in increasing order + * @param n number of points required for interpolation, i.e. 2 for linear, 3 + * for quadratic... + * @exception MathIllegalArgumentException if grid size is smaller than {@code n} + * or if the grid is not sorted in strict increasing order + */ + public FieldGridAxis(final T[] grid, final int n) + throws MathIllegalArgumentException { + + // safety checks + if (grid.length < n) { + throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DIMENSION, + grid.length, n); + } + MathArrays.checkOrder(grid); + + this.grid = grid.clone(); + this.n = n; + this.cache = new AtomicInteger(0); + + } + + /** Get the number of points of the grid. + * @return number of points of the grid + */ + public int size() { + return grid.length; + } + + /** Get the number of points required for interpolation. + * @return number of points required for interpolation + */ + public int getN() { + return n; + } + + /** Get the interpolation node at specified index. + * @param index node index + * @return coordinate of the node at specified index + */ + public T node(final int index) { + return grid[index]; + } + + /** Get the index of the first interpolation node for some coordinate along the grid. + *

+ * The index return is the one for the lowest interpolation node suitable for + * {@code t}. This means that if {@code i} is returned the nodes to use for + * interpolation at coordinate {@code t} are at indices {@code i}, {@code i+1}, + * ..., {@code i+n-1}, where {@code n} is the number of points required for + * interpolation passed at construction. + *

+ *

+ * The index is selected in order to have the subset of nodes from {@code i} to + * {@code i+n-1} as balanced as possible around {@code t}: + *

+ *
    + *
  • + * if {@code t} is inside the grid and sufficiently far from the endpoints + *
      + *
    • + * if {@code n} is even, the returned nodes will be perfectly balanced: + * there will be {@code n/2} nodes smaller than {@code t} and {@code n/2} + * nodes larger than {@code t} + *
      • + *
    • + * if {@code n} is odd, the returned nodes will be slightly unbalanced by + * one point: there will be {@code (n+1)/2} nodes smaller than {@code t} + * and {@code (n-1)/2} nodes larger than {@code t} + *
      • + *
    + *
    • + *
  • + * if {@code t} is inside the grid and close to endpoints, the returned nodes + * will be unbalanced: there will be less nodes on the endpoints side and + * more nodes on the interior side + *
    • + *
  • + * if {@code t} is outside of the grid, the returned nodes will completely + * off balance: all nodes will be on the same side with respect to {@code t} + *
    • + *
+ *

+ * It is not an error to call this method with {@code t} outside of the grid, + * it simply implies that the interpolation will become an extrapolation and accuracy + * will decrease as {@code t} goes farther from the grid points. This is intended so + * interpolation does not fail near the end of the grid. + *

+ * @param t coordinate of the point to interpolate + * @return index {@code i} such {@link #node(int) node(i)}, {@link #node(int) node(i+1)}, + * ... {@link #node(int) node(i+n-1)} can be used for interpolating a value at + * coordinate {@code t} + * @since 1.4 + */ + public int interpolationIndex(final T t) { + + final int middleOffset = (n - 1) / 2; + int iInf = middleOffset; + int iSup = grid.length - (n - 1) + middleOffset; + + // first try to simply reuse the cached index, + // for faster return in a common case + final int cached = cache.get(); + final int middle = cached + middleOffset; + final T aMid0 = grid[middle]; + final T aMid1 = grid[middle + 1]; + if (t.getReal() < aMid0.getReal()) { + if (middle == iInf) { + // we are in the unbalanced low area + return cached; + } + } else if (t.getReal() < aMid1.getReal()) { + // we are in the balanced middle area + return cached; + } else { + if (middle == iSup - 1) { + // we are in the unbalanced high area + return cached; + } + } + + // we need to find a new index + T aInf = grid[iInf]; + T aSup = grid[iSup]; + while (iSup - iInf > 1) { + final int iInterp = (int) ((iInf * (aSup.getReal() - t.getReal()) + iSup * (t.getReal() - aInf.getReal())) / (aSup.getReal() - aInf.getReal())); + final int iMed = FastMath.max(iInf + 1, FastMath.min(iInterp, iSup - 1)); + if (t.getReal() < grid[iMed].getReal()) { + // keeps looking in the lower part of the grid + iSup = iMed; + aSup = grid[iSup]; + } else { + // keeps looking in the upper part of the grid + iInf = iMed; + aInf = grid[iInf]; + } + } + + final int newCached = iInf - middleOffset; + cache.compareAndSet(cached, newCached); + return newCached; + + } + +} diff --git a/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/GridAxis.java b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/GridAxis.java index 7eb1e3155..7f5323396 100644 --- a/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/GridAxis.java +++ b/hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/GridAxis.java @@ -36,7 +36,7 @@ * The method uses linear interpolation to select the nodes indices. * It should be O(1) for sufficiently regular data, therefore much faster * than bisection. It also features caching, which improves speed when - * interpolating several points in raw in the close locations, i.e. when + * interpolating several points in row in the close locations, i.e. when * successive calls have a high probability to return the same interpolation * nodes. This occurs for example when scanning with small steps a loose * grid. The method also works on non-regular grids, but may be slower in @@ -56,7 +56,7 @@ public class GridAxis implements Serializable { private final double[] grid; /** Number of points required for interpolation. */ - private int n; + private final int n; /** Cached value of last x index. */ private final AtomicInteger cache; diff --git a/hipparchus-core/src/test/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolatingFunctionTest.java b/hipparchus-core/src/test/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolatingFunctionTest.java new file mode 100644 index 000000000..8cd4ce49f --- /dev/null +++ b/hipparchus-core/src/test/java/org/hipparchus/analysis/interpolation/FieldBilinearInterpolatingFunctionTest.java @@ -0,0 +1,215 @@ +/* + * Licensed to the Hipparchus project under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The Hipparchus project 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 + * + * https://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.hipparchus.analysis.interpolation; + +import org.hipparchus.analysis.CalculusFieldBivariateFunction; +import org.hipparchus.analysis.differentiation.UnivariateDerivative1; +import org.hipparchus.random.RandomVectorGenerator; +import org.hipparchus.random.SobolSequenceGenerator; +import org.hipparchus.util.FastMath; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * Test case for the field bilinear function. + */ +final class FieldBilinearInterpolatingFunctionTest { + + @Test + void testConstant() { + + double xMin = 0.0; + double xMax = 7.0; + int nx = 15; + UnivariateDerivative1[] xVal = createLinearGrid(xMin, xMax, nx, 1.0); + + double yMin = -5.0; + double yMax = 5.0; + int ny = 11; + UnivariateDerivative1[] yVal = createLinearGrid(yMin, yMax, ny, 2.0); + + CalculusFieldBivariateFunction f = (x, y) -> new UnivariateDerivative1(3.5, 1.0); + FieldBilinearInterpolatingFunction bif = createInterpolatingFunction(xVal, yVal, f); + + assertEquals(xMin, bif.getXInf().getValue(), 1.0e-15); + assertEquals(1.0, bif.getXInf().getFirstDerivative(), 1.0e-15); + assertEquals(xMax, bif.getXSup().getValue(), 1.0e-15); + assertEquals(1.0, bif.getXSup().getFirstDerivative(), 1.0e-15); + assertEquals(yMin, bif.getYInf().getValue(), 1.0e-15); + assertEquals(2.0, bif.getYInf().getFirstDerivative(), 1.0e-15); + assertEquals(yMax, bif.getYSup().getValue(), 1.0e-15); + assertEquals(2.0, bif.getYSup().getFirstDerivative(), 1.0e-15); + + checkInterpolationAtNodes(xVal, yVal, bif, f, 1.0e-15, 1.0e-15); + checkInterpolationRandom(new SobolSequenceGenerator(2), xMin, xMax, yMin, yMax, bif, f, 1.0e-15, 1.0e-15); + + } + + @Test + void testLinear() { + + double xMin = -5.0; + double xMax = 5.0; + int nx = 11; + UnivariateDerivative1[] xVal = createLinearGrid(xMin, xMax, nx, 1.0); + + double yMin = 0.0; + double yMax = 7.0; + int ny = 15; + UnivariateDerivative1[] yVal = createLinearGrid(yMin, yMax, ny, 2.0); + + CalculusFieldBivariateFunction f = (x, y) -> x.multiply(2).subtract(y); + FieldBilinearInterpolatingFunction bif = createInterpolatingFunction(xVal, yVal, f); + + assertEquals(xMin, bif.getXInf().getValue(), 1.0e-15); + assertEquals(1.0, bif.getXInf().getFirstDerivative(), 1.0e-15); + assertEquals(xMax, bif.getXSup().getValue(), 1.0e-15); + assertEquals(1.0, bif.getXSup().getFirstDerivative(), 1.0e-15); + assertEquals(yMin, bif.getYInf().getValue(), 1.0e-15); + assertEquals(2.0, bif.getYInf().getFirstDerivative(), 1.0e-15); + assertEquals(yMax, bif.getYSup().getValue(), 1.0e-15); + assertEquals(2.0, bif.getYSup().getFirstDerivative(), 1.0e-15); + + checkInterpolationAtNodes(xVal, yVal, bif, f, 1.0e-15, 1.0e-15); + checkInterpolationRandom(new SobolSequenceGenerator(2), xMin, xMax, yMin, yMax, bif, f, 1.0e-15, 1.0e-15); + + } + + @Test + void testQuadratic() { + + double xMin = -5.0; + double xMax = 5.0; + int nx = 11; + UnivariateDerivative1[] xVal = createLinearGrid(xMin, xMax, nx, 1.0); + + double yMin = 0.0; + double yMax = 7.0; + int ny = 15; + UnivariateDerivative1[] yVal = createLinearGrid(yMin, yMax, ny, 2.0); + + CalculusFieldBivariateFunction f = (x, y) -> x.multiply(3).subtract(2).multiply(y.multiply(-0.5).add(6)); + FieldBilinearInterpolatingFunction bif = createInterpolatingFunction(xVal, yVal, f); + + assertEquals(xMin, bif.getXInf().getValue(), 1.0e-15); + assertEquals(1.0, bif.getXInf().getFirstDerivative(), 1.0e-15); + assertEquals(xMax, bif.getXSup().getValue(), 1.0e-15); + assertEquals(1.0, bif.getXSup().getFirstDerivative(), 1.0e-15); + assertEquals(yMin, bif.getYInf().getValue(), 1.0e-15); + assertEquals(2.0, bif.getYInf().getFirstDerivative(), 1.0e-15); + assertEquals(yMax, bif.getYSup().getValue(), 1.0e-15); + assertEquals(2.0, bif.getYSup().getFirstDerivative(), 1.0e-15); + + checkInterpolationAtNodes(xVal, yVal, bif, f, 1.0e-15, 1.0e-15); + checkInterpolationRandom(new SobolSequenceGenerator(2), xMin, xMax, yMin, yMax, bif, f, 1.0e-15, 1.0e-15); + + } + + @Test + void testSinCos() { + doTestSinCos( 10, 10, 1.8e-2, 8.3e-2); + doTestSinCos( 100, 100, 1.5e-4, 7.6e-3); + doTestSinCos(1000, 1000, 1.4e-6, 7.6e-4); + } + + private void doTestSinCos(final int nx, final int ny, final double tol0, final double tol1) { + double xMin = -1.0; + double xMax = 2.0; + UnivariateDerivative1[] xVal = createLinearGrid(xMin, xMax, nx, 1.0); + + double yMin = 0.0; + double yMax = 1.5; + UnivariateDerivative1[] yVal = createLinearGrid(yMin, yMax, ny, 2.0); + + CalculusFieldBivariateFunction f = (x, y) -> FastMath.sin(x).multiply(FastMath.cos(y)); + FieldBilinearInterpolatingFunction bif = createInterpolatingFunction(xVal, yVal, f); + + assertEquals(xMin, bif.getXInf().getValue(), 1.0e-15); + assertEquals(1.0, bif.getXInf().getFirstDerivative(), 1.0e-15); + assertEquals(xMax, bif.getXSup().getValue(), 1.0e-15); + assertEquals(1.0, bif.getXSup().getFirstDerivative(), 1.0e-15); + assertEquals(yMin, bif.getYInf().getValue(), 1.0e-15); + assertEquals(2.0, bif.getYInf().getFirstDerivative(), 1.0e-15); + assertEquals(yMax, bif.getYSup().getValue(), 1.0e-15); + assertEquals(2.0, bif.getYSup().getFirstDerivative(), 1.0e-15); + + checkInterpolationAtNodes(xVal, yVal, bif, f, 1.0e-15, 1.0e-15); + checkInterpolationRandom(new SobolSequenceGenerator(2), xMin, xMax, yMin, yMax, bif, f, tol0, tol1); + + } + + private UnivariateDerivative1[] createLinearGrid(final double min, final double max, final int n, + final double derivative) { + final UnivariateDerivative1[] grid = new UnivariateDerivative1[n]; + for (int i = 0; i < n; ++i) { + grid[i] = new UnivariateDerivative1(((n - 1 - i) * min + i * max) / (n - 1), derivative); + } + return grid; + } + + private FieldBilinearInterpolatingFunction + createInterpolatingFunction(UnivariateDerivative1[] xVal, UnivariateDerivative1[] yVal, + CalculusFieldBivariateFunction f) { + final UnivariateDerivative1[][] fVal = new UnivariateDerivative1[xVal.length][yVal.length]; + for (int i = 0; i < xVal.length; ++i) { + for (int j = 0; j < yVal.length; ++j) { + fVal[i][j] = f.value(xVal[i], yVal[j]); + } + } + return new FieldBilinearInterpolator().interpolate(xVal, yVal, fVal); + } + + private void checkInterpolationAtNodes(final UnivariateDerivative1[] xVal, + final UnivariateDerivative1[] yVal, + final FieldBilinearInterpolatingFunction bif, + final CalculusFieldBivariateFunction f, + final double tol0, final double tol1) { + for (final UnivariateDerivative1 x : xVal) { + for (final UnivariateDerivative1 y : yVal) { + assertEquals(f.value(x, y).getValue(), bif.value(x, y).getValue(), tol0); + assertEquals(f.value(x, y).getFirstDerivative(), bif.value(x, y).getFirstDerivative(), tol1); + } + } + } + + private void checkInterpolationRandom(final RandomVectorGenerator random, + final double xMin, final double xMax, + final double yMin, final double yMax, + final FieldBilinearInterpolatingFunction bif, + final CalculusFieldBivariateFunction f, + final double tol0, final double tol1) { + double maxError0 = 0.0; + double maxError1 = 0.0; + for (int i = 0; i < 10000; ++i) { + + final double[] v = random.nextVector(); + + final UnivariateDerivative1 x = new UnivariateDerivative1(xMin + v[0] * (xMax - xMin), 1.0); + final UnivariateDerivative1 y = new UnivariateDerivative1(yMin + v[1] * (yMax - yMin), 1.0); + final UnivariateDerivative1 delta = f.value(x, y).subtract(bif.value(x, y)); + maxError0 = FastMath.max(maxError0, FastMath.abs(delta.getValue())); + maxError1 = FastMath.max(maxError1, FastMath.abs(delta.getFirstDerivative())); + + } + + assertEquals(0.0, maxError0, tol0); + assertEquals(0.0, maxError1, tol1); + + } + +}