FunctionOfMultipleMeansTallyWithCV.java

/*
 * Class:        FunctionOfMultipleMeansTallyWithCV
 * Description:  function of multiple means tally with control variables
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Éric Buist
 * @since        August 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.stat.list.lincv;
 
import cern.colt.matrix.DoubleMatrix2D;
 
import umontreal.ssj.util.MultivariateFunction;
import umontreal.ssj.stat.FunctionOfMultipleMeansTally;
import umontreal.ssj.stat.Tally;

public class FunctionOfMultipleMeansTallyWithCV extends FunctionOfMultipleMeansTally {
   private MultivariateFunction funcNoCV;
   private double[] beta;

   public FunctionOfMultipleMeansTallyWithCV(MultivariateFunction funcNoCV, int p, int q) {
      // Assumes that the superclass' constructor do not evaluate the function
      super(new LinCVFunction(p + q), ListOfTalliesWithCV.createWithTally(p, q));
      this.funcNoCV = funcNoCV;
      beta = new double[q];
      ((LinCVFunction) getFunction()).initFunctionOfMultipleMeansTallyWithCV(this);
   }

   public FunctionOfMultipleMeansTallyWithCV(MultivariateFunction funcNoCV, ListOfTalliesWithCV<Tally> l) {
      // Assumes that the superclass' constructor do not evaluate the function
      super(new LinCVFunction(l.size()), l);
      this.funcNoCV = funcNoCV;
      beta = new double[l.getNumControlVariables()];
      ((LinCVFunction) getFunction()).initFunctionOfMultipleMeansTallyWithCV(this);
   }

   public MultivariateFunction getFunctionWithoutCV() {
      return funcNoCV;
   }

   public int getNumControlVariables() {
      return beta.length;
   }

   public int getDimensionWithoutCV() {
      return getDimension() - beta.length;
   }

   public double getBeta(int i) {
      return beta[i];
   }

   public void setBeta(int i, double b) {
      beta[i] = b;
   }

   public double[] getBeta() {
      return beta;
   }

   public void setBeta(double[] beta) {
      if (beta.length != this.beta.length)
         throw new IllegalArgumentException("Invalid length of beta");
      this.beta = beta;
   }

   public ListOfTalliesWithCV<Tally> getListOfTalliesWithCV() {
      return (ListOfTalliesWithCV<Tally>) super.getListOfTallies();
   }

   public double getExpectedValue(int i) {
      return getListOfTalliesWithCV().getExpectedValue(i);
   }

   public void setExpectedValue(int i, double e) {
      getListOfTalliesWithCV().setExpectedValue(i, e);
   }

   public double[] getExpectedValues() {
      return getListOfTalliesWithCV().getExpectedValues();
   }

   public void setExpectedValues(double[] exp) {
      getListOfTalliesWithCV().setExpectedValues(exp);
   }

   public void estimateBeta() {
      ListOfTalliesWithCV a = getListOfTalliesWithCV();
      a.estimateBeta();
      estimateBetaFromMatrix(a.getBeta());
   }

   public void estimateBetaFromMatrix(DoubleMatrix2D mbeta) {
      ListOfTalliesWithCV a = getListOfTalliesWithCV();
      int p = getDimension() - beta.length;
      int q = beta.length;
      double[] avg = new double[p];
      double[] gradient = new double[p];
      for (int i = 0; i < p; i++)
         avg[i] = a.get(i).average();
      for (int i = 0; i < p; i++)
         gradient[i] = funcNoCV.evaluateGradient(i, avg);
      for (int i = 0; i < q; i++) {
         beta[i] = 0;
         for (int j = 0; j < p; j++)
            beta[i] += mbeta.getQuick(i, j) * gradient[j];
      }
   }

   public FunctionOfMultipleMeansTallyWithCV clone() {
      FunctionOfMultipleMeansTallyWithCV mta = (FunctionOfMultipleMeansTallyWithCV) super.clone();
      mta.beta = beta.clone();
 
      LinCVFunction fct = new LinCVFunction(getDimension());
      mta.func = fct;
      fct.initFunctionOfMultipleMeansTallyWithCV(mta);
 
      return mta;
   }
 
   private static class LinCVFunction implements MultivariateFunction {
      private FunctionOfMultipleMeansTallyWithCV fcv;
      private int pplusq;
      private double[] tmp;
 
      public LinCVFunction(int pplusq) {
         this.pplusq = pplusq;
      }
 
      public void initFunctionOfMultipleMeansTallyWithCV(FunctionOfMultipleMeansTallyWithCV fcv) {
         this.fcv = fcv;
         tmp = new double[fcv.getListOfTalliesWithCV().sizeWithoutCV()];
      }
 
      public int getDimension() {
         return pplusq;
      }
 
      public double evaluate(double... x) {
         MultivariateFunction funcNoCV = fcv.getFunctionWithoutCV();
 
         if (x.length != getDimension())
            throw new IllegalArgumentException(
                  "x has length " + x.length + ", which differs from the dimension " + getDimension());
         int d = getDimension() - fcv.beta.length;
         System.arraycopy(x, 0, tmp, 0, d);
         double gValue = funcNoCV.evaluate(tmp);
         ListOfTalliesWithCV<?> a = fcv.getListOfTalliesWithCV();
         for (int i = 0; i < fcv.beta.length; i++)
            gValue -= fcv.beta[i] * (x[d + i] - a.getExpectedValue(i));
         return gValue;
      }
 
      public double evaluateGradient(int i, double... x) {
         MultivariateFunction funcNoCV = fcv.getFunctionWithoutCV();
 
         if (x.length != getDimension())
            throw new IllegalArgumentException(
                  "x has length " + x.length + ", which differs from the dimension " + getDimension());
         int d = getDimension() - fcv.beta.length;
         if (i < d) {
            System.arraycopy(x, 0, tmp, 0, d);
            return funcNoCV.evaluateGradient(i, tmp);
         }
         return -fcv.beta[i - d];
      }
   }
}