C# Quasi Random Example

using System;
using System.Collections;

using CenterSpace.NMath.Core;

namespace CenterSpace.NMath.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# showing how to generate a sequence of quasirandom points.
  /// </summary>
  class QuasiRandomExample
  {
    static void Main( string[] args )
    {
      // 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).
      int dim = 1;
      var nqrg = new NiederreiterQuasiRandomGenerator( dim );

      // 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.
      int numPts = 100;
      var A = 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.
      var h = 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.
      int seed = 0x345;
      var stream = new RandomNumberStream( seed, RandomNumberStream.BasicRandGenType.MersenneTwister );
      var unifDist = new DoubleRandomUniformDistribution();
      var pseudoRandom = new DoubleVector( numPts, stream, unifDist );
      h.Reset();
      h.AddData( pseudoRandom );
      Console.WriteLine( "\n{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.
      dim = 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.
      var primitivePolys = new BitArray[dim];
      primitivePolys[0] = new BitArray( new bool[] { true, true } ); // x + 1
      primitivePolys[1] = new BitArray( new bool[] { true, true, true } ); // x^2 + x + 1
      primitivePolys[2] = new BitArray( new bool[] { true, false, true, true } ); // x^3 + x + 1
      primitivePolys[3] = new BitArray( new bool[] { true, true, false, true } ); // x^3 + x^2 + 1
      primitivePolys[4] = new BitArray( new bool[] { true, false, false, true, true } ); // x^4 + x + 1
      primitivePolys[5] = new BitArray( new bool[] { true, true, false, false, true } ); // x^4 + x^3 + 1

      var sobol = new SobolQuasiRandomGenerator( dim, primitivePolys );
      int numPoints = 180000;
      var points = new DoubleMatrix( dim, numPoints );
      sobol.Fill( points );
      double sum = 0;
      for ( int i = 0; i < numPoints; i++ )
      {
        sum += F( points.Col( i ) );
      }
      sum /= numPoints;
      double actualIntegralValue = -0.022;
      Console.WriteLine( "Actual integral value = " + actualIntegralValue );
      Console.WriteLine( "\nMonte-Carlo approximated integral value = " + sum );

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

    /// <summary>
    /// F(x) = cos(x[0])*2cos(2x[1])*...*6cos(6*x[5])
    /// </summary>
    /// <param name="x"></param>
    /// <returns></returns>
    static double F( DoubleVector x )
    {
      if ( x.Length != 6 )
        throw new InvalidArgumentException( "x must have length 6!" );

      double y = 1;
      for ( int i = 1; i <= x.Length; i++ )
      {
        y *= i * Math.Cos( i * x[i - 1] );
      }
      return y;
    }

  } // class

} // namespace