1 /*
2 * Copyright (C) 2015 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical;
17
18 import com.irurueta.algebra.Matrix;
19 import com.irurueta.algebra.WrongSizeException;
20
21 /**
22 * Class to estimate the Jacobian of a multi variate and multidimensional
23 * function.
24 * This class evaluates a function at very close locations of a given input
25 * point in order to determine the Jacobian at such point.
26 */
27 public class JacobianEstimator {
28
29 /**
30 * Constant considered as machine precision.
31 */
32 public static final double EPS = 1e-8;
33
34 /**
35 * Listener to evaluate a multivariate function.
36 */
37 private final MultiVariateFunctionEvaluatorListener listener;
38
39 /**
40 * Internal array to hold input parameter's values.
41 */
42 private double[] xh;
43
44 /**
45 * Constructor.
46 *
47 * @param listener listener to evaluate a multivariate function.
48 */
49 public JacobianEstimator(final MultiVariateFunctionEvaluatorListener listener) {
50 this.listener = listener;
51 }
52
53 /**
54 * Returns the Jacobian of a multivariate function at provided point.
55 *
56 * @param point input point.
57 * @return jacobian.
58 * @throws EvaluationException raised if function cannot be evaluated.
59 */
60 public Matrix jacobian(final double[] point) throws EvaluationException {
61 try {
62 final var result = new Matrix(listener.getNumberOfVariables(), point.length);
63 jacobian(point, result);
64 return result;
65 } catch (final WrongSizeException e) {
66 throw new EvaluationException(e);
67 }
68 }
69
70 /**
71 * Sets estimated jacobian in provided result matrix of a multivariate
72 * function at provided point.
73 * This method is preferred respect to jacobian(double[]) because result
74 * matrix can be reused and hence is more memory efficient.
75 *
76 * @param point input point.
77 * @param result output matrix containing estimated jacobian. This parameter
78 * must be a matrix having the same number of columns as the length of the
79 * point and the same number of rows as the number of variables returned by
80 * function evaluations.
81 * @throws EvaluationException raised if function cannot be evaluated.
82 * @throws IllegalArgumentException if size of result is not valid.
83 */
84 public void jacobian(final double[] point, final Matrix result) throws EvaluationException {
85 final var numdims = point.length;
86 final var numvars = listener.getNumberOfVariables();
87 if (result.getColumns() != numdims) {
88 throw new IllegalArgumentException();
89 }
90 if (result.getRows() != numvars) {
91 throw new IllegalArgumentException();
92 }
93
94 final var fold = new double[numvars];
95 final var fh = new double[numvars];
96 if (xh == null || xh.length != numdims) {
97 xh = new double[numdims];
98 }
99 System.arraycopy(point, 0, xh, 0, numdims);
100
101 double temp;
102 double h;
103 listener.evaluate(point, fold);
104 for (var j = 0; j < numdims; j++) {
105 temp = point[j];
106 h = EPS * Math.abs(temp);
107 if (h == 0.0) {
108 h = EPS;
109 }
110 xh[j] = temp + h;
111 h = xh[j] - temp;
112 listener.evaluate(xh, fh);
113 xh[j] = temp;
114 for (var i = 0; i < numvars; i++) {
115 result.setElementAt(i, j, (fh[i] - fold[i]) / h);
116 }
117 }
118 }
119 }