HaltonSequence.java

/*
 * Class:        HaltonSequence
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since
 *
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
package umontreal.ssj.hups;

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);

   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);
   }

   public void setStart(double[] x0) {
      for (int i = 0; i < dim; i++)
         start[i] = RadicalInverse.radicalInverseInteger(base[i], x0[i]);
   }

   public void init(double[] x0) {
      radinv = new RadicalInverse[dim];
      for (int i = 0; i < dim; i++)
         radinv[i] = new RadicalInverse(base[i], x0[i]);
   }

   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;
   }

   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;
   }

   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);
      }
   }
}