VB.NET Multi Variable Curve Fitting Example

Imports System

Imports CenterSpace.NMath.Core
Imports CenterSpace.NMath.Analysis

Namespace CenterSpace.NMath.Analysis.Examples.VisualBasic

  Module MultiVariableCurveFittingExample

    ' The MultiVariableFunctionFitter&ltT&gt Needs a parameterized function
    ' and a set of data points. One way to specify the parameterized function,
    ' and optionally it's gradient with respect to the parameters, is to
    ' implement an instance of the abstract class DoubleParameterizedFunctional.
    ' You must overwrite the Evaluate() method which computes and returns the
    ' parameterized function value at a specified set of parameters and 
    ' point. It is optional to overwrite the GradientWithRespectToParams() method.
    ' If you do not overwrite it a numerical approximation using fintite differences
    ' will be used to approximate the gradient if it is needed.
    ' 
    ' Here the parameterized function we are defining is a real valued function
    ' of two variables, x0 and x1, and three parameters, p0, p1, and p2, defined
    ' by the formula:
    ' 
    ' p0*x1*x0^2 + p1*sin(x0) + p2*x1^3
    ' 
    Private Class ParameterizedFunction
      Inherits DoubleParameterizedFunctional

      ' Creates an instance of our parameterized function. We must 
      ' initialize the base class with the dimension of our functions
      ' domain. Since our function is a function of two variables 
      ' we initialize the base class with 2.

      Public Sub New()
        MyBase.New(2)
      End Sub

      ' Override the abstract evaluate function.
      ' <param name="parameters">The parameter values.</param>
      ' <param name="x">The point to evaluate at.</param>
      ' <returns>The value of the parameterized function at the given
      ' point and parameters.</returns>
      Public Overrides Function Evaluate(ByVal Parameters As DoubleVector, ByVal X As DoubleVector) As Double

        Return Parameters(0) * X(1) * Math.Pow(X(0), 2.0) + Parameters(1) * Math.Sin(X(0)) + Parameters(2) * Math.Pow(X(1), 3.0)

      End Function

      ' Since the gradient of our function is rather easy to derive, we will
      ' override the GradientWithRespectToParams() function. Remember, this is
      ' the vector of partials with respect to the parameters, NOT the variables.

      ' <param name="parameters">Evaluate the gradient at these parameter values.</param>
      ' <param name="x">Evaluate the gradient at this point.</param>
      ' <param name="grad">Place the value of the gradient in this vector.</param>
      ' <remarks>Note how this function does not return the gradient as a new
      ' vector, but places the gradient value in a vector supplied by the 
      ' calling routine. This is for optimization purposes. The curve fitter uses 
      ' a optimization algorithm that will most likely be iterative, and thus may 
      ' need to evaluate the gradient many times. Having the vector 
      ' passed in to the routine allows the calling code to allocate space for the 
      ' gradient once and reuse it on successive calls, thus avoiding the potential 
      ' of allocating a large number of small objects on the managed heap.</remarks>
      Public Overrides Sub GradientWithRespectToParams(ByVal Parameters As DoubleVector, ByVal X As DoubleVector, ByRef Grad As DoubleVector)
        Grad(0) = X(0) * X(0) * X(1)
        Grad(1) = Math.Sin(X(0))
        Grad(2) = Math.Pow(X(1), 3)
      End Sub

    End Class

    ' A .NET example in C# showing how to fit a generalized multivariable function to a set 
    ' of points.
    ' Uses the trust-region algorithm.
    Sub Main()

      Console.WriteLine()

      ' Class MultiVariableFunctionFitter fits a parameterized multivariable function to a
      ' set of points. In the space of the function parameters, begining at a specified
      ' starting point, the Fit() method finds a minimum (possibly local) in the sum of
      ' the squared residuals with respect to the data. Fit() uses a nonlinear least
      ' squares minimizer specified as a generic.

      ' For example, here is dataset from the Matlab docs, which fits a function
      ' z = f(x, y) to three-dimensional data describing a surface  
      ' http:'www.mathworks.com/support/solutions/data/1-17YMU.html?solution=1-17YMU

      ' Since the domain of the function has two dimensions, we use a two-column matrix to
      ' hold the x,y data.
      Dim XYValues As New DoubleMatrix(10, 2)
      XYValues(Slice.All, 0) = New DoubleVector("3.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 10.0 5.4")
      XYValues(Slice.All, 1) = New DoubleVector("16.5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3")

      Dim ZValues As New DoubleVector("95.09 23.11 60.63 48.59 89.12 76.97 45.68 1.84 82.17 44.47")

      ' Published starting guess in the space of the function parameters.
      Dim Start As New DoubleVector("10 10 10")

      ' Construct a curve fitting object for our function, then perform the fit. We will use the
      ' TrustRegionMinimizer implementation of the non-linear least squares minimizer to find the optimal
      ' set of parameters. 
      Dim F As New ParameterizedFunction()
      Dim Fitter As New MultiVariableFunctionFitter(Of TrustRegionMinimizer)(F)
      Dim Solution As DoubleVector = Fitter.Fit(XYValues, ZValues, Start)

      ' Display the results
      Console.WriteLine("Fit #1")
      Console.WriteLine("Matlab solution: " & New DoubleVector("0.0074 -19.9749 -0.0000").ToString())
      Console.WriteLine("NMath solution: " & Solution.ToString())
      Console.WriteLine("NMath residual: " & Fitter.Minimizer.FinalResidual)
      Console.WriteLine()

      ' The parameterized function used by the fitter may also be specified using a delegate.
      ' here we define a delegate for the same function
      ' p0*x1*x0^2 + p1*sin(x0) + p2*x1^3
      Dim XDimension As Integer = 2 ' The dimension of the domain of f.

      Dim FDelegate As Func(Of DoubleVector, DoubleVector, Double) = AddressOf Foo

      ' The delegate for the parameterized function may be used directly in MultiVariableFunctionFitter
      ' constructors, or may be wrapped by the DoubleVectorParameterizedDelegate, which implements 
      ' DoubleParameterizedFunctional. Here we do the latter.
      ' Note that we do not supply the gradient with respect
      ' to parameters here. The gradient will be computed using a finite difference algorithm if
      ' needed.
      Fitter.Function = New DoubleVectorParameterizedDelegate(fdelegate, XDimension)
      ' Perform the fit and display the results
      Solution = Fitter.Fit(XYValues, ZValues, Start)
      Console.WriteLine("Fit #1 (Repeated without user specified Partial Derivatives)")
      Console.WriteLine("NMath solution: " & Solution.ToString())
      Console.WriteLine("NMath residual: " & Fitter.Minimizer.FinalResidual)
      Console.WriteLine()

      ' Now let's perform the fit again using some random data. First we generate
      ' 50 random x,y points in range (0,10).
      XYValues = New DoubleMatrix(50, 2, New RandGenUniform(0, 10))

      ' The target solution.
      Dim Target As New DoubleVector("1 2 3")

      ' When caculating the z values, we add some noise, so the points
      ' don't lie exactly on the target surface.
      ZValues = New DoubleVector(50)
      Dim Rnd As New RandGenUniform(-1, 1)

      Dim I As Integer
      For I = 0 To ZValues.Length - 1
        ZValues(I) = fdelegate(Target, XYValues.Row(I)) + Rnd.Next()
      Next

      ' Perform the fit and display the results
      Solution = Fitter.Fit(XYValues, ZValues, Start)
      Console.WriteLine("Fit #2")
      Console.WriteLine("Target solution: " & Target.ToString())
      Console.WriteLine("Actual solution: " & Solution.ToString())
      Console.WriteLine("Residual: " & Fitter.Minimizer.FinalResidual)
      Console.WriteLine()

      Console.WriteLine()
      Console.WriteLine("Press Enter Key")
      Console.Read()

    End Sub

    Private Function Foo(ByVal P As DoubleVector, ByVal X As DoubleVector) As Double
      Return P(0) * X(1) * Math.Pow(X(0), 2.0) + P(1) * Math.Sin(X(0)) + P(2) * Math.Pow(X(1), 3.0)
    End Function

  End Module

End Namespace