Computation assignment- Numerical Integration- Niederreiter sequence

[navneeraj]

Niederreiter sequence is a generalization of the Sobol sequence for any natural number base. Sobol Sequence is defined only for base 2. This is the Wikipedia reference:

[https://en.wikipedia.org/wiki/Sobol_sequence]

This is the relevant excerpt from the Wikipedia page:-

"In his article, Sobol described Πτ-meshes and LPτ sequences, which are (t,m,s)-nets and (t,s)-sequences in base 2 respectively. The terms (t,m,s)-nets and (t,s)-sequences in base b (also called Niederreiter sequences) were coined in 1988 by Harald Niederreiter.[2] The term Sobol sequences was introduced in late English-speaking papers in comparison with Halton, Faure and other low-discrepancy sequences."

We(@rebekahanne, @mattdbrown17,@navneeraj) were also able to find the paper, which has Figure 14.2 of the Numerical Integration problem set. This is the link:

https://cseweb.ucsd.edu/~dstefan/pubs/dalal:2008:low.pdf

I think for both Baker and Niederreiter sequences, we have to use something known as Gray Code to implement them. The reference on the Gray Code is here:

https://en.wikipedia.org/wiki/Gray_code

Someone has implemented the Sobol sequence in Python. The link is here:-

https://github.com/naught101/sobol_seq/blob/master/sobol_seq/sobol_seq.py

I am still trying to understand it and may try to implement it for calculating pi for nth time.:-)

[rickecon]

Here is some help on your equidistributed sequences problem. The third function of this code below is adapted from the code of one of my former students @jmbejara. I have created a module that has the following four functions for generating equidistributed sequences. The first function checks if a number n is a prime number.

def isPrime(n):
    '''
    --------------------------------------------------------------------
    This function returns a boolean indicating whether an integer n is a
    prime number
    --------------------------------------------------------------------
    INPUTS:
    n = scalar, any scalar value

    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None

    OBJECTS CREATED WITHIN FUNCTION:
    i = integer in [2, sqrt(n)]

    FILES CREATED BY THIS FUNCTION: None

    RETURN: boolean
    --------------------------------------------------------------------
    '''
    for i in range(2, int(np.sqrt(n) + 1)):
        if n % i == 0:
            return False

    return True

The second function creates an ascending list of consecutive prime numbers of length N. And it calls the first function isPrime().

def primes_ascend(N, min_val=2):
    '''
    --------------------------------------------------------------------
    This function generates an ordered sequence of N consecutive prime
    numbers, the smallest of which is greater than or equal to 1 using
    the Sieve of Eratosthenes algorithm.
    (https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes)
    --------------------------------------------------------------------
    INPUTS:
    N       = integer, number of elements in sequence of consecutive
              prime numbers
    min_val = scalar >= 2, the smallest prime number in the consecutive
              sequence must be greater-than-or-equal-to this value

    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:
        isPrime()

    OBJECTS CREATED WITHIN FUNCTION:
    primes_vec     = (N,) vector, consecutive prime numbers greater than
                     min_val
    MinIsEven      = boolean, =True if min_val is even, =False otherwise
    MinIsGrtrThn2  = boolean, =True if min_val is
                     greater-than-or-equal-to 2, =False otherwise
    curr_prime_ind = integer >= 0, running count of prime numbers found

    FILES CREATED BY THIS FUNCTION: None

    RETURN: primes_vec
    --------------------------------------------------------------------
    '''
    primes_vec = np.zeros(N, dtype=int)
    MinIsEven = 1 - min_val % 2
    MinIsGrtrThn2 = min_val > 2
    curr_prime_ind = 0
    if not MinIsGrtrThn2:
        i = 2
        curr_prime_ind += 1
        primes_vec[0] = i
    i = min(3, min_val + (MinIsEven * 1))
    while curr_prime_ind < N:
        if isPrime(i):
            curr_prime_ind += 1
            primes_vec[curr_prime_ind - 1] = i
        i += 2

    return primes_vec

The third function creates the nth element of an equidistributed sequence. Notice that I have removed the "weyl", "haber", and "baker" sequence element code because I want you to write that one yourselves. But I left in the "niederreiter" sequence code.

def equidistr_nth(n, d, seq_type='niederreiter'):
    '''
    --------------------------------------------------------------------
    This function returns the nth element of a d-dimensional
    equidistributed sequence. Sequence types available to this function
    are Weyl, Haber, Niederreiter, and Baker. This function follows the
    exposition in Judd (1998, chap. 9).
    --------------------------------------------------------------------
    INPUTS:
    n        = integer >= 1, index of nth value of equidistributed
               sequence where n=1 is the first element
    d        = integer >= 1, number of dimensions in each point of the
               sequence
    seq_type = string, sequence type: "weyl", "haber", "niederreiter",
               or "baker"

    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:
        primes_ascend()

    OBJECTS CREATED WITHIN FUNCTION:
    seq_elem  = (d,) vector, coordinates of element of eq'distr seq
    prime_vec = (d,) vector, consecutive prime numbers
    x         = (d,) vector, whole and fractional part of
                equidistributed sequence coordinates
    x_floor   = (d,) vector, whole part of equidistributed sequence
                coordinates
    ar        = (d,) vector, ascending integer values from 1 to d
    error_str = string, error message for unrecognized name error

    FILES CREATED BY THIS FUNCTION: None

    RETURN: seq_elem
    --------------------------------------------------------------------
    '''
    seq_elem = np.empty(d)
    prime_vec = primes_ascend(d)

    if seq_type == 'weyl':
        # I have removed "weyl" because I want you to do it yourself.
        seq_elem = np.ones(d)

    elif seq_type == 'haber':
        # I have removed "haber" because I want you to do it yourself.
        seq_elem = np.ones(d)

    elif seq_type == 'niederreiter':
        ar = np.arange(1, d + 1)
        x = n * (2 ** (ar / (d + 1)))
        x_floor = np.floor(x)
        seq_elem = x - x_floor

    elif seq_type == 'baker':
        # I have removed "baker" because I want you to do it yourself.
        seq_elem = np.ones(d)

    else:
        error_str = ('Equidistributed sequence name in seq_type not ' +
                     'recognized.')
        raise NameError(error_str)

    return seq_elem

And the fourth function creates a sequence of N consecutive elements of an equidistributed sequence n=1,2,...N.

def equidistr_seq(N, d, seq_type):
    '''
    --------------------------------------------------------------------
    This function returns a vector of N d-dimensional elements of an
    equidistributed sequence. Sequence types available to this function
    are Weyl, Haber, Niederreiter, and Baker.
    --------------------------------------------------------------------
    INPUTS:
    N        = integer >= 1, number of consecutive elements of
               equidistributed sequence to compute, n=1,2,...N
    d        = integer >= 1, number of dimensions in each point of the
               sequence
    seq_type = string, sequence type: "weyl", "haber", "niederreiter",
               or "baker"

    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:
        equidistr_nth()

    OBJECTS CREATED WITHIN FUNCTION:
    eq_seq = (N, d) matrix, N elements of d-dimensional equidistributed
             sequence
    n_val  = integer in [1, N], index of nth value of equidistributed
             sequence

    FILES CREATED BY THIS FUNCTION: None

    RETURN: eq_seq
    --------------------------------------------------------------------
    '''
    eq_seq = np.zeros((N, d))
    for n_val in range(1, N + 1):
        eq_seq[n_val - 1, :] = equidistr_nth(n_val, d, seq_type)

    return eq_seq

You can test this out by writing these functions and saving them in a script named something like equidistributed.py. Because this script will only have import commands (all you need is import numpy as np) and function definitions, this script is called a module. You can import a module in the following way.

import equidistributed as eqd

Then test out your module by doing commands like the following.

nied_500 = eqd.equidistr_seq(500, 2, 'niederreiter')
weyl_1000 = eqd.equidistr_seq(1000, 2, 'weyl') 

You can also replicate the plots in your problem set by the following commands.

import matplotlib.pyplot as plt

plt.scatter(nied_500[:,0], nied_500[:, 1], s=2)
plt.show()

plt.scatter(weyl_1000[:,0], weyl_1000[:, 1], s=1)
plt.show()

Enjoy!

@rebekahanne, @mattdbrown17, @navneeraj, @zeshunzong

[rickecon]

By the way, that is what clean, well-documented, PEP8-compliant code looks like.