ningshuxia
2022-12-01 ad494f13d2ddf31f142cf7fb908b3a6e90395a1a
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
// <copyright file="NumericalJacobian.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
 
using System;
 
namespace IStation.Numerics.Differentiation
{
    /// <summary>
    /// Class for evaluating the Jacobian of a function using finite differences.
    /// By default, a central 3-point method is used.
    /// </summary>
    public class NumericalJacobian
    {
        /// <summary>
        /// Number of function evaluations.
        /// </summary>
        public int FunctionEvaluations => _df.Evaluations;
 
        private readonly NumericalDerivative _df;
 
        /// <summary>
        /// Creates a numerical Jacobian object with a three point central difference method.
        /// </summary>
        public NumericalJacobian() : this(3, 1) { }
 
        /// <summary>
        /// Creates a numerical Jacobian with a specified differentiation scheme.
        /// </summary>
        /// <param name="points">Number of points for Jacobian evaluation.</param>
        /// <param name="center">Center point for differentiation.</param>
        public NumericalJacobian(int points, int center)
        {
            _df = new NumericalDerivative(points, center);
        }
 
        /// <summary>
        /// Evaluates the Jacobian of scalar univariate function f at point x.
        /// </summary>
        /// <param name="f">Scalar univariate function handle.</param>
        /// <param name="x">Point at which to evaluate Jacobian.</param>
        /// <returns>Jacobian vector.</returns>
        public double[] Evaluate(Func<double, double> f, double x)
        {
            return new[] { _df.EvaluateDerivative(f, x, 1) };
        }
 
        /// <summary>
        /// Evaluates the Jacobian of a multivariate function f at vector x.
        /// </summary>
        /// <remarks>
        /// This function assumes that the length of vector x consistent with the argument count of f.
        /// </remarks>
        /// <param name="f">Multivariate function handle.</param>
        /// <param name="x">Points at which to evaluate Jacobian.</param>
        /// <returns>Jacobian vector.</returns>
        public double[] Evaluate(Func<double[], double> f, double[] x)
        {
            var jacobian = new double[x.Length];
 
            for (var i = 0; i < jacobian.Length; i++)
                jacobian[i] = _df.EvaluatePartialDerivative(f, x, i, 1);
 
            return jacobian;
        }
 
        /// <summary>
        /// Evaluates the Jacobian of a multivariate function f at vector x given a current function value.
        /// </summary>
        /// <remarks>
        /// To minimize the number of function evaluations, a user can supply the current value of the function
        /// to be used in computing the Jacobian. This value must correspond to the "center" location for the
        /// finite differencing. If a scheme is used where the center value is not evaluated, this will provide no
        /// added efficiency. This method also assumes that the length of vector x consistent with the argument count of f.
        /// </remarks>
        /// <param name="f">Multivariate function handle.</param>
        /// <param name="x">Points at which to evaluate Jacobian.</param>
        /// <param name="currentValue">Current function value at finite difference center.</param>
        /// <returns>Jacobian vector.</returns>
        public double[] Evaluate(Func<double[], double> f, double[] x, double currentValue)
        {
            var jacobian = new double[x.Length];
 
            for (var i = 0; i < jacobian.Length; i++)
                jacobian[i] = _df.EvaluatePartialDerivative(f, x, i, 1, currentValue);
 
            return jacobian;
        }
 
        /// <summary>
        /// Evaluates the Jacobian of a multivariate function array f at vector x.
        /// </summary>
        /// <param name="f">Multivariate function array handle.</param>
        /// <param name="x">Vector at which to evaluate Jacobian.</param>
        /// <returns>Jacobian matrix.</returns>
        public double[,] Evaluate(Func<double[], double>[] f, double[] x)
        {
            var jacobian = new double[f.Length, x.Length];
            for (int i = 0; i < f.Length; i++)
            {
                var gradient = Evaluate(f[i], x);
                for (int j = 0; j < gradient.Length; j++)
                    jacobian[i, j] = gradient[j];
            }
 
            return jacobian;
        }
 
        /// <summary>
        /// Evaluates the Jacobian of a multivariate function array f at vector x given a vector of current function values.
        /// </summary>
        /// <remarks>
        /// To minimize the number of function evaluations, a user can supply a vector of current values of the functions
        /// to be used in computing the Jacobian. These value must correspond to the "center" location for the
        /// finite differencing. If a scheme is used where the center value is not evaluated, this will provide no
        /// added efficiency. This method also assumes that the length of vector x consistent with the argument count of f.
        /// </remarks>
        /// <param name="f">Multivariate function array handle.</param>
        /// <param name="x">Vector at which to evaluate Jacobian.</param>
        /// <param name="currentValues">Vector of current function values.</param>
        /// <returns>Jacobian matrix.</returns>
        public double[,] Evaluate(Func<double[], double>[] f, double[] x, double[] currentValues)
        {
            var jacobian = new double[f.Length, x.Length];
            for (int i = 0; i < f.Length; i++)
            {
                var gradient = Evaluate(f[i], x, currentValues[i]);
                for (int j = 0; j < gradient.Length; j++)
                    jacobian[i, j] = gradient[j];
            }
 
            return jacobian;
        }
 
        /// <summary>
        /// Resets the function evaluation counter for the Jacobian.
        /// </summary>
        public void ResetFunctionEvaluations()
        {
            _df.ResetEvaluations();
        }
    }
}