SumToOne { public double[][] jacobian(final double[] fitParms) { ArgumentChecker.isTrue(fitParms.length == _n - 1, "length of fitParms is {}, but must be {} ", fitParms.length, _n - 1); final double[] sin = new double[_n - 1]; final double[] cos = new double[_n - 1]; for (int j = 0; j < _n - 1; j++) { sin[j] = Math.sin(fitParms[j]); cos[j] = Math.cos(fitParms[j]); } final double[] a = new double[_n]; for (int i = 0; i < _n; i++) { double prod = 1.0; for (int j = 0; j < _n - 1; j++) { if (_set[i][j] == 1) { prod *= sin[j]; } else if (_set[i][j] == -1) { prod *= cos[j]; } } a[i] = 2 * prod * prod; } final double[][] res = new double[_n][_n - 1]; for (int i = 0; i < _n; i++) { for (int j = 0; j < _n - 1; j++) { if (_set[i][j] == 1 && a[i] != 0.0) { res[i][j] = a[i] * cos[j] / sin[j]; } else if (_set[i][j] == -1 && a[i] != 0.0) { res[i][j] = -a[i] * sin[j] / cos[j]; } } } return res; } SumToOne(final int n); double[] transform(final double[] fitParms); DoubleMatrix1D transform(final DoubleMatrix1D fitParms); double[] inverseTransform(final double[] modelParms); DoubleMatrix1D inverseTransform(final DoubleMatrix1D modelParms); double[][] jacobian(final double[] fitParms); DoubleMatrix2D jacobian(final DoubleMatrix1D fitParms);  }
@Test public void jacobianTest() { final int n = 5; final SumToOne trans = new SumToOne(n); Function1D<DoubleMatrix1D, DoubleMatrix1D> func = new Function1D<DoubleMatrix1D, DoubleMatrix1D>() { @Override public DoubleMatrix1D evaluate(DoubleMatrix1D theta) { return trans.transform(theta); } }; Function1D<DoubleMatrix1D, DoubleMatrix2D> jacFunc = new Function1D<DoubleMatrix1D, DoubleMatrix2D>() { @Override public DoubleMatrix2D evaluate(DoubleMatrix1D theta) { return trans.jacobian(theta); } }; Function1D<DoubleMatrix1D, DoubleMatrix2D> fdJacFunc = DIFFER.differentiate(func); double[] theta = new double[n - 1]; for (int tries = 0; tries < 10; tries++) { for (int i = 0; i < n - 1; i++) { theta[i] = RANDOM.nextDouble(); } DoubleMatrix1D vTheta = new DoubleMatrix1D(theta); DoubleMatrix2D jac = jacFunc.evaluate(vTheta); DoubleMatrix2D fdJac = fdJacFunc.evaluate(vTheta); for (int j = 0; j < n - 1; j++) { double sum = 0.0; for (int i = 0; i < n; i++) { sum += jac.getEntry(i, j); assertEquals("element " + i + " " + j, fdJac.getEntry(i, j), jac.getEntry(i, j), 1e-6); } assertEquals("wrong sum of sensitivities", 0.0, sum, 1e-15); } } }