VB.NET Quasi Random Example

Imports System
Imports System.Collections

Imports CenterSpace.NMath.Core

Namespace CenterSpace.NMath.Core.Examples.VisualBasic

  ' <summary>
  ' A .NET example in C# showing how to generate a sequence of quasirandom points.
  ' </summary>
  Module QuasiRandomExample

    Sub Main()

    ' A quasirandom sequence is a sequence of n-tuples that fills n-space
    ' more uniformly than uncorrelated random points. NMath provides both
    ' Sobol and Niederreiter quasirandom number generators.

    ' Create a Niederreiter quasirandom number generator object with dimension 1
    ' to generate quasirandom numbers (not vectors).
      Dim Dimensions As Integer = 1
      Dim nqrg As New NiederreiterQuasiRandomGenerator(Dimensions)

    ' You can fill an existing matrix or array. (The points are the columns of the matrix,
    ' so the number of rows in the given matrix must be equal to the Dimension of the
    ' quasirandom number generator.)

    ' Here we create and fill a 1 x numPts matrix with quasirandom numbers from a uniform
    ' (0,1) distribution. The resulting numbers will occupy the first row of A.
      Dim NumPts As Integer = 100
      Dim A As New DoubleMatrix(nqrg.Dimension, NumPts)
      nqrg.Fill(A)

    ' Create a histogram to verify that the quasirandom numbers are uniformly
    ' distributed in the interval (0,1). We create 10 bins expecting that each
    ' bin will contain approximately numPts/10 points.
      Dim H As New Histogram(10, 0, 1)
      h.AddData(A.Row(0))
      Console.WriteLine("{0} uniform (0,1) quasirandom points:", numPts)
      Console.WriteLine(h.StemLeaf())

    ' Compare the above results with those from a pseudo-random generator. One expects
    ' a not-so uniform distribution in space.
      Dim Seed As Integer = &H345
      Dim Stream As New RandomNumberStream(Seed, RandomNumberStream.BasicRandGenType.MersenneTwister)
      Dim UnifDist As New DoubleRandomUniformDistribution()
      Dim pseudoRandom As New DoubleVector(NumPts, Stream, UnifDist)
      h.Reset()
      H.AddData(pseudoRandom)
      Console.WriteLine()
      Console.WriteLine("{0} uniform (0,1) pseudo-random points:", NumPts)
      Console.WriteLine(h.StemLeaf())

    ' Quasi-random numbers are useful for evaluating the integral of a function in the 
    ' unit n-dimensional cube - essentially calculating the average of the function
    ' at a set of randomly selected points. 

    ' In the following example we use a Sobol generator to approximate the integral
    ' of a function F over the 6-dimensional unit cube. The function is defined by
    ' F(x) = 1*cos(1*x1) * 2*cos(2*x2) * 3*cos(3*x3) *...* 6*cos(6*x6)
    ' where x = (x1, x2,..., x6) is a point in 6-dimensional space.
      Dimensions = 6

    ' In this example we provide our own primitive polynomials for initializing
    ' the generator. Primitive polynomials have coefficients in {0, 1} and are
    ' specified by BitArray's containing the polynomial coefficients with the
    ' leading coefficient at index 0.

    ' We can supply either dim polynomials, or dim - 1 polynomials. If we specify
    ' dim - 1 polynomials the primitive polynomial for the first dimension will 
    ' be initialized with a default value.
      Dim primitivePolys(Dimensions - 1) As BitArray
      primitivePolys(0) = New BitArray(New Boolean() {True, True}) ' x + 1
      primitivePolys(1) = New BitArray(New Boolean() {True, True, True}) ' x^2 + x + 1
      primitivePolys(2) = New BitArray(New Boolean() {True, False, True, True}) ' x^3 + x + 1
      primitivePolys(3) = New BitArray(New Boolean() {True, True, False, True}) ' x^3 + x^2 + 1
      primitivePolys(4) = New BitArray(New Boolean() {True, False, False, True, True}) ' x^4 + x + 1
      primitivePolys(5) = New BitArray(New Boolean() {True, True, False, False, True}) ' x^4 + x^3 + 1

      Dim Sobol As New SobolQuasiRandomGenerator(Dimensions, primitivePolys)
      Dim NumPoints As Integer = 180000
      Dim Points As New DoubleMatrix(Dimensions, numPoints)
      sobol.Fill(points)
      Dim Sum As Double = 0

      For I = 0 To NumPoints - 1
          Sum += F(Points.Col(I))
      Next

      Sum /= NumPoints
      Dim actualIntegralValue As Double = -0.022
      Console.WriteLine("Actual integral value = " & actualIntegralValue)
      Console.WriteLine()
      Console.WriteLine("Monte-Carlo approximated integral value = " & Sum)

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

    End Sub

    ' <summary>
    ' F(x) = cos(x(0))*2cos(2x(1))*...*6cos(6*x(5))
    ' </summary>
    ' <param name="x"></param>
    ' <returns></returns>
    Function F(ByVal X As DoubleVector) As Double

      If (X.Length <> 6) Then
        Throw New InvalidArgumentException("x must have length 6!")
      End If

      Dim Y As Double = 1

      For i = 1 To X.Length
        Y *= i * Math.Cos(i * X(i - 1))
      Next
      Return Y
    End Function

  End Module
End Namespace