From 3ed9bf97984587f0d815496de1b37190bd9e21ab Mon Sep 17 00:00:00 2001 From: Hythem Sidky Date: Fri, 9 Jan 2015 21:18:33 -0500 Subject: [PATCH 1/2] Added numerical differentiation classes. --- .../FiniteDifferenceCoefficients.cs | 140 +++++ .../Differentiation/NumericalDerivative.cs | 483 ++++++++++++++++++ src/Numerics/Numerics.csproj | 2 + .../FiniteDifferenceCoefficientsTests.cs | 77 +++ .../NumericalDerivativeTests.cs | 202 ++++++++ src/UnitTests/UnitTests.csproj | 2 + 6 files changed, 906 insertions(+) create mode 100644 src/Numerics/Differentiation/FiniteDifferenceCoefficients.cs create mode 100644 src/Numerics/Differentiation/NumericalDerivative.cs create mode 100644 src/UnitTests/DifferentiationTests/FiniteDifferenceCoefficientsTests.cs create mode 100644 src/UnitTests/DifferentiationTests/NumericalDerivativeTests.cs diff --git a/src/Numerics/Differentiation/FiniteDifferenceCoefficients.cs b/src/Numerics/Differentiation/FiniteDifferenceCoefficients.cs new file mode 100644 index 00000000..6c5de78f --- /dev/null +++ b/src/Numerics/Differentiation/FiniteDifferenceCoefficients.cs @@ -0,0 +1,140 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// 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. +// + +using System; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.Differentiation +{ + /// + /// Class to calculate finite difference coefficients using Taylor series expansion method. + /// + /// + /// For n points, coefficients are calculated up to the maximum derivative order possible (n-1). + /// The current function value position specifies the "center" for surrounding coefficients. + /// Selecting the first, middle or last positions represent forward, backwards and central difference methods. + /// + /// + /// + public class FiniteDifferenceCoefficients + { + /// + /// Number of points for finite difference coefficients. Changing this value recalculates the coefficients table. + /// + public int Points + { + get { return _points; } + set + { + CalculateCoefficients(value); + _points = value; + } + } + + private double[][,] _coefficients; + private int _points; + + /// + /// Initializes a new instance of the class. + /// + /// Number of finite difference coefficients. + public FiniteDifferenceCoefficients(int points) + { + Points = points; + CalculateCoefficients(Points); + } + + /// + /// Gets the finite difference coefficients for a specified center and order. + /// + /// Current function position with respect to coefficients. Must be within point range. + /// Order of finite difference coefficients. + /// Vector of finite difference coefficients. + public double[] GetCoefficients(int center, int order) + { + if (center >= _coefficients.Length) + throw new ArgumentOutOfRangeException("center", "Center position must be within the point range."); + if (order >= _coefficients.Length) + throw new ArgumentOutOfRangeException("order", "Maximum difference order is points-1."); + + // Return proper row + var columns = _coefficients[center].GetLength(1); + var array = new double[columns]; + for (int i = 0; i < columns; ++i) + array[i] = _coefficients[center][order, i]; + return array; + } + + /// + /// Gets the finite difference coefficients for all orders at a specified center. + /// + /// Current function position with respect to coefficients. Must be within point range. + /// Rectangular array of coefficients, with columns specifing order. + public double[,] GetCoefficientsForAllOrders(int center) + { + if (center >= _coefficients.Length) + throw new ArgumentOutOfRangeException("center", "Center position must be within the point range."); + + return _coefficients[center]; + } + + private void CalculateCoefficients(int points) + { + var c = new double[points][,]; + + // For ever possible center given the number of points, compute ever possible coefficeint for all possible orders. + for (int center = 0; center < points; center++) + { + // Deltas matrix for center located at 'center'. + var A = new DenseMatrix(points); + var l = points - center - 1; + for (int row = points - 1; row >= 0; row--) + { + A[row, 0] = 1.0; + for (int col = 1; col < points; col++) + { + A[row, col] = A[row, col - 1] * l / col; + } + l -= 1; + } + + c[center] = A.Inverse().ToArray(); + + // "Polish" results by rounding. + var fac = SpecialFunctions.Factorial(points); + for (int j = 0; j < points; j++) + for (int k = 0; k < points; k++) + c[center][j, k] = (Math.Round(c[center][j, k] * fac, MidpointRounding.AwayFromZero)) / fac; + } + + _coefficients = c; + } + } +} diff --git a/src/Numerics/Differentiation/NumericalDerivative.cs b/src/Numerics/Differentiation/NumericalDerivative.cs new file mode 100644 index 00000000..a0bfbf9a --- /dev/null +++ b/src/Numerics/Differentiation/NumericalDerivative.cs @@ -0,0 +1,483 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// 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. +// + +using System; +using System.Linq; + + +namespace MathNet.Numerics.Differentiation +{ + /// + /// Type of finite different step size. + /// + public enum StepType + { + /// + /// The absolute step size value will be used in numerical derivatives, regardless of order or function parameters. + /// + Absolute, + + /// + /// A base step size value, h, will be scaled according to the function input parameter. A common example is hx = h*(1+abs(x)), however + /// this may vary depending on implementation. This definition only guarantees that the only scaling will be relative to the + /// function input parameter and not the order of the finite difference derivative. + /// + RelativeX, + + /// + /// A base step size value, eps (typically machine precision), is scaled according to the finite difference coefficient order + /// and function input parameter. The initial scaling according to finite different coefficient order can be thought of as producing a + /// base step size, h, that is equivalent to scaling. This stepsize is then scaled according to the function + /// input parameter. Although implementation may vary, an example of second order accurate scaling may be (eps)^(1/3)*(1+abs(x)). + /// + Relative + }; + + /// + /// Class to evaluate the numerical derivative of a function using finite difference approximations. + /// Variable point and center methods can be initialized . + /// This class can also be used to return function handles (delagates) for a fixed derivative order and variable. + /// It is possible to evaluate the derivative and partial derivative of univariate and multivariate functions respectively. + /// + public class NumericalDerivative + { + /// + /// Sets and gets the finite difference step size. This value is for each function evaluation if relative stepsize types are used. + /// If the base step size used in scaling is desired, see . + /// + /// + /// Setting then getting the StepSize may return a different value. This is not unusual since a user-defined step size is converted to a + /// base-2 representable number to improve finite difference accuracy. + /// + public double StepSize + { + get { return _stepSize; } + set + { + //Base 2 yields more accurate results... + var p = Math.Log(Math.Abs(value))/Math.Log(2); + _stepSize = Math.Pow(2, Math.Round(p)); + } + } + + /// + /// Sets and gets the base fininte difference step size. This assigned value to this parameter is only used if is set to RelativeX. + /// However, if the StepType is Relative, it will contain the base step size computed from based on the finite difference order. + /// + public double BaseStepSize + { + get { return _baseStepSize; } + set + { + //Base 2 yields more accurate results... + var p = Math.Log(Math.Abs(value)) / Math.Log(2); + _baseStepSize = Math.Pow(2, Math.Round(p)); + } + } + + /// + /// Sets and gets the base finite difference step size. This parameter is only used if is set to Relative. + /// By default this is set to machine epsilon, from which is computed. + /// + public double Epsilon + { + get { return _epsilon; } + set + { + //Base 2 yields more accurate results... + var p = Math.Log(Math.Abs(value)) / Math.Log(2); + _epsilon = Math.Pow(2, Math.Round(p)); + } + } + + /// + /// Sets and gets the location of the center point for the finite difference derivative. + /// + public int Center + { + get { return _center; } + set + { + if (value >= _points || value < 0) + throw new ArgumentOutOfRangeException("value", "Center must lie between 0 and points -1"); + _center = value; + } + } + + /// + /// Number of times a function is evaluated for numerical derivatives. + /// + public int Evaluations { get; private set; } + + /// + /// Type of step size for computing finite differences. If set to absolute, dx = h. + /// If set to relative, dx = (1+abs(x))*h^(2/(order+1)). This provides accurate results when + /// h is approximately equal to the square-root of machine accuracy, epsilon. + /// + public StepType StepType + { + get { return _stepType; } + set { _stepType = value; } + } + + private readonly int _points; + private int _center; + private double _stepSize = Math.Pow(2, -10); + private double _epsilon = Math.Pow(2, -52); + private double _baseStepSize = Math.Pow(2, -26); + private StepType _stepType = StepType.Relative; + private readonly FiniteDifferenceCoefficients _coefficients; + + /// + /// Initializes a NumericalDerivative class with the default 3 point center difference method. + /// + public NumericalDerivative() : this(3, 1) + { + } + + /// + /// Initialized a NumericalDerivative class. + /// + /// Number of points for finite difference derivatives. + /// Location of the center with respect to other points. Value ranges from zero to points-1. + public NumericalDerivative(int points, int center) + { + _center = center; + if (points < 2) + throw new ArgumentOutOfRangeException("points", "Points must be two or greater."); + _points = points; + Center = center; + _epsilon = CalculateMachineEpsilon(); + _coefficients = new FiniteDifferenceCoefficients(points); + } + + /// + /// Evaluates the derivative of equidistant points using the finite difference method. + /// + /// Vector of points StepSize apart. + /// Derivative order. + /// Finite difference step size. + /// Derivative of points of the specified order. + public double EvaluateDerivative(double[] points, int order, double stepSize) + { + if (points == null) + throw new ArgumentNullException("points"); + + if (order >= _points || order < 0) + throw new ArgumentOutOfRangeException("order", "Order must be between zero and points-1."); + + var c = _coefficients.GetCoefficients(Center, order); + var result = c.Select((t, i) => t*points[i]).Sum(); + result /= Math.Pow(stepSize, order); + return result; + } + + /// + /// Evaluates the derivative of a scalar univariate function. + /// + /// + /// Supplying the optional argument currentValue will reduce the number of function evaluations + /// required to calculate the finite difference derivative. + /// + /// Function handle. + /// Point at which to compute the derivative. + /// Derivative order. + /// Current function value at center. + /// Function derivative at x of the specified order. + public double EvaluateDerivative(Func f, double x, int order, double? currentValue = null) + { + var c = _coefficients.GetCoefficients(Center, order); + var h = CalculateStepSize(_points, x, order); + + var points = new double[_points]; + for (int i = 0; i < _points; i++) + { + if (i == Center && currentValue.HasValue) + points[i] = currentValue.Value; + else if(c[i] != 0) // Only evaluate function if it will actually be used. + { + points[i] = f(x + (i - Center) * h); + Evaluations++; + } + } + + return EvaluateDerivative(points, order, h); + } + + /// + /// Creates a function handle for the derivative of a scalar univariate function. + /// + /// Input function handle. + /// Derivative order. + /// Function handle that evaluates the derivative of input function at a fixed order. + public Func CreateDerivativeFunctionHandle(Func f, int order) + { + return x => EvaluateDerivative(f, x, order); + } + + /// + /// Evaluates the partial derivative of a multivariate function. + /// + /// Multivariate function handle. + /// Vector at which to evaluate the derivative. + /// Index of independent variable for partial derivative. + /// Derivative order. + /// Current function value at center. + /// Function partial derivative at x of the specified order. + public double EvaluatePartialDerivative(Func f, double[] x, int parameterIndex, int order, double? currentValue = null) + { + var xi = x[parameterIndex]; + var c = _coefficients.GetCoefficients(Center, order); + var h = CalculateStepSize(_points, x[parameterIndex], order); + + var points = new double[_points]; + for (int i = 0; i < _points; i++) + { + if (i == Center && currentValue.HasValue) + points[i] = currentValue.Value; + else if(c[i] != 0) // Only evaluate function if it will actually be used. + { + x[parameterIndex] = xi + (i - Center) * h; + points[i] = f(x); + Evaluations++; + } + } + + //restore original value + x[parameterIndex] = xi; + return EvaluateDerivative(points, order, h); + } + + /// + /// Evaluates the partial derivatives of a multivariate function array. + /// + /// + /// This function assumes the input vector x is of the correct length for f. + /// + /// Multivariate vector function array handle. + /// Vector at which to evaluate the derivatives. + /// Index of independent variable for partial derivative. + /// Derivative order. + /// Current function value at center. + /// Vector of functions partial derivatives at x of the specified order. + public double[] EvaluatePartialDerivative(Func[] f, double[] x, int parameterIndex, int order, double?[] currentValue = null) + { + var df = new double[f.Length]; + for (int i = 0; i < f.Length; i++) + { + if(currentValue != null && currentValue[i].HasValue) + df[i] = EvaluatePartialDerivative(f[i], x, parameterIndex, order, currentValue[i].Value); + else + df[i] = EvaluatePartialDerivative(f[i], x, parameterIndex, order); + } + + return df; + } + + /// + /// Creates a function handle for the partial derivative of a multivariate function. + /// + /// Input function handle. + /// Index of the independent variable for partial derivative. + /// Derivative order. + /// Function handle that evaluates partial derivative of input function at a fixed order. + public Func CreatePartialDerivativeFunctionHandle(Func f, int parameterIndex, + int order) + { + return x => EvaluatePartialDerivative(f, x, parameterIndex, order); + } + + /// + /// Creates a function handle for the partial derivative of a vector multivariate function. + /// + /// Input function handle. + /// Index of the independent variable for partial derivative. + /// Derivative order. + /// Function handle that evaluates partial derivative of input function at fixed order. + public Func CreatePartialDerivativeFunctionHandle(Func[] f, + int parameterIndex, + int order) + { + return x => EvaluatePartialDerivative(f, x, parameterIndex, order); + } + + /// + /// Evaluates the mixed partial derivative of variable order for multivariate functions. + /// + /// + /// This function recursively uses to evaluate mixed partial derivative. + /// Therefore, it is more efficient to call for higher order derivatives of + /// a single independent variable. + /// + /// Multivariate function handle. + /// Points at which to evaluate the derivative. + /// Vector of indices for the independent variables at descending derivative orders. + /// Highest order of differentiation. + /// Current function value at center. + /// Function mixed partial derivative at x of the specified order. + public double EvaluateMixedPartialDerivative(Func f, double[] x, int[] parameterIndex, + int order, double? currentValue = null) + { + if (parameterIndex.Length != order) + throw new ArgumentOutOfRangeException("parameterIndex", + "The number of parameters must match derivative order."); + + if (order == 1) + return EvaluatePartialDerivative(f, x, parameterIndex[0], order, currentValue); + + int reducedOrder = order - 1; + var reducedParameterIndex = new int[reducedOrder]; + Array.Copy(parameterIndex, 0, reducedParameterIndex, 0, reducedOrder); + + var points = new double[_points]; + var currentParameterIndex = parameterIndex[order - 1]; + var h = CalculateStepSize(_points, x[currentParameterIndex], order); + + var xi = x[currentParameterIndex]; + for (int i = 0; i < _points; i++) + { + x[currentParameterIndex] = xi + (i - Center)*h; + points[i] = EvaluateMixedPartialDerivative(f, x, reducedParameterIndex, reducedOrder); + } + + // restore original value + x[currentParameterIndex] = xi; + + // This will always be to the first order + return EvaluateDerivative(points, 1, h); + } + + /// + /// Evaluates the mixed partial derivative of variable order for multivariate function arrays. + /// + /// + /// This function recursively uses to evaluate mixed partial derivative. + /// Therefore, it is more efficient to call for higher order derivatives of + /// a single independent variable. + /// + /// Multivariate function array handle. + /// Vector at which to evaluate the derivative. + /// Vector of indices for the independent variables at descending derivative orders. + /// Highest order of differentiation. + /// Current function value at center. + /// Function mixed partial derivatives at x of the specified order. + public double[] EvaluateMixedPartialDerivative(Func[] f, double[] x, int[] parameterIndex, + int order, double?[] currentValue = null) + { + var df = new double[f.Length]; + for (int i = 0; i < f.Length; i++) + { + if(currentValue != null && currentValue[i].HasValue) + df[i] = EvaluateMixedPartialDerivative(f[i], x, parameterIndex, order, currentValue[i].Value); + else + df[i] = EvaluateMixedPartialDerivative(f[i], x, parameterIndex, order); + } + + return df; + } + + /// + /// Creates a function handle for the mixed partial derivative of a multivariate function. + /// + /// Input function handle. + /// Vector of indices for the independent variables at descending derivative orders. + /// Highest derivative order. + /// Function handle that evaluates the fixed mixed partial derivative of input function at fixed order. + public Func CreateMixedPartialDerivativeFunctionHandle(Func f, + int[] parameterIndex, int order) + { + return x => EvaluateMixedPartialDerivative(f, x, parameterIndex, order); + } + + /// + /// Creates a function handle for the mixed partial derivative of a multivariate vector function. + /// + /// Input vector function handle. + /// Vector of indices for the independent variables at descending derivative orders. + /// Highest derivative order. + /// Function handle that evaluates the fixed mixed partial derivative of input function at fixed order. + public Func CreateMixedPartialDerivativeFunctionHandle(Func[] f, + int[] parameterIndex, int order) + { + return x => EvaluateMixedPartialDerivative(f, x, parameterIndex, order); + } + + /// + /// Resets the evaluation counter. + /// + public void ResetEvaluations() + { + Evaluations = 0; + } + + /// + /// Calculates machine epsilon - the smallest number that can be added to 1, yeilding a results different than 1. + /// This is also known as roundoff error. + /// + /// Machine epislon + public static double CalculateMachineEpsilon() + { + double eps = 1; + + while ((1.0d + (eps / 2.0d)) > 1.0d) + eps /= 2.0d; + + return eps; + } + + private double[] CalculateStepSize(int points, double[] x, double order) + { + var h = new double[x.Length]; + for (int i = 1; i < h.Length; i++) + h[i] = CalculateStepSize(points, x[i], order); + + return h; + } + + private double CalculateStepSize(int points, double x, double order) + { + // Step size relative to function input parameter + if (StepType == StepType.RelativeX) + { + StepSize = BaseStepSize*(1 + Math.Abs(x)); + } + // Step size relative to function input parameter and order + else if (StepType == StepType.Relative) + { + var accuracy = points - order; + BaseStepSize = Math.Pow(Epsilon,(1/(accuracy + order))); + StepSize = BaseStepSize*(1 + Math.Abs(x)); + } + // Do nothing for absolute step size. + + return StepSize; + } + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index d88a75ae..166b4258 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -80,6 +80,8 @@ + + diff --git a/src/UnitTests/DifferentiationTests/FiniteDifferenceCoefficientsTests.cs b/src/UnitTests/DifferentiationTests/FiniteDifferenceCoefficientsTests.cs new file mode 100644 index 00000000..dfeb0006 --- /dev/null +++ b/src/UnitTests/DifferentiationTests/FiniteDifferenceCoefficientsTests.cs @@ -0,0 +1,77 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// 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. +// + + +namespace MathNet.Numerics.UnitTests.DifferentiationTests +{ + using Differentiation; + using NUnit.Framework; + + [TestFixture, Category("Differentiation")] + public class FiniteDifferenceCoefficientsTests + { + [Test] + public void CentralDifferenceFirstOrderThreePointTest() + { + double[] results = { -0.5, 0, 0.5 }; + var finite = new FiniteDifferenceCoefficients(3); + var coeff = finite.GetCoefficients(1, 1); + Assert.AreEqual(results, coeff); + } + + [Test] + public void CentralDifferenceSecondOrderFivePointsTest() + { + double[] results = { (double)-1 / 12, (double)4 / 3, (double)-5 / 2, (double)4 / 3, (double)-1 / 12 }; + var finite = new FiniteDifferenceCoefficients(5); + var coeff = finite.GetCoefficients(2, 2); + for (int i = 0; i < coeff.Length; i++) + Assert.AreEqual(results[i], coeff[i]); + } + + [Test] + public void ForwardDifferenceThirdOrderEightPointsTest() + { + double[] results = { (double)-967 / 120, (double)638 / 15, (double)-3929 / 40, (double)389 / 3, + (double)-2545 / 24, (double)268 / 5, (double)-1849 / 120, (double)29 / 15 }; + var finite = new FiniteDifferenceCoefficients(8); + var coeff = finite.GetCoefficients(0, 3); + for (int i = 0; i < coeff.Length; i++) + Assert.AreEqual(results[i], coeff[i]); + } + + [Test] + public void BackwardDifferenceThirdOrderFourPointsTest() + { + double[] results = { -1, 3, -3, 1 }; + var finite = new FiniteDifferenceCoefficients(4); + var coeff = finite.GetCoefficients(3, 3); + for (int i = 0; i < coeff.Length; i++) + Assert.AreEqual(results[i], coeff[i]); + + } + } +} diff --git a/src/UnitTests/DifferentiationTests/NumericalDerivativeTests.cs b/src/UnitTests/DifferentiationTests/NumericalDerivativeTests.cs new file mode 100644 index 00000000..c8bf2c59 --- /dev/null +++ b/src/UnitTests/DifferentiationTests/NumericalDerivativeTests.cs @@ -0,0 +1,202 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// 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. +// + +namespace MathNet.Numerics.UnitTests.DifferentiationTests +{ + using System; + using Differentiation; + using NUnit.Framework; + + [TestFixture, Category("Differentiation")] + class NumericalDerivativeTests + { + [Test] + public void SinFirstDerivativeAtZeroTest() + { + Func f = Math.Sin; + var df = new NumericalDerivative(); + Assert.AreEqual(1, df.EvaluateDerivative(f, 0, 1), 1e-10); + } + + [Test] + public void CubicPolynomialThirdDerivativeAtAnyValTest() + { + Func f = x => 3 * Math.Pow(x, 3) + 2 * x - 6; + var df = new NumericalDerivative(5, 2); + Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3)); + Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3)); + df.Center = 0; + Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3)); + Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3)); + df.Center = 1; + Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3)); + Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3)); + df.Center = 2; + Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3)); + Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3)); + df.Center = 3; + Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3)); + Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3)); + df.Center = 4; + Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3)); + Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3)); + } + + [Test] + public void CubicPolynomialFunctionValueTest() + { + Func f = x => 3 * Math.Pow(x, 3) + 2 * x - 6; + var current = f(2); + var df = new NumericalDerivative(3, 0); + Assert.AreEqual(38, df.EvaluateDerivative(f, 2, 1, current), 1e-8); + } + + [Test] + public void CreateDerivativeFunctionHandleTest() + { + Func f = x => 3 * Math.Pow(x, 3) + 2 * x - 6; + var nd = new NumericalDerivative(5, 2); + var df = nd.CreateDerivativeFunctionHandle(f, 3); + + Assert.AreEqual(18, df(0)); + + // Test new function with same nd class + Func f2 = x => 2 * Math.Pow(x, 3) + 2 * x - 6; + var df2 = nd.CreateDerivativeFunctionHandle(f2, 3); + + Assert.AreEqual(12, df2(0)); + + // Original delegate not changed + Assert.AreEqual(18, df(0)); + + } + + [Test] + public void ExponentialFunctionPartialDerivativeTest() + { + //Test Function + Func f = (x) => Math.Sin(x[0] * x[1]) + Math.Exp(-x[0] / 2) + x[1] / x[0]; + + //Analytical partial dfdx + Func dfdx = + (x) => Math.Cos(x[0] * x[1]) * x[1] - Math.Exp(-x[0] / 2) / 2 - x[1] / Math.Pow(x[0], 2); + + //Analytical partial dfdy + Func dfdy = (x) => Math.Cos(x[0] * x[1]) * x[0] + 1 / x[0]; + + var df = new NumericalDerivative(3, 1); + var x1 = new double[] { 3, 3 }; + Assert.AreEqual(dfdx(x1), df.EvaluatePartialDerivative(f, x1, 0, 1), 1e-8); + + Assert.AreEqual(dfdy(x1), df.EvaluatePartialDerivative(f, x1, 1, 1), 1e-8); + + var x2 = new double[] { 300, -50 }; + df.StepType = StepType.Absolute; + Assert.AreEqual(dfdx(x2), df.EvaluatePartialDerivative(f, x2, 0, 1), 1e-5); + Assert.AreEqual(dfdy(x2), df.EvaluatePartialDerivative(f, x2, 1, 1), 1e-2); + } + + [Test] + public void ExponentialFunctionPartialDerivativeCurrentValueTest() + { + //Test Function + Func f = (x) => Math.Sin(x[0] * x[1]) + Math.Exp(-x[0] / 2) + x[1] / x[0]; + + //Analytical partial dfdx + Func dfdx = + (x) => Math.Cos(x[0] * x[1]) * x[1] - Math.Exp(-x[0] / 2) / 2 - x[1] / Math.Pow(x[0], 2); + + //Analytical partial dfdy + Func dfdy = (x) => Math.Cos(x[0] * x[1]) * x[0] + 1 / x[0]; + + // Current value + var x1 = new double[] { 3, 3 }; + var current = f(x1); + + var df = new NumericalDerivative(5, 2); + Assert.AreEqual(dfdx(x1), df.EvaluatePartialDerivative(f, x1, 0, 1, current), 1e-8); + + Assert.AreEqual(dfdy(x1), df.EvaluatePartialDerivative(f, x1, 1, 1, current), 1e-8); + } + + [Test] + public void RosenbrockFunctionMixedDerivativeOneVariableSecondOrderTest() + { + Func f = x => Math.Pow(1 - x[0], 2) + 100 * Math.Pow(x[1] - Math.Pow(x[0], 2), 2); + var df = new NumericalDerivative(); + var x0 = new double[] { 2, 2 }; + var parameterindex = new int[] { 0, 0 }; + Assert.AreEqual(1602, df.EvaluatePartialDerivative(f, x0, 0, 1), 1e-6); + Assert.AreEqual(4002, df.EvaluateMixedPartialDerivative(f, x0, parameterindex, 2)); + } + + [Test] + public void RosenbrockFunctionMixedDerivativeTwoVariableSecondOrderTest() + { + Func f = x => Math.Pow(1 - x[0], 2) + 100 * Math.Pow(x[1] - Math.Pow(x[0], 2), 2); + var df = new NumericalDerivative(); + var x0 = new double[] { 2, 2 }; + var parameterIndex = new[] { 0, 1 }; + + Assert.AreEqual(-800, df.EvaluateMixedPartialDerivative(f, x0, parameterIndex, 2)); + } + + [Test] + public void VectorFunction1PartialDerivativeTest() + { + Func[] f = + { + (x) => Math.Pow(x[0],2) - 3*x[1], + (x) => x[1]*x[1] + 2*x[0]*x[1] + }; + + var x0 = new double[] { 2, 2 }; + var g = new double[] { 4, 4 }; + + var df = new NumericalDerivative(); + Assert.AreEqual(g, df.EvaluatePartialDerivative(f, x0, 0, 1)); + Assert.AreEqual(new double[] { 2, 0 }, df.EvaluatePartialDerivative(f, x0, 0, 2)); + Assert.AreEqual(new double[] { -3, 8 }, df.EvaluatePartialDerivative(f, x0, 1, 1)); + } + + [Test] + public void VectorFunctionMixedPartialDerivativeTest() + { + Func[] f = + { + (x) => Math.Pow(x[0],2) - 3*x[1], + (x) => x[1]*x[1] + 2*x[0]*x[1] + }; + + var x0 = new double[] { 2, 2 }; + + var df = new NumericalDerivative(); + Assert.AreEqual(new double[] { 0, 2 }, df.EvaluateMixedPartialDerivative(f, x0, new int[] { 0, 1 }, 2)); + Assert.AreEqual(new double[] { 0, 2 }, df.EvaluateMixedPartialDerivative(f, x0, new int[] { 1, 0 }, 2)); + + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 1836e16d..1ac8eb41 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -77,6 +77,8 @@ + + From 36050ed83d192f90ef83cd5416ed9675af28e7ee Mon Sep 17 00:00:00 2001 From: Hythem Sidky Date: Fri, 9 Jan 2015 21:24:49 -0500 Subject: [PATCH 2/2] Numerical Jacobian and Hessian wrappers --- .../Differentiation/NumericalHessian.cs | 119 ++++++++++++ .../Differentiation/NumericalJacobian.cs | 172 ++++++++++++++++++ src/Numerics/Numerics.csproj | 2 + .../NumericalHessianTests.cs | 76 ++++++++ .../NumericalJacobianTests.cs | 75 ++++++++ src/UnitTests/UnitTests.csproj | 2 + 6 files changed, 446 insertions(+) create mode 100644 src/Numerics/Differentiation/NumericalHessian.cs create mode 100644 src/Numerics/Differentiation/NumericalJacobian.cs create mode 100644 src/UnitTests/DifferentiationTests/NumericalHessianTests.cs create mode 100644 src/UnitTests/DifferentiationTests/NumericalJacobianTests.cs diff --git a/src/Numerics/Differentiation/NumericalHessian.cs b/src/Numerics/Differentiation/NumericalHessian.cs new file mode 100644 index 00000000..6edf6ba9 --- /dev/null +++ b/src/Numerics/Differentiation/NumericalHessian.cs @@ -0,0 +1,119 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// 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. +// + +using System; + +namespace MathNet.Numerics.Differentiation +{ + /// + /// Class for evaluating the Hessian of a smooth continuously differentiable function using finite differences. + /// By default, a central 3-point method is used. + /// + public class NumericalHessian + { + /// + /// Number of function evaluations. + /// + public int FunctionEvaluations + { + get { return _df.Evaluations; } + } + + private readonly NumericalDerivative _df; + + /// + /// Creates a numerical Hessian object with a three point central difference method. + /// + public NumericalHessian() : this(3, 1) { } + + /// + /// Creates a numerical Hessian with a specified differentiation scheme. + /// + /// Number of points for Hessian evaluation. + /// Center point for differentiation. + public NumericalHessian(int points, int center) + { + _df = new NumericalDerivative(points, center); + } + + /// + /// Evaluates the Hessian of the scalar univariate function f at points x. + /// + /// Scalar univariate function handle. + /// Point at which to evaluate Hessian. + /// Hessian tensor. + public double[] Evaluate(Func f, double x) + { + return new[] { _df.EvaluateDerivative(f, x, 2) }; + } + + /// + /// Evaluates the Hessian of a multivariate function f at points x. + /// + /// + /// This method of computing the Hessian is only vaid for Lipschitz continuous functions. + /// The function mirrors the Hessian along the diagonal since d2f/dxdy = d2f/dydx for continuously differentiable functions. + /// + /// Multivariate function handle.> + /// Points at which to evaluate Hessian.> + /// Hessian tensor. + public double[,] Evaluate(Func f, double[] x) + { + var hessian = new double[x.Length, x.Length]; + + // Compute diagonal elements + for (var row = 0; row < x.Length; row++) + { + hessian[row, row] = _df.EvaluatePartialDerivative(f, x, row, 2); + } + + // Compute non-diagonal elements + for (var row = 0; row < x.Length; row++) + { + for (var col = 0; col < row; col++) + { + var mixedPartial = _df.EvaluateMixedPartialDerivative(f, x, new[] { row, col }, 2); + + hessian[row, col] = mixedPartial; + hessian[col, row] = mixedPartial; + } + } + return hessian; + } + + /// + /// Resets the function evaluation counter for the Hessian. + /// + public void ResetFunctionEvaluations() + { + _df.ResetEvaluations(); + } + } +} diff --git a/src/Numerics/Differentiation/NumericalJacobian.cs b/src/Numerics/Differentiation/NumericalJacobian.cs new file mode 100644 index 00000000..ec04debc --- /dev/null +++ b/src/Numerics/Differentiation/NumericalJacobian.cs @@ -0,0 +1,172 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// 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. +// + +using System; + +namespace MathNet.Numerics.Differentiation +{ + /// + /// Class for evaluating the Jacobian of a function using finite differences. + /// By default, a central 3-point method is used. + /// + public class NumericalJacobian + { + /// + /// Number of function evaluations. + /// + public int FunctionEvaluations + { + get { return _df.Evaluations; } + } + + private readonly NumericalDerivative _df; + + /// + /// Creates a numerical Jacobian object with a three point central difference method. + /// + public NumericalJacobian() : this(3, 1) { } + + /// + /// Creates a numerical Jacobian with a specified differentiation scheme. + /// + /// Number of points for Jacobian evaluation. + /// Center point for differentation. + public NumericalJacobian(int points, int center) + { + _df = new NumericalDerivative(points, center); + } + + /// + /// Evaluates the Jacbian of scalar univariate function f at point x. + /// + /// Scalar univariate function handle. + /// Point at which to evaluate Jacobian. + /// Jacobian vector. + public double[] Evaluate(Func f, double x) + { + return new[] { _df.EvaluateDerivative(f, x, 1) }; + } + + /// + /// Evaluates the Jacobian of a multivariate function f at vector x. + /// + /// + /// This function assumes that the length of vector x consistent with the argument count of f. + /// + /// Multivariate function handle. + /// Points at which to evaluate Jacobian. + /// Jacobian vector. + public double[] Evaluate(Func 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; + } + + /// + /// Evaluates the Jacobian of a multivariate function f at vector x given a current function value. + /// + /// + /// 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. + /// + /// Multivariate function handle. + /// Points at which to evaluate Jacobian. + /// Current function value at finite difference center. + /// Jacobian vector. + public double[] Evaluate(Func 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; + } + + /// + /// Evaluates the Jacobian of a multivariate function array f at vector x. + /// + /// Multivariate function array handle. + /// Vector at which to evaluate Jacobian. + /// Jacobian matrix. + public double[,] Evaluate(Func[] 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; + } + + /// + /// Evaluates the Jacobian of a multivariate function array f at vector x given a vector of current function values. + /// + /// + /// 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. + /// + /// Multivariate function array handle. + /// Vector at which to evaluate Jacobian. + /// Vector of current function values. + /// Jacobian matrix. + public double[,] Evaluate(Func[] 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; + } + + /// + /// Resets the function evaluation counter for the Jacobian. + /// + public void ResetFunctionEvaluations() + { + _df.ResetEvaluations(); + } + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 166b4258..2172e7d4 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -82,6 +82,8 @@ + + diff --git a/src/UnitTests/DifferentiationTests/NumericalHessianTests.cs b/src/UnitTests/DifferentiationTests/NumericalHessianTests.cs new file mode 100644 index 00000000..a833a3a3 --- /dev/null +++ b/src/UnitTests/DifferentiationTests/NumericalHessianTests.cs @@ -0,0 +1,76 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// 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. +// + +namespace MathNet.Numerics.UnitTests.DifferentiationTests +{ + using System; + using Differentiation; + using NUnit.Framework; + + [TestFixture, Category("Differentiation")] + class NumericalHessianTests + { + [Test] + public void ScalarFunctonSecondDerivativeTest() + { + // d2/dx2((2x)^2) = 8 + Func f = x => Math.Pow(2 * x, 2); + var J = new NumericalHessian(); + var jeval = J.Evaluate(f, 2); + Assert.AreEqual((double)8, jeval[0], 1e-7); + } + + [Test] + public void OneEquationTwoVarVectorFunctionTest() + { + Func f = x => 4 * Math.Pow(x[0], 2) + 3 * Math.Pow(x[1], 3); + var H = new NumericalHessian(); + var heval = H.Evaluate(f, new double[] { 1, 1 }); + var solution = new double[,] { { 8, 0 }, { 0, 18 } }; + Assert.AreEqual(solution, heval); + Assert.AreEqual(12, H.FunctionEvaluations); + H.ResetFunctionEvaluations(); + Assert.AreEqual(0, H.FunctionEvaluations); + } + + [Test] + public void RosenbrockFunctionHessianTest() + { + Func f = x => Math.Pow((1 - x[0]), 2) + 100 * Math.Pow(x[1] - Math.Pow(x[0], 2), 2); + var H = new NumericalHessian(5, 2); + var heval = H.Evaluate(f, new double[] { 2, 3 }); + var solution = new double[,] { { 3602, -800 }, { -800, 200 } }; + for (int row = 0; row < solution.Rank; row++) + { + for (int col = 0; col < solution.Rank; col++) + { + Assert.AreEqual(solution[row, col], heval[row, col], 1e-7); + } + } + + } + } +} diff --git a/src/UnitTests/DifferentiationTests/NumericalJacobianTests.cs b/src/UnitTests/DifferentiationTests/NumericalJacobianTests.cs new file mode 100644 index 00000000..c9833070 --- /dev/null +++ b/src/UnitTests/DifferentiationTests/NumericalJacobianTests.cs @@ -0,0 +1,75 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// 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. +// + +namespace MathNet.Numerics.UnitTests.DifferentiationTests +{ + using System; + using Differentiation; + using NUnit.Framework; + + [TestFixture, Category("Differentiation")] + class NumericalJacobianTests + { + [Test] + public void OneEquationVectorFunctionTest() + { + Func f = x => Math.Sin(x[0] * x[1]) * Math.Exp(-x[0] / 2) + x[1] / x[0]; + var J = new NumericalJacobian(); + var jeval = J.Evaluate(f, new double[] { 1, 1 }); + Assert.AreEqual(-0.92747906175, jeval[0], 1e-9); + Assert.AreEqual(1.32770991402, jeval[1], 1e-9); + Assert.AreEqual(4, J.FunctionEvaluations); + J.ResetFunctionEvaluations(); + Assert.AreEqual(0, J.FunctionEvaluations); + } + + [Test] + public void ScalarFunctonDerivativeTest() + { + Func f = x => 1 / x; + var J = new NumericalJacobian(); + var jeval = J.Evaluate(f, 2); + Assert.AreEqual((double)-1 / 4, jeval[0], 1e-7); + } + + [Test] + public void VectorFunctionJacobianTest() + { + Func[] f = + { + (x) => Math.Pow(x[0], 2)*x[1], + (x) => 5*x[0] + Math.Sin(x[1]) + }; + var j = new NumericalJacobian(); + var jeval = j.Evaluate(f, new double[] { 1, 1 }); + Assert.AreEqual(2, jeval[0, 0]); + Assert.AreEqual(1, jeval[0, 1]); + Assert.AreEqual(5, jeval[1, 0]); + Assert.AreEqual(Math.Cos(1), jeval[1, 1], 1e-7); + + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 1ac8eb41..54482cd3 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -79,6 +79,8 @@ + +