package fig.basic;

import java.util.*;

public class NumUtils {
  // This random stuff should be deprecated.  DON'T USE IT!
  @Deprecated
  private static Random random = new Random();
  @Deprecated
  public static Random getRandom() { return random; }
  @Deprecated
  public static void setRandom(long seed) { setRandom(new Random(seed)); }
  @Deprecated
  public static void setRandom(Random random) { NumUtils.random = random; }
  @Deprecated
  public static double randDouble() { return random.nextDouble(); }
  @Deprecated
  public static int randInt(int n) { return random.nextInt(n); }
  @Deprecated
  public static boolean randBernoulli(double p) { return random.nextDouble() < p; }
  @Deprecated
  public static int randMultinomial(double[] probs, Random random) {
    double v = random.nextDouble();
    double sum = 0;
    for(int i = 0; i < probs.length; i++) {
      sum += probs[i];
      if(v < sum) return i;
    }
    throw new RuntimeException(sum + " < " + v);
  }
  @Deprecated
  public static int randMultinomial(double[] probs) {
    return randMultinomial(probs, random);
  }

  public static boolean isFinite(double x) { return !Double.isNaN(x) && !Double.isInfinite(x); }
  public static void assertIsFinite(double x) { assert isFinite(x): "Not finite: " + x; }
  public static void assertIsFinite(double[] xs) { for(double x : xs) assert isFinite(x): "Not finite: " + Fmt.D(xs); }
  public static void assertIsFinite(double[][] xss) { for(double[] xs : xss) for(double x : xs) assert isFinite(x): "Not finite: " + Fmt.D(xss); }
  public static boolean isProb(double x) { return x >= 0 && x <= 1 && !Double.isNaN(x); }
  public static void assertIsProb(double x) { assert isProb(x): "Not a probability [0, 1]: " + x; }
  public static void assertEquals(double x, double y) { assertEquals(x, y, 1e-10); }
  public static void assertEquals(double x, double y, double tol) { assert Math.abs(x-y)<tol: x + " != " + y; }
  public static void assertNormalized(double[] p) { assertEquals(ListUtils.sum(p), 1); }
  public static void assertNormalized(double[] p, double tol) { assertEquals(ListUtils.sum(p), 1, tol); }

  // Vector, matrix operations {
  public static boolean normalize(float[] data) {
    float sum = 0;
    for(float x : data) sum += x;
    if(sum == 0) return false;
    for(int i = 0; i < data.length; i++) data[i] /= sum;
    return true;
  }
  public static boolean normalize(double[] data) {
    double sum = 0;
    for(double x : data) sum += x;
    if(sum == 0) return false;
    for(int i = 0; i < data.length; i++) data[i] /= sum;
    return true;
  }
  public static boolean normalize(double[][] data) {
    double sum = 0;
    for(double[] v : data)
      for(double x : v) sum += x;
    if(sum == 0) return false;
    for(double[] v : data)
      for(int i = 0; i < v.length; i++) v[i] /= sum;
    return true;
  }
  public static boolean normalizeEachRow(double[][] data) {
    boolean allRowsOkay = true;
    for (double[] row: data) {
      if(!NumUtils.normalize(row)) allRowsOkay = false;
    }
    return allRowsOkay;
  }
  public static boolean normalize(double[][][] data) {
    double sum = 0;
    for(double[][] m : data)
      for(double[] v : m)
        for(double x : v) sum += x;
    if(sum == 0) return false;
    for(double[][] m : data)
      for(double[] v : m)
        for(int i = 0; i < v.length; i++) v[i] /= sum;
    return true;
  }

  public static boolean expNormalize(double[] probs) {
    // Input: log probabilities (unnormalized too)
    // Output: normalized probabilities
    // probs actually contains log probabilities; so we can add an arbitrary constant to make
    // the largest log prob 0 to prevent overflow problems
    double max = Double.NEGATIVE_INFINITY;
    for(int i = 0; i < probs.length; i++)
      max = Math.max(max, probs[i]);
    for(int i = 0; i < probs.length; i++)
      probs[i] = Math.exp(probs[i]-max);
    return normalize(probs);
  }

  public static boolean expNormalize(double[][] probs) {
    double max = Double.NEGATIVE_INFINITY;
    for(double[] v : probs)
      for(int i = 0; i < v.length; i++)
        max = Math.max(max, v[i]);
    for(double[] v : probs)
      for(int i = 0; i < v.length; i++)
        v[i] = Math.exp(v[i]-max);
    return normalize(probs);
  }

  public static boolean expNormalize(double[][][] probs) {
    double max = Double.NEGATIVE_INFINITY;
    for(double[][] m : probs)
      for(double[] v : m)
        for(int i = 0; i < v.length; i++)
          max = Math.max(max, v[i]);
    for(double[][] m : probs)
      for(double[] v : m)
        for(int i = 0; i < v.length; i++)
          v[i] = Math.exp(v[i]-max);
    return normalize(probs);
  }
  
  public static int[][] toInt(double[][] data) {
    int[][] newdata = new int[data.length][];
    for(int r = 0; r < data.length; r++) {
      newdata[r] = new int[data[r].length];
      for(int c = 0; c < data[r].length; c++)
        newdata[r][c] = (int)data[r][c];
    }
    return newdata;
  }

  public static double l1Dist(double[] x, double[] y) {
    double sum = 0;
    for(int i = 0; i < x.length; i++)
      sum += Math.abs(x[i]-y[i]);
    return sum;
  }

  public static double lInfDist(double[] x, double[] y) {
    double max = 0;
    for(int i = 0; i < x.length; i++)
      max = Math.max(max, Math.abs(x[i]-y[i]));
    return max;
  }

  public static double l2Dist(double[] x, double[] y) {
    return Math.sqrt(l2DistSquared(x, y));
  }
  public static double l2DistSquared(double[] x, double[] y) {
    double sum = 0;
    for(int i = 0; i < x.length; i++)
      sum += (x[i]-y[i])*(x[i]-y[i]);
    return sum;
  }

  public static double l2Norm(double[] x) { return Math.sqrt(l2NormSquared(x)); }
  public static double l2NormSquared(double[] x) {
    double sum = 0;
    for(int i = 0; i < x.length; i++)
      sum += x[i]*x[i];
    return sum;
  }

  public static double[] l2NormalizedMut(double[] x) {
    double norm = l2Norm(x);
    if(norm > 0) ListUtils.multMut(x, 1.0/norm);
    return x;
  }

  // If sum is 0, set to uniform
  // Return false if we had to set to uniform
  public static boolean normalizeForce(double[] data) {
    double sum = 0;
    for(double x : data) sum += x;
    if(sum == 0) {
      for(int i = 0; i < data.length; i++) data[i] = 1.0/data.length;
      return false;
    }
    else {
      for(int i = 0; i < data.length; i++) data[i] /= sum;
      return true;
    }
  }

  public static double[][] transpose(double[][] mat) {
    int m = mat.length, n = mat[0].length;
    double[][] newMat = new double[n][m];
    for(int r = 0; r < m; r++)
      for(int c = 0; c < n; c++)
        newMat[c][r] = mat[r][c];
    return newMat;
  }

  public static double[][] elementWiseMult(double[][] mat1, double[][] mat2) {
    int m = mat1.length, n = mat1[0].length;
    double[][] newMat = new double[m][n];
    for(int r = 0; r < m; r++)
      for(int c = 0; c < n; c++)
        newMat[r][c] = mat1[r][c] * mat2[r][c];
    return newMat;
  }

  public static void scalarMult(double[][] mat, double x) {
    int m = mat.length, n = mat[0].length;
    for(int r = 0; r < m; r++)
      for(int c = 0; c < n; c++)
        mat[r][c] *= x;
  }

  public static double[][] copy(double[][] mat) {
    int m = mat.length;
    double[][] newMat = new double[m][];
    for (int r = 0; r < m; r++) {
      int n = mat[r].length;
      newMat[r] = new double[n];
      for (int c = 0; c < n; c++)
        newMat[r][c] = mat[r][c];
    }
    return newMat;
  }

  public static boolean equals(double x, double y) {
    return Math.abs(x-y) < 1e-10;
  }
  public static boolean equals(double x, double y, double tol) {
    return Math.abs(x-y) < tol;
  }

  public static double round(double x, int numPlaces) {
    double scale = Math.pow(10, numPlaces);
    return Math.round(x * scale)/scale;
  }

  public static double[] round(double[] vec, int numPlaces) {
    double[] newVec = new double[vec.length];
    double scale = Math.pow(10, numPlaces);
    for(int i = 0; i < vec.length; i++)
      newVec[i] = Math.round(vec[i] * scale)/scale;
    return newVec;
  }

  public static double bound(double x, double lower, double upper) {
    if(x < lower) return lower;
    if(x > upper) return upper;
    return x;
  }
  public static int bound(int x, int lower, int upper) {
    if(x < lower) return lower;
    if(x > upper) return upper;
    return x;
  }
  // }

  public static double entropy(double[] probs) {
    double e = 0;
    for(double p : probs) {
      if(p > 0) e += -p * Math.log(p);
    }
    return e;
  }

  // Suppose we have a joint probability distribution p(X, Y)
  // and we want to calculate the conditional entropy
  // H(Y | X) = -\sum_x \sum_y p(x, y) \log p(y | x)
  // For a fixed x, we input y -> p(X=x, Y=y)
  // and get out the contribution to H(Y|X).
  // If we sum over all values of X, we get the desired result.
  public static double condEntropy(double[] probs) {
    double sum = ListUtils.sum(probs); // This is p(X) = \sum_y P(X, Y=y)
    double e = 0;
    for(double p : probs) {
      if(p > 0) e += -p * Math.log(p/sum);
    }
    assertIsFinite(sum);
    return e;
  }
  public static double condEntropy(double[][] probs) {
    double sum = 0;
    for(int i = 0; i < probs.length; i++)
      sum += condEntropy(probs[i]);
    return sum;
  }

  // Return log(gamma(xx))
  private static double[] logGammaCoeff =
  {
    76.18009172947146, -86.50532032941677,
    24.01409824083091, -1.231739572450155,
    0.1208650973866179e-2, -0.5395239384953e-5
  };
  public static double logGamma(double xx)
  {
    double x, y, tmp, ser;
    y = x = xx;
    tmp = x + 5.5;
    tmp -= (x+0.5) * Math.log(tmp);
    ser = 1.000000000190015;
    for(int j = 0; j <= 5; j++) ser += logGammaCoeff[j]/++y;
    return -tmp + Math.log(2.5066282746310005*ser/x);
  }

  /**
  * Stolen from Radford Neal's fbm package.
  * digamma(x) is defined as (d/dx) log Gamma(x).  It is computed here
  * using an asymptotic expansion when x>5.  For x<=5, the recurrence
  * relation digamma(x) = digamma(x+1) - 1/x is used repeatedly.  See
  * Venables & Ripley, Modern Applied Statistics with S-Plus, pp. 151-152.
  * COMPUTE THE DIGAMMA FUNCTION.  Returns -inf if the argument is an integer
  * less than or equal to zero.
  */
  public static double digamma(double x) {
    assert x > 0 : x;
    double r, f, t;
    r = 0;
    while (x<=5) {
      r -= 1/x;
      x += 1;
    }
    f = 1/(x*x);
    t = f*(-1/12.0 + f*(1/120.0 + f*(-1/252.0 + f*(1/240.0 + f*(-1/132.0
        + f*(691/32760.0 + f*(-1/12.0 + f*3617/8160.0)))))));
    return r + Math.log(x) - 0.5/x + t;
  }

  // Return log(exp(a)+exp(b))
  private static double logMaxValue = Math.log(Double.MAX_VALUE);
  public static double logAdd(double a, double b)
  {
    if(a > b) {
      if(Double.isInfinite(b) || a-b > logMaxValue) return a;
      return b+Math.log(1+Math.exp(a-b));
    }
    else {
      if(Double.isInfinite(a) || b-a > logMaxValue) return b;
      return a+Math.log(1+Math.exp(b-a));
    }
  }

  // Fast exponential
  // http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/
  public static double fastExp(double val) {
    final long tmp = (long) (1512775 * val + (1072693248 - 60801));
    return Double.longBitsToDouble(tmp << 32);
  }
  public static double fastLog(double val) {
    final double x = (Double.doubleToLongBits(val) >> 32);
    return (x - 1072632447) / 1512775;
  }
}
