Add non-linear least squares MKL wrapper

Add non-linear least squares MKL wrapper and suggested form for
Optimization routines.

src/Numerics/Numerics.csproj

@@ -87,6 +87,9 @@
     <Reference Include="System.Xml" />
   </ItemGroup>
   <ItemGroup>
+    <Compile Include="Optimization\BrentMinimizer.cs" />
+    <Compile Include="Optimization\NonLinearLeastSquaresMinimizer.cs" />
+    <Compile Include="Optimization\PowellMinimizer.cs" />
     <Compile Include="Precision.Comparison.cs" />
     <Compile Include="Precision.Equality.cs" />
     <Compile Include="Distributions\Bernoulli.cs" />
@@ -178,6 +181,9 @@
     <Compile Include="LinearAlgebra\Vector.BCL.cs" />
     <Compile Include="LinearAlgebra\Vector.Operators.cs" />
     <Compile Include="NonConvergenceException.cs" />
+    <Compile Include="Providers\Optimization\IOptimizationProvider.cs" />
+    <Compile Include="Providers\Optimization\Mkl\MklOptimizationProvider.cs" />
+    <Compile Include="Providers\Optimization\Mkl\SafeNativeMethods.cs" />
     <Compile Include="Random\RandomSeed.cs" />
     <Compile Include="RootFinding\Broyden.cs" />
     <Compile Include="RootFinding\NewtonRaphson.cs" />

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
+    {
+    }
+}

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
+{   
+    /// <summary>
+    /// 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.  
+    /// </summary>
+    public class NonLinearLeastSquaresMinimizer
+    {
+        /// <summary>
+        /// Criterion0: Δ &lt; eps(0) (trust region solvers only)
+        /// Criterion1: ||F(x)||2 &lt; eps(1)
+        /// Criterion2: The Jacobian matrix is singular.||J(x)(1:m,j)||2 &lt; eps(2), j = 1, ..., n
+        /// Criterion3: ||s||2 &lt; eps(3)
+        /// Criterion4: ||F(x)||2 - ||F(x) - J(x)s||2 &lt; eps(4)
+        /// </summary>
+        public enum ConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, SingularJacobian, Error };
+        
+        public class Result
+        {
+            public int NumberOfIterations;
+
+            public ConvergenceType ConvergenceType;
+        }
+        
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <param name="x"></param>
+        /// <param name="y"></param>
+        /// <param name="f"></param>
+        /// <param name="pStart"></param>
+        /// <param name="jacobian">jac_j(x, p) = df / dp_j</param>
+        /// <returns></returns>
+        public static double[] CurveFit(double[] x, double[] y, Func<double, double[], double> f,
+            double[] pStart, Func<double, double[], double[]> 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;
+        }
+    }
+}

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
+    {
+
+    }
+}

src/Numerics/Providers/Optimization/IOptimizationProvider.cs

@@ -0,0 +1,63 @@
+// <copyright file="IOptimizationProvider.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+using MathNet.Numerics.Optimization;
+
+namespace MathNet.Numerics.Providers.Optimization
+{
+    /// <summary>
+    /// Function specifying the model. This takes in model parameters, calculates residuals
+    /// and updates the residuals array with these.
+    /// </summary>
+    /// <param name="parameters">The model parameters. The function must not change these.</param>
+    /// <param name="r">The residuals to be updated. The existing array should be updated.</param>
+    /// <returns></returns>
+    public delegate void LeastSquaresForwardModel(double[] p, double[] r);
+
+    /// <summary>
+    /// 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
+    /// </summary>
+    /// <param name="p">THe model parameters. The function must not change these.</param>
+    /// <param name="jacobian">THe Jacobian matrix in column-major format.</param>
+    public delegate void Jacobian(double[] p, double[] jacobian);
+    
+    /// <summary>
+    /// Interface to linear algebra algorithms that work off 1-D arrays.
+    /// </summary>
+    /// <typeparam name="T">Supported data types are Double, Single, Complex, and Complex32.</typeparam>
+    public interface IOptimizationProvider<T>
+        where T : struct
+    {
+        NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize(
+            int residualsLength, T[] initialGuess, LeastSquaresForwardModel function, 
+            out T[] parameters, Jacobian jacobianFunction = null);
+    }
+}

src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs

@@ -0,0 +1,206 @@
+// <copyright file="MklOptimizationProvider.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+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<double>
+    {
+        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

src/Numerics/Providers/Optimization/Mkl/SafeNativeMethods.cs

@@ -0,0 +1,84 @@
+// <copyright file="SafeNativeMethods.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+#if NATIVEMKL
+
+using System.Numerics;
+using System.Runtime.InteropServices;
+using System.Security;
+using System;
+
+namespace MathNet.Numerics.Providers.Optimization.Mkl
+{
+    /// <summary>
+    /// P/Invoke methods to the native math libraries.
+    /// </summary>
+    [SuppressUnmanagedCodeSecurity]
+    [SecurityCritical]
+    internal static class SafeNativeMethods
+    {
+        /// <summary>
+        /// Name of the native DLL.
+        /// </summary>
+        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

src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs

@@ -0,0 +1,63 @@
+// <copyright file="NonLinearLeastSquaresTest.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+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<double, double[], double> function = (x, p) => p[0] * (1 - Math.Exp(-p[1] * x));
+            Func<double, double[], double[]> 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);
+        }
+
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -353,6 +353,7 @@
     <Compile Include="NumberTheoryTests\GcdRelatedTest.cs" />
     <Compile Include="NumberTheoryTests\GcdRelatedTestBigInteger.cs" />
     <Compile Include="NumberTheoryTests\IntegerTheoryTest.cs" />
+    <Compile Include="OptimizationTests\NonLinearLeastSquaresTest.cs" />
     <Compile Include="RootFindingTests\BisectionTest.cs" />
     <Compile Include="PermutationTest.cs" />
     <Compile Include="PrecisionTest.cs" />