RQMCExperimentSeries.java

/*
 * Class:        RQMCExperimentSeries
 * Description:  randomized quasi-Monte Carlo simulations
 * 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.mcqmctools;
 
import umontreal.ssj.functionfit.LeastSquares;
import umontreal.ssj.hups.*;
import umontreal.ssj.stat.PgfDataTable;
import umontreal.ssj.stat.Tally;
import umontreal.ssj.stat.list.lincv.ListOfTalliesWithCV;
import umontreal.ssj.util.Chrono;
import umontreal.ssj.util.PrintfFormat;
import java.util.ArrayList;

 
public class RQMCExperimentSeries {
   int numSets = 0; // Number of point sets in the series.
   RQMCPointSet[] theSets;
   double base = 2.0; // Base for the logs (in base 2 by default)
   double logOfBase; // Math.log(base)
   double[] size = new double[numSets]; // values of n
   double[] mean = new double[numSets]; // average performance for each point set
   double[] variance = new double[numSets]; // variance of the average for each point set
   double[] logn = new double[numSets]; // log_base n
   double[] logVar = new double[numSets]; // log_base (variance)
   String[] tableFields = { "n", "mean", "variance", "log(n)", "log(variance)" };
   // Names of fields for table.
   boolean displayExec = false; // When true, prints a display of execution in real time
   int numReplicates; // last value of m
   MonteCarloModelDouble model;
   // int numSkipRegression = 0; // Number of values of n that are skipped for the
   // regression
   String cpuTime; // time for last experiment\
   String title;

   public RQMCExperimentSeries() {};

   public RQMCExperimentSeries(RQMCPointSet[] theSets, double base) {
      init(theSets, base);
   }

   public void init(RQMCPointSet[] theSets, double base) {
      this.base = base;
      this.logOfBase = Math.log(base);
      this.numSets = theSets.length;
      this.theSets = theSets;
      size = new double[numSets]; // n for each point set
      mean = new double[numSets]; // average for each point set
      variance = new double[numSets]; // variance for each point set
      logn = new double[numSets]; // log n
      logVar = new double[numSets]; // log (variance)
   }

   public void setExecutionDisplay(boolean display) {
      this.displayExec = display;
   }

   public void setBase(double b) {
      base = b;
      logOfBase = Math.log(base);
   }

   public double getBase() {
      return base;
   }

   public RQMCPointSet getSet(int i) {
      return theSets[i];
   }

   public double[] getValuesn() {
      return size;
   }

   public double[] getLogn() {
      return logn;
   }

   public double[] getMeans() {
      return mean;
   }

   public double[] getVariances() {
      return variance;
   }

   public double[] getLogVariances() {
      return logVar;
   }

   public void testVarianceRate(MonteCarloModelDouble model, int m) {
      int n;
      Tally statReps = new Tally();
      Chrono timer = new Chrono();
      numReplicates = m;
      this.model = model;
      if (displayExec) {
         System.out.println("\n ============================================= ");
         System.out.println("RQMC simulation for mean estimation:  ");
         System.out.println("Model: " + model.toString());
         System.out.println(" Number of indep copies m  = " + m);
         System.out.println(" Point sets: " + theSets[0].toString() + "\n");
         System.out.println("    n     CPU time         mean      log(var) ");
      }
      for (int s = 0; s < numSets; s++) { // For each cardinality n
         n = theSets[s].getNumPoints();
         size[s] = n;
         logn[s] = Math.log(n) / logOfBase;
         // System.out.println(" n = " + n + ", log n = " + logn[s] + "\n"); // ****
         // System.out.println(" " + n + " " + timer.format());
         RQMCExperiment.simulReplicatesRQMC(model, theSets[s], m, statReps);
         mean[s] = statReps.average();
         variance[s] = statReps.variance();
         logVar[s] = Math.log(variance[s]) / logOfBase;
         if (displayExec) {
            System.out.println("  " + n + "     " + timer.format() + "   " + PrintfFormat.f(10, 5, mean[s]) + "   "
                  + PrintfFormat.f(7, 2, logVar[s]));
         }
      }
      cpuTime = timer.format();
   }

   public void testVarianceRateCV(MonteCarloModelCV model, int m) {
      int numCV = model.getNumberCV();
      ListOfTalliesWithCV<Tally> statWithCV = ListOfTalliesWithCV.createWithTally(1, numCV);
      statWithCV.setExpectedValue(0, 0.0); // The CV is centered to 0.
      int n;
      Chrono timer = new Chrono();
      numReplicates = m;
      this.model = model;
      if (displayExec) {
         System.out.println("\n ============================================= ");
         System.out.println("RQMC simulation for mean estimation with control variates:  ");
         System.out.println("Model: " + model.toString());
         System.out.println(" Number of indep copies m  = " + m);
         System.out.println(" Point sets: " + theSets[0].toString() + "\n");
         System.out.println("    n     CPU time         mean     log(var) ");
      }
      for (int s = 0; s < numSets; s++) { // For each cardinality n
         n = theSets[s].getNumPoints();
         size[s] = n;
         logn[s] = Math.log(n) / logOfBase;
         // System.out.println(" n = " + n + ", log n = " + logn[s] + "\n"); // ****
         // System.out.println(" " + n + " " + timer.format());
         RQMCExperiment.simulReplicatesRQMCCV(model, theSets[s], m, statWithCV);
         statWithCV.estimateBeta(); // This is where the var. and covar. are computed!
         mean[s] = statWithCV.averageWithCV(0);
         variance[s] = statWithCV.covarianceWithCV(0, 0);
         logVar[s] = Math.log(variance[s]) / logOfBase;
         if (displayExec) {
            System.out.println("  " + n + "     " + timer.format() + "   " + PrintfFormat.f(10, 5, mean[s]) + "   "
                  + PrintfFormat.f(7, 2, logVar[s]));
         }
      }
      cpuTime = timer.format();
   }

   public double[] regressionLogVariance(int numSkip) {
      double[] x2 = new double[numSets - numSkip], y2 = new double[numSets - numSkip];
      for (int i = 0; i < numSets - numSkip; ++i) {
         x2[i] = logn[i + numSkip];
         y2[i] = logVar[i + numSkip];
      }
      return LeastSquares.calcCoefficients(x2, y2, 1);
   }

   public String formatRegression(double[] regCoeff) {
      StringBuffer sb = new StringBuffer("");
      // double[] regCoeff = regressionLogVariance (numSkipRegression);
      sb.append("  Slope of log(var) = " + PrintfFormat.f(8, 5, regCoeff[1]) + "\n");
      sb.append("    constant term      = " + PrintfFormat.f(8, 5, regCoeff[0]) + "\n\n");
      return sb.toString();
   }

   public String reportVarianceRate(int numSkip, boolean details) {
      StringBuffer sb = new StringBuffer("");
      sb.append("\n ============================================= \n");
      sb.append("RQMC simulation for mean estimation: \n ");
      sb.append("Model: " + model.toString() + "\n");
      sb.append(" Number of indep copies m  = " + numReplicates + "\n");
      sb.append(" RQMC point sets: " + theSets[0].toString() + "\n\n");
      sb.append("RQMC variance \n");
      if (details)
         sb.append(dataLogForPlot());
      sb.append(formatRegression(regressionLogVariance(numSkip)));
      // sb.append(" Slope of log(var) = " + PrintfFormat.f(8, 5, regCoeff[1]) +
      // "\n");
      // sb.append(" constant term = " + PrintfFormat.f(8, 5, regCoeff[0]) + "\n\n");
      sb.append("  Total CPU Time = " + cpuTime + "\n");
      sb.append("-----------------------------------------------------\n");
      return sb.toString();
   }

   public PgfDataTable toPgfDataTable(String tableName, String tableLabel) {
      double[][] data = new double[numSets][5];
      for (int s = 0; s < numSets; s++) { // For each cardinality n
         data[s][0] = size[s];
         data[s][1] = mean[s];
         data[s][2] = variance[s];
         data[s][3] = logn[s];
         data[s][4] = logVar[s];
      }
      return new PgfDataTable(tableName, tableLabel, tableFields, data);
   }
 
   public PgfDataTable toPgfDataTable(String tableLabel) {
      return toPgfDataTable(title, tableLabel);
   }

   public String dataForPlot() {
      StringBuffer sb = new StringBuffer("");
      sb.append("    n      mean       variance \n");
      for (int s = 0; s < numSets; s++) // For each cardinality n
         sb.append(
               " " + size[s] + " " + PrintfFormat.f(10, 5, mean[s]) + " " + PrintfFormat.f(10, 5, variance[s]) + "\n");
      return sb.toString();
   }

   public String dataLogForPlot() {
      StringBuffer sb = new StringBuffer("");
      sb.append("   log(n)      mean       log(variance) \n");
      for (int s = 0; s < numSets; s++) // For each cardinality n
         sb.append(
               " " + logn[s] + " " + PrintfFormat.f(10, 5, mean[s]) + " " + PrintfFormat.f(10, 5, logVar[s]) + "\n");
      return sb.toString();
   }
 
   /*
    * Returns the data on the mean and variance for each n, in an appropriate
    * format to produce a plot with the pgfplot package. This is OBSOLETE.
    * 
    * @return Report as a string.
    * 
    * public String XformatPgfCurve (String curveName) { StringBuffer sb = new
    * StringBuffer("");
    * sb.append("      \\addplot+[no marks] table[x=n,y=variance] {" + "\n");
    * sb.append( dataForPlot() + " } \n"); sb.append("      \\addlegendentry{" +
    * curveName + "}\n"); sb.append("      % \n"); return sb.toString(); }
    */

   public String testVarianceRateManyPointTypes(MonteCarloModelDouble model, ArrayList<RQMCPointSet[]> list, int m,
         int numSkip, boolean makePgfTable, boolean printReport, boolean details, ArrayList<PgfDataTable> listCurves) {
      StringBuffer sb = new StringBuffer("");
      // if (makePgfTable)
      // listCurves = new ArrayList<PgfDataTable>();
      for (RQMCPointSet[] ptSeries : list) {
         init(ptSeries, base);
         testVarianceRate(model, m);
         if (printReport)
            System.out.println(reportVarianceRate(numSkip, details));
         if (printReport)
            sb.append(reportVarianceRate(numSkip, details));
         if (makePgfTable == true)
            listCurves.add(toPgfDataTable(ptSeries[0].getLabel()));
      }
      return sb.toString();
   }

   public String toString() {
      return ("RQMC Experiment:" + title); // theSets[0].toString();
   }
}