/*

* Copyright (c) NM LTD.

* https://nm.dev/

*

* THIS SOFTWARE IS LICENSED, NOT SOLD.

*

* YOU MAY USE THIS SOFTWARE ONLY AS DESCRIBED IN THE LICENSE.

* IF YOU ARE NOT AWARE OF AND/OR DO NOT AGREE TO THE TERMS OF THE LICENSE,

* DO NOT USE THIS SOFTWARE.

*

* THE SOFTWARE IS PROVIDED "AS IS", WITH NO WARRANTY WHATSOEVER,

* EITHER EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION,

* ANY WARRANTIES OF ACCURACY, ACCESSIBILITY, COMPLETENESS,

* FITNESS FOR A PARTICULAR PURPOSE, MERCHANTABILITY, NON-INFRINGEMENT,

* TITLE AND USEFULNESS.

*

* IN NO EVENT AND UNDER NO LEGAL THEORY,

* WHETHER IN ACTION, CONTRACT, NEGLIGENCE, TORT, OR OTHERWISE,

* SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR

* ANY CLAIMS, DAMAGES OR OTHER LIABILITIES,

* ARISING AS A RESULT OF USING OR OTHER DEALINGS IN THE SOFTWARE.

*/

package dev.nm.nmj;

import java.io.BufferedWriter;

import java.io.FileWriter;

import java.io.IOException;

import dev.nm.algebra.linear.vector.doubles.dense.DenseVector;

import dev.nm.analysis.curvefit.LeastSquares;

import dev.nm.analysis.curvefit.interpolation.bivariate.BicubicInterpolation;

import dev.nm.analysis.curvefit.interpolation.bivariate.BicubicSpline;

import dev.nm.analysis.curvefit.interpolation.bivariate.BilinearInterpolation;

import dev.nm.analysis.curvefit.interpolation.bivariate.BivariateArrayGrid;

import dev.nm.analysis.curvefit.interpolation.bivariate.BivariateGrid;

import dev.nm.analysis.curvefit.interpolation.bivariate.BivariateGridInterpolation;

import dev.nm.analysis.curvefit.interpolation.bivariate.BivariateRegularGrid;

import dev.nm.analysis.curvefit.interpolation.multivariate.MultivariateArrayGrid;

import dev.nm.analysis.curvefit.interpolation.multivariate.RecursiveGridInterpolation;

import dev.nm.analysis.curvefit.interpolation.univariate.CubicHermite;

import dev.nm.analysis.curvefit.interpolation.univariate.CubicSpline;

import dev.nm.analysis.curvefit.interpolation.univariate.Interpolation;

import dev.nm.analysis.curvefit.interpolation.univariate.LinearInterpolation;

import dev.nm.analysis.curvefit.interpolation.univariate.NewtonPolynomial;

import dev.nm.analysis.function.rn2r1.RealScalarFunction;

import dev.nm.analysis.function.rn2r1.univariate.UnivariateRealFunction;

import dev.nm.analysis.function.tuple.OrderedPairs;

import dev.nm.analysis.function.tuple.SortedOrderedPairs;

import dev.nm.misc.datastructure.MultiDimensionalArray;

import java.io.File;

import static java.lang.Math.log;

/**

* Numerical Methods Using Java: For Data Science, Analysis, and Engineering

*

* @author haksunli

* @see

* https://www.amazon.com/Numerical-Methods-Using-Java-Engineering/dp/1484267966

* https://nm.dev/

*/

public class Chapter5 {

public static void main(String[] args) throws IOException {

System.out.println("Chapter 5 demos");

Chapter5 chapter5 = new Chapter5();

chapter5.least_square_curve_fitting();

chapter5.linear_interpolation();

chapter5.cubic_Hermite_interpolation();

chapter5.cubic_spline_interpolation();

chapter5.newton_polynomial_interpolation();

chapter5.bivariate_interpolation();

chapter5.bivariate_interpolation_using_derivatives();

chapter5.multivariate_interpolation();

}

public void least_square_curve_fitting() {

System.out.println("least square curve fitting");

// the data set

OrderedPairs data = new SortedOrderedPairs(

new double[]{0., 1., 2., 3., 4., 5.},

new double[]{0., 1., 1.414, 1.732, 2., 2.236}

);

LeastSquares ls = new LeastSquares(2);

UnivariateRealFunction f = ls.fit(data);

System.out.println(String.format("f(%.0f)=%f", 0., f.evaluate(0.))); // f(0) = 0.09

System.out.println(String.format("f(%.0f)=%f", 1., f.evaluate(1.))); // f(1) = 0.82

System.out.println(String.format("f(%.0f)=%f", 2., f.evaluate(2.))); // f(2) = 1.39

System.out.println(String.format("f(%.0f)=%f", 3., f.evaluate(3.))); // f(3) = 1.81

System.out.println(String.format("f(%.0f)=%f", 4., f.evaluate(4.))); // f(4) = 2.07

System.out.println(String.format("f(%.0f)=%f", 5., f.evaluate(5.))); // f(5) = 2.17

}

public void linear_interpolation() throws IOException {

System.out.println("linear interpolation");

// the data set

OrderedPairs data = new SortedOrderedPairs(

new double[]{0., 0.7, 1.4, 2.1, 2.8, 3.5, 4.2, 4.9, 5.6, 6.3},

new double[]{0., 0.644218, 0.98545, 0.863209, 0.334988, -0.350783, -0.871576, -0.982453, -0.631267, 0.0168139}

);

LinearInterpolation li = new LinearInterpolation();

UnivariateRealFunction f = li.fit(data);

System.out.println(f.evaluate(2)); // f(2) = 0.880672

System.out.println(f.evaluate(3)); // f(3) = 0.139053

plot(f, 100, 0, 6.5, "./plots/chapter5/figure_5_6/f_sample.txt");

}

public void cubic_Hermite_interpolation() throws IOException {

System.out.println("cubic_Hermite interpolation");

// the data set

OrderedPairs data = new SortedOrderedPairs(

new double[]{0., 0.7, 1.4, 2.1, 2.8, 3.5, 4.2, 4.9, 5.6, 6.3},

new double[]{0., 0.644218, 0.98545, 0.863209, 0.334988, -0.350783, -0.871576, -0.982453, -0.631267, 0.0168139}

);

CubicHermite spline = new CubicHermite(CubicHermite.Tangents.CATMULL_ROM);

// CubicHermite spline = new CubicHermite(CubicHermite.Tangents.FINITE_DIFFERENCE);

UnivariateRealFunction f = spline.fit(data);

System.out.println(f.evaluate(2)); // f(2) = 0.906030

System.out.println(f.evaluate(3)); // f(3) = 0.145727

plot(f, 100, 0, 6.3, "./plots/chapter5/figure_5_7/f_sample.txt");

plot(f, 100, 0.7, 2.1, "./plots/chapter5/figure_5_8/f_sample.txt");

}

public void cubic_spline_interpolation() throws IOException {

System.out.println("cubic spline interpolation");

// the data set

OrderedPairs data = new SortedOrderedPairs(

new double[]{0., 1., 2., 3., 4., 5.},

new double[]{0., 3.5, 5., 3., 1., 4.}

);

CubicSpline cs1 = CubicSpline.natural();

UnivariateRealFunction f1 = cs1.fit(data);

plot(f1, 100, 0, 5, "./plots/chapter5/figure_5_9/natural_cspline_sample.txt");

CubicSpline cs2 = CubicSpline.clamped();

UnivariateRealFunction f2 = cs2.fit(data);

plot(f2, 100, 0, 5, "./plots/chapter5/figure_5_9/clamped_cspline_sample.txt");

CubicSpline cs3 = CubicSpline.notAKnot();

UnivariateRealFunction f3 = cs3.fit(data);

plot(f3, 100, 0, 5, "./plots/chapter5/figure_5_9/notaknot_cspline_sample.txt");

}

public void newton_polynomial_interpolation() throws IOException {

System.out.println("Newton polynomial interpolation");

// 2 data points, linear form

OrderedPairs data1 = new SortedOrderedPairs(

new double[]{1., 3.},

new double[]{log(1.), log(3.)}

);

Interpolation np1 = new NewtonPolynomial();

UnivariateRealFunction f1 = np1.fit(data1);

plot(f1, 100, 1, 3, "./plots/chapter5/figure_5_10/linear_newton_sample.txt");

// 3 data points, quadratic form

OrderedPairs data2 = new SortedOrderedPairs(

new double[]{1., 2., 3.},

new double[]{log(1.), log(2.), log(3.)}

);

Interpolation np2 = new NewtonPolynomial();

UnivariateRealFunction f2 = np2.fit(data2);

plot(f2, 100, 1, 3, "./plots/chapter5/figure_5_11/quadratic_newton_sample.txt");

// comparison between Newton polynomial and cubic spline

OrderedPairs data3 = new SortedOrderedPairs(

new double[]{1., 2., 3., 4., 5., 6., 7.},

new double[]{3., 4., 2., 5., 4., 3., 6.}

);

Interpolation np3 = new NewtonPolynomial();

UnivariateRealFunction f3_1 = np3.fit(data3);

plot(f3_1, 500, 1, 7, "./plots/chapter5/figure_5_12/newton_sample.txt");

Interpolation cs = CubicSpline.natural();

UnivariateRealFunction f3_2 = cs.fit(data3);

plot(f3_2, 500, 1, 7, "./plots/chapter5/figure_5_12/cspline_sample.txt");

}

public void bivariate_interpolation() throws IOException {

System.out.println("bivariate interpolation");

BivariateGrid grids = new BivariateArrayGrid(

new double[][]{

{1, 1, 1}, // z(1, 1) = 1, z(1, 2) = 1, z(1, 3) = 1

{2, 4, 8}, // z(2, 1) = 2, z(2, 2) = 4, z(2, 3) = 8

{3, 9, 27} // z(3, 1) = 3, z(3, 2) = 9, z(3, 3) = 27

},

new double[]{1, 2, 3}, // x

new double[]{1, 2, 3} // y

);

BivariateGridInterpolation bl = new BilinearInterpolation();

RealScalarFunction f1 = bl.interpolate(grids); // f3(1.5, 1.5) = 2.0

System.out.println(f1.evaluate(new DenseVector(new double[]{1.5, 1.5})));

BivariateGridInterpolation bs = new BicubicSpline();

RealScalarFunction f2 = bs.interpolate(grids); // f2(1.5, 1.5) = 1.8828125

System.out.println(f2.evaluate(new DenseVector(new double[]{1.5, 1.5})));

BivariateGridInterpolation bi = new BicubicInterpolation();

RealScalarFunction f3 = bi.interpolate(grids); // f1(1.5, 1.5) = 1.90625

System.out.println(f3.evaluate(new DenseVector(new double[]{1.5, 1.5})));

plot(f3, 30, new double[]{1, 3}, new double[]{1, 3}, "./plots/chapter5/figure_5_13/bicspline_sample.txt");

}

public void bivariate_interpolation_using_derivatives() throws IOException {

System.out.println("bivariate interpolation using derivatives");

// derivatives and answers from Michael Flanagan's library

double[][] z = new double[][]{

{1.0, 3.0, 5.0},

{2.0, 4.0, 8.0},

{9.0, 10.0, 11.0}

};

final double[][] dx = new double[][]{

{6.0, 2.0, 2.0},

{6.0, 7.0, 8.0},

{6.0, 12.0, 14.0}

};

final double[][] dy = new double[][]{

{8.0, 8.0, 8.0},

{16.0, 12.0, 8.0},

{4.0, 4.0, 4.0}

};

final double[][] dxdy = new double[][]{

{16.0, 8.0, 0.0},

{-4.0, -4.0, -4.0},

{-24.0, -16.0, -8.0}

};

BicubicInterpolation.PartialDerivatives deriv

= new BicubicInterpolation.PartialDerivatives() {

@Override

public double dx(BivariateGrid grid, int i, int j) {

return getDeriv(dx, i, j); // for some reason the y-axis is written in reverse...

}

@Override

public double dy(BivariateGrid grid, int i, int j) {

return getDeriv(dy, i, j);

}

@Override

public double dxdy(BivariateGrid grid, int i, int j) {

return getDeriv(dxdy, i, j);

}

private double getDeriv(double[][] dx, int i, int j) {

return dx[i][2 - j];

}

};

BivariateGridInterpolation interpolation = new BicubicInterpolation(deriv);

BivariateGrid grid = new BivariateRegularGrid(z, 0.0, 0.0, 0.5, 0.25);

RealScalarFunction f = interpolation.interpolate(grid);

System.out.println(f.evaluate(new DenseVector(0.0, 0.0))); // 1.0

System.out.println(f.evaluate(new DenseVector(0.0, 0.125))); // 2.0

System.out.println(f.evaluate(new DenseVector(0.0, 0.25))); // 3.0

System.out.println(f.evaluate(new DenseVector(0.0, 0.375))); // 4.0

System.out.println(f.evaluate(new DenseVector(0.0, 0.5))); // 5.0

System.out.println(f.evaluate(new DenseVector(0.25, 0.0))); // 1.125

System.out.println(f.evaluate(new DenseVector(0.25, 0.125))); // 2.078125

System.out.println(f.evaluate(new DenseVector(0.25, 0.25))); // 3.1875

System.out.println(f.evaluate(new DenseVector(0.25, 0.375))); // 4.765625

System.out.println(f.evaluate(new DenseVector(0.25, 0.5))); // 6.5

System.out.println(f.evaluate(new DenseVector(0.5, 0.0))); // 2.0

System.out.println(f.evaluate(new DenseVector(0.5, 0.125))); // 2.875

System.out.println(f.evaluate(new DenseVector(0.5, 0.25))); // 4.0

System.out.println(f.evaluate(new DenseVector(0.5, 0.375))); // 5.875

System.out.println(f.evaluate(new DenseVector(0.5, 0.5))); // 8.0

System.out.println(f.evaluate(new DenseVector(0.75, 0.0))); // 5.125

System.out.println(f.evaluate(new DenseVector(0.75, 0.125))); // 5.828125

System.out.println(f.evaluate(new DenseVector(0.75, 0.25))); // 6.6875

System.out.println(f.evaluate(new DenseVector(0.75, 0.375))); // 8.015625

System.out.println(f.evaluate(new DenseVector(0.75, 0.5))); // 9.5

System.out.println(f.evaluate(new DenseVector(1.0, 0.0))); // 9.0

System.out.println(f.evaluate(new DenseVector(1.0, 0.125))); // 9.5

System.out.println(f.evaluate(new DenseVector(1.0, 0.25))); // 10.0

System.out.println(f.evaluate(new DenseVector(1.0, 0.375))); // 10.5

System.out.println(f.evaluate(new DenseVector(1.0, 0.5))); // 11.0

}

public void multivariate_interpolation() {

// the data set

MultiDimensionalArray<Double> mda

= new MultiDimensionalArray<>(2, 2, 2);

mda.set(1., 0, 0, 0); // mda[0][0][0] = 1.

mda.set(2., 1, 0, 0);

mda.set(3., 0, 1, 0);

mda.set(4., 0, 0, 1);

mda.set(5., 1, 1, 0);

mda.set(6., 1, 0, 1);

mda.set(7., 0, 1, 1);

mda.set(8., 1, 1, 1);

MultivariateArrayGrid mvGrid = new MultivariateArrayGrid(

mda,

new double[]{1, 2},

new double[]{1, 2},

new double[]{1, 2}

);

RecursiveGridInterpolation rgi

= new RecursiveGridInterpolation(new LinearInterpolation());

RealScalarFunction f = rgi.interpolate(mvGrid);

System.out.println(f.evaluate(new DenseVector(new double[]{1.5, 1.5, 1.5}))); // f(1.5, 1.5, 1.5) = 4.5

}

/**

* Export sample points to plot a function curve.

*

* @param f the function to plot

* @param nSamples number of sample points

* @param rangeStart start of the plotting range

* @param rangeEnd end of the plotting range

* @param filename export file name

* @throws IOException

*/

private void plot(

UnivariateRealFunction f,

int nSamples,

double rangeStart,

double rangeEnd,

String filename

) throws IOException {

File file = new File(filename);

if (!file.getParentFile().exists()) {

file.getParentFile().mkdirs();

}

file.createNewFile();

try (BufferedWriter writer = new BufferedWriter(new FileWriter(file))) {

double gridSize = (rangeEnd - rangeStart) / (nSamples - 1);

double x = rangeStart;

for (int i = 0; i < nSamples; i++) {

writer.write(String.format("%f %f\n", x, f.evaluate(x)));

// writer.write(new StringBuilder().append(x).append(" ").append(f.evaluate(x)).toString());

// writer.newLine();

x += gridSize;

}

}

}

/**

* Export sample points to plot a bivariate function surface.

*

* @param f the function to plot

* @param nSamples number of sample points

* @param rangeX plotting range of x

* @param rangeY plotting range of y

* @param filename export file name

* @throws IOException

*/

private void plot(

RealScalarFunction f,

int nSamples,

double[] rangeX,

double[] rangeY,

String filename

) throws IOException {

File file = new File(filename);

if (!file.getParentFile().exists()) {

file.getParentFile().mkdirs();

}

file.createNewFile();

try (BufferedWriter writer = new BufferedWriter(new FileWriter(file))) {

double gridSizeX = (rangeX[1] - rangeX[0]) / (nSamples - 1);

double gridSizeY = (rangeY[1] - rangeY[0]) / (nSamples - 1);

double x = rangeX[0];

for (int i = 0; i < nSamples; i++) {

double y = rangeY[0];

for (int j = 0; j < nSamples; j++) {

writer.write(String.format("%f %f %s\n", x, y, f.evaluate(new DenseVector(x, y)).toString()));

// writer.write(new StringBuilder().append(x).append(" ").append(y).append(" ").append(f.evaluate(new DenseVector(x, y))).toString());

// writer.newLine();

y += gridSizeY;

}

x += gridSizeX;

}

}

}

}