BatchMeansSim.java

/*
 * Class:        BatchMeansSim
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Éric Buist
 * @since        2007
 
 * 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
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.simexp;
 
import java.lang.ref.SoftReference;
import cern.colt.list.DoubleArrayList;
import cern.colt.matrix.DoubleMatrix1D;
import cern.colt.matrix.DoubleMatrix2D;
import umontreal.ssj.simevents.Event;
import umontreal.ssj.simevents.Simulator;
import umontreal.ssj.util.Misc;

public abstract class BatchMeansSim extends SimExp {
   private DoubleArrayList batchTimes = new DoubleArrayList();
   private boolean batchLengths;
   private boolean aggregation;
   private int minBatches;
   private int maxBatches;
   private double warmupTime;
   private EndSimEvent endSimEvent;
 
   private boolean warmupDone;
   private int targetBatches;
   private int doneBatches;
   private int droppedBatches;
   private double batchFraction;
   private double batchSizeMultiplier;
   private boolean aggregateUpdated;
   private double batchSize;
 
   private int nAgr;

   public BatchMeansSim(int minBatches, double batchSize, double warmupTime) {
      this(minBatches, Integer.MAX_VALUE, batchSize, warmupTime);
   }

   public BatchMeansSim(int minBatches, int maxBatches, double batchSize, double warmupTime) {
      this(Simulator.getDefaultSimulator(), minBatches, maxBatches, batchSize, warmupTime);
   }

   public BatchMeansSim(Simulator sim, int minBatches, double batchSize, double warmupTime) {
      this(sim, minBatches, Integer.MAX_VALUE, batchSize, warmupTime);
   }

   public BatchMeansSim(Simulator sim, int minBatches, int maxBatches, double batchSize, double warmupTime) {
      super(sim);
      if (minBatches <= 0)
         throw new IllegalArgumentException("minBatches <= 0");
      if (maxBatches < minBatches)
         throw new IllegalArgumentException("maxBatches < minBatches");
      if (warmupTime < 0)
         throw new IllegalArgumentException("Warmup time must not be negative");
      if (batchSize <= 0)
         throw new IllegalArgumentException("Batch size must not be 0 or negative");
      this.minBatches = minBatches;
      this.maxBatches = maxBatches;
      this.warmupTime = warmupTime;
      this.batchSize = batchSize;
      targetBatches = minBatches;
   }

   public boolean getBatchAggregation() {
      return aggregation;
   }

   public void setBatchAggregation(boolean a) {
      if (warmupDone && simulating)
         throw new IllegalStateException("Cannot change the aggregation status during simulation");
      aggregation = a;
   }

   public boolean getBatchLengthsKeeping() {
      return batchLengths;
   }

   public void setBatchLengthsKeeping(boolean b) {
      if (warmupDone && simulating)
         throw new IllegalStateException("Cannot change the aggregation status during simulation");
      batchLengths = b;
   }

   public int getMinBatches() {
      return minBatches;
   }

   public void setMinBatches(int minBatches) {
      if (minBatches <= 0)
         throw new IllegalArgumentException("minBatches <= 0");
      this.minBatches = minBatches;
      if (maxBatches < minBatches)
         maxBatches = minBatches;
   }

   public int getMaxBatches() {
      return maxBatches;
   }

   public void setMaxBatches(int maxBatches) {
      if (maxBatches < minBatches)
         throw new IllegalArgumentException("maxBatches < minBatches");
      this.maxBatches = maxBatches;
   }

   public double getBatchSize() {
      return batchSize;
   }

   public void setBatchSize(double batchSize) {
      if (batchSize <= 0)
         throw new IllegalArgumentException("batchSize <= 0");
      this.batchSize = batchSize;
   }

   public double getWarmupTime() {
      return warmupTime;
   }

   public void setWarmupTime(double warmupTime) {
      if (warmupTime < 0)
         throw new IllegalArgumentException("warmupTime < 0");
      this.warmupTime = warmupTime;
   }

   public double getBatchFraction() {
      return batchFraction;
   }

   public double getBatchSizeMultiplier() {
      return batchSizeMultiplier;
   }

   public int getTargetBatches() {
      return targetBatches;
   }

   public void setTargetBatches(int targetBatches) {
      if (targetBatches < minBatches)
         throw new IllegalArgumentException("Target number of batches too small");
      if (targetBatches > maxBatches)
         throw new IllegalArgumentException("Target number of batches too large");
      this.targetBatches = targetBatches;
   }

   public int getCompletedRealBatches() {
      return doneBatches;
   }

   public int getDroppedRealBatches() {
      return droppedBatches;
   }

   public void dropFirstRealBatches(int n) {
      if (n < 0 || n > doneBatches - droppedBatches)
         throw new IllegalArgumentException(
               "Cannot drop less than 0 or more than " + (doneBatches - droppedBatches) + " real batches");
      if (batchLengths || aggregation)
         batchTimes.removeFromTo(0, n - 1);
      droppedBatches += n;
   }

   public int getBatch(double time) {
      if (!warmupDone)
         return -1;
      if (!batchLengths && !aggregation)
         return -1;
      if (batchTimes.size() == 0)
         return -1;
      return Misc.getTimeInterval(batchTimes.elements(), 0, doneBatches, time);
   }

   public boolean isWarmupDone() {
      return warmupDone;
   }

   public int getNumAggregates() {
      if (nAgr == 0)
         throw new IllegalStateException("Number of real batches not available");
      return nAgr;
   }

   public double getRealBatchLength(int batch) {
      if (batchLengths || aggregation)
         // The first get will fail if batch + 1 is out of bounds.
         // If the first get succeeds, the second will, so we can
         // use getQuick to eliminate a check.
         return batchTimes.get(batch + 1 - droppedBatches) - batchTimes.getQuick(batch - droppedBatches);
      else if (batch != doneBatches - 1)
         throw new IllegalArgumentException("Unavailable batch length");
      else
         return batchTimes.getQuick(1) - batchTimes.getQuick(0);
   }

   public double getRealBatchStartingTime(int batch) {
      if (batchLengths || aggregation)
         return batchTimes.get(batch);
      else if (batch != doneBatches - 1)
         throw new IllegalArgumentException("Unavailable batch time");
      else
         return batchTimes.getQuick(0);
   }

   public double getRealBatchEndingTime(int batch) {
      if (batchLengths || aggregation)
         return batchTimes.get(batch + 1);
      else if (batch != doneBatches - 1)
         throw new IllegalArgumentException("Unavailable batch time");
      else
         return batchTimes.getQuick(1);
   }

   public void allocateCapacity(int capacity) {
      throw new UnsupportedOperationException();
   }

   public void regroupRealBatches(int x) {
      throw new UnsupportedOperationException();
   }
 
   private final void allocateCapacity() {
      // Preallocate the necessary memory to store
      // all the observations when using fixed
      // number of batches.
      // A maximal number of batches is used
      // to prevent memory overflow problems.
      // Unfortunately, in Java, we can only
      // use heuristics to prevent the memory
      // problem. In addition to a maximum
      // number of batches, we use a soft reference
      // that will be cleared before the OutOfMemoryError
      // is thrown.
      // The size of the referent can be tuned; it is only
      // an heuristic.
      SoftReference<?> softRef = new SoftReference<double[]>(new double[50000]);
      batchFraction = 1.0;
      batchSizeMultiplier = 1.0;
      int currentCapacity = doneBatches;
      while (currentCapacity < targetBatches && softRef.get() != null) {
         // Try to double the capacity
         int newCapacity = Math.min(2 * currentCapacity + 1, targetBatches);
         // Since we want all collectors to have the same capacity
         // we must not stop the for loop if the softRef
         // is cleared.
         if (aggregation || batchLengths)
            batchTimes.ensureCapacity(newCapacity + 1);
         try {
            allocateCapacity(newCapacity);
         } catch (UnsupportedOperationException use) {
            return;
         }
         currentCapacity = newCapacity;
      }
      if (currentCapacity < targetBatches && softRef.get() == null) {
         if (!aggregation)
            throw new IllegalStateException("Cannot increase batch size");
         // The soft reference was cleared.
         // The memory taken by the simulation
         // could crash the program.
         // To try avoiding that, we will stop taking memory
         // by increasing the batch size.
         double ftargetBatches = targetBatches;
         double fdoneBatches = doneBatches;
         while (targetBatches > currentCapacity) {
            batchSizeMultiplier *= 2;
            // From now, we can have a fraction of batches done
            fdoneBatches /= 2;
            ftargetBatches /= 2;
            targetBatches = (int) ftargetBatches;
            doneBatches = (int) fdoneBatches;
            if (aggregation || batchLengths) {
               for (int i = 0; i < batchTimes.size() / 2; i++)
                  batchTimes.setQuick(i, batchTimes.getQuick(2 * i));
               batchTimes.setSize(batchTimes.size() / 2);
            }
            regroupRealBatches(2);
            // Compute the remaining fraction of batch to simulate
            batchFraction = 1.0 - (fdoneBatches - doneBatches);
         }
         // Round the new total number of batches to a multiple
         // of minBatches.
         if (targetBatches % minBatches != 0) {
            targetBatches /= minBatches;
            targetBatches *= minBatches;
            targetBatches += minBatches;
         }
      }
   }

   public abstract void initSimulation();

   public abstract void initBatchStat();

   public abstract void initRealBatchProbes();

   public abstract void initEffectiveBatchProbes();

   public abstract void addRealBatchObs();

   public abstract void addEffectiveBatchObs(int s, int h, double l);

   public int getRequiredNewBatches() {
      return 0;
   }

   public void init() {
      if (simulating)
         throw new IllegalStateException("Already simulating");
      nAgr = 0;
      warmupDone = false;
      doneBatches = 0;
      droppedBatches = 0;
      if (targetBatches < minBatches)
         targetBatches = minBatches;
      if (aggregation && targetBatches % minBatches != 0) {
         targetBatches /= minBatches;
         targetBatches += minBatches;
         if (targetBatches + minBatches <= maxBatches)
            targetBatches += minBatches;
      }
      aggregateUpdated = false;
      simulator().init();
      endSimEvent = new EndSimEvent(simulator());
      initSimulation();
   }

   public Event getEndSimEvent() {
      return endSimEvent;
   }

   public void warmup() {
      warmup(warmupTime);
   }

   public void warmup(double warmupTime) {
      if (warmupDone)
         throw new IllegalStateException("Warmup already done");
      if (!Double.isInfinite(warmupTime) && !Double.isNaN(warmupTime))
         endSimEvent.schedule(warmupTime);
      simulator().start();
      endSimEvent.cancel();
      initRealBatchProbes();
      if (!aggregation)
         initEffectiveBatchProbes();
      if (aggregation || batchLengths) {
         batchTimes.setSize(0);
         batchTimes.add(simulator().time());
      } else {
         batchTimes.setSize(2);
         batchTimes.trimToSize();
         batchTimes.setQuick(0, simulator().time());
         batchTimes.setQuick(1, simulator().time());
      }
      warmupDone = true;
   }

   public void simulateBatch() {
      double f = getBatchSizeMultiplier();
      if (f > 1)
         batchSize *= f;
      simulateBatch(batchSize * getBatchFraction());
   }

   public void simulateBatch(double batchLength) {
      initBatchStat();
      endSimEvent.schedule(batchLength);
      simulator().start();
      endSimEvent.cancel();
      batchFraction = 1.0;
      batchSizeMultiplier = 1.0;
      aggregateUpdated = false;
      ++doneBatches;
      if (batchLengths || aggregation)
         batchTimes.add(simulator().time());
      else {
         batchTimes.setQuick(0, batchTimes.getQuick(1));
         batchTimes.setQuick(1, simulator().time());
      }
      addRealBatchObs();
      if (!aggregation) {
         nAgr = 1;
         addEffectiveBatchObs(doneBatches - 1, 1, getRealBatchLength(doneBatches - 1));
      }
   }

   public void adjustTargetBatches(int numNewBatches) {
      if (doneBatches + numNewBatches > maxBatches) {
         numNewBatches = maxBatches - doneBatches;
         if (aggregation && numNewBatches < minBatches)
            // We must simulate at least minBatches real batches, so
            // do nothing
            return;
      }
      if (aggregation) {
         // For effective batches to change,
         // we need to simulate at least minBatches additional real batches.
         // The number of new batches is therefore rounded to
         // the smallest multiple greater than or
         // equal to minBatches.
         if (numNewBatches % minBatches != 0) {
            numNewBatches /= minBatches;
            numNewBatches *= minBatches;
            if (doneBatches + numNewBatches + minBatches <= maxBatches)
               numNewBatches += minBatches;
         }
      }
      targetBatches = doneBatches + numNewBatches;
   }

   public void simulateBatches() {
      if (doneBatches < targetBatches)
         // Increased the target number of required batches
         // and allocate memory to store the observations
         // if this is required.
         allocateCapacity();
      while (doneBatches < targetBatches)
         // targetBatches may be changed by a simulation event.
         // This can be used to stop the simulation before the
         // desired number of batches is obtained, e.g., when
         // the simulated system seems unstable.
         simulateBatch();
      if (aggregation)
         makeAggregateBatches();
   }

   public void simulate() {
      init();
      simulating = true;
      try {
         warmup(warmupTime);
         while (doneBatches < targetBatches) {
            simulateBatches();
            adjustTargetBatches(getRequiredNewBatches());
         }
         if (aggregation && !aggregateUpdated)
            makeAggregateBatches();
      } finally {
         simulating = false;
      }
   }
 
   private final void makeAggregateBatches() {
      // When using fixed number of batches, we must
      // add the observations to statistical collectors
      // only at the end of the simulation.
      initEffectiveBatchProbes();
      // nAgr gives the number of real batches in each effective batch
      nAgr = (doneBatches - droppedBatches) / minBatches;
      // tnb gives the number of observations in the statistical counters
      int tnb = minBatches;
      if (nAgr == 0) {
         // doneBatches is smaller than minBatches, but we would like
         // observations to be collected and displayed.
         // We will have less than minBatches observations.
         nAgr = 1;
         tnb = doneBatches - droppedBatches;
      }
      for (int i = 0; i < tnb; i++) {
         int startAgr = i * nAgr + droppedBatches;
         addEffectiveBatchObs(startAgr, nAgr,
               batchTimes.getQuick(startAgr - droppedBatches + nAgr) - batchTimes.getQuick(startAgr - droppedBatches));
      }
      aggregateUpdated = true;
   }
 
   private static final class EndSimEvent extends Event {
      public EndSimEvent(Simulator sim) {
         super(sim);
      }
 
      public void actions() {
         simulator().stop();
      }
   }
 
   public String toString() {
      StringBuffer sb = new StringBuffer(getClass().getName());
      sb.append('[');
      sb.append("minimal number of batches: ").append(minBatches);
      if (maxBatches < Integer.MAX_VALUE)
         sb.append(", maximal number of batches: ").append(maxBatches);
      sb.append(", target number of batches: ").append(targetBatches);
      sb.append(", warmup time: ").append(warmupTime);
      if (simulating)
         sb.append(", simulation in progress");
      else
         sb.append(", simulation stopped");
      if (warmupDone)
         sb.append(", number of completed batches: ").append(doneBatches);
      else
         sb.append(", warmup not completed");
      if (aggregation)
         sb.append(", batch aggregation turned ON");
      else {
         sb.append(", batch aggregation turned OFF");
         if (batchLengths)
            sb.append(", batch lengths are kept");
         else
            sb.append(", batch lengths are not kept");
      }
      sb.append(']');
      return sb.toString();
   }

   public static double getSum(double[] a, int start, int length) {
      if (start < 0)
         throw new IndexOutOfBoundsException("start must not be negative");
      if (length < 0)
         throw new IllegalArgumentException("length is negative");
      if (start + length > a.length)
         throw new IndexOutOfBoundsException("Not enough elements in the array");
      double s = 0;
      for (int i = start; i < start + length; i++)
         s += a[i];
      return s;
   }

   public static double getSum(DoubleArrayList l, int start, int length) {
      if (start < 0)
         throw new IndexOutOfBoundsException("start must not be negative");
      if (length < 0)
         throw new IllegalArgumentException("length is negative");
      if (start + length > l.size())
         throw new IndexOutOfBoundsException("Not enough elements in the array list");
      double s = 0;
      for (int i = start; i < start + length; i++)
         s += l.getQuick(i);
      return s;
   }

   public static double getSum(DoubleMatrix1D m, int start, int length) {
      if (start < 0)
         throw new IndexOutOfBoundsException("start must not be negative");
      if (length < 0)
         throw new IllegalArgumentException("length is negative");
      if (start + length > m.size())
         throw new IndexOutOfBoundsException("Not enough elements in the matrix");
      double s = 0;
      for (int i = start; i < start + length; i++)
         s += m.getQuick(i);
      return s;
   }

   public static double[] getSum(double[][] a, int startColumn, int numColumns) {
      if (a.length == 0)
         return new double[0];
      if (startColumn < 0)
         throw new IndexOutOfBoundsException("startColumn must not be negative");
      if (numColumns < 0)
         throw new IllegalArgumentException("numColumns is negative");
      if (startColumn + numColumns > a[0].length)
         throw new IndexOutOfBoundsException("Not enough elements in the array");
      double[] res = new double[a.length];
      for (int i = 0; i < res.length; i++)
         for (int j = startColumn; j < startColumn + numColumns; j++)
            res[i] += a[i][j];
      return res;
   }

   public static double[] getSum(DoubleMatrix2D m, int startColumn, int numColumns) {
      if (startColumn < 0)
         throw new IndexOutOfBoundsException("startColumn must not be negative");
      if (numColumns < 0)
         throw new IllegalArgumentException("numColumns is negative");
      if (startColumn + numColumns > m.columns())
         throw new IndexOutOfBoundsException("Not enough columns in the matrix");
      double[] res = new double[m.rows()];
      for (int i = 0; i < res.length; i++)
         for (int j = startColumn; j < startColumn + numColumns; j++)
            res[i] += m.getQuick(i, j);
      return res;
   }

   public static void regroupElements(double[] a, int x) {
      if (x < 1)
         throw new IllegalArgumentException("x < 1");
      int gs = a.length / x;
      for (int i = 0; i < gs; i++) {
         double o = 0;
         for (int j = 0; j < x; j++)
            o += a[i * x + j];
         a[i] = o;
      }
      int m = a.length % x;
      if (m > 0) {
         double o = 0;
         int r = a.length - x * gs;
         for (int j = 0; j < r; j++)
            o += a[gs * x + j];
         a[gs] = o;
         ++gs;
      }
      for (int i = gs; i < a.length; i++)
         a[i] = 0;
   }

   public static void regroupElements(DoubleArrayList l, int x) {
      if (x < 1)
         throw new IllegalArgumentException("x < 1");
      int gs = l.size() / x;
      for (int i = 0; i < gs; i++) {
         double o = 0;
         for (int j = 0; j < x; j++)
            o += l.getQuick(i * x + j);
         l.setQuick(i, o);
      }
      int m = l.size() % x;
      if (m > 0) {
         double o = 0;
         int r = l.size() - x * gs;
         for (int j = 0; j < r; j++)
            o += l.getQuick(gs * x + j);
         l.setQuick(gs, o);
         ++gs;
      }
      l.setSize(gs);
   }

   public static void regroupElements(DoubleMatrix1D mat, int x) {
      if (x < 1)
         throw new IllegalArgumentException("x < 1");
      int gs = mat.size() / x;
      for (int i = 0; i < gs; i++) {
         double o = 0;
         for (int j = 0; j < x; j++)
            o += mat.getQuick(i * x + j);
         mat.setQuick(i, o);
      }
      int m = mat.size() % x;
      if (m > 0) {
         double o = 0;
         int r = mat.size() - x * gs;
         for (int j = 0; j < r; j++)
            o += mat.getQuick(gs * x + j);
         mat.setQuick(gs, o);
         ++gs;
      }
      for (int i = gs; i < mat.size(); i++)
         mat.setQuick(i, 0);
   }

   public static void regroupElements(DoubleMatrix2D mat, int x) {
      if (x < 1)
         throw new IllegalArgumentException("x < 1");
      int gs = mat.columns() / x;
      for (int i = 0; i < gs; i++) {
         for (int r = 0; r < mat.rows(); r++) {
            double o = 0;
            for (int j = 0; j < x; j++)
               o += mat.getQuick(r, i * x + j);
            mat.setQuick(r, i, o);
         }
      }
      int m = mat.columns() % x;
      if (m > 0) {
         int r = mat.columns() - x * gs;
         for (int row = 0; row < mat.rows(); row++) {
            double o = 0;
            for (int j = 0; j < r; j++)
               o += mat.getQuick(row, gs * x + j);
            mat.setQuick(row, gs, o);
         }
         ++gs;
      }
      for (int row = 0; row < mat.rows(); row++)
         for (int i = gs; i < mat.columns(); i++)
            mat.setQuick(row, i, 0);
   }
}