#! /usr/bin/env python # def i4_to_van_der_corput ( seed, base ): #*****************************************************************************80 # ## I4_TO_VAN_DER_CORPUT computes an element of a van der Corput sequence. # # Discussion: # # The van der Corput sequence is often used to generate a "subrandom" # sequence of points which have a better covering property # than pseudorandom points. # # The van der Corput sequence generates a sequence of points in [0,1] # which (theoretically) never repeats. Except for SEED = 0, the # elements of the van der Corput sequence are strictly between 0 and 1. # # The van der Corput sequence writes an integer in a given base B, # and then its digits are "reflected" about the decimal point. # This maps the numbers from 1 to N into a set of numbers in [0,1], # which are especially nicely distributed if N is one less # than a power of the base. # # Hammersley suggested generating a set of N nicely distributed # points in two dimensions by setting the first component of the # Ith point to I/N, and the second to the van der Corput # value of I in base 2. # # Halton suggested that in many cases, you might not know the number # of points you were generating, so Hammersley's formulation was # not ideal. Instead, he suggested that to generated a nicely # distributed sequence of points in M dimensions, you simply # choose the first M primes, P(1:M), and then for the J-th component of # the I-th point in the sequence, you compute the van der Corput # value of I in base P(J). # # Thus, to generate a Halton sequence in a 2 dimensional space, # it is typical practice to generate a pair of van der Corput sequences, # the first with prime base 2, the second with prime base 3. # Similarly, by using the first K primes, a suitable sequence # in K-dimensional space can be generated. # # The generation is quite simple. Given an integer SEED, the expansion # of SEED in base BASE is generated. Then, essentially, the result R # is generated by writing a decimal point followed by the digits of # the expansion of SEED, in reverse order. This decimal value is actually # still in base BASE, so it must be properly interpreted to generate # a usable value. # # Example: # # BASE = 2 # # SEED SEED van der Corput # decimal binary binary decimal # ------- ------ ------ ------- # 0 = 0 => .0 = 0.0 # 1 = 1 => .1 = 0.5 # 2 = 10 => .01 = 0.25 # 3 = 11 => .11 = 0.75 # 4 = 100 => .001 = 0.125 # 5 = 101 => .101 = 0.625 # 6 = 110 => .011 = 0.375 # 7 = 111 => .111 = 0.875 # 8 = 1000 => .0001 = 0.0625 # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 12 May 2015 # # Author: # # John Burkardt # # Reference: # # J H Halton, # On the efficiency of certain quasi-random sequences of points # in evaluating multi-dimensional integrals, # Numerische Mathematik, # Volume 2, pages 84-90, 1960. # # J M Hammersley, # Monte Carlo methods for solving multivariable problems, # Proceedings of the New York Academy of Science, # Volume 86, pages 844-874, 1960. # # J G van der Corput, # Verteilungsfunktionen I & II, # Nederl. Akad. Wetensch. Proc., # Volume 38, 1935, pages 813-820, pages 1058-1066. # # Parameters: # # Input, integer SEED, the seed or index of the desired element. # SEED should be nonnegative. Only the integer part of SEED is used. # SEED = 0 is allowed, and returns R = 0. # # Input, integer BASE, the van der Corput base, which is typically # a prime number. Only the integer part of BASE is used. # BASE must be greater than 1. # # Output, real R, the SEED-th element of the van der Corput sequence # for base BASE. # from sys import exit r = 0.0 # # Ensure that BASE is an integer, and acceptable. # base = int ( base ) if ( base <= 1 ): print ( '' ) print ( 'I4_TO_VAN_DER_CORPUT - Fatal error!' ) print ( ' The input base BASE is <= 1!' ) print ( ' BASE = %d' % ( base ) ) exit ( 'I4_TO_VAN_DER_CORPUT - Fatal error!' ) # # Ensure that SEED is an integer, and acceptable. # seed = int ( seed ) if ( seed < 0 ): print ( '' ) print ( 'I4_TO_VAN_DER_CORPUT - Fatal error!' ) print ( ' The input base SEED is < 0!' ) print ( ' SEED = %d' % ( seed ) ) exit ( 'I4_TO_VAN_DER_CORPUT - Fatal error!' ); # # Carry out the computation. # base_inv = 1.0 / float ( base ) while ( seed != 0 ): digit = ( seed % base ) r = r + float ( digit ) * base_inv base_inv = base_inv / float ( base ) seed = ( seed // base ) return r def i4_to_van_der_corput_test ( ): #*****************************************************************************80 # ## I4_TO_VAN_DER_CORPUT_TEST tests I4_TO_VAN_DER_CORPUT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 12 May 2015 # # Author: # # John Burkardt # import numpy as np import platform from prime import prime n_prime = 5 n_test = 10 print ( '' ) print ( 'I4_TO_VAN_DER_CORPUT_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' I4_TO_VAN_DER_CORPUT computes the elements' ) print ( ' of a van der Corput sequence.' ) print ( ' The sequence depends on the prime numbers used' ) print ( ' as a base.' ) print ( '' ) print ( ' Bases:' ) print ( '' ) print ( '' ) for j in range ( 1, n_prime + 1 ): print ( ' %6d' % prime ( j ) ), print ( '' ) print ( '' ) h = np.zeros ( n_prime ) for i in range ( 0, n_test ): for j in range ( 1, n_prime + 1 ): jm1 = j - 1 p = prime ( j ) h[jm1] = i4_to_van_der_corput ( i, p ) print ( ' %2d' % ( i ) ), for j in range ( 1, n_prime + 1 ): jm1 = j - 1 print ( ' %12f' % ( h[jm1] ) ), print ( '' ) # # Terminate. # print ( '' ) print ( 'I4_TO_VAN_DER_CORPUT_TEST' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) i4_to_van_der_corput_test ( ) timestamp ( )