@Test public void testPolyFit() throws Exception { final int maxDegree = 6; double[] x = MathUtils.increment(11, 0, 0.1); double[] y = new double[x.length]; for (int i = 0; i < x.length; i++) { y[i] = erf(x[i]); } ArrayList<double[]> coeff_EXPECTED = new ArrayList<double[]>(); coeff_EXPECTED.add(new double[]{0.0531277556004621, 0.853026625583286}); coeff_EXPECTED.add(new double[]{-0.00634545069301346, 1.24951466753979, -0.396488041956504}); coeff_EXPECTED.add(new double[]{-0.00117120913502653, 1.16730171834066, -0.180894643707045, -0.143728932166306}); coeff_EXPECTED.add(new double[]{1.87578727256749e-05, 1.12598341946038, 0.0256968506943939, -0.474275323208611, 0.165273195521153}); coeff_EXPECTED.add(new double[]{7.50968728582951e-06, 1.12686515221902, 0.0184323975978185, -0.453809874699964, 0.141839475854753, 0.00937348786656298}); coeff_EXPECTED.add(new double[]{6.40296881978214e-08, 1.12837041599673, -3.2419597312071e-05, -0.374389526989057, -0.0116238001908376, 0.145877210494725, -0.0455012408760536}); for (int degree = 0; degree < maxDegree; degree++) { double[] coeff_ACTUAL = PolyUtils.polyFit(new DoubleMatrix(x), new DoubleMatrix(y), degree + 1); Assert.assertArrayEquals(coeff_EXPECTED.get(degree), coeff_ACTUAL, DELTA_03); } }
public static double[] polyFit(DoubleMatrix t, DoubleMatrix y, final int degree) throws IllegalArgumentException { setLoggerLevel(); if (t.length != y.length || !t.isVector() || !y.isVector()) { logger.error("polyfit: require same size vectors."); throw new IllegalArgumentException("polyfit: require same size vectors."); } final int numOfPoints = t.length; final int numOfUnknowns = degree + 1; logger.debug("Degree of interpolating polynomial: {}", degree); logger.debug("Number of unknowns: {}", numOfUnknowns); logger.debug("Number of data points: {}", numOfPoints); if (numOfPoints < numOfUnknowns) { logger.error("Number of points is smaller than parameters solved for."); throw new IllegalArgumentException("Number of points is smaller than parameters solved for."); } logger.debug("Setting up linear system of equations"); DoubleMatrix A = new DoubleMatrix(numOfPoints, numOfUnknowns); for (int j = 0; j <= degree; j++) { A.putColumn(j, pow(t, j)); } logger.debug("Solving lin. system of equations with Cholesky."); DoubleMatrix N = A.transpose().mmul(A); DoubleMatrix rhs = A.transpose().mmul(y); DoubleMatrix x = Solve.solveSymmetric(N, rhs); DoubleMatrix Qx_hat = Solve.solveSymmetric(N, DoubleMatrix.eye(N.getRows())); double maxDeviation = (N.mmul(Qx_hat).sub(DoubleMatrix.eye(Qx_hat.rows))).normmax(); logger.debug("polyfit orbit: max(abs(N*inv(N)-I)) = " + maxDeviation); if (maxDeviation > 1e-6) { logger.warn("polyfit orbit: max(abs(N*inv(N)-I)) = {}", maxDeviation); logger.warn("polyfit orbit interpolation unstable!"); } DoubleMatrix y_hat = A.mmul(x); DoubleMatrix e_hat = y.sub(y_hat); if (e_hat.normmax() > 0.02) { logger.warn("WARNING: Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } else { logger.debug("Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } if (logger.isDebugEnabled()) { logger.debug("REPORTING POLYFIT LEAST SQUARES ERRORS"); logger.debug(" time \t\t\t y \t\t\t yhat \t\t\t ehat"); for (int i = 0; i < numOfPoints; i++) { logger.debug(" " + t.get(i) + "\t" + y.get(i) + "\t" + y_hat.get(i) + "\t" + e_hat.get(i)); } for (int i = 0; i < numOfPoints - 1; i++) { double dt = t.get(i + 1) - t.get(i); logger.debug("Time step between point " + i + 1 + " and " + i + "= " + dt); if (Math.abs(dt - (t.get(1) - t.get(0))) > 0.001) logger.warn("WARNING: Orbit: data does not have equidistant time interval?"); } } return x.toArray(); }
PolyUtils { public static double[] polyFit(DoubleMatrix t, DoubleMatrix y, final int degree) throws IllegalArgumentException { setLoggerLevel(); if (t.length != y.length || !t.isVector() || !y.isVector()) { logger.error("polyfit: require same size vectors."); throw new IllegalArgumentException("polyfit: require same size vectors."); } final int numOfPoints = t.length; final int numOfUnknowns = degree + 1; logger.debug("Degree of interpolating polynomial: {}", degree); logger.debug("Number of unknowns: {}", numOfUnknowns); logger.debug("Number of data points: {}", numOfPoints); if (numOfPoints < numOfUnknowns) { logger.error("Number of points is smaller than parameters solved for."); throw new IllegalArgumentException("Number of points is smaller than parameters solved for."); } logger.debug("Setting up linear system of equations"); DoubleMatrix A = new DoubleMatrix(numOfPoints, numOfUnknowns); for (int j = 0; j <= degree; j++) { A.putColumn(j, pow(t, j)); } logger.debug("Solving lin. system of equations with Cholesky."); DoubleMatrix N = A.transpose().mmul(A); DoubleMatrix rhs = A.transpose().mmul(y); DoubleMatrix x = Solve.solveSymmetric(N, rhs); DoubleMatrix Qx_hat = Solve.solveSymmetric(N, DoubleMatrix.eye(N.getRows())); double maxDeviation = (N.mmul(Qx_hat).sub(DoubleMatrix.eye(Qx_hat.rows))).normmax(); logger.debug("polyfit orbit: max(abs(N*inv(N)-I)) = " + maxDeviation); if (maxDeviation > 1e-6) { logger.warn("polyfit orbit: max(abs(N*inv(N)-I)) = {}", maxDeviation); logger.warn("polyfit orbit interpolation unstable!"); } DoubleMatrix y_hat = A.mmul(x); DoubleMatrix e_hat = y.sub(y_hat); if (e_hat.normmax() > 0.02) { logger.warn("WARNING: Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } else { logger.debug("Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } if (logger.isDebugEnabled()) { logger.debug("REPORTING POLYFIT LEAST SQUARES ERRORS"); logger.debug(" time \t\t\t y \t\t\t yhat \t\t\t ehat"); for (int i = 0; i < numOfPoints; i++) { logger.debug(" " + t.get(i) + "\t" + y.get(i) + "\t" + y_hat.get(i) + "\t" + e_hat.get(i)); } for (int i = 0; i < numOfPoints - 1; i++) { double dt = t.get(i + 1) - t.get(i); logger.debug("Time step between point " + i + 1 + " and " + i + "= " + dt); if (Math.abs(dt - (t.get(1) - t.get(0))) > 0.001) logger.warn("WARNING: Orbit: data does not have equidistant time interval?"); } } return x.toArray(); } }
PolyUtils { public static double[] polyFit(DoubleMatrix t, DoubleMatrix y, final int degree) throws IllegalArgumentException { setLoggerLevel(); if (t.length != y.length || !t.isVector() || !y.isVector()) { logger.error("polyfit: require same size vectors."); throw new IllegalArgumentException("polyfit: require same size vectors."); } final int numOfPoints = t.length; final int numOfUnknowns = degree + 1; logger.debug("Degree of interpolating polynomial: {}", degree); logger.debug("Number of unknowns: {}", numOfUnknowns); logger.debug("Number of data points: {}", numOfPoints); if (numOfPoints < numOfUnknowns) { logger.error("Number of points is smaller than parameters solved for."); throw new IllegalArgumentException("Number of points is smaller than parameters solved for."); } logger.debug("Setting up linear system of equations"); DoubleMatrix A = new DoubleMatrix(numOfPoints, numOfUnknowns); for (int j = 0; j <= degree; j++) { A.putColumn(j, pow(t, j)); } logger.debug("Solving lin. system of equations with Cholesky."); DoubleMatrix N = A.transpose().mmul(A); DoubleMatrix rhs = A.transpose().mmul(y); DoubleMatrix x = Solve.solveSymmetric(N, rhs); DoubleMatrix Qx_hat = Solve.solveSymmetric(N, DoubleMatrix.eye(N.getRows())); double maxDeviation = (N.mmul(Qx_hat).sub(DoubleMatrix.eye(Qx_hat.rows))).normmax(); logger.debug("polyfit orbit: max(abs(N*inv(N)-I)) = " + maxDeviation); if (maxDeviation > 1e-6) { logger.warn("polyfit orbit: max(abs(N*inv(N)-I)) = {}", maxDeviation); logger.warn("polyfit orbit interpolation unstable!"); } DoubleMatrix y_hat = A.mmul(x); DoubleMatrix e_hat = y.sub(y_hat); if (e_hat.normmax() > 0.02) { logger.warn("WARNING: Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } else { logger.debug("Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } if (logger.isDebugEnabled()) { logger.debug("REPORTING POLYFIT LEAST SQUARES ERRORS"); logger.debug(" time \t\t\t y \t\t\t yhat \t\t\t ehat"); for (int i = 0; i < numOfPoints; i++) { logger.debug(" " + t.get(i) + "\t" + y.get(i) + "\t" + y_hat.get(i) + "\t" + e_hat.get(i)); } for (int i = 0; i < numOfPoints - 1; i++) { double dt = t.get(i + 1) - t.get(i); logger.debug("Time step between point " + i + 1 + " and " + i + "= " + dt); if (Math.abs(dt - (t.get(1) - t.get(0))) > 0.001) logger.warn("WARNING: Orbit: data does not have equidistant time interval?"); } } return x.toArray(); }  }
PolyUtils { public static double[] polyFit(DoubleMatrix t, DoubleMatrix y, final int degree) throws IllegalArgumentException { setLoggerLevel(); if (t.length != y.length || !t.isVector() || !y.isVector()) { logger.error("polyfit: require same size vectors."); throw new IllegalArgumentException("polyfit: require same size vectors."); } final int numOfPoints = t.length; final int numOfUnknowns = degree + 1; logger.debug("Degree of interpolating polynomial: {}", degree); logger.debug("Number of unknowns: {}", numOfUnknowns); logger.debug("Number of data points: {}", numOfPoints); if (numOfPoints < numOfUnknowns) { logger.error("Number of points is smaller than parameters solved for."); throw new IllegalArgumentException("Number of points is smaller than parameters solved for."); } logger.debug("Setting up linear system of equations"); DoubleMatrix A = new DoubleMatrix(numOfPoints, numOfUnknowns); for (int j = 0; j <= degree; j++) { A.putColumn(j, pow(t, j)); } logger.debug("Solving lin. system of equations with Cholesky."); DoubleMatrix N = A.transpose().mmul(A); DoubleMatrix rhs = A.transpose().mmul(y); DoubleMatrix x = Solve.solveSymmetric(N, rhs); DoubleMatrix Qx_hat = Solve.solveSymmetric(N, DoubleMatrix.eye(N.getRows())); double maxDeviation = (N.mmul(Qx_hat).sub(DoubleMatrix.eye(Qx_hat.rows))).normmax(); logger.debug("polyfit orbit: max(abs(N*inv(N)-I)) = " + maxDeviation); if (maxDeviation > 1e-6) { logger.warn("polyfit orbit: max(abs(N*inv(N)-I)) = {}", maxDeviation); logger.warn("polyfit orbit interpolation unstable!"); } DoubleMatrix y_hat = A.mmul(x); DoubleMatrix e_hat = y.sub(y_hat); if (e_hat.normmax() > 0.02) { logger.warn("WARNING: Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } else { logger.debug("Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } if (logger.isDebugEnabled()) { logger.debug("REPORTING POLYFIT LEAST SQUARES ERRORS"); logger.debug(" time \t\t\t y \t\t\t yhat \t\t\t ehat"); for (int i = 0; i < numOfPoints; i++) { logger.debug(" " + t.get(i) + "\t" + y.get(i) + "\t" + y_hat.get(i) + "\t" + e_hat.get(i)); } for (int i = 0; i < numOfPoints - 1; i++) { double dt = t.get(i + 1) - t.get(i); logger.debug("Time step between point " + i + 1 + " and " + i + "= " + dt); if (Math.abs(dt - (t.get(1) - t.get(0))) > 0.001) logger.warn("WARNING: Orbit: data does not have equidistant time interval?"); } } return x.toArray(); }  static double normalize2(double data, final int min, final int max); static double normalize2(double data, final double min, final double max); static DoubleMatrix normalize(DoubleMatrix t); static int degreeFromCoefficients(int numOfCoefficients); static int numberOfCoefficients(final int degree); static double[] polyFitNormalized(DoubleMatrix t, DoubleMatrix y, final int degree); static double[] polyFit2D(final DoubleMatrix x, final DoubleMatrix y, final DoubleMatrix z, final int degree); static double[] polyFit(DoubleMatrix t, DoubleMatrix y, final int degree); static double polyVal1D(double x, double[] coeffs); static double[][] polyval(final double[] x, final double[] y, final double coeff[], int degree); static DoubleMatrix polyval(final DoubleMatrix x, final DoubleMatrix y, final DoubleMatrix coeff, int degree); static double polyval(final double x, final double y, final DoubleMatrix coeff, int degree); static double polyval(final double x, final double y, final double[] coeff); static double polyval(final double x, final double y, final double[] coeff, int degree); }
PolyUtils { public static double[] polyFit(DoubleMatrix t, DoubleMatrix y, final int degree) throws IllegalArgumentException { setLoggerLevel(); if (t.length != y.length || !t.isVector() || !y.isVector()) { logger.error("polyfit: require same size vectors."); throw new IllegalArgumentException("polyfit: require same size vectors."); } final int numOfPoints = t.length; final int numOfUnknowns = degree + 1; logger.debug("Degree of interpolating polynomial: {}", degree); logger.debug("Number of unknowns: {}", numOfUnknowns); logger.debug("Number of data points: {}", numOfPoints); if (numOfPoints < numOfUnknowns) { logger.error("Number of points is smaller than parameters solved for."); throw new IllegalArgumentException("Number of points is smaller than parameters solved for."); } logger.debug("Setting up linear system of equations"); DoubleMatrix A = new DoubleMatrix(numOfPoints, numOfUnknowns); for (int j = 0; j <= degree; j++) { A.putColumn(j, pow(t, j)); } logger.debug("Solving lin. system of equations with Cholesky."); DoubleMatrix N = A.transpose().mmul(A); DoubleMatrix rhs = A.transpose().mmul(y); DoubleMatrix x = Solve.solveSymmetric(N, rhs); DoubleMatrix Qx_hat = Solve.solveSymmetric(N, DoubleMatrix.eye(N.getRows())); double maxDeviation = (N.mmul(Qx_hat).sub(DoubleMatrix.eye(Qx_hat.rows))).normmax(); logger.debug("polyfit orbit: max(abs(N*inv(N)-I)) = " + maxDeviation); if (maxDeviation > 1e-6) { logger.warn("polyfit orbit: max(abs(N*inv(N)-I)) = {}", maxDeviation); logger.warn("polyfit orbit interpolation unstable!"); } DoubleMatrix y_hat = A.mmul(x); DoubleMatrix e_hat = y.sub(y_hat); if (e_hat.normmax() > 0.02) { logger.warn("WARNING: Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } else { logger.debug("Max. approximation error at datapoints (x,y,or z?): {}", e_hat.normmax()); } if (logger.isDebugEnabled()) { logger.debug("REPORTING POLYFIT LEAST SQUARES ERRORS"); logger.debug(" time \t\t\t y \t\t\t yhat \t\t\t ehat"); for (int i = 0; i < numOfPoints; i++) { logger.debug(" " + t.get(i) + "\t" + y.get(i) + "\t" + y_hat.get(i) + "\t" + e_hat.get(i)); } for (int i = 0; i < numOfPoints - 1; i++) { double dt = t.get(i + 1) - t.get(i); logger.debug("Time step between point " + i + 1 + " and " + i + "= " + dt); if (Math.abs(dt - (t.get(1) - t.get(0))) > 0.001) logger.warn("WARNING: Orbit: data does not have equidistant time interval?"); } } return x.toArray(); }  static double normalize2(double data, final int min, final int max); static double normalize2(double data, final double min, final double max); static DoubleMatrix normalize(DoubleMatrix t); static int degreeFromCoefficients(int numOfCoefficients); static int numberOfCoefficients(final int degree); static double[] polyFitNormalized(DoubleMatrix t, DoubleMatrix y, final int degree); static double[] polyFit2D(final DoubleMatrix x, final DoubleMatrix y, final DoubleMatrix z, final int degree); static double[] polyFit(DoubleMatrix t, DoubleMatrix y, final int degree); static double polyVal1D(double x, double[] coeffs); static double[][] polyval(final double[] x, final double[] y, final double coeff[], int degree); static DoubleMatrix polyval(final DoubleMatrix x, final DoubleMatrix y, final DoubleMatrix coeff, int degree); static double polyval(final double x, final double y, final DoubleMatrix coeff, int degree); static double polyval(final double x, final double y, final double[] coeff); static double polyval(final double x, final double y, final double[] coeff, int degree);  }