package nmrsvd;

/**
 * <p>Title: </p>
 *
 * <p>Description: </p>
 *
 * <p>Copyright: Copyright (c) 2012</p>
 *
 * <p>Company: </p>
 *
 * @author Pascal P. Man
 * @version 2.0
 */

import java.util.Observable;
import java.io.*;
import com.sun.jna.Library;
import com.sun.jna.Native;
import com.sun.jna.NativeLibrary;

public class SVDMVCModelGPU extends Observable implements Library, Serializable {
	
  public static final String JNA_LIBRARY_NAME = "svdComplexDevice1";
  public static final NativeLibrary JNA_NATIVE_LIB = 
	NativeLibrary.getInstance(SVDMVCModelGPU.JNA_LIBRARY_NAME);
  static {
	Native.register(SVDMVCModelGPU.JNA_LIBRARY_NAME);
  }
  /**
	* Function declarations<br>
	* Original signature : <code>void selfSVD(const int, const int, 
	* float*, float*, float*, float*)</code><br>
	* <i>native declaration : line 5</i>
	*/
  public static native void selfSVD(int M, int N, 
	float[] matrixA, float[] matrixS, 
	float[] matrixR, float[] matrixT);
	
  /**
   * Function declarations<br>
   * Original signature : <code>void selfgemm(int, int, int, 
   * float*, float*, float*)</code><br>
   * <i>native declaration : line 8</i>
   */
  public static native void selfgemm(int M, int N, int K, 
	float[] matrixA, float[] matrixB, float[] matrixC);
	
  private static final long serialVersionUID = 1L;

  private double[] fidRe,
                   fidIm;

  private int maxNbComplex,
              td,    //initial TD/2 value
              tdeff, //TDeff/2 value
              n,     //number of lines of Hankel matrix
              m,     //number of columns of Hankel matrix
              minnm; //min(n,m)
  
  private double[] eigenvalue;  //1D singular value matrix
  
  private float[] MatrixA;   //1D Hankel matrix
  private float[] MatrixS;   //1D singular value matrix
  private float[] MatrixR;   //1D U matrix
  private float[] MatrixT;   //1D VT matrix
  
  public SVDMVCModelGPU(){
  }

  //---------------------//
  // Setters and getters //
  //---------------------//

  public int getmaxNbComplex() {return maxNbComplex;}
  public int setmaxNbComplex(int x) {return maxNbComplex = x;}
  
  public int getmvalue()      {return m;}
  public int setmvalue(int x) {return m = x;}

  public int getnvalue()      {return n;}
  public int setnvalue(int x) {return n = maxNbComplex - x + 1;}

  public int getTD()      {return td;}     //modified 25 march 2012
  public int setTD(int x) {return td = x;}

  public int getTDeff()      {return tdeff;}  //modified 25 march 2012
  public int setTDeff(int x) {return tdeff = x;}

  public double[] getfidRe()           {return fidRe;}
  public double[] setfidRe(double[] x) {return fidRe = x;}

  public double[] getfidIm()           {return fidIm;}
  public double[] setfidIm(double[] x) {return fidIm = x;}

  public double[] getSVD()           {return eigenvalue;} 
  public double[] setSVD(double[] x) {return eigenvalue = x;}

  public void essaiZsvd() {
	minnm = Math.min(n, m);
	eigenvalue = new double[minnm];

	MatrixA = new float[2*n*m];
	MatrixS = new float[minnm];
	MatrixR = new float[2*n*n];
	MatrixT = new float[2*m*m];	    
	    
	int p = 0;
	float ary[][] = new float[2*n][m];

    //first row
    for (int i = 0; i < m; i++) {
      ary[0][i] = (float)fidRe[p];
      ary[1][i] = (float)fidIm[p];
      p++;
    }
    
    //last column
    for (int i = 2; i < 2*n; i =i+2) {
      ary[i]  [m - 1] = (float)fidRe[p];
      ary[i+1][m - 1] = (float)fidIm[p];
      p++;
    }

    //remaining Hankel matrix elements
    for (int i = 2; i < 2*n; i=i+2) {
      for (int j = 0; j < m - 1; j++) {
        ary[i]  [j] = ary[i-2][j+1];
        ary[i+1][j] = ary[i-1][j+1];
      }
    }
    
    //input Hankel as vector for selfSVD
    int q = 0;
    for (int i = 0; i < m; i++) {  //column
      for (int j = 0; j < 2*n; j++) {  //line
        MatrixA[q] = ary[j][i];
        q++;
      }
    }
    
    selfSVD(n, m, MatrixA, MatrixS, MatrixR, MatrixT);

	for (int i = 0; i < minnm; i++){ 
	  eigenvalue[i] = MatrixS[i]; //conversion from float to double
	}
  }//end of essaiZsvd

  public void createFID() {
	
	float[] MatrixS2D =         new float[2*n*m];
	float[] MatrixTmp =         new float[2*n*m];	
	float[] MatrixHankel =      new float[2*n*m];
	double[][] hankelreDouble = new double[n][m];
	double[][] hankelimDouble = new double[n][m];
		
	//Initialization of MatrixS2D to zero
	for (int i = 0; i < 2*n*m; i++) {
	  MatrixS2D[i] = 0;
	}
    
	//conversion from 1D double to 2D float of singular value matrix
	int q = 0;
	for (int p = 0; p < minnm; p++) {
	  MatrixS2D[q] = (float)eigenvalue[p];
	  q = q + 2*(n+1);
	}

	selfgemm(n, m, m, MatrixS2D, MatrixT, MatrixTmp);		
	selfgemm(n, m, n, MatrixR, MatrixTmp, MatrixHankel);

    //Hankel matrix transformation from 1D vector to 2D matrix
	int p = 0;	
	for (int i = 0; i < m; i++) {  //column
	  for (int j = 0; j < n; j++) {   //row
		hankelreDouble[j][i] = MatrixHankel[p];	
		hankelimDouble[j][i] = MatrixHankel[p+1];
		p = p + 2;	
	  }
	}

    //clear
    for (int i = 0; i < maxNbComplex; i++) {
      fidRe[i] = 0;
      fidIm[i] = 0;
    }

    //first row, upper part
    for (int i = 0; i < m - 1; i++) {
      for (int j = 0; j <= i; j++) {
        fidRe[i] += hankelreDouble[i - j][j];
        fidIm[i] += hankelimDouble[i - j][j];
      }
      fidRe[i] /= i + 1;
      fidIm[i] /= i + 1;
    }

    //last column, middle part
    for (int i = 0; i <= n - m; i++) {
      for (int j = 0; j < m; j++) {
        fidRe[m - 1 + i] += hankelreDouble[m - 1 + i - j][j];
        fidIm[m - 1 + i] += hankelimDouble[m - 1 + i - j][j];
      }
      fidRe[m - 1 + i] /= m;
      fidIm[m - 1 + i] /= m;
    }

    //last column, lower part
    for (int i = 0; i < m - 1; i++) {
      for (int j = 0; j < m - 1 - i; j++) {
        fidRe[n + i] += hankelreDouble[n - m + 1 + i + j][m - 1 - j];
        fidIm[n + i] += hankelimDouble[n - m + 1 + i + j][m - 1 - j];
      }
      fidRe[n + i] /= m - 1 - i;
      fidIm[n + i] /= m - 1 - i;
    }
  }//end of createFID

}//end of class SVDMVCModelGPU