ArrayOfComparableChains.java

package umontreal.ssj.markovchainrqmc;
 
import umontreal.ssj.stat.PgfDataTable;
import umontreal.ssj.stat.Tally;
import umontreal.ssj.util.*;
import umontreal.ssj.util.multidimsort.*;
import umontreal.ssj.hups.*;
import umontreal.ssj.charts.*;
import umontreal.ssj.functionfit.LeastSquares;
import java.io.IOException;
import java.io.Writer;
import java.io.FileWriter;
import java.lang.reflect.Array;

public class ArrayOfComparableChains<T extends MarkovChainComparable> {
   protected T baseChain; // The base chain.
   protected int n; // Current number of chains.
   // protected int stateDim; // Dimension of the chain state.
   protected T[] chains; // Array of n comparable chains.
   protected double[] performances; // Performances for the n chains.
   protected PointSetRandomization randomization;
   protected MultiDimSort<T> savedSort;
   protected int sortCoordPts = 0; // Point coordinates used to sort points.

   public Tally[] performancePerRun;

   public ArrayOfComparableChains(T baseChain) {
      this.baseChain = baseChain;
      // stateDim = baseChain.stateDim;
   }

   public ArrayOfComparableChains(T baseChain, PointSetRandomization rand, MultiDimSort<T> sort) {
      this.baseChain = baseChain;
      // stateDim = baseChain.stateDim;
      randomization = rand;
      savedSort = sort;
   }

   @SuppressWarnings("unchecked")
   public void makeCopies(int n) {
 
      final T[] c = (T[]) Array.newInstance(baseChain.getClass(), n);
      chains = c; // Array of Markov chains.
 
      this.n = n;
      performances = new double[n]; // Array to store the performances.
      for (int i = 0; i < n; i++) {
         try {
            chains[i] = (T) baseChain.clone();
         } catch (CloneNotSupportedException e) {
            System.err.println("ArrayOfComparableChains:");
            e.printStackTrace();
         }
      }
   }

   public void initialStates() {
      int i = 0;
      for (T mc : chains) {
         mc.initialState();
         performances[i] = mc.getPerformance(); // Needed? Why do this?
         ++i;
      }
   }

   public int getN() {
      return n;
   }

   public T[] getChains() {
      return chains;
   }

   public void setRandomization(PointSetRandomization rand) {
      randomization = rand;
   }

   public PointSetRandomization getRandomization() {
      return randomization;
   }

   public void setSort(MultiDimSort<T> sort) {
      savedSort = sort;
   }

   public MultiDimSort<T> getSort() {
      return savedSort;
   }

   public int simulOneStepArrayRQMC(PointSet p, PointSetRandomization rand, MultiDimSort<T> sort, int sortCoordPts) {
      int nStopped = 0;
      p.randomize(rand); // Randomize point set.
      if (sortCoordPts > 0) {
         if (!(p instanceof CachedPointSet))
            throw new IllegalArgumentException("p is not a CachedPointSet.");
         if (sortCoordPts > 1)
            ((CachedPointSet) p).sort(sort); // Sort points using first sortCoordPts coordinates.
         else
            ((CachedPointSet) p).sortByCoordinate(0); // Sort by first coordinate.
      }
      PointSetIterator stream = p.iterator();
      stream.resetCurPointIndex(); // Go to first point.
      int i = 0;
      for (T mc : chains) { // Assume the chains are sorted
         if (mc.hasStopped()) {
            ++nStopped;
         } else {
            stream.setCurCoordIndex(sortCoordPts); // Skip first sortCoordPts coord.
            mc.nextStep(stream); // simulate next step of the chain.
            stream.resetNextSubstream(); // Go to next point.
            if (mc.hasStopped())
               ++nStopped;
         }
         performances[i] = mc.getPerformance();
         ++i;
      }
      return n - nStopped;
   }

   public int simulOneStepArrayRQMC(PointSet p) {
      return simulOneStepArrayRQMC(p, randomization, savedSort, 0);
   }

   public double simulArrayRQMC(PointSet p, PointSetRandomization rand, MultiDimSort<T> sort, int sortCoordPts,
         int numSteps) {
      int numNotStopped = n;
      initialStates();
      int step = 0;
      while (step < numSteps && numNotStopped > 0) {
         if (numNotStopped == n)
            sort.sort(chains, 0, n);
         else
            sortNotStoppedChains(sort); // Sort the numNotStopped first chains.
         p.randomize(rand); // Randomize the point set.
         if (sortCoordPts > 0) {
            if (!(p instanceof CachedPointSet))
               throw new IllegalArgumentException("p is not a CachedPointSet.");
            if (sortCoordPts > 1)
               ((CachedPointSet) p).sort(sort); // Sort points using first sortCoordPts coordinates.
            else
               ((CachedPointSet) p).sortByCoordinate(0); // Sort by first coordinate.
         }
         PointSetIterator stream = p.iterator();
         stream.resetCurPointIndex(); // Go to first point.
         int i = 0;
         for (T mc : chains) { // Assume the chains are sorted
            if (mc.hasStopped()) {
               numNotStopped--;
            } else {
               stream.setCurCoordIndex(sortCoordPts); // Skip first sortCoordPts coord.
               mc.nextStep(stream); // simulate next step of the chain.
               stream.resetNextSubstream(); // Go to next point.
               if (mc.hasStopped())
                  numNotStopped--;
            }
            performances[i] = mc.getPerformance();
            ++i;
         }
         ++step;
      }
      return calcMeanPerf();
   }

   public double simulArrayRQMC(PointSet p, PointSetRandomization rand, MultiDimSort<T> sort, int numSteps) {
      return simulArrayRQMC(p, rand, sort, 0, numSteps);
   }

   public double simulArrayRQMC(PointSet p, int numSteps) {
      return simulArrayRQMC(p, randomization, savedSort, 0, numSteps);
   }

   public double[] getPerformances() {
      return performances;
   }

   public double calcMeanPerf() {
      double sumPerf = 0.0; // Sum of performances.
      for (int i = 0; i < n; ++i) {
         sumPerf += performances[i];
      }
      return sumPerf / n;
   }

   public void simulReplicatesArrayRQMC(PointSet p, PointSetRandomization rand, MultiDimSort<T> sort, int sortCoordPts,
         int numSteps, int m, Tally statReps) {
      makeCopies(p.getNumPoints());
      statReps.init();
      for (int rep = 0; rep < m; rep++) {
         statReps.add(simulArrayRQMC(p, rand, sort, sortCoordPts, numSteps));
      }
   }

   public String simulReplicatesArrayRQMCFormat(PointSet p, PointSetRandomization rand, MultiDimSort<T> sort,
         int sortCoordPts, int numSteps, int m, Tally statReps) {
      Chrono timer = Chrono.createForSingleThread();
      makeCopies(p.getNumPoints());
      timer.init();
      statReps.init();
      for (int rep = 0; rep < m; rep++) {
         statReps.add(simulArrayRQMC(p, rand, sort, sortCoordPts, numSteps));
      }
      StringBuffer sb = new StringBuffer("----------------------------------------------" + PrintfFormat.NEWLINE);
      sb.append("Array-RQMC simulations:" + PrintfFormat.NEWLINE);
      sb.append(PrintfFormat.NEWLINE + p.toString() + ":" + PrintfFormat.NEWLINE);
      sb.append(" Number of indep copies m  = " + m);
      sb.append(PrintfFormat.NEWLINE + " Number of points n        = " + n + PrintfFormat.NEWLINE);
      sb.append(baseChain.formatResultsRQMC(statReps, n));
      sb.append(" CPU Time = " + timer.format() + PrintfFormat.NEWLINE);
      return sb.toString();
   }

   public String varianceImprovementFormat(double varRQMC, double varMC) {
      // double varRQMC = p.getNumPoints() * statReps.variance();
      StringBuffer sb = new StringBuffer(
            " Variance ratio MC / RQMC: " + PrintfFormat.format(15, 10, 4, varMC / varRQMC) + PrintfFormat.NEWLINE);
      return sb.toString();
   }

   public String testVarianceRateFormat(PointSet[] pointSets, PointSetRandomization rand, MultiDimSort<T> sort,
         int sortCoordPts, int numSteps, int m, double varMC, String filenamePlot, String methodLabel) {
 
      int numSets = pointSets.length; // Number of point sets.
      Tally statPerf = new Tally("Performance");
      double[] logn = new double[numSets];
      double[] variance = new double[numSets];
      double[] logVariance = new double[numSets];
      long initTime; // For timings.
 
      StringBuffer str = new StringBuffer("\n\n --------------------------");
      str.append(methodLabel + "\n  MC Variance : " + varMC + "\n\n");
 
      // Array-RQMC experiment with each pointSet.
      for (int i = 0; i < numSets; ++i) {
         initTime = System.currentTimeMillis();
         n = pointSets[i].getNumPoints();
         str.append("n = " + n + "\n");
         simulReplicatesArrayRQMC(pointSets[i], rand, sort, sortCoordPts, numSteps, m, statPerf);
         logn[i] = Num.log2(n);
         variance[i] = statPerf.variance();
         logVariance[i] = Num.log2(variance[i]);
         str.append("  Average = " + statPerf.average() + "\n");
         str.append("  VRF =  " + varMC / (n * variance[i]) + "\n");
         str.append(formatTime((System.currentTimeMillis() - initTime) / 1000.) + "\n");
      }
      // Estimate regression slope and print plot and overall results.
      double regSlope = slope(logn, logVariance, numSets);
      str.append("Regression slope (log) for variance = " + regSlope + "\n\n");
 
      // Print plot and overall results in files.
      if (filenamePlot != null)
         try {
            Writer file = new FileWriter(filenamePlot + ".tex");
            XYLineChart chart = new XYLineChart();
            // ("title", "$log_2(n)$", "$log_2 Var[hat mu_{rqmc,s,n}]$");
            chart.add(logn, logVariance);
            file.write(chart.toLatex(12, 8));
            file.close();
         } catch (IOException e) {
            e.printStackTrace();
         }
      return str.toString();
   }

   public String testVarianceRateFormat(RQMCPointSet[] rqmcPts, MultiDimSort<T> sort, int sortCoordPts, int numSteps,
         int m, double varMC, String filenamePlot, String methodLabel) {
      int numSets = rqmcPts.length; // Number of point sets.
      Tally statPerf = new Tally("Performance");
      double[] logn = new double[numSets];
      double[] variance = new double[numSets];
      double[] logVariance = new double[numSets];
      long initTime; // For timings.
 
      StringBuffer str = new StringBuffer("\n\n --------------------------");
      str.append(methodLabel + "\n  MC Variance per Run : " + varMC / (double) n + "\n\n");
 
      // Array-RQMC experiment with each pointSet.
      for (int i = 0; i < numSets; ++i) {
 
//          performancePerRun = new Tally[numSteps]; 
         for (int j = 0; j < numSteps; ++j)
            performancePerRun[j] = new Tally();
         initTime = System.currentTimeMillis();
         n = rqmcPts[i].getNumPoints();
         str.append("n = " + n + "\n");
         simulReplicatesArrayRQMC(rqmcPts[i].getPointSet(), rqmcPts[i].getRandomization(), sort, sortCoordPts, numSteps,
               m, statPerf);
         logn[i] = Num.log2(n);
         variance[i] = statPerf.variance();
         logVariance[i] = Num.log2(variance[i]);
         str.append("  Average = " + statPerf.average() + "\n");
         str.append(" RQMC Variance : " + variance[i] + "\n\n");
         str.append("  VRF =  " + varMC / (n * variance[i]) + "\n");
         str.append(formatTime((System.currentTimeMillis() - initTime) / 1000.) + "\n");
      }
      // Estimate regression slope and print plot and overall results.
      double regSlope = slope(logn, logVariance, numSets);
      str.append("Regression slope (log) for variance = " + regSlope + "\n\n");
 
      String[] tableField = { "log(n)", "log(Var)" };
      double[][] data = new double[numSets][2];
      for (int s = 0; s < numSets; s++) { // For each cardinality n
         data[s][0] = logn[s];
         data[s][1] = logVariance[s];
      }
 
      // Print plot and overall results in files.
      if (filenamePlot != null)
         try {
 
            PgfDataTable pgf = new PgfDataTable(filenamePlot, rqmcPts[0].getLabel(), tableField, data);
            String pVar = pgf.drawPgfPlotSingleCurve(filenamePlot, "axis", 0, 1, 2, "", "");
            String plotIV = (PgfDataTable.pgfplotFileHeader() + pVar + PgfDataTable.pgfplotEndDocument());
 
            FileWriter fileIV = new FileWriter(filenamePlot + "_" + "VAr.tex");
            fileIV.write(plotIV);
            fileIV.close();
         } catch (IOException e) {
            e.printStackTrace();
         }
 
      tableField = new String[] { "step", "average", "variance" };
      data = new double[numSteps][3];
      for (int s = 0; s < numSteps; s++) { // For each cardinality n
         data[s][0] = s + 1;
         data[s][1] = performancePerRun[s].average();
         data[s][2] = performancePerRun[s].variance();
      }
 
      if (filenamePlot != null)
         try {
 
            PgfDataTable pgf = new PgfDataTable(filenamePlot, rqmcPts[0].getLabel(), tableField, data);
            String pMean = pgf.drawPgfPlotSingleCurve(filenamePlot, "axis", 0, 1, 2, "", "");
            String pVar = pgf.drawPgfPlotSingleCurve(filenamePlot, "axis", 0, 2, 2, "", "");
 
            String plotMean = (PgfDataTable.pgfplotFileHeader() + pMean + PgfDataTable.pgfplotEndDocument());
            String plotVar = (PgfDataTable.pgfplotFileHeader() + pVar + PgfDataTable.pgfplotEndDocument());
 
            FileWriter fileIV = new FileWriter(filenamePlot + "_" + "MeanPerStep.tex");
            fileIV.write(plotMean);
            fileIV.close();
 
            fileIV = new FileWriter(filenamePlot + "_" + "VariancePerStep.tex");
            fileIV.write(plotVar);
            fileIV.close();
         } catch (IOException e) {
            e.printStackTrace();
         }
 
      return str.toString();
   }

   public void sortNotStoppedChains(MultiDimSort<T> sort) {
      int j = n - 1;
      int i = 0;
      T mc;
      while (j >= 0 && chains[j].hasStopped())
         --j;
      while (i < n && !chains[i].hasStopped())
         ++i;
      while (i < j) {
         while (!chains[i].hasStopped())
            ++i;
         while (chains[j].hasStopped())
            --j;
         mc = chains[i];
         chains[i] = chains[j];
         chains[j] = mc;
      }
      sort.sort(chains, 0, i);
   }

   public void sortChains() {
      savedSort.sort(chains, 0, n);
   }
 
   // Takes time in seconds and formats it.
   public String formatTime(double time) {
      int second, hour, min, centieme;
      hour = (int) (time / 3600.0);
      if (hour > 0) {
         time -= ((double) hour * 3600.0);
      }
      min = (int) (time / 60.0);
      if (min > 0) {
         time -= ((double) min * 60.0);
      }
      second = (int) time;
      centieme = (int) (100.0 * (time - (double) second) + 0.5);
      return String.valueOf(hour) + ":" + min + ":" + second + "." + centieme;
   }
 
   // Compute slope of linear regression of y on x, using only first n
   // observations.
   public double slope(double[] x, double[] y, int n) {
      if (n < 2) {
         return 0.0;
      } else {
         double[] x2 = new double[n], y2 = new double[n];
         for (int i = 0; i < n; ++i) {
            x2[i] = x[i];
            y2[i] = y[i];
         }
         return LeastSquares.calcCoefficients(x2, y2, 1)[1];
      }
   }
}