// Copyright (c) 2016 - Patrick Schäfer (patrick.schaefer@hu-berlin.de)
// Distributed under the GLP 3.0 (See accompanying file LICENSE)
package dotnet54.classifiers.tschief.splitters.boss;

import dotnet54.tscore.data.TimeSeries;
import org.jtransforms.fft.DoubleFFT_1D;

import java.io.IOException;
import java.io.Serializable;

/**
 * The Momentary Fourier Transform is alternative algorithm of
 * the Discrete Fourier Transform for overlapping windows. It has
 * a constant computational complexity for in the window queryLength n as
 * opposed to O(n log n) for the Fast Fourier Transform algorithm.
 * <p>
 * It was first published in:
 * Albrecht, S., Cumming, I., Dudas, J.: The momentary fourier transformation
 * derived from recursive matrix transformations. In: Digital Signal Processing
 * Proceedings, 1997., IEEE (1997)
 *
 */
public class MFT implements Serializable {
  private static final long serialVersionUID = 8508604292241736378L;

  private int windowSize = 0;
  private int startOffset = 0;
  private double norm = 0;

  private transient DoubleFFT_1D fft = null;
  private boolean useMaxOrMin = false;

  public MFT() {
  }

  public MFT(int windowSize, boolean normMean, boolean lowerBounding) {
    this(windowSize, normMean, lowerBounding, false);
  }

  public MFT(int windowSize, boolean normMean, boolean lowerBounding, boolean useMinOrMax) {
    this.windowSize = windowSize;
    this.useMaxOrMin = useMinOrMax;

    initFFT();

    // ignore DC value?
    this.startOffset = normMean ? 2 : 0;
    this.norm = lowerBounding ? 1.0 / Math.sqrt(windowSize) : 1.0;
  }

  /**
   * Transforms a time series using the *discrete* fourier transform. Results in
   * a single Fourier transform of the time series.
   *
   * @param series the time series to be transformed
   * @param l          the number of Fourier values to keep
   * @return the first l Fourier values
   */
  public double[] transform(TimeSeries series, int l) {
    double[] data = new double[this.windowSize];
    //TODO using dimension 1 only
    System.arraycopy(series.data()[0], 0, data, 0, Math.min(this.windowSize, series.length()));
    this.fft.realForward(data);
    data[1] = 0; // DC-coefficient imaginary part

    // make it even length for uneven windowSize
    double[] copy = new double[l];
    int length = Math.min(this.windowSize - this.startOffset, l);
    System.arraycopy(data, this.startOffset, copy, 0, length);

    // norming
    int sign = 1;
    for (int i = 0; i < copy.length; i++) {
      copy[i] *= this.norm * sign;
      sign *= -1;
    }

    return copy;
  }

  /**
   * Transforms a time series, extracting windows and using *momentary* fourier
   * transform for each window. Results in one Fourier transform for each
   * window. Returns only the first l/2 Fourier coefficients for each window.
   *
   * @param series the time series to be transformed
   * @param l          the number of Fourier values to use (equal to l/2 Fourier
   *                   coefficients). If l is uneven, l+1 Fourier values are returned. If
   *                   windowSize is smaller than l, only the first windowSize Fourier
   *                   values are set.
   * @return returns only the first l/2 Fourier coefficients for each window.
   */
  public double[][] transformWindowing(TimeSeries series, int l) {
    int wordLength = useMaxOrMin ?
        Math.max(windowSize, l + this.startOffset) : // MUSE uses 'max'
        Math.min(windowSize, l + this.startOffset); // WEASEL uses 'min'
    wordLength += wordLength%2; // make it even
    double[] phis = new double[wordLength];

    for (int u = 0; u < phis.length; u += 2) {
      double uHalve = -u / 2;
      phis[u] = realPartEPhi(uHalve, this.windowSize);
      phis[u + 1] = complexPartEPhi(uHalve, this.windowSize);
    }

    // means and stddev for each sliding window
    int end = Math.max(1, series.length() - this.windowSize + 1);
    double[] means = new double[end];
    double[] stds = new double[end];
    //TODO using dimension 1 only
    series.calcIncrementalMeanStddev(this.windowSize, series.data()[0], means, stds);

    double[][] transformed = new double[end][];

    // holds the DFT of each sliding window
    double[] mftData = new double[wordLength];
    double[] data = series.data(0);   //TODO add support for multivariate

    for (int t = 0; t < end; t++) {
      // use the MFT
      if (t > 0) {
        for (int k = 0; k < wordLength; k += 2) {
          double real1 = (mftData[k] + data[t + this.windowSize - 1] - data[t - 1]);
          double imag1 = (mftData[k + 1]);

          double real = complexMultiplyRealPart(real1, imag1, phis[k], phis[k + 1]);
          double imag = complexMultiplyImagPart(real1, imag1, phis[k], phis[k + 1]);

          mftData[k] = real;
          mftData[k + 1] = imag;
        }
      }
      // use the DFT for the first offset
      else {
        double[] dft = new double[this.windowSize];
        //TODO add support for multivariate
        System.arraycopy(series.data(0), 0, dft, 0, Math.min(this.windowSize, series.length()));

        this.fft.realForward(dft);
        dft[1] = 0; // DC-coefficient imag part

        // if windowSize > mftData.queryLength, the remaining  should be 0 now.
        System.arraycopy(dft, 0, mftData, 0, Math.min(mftData.length, dft.length));
      }

      // normalization for lower bounding
      double[] copy = new double[l];
      System.arraycopy(mftData, this.startOffset, copy, 0, Math.min(l, mftData.length-this.startOffset));

      transformed[t] = normalizeFT(copy, stds[t]);
    }

    return transformed;
  }

  /**
   * Calculate the real part of a multiplication of two complex numbers
   */
  private static double complexMultiplyRealPart(double r1, double im1, double r2, double im2) {
    return r1 * r2 - im1 * im2;
  }

  /**
   * Caluculate the imaginary part of a multiplication of two complex numbers
   */
  private static double complexMultiplyImagPart(double r1, double im1, double r2, double im2) {
    return r1 * im2 + r2 * im1;
  }

  /**
   * Real part of e^(2*pi*u/M)
   */
  private static double realPartEPhi(double u, double M) {
    return Math.cos(2 * Math.PI * u / M);
  }

  /**
   * Imaginary part of e^(2*pi*u/M)
   */
  private static double complexPartEPhi(double u, double M) {
    return -Math.sin(2 * Math.PI * u / M);
  }

  /**
   * Apply normalization to the Fourier coefficients to allow lower bounding in Euclidean space
   */
  private double[] normalizeFT(double[] copy, double std) {
    double normalisingFactor = (TimeSeries.zNormalize && std > 0 ? 1.0 / std : 1.0) * this.norm;
    int sign = 1;
    for (int i = 0; i < copy.length; i++) {
      copy[i] *= sign * normalisingFactor;
      sign *= -1;
    }
    return copy;
  }

  public int getStartOffset() {
    return startOffset;
  }

  private void readObject(java.io.ObjectInputStream in) throws IOException, ClassNotFoundException {
    in.defaultReadObject();
    initFFT();
  }

  private void initFFT() {
    this.fft = new DoubleFFT_1D(this.windowSize);
  }
}
