umontreal.iro.lecuyer.hups.HaltonSequence Maven / Gradle / Ivy
Show all versions of ssj Show documentation
/*
* Class: HaltonSequence
* Description:
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Université de Montréal
* Organization: DIRO, Université de Montréal
* @author
* @since
* SSJ is free software: you can redistribute it and/or modify it under
* the terms of the GNU General Public License (GPL) as published by the
* Free Software Foundation, either version 3 of the License, or
* any later version.
* SSJ is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* A copy of the GNU General Public License is available at
GPL licence site.
*/
package umontreal.iro.lecuyer.hups;
/**
* 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.
* The ith point in s dimensions is
*
*
*
* 
i = (ψb1(i), ψb2(i),..., ψbs(i)),
*
* for
* i = 0, 1, 2,..., where ψb is the radical inverse function
* in base b, defined in class {@link RadicalInverse}, and where
*
* 2 = b1 < ... < bs are the s smallest prime numbers in
* increasing order.
*
*
* A fast method is implemented to generate randomized Halton sequences, starting from an arbitrary point x0.
*
*
* The points can be ``scrambled'' by applying a permutation to the
* digits of i before computing each coordinate, in the same way as for the class
* {@link HammersleyPointSet}, for all coordinates j >= 0.
*
*/
public class HaltonSequence extends PointSet {
private int[] base; // Vector of prime bases.
private int[][] permutation; // Digits permutation, for each dimension.
private boolean permuted; // Permute digits?
private RadicalInverse[] radinv; // Vector of RadicalInverse's.
private int[] start; // starting indices
private final static int positiveBitMask = ~Integer.reverse(1);
/**
* Constructs a new Halton sequence
* in dim dimensions.
*
* @param dim dimension
*
*/
public HaltonSequence (int dim) {
if (dim < 1)
throw new IllegalArgumentException
("Halton sequence must have positive dimension dim");
this.dim = dim;
numPoints = Integer.MAX_VALUE;
base = RadicalInverse.getPrimes (dim);
start = new int[dim];
java.util.Arrays.fill(start, 0);
}
/**
* Initializes the Halton sequence starting at point x0.
* For each coordinate j, the sequence starts at index ij such that
* x0[j] is the radical inverse of ij.
* The dimension of x0 must be at least as large as the dimension
* of this object.
*
* @param x0 starting point of the Halton sequence
*
*
*/
public void setStart (double[] x0) {
for (int i = 0; i < dim; i++)
start[i] = RadicalInverse.radicalInverseInteger(base[i], x0[i]);
}
/**
* Initializes the Halton sequence starting at point x0.
* The dimension of x0 must be at least as large as the dimension
* of this object.
*
* @param x0 starting point of the Halton sequence
*
*
*/
public void init (double[] x0) {
radinv = new RadicalInverse[dim];
for (int i = 0; i < dim; i++)
radinv[i] = new RadicalInverse (base[i], x0[i]);
}
/**
* Permutes the digits using permutations from for all coordinates.
* After the method is called, the coordinates ui, j are generated via
*
*
*
* ui, j = ∑r=0k-1πj[ar]bj-r-1,
*
* for
* j = 0,..., s - 1,
* where πj is the Faure-Lemieux (2008) permutation of
* {0,..., bj -1}.
*
*/
public void addFaureLemieuxPermutations() {
permutation = new int[dim][];
for (int i = 0; i < dim; i++) {
permutation[i] = new int[base[i]];
RadicalInverse.getFaureLemieuxPermutation (i, permutation[i]);
}
permuted = true;
}
/**
* Permutes the digits using Faure permutations for all coordinates.
* After the method is called, the coordinates ui, j are generated via
*
*
*
* ui, j = ∑r=0k-1πj[ar]bj-r-1,
*
* for
* j = 0,..., s - 1,
* where πj is the Faure permutation of
* {0,..., bj -1}.
*
*/
public void addFaurePermutations() {
permutation = new int[dim][];
for (int i = 0; i < dim; i++) {
permutation[i] = new int[base[i]];
RadicalInverse.getFaurePermutation (base[i], permutation[i]);
}
permuted = true;
}
/**
* Erases the permutations: from now on, the digits will not be
* permuted.
*
*/
public void ErasePermutations() {
permuted = false;
permutation = null;
}
public int getNumPoints () {
return Integer.MAX_VALUE;
}
public double getCoordinate (int i, int j) {
if (radinv != null) {
if (!permuted) {
return radinv[j].nextRadicalInverse ();
} else {
throw new UnsupportedOperationException (
"Fast radical inverse is not implemented in case of permutation");
}
} else {
int k = start[j] + i;
// if overflow, restart at first nonzero point
// (Struckmeier restarts at zero)
if (k < 0)
k = (k & positiveBitMask) + 1;
if (permuted)
return RadicalInverse.permutedRadicalInverse
(base[j], permutation[j], k);
else
return RadicalInverse.radicalInverse (base[j], k);
}
}
}