Misc.java

/*
 * Class:        Misc
 * Description:  Miscellaneous functions that are hard to classify.
 * 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.util;

public class Misc {
   private Misc() {
   }

   public static double quickSelect(double[] A, int n, int k) {
      double[] U = new double[n];
      double[] V = new double[n];
      double p = A[k - 1];
      int u = 0;
      int v = 0;
      int indV = 0;
 
      for (int i = 0; i < n; i++) {
         if (A[i] <= p) {
            v++;
            if (A[i] != p) {
               U[u++] = A[i];
            }
         } else
            V[indV++] = A[i];
      }
 
      if (k <= u)
         return quickSelect(U, u, k);
      else if (k > v)
         return quickSelect(V, indV, k - v);
      else
         return p;
   }

   public static int quickSelect(int[] A, int n, int k) {
      int[] U = new int[n];
      int[] V = new int[n];
      int p = A[k - 1];
      int u = 0;
      int v = 0;
      int indV = 0;
 
      for (int i = 0; i < n; i++) {
         if (A[i] <= p) {
            v++;
            if (A[i] != p) {
               U[u++] = A[i];
            }
         } else
            V[indV++] = A[i];
      }
 
      if (k <= u)
         return quickSelect(U, u, k);
      else if (k > v)
         return quickSelect(V, indV, k - v);
      else
         return p;
   }

   public static double getMedian(double[] A, int n) {
      int k = (n + 1) / 2; // median index
      double med = quickSelect(A, n, k);
      double y;
      if ((n & 1) == 0) {
         y = quickSelect(A, n, k + 1);
         med = (med + y) / 2.0;
      }
      return med;
   }

   public static double getMedian(int[] A, int n) {
      int k = (n + 1) / 2; // median index
      double med = quickSelect(A, n, k);
      double y;
      if ((n & 1) == 0) {
         y = quickSelect(A, n, k + 1);
         med = (med + y) / 2.0;
      }
      return med;
   }

   public static int getTimeInterval(double[] times, int start, int end, double t) {
      if (start < 0)
         throw new IllegalArgumentException("The starting index must not be negative");
      int n = end - start;
      if (n < 0)
         throw new IllegalArgumentException("The ending index must be greater than or equal to the starting index");
      if (t < times[start])
         return -1;
      if (t >= times[end])
         return n;
 
      int start0 = start;
      // Perform binary search to find the interval index
      int mid = (start + end) / 2;
      // Test if t is inside the interval mid.
      // The interval mid starts at times[mid],
      // and the interval mid+1 starts at times[mid + 1].
      while (t < times[mid] || t >= times[mid + 1]) {
         if (start == end)
            // Should not happen, safety check to avoid infinite loops.
            throw new IllegalStateException();
         if (t < times[mid])
            // time corresponds to an interval before mid.
            end = mid - 1;
         else
            // time corresponds to an interval after mid.
            start = mid + 1;
         mid = (start + end) / 2;
      }
      return mid - start0;
   }

   public static void interpol(int n, double[] X, double[] Y, double[] C) {
      int j;
      // Compute divided differences for the Newton interpolating polynomial
      for (j = 0; j <= n; ++j)
         C[j] = Y[j];
      for (int i = 1; i <= n; ++i)
         for (j = n; j >= i; --j) {
            if (X[j] == X[j - i])
               C[j] = 0;
            else
               C[j] = (C[j] - C[j - 1]) / (X[j] - X[j - i]);
         }
   }

   public static double evalPoly(int n, double[] X, double[] C, double z) {
      double v = C[n];
      for (int j = n - 1; j >= 0; --j)
         v = v * (z - X[j]) + C[j];
      return v;
   }

   public static double evalPoly(double[] C, int n, double x) {
      double v = C[n];
      for (int j = n - 1; j >= 0; --j)
         v = v * x + C[j];
      return v;
   }
 
}