package brn.analysis;

import brn.sim.data.DiagramData;
import brn.sim.data.Line;
import brn.swans.radio.biterrormodels.BitErrorMask.SpecifiedErrors;

public class Statistics {

  public static Line autoCorrelate(DiagramData lInput, String title, int length, double limit) {
    if (null == lInput)
      return null;

    double[] correlation = new double[length];
    Statistics.autoCorrelate(lInput.getY(), correlation, limit);
    Line line = new Line(title);
    for (int x = 0; x < correlation.length; x++)
      line.add(x, correlation[x]);
    return line;
  }

  /**
   * Calculates the autocorrelation of the input array.
   *
   * @param v
   *          the input array
   * @param ret
   *          array to store the results
   * @param limit
   *          min autocorrelation
   */
  public static void autoCorrelate(double[] v, double[] ret, double limit) {

    // The devMean array will contain the deviation from the mean.
    // That is, v[i] - mean(v).
    double[] devMean = new double[v.length];

    double sum;
    int i;

    // Calculate the mean and copy the vector v, into devMean
    sum = 0;
    for (i = 0; i < v.length; i++)
      sum += v[i];
    double mean = sum / (double) v.length;

    // Compute the values for the devMean array. Also calculate the
    // denominator for the autocorrelation function.
    sum = 0;
    for (i = 0; i < v.length; i++) {
      devMean[i] = v[i] - mean;
      sum += devMean[i] * devMean[i];
    }
    double denom = sum / (double) v.length;

    // Calculate a std::vector of values which will make up
    // the autocorrelation function.
    double cor = 1.0;
    double numerator;
    int n;
    for (int shift = 1; shift <= ret.length && Math.abs(cor) > limit; shift++) {
      ret[shift - 1] = cor;
      n = v.length - shift;
      sum = 0.0;
      for (i = 0; i < n; i++)
        sum += devMean[i] * devMean[i + shift];
      numerator = sum / (double) n;
      cor = numerator / denom;
    }
  }

  /**
   * Calculates the autocorrelation of the input array.
   *
   * @param v1
   *          the input array 1
   * @param v2
   *          the input array 2
   * @param ret
   *          array to store the results
   * @param limit
   *          min autocorrelation
   */
  public static void correlate(double[] v1, double[] v2, double[] ret, double limit) {

    // The devMean array will contain the deviation from the mean.
    // That is, v[i] - mean(v).
    int N = Math.min(v1.length, v2.length);

    double[] devMean1 = new double[v1.length];
    double[] devMean2 = new double[v2.length];

    // Calculate the mean and copy the vector v, into devMean
    double sum1 = 0, sum2 = 0;
    for (int i = 0; i < v1.length; i++)
      sum1 += v1[i];
    double mean1 = sum1 / (double) v1.length;

    for (int i = 0; i < v2.length; i++)
      sum2 += v2[i];
    double mean2 = sum2 / (double) v2.length;

    // Compute the values for the devMean array. Also calculate the
    // denominator for the autocorrelation function.
    sum1 = 0;
    for (int i = 0; i < v1.length; i++) {
      devMean1[i] = v1[i] - mean1;
      sum1 += devMean1[i] * devMean1[i];
    }
    double dev1 = Math.sqrt(sum1 / (double) v1.length);

    sum2 = 0;
    for (int i = 0; i < v2.length; i++) {
      devMean2[i] = v2[i] - mean2;
      sum2 += devMean2[i] * devMean2[i];
    }
    double dev2 = Math.sqrt(sum2 / (double) v2.length);

    // Calculate a std::vector of values which will make up
    // the autocorrelation function.
    double cor = 1.0;
    double numerator;
    int n;
    for (int shift = 0; shift < ret.length && Math.abs(cor) > limit; shift++) {
      n = N - shift;
      sum1 = 0.0;
      for (int i = 0; i < n; i++)
        sum1 += devMean1[i] * devMean2[i + shift];
      numerator = sum1 / (double) n;
      cor = numerator / dev1 / dev2;
      ret[shift] = cor;
    }
  }

  // public double crossCorrelate(double[] v1, double[] v2) {
  // double mean1 = 0;
  // double mean2 = 0;
  // double variance = 0;
  // double std1 = 0;
  // double std2 = 0;
  //
  // for (int i = 0; i < v1.length; i++) {
  // mean1 += v1[i];
  // mean2 += v2[i];
  // }
  // double temp1, temp2;
  // mean1 /= (double)v1.length;
  // mean2 /= (double)v2.length;
  //
  // for (int i = 0; i < v1.length; i++) {
  // temp1 = v1[i];
  // temp2 = v2[i];
  // variance += (temp1-mean1)*(temp2-mean2);
  // std1 += (temp1-mean1)*(temp1-mean1);
  // std2 += (temp2-mean2)*(temp2-mean2);
  // }
  //
  // return variance/std1/std2;
  // }

  /**
   * Calculates Pearson correlation coefficient, r, from two double[]s
   */
  public static double correlationCoefficient(double[] x, double[] y) {
    double N = Math.min(x.length, y.length);
    double xTot = 0;
    double yTot = 0;
    double sqrXTot = 0;
    double sqrYTot = 0;
    double xYTot = 0;
    for (int i = 0; i < N; i++) {
      xTot += x[i];
      yTot += y[i];
      sqrXTot += Math.pow(x[i], 2);
      sqrYTot += Math.pow(y[i], 2);
      xYTot += (x[i] * y[i]);
    }
    double top = (N * xYTot) - (xTot * yTot);
    double botLeft = Math.sqrt((N * sqrXTot) - Math.pow(xTot, 2));
    double botRight = Math.sqrt((N * sqrYTot) - Math.pow(yTot, 2));
    return top / (botLeft * botRight);
  }


  public static double coVariance(double[] x, double[] y) {
    double N = Math.min(x.length, y.length);
    double xTot = 0;
    double yTot = 0;
    double sqrXTot = 0;
    double sqrYTot = 0;
    double xYTot = 0;
    for (int i = 0; i < N; i++) {
      xTot += x[i];
      yTot += y[i];
      sqrXTot += Math.pow(x[i], 2);
      sqrYTot += Math.pow(y[i], 2);
      xYTot += (x[i] * y[i]);
    }
    return (N* xYTot - xTot * yTot) / N / (N-1);
  }


  public static double[] timeCorrelation(double y1[], double y2[]) {
    int N = Math.min(y1.length, y2.length);
    int M = (int) Math.floor(1 + Math.log((N - 1.0) / 3.0) / Math.log(2.0));

    double[] ret = new double[M];
    double[] x1 = y1;
    double[] x2 = y2;

    int m = 1; // m = 2^j
    for (int j = 0; j < M; j++) {
      // calc current covariance
      ret[j] = correlationCoefficient(x1, x2);
      if (Double.isNaN(ret[j]))
        ret[j] = .0;

      // prepare next run, aggregate data
      int n = (int) Math.floor((N+1)/2);
      double[] w1 = new double[n];
      double[] w2 = new double[n];

      for (int i = 0; i < n-1; i++) {
        w1[i] = (x1[2*i] + x1[2*i + 1]) / 2;
        w2[i] = (x2[2*i] + x2[2*i + 1]) / 2;
      }
      w1[n-1] = (x1[2*n-2] + x1[N < 2*n ? 2*n-2 : 2*n-1]) / 2;
      w2[n-1] = (x2[2*n-2] + x2[N < 2*n ? 2*n-2 : 2*n-1]) / 2;

      N = n;
      x1 = w1;
      x2 = w2;
      m *= 2; // m = 2^j
    }

    return ret;
  }


  public static double[] timeCovariance(double y1[], double y2[]) {
    int N = Math.min(y1.length, y2.length);
    int M = (int) Math.floor(1 + Math.log((N - 1.0) / 3.0) / Math.log(2.0));

    double[] ret = new double[M];
    double[] x1 = y1;
    double[] x2 = y2;

    int m = 1; // m = 2^j
    for (int j = 0; j < M; j++) {
      // calc current covariance
      ret[j] = coVariance(x1, x2);

      // prepare next run, aggregate data
      int n = (int) Math.floor((N+1)/2);
      double[] w1 = new double[n];
      double[] w2 = new double[n];

      for (int i = 0; i < n-1; i++) {
        w1[i] = (x1[2*i] + x1[2*i + 1]) / 2;
        w2[i] = (x2[2*i] + x2[2*i + 1]) / 2;
      }
      w1[n-1] = (x1[2*n-2] + x1[N < 2*n ? 2*n-2 : 2*n-1]) / 2;
      w2[n-1] = (x2[2*n-2] + x2[N < 2*n ? 2*n-2 : 2*n-1]) / 2;

      N = n;
      x1 = w1;
      x2 = w2;
      m *= 2; // m = 2^j
    }

    return ret;
  }

  /**
   * computes the ALLAN DEVIATION
   *
   * @param y
   *          pointer to data, already allocated
   * @param tau0
   *          sampling periode in seconds
   *
   * @return ret[0]
   *          tau
   *         ret[1]
   *          overlapped Allan STD DEV (ADEV)
   *         ret[2]
   *          modified Allan STD DEV (MDEV)
   *         ret[3]
   *          Time Allan STD DEV (TDEV)
   *
   * Thanks to Alaa MAKDISSI.
   */
  public static double[][] allanDeviation(double y[], double tau0) {
    int k, m = 1, n = y.length; // m = 2^j ; j=0 ==> tau=tau0
    double s, v;

    // nb of points for the results
    int N = (int) Math.floor(1 + Math.log((n - 1.0) / 3.0) / Math.log(2.0));
    double[] D = new double[n];
    double[] tau = new double[N];
    double[] ADEV = new double[N];
    double[] MDEV = new double[N];
    double[] TDEV = new double[N];

    for (int j = 0; j < N; j++) {
      tau[j] = m * tau0;
      s = 0;
      v = 0; // s =Y(0) and v =Y(m)
      for (k = 0; k < m; k++) {
        s += y[k];
        v += y[k + m];
      }
      D[0] = (v - s) / m;

      // D(k)
      for (k = 1; k < n + 1 - 2 * m; k++)
        D[k] = D[k - 1] + (y[k - 1] - 2 * y[k + m - 1] + y[k + 2 * m - 1]) / m;

      // ADEV
      s = 0;
      for (k = 0; k < n + 1 - 2 * m; k++) {
        s = s + D[k] * D[k];
      }
      ADEV[j] = Math.sqrt(s / (n + 1 - 2 * m) / 2);

      // MDEV
      v = 0;
      for (int l = 0; l < n + 2 - (3 * m); l++) {
        s = 0;
        for (k = 0; k < m; k++) {
          s = s + D[l + k];
        }
        v = v + s * s;
      }
      v = v / (n + 2 - 3 * m);
      MDEV[j] = Math.sqrt(v / 2) / m;

      // TDEV
      TDEV[j] = MDEV[j] * m * tau0 / Math.sqrt(3.0);
      m = m * 2; // m = 2^j
    }

    return new double[][] { tau, ADEV, MDEV, TDEV };
  }

  public static double[] allanVariance(double y[]) {
    int k, m = 1, n = y.length; // m = 2^j ; j=0 ==> tau=tau0
    double s, v;

    // nb of points for the results
    int N = (int) Math.floor(1 + Math.log((n - 1.0) / 3.0) / Math.log(2.0));
    double[] D = new double[n];
    double[] AVAR = new double[N];

    for (int j = 0; j < N; j++) {
      s = 0;
      v = 0; // s =Y(0) and v =Y(m)
      for (k = 0; k < m; k++) {
        s += y[k];
        v += y[k + m];
      }
      D[0] = (v - s) / m;

      // D(k)
      for (k = 1; k < n + 1 - 2 * m; k++)
        D[k] = D[k - 1] + (y[k - 1] - 2 * y[k + m - 1] + y[k + 2 * m - 1]) / m;

      // AVAR
      s = 0;
      for (k = 0; k < n + 1 - 2 * m; k++) {
        s = s + D[k] * D[k];
      }
      AVAR[j] = s / (n + 1 - 2 * m) / 2;
      m = m * 2; // m = 2^j
    }

    return AVAR;
  }
  
  /**
   * Calculates allan variance from a bit error mask instead of an array.
   * This is incredibly slow because of all the map lookups. Don't use it
   * @param y
   * @param length
   * @return
   */
  public static double[] allanVariance(SpecifiedErrors y, int length) {
    int k, m = 1, n = length; // m = 2^j ; j=0 ==> tau=tau0
    double s, v;

    // nb of points for the results
    int N = (int) Math.floor(1 + Math.log((n - 1.0) / 3.0) / Math.log(2.0));
    double[] D = new double[n];
    double[] AVAR = new double[N];

    for (int j = 0; j < N; j++) {
      s = 0;
      v = 0; // s =Y(0) and v =Y(m)
      for (k = 0; k < m; k++) {
        s += y.errVal(k);
        v += y.errVal(k + m);
      }
      D[0] = (v - s) / m;

      // D(k)
      for (k = 1; k < n + 1 - 2 * m; k++)
        D[k] = D[k - 1] + (y.errVal(k - 1) - 2 * y.errVal(k + m - 1) + y.errVal(k + 2 * m - 1)) / m;

      // AVAR
      s = 0;
      for (k = 0; k < n + 1 - 2 * m; k++) {
        s = s + D[k] * D[k];
      }
      AVAR[j] = s / (n + 1 - 2 * m) / 2;
      m = m * 2; // m = 2^j
    }

    return AVAR;
  }

  
  
  public static double[] allanCovariance(double y1[], double y2[]) {
    int k, m = 1, n = Math.min(y1.length, y2.length); // m = 2^j ; j=0 ==> tau=tau0
    double s1, s2, v1, v2;

    // nb of points for the results
    int N = (int) Math.floor(1 + Math.log((n - 1.0) / 3.0) / Math.log(2.0));
    double[] D1 = new double[n];
    double[] D2 = new double[n];
    double[] AVAR = new double[N];

    for (int j = 0; j < N; j++) {
      s1 = s2 = 0;
      v1 = v2 = 0; // s =Y(0) and v =Y(m)
      for (k = 0; k < m; k++) {
        s1 += y1[k];
        s2 += y2[k];
        v1 += y1[k + m];
        v2 += y2[k + m];
      }
      D1[0] = (v1 - s1) / m;
      D2[0] = (v2 - s2) / m;

      // D(k)
      for (k = 1; k < n + 1 - 2 * m; k++) {
        D1[k] = D1[k - 1] + (y1[k - 1] - 2 * y1[k + m - 1] + y1[k + 2 * m - 1]) / m;
        D2[k] = D2[k - 1] + (y2[k - 1] - 2 * y2[k + m - 1] + y2[k + 2 * m - 1]) / m;
      }

      // AVAR
      s1 = 0;
      for (k = 0; k < n + 1 - 2 * m; k++) {
        s1 = s1 + D1[k] * D2[k];
      }
      AVAR[j] = s1 / (n + 1 - 2 * m) / 2;
      m = m * 2; // m = 2^j
    }

    return AVAR;
  }

  public static double[] average(double[] x1, int level) {
    if (null == x1 || 0 == x1.length)
      return null;

    int N = x1.length;
    int n = N / level;

    // aggregate data
    double[] w1 = new double[n];

    for (int i = 0; i < n-1; i++) {
      w1[i] = 0;
      for (int j = 0; j < level; j++)
        w1[i] += x1[i*level + j];
      w1[i] /= (double)level;
    }

    return w1;
  }

  /**
   * Returns the mean of an array of numbers.
   *
   * @param values  the values (<code>null</code> permitted, returns <code>Double.NaN</code>).
   *
   * @return The mean.
   */
  public static double mean(double[] values) {
      double result = Double.NaN;
      if (values != null && values.length > 0) {
          double sum = 0.0;
          int counter = 0;
          for (; counter < values.length; counter++) {
              sum = sum + values[counter];
          }
          result = (sum / counter);
      }
      return result;
  }

  /**
   * Returns the standard deviation of a set of numbers.
   *
   * @param data  the data.
   *
   * @return the standard deviation of a set of numbers.
   */
  public static double stddev(double[] data) {
      final double avg = mean(data);
      double sum = 0.0;

      for (int counter = 0; counter < data.length; counter++) {
          final double diff = data[counter] - avg;
          sum = sum + diff * diff;
      }
      return Math.sqrt(sum / (data.length - 1));
  }

  /**
   * Calculates autocorrelation from a bit error mask instead of an array 
   * this is incredibly slow because of all the map lookups. Don't use it
   * @param v
   *          the input bit error mask
   * @param length
   *          length of the mask
   * @param ret
   *          array to store the results
   * @param limit
   *          min autocorrelation
   */
  public static void autoCorrelate(SpecifiedErrors v, int length, double[] ret, double limit) {

    // The devMean array will contain the deviation from the mean.
    // That is, v[i] - mean(v).
    double[] devMean = new double[length];

    double sum;
    int i;

    // Calculate the mean and copy the vector v, into devMean
    sum = 0;
    
    sum += v.getNumErrors();
    double mean = sum / (double) length;

    // Compute the values for the devMean array. Also calculate the
    // denominator for the autocorrelation function.
    sum = 0;
    for (i = 0; i < length; i++) {
      devMean[i] = v.errVal(i) - mean;
      sum += devMean[i] * devMean[i];
    }
    double denom = sum / (double) length;

    // Calculate a std::vector of values which will make up
    // the autocorrelation function.
    double cor = 1.0;
    double numerator;
    int n;
    for (int shift = 1; shift <= ret.length && Math.abs(cor) > limit; shift++) {
      ret[shift - 1] = cor;
      n = length - shift;
      sum = 0.0;
      for (i = 0; i < n; i++)
        sum += devMean[i] * devMean[i + shift];
      numerator = sum / (double) n;
      cor = numerator / denom;
    }
    
  }

//  public static void main(String[] args) {
//    SecureRandom rand = new SecureRandom();
//
//    double[] v1 = new double[400000];
//    double[] v2 = new double[v1.length];
//    double[] v3 = new double[v1.length];
//
//    double amp = 1.;
//    double period = 512;
//
//    v1[0] = amp * Math.sin(0 * 2*Math.PI/period) + rand.nextDouble() - .5;
//    for (int i = 1; i < v1.length; i++) {
//      v1[i] = //amp * Math.signum(Math.sin(i * 2*Math.PI/period))
//        //amp * Math.sin(i * 2*Math.PI/period)
//        2*(rand.nextDouble() - .5)
//        ; //+ amp/4 * Math.sin(i * 2*Math.PI/(period*16));
//      v2[i - 1] = v1[i] + 2*(rand.nextDouble() - .5);
//    }
//    v2[v1.length - 1] = rand.nextDouble();
//    for (int i = 1; i < v3.length; i++)
//      v3[i] = 2*(rand.nextDouble() - .5);
//
//    double[] corr1 = new double[50];
//    autoCorrelate(v1, corr1, 0.0000001);
//
//    double[] corr2 = new double[50];
//    correlate(v1, v1, corr2, 0.0000001);
//    for (int i = 0; i < corr1.length; i++) {
//      if (Math.abs(corr1[i] - corr2[i]) > 0.0001)
//        System.out.println("error at " + i);
//    }
//
//    double x = correlationCoefficient(v1, v1);
//    double y = correlationCoefficient(v1, v2);
//
//    double[] cov1 = timeCorrelation(v1, v1);
//    double[] cov2 = timeCorrelation(v1, v2);
//    double[] cov3 = timeCorrelation(v1, v3);
//
//    double[] AVAR = allanVariance(v1);
//    double[] AAVAR = allanCovariance(v1, v1);
//    double[] ADEV = allanDeviation(v1, 1)[1];
//    double[] MDEV = allanDeviation(v1, 1)[2];
//    for (int i = 0; i < AVAR.length; i++) {
//      AVAR[i] = Math.log(Math.sqrt(AVAR[i])) / Math.log(2);
//      AAVAR[i] = Math.log(Math.sqrt(AAVAR[i])) / Math.log(2);
//      ADEV[i] = Math.log(ADEV[i]) / Math.log(2);
//      MDEV[i] = Math.log(MDEV[i]) / Math.log(2);
//      if (Math.abs(ADEV[i] - AVAR[i]) > 0.0001)
//        System.out.println("error at " + i);
//      if (Math.abs(ADEV[i] - AAVAR[i]) > 0.0001)
//        System.out.println("error at " + i);
//    }
//
//    System.out.println();
//  }

}
