VB Quasi Random Example

Imports System
Imports System.Collections

Imports CenterSpace.NMath.Core

Namespace CenterSpace.NMath.Examples.VisualBasic

  <summary>
  A .NET example in Visual Basic 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()

      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 BitArrays 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