Discrepancy.java

/*
 * Class:        Discrepancy
 * Description:  Base class of all discrepancies
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Richard Simard
 * @since        January 2009
 
 * 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.discrepancy;
 
import umontreal.ssj.hups.*;
import cern.colt.list.DoubleArrayList;
import umontreal.ssj.util.PrintfFormat;

public abstract class Discrepancy {
   protected static double[] ONES = { 1 }; // vector of 1, dimension dim
   protected double[] gamma; // gamma[]: weights, dimension dim
   protected double[][] Points; // Points[n][s]: n points in s dimensions
   protected int dim; // Dimension of points
   protected int numPoints; // Number of points
//   protected boolean primeF;                  // true if numPoints is prime
//   protected boolean power2F;                 // true if numPoints is a power of 2
 
   static final double UNSIX = 1.0 / 6.0;
   static final double QUARAN = 1.0 / 42.0;
   static final double UNTRENTE = 1.0 / 30.0;
   static final double DTIERS = 2.0 / 3.0;
   static final double STIERS = 7.0 / 3.0;
   static final double QTIERS = 14.0 / 3.0;
 
   static protected void setONES(int s) {
      if (s < ONES.length)
         return;
      ONES = new double[s + 1];
      for (int i = 0; i <= s; i++)
         ONES[i] = 1.0;
   }
 
   protected void appendGamma(StringBuffer sb, double[] gamma, int s) {
      // append the first s components of gamma to sb
      if (gamma != null) {
         for (int j = 0; j < s; ++j)
            sb.append("  " + gamma[j]);
      }
   }
 
   private void initD(double[][] points, int n, int s, double[] gamma) {
      setONES(s);
      if (gamma == null)
         setGamma(ONES, s);
      else
         setGamma(gamma, s);
      dim = s;
      numPoints = n;
      this.Points = points;
   }

   public Discrepancy(double[][] points, int n, int s) {
      initD(points, n, s, null);
   }

   public Discrepancy(double[][] points, int n, int s, double[] gamma) {
      initD(points, n, s, gamma);
   }

   public Discrepancy(int n, int s, double[] gamma) {
      initD(null, n, s, gamma);
   }

   public Discrepancy(PointSet set) {
      numPoints = set.getNumPoints();
      dim = set.getDimension();
      Points = toArray(set);
      initD(Points, numPoints, dim, null);
   }

   public Discrepancy() {
   }

   public double compute() {
      return compute(Points, numPoints, dim, gamma);
   }

   public double compute(int s) {
      return compute(Points, numPoints, s, gamma);
   }

   public double compute(double[][] points, int n, int s, double[] gamma) {
      return -1;
   }

   public abstract double compute(double[][] points, int n, int s);

   public double compute(double[][] points) {
      return compute(points, points.length, points[0].length, gamma);
   }

   public double compute(double[] T, int n) {
      double[][] p = new double[n][1];
      for (int i = 0; i < n; ++i)
         p[i][0] = T[i];
      return compute(p, n, 1);
   }

   public double compute(double[] T) {
      return compute(T, T.length, gamma[0]);
   }

   public double compute(double[] T, int n, double gamma) {
      return -1;
   }

   public double compute(PointSet set, double[] gamma) {
      int n = set.getNumPoints();
      int dim = set.getDimension();
      if (dim > 1) {
         double[][] points = new double[n][dim];
         PointSetIterator it = set.iterator();
         for (int i = 0; i < n; ++i) {
            it.nextPoint(points[i], dim);
         }
         return compute(points, n, dim, gamma);
      } else {
         double[] points = new double[n];
         PointSetIterator it = set.iterator();
         for (int i = 0; i < n; ++i) {
            points[i] = it.nextCoordinate();
            it.resetToNextPoint();
         }
         return compute(points, n, gamma[0]);
      }
   }

   public double compute(PointSet set) {
      int dim = set.getDimension();
      setONES(dim);
      return compute(set, ONES);
   }

   public int getNumPoints() {
      return numPoints;
   }

   public int getDimension() {
      return dim;
   }

   public void setPoints(double[][] points, int n, int s) {
      initD(points, n, s, null);
   }

   public void setPoints(double[][] points) {
      setPoints(points, points.length, points[0].length);
   }

   public void setGamma(double[] gam, int s) {
      gamma = gam;
   }

   public double[] getGamma() {
      return gamma;
   }

   public static double[][] toArray(PointSet set) {
      int n = set.getNumPoints();
      int dim = set.getDimension();
      if (dim > 1) {
         double[][] points = new double[n][dim];
         PointSetIterator it = set.iterator();
         for (int i = 0; i < n; ++i)
            it.nextPoint(points[i], dim);
         return points;
 
      } else {
         double[][] po1 = new double[n][1];
         PointSetIterator it = set.iterator();
         for (int i = 0; i < n; ++i) {
            po1[i][0] = it.nextCoordinate();
            it.resetToNextPoint();
         }
         return po1;
      }
   }

   public static DoubleArrayList sort(double[] T, int n) {
      DoubleArrayList sortedU = new DoubleArrayList(T);
      sortedU.sortFromTo(0, n - 1);
      return sortedU;
   }

   public String toString() {
      StringBuffer sb = new StringBuffer(getName() + ":" + PrintfFormat.NEWLINE);
      sb.append("n = " + numPoints + ",   dim = " + dim + PrintfFormat.NEWLINE);
      sb.append("gamma = [");
      appendGamma(sb, gamma, dim);
      sb.append(" ]" + PrintfFormat.NEWLINE);
      return sb.toString();
   }

   public String formatPoints() {
      StringBuffer sb = new StringBuffer("Points = [" + PrintfFormat.NEWLINE);
      if (Points != null) {
         for (int i = 0; i < Points.length; ++i) {
            sb.append(" [ " + Points[i][0]);
            for (int j = 1; j < Points[0].length; ++j)
               sb.append(", " + Points[i][j]);
            sb.append(" ]" + PrintfFormat.NEWLINE);
         }
      }
      sb.append(" ]" + PrintfFormat.NEWLINE);
      return sb.toString();
   }

   public String getName() {
      return getClass().getSimpleName();
   }
 
}