import SVDC_j.*;
import Blas_j.*;
import Blas_j2.*;
import java.lang.*;
import java.util.*;
//import gt2.*;

public class gt_util9 extends Object
{
  static double pdv, pdelta;
  static int pspandim, pds;
 
  public gt_util9() 
  {
    double pdv = 0.0;
    double pdelta = 0.0;
    int pspandim = 0;
    int pds = 0;
  }

  public static int gen_norm_variates(int n, int p, double vars[]) 
  {
    int check =1;
    double frunif[];
    frunif = new double[2];
    double r, fac, frnorm[];
    frnorm = new double[2];
    int np; 
    np=n*p+1;
    np=np/2;
    Random rand = new Random();
//  System.out.println("n = "+n+" p = "+p);

    for (int i=0; i<np; i++) 
    {
      while (check>0) 
      {
        frunif[0] = 2.0*rand.nextDouble()-1.0;
        frunif[1] = 2.0*rand.nextDouble()-1.0;
        r = frunif[0]*frunif[0] + frunif[1]*frunif[1];
        if (r<1) 
        {
          check = 0;
          fac = Math.sqrt(-2.0*Math.log(r)/r);
          frnorm[0] = frunif[0]*fac;
          frnorm[1] = frunif[1]*fac;
        }
      } 
      check = 1;
      vars[2*i] = (double) (frnorm[0]);
// Following one line was added on 10/30/98.
      if ((2*i+1) <= (n*p-1)) 
      vars[2*i+1] = (double) (frnorm[1]);
    }
    return 1;
  }

// Generate two random p dimensional vectors to form new ending basis

  public static void new_random_basis(double basis[][], int ncols,
       int ntour_vars, int tour_vars[], int ndim)
  {

// INPUT : ncols, ntour_vars, tour_vars[], ndim
// OUTPUT : basis
    Blas_j2 Bj2;
    Bj2 = new Blas_j2();
    double tmp_basis[];
    tmp_basis = new double[ndim*ncols];

    gen_norm_variates(ndim,ntour_vars,tmp_basis);

// Zero Basis before Filling

    for (int k=0; k<ndim; k++) 
    {
      for (int j=0; j<ncols; j++) 
      {
        basis[k][j] = 0.0;
      }
    }

    for (int k=0; k<ndim; k++) 
    {
      for (int j=0; j<ntour_vars; j++) 
      {
        basis[k][tour_vars[j]] = tmp_basis[k*ntour_vars+j];
      }
    }

    for (int k=0; k<ndim; k++) 
    {
      Bj2.norm(ncols,basis[k]);
    }

    if (ndim > 1)
    {
      for (int k=0; k<(ndim-1); k++) 
      {
        for (int j=k+1; j<ndim; j++) 
        { 
          Bj2.gram_schmidt(basis[k], basis[j], ncols);
        }
      }
    }
  }
                                                                           
  static void path(
             double u0[][], double u1[][],  // U0 : Fa, U1 : Fz
             double v0[][], double v1[][], double b[][],
             double w0[][], double w1[][], int ncols, int tour_vars[], 
             int ndim, double tau[], double tv[][], 
             double step) throws SVDCException 
  {

// INPUT : u0, u1, ndim, ncols, step
// OUTPUT : tau, tv, v0, v1, w0, w1, b
// int ncols;
    double tmpf;
    double tol = 0.000001;   // original tol=0.000001
// PROJ_DIM =2
//    double costf[PROJ_DIM], sintf[PROJ_DIM];
    double costf[], sintf[];
    costf = new double[ndim];
    sintf = new double[ndim];
    double dv, delta;
    double sorttau[];
    sorttau = new double[3];
//    sorttau = new double[ndim];
    int tauindx[];
    tauindx = new int[3];
//    tauindx = new int[ndim];
    int spandim, ds;
    spandim = ndim;
// Variables for SVDC
    double uu[][];
    double vv[][];
    double work[];
    uu = new double[ndim][ndim];
    vv = new double[ndim][ndim];    
    work = new double[ndim];

// Form U0'U1
    for (int i = 0; i < ndim; i++) 
    {
      for (int j=0; j<ndim; j++) 
      {
        tmpf=0.0;
        for (int k = 0; k < ncols; k++) 
        {
          tmpf += (u0[i][k]*u1[j][k]);
        }
        tv[i][j] = tmpf;
      }
    }         

// Singular Value Decomposition
// tau[] <- s[] in SVDC_j.dsvdc_j    

    SVDC_j svd1 = new SVDC_j();
    double e[];
    e = new double[ndim];

/* SVD
   INPUT :  tv : Matrix to be SVDecomposed
   OUTPUT:  tau : Singular Values of tv
            uu : left singular vectors
            vv : right singular vectors*/ 


    svd1.dsvdc_j(tv, ndim, ndim, tau, e, uu, vv, work, 11);


// Sort Singular Values and TV, V1
    int k;
    if (ndim > 1)
    {
      for (int i=0; i<ndim-1; i++) 
      {
        tmpf = tau[i];
        k = i;
        for (int j=i+1; j<ndim; j++) 
        {
          if (tau[j] >= tmpf) 
          {
            tmpf = tau[j]; k = j;
          }
        } 
        if (k != i) 
        {
          tau[k] = tau[i];
          tau[i] = tmpf;
          for (int j=0; j<ndim; j++) 
          {
            tmpf = uu[j][i]; 
            uu[j][i] = uu[j][k]; 
            uu[j][k] = tmpf;
          }
          for (int j=0; j<ndim; j++) 
          {
            tmpf = vv[j][i]; 
            vv[j][i] = vv[j][k]; 
            vv[j][k] = tmpf;
          }    
	}
      }        
    }            

// Check if any singular values too close to 1
//   System.out.println("singular value = "+tau[0]+" "+tau[1]);                                                                          
    for (int kk=0; kk<ndim; kk++) 
    {
      tmpf = tau[kk]+tol;
/////////////// May 2, 1998

      if (tmpf > 1.0) 
      {
        spandim = spandim-1; 
//      System.out.println("kk is "+kk+" Spandim is "+spandim+" tau is "+tau[kk]);
      }

/////////////// I put "//" on lineno 166. 
    }
                                                                             
// Generate Principal Directions

    Blas_j2.mat_mult(u0,uu,v0,ncols,ndim,ndim);
    for (int ku=0; ku<ndim; ku++) 
    {
      for(int i=0; i<ndim; i++) 
      {
        uu[i][ku] = vv[i][ku];
        vv[i][ku] = 0.0;
      }
    }

    Blas_j2.mat_mult(u1,uu,v1,ncols,ndim,ndim);
    // for (int i=0; i<ndim; i++) {
   //  System.out.println("v0: " +v0[i][0]+" "+v0[i][1]+" "+v0[i][2]);
   // }
// for (int i=0; i<ndim; i++) {
//  System.out.println("v1: " +v1[i][0]+" "+v1[i][1]+" "+v1[i][2]);
// }
//    for (int kk=ndim-1; kk>ndim-spandim-1; kk--) 
    for (int kk=ndim-1; kk>ndim-spandim-1; kk--) 
    {
      Blas_j2.gram_schmidt(v0[kk], v1[kk], ncols);
    }
// System.out.println("         ");
//    for (int i=0; i<ndim; i++) 
//    {
//  System.out.println("v1: " +v1[i][0]+" "+v1[i][1]+" "+v1[i][ndim]);
//    }

    for (int i=0; i<ndim; i++) 
    {
      if (tau[i] > 1.0) 
        tau[i]=1.0;
      if (tau[i] < -1.0) 
        tau[i]=-1.0;
    }  
// Compute Principal Angles

    for (int kk=0; kk<spandim; kk++) 
    {
      uu[0][kk] = Math.acos(tau[ndim-kk-1]);
    }
    for (int kk=spandim; kk<ndim; kk++) 
    {
      uu[0][kk] = 0.0;
    }
    for (int kk=0; kk<ndim; kk++) 
    {
      tau[kk] = uu[0][kk];
    }
                                                                            
// Construct Preprojection Basis
    for (int kk=0; kk<spandim; kk++) 
    {
      for (int j=0; j<ncols; j++) 
      {
        b[2*kk][j] = v0[ndim-kk-1][j];
        b[2*kk+1][j] = v1[ndim-kk-1][j];
//       System.out.println("b["+Integer.toString(2*kk)+"]["
//                          +Integer.toString(j)+"]= "+b[2*kk][j]);
//       System.out.println("b["+Integer.toString(2*kk+1)+"]["
//                          +Integer.toString(j)+"]= "+b[2*kk+1][j]);
      }
    }

//    for (int kk=2; kk<ndim; kk++) 
    for (int kk=spandim; kk<ndim; kk++) 
    {
      for (int j=0; j<ncols; j++) 
      {
        b[spandim+kk][j] = v0[ndim-kk-1][j];
      }
    }

    ds= spandim+ndim;
                                                                     
// Construct Frames
                                                                           
    for (int i=0; i<ndim; i++) 
    {
      for (int kkk=0; kkk<ds; kkk++) 
      {
        tmpf = 0.0;
        for (int j=0; j<ncols; j++) 
        {
          tmpf += (b[kkk][j]*u0[i][j]);
        }
        w0[kkk][i] = tmpf;
        tmpf=0.0;
        for (int j=0; j<ncols; j++) 
        {
          tmpf += (b[kkk][j]*u1[i][j]);
        }
        w1[kkk][i] = tmpf; 
//       System.out.println("wo="+w0[kkk][i]); 
//       System.out.println("w1="+w1[kkk][i]);       
      }
    }   

    for (int i=0; i<ndim; i++) 
    {
      costf[i] = Math.cos(tau[i]);
      sintf[i] = Math.sin(tau[i]);
    }

    tmpf = 0.0;

    for (int kj=0; kj<spandim; kj++) 
    {
      tmpf += (tau[kj]*tau[kj]);
    }
    pdv = Math.sqrt(tmpf);
    pdelta = step/pdv;
    pspandim = spandim;
    pds = ds;                     
  }


// Checks if the current position is within an increment of
// the Target Plane
                                                                           
  public static int reached_target_plane(int ntour_vars, double tinc[],
                       double tau[], double delta, int spandim, 
                       int got_new_target_plane) 
  {

    int ret_val=0;
    int all_equal=1;
    double tol=0.000001;  //original tol=0.000001
    for (int k=0; k<spandim; k++) 
    {
      if (Math.abs(tinc[k]-tau[k]) > tol) 
      { 
        all_equal=0;
      }
    }

    if (all_equal >0) 
      ret_val=1;
    else 
    {
      for (int k=0; k<spandim; k++) 
      {
        tinc[k] += (delta*tau[k]);
//   System.out.println("tinc="+tinc[k]+" tau="+tau[k]);
      }
      for (int k=0; k<spandim; k++) 
      {
        if (tinc[k] > tau[k]) 
          ret_val=1;
      }
    } 

    return ret_val;

  }

/* MOVE these two methods to gt2.java

  public void increment_tour() 
  {
    gt2 GT1;
    GT1 = new gt2();
    GT1.tour_reproject(); 
  }

  public static void finishing_step(double tinc[], double tau[], int spandim) 
  {
    gt2 GT1;
    GT1 = new gt2();
    int do_final_step=0;

    for (int k=0; k<spandim; k++) 
    {
      if (tinc[k] != tau[k]) 
      {
        tinc[k]=tau[k]; 
        do_final_step=1;
      }
      if (do_final_step>0) 
      {
        GT1.tour_reproject(); 
      }
    }
  }
*/

// Find U0 basis

  public static void u0_gets_current_pos(double u0[][], double u[][],
          int ncols, int ntour_vars, int ndim) 
  {
// INPUT : u, ncols, ntour_vars, ndim
// OUTPUT : u0

    Blas_j2 Bj2;
    Bj2 = new Blas_j2();

    for (int k=0; k<ndim; k++) 
    {
      for (int j=0; j<ncols; j++) 
      {
        u0[k][j] = u[k][j];
      }
    }

    for (int k=0; k<ndim; k++) 
    {
      Bj2.norm(ncols,u0[k]); 
    }

    if (ndim > 1)
    {
      for (int k=0; k<(ndim-1); k++) 
      {
        for (int j=k+1; j<ndim; j++) 
        {
          Bj2.gram_schmidt(u0[k], u0[j], ncols);
        }
      }
    }
  }

}