PenaltyMatrixGenerator { public static DoubleMatrix2D getDerivativeMatrix(double[] x, int k, boolean includeEnds) { ArgumentChecker.notEmpty(x, "x"); ArgumentChecker.notNegative(k, "k"); int size = x.length; ArgumentChecker.isTrue(k < size, "order too high. Length of x is {}, and k is {}", size, k); if (k == 0) { return new IdentityMatrix(size); } else if (k > 2) { throw new NotImplementedException("cannot handle order (k) > 2"); } ArgumentChecker.isTrue(size > 2, "Need at least 3 points for a three point estimate"); double[] dx = new double[size - 1]; double[] dx2 = new double[size - 1]; for (int i = 0; i < (size - 1); i++) { double temp = x[i + 1] - x[i]; ArgumentChecker.isTrue(temp > 0.0, "x not in ascending order, or two identical points"); dx[i] = temp; dx2[i] = temp * temp; } double[] w = new double[size - 2]; for (int i = 0; i < (size - 2); i++) { w[i] = 1.0 / (dx[i] * dx[i + 1] * (dx[i] + dx[i + 1])); } DoubleMatrix2D res = new DoubleMatrix2D(size, size); double[][] data = res.getData(); if (k == 1) { for (int i = 1; i < (size - 1); i++) { data[i][i - 1] = -w[i - 1] * dx2[i]; data[i][i] = w[i - 1] * (dx2[i] - dx2[i - 1]); data[i][i + 1] = w[i - 1] * dx2[i - 1]; } if (includeEnds) { data[0][0] = -w[0] * dx[1] * (2 * dx[0] + dx[1]); data[0][1] = w[0] * (dx2[0] + dx2[1] + 2 * dx[0] * dx[1]); data[0][2] = -w[0] * dx2[0]; data[size - 1][size - 3] = w[size - 3] * dx2[size - 2]; data[size - 1][size - 2] = -w[size - 3] * (dx2[size - 3] + dx2[size - 2] + 2 * dx[size - 2] * dx[size - 3]); data[size - 1][size - 1] = w[size - 3] * dx[size - 3] * (2 * dx[size - 2] + dx[size - 3]); } return res; } else { for (int i = 1; i < (size - 1); i++) { double tmp = 2 * w[i - 1]; data[i][i - 1] = tmp * dx[i]; data[i][i] = -tmp * (dx[i] + dx[i - 1]); data[i][i + 1] = tmp * dx[i - 1]; } if (includeEnds) { data[0] = data[1]; data[size - 1] = data[size - 2]; } return res; } }  static DoubleMatrix2D getDifferenceMatrix(int m, int k); static DoubleMatrix2D getPenaltyMatrix(int m, int k); static DoubleMatrix2D getPenaltyMatrix(int[] numElements, int k, int index); static DoubleMatrix2D getPenaltyMatrix(int[] numElements, int[] k, double[] lambda); static DoubleMatrix2D getDerivativeMatrix(double[] x, int k, boolean includeEnds); static DoubleMatrix2D getPenaltyMatrix(double[] x, int k); static DoubleMatrix2D getPenaltyMatrix(double[][] x, int k, int index); static DoubleMatrix2D getPenaltyMatrix(double[][] x, int[] k, double[] lambda); static DoubleMatrix1D flattenMatrix(DoubleMatrix2D aMatrix); static DoubleMatrix2D getMatrixForFlattened(int[] numElements, DoubleMatrix2D m, int index);  }
@Test(dataProvider = "data") public void derivativeMatrix1DTest(double[] x, DoubleMatrix1D y, DoubleMatrix1D dydx, DoubleMatrix1D d2ydx2) { int n = x.length; DoubleMatrix2D d0 = getDerivativeMatrix(x, 0, true); AssertMatrix.assertEqualsMatrix(new IdentityMatrix(n), d0, 1e-14); DoubleMatrix2D d1 = getDerivativeMatrix(x, 1, true); DoubleMatrix2D d2 = getDerivativeMatrix(x, 2, true); DoubleMatrix1D d1y = (DoubleMatrix1D) MA.multiply(d1, y); DoubleMatrix1D d2y = (DoubleMatrix1D) MA.multiply(d2, y); AssertMatrix.assertEqualsVectors(dydx, d1y, 1e-13); AssertMatrix.assertEqualsVectors(d2ydx2, d2y, 1e-13); }
@Test(expectedExceptions = NotImplementedException.class) public void orderG2Test() { double[] x = new double[] {0.0, 0.3, 0.7, 0.8, 1.2, 2.0 }; @SuppressWarnings("unused") DoubleMatrix2D d3 = getDerivativeMatrix(x, 3, true); }
@Test(expectedExceptions = IllegalArgumentException.class) public void twoPointsTest() { double[] x = new double[] {0.2, 0.5 }; @SuppressWarnings("unused") DoubleMatrix2D d1 = getDerivativeMatrix(x, 1, true); }