diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 00054cff..32c01fa0 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -87,6 +87,9 @@ + + + @@ -178,6 +181,9 @@ + + + diff --git a/src/Numerics/Optimization/BrentMinimizer.cs b/src/Numerics/Optimization/BrentMinimizer.cs new file mode 100644 index 00000000..6b6de901 --- /dev/null +++ b/src/Numerics/Optimization/BrentMinimizer.cs @@ -0,0 +1,11 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public class BrentMinimizer + { + } +} diff --git a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs new file mode 100644 index 00000000..533869b5 --- /dev/null +++ b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs @@ -0,0 +1,71 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using MathNet.Numerics.Providers.Optimization; +using MathNet.Numerics.Providers.Optimization.Mkl; + +namespace MathNet.Numerics.Optimization +{ + /// + /// This class is a special function minimizer that minimizes functions of the form + /// f(p) = |r(p)|^2 where r is a vector of residuals and p is a vector of model parameters. + /// + public class NonLinearLeastSquaresMinimizer + { + /// + /// Criterion0: Δ < eps(0) (trust region solvers only) + /// Criterion1: ||F(x)||2 < eps(1) + /// Criterion2: The Jacobian matrix is singular.||J(x)(1:m,j)||2 < eps(2), j = 1, ..., n + /// Criterion3: ||s||2 < eps(3) + /// Criterion4: ||F(x)||2 - ||F(x) - J(x)s||2 < eps(4) + /// + public enum ConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, SingularJacobian, Error }; + + public class Result + { + public int NumberOfIterations; + + public ConvergenceType ConvergenceType; + } + + /// + /// Non-Linear Least-Squares fitting the points (x,y) to a specified function of y : x -> f(x, p), p being a vector of parameters. + /// returning its best fitting parameters p. + /// + /// + /// + /// + /// + /// jac_j(x, p) = df / dp_j + /// + public static double[] CurveFit(double[] x, double[] y, Func f, + double[] pStart, Func jacobian = null) + { + if (x.Length != y.Length) throw new ArgumentException("x and y lengths different"); + var provider = new MklOptimizationProvider(); + LeastSquaresForwardModel function = (p, r) => + { + for (int i = 0; i < r.Length; ++i) + r[i] = y[i] - f(x[i], p); + }; + + // jac is df_i / dp_j + + Jacobian jacobianFunction = null; + if (jacobian != null) jacobianFunction = (p, jac) => + { + for (int i = 0; i < y.Length; ++i) + { + double[] values = jacobian(x[i], p); + for (int j = 0; j < values.Length; ++j) + jac[j * y.Length + i] = -values[j]; + } + }; + + double[] parameters; + Result result = provider.NonLinearLeastSquaresUnboundedMinimize(y.Length, pStart, function, out parameters, jacobianFunction); + return parameters; + } + } +} diff --git a/src/Numerics/Optimization/PowellMinimizer.cs b/src/Numerics/Optimization/PowellMinimizer.cs new file mode 100644 index 00000000..a1377516 --- /dev/null +++ b/src/Numerics/Optimization/PowellMinimizer.cs @@ -0,0 +1,12 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public class PowellSolver + { + + } +} diff --git a/src/Numerics/Providers/Optimization/IOptimizationProvider.cs b/src/Numerics/Providers/Optimization/IOptimizationProvider.cs new file mode 100644 index 00000000..b5463849 --- /dev/null +++ b/src/Numerics/Providers/Optimization/IOptimizationProvider.cs @@ -0,0 +1,63 @@ +// +// 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-2013 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 MathNet.Numerics.Optimization; + +namespace MathNet.Numerics.Providers.Optimization +{ + /// + /// Function specifying the model. This takes in model parameters, calculates residuals + /// and updates the residuals array with these. + /// + /// The model parameters. The function must not change these. + /// The residuals to be updated. The existing array should be updated. + /// + public delegate void LeastSquaresForwardModel(double[] p, double[] r); + + /// + /// Function providing the Jacobian matrix in column-major format for a set of model parameter values. + /// Jacobian is dr_i / dp_i, r being the residuals vector and p the vector of parameters + /// + /// THe model parameters. The function must not change these. + /// THe Jacobian matrix in column-major format. + public delegate void Jacobian(double[] p, double[] jacobian); + + /// + /// Interface to linear algebra algorithms that work off 1-D arrays. + /// + /// Supported data types are Double, Single, Complex, and Complex32. + public interface IOptimizationProvider + where T : struct + { + NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize( + int residualsLength, T[] initialGuess, LeastSquaresForwardModel function, + out T[] parameters, Jacobian jacobianFunction = null); + } +} diff --git a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs new file mode 100644 index 00000000..571ca076 --- /dev/null +++ b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs @@ -0,0 +1,206 @@ +// +// 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-2013 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 MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; +using MathNet.Numerics.Threading; +using System; +using MathNet.Numerics.Optimization; + +#if NATIVEMKL + +namespace MathNet.Numerics.Providers.Optimization.Mkl +{ + public class MklOptimizationProvider : IOptimizationProvider + { + const int TR_SUCCESS = 1501; + + public NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize(int residualsLength, double[] initialGuess, LeastSquaresForwardModel function, out double[] parameters, Jacobian jacobianFunction = null) + { + bool analyticJacobian = jacobianFunction != null; + double[] residuals = new double[residualsLength]; + double[] residualsMinus = new double[residualsLength]; + double[] residualsPlus = new double[residualsLength]; + double[] jacobian = new double[residualsLength * initialGuess.Length]; + parameters = new double[initialGuess.Length]; + + double[] eps = new double[6]; // stop criteria + int i; + for (i = 0; i < 6; i++) + eps[i] = 1e-8; + + for (i = 0; i < initialGuess.Length; i++) + parameters[i] = initialGuess[i]; + + int successful; + + int maxIterations = 1000, maxTrialStepIterations = 100; + + IntPtr solverHandle = IntPtr.Zero; + IntPtr jacobianHandle = IntPtr.Zero; + + int[] info = new int[6]; // for parameter checking + + double initialStepBound = 0.0; + + double jacobianPrecision = 1e-8; + + // zero initial values: + for (i = 0; i < residuals.Length; i++) + residuals[i] = 0.0; + for (i = 0; i < residuals.Length * parameters.Length; i++) + jacobian[i] = 0.0; + + if (SafeNativeMethods.unbound_nonlinearleastsq_init(ref solverHandle, parameters.Length, residualsLength, parameters, eps, maxIterations, maxTrialStepIterations, initialStepBound) != + TR_SUCCESS) + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + + if (SafeNativeMethods.unbound_nonlinearleastsq_check(ref solverHandle, parameters.Length, residualsLength, jacobian, residuals, eps, info) != TR_SUCCESS) + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + else + { + if (info[0] != 0 || // Handle invalid + info[1] != 0 || // Jacobian array not valid + info[2] != 0 || // Parameters array not valid + info[3] != 0) // Eps array not valid + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + } + + if (SafeNativeMethods.jacobi_init(ref jacobianHandle, parameters.Length, residuals.Length, parameters, jacobian, jacobianPrecision) != TR_SUCCESS) + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + + int rciRequest = 0; + successful = 0; + while (successful == 0) + { + if (SafeNativeMethods.unbound_nonlinearleastsq_solve(ref solverHandle, residuals, jacobian, ref rciRequest) != TR_SUCCESS) + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + if (rciRequest == -1 || rciRequest == -2 || rciRequest == -3 || + rciRequest == -4 || rciRequest == -5 || rciRequest == -6) + successful = 1; + if (rciRequest == 1) // recalculate function to update parameters + { + function(parameters, residuals); + } + if (rciRequest == 2) + { + if (analyticJacobian) + jacobianFunction(parameters, jacobian); + else + { + // calculate by central differences: + int rciRequestJacobian = 0; + int jacobianSuccessful = 0; + + // update Jacobian matrix: + while (jacobianSuccessful == 0) + { + if (SafeNativeMethods.jacobi_solve(ref jacobianHandle, residualsPlus, residualsMinus, ref rciRequestJacobian) != TR_SUCCESS) + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + if (rciRequestJacobian == 1) + function(parameters, residualsPlus); + else if (rciRequestJacobian == 2) + function(parameters, residualsMinus); + else if (rciRequestJacobian == 0) + jacobianSuccessful = 1; + } + } + } + } + + int stopCriterionNumber = 0, iterations = 0; + double initialResidual = 0, finalResidual = 0; + if (SafeNativeMethods.unbound_nonlinearleastsq_get(ref solverHandle, ref iterations, ref stopCriterionNumber, ref initialResidual, ref finalResidual) != TR_SUCCESS) + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + + if (SafeNativeMethods.unbound_nonlinearleastsq_delete(ref solverHandle) != TR_SUCCESS) + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + + if (SafeNativeMethods.jacobi_delete(ref jacobianHandle) != TR_SUCCESS) + { + SafeNativeMethods.FreeBuffers(); + return ErrorResult(); + } + + SafeNativeMethods.FreeBuffers(); + + NonLinearLeastSquaresMinimizer.ConvergenceType convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Error; + switch (rciRequest) + { + case -1: + convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.NoneMaxIterationExceeded; break; + case -2: + convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion0; break; + case -3: + convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion1; break; + case -4: + convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.SingularJacobian; break; + case -5: + convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion3; break; + case -6: + convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion4; break; + } + + // no errors, find reason for stopping; + return new NonLinearLeastSquaresMinimizer.Result() { ConvergenceType = convergenceType }; + } + + public static NonLinearLeastSquaresMinimizer.Result ErrorResult() + { + return new NonLinearLeastSquaresMinimizer.Result() { ConvergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Error }; + } + } +} + +#endif \ No newline at end of file diff --git a/src/Numerics/Providers/Optimization/Mkl/SafeNativeMethods.cs b/src/Numerics/Providers/Optimization/Mkl/SafeNativeMethods.cs new file mode 100644 index 00000000..7e55cb60 --- /dev/null +++ b/src/Numerics/Providers/Optimization/Mkl/SafeNativeMethods.cs @@ -0,0 +1,84 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009-2013 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. +// + +#if NATIVEMKL + +using System.Numerics; +using System.Runtime.InteropServices; +using System.Security; +using System; + +namespace MathNet.Numerics.Providers.Optimization.Mkl +{ + /// + /// P/Invoke methods to the native math libraries. + /// + [SuppressUnmanagedCodeSecurity] + [SecurityCritical] + internal static class SafeNativeMethods + { + /// + /// Name of the native DLL. + /// + const string DllName = "MathNet.Numerics.MKL.dll"; + + #region Non-Linear Least Squares Unbounded + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int unbound_nonlinearleastsq_init(ref IntPtr handle, int n, int m, double[] x, double[] eps, int iter1, int iter2, double rs); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int unbound_nonlinearleastsq_check(ref IntPtr handle, int n, int m, double[] fjac, double[] fvec, double[] eps, int[] info); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int unbound_nonlinearleastsq_solve(ref IntPtr handle, double[] fvec, double[] fjac, ref int RCI_Request); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int unbound_nonlinearleastsq_get(ref IntPtr handle, ref int iter, ref int st_cr, ref double r1, ref double r2); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int unbound_nonlinearleastsq_delete(ref IntPtr handle); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int jacobi_init(ref IntPtr handle, int n, int m, double[] x, double[] fjac, double eps); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int jacobi_solve(ref IntPtr handle, double[] f1, double[] f2, ref int RCI_Request); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int jacobi_delete(ref IntPtr handle); + + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int FreeBuffers(); + + #endregion + + } +} + +#endif diff --git a/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs new file mode 100644 index 00000000..17cf12d7 --- /dev/null +++ b/src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs @@ -0,0 +1,63 @@ +// +// 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-2013 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.Optimization; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class NonLinearLeastSquaresTest + { + [Test] + public void CurveFit() + { + // y = b1*(1-exp[-b2*x]) + e + var xin = new double[] { 1, 2, 3, 5, 7, 10 }; + var yin = new double[] { 109, 149, 149, 191, 213, 224 }; + var popt = NonLinearLeastSquaresMinimizer.CurveFit(xin, yin, (x, p) => p[0] * (1 - Math.Exp(-p[1] * x)), new double[] { 1, 1 }); + + Func function = (x, p) => p[0] * (1 - Math.Exp(-p[1] * x)); + Func jacobian = (x, p) => new double[] { + 1 - Math.Exp(-p[1] * x), + p[0] * x * Math.Exp(-p[1] * x) }; + + popt = NonLinearLeastSquaresMinimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }, jacobian); // 100, 0.75 + + double[] expected = new double[] { 2.1380940889E+02, 5.4723748542E-01 }; + + double residual = 0; + for (int i = 0; i < yin.Length; ++i) residual += (yin[i] - function(xin[i], popt)) * (yin[i] - function(xin[i], popt)); + //Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001); + } + + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 271014d2..0776c7da 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -353,6 +353,7 @@ +