
SSJ V. 2.6. 

PREV PACKAGE NEXT PACKAGE  FRAMES NO FRAMES 
See:
Description
Interface Summary  

PointSetIterator  Objects of classes that implement this interface are iterators that permit one to enumerate (or observe) the successive points of a point set and the successive coordinates of these points. 
PointSetRandomization  This interface is used to randomize a
PointSet . 
Class Summary  

AntitheticPointSet  This container class provides antithetic points. 
BakerTransformedPointSet  This container class embodies a point set to which a Baker transformation is applied. 
CachedPointSet  This container class caches a point set by precomputing and storing its points locally in an array. 
ContainerPointSet  This acts as a generic base class for all container classes that contain a point set and apply some kind of transformation to the coordinates to define a new point set. 
CycleBasedPointSet  This abstract class provides the basic structures for storing and manipulating a highly uniform point set defined by a set of cycles. 
CycleBasedPointSetBase2  Similar to CycleBasedPointSet , except that the successive
values in the cycles are stored as integers in the range
{0,..., 2^{k} 1}, where
1 <= k <= 31. 
DigitalNet  This class provides the basic structures for storing and manipulating linear digital nets in base b, for an arbitrary base b >= 2. 
DigitalNetBase2  A special case of DigitalNet for the base b = 2. 
DigitalNetBase2FromFile  This class allows us to read the parameters defining a digital net in base 2 either from a file, or from a URL address on the World Wide Web. 
DigitalNetFromFile  This class allows us to read the parameters defining a digital net either from a file, or from a URL address on the World Wide Web. 
DigitalSequence  This abstract class describes methods specific to digital sequences. 
DigitalSequenceBase2  This abstract class describes methods specific to digital sequences in base 2. 
EmptyRandomization  This class implements an empty
PointSetRandomization . 
F2wCycleBasedLFSR  This class creates a point set based upon a linear feedback shift register sequence. 
F2wCycleBasedPolyLCG  This class creates a point set based upon a linear congruential sequence in the finite field F_{2w}[z]/P(z). 
F2wNetLFSR  This class implements a digital net in base 2 starting from a linear feedback shift register generator. 
F2wNetPolyLCG  This class implements a digital net in base 2 starting from a polynomial LCG in F_{2w}[z]/P(z). 
F2wStructure  This class implements methods and fields needed by the classes
F2wNetLFSR ,
F2wNetPolyLCG ,
F2wCycleBasedLFSR and
F2wCycleBasedPolyLCG . 
FaureSequence  This class implements digital nets or digital sequences formed by the first n = b^{k} points of the Faure sequence in base b. 
HaltonSequence  This class implements the sequence of Halton, which is essentially a modification of Hammersley nets for producing an infinite sequence of points having low discrepancy. 
HammersleyPointSet  This class implements Hammersley point sets, which are defined as follows. 
KorobovLattice  This class implements Korobov lattices, which represents the same point
sets as in class LCGPointSet , but implemented differently. 
KorobovLatticeSequence  This class implements Korobov lattice sequences, defined as follows. 
LCGPointSet  Implements a recurrencebased point set defined via a linear congruential recurrence of the form x_{i} = ax_{i1}mod n and u_{i} = x_{i}/n. 
LMScrambleShift  This class implements a
PointSetRandomization
that performs a left matrix scrambling and adds a random digital
shift. 
NiedSequenceBase2  This class implements digital sequences constructed from the Niederreiter sequence in base 2. 
NiedXingSequenceBase2  This class implements digital sequences based on the NiederreiterXing sequence in base 2. 
PaddedPointSet  This container class realizes padded point sets, constructed by taking some coordinates from a point set P_{1}, other coordinates from a point set P_{2}, and so on. 
PointSet  This abstract class defines the basic methods for accessing and manipulating point sets. 
RadicalInverse  This class implements basic methods for working with radical inverses of integers in an arbitrary basis b. 
RandomShift  This class implements a
PointSetRandomization . 
RandomStart  This class implements a
PointSetRandomization
that randomizes a sequence with a random starting point. 
RandShiftedPointSet  This container class embodies a point set to which a random shift modulo 1 is applied (i.e., a single uniform random point is added to all points, modulo 1, to randomize the inner point set). 
Rank1Lattice  This class implements point sets specified by integration lattices of rank 1. 
RQMCPointSet  This class is used for randomized quasiMonte Carlo (RQMC) simulations. 
SMScrambleShift  This class implements a
PointSetRandomization
that performs a striped matrix scrambling and adds a random
digital shift. 
SobolSequence  This class implements digital nets or digital sequences in base 2 formed by the first n = 2^{k} points of a Sobol' sequence. 
SubsetOfPointSet  This container class permits one to select a subset of a point set. 
This package provides classes implementing highly uniform point sets (HUPS) over the sdimensional unit hypercube [0, 1)^{s}, and tools for their randomization. The terminology lowdiscrepancy sequence (LDS) is often used for infinite sequences of points such that the discrepancy between the distribution of the first n points of the sequence and the uniform distribution converges to zero at a certain rate when n > ∞. HUPS and LDS are used for quasiMonte Carlo integration, as we now briefly explain. See, e.g., for further details.
Suppose we want to estimate the integral of a function f defined over the sdimensional unit hypercube,
Practically any mathematical expectation that can be estimated by simulation can be written in this way, usually for a very complicated f and sometimes for s = ∞. Indeed, the source of randomness of stochastic simulations is usually a stream of real numbers u = (u_{0}, u_{1}, u_{2},...) whose purpose is to imitate i.i.d. U(0, 1) random variables. These real numbers are transformed in complicated ways to produce the estimator. Thus, the dimension s of the integral represents the number of calls to the uniform random number generator if that number is deterministic. If it is random and unbounded, we take s = ∞. In the latter case, however, we can assume that the actual number of calls is finite with probability one (otherwise the simulation may never end).We consider an estimator of μ of the form
which is the average of f over the point set P_{n} = {u_{0},...,u_{n1}}⊂[0, 1)^{s}.With the Monte Carlo (MC) method, the u_{i}'s are i.i.d. random vectors uniformly distributed over [0, 1)^{s}. Then, Q_{n} is an unbiased estimator of μ with variance σ^{2}/n, where
QuasiMonte Carlo (QMC) methods use point sets P_{n} that are more evenly distributed over the unit hypercube than typical random points. We call them highly uniform point sets (HUPS). The aim is to reduce the size of the integration error Q_{n}  μ. Two important classes of methods for constructing such point sets are digital nets and integration lattices. Both are implemented in this package, in various flavors.
To give an idea of how HUPS and LDS can be constructed, we start with a simple onedimensional example. If s = 1 and n is fixed, very simple highly uniform constructions are the point sets P_{n} = {0, 1/n, ...,(n  1)/n} and the shifted version P'_{n} = {1/(2n), 3/(2n), ...,(2n  1)/(2n)}.
In s > 1 dimensions, the simplest extensions would be as follows. Let n = d^{s} for some integer d and define P_{n} as the Cartesian product of s copies of the onedimensional sets P_{d}; that is, P_{n} = {(u_{0},..., u_{s1}) : u_{j}∈{0, 1/d, ...,(d  1)/d} for each j}, and similarly for P'_{n}. The point sets thus obtained are regular rectangular grids. Unfortunately, this approach breaks down rapidly when s gets large, because n must increase exponentially fast with s for fixed d. Another important drawback is that when P_{n} is projected over lowerdimensional subspaces, several points are projected onto each other and become redundant.
A better idea is to construct a point set P_{n} in s dimensions such that each onedimensional projection of P_{n} is the set of values {0, 1/n, ...,(n  1)/n}. Of course, these values should not be visited in the same order for all coordinates, because otherwise all the points would lie on the diagonal line going from (0,..., 0) to (1,..., 1). In other words, for each coordinate j, 0 <= j < s, we must define a different permutation of the integers {0,..., n  1} and visit the values {0, 1/n, ...,(n  1)/n} in the order determined by that permutation. The trick is to select those permutations in a way that P_{n} itself is highly uniform over [0, 1)^{s} in a welldefined sense (to be determined). This is what most construction methods attempt to achieve. Before looking at concrete ways of defining such permutations, we introduce a related issue: what to do if n is not fixed.
For s = 1, a simple way of filling up the unit interval [0, 1) uniformly is via the lowdiscrepancy sequence 0, 1/2, 1/4, 3/4, 1/8, 5/8, 3/8, 7/8, 1/16, 9/16, ..., called the van der Corput sequence in base 2. More generally, select an integer b >= 2, called the base. The radical inverse function in base b, ψ_{b} : N > [0, 1), is defined as follows. If i is a kdigit integer in base b with digital bary expansion
For s > 1, one could either take different (relatively prime) bases for the different coordinates, or take the same basis b but permute the successive values using a different permutation for each coordinate. These permutations are usually selected in a way that for every integer k, the first b^{k} values that are enumerated remain the same (they are the values of ψ_{b}(i) for i = 0,..., b^{k}  1), but they are enumerated in a different order. Several digital net constructions (to be defined later) fit this framework.
If we decide to take different bases, the most natural choice is to take the jth smallest prime, b_{j}, as a base for coordinate j  1; that is, base 2 for coordinate 0, base 3 for coordinate 1, base 5 for coordinate 2, and so on. The infinite sequence thus defined, where point i is
In the case where n is fixed, we can always take i/n as the first coordinate of point i. In particular, the Hammersley point set with n points in s dimensions contains the points
Digital nets and sequences are an important class of HUPS and LDS constructions. Most concrete implementations, e.g., those proposed by Sobol', Faure, Niederreiter, and Niederreiter and Xing, are linear digital nets and sequences, defined as follows (see also).
Let b >= 2 be an arbitrary integer (usually a prime number), called
the base.
A net that contains n = b^{k} points in s dimensions is defined via
s generator matrices
C_{0},...,C_{s1}, which are
(in theory)
∞×k matrices whose elements are in
Z_{b} = {0,..., b  1}.
The matrix C_{j} is used for coordinate j of all the points, for j >= 0.
To define the ith point
u_{i}, for
i = 0,..., b^{k}  1, write
the digital expansion of i in base b and multiply the vector of its
digits by C_{j} to obtain the digits of the expansion of u_{i, j},
the jth coordinate of
u_{i}. That is,
i  =  ∑_{=0}^{k1}a_{i,[tex2html_wrap_indisplay1174]}b^{[tex2html_wrap_indisplay1175]},  
u_{i, j, 1}u_{i, j, 2}[tex2html_wrap_indisplay1179] [tex2html_wrap_indisplay1181]  =  C_{j}[tex2html_wrap_indisplay1184]a_{i, 0}[tex2html_wrap_indisplay1185]a_{i, 1}[tex2html_wrap_indisplay1186] [tex2html_wrap_indisplay1187] [tex2html_wrap_indisplay1188]a_{i, k1}[tex2html_wrap_indisplay1189],  
u_{i, j}  =  ∑_{[tex2html_wrap_indisplay1193]=1}^{∞}u_{i, j,[tex2html_wrap_indisplay1194]}b^{[tex2html_wrap_indisplay1195]},  
[tex2html_wrap_indisplay1197]_{i}  =  (u_{i, 0},..., u_{i, s1}). 
Usually, the first k lines of each C_{j} form a nonsingular k×k matrix. Then, the n output values for coordinate j, u_{0, j},..., u_{n1, j}, when truncated to their first k fractional digits in base b, are a permutation of the numbers 0, 1/n,...,(n  1)/n. Different coordinates simply use different permutations, implemented via the matrices C_{j}.
When the first k lines of C_{j} form the identity and the other lines are zero, the first n output values are the first n elements of the van der Corput sequence in base b. If we reverse the order of the columns of that matrix C_{j} (i.e., column c will contain a one in line k  c + 1 and zeros elsewhere, for 0 <= c < k), we obtain the output values 0, 1/n,...,(n  1)/n in that order. With a slight abuse of language, we shall call this first matrix (with the identity followed by lines of zeros) the identity and the second one (with the columns in reverse order) the reflected identity. It is customary to take C_{0} as the identity for digital sequences, and often for digital nets as well. But for digital nets (where n is fixed in advance), one can take C_{0} as the reflected identity instead, then C_{1} as the identity, and so on. That is, the matrix C_{j} for the digital net is taken as the matrix C_{j1} of the digital sequence. Our package often gives the choice.
For digital sequences, the matrices C_{j} actually have an infinite number of columns, although only the first k columns are needed to generate the first b^{k} points. So in practice, we never need to store more than a finite number of columns at a time. Whenever we find that we need more than b^{k} points for the current value of k, we can simply increase k and add the corresponding columns to the matrices C_{j}.
The classes DigitalNet
and DigitalSequence
implement generic digital nets and sequences.
Specific instances are constructed in subclasses of these two classes.
An integration lattice is a discrete (but infinite) subset of R^{s} of the form
Let V be the matrix whose rows are the basis vectors v_{1},^{ ... },v_{s} and V^{1} its inverse. One has Z^{s}⊆L_{s} if and only if all entries of V^{1} are integer. When this holds, n = det(V^{1}) and all entries of V are multiples of 1/n.
The rank of the lattice is the smallest r such that one can find a basis of the form v_{1},...,v_{r},e_{r+1},^{ ... },e_{s}, where e_{j} is the jth unit vector in s dimensions. In particular, a lattice rule of rank 1 has a basis of the form v_{1} = (a_{1},..., a_{s})/n and v_{j} = e_{j} for j > 1, where a_{j}∈Z_{n} for each j. It is a Korobov rule if v_{1} has the special form v_{1} = (1, a, a^{2}mod n, ..., a^{s1}mod n)/n for some a∈Z_{n}. The point set P_{n} of a Korobov lattice rule can also be written as P_{n} = {(x_{1},..., x_{s})/n such that x_{1}∈Z_{n} and x_{j} = ax_{j1}mod n for all j > 1}. This is the set of all vectors of successive values produced by a linear congruential generator (LCG) with modulus n and multiplier a, from all possible initial states, including 0. In this case, the points are easy to enumerate by using the recurrence.
Certain types of point sets are defined pretty much like random number generators: choose a finite state space S, a transition function f : S > S, an output function g : S > [0, 1), and define
Let us assume that n is finite and that for each s_{0}∈S, the recurrence s_{j} = f (s_{j1}) is purely periodic, i.e., there is always an integer j such that s_{j} = s_{0}. The smallest such j, called the period length, depends in general on s_{0}. Thus, the state space S is partitioned into a finite number of cycles. The successive coordinates of any point u∈P_{n} are periodic with period length equal to the length of the cycle that contains s_{0} (and the following s_{j}'s).
One way of implementing such a point set while avoiding to recompute
f and g each time a coordinate is needed is to store explicitly
all the cycles of the recurrence, in the form of a list of cycles.
We can store either the successive u_{j}'s directly, or the successive
s_{j}'s, over each cycle.
The class CycleBasedPointSet
provides the framework for doing that.
For example, a Korobov lattice point set is defined via the recurrence
x_{j} = ax_{j1}mod n and output function
u_{j} = x_{j}/n.
If n is prime and a is a primitive element modulo n, then there
are two cycles: one of period 1 that contains only 0, and
the other of period n  1.
For more general n and a, there will be more cycles.
The class LCGPointSet
constructs this type of point set and
stores explicitly the successive values of u_{j} over the different cycles.
There are cases where n is a power of two, say n = 2^{k},
and where the state s_{j} is represented as a kbit string.
In that context, it is often more convenient to store the successive
s_{j}'s instead of the successive u_{j}'s, over the set of cycles
(e.g., if a random digital shift in base 2 is to be applied to
randomize the points, it can be performed by applying a bitwise xor
directly to s_{j}).
When generating the coordinates, the s_{j}'s can be interpreted as
2^{k}bit integers and multiplied by 2^{k} to produce the output.
This is supported by the class
CycleBasedPointSetBase2
.
Special instances of this class are usually based on
linear recurrences modulo 2 and they include the Korobovtype
polynomial lattice rules.
In their original versions, these HUPS are deterministic, and the corresponding QMC methods give a deterministic integration error that is difficult to estimate. In randomized QMC methods, P_{n} is randomized, preferably in a way that it retains its high uniformity over [0, 1)^{s} when taken as a set, while each of its points has the uniform distribution over [0, 1)^{s} when taken individually. Then, Q_{n} becomes an unbiased estimator of μ, hopefully with smaller variance than the standard MC estimator. To estimate the variance and compute a confidence interval on μ, one can apply m independent randomizations to the same P_{n}, and compute bar(X)_{m} and S_{m, x}^{2}, the sample mean and sample variance of the m corresponding (independent) copies of Q_{n}. Then, E[bar(X)_{m}] = μ and E[S_{m, x}^{2}] = Var[Q_{n}] = mVar[bar(X)_{m}].
Two examples of such randomizations are the random shift modulo 1,
proposed in and implemented in class
RandShiftedPointSet
,
and the random digital shift in base b, described and
implemented in class DigitalNet
.
These randomizations are also incorporated directly in certain types
of point sets such as
CycleBasedPointSet
,
CycleBasedPointSetBase2
, etc.
In the random shift modulo 1, we generate a single point u uniformly over [0, 1)^{s} and add it to each point of P_{n}, coordinatewise, modulo 1. Since all points of P_{n} are shifted by the same amount, the set retains most of its structure and uniformity.
For the random digital shift in base b, we generate again a single u = (u_{1},..., u_{s}) uniformly over [0, 1)^{s}, write the digital expansion in base b of each of its coordinates, say u_{j} = ∑_{[tex2html_wrap_inline1364]=1}^{∞}d_{j,[tex2html_wrap_inline1365]}b^{[tex2html_wrap_inline1366]}, then add d_{j,[tex2html_wrap_inline1368]} modulo b to the th digit of the digital expansion in base b of the jth coordinate of each point u_{i}∈P_{n}. For b = 2, the digitwise addition modulo b becomes a bitwise exclusiveor, which is fast to perform on a computer.
An interesting property of the digital shift in base b is that
if the hypercube [0, 1)^{s} is partitioned into
b^{q1+ ... +qs}
rectangular boxes
of the same size by partitioning the jth axis into b^{qj} equal
parts for each j, for some integers q_{j} >= 0
(such a partition is called a qequidissection in base b
of the unit hypercube, where
q = (q_{1},..., q_{s})), then the number of
boxes that contain m points, for each integer m, is unchanged by
the randomization. In particular, if each box contains the same number
of points of P_{n} before the randomization, then it also does after the
randomization.
In this case, we say that P_{n} is qequidistributed in base b.
Several other randomization methods exist and most are adapted to
special types of point sets; see, e.g.,
DigitalNet
and.
Let
u_{i} = (u_{i, 0}, u_{i, 1},..., u_{i, s1})
be the elements of the point set P_{n}, for
i = 0,..., n  1.
Both the number of points n and the dimension s can be finite
or infinite. The point set can be viewed as a twodimensional array
whose element (i, j) contains u_{i, j}, the coordinate j of point i.
In the implementations of typical point sets, the values u_{i, j} are
not stored explicitly in a twodimensional array, but pertinent
information is organized so that the points and their coordinates can be
generated efficiently. The base class for point sets is the abstract
class PointSet
.
To enumerate the successive points or the successive coordinates of a
given point, we use point set iterators,
which resemble the iterators defined
in Java collections, except that they loop over bidimensional sets.
Their general behavior is defined in the interface
PointSetIterator
.
Several independent iterators can coexist at any given time for the same
point set. Each one maintains a current point index and a current
coordinate index, which are incremented by one when the iterator advances
to the next point or the next coordinate. Both are initialized to 0.
Each subclass of PointSet
has its own implementation of
PointSetIterator
and has
a method iterator that
creates and returns a new point set iterator of the correct type.
An important feature of the
PointSetIterator
interface is that
it extends the RandomStream
interface. This means that any
point set iterator can be used in place of a random stream that is
supposed to generate i.i.d. U(0, 1) random variables, anywhere in a
simulation program. It then becomes very easy to replace the
(pseudo)random numbers by the coordinates u_{i, j} of a randomized HUPS
without changing the internal code of the simulation program.
HUPS are often transformed either deterministically or randomly. Deterministic transformations can be applied to improve the uniformity, or to eliminate some points or coordinates (i.e., selecting subsets), or to concatenate point sets (padding), or to take an antithetic version of a point set, etc. Random transformations are used for randomized QMC. They are also useful when the average of a uniformity measure of interest over the outcomes of a certain type of randomization is much better than the worst case and may be better than the uniformity measure of the original point set. When a point set is transformed, we often want to keep the original as well, and we may want to apply different types of transformations or different independent randomizations to the same point set.
This can be achieved via container point sets, which are defined
in terms of another point set to which they keep a reference
and apply certain transformations.
ContainerPointSet
is the base
class for such containers.
One example is RandShiftedPointSet
,
which applies a random shift
modulo 1 to the point set that it contains.
Of course, the contained point set can be a container itself and this
can be done recursively, but too many levels of recursiveness may impair
the performance (speed).
To be done...

SSJ V. 2.6. 

PREV PACKAGE NEXT PACKAGE  FRAMES NO FRAMES 