#include <iostream.h>
#include <stdio.h>
#include <memory.h>
#include <stdlib.h>
#include <Console.h>

#include <VectorDouble.h>
#include <Sof.h>
#include "Kpca.h"

// 0: data set 1
// 1: data set 2
//
long EVAL_ID = 0;

// 11: data set 1
// 5: data set 2
//
long CLASS_NUM = 0;

#define TRAIN_NUM      999
long PATTERN_LEN = 0;

#define RBF_VAR0        0.015
#define RBF_VAR1        0.25

static long MAX_EIGEN = 0;
#define EIGEN_RATIO    0.99999

#define MARGIN_H         1
#define MARGIN_L         -1

long KERNEL_TYPE = 0;
double POLY_DEGREE = 5;

unichar KFILE[] = L"k.data";
char OUTFILE[] = "class.data";

Sof k_file;
String k_filename;

// merge vectors
//
boolean merge(Vector<VectorDouble>& patterns_a) {

  long len = patterns_a.length();
  if (len != ((long)(len / 5)) * 5) {
    printf("error in pattern length!\n");
    return false;
  }

  long merge_len = len / 5;
  VectorDouble vec_tmp;
  vec_tmp.setLength(patterns_a(0).length() * 5);
  for (long i = 0; i < merge_len; i++) {
    long start_pos = i * 5;
    if (start_pos != i) {
      patterns_a(i).assign(patterns_a(start_pos));
    }
    patterns_a(i).concat(patterns_a(start_pos + 1));
    patterns_a(i).concat(patterns_a(start_pos + 2));
    patterns_a(i).concat(patterns_a(start_pos + 3));
    patterns_a(i).concat(patterns_a(start_pos + 4));
  }

  return patterns_a.setLength(merge_len, true);
}
  
    
		   
// kernel method:
//
boolean kernel(MatrixDouble& output_a,
	       const Vector<VectorDouble>& patterns_a,
	       const Vector<VectorDouble>& test_patterns_a, boolean self_a) {

  // check the arguments
  //
  if (test_patterns_a.length() <= 0) {
    printf("empty test patterns!\n");
    return false;
  }

  double RBF_VAR;
  if (EVAL_ID == 0) {
    RBF_VAR = RBF_VAR0;
  }
  else {
    RBF_VAR = RBF_VAR1;
  }
    
  if (KERNEL_TYPE == 0) {

    if (self_a) {

      // loop through the input test feature vectors and compute the K matrix
      //
      long len = patterns_a.length();
      output_a.setDimensions(len, len, false, Integral::SYMMETRIC);
      VectorDouble vec_tmp;
      Double value;
    
      for (long i = 0; i < len; i++) {
	for (long j = 0; j < i; j++) {
	  vec_tmp.sub(patterns_a(i), patterns_a(j));
	  value = vec_tmp.sumSquare();
	  value /= -RBF_VAR;
	  value.exp();
	  output_a.setValue(i, j, value);
	}
	output_a.setValue(i, i, 1);
      }
    }
    
    else {
      
      // loop through the input test feature vectors and compute the
      // K_test matrix
      //
      long len = patterns_a.length();
      long test_num = test_patterns_a.length();
      output_a.setDimensions(test_num, len, false, Integral::FULL);
      VectorDouble vec_tmp;
      Double value;

      for (long i = 0; i < test_num; i++) {
	for (long j = 0; j < len; j++) {
	  vec_tmp.sub(test_patterns_a(i), patterns_a(j));
	  value = vec_tmp.sumSquare();
	  value /= -RBF_VAR;
	  value.exp();
	  output_a.setValue(i, j, value);
	}
      }
    }
  }

  else if (KERNEL_TYPE == 1) {

    // polynomial kernel with degree
    //
    if (self_a) {
      
      // loop through the input test feature vectors and compute the K matrix
      //
      long len = patterns_a.length();
      output_a.setDimensions(len, len, false, Integral::SYMMETRIC);
      VectorDouble vec_tmp;
      Double value;
    
      for (long i = 0; i < len; i++) {
	for (long j = 0; j <= i; j++) {
	  value = patterns_a(i).dotProduct(patterns_a(j)) + 1;
	  value.pow((double)POLY_DEGREE);
	  output_a.setValue(i, j, value);
	}
      }
    }
    
    else {
      
      // loop through the input test feature vectors and compute the
      // K_test matrix
      //
      long len = patterns_a.length();
      long test_num = test_patterns_a.length();
      output_a.setDimensions(test_num, len, false, Integral::FULL);
      VectorDouble vec_tmp;
      Double value;

      for (long i = 0; i < test_num; i++) {
	for (long j = 0; j < len; j++) {
	  value = test_patterns_a(i).dotProduct(patterns_a(j)) + 1;
	  value.pow((double)POLY_DEGREE);
	  output_a.setValue(i, j, value);
	}
      }
    }
  }
  
  // exit gracefully
  //
  return true;
}


// method: eval
//
//  MatrixDouble& output: (output) output data
//  const Vector<VectorDouble>& input: (input) feature data
//
// return: a boolean value indicating status
//
// this method performs kernel conversion on input feature data and output
// pricinple components of the extended space
//
boolean eval(const char* name_a) {

  // read test patterns
  //
  Vector<VectorDouble> test_patterns;
  FILE* fp_ref = fopen(name_a, "r");

  char* buf_ref;
  buf_ref = new char[999];
  if (!buf_ref) {
    printf("error in memory!");
    return false;
  }

  memset(buf_ref, 999, 0);
  char buf_num[8];

  char* pref;
  char* pnum;

  int num_data = 0;
  long num_colon;

  while (fgets(buf_ref, 999, fp_ref) != (char*)NULL) {

    // get the chars before ':'
    //
    pref = buf_ref;
    pnum = buf_num;
    num_colon = 7;
    while ((num_colon > 0) && (*pref != ':')) {
      num_colon--;
      *pnum++ = *pref++;
      
    }
    *pnum = 0;
    pref++;

    if (num_colon <= 0) {

      // no ':'
      //
      pref = buf_ref;
    }
    
    // try to find out the PATTERN_LEN
    //
    if (PATTERN_LEN <= 0) {
      char* pCount = pref;
      PATTERN_LEN = 0;
      while(*pCount) {
	if (*pCount++ == '.') {
	  PATTERN_LEN++;
	}
      }

      // output the pattern length
      //
      printf("The pattern length = %ld\n", PATTERN_LEN);

      if (PATTERN_LEN <=0) {
	delete buf_ref;
	printf("data file %s is wrong\n", name_a);
	return false;
      }

      else if (PATTERN_LEN == 10) {
	EVAL_ID = 0;
      }

      else if (PATTERN_LEN == 39) {
	EVAL_ID = 1;
      }
    }

    // resize the test patterns buffer
    //
    long len = test_patterns.length();
    test_patterns.setLength(len + 1, true);
    test_patterns(len).setLength(PATTERN_LEN);

    double temp;
    for (long i = 0; i < PATTERN_LEN; i++) {
      sscanf(pref, " %lf", &temp);
      test_patterns(len)(i) = temp;
      while(*pref == ' ') {
	pref++;
      }
      while(*pref != ' ') {
	pref++;
      }
    }

    num_data++;
    memset(buf_ref, 999, 0);
  }

  if (EVAL_ID == 0) {
    CLASS_NUM = 11;
  }
  else {
    CLASS_NUM = 5;
  }
  
  if (EVAL_ID == 1) {
    num_data /= 5;
    
    if (!merge(test_patterns)) {
      printf("error in merge!");
    }
  }

  delete buf_ref;
  buf_ref = 0;
  printf("Total %d testing tokens\n", num_data);
  fclose(fp_ref);

  // read parameters
  //
  k_filename.assign(KFILE);
  k_file.open(k_filename, File::READ_ONLY, File::BINARY);

  Long len;
  Vector<VectorDouble> patterns;
  MatrixDouble K;
  MatrixDouble eigvecs_trim;
  MatrixDouble a;
  Double norm_val;
  
  len.read(k_file, (long)0);
  patterns.read(k_file, (long)1);
  K.read(k_file, (long)2);
  eigvecs_trim.read(k_file, (long)3);
  a.read(k_file, (long)4);
  norm_val.read(k_file, (long)5);

  k_file.close();
  
  printf("extracting testing features...\n");

  // normalize test patterns
  //
  for (long i = 0; i < test_patterns.length(); i++) {
    test_patterns(i).mult(norm_val);
  }
  
  // check whether we need extract test features
  //
  if (test_patterns.length() <= 0) {
    printf("empty test patterns!\n");
    return false;
  }

  // loop through the input test feature vectors and compute the K_test matrix
  //
  long test_num = test_patterns.length();
  MatrixDouble K_test(test_num, len, Integral::FULL);
  if (!kernel(K_test, patterns, test_patterns)) {
    return false;
  }
  
  // centering matrix K_test in the feature space
  //
  MatrixDouble unit(len, len, Integral::FULL);
  unit.assign((double)1 / len);
  MatrixDouble mat_out_test(test_num, len, Integral::FULL);
  MatrixDouble mat_tmp_test(test_num, len, Integral::FULL);
  MatrixDouble unit_test(test_num, len, Integral::FULL);
  unit_test.assign((double)1 / len);

  // K_test * unit
  //
  mat_tmp_test.mult(K_test, unit);

  // K_test_n = K_test - K_test * unit
  //
  mat_out_test.sub(K_test, mat_tmp_test);

  // unit_test * K
  //
  mat_tmp_test.mult(unit_test, K);

  // K_test_n = K_test - unit_test * K - K_test * unit
  //
  mat_out_test.sub(mat_tmp_test);

  // unit_test * K * unit
  //
  mat_tmp_test.mult(unit);

  // K_test_n = K_test - unit_test * K - K_test * unit + unit_test * K * unit
  //
  mat_out_test.add(mat_tmp_test);

  // output K_test_n
  //
  MatrixDouble K_test_n(mat_out_test);

  // extract new features for test set
  //
  MatrixDouble test_features;
  test_features.mult(K_test_n, eigvecs_trim);

  // compute the scores for test features
  //
  MatrixDouble Yt;
  VectorDouble vec_tmp0;
  long feature_len = eigvecs_trim.getNumColumns();
  long test_feature_num = test_features.getNumRows();
  Yt.setDimensions(test_feature_num, feature_len + 1, false);
  vec_tmp0.setLength(test_feature_num);
  vec_tmp0.assign((double)1);
  Yt.setColumn(0, vec_tmp0);

  // copy test features to Yt buffer
  //
  for (long i = 0; i < feature_len; i++) {
    test_features.getColumn(vec_tmp0, i);
    Yt.setColumn(i + 1, vec_tmp0);
  }

  // compute scores
  //
  Yt.mult(a);

  // count the highest score to determine the class
  //
  VectorLong vec_class(test_feature_num);
  long max_index;
  double max_val;
  double cur_val;
  for (long i = 0; i < test_feature_num; i++) {
    max_index = 0;
    max_val = Yt(i, 0);
    for (long j = 1; j < CLASS_NUM; j++) {
      cur_val = Yt(i, j);
      if (max_val < cur_val) {
	max_val = cur_val;
	max_index = j;
      }
    }
    vec_class(i) = max_index + 1;
  }

  printf("successful classification...\n");

  // write the result to data file
  //
  fp_ref = fopen(OUTFILE, "w");

  if (EVAL_ID == 0) {
    for (long i = 0; i < test_feature_num; i++) {
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
    }
  }
  else {
    for (long i = 0; i < test_feature_num; i++) {
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
    }
  }
  fclose(fp_ref);

  return true;
}
  
// method: compute
//
//  MatrixDouble& output: (output) output data
//  const Vector<VectorDouble>& input: (input) feature data
//
// return: a boolean value indicating status
//
// this method performs kernel conversion on input feature data and output
// pricinple components of the extended space
//
boolean compute(MatrixDouble& features_a, MatrixDouble& test_features_a,
		const Vector<VectorDouble>& patterns_a,
		const Vector<VectorDouble>& test_patterns_a) {

  // declare local variable
  //
  long len = patterns_a.length();
  VectorDouble vec_tmp;
  Double value;

  printf("extracting training features...\n");

  // set space for output matrix
  //
  MatrixDouble K(len, len, Integral::SYMMETRIC);

  // loop through the input feature vectors and compute the K matrix
  //
  if (!kernel(K, patterns_a, patterns_a, true)) {
    Integral::exit();
  }
  
  // centering matrix K in the feature space
  //
  MatrixDouble mat_out(len, len, Integral::FULL);
  MatrixDouble mat_tmp(len, len, Integral::FULL);
  MatrixDouble unit(len, len, Integral::FULL);
  unit.assign((double)1 / len);

  // K * unit
  //
  mat_tmp.mult(K, unit);

  // K_n = K - K * unit
  //
  mat_out.sub(K, mat_tmp);

  // unit * K
  //
  mat_tmp.mult(unit, K);

  // K_n = K - unit * K - K * unit
  //
  mat_out.sub(mat_tmp);

  // unit * K * unit
  //
  mat_tmp.mult(unit);

  // K_n = K - unit * K - K * unit + unit * K * unit
  //
  mat_out.add(mat_tmp);

  // output mat_out
  //
  MatrixDouble K_n(mat_out);

  // compute the eigen values and eigenvectors
  //
  VectorDouble eigvals;
  MatrixDouble eigvecs;
  mat_out.mult((double) 1 / len);

  // try to convert the matrix to symmetric matrix
  //
  if (mat_out.changeType(Integral::SYMMETRIC)) {
    printf("Symmetric K matrix\n");
  }
  
  //  if (mat_out.eigen(eigvals, eigvecs)) {
  if (!computeEigen(eigvals, eigvecs, mat_out)) {
    return false;
  }

  // choose the number of eigenvalus which are worth of keeping
  //
  MAX_EIGEN = chooseEigenNum(eigvals);
    
  // normalize eigenvectors
  //
  for (long i = 0; i < MAX_EIGEN; i++) {
      
    // compute the normalize factor
    //
    value = (double) 1 / eigvals(i);
    value.sqrt();

    // normalize the eigenvector
    //
    eigvecs.getRow(vec_tmp, i);
    vec_tmp.mult(value);
    eigvecs.setRow(i, vec_tmp);
  }
      
  // outupt eigenvalues/vectors
  //
  printf("Choosed number for eigen values: %ld\n", MAX_EIGEN);
  //eigvals.debug(L"eigenvalues");

  // trim the eigvecs and swap it to be column-based eigen vectors
  //
  MatrixDouble eigvecs_trim(len, MAX_EIGEN, Integral::FULL);
  for (long i = 0; i < MAX_EIGEN; i++) {
    eigvecs.getRow(vec_tmp, i);
    eigvecs_trim.setColumn(i, vec_tmp);
  }
  
  // extract new features for training set
  //
  features_a.mult(K_n, eigvecs_trim);

  printf("extracting testing features...\n");

  // check whether we need extract test features
  //
  if (test_patterns_a.length() <= 0) {
    return test_features_a.setDimensions(0, 0, false);
  }

  // loop through the input test feature vectors and compute the K_test matrix
  //
  long test_num = test_patterns_a.length();
  MatrixDouble K_test(test_num, len, Integral::FULL);

  if (!kernel(K_test, patterns_a, test_patterns_a)) {
    Integral::exit();
  }

  
  // centering matrix K_test in the feature space
  //
  MatrixDouble mat_out_test(test_num, len, Integral::FULL);
  MatrixDouble mat_tmp_test(test_num, len, Integral::FULL);
  MatrixDouble unit_test(test_num, len, Integral::FULL);
  unit_test.assign((double)1 / len);

  // K_test * unit
  //
  mat_tmp_test.mult(K_test, unit);

  // K_test_n = K_test - K_test * unit
  //
  mat_out_test.sub(K_test, mat_tmp_test);

  // unit_test * K
  //
  mat_tmp_test.mult(unit_test, K);

  // K_test_n = K_test - unit_test * K - K_test * unit
  //
  mat_out_test.sub(mat_tmp_test);

  // unit_test * K * unit
  //
  mat_tmp_test.mult(unit);

  // K_test_n = K_test - unit_test * K - K_test * unit + unit_test * K * unit
  //
  mat_out_test.add(mat_tmp_test);

  // output K_test_n
  //
  MatrixDouble K_test_n(mat_out_test);

  // extract new features for test set
  //
  test_features_a.mult(K_test_n, eigvecs_trim);
  
  // write features to data file for evaluation
  //
  k_filename.assign(KFILE);
  k_file.open(k_filename, File::WRITE_ONLY, File::BINARY);

  Long temp;
  temp = len;
  temp.write(k_file, (long)0);
  patterns_a.write(k_file, (long)1);
  K.write(k_file, (long)2);
  eigvecs_trim.write(k_file, (long)3);
  
  // exit gracefully
  //
  return true;
}

// method: computeEigen
//
//  VectorDouble& eigenvals: (output) output eigenvalues
//  MatrixDouble& eigenvecs: (output) output eigenvectors
//  const MatrixDouble& input: (input) symmetric source matrix
//
// return: a boolean value indicating status
//
// this method performs Jacobi eigenvalue/vector computation on symmetric
// matrix
//
boolean computeEigen(VectorDouble& eigenvals_a, MatrixDouble& eigenvecs_a,
		     const MatrixDouble& input_a) {

  // declare local variables
  //
  long n = input_a.getNumRows();
  double **a;
  double *d;
  double **v;
  long nrot = 0;

  a = new double*[n + 1];
  for (long i = 0; i <= n; i++) {
    a[i] = new double[n + 1];
  }

  d = new double[n + 1];

  v = new double*[n + 1];
  for (long i = 0; i <= n; i++) {
    v[i] = new double[n + 1];
  }

  for (long i = 0; i < n; i++) {
    for (long j = 0; j < n; j++) {
      a[i + 1][j + 1] = input_a(i, j);
    }
  }

#define ROTATE(a, i, j, k, l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
a[k][l]=h+s*(g-h*tau);

  long j, iq, ip, i;
  double tresh, theta, tau, t, sm, s, h, g, c, *b, *z;

  b = new double[n + 1];
  z = new double[n + 1];

  for (ip = 1; ip <= n; ip++) {
    for (iq = 1; iq <= n; iq++) {
      v[ip][iq] = 0;
    }
    v[ip][ip] = 1;
  }

  for (ip = 1; ip <= n; ip++) {
    b[ip] = d[ip] = a[ip][ip];
    z[ip] = 0;
  }

  nrot = 0;

  for (i = 1; i <= 50; i++) {
    sm = 0;
    for (ip = 1; ip <= n - 1; ip++) {
      for (iq = ip + 1; iq <= n; iq++) {
	sm += Integral::abs(a[ip][iq]);
      }
    }
    if (sm == 0) {

      // abnormal exit
      //
      goto ab_exit;
    }
    if (i < 4) {
      tresh = 0.2 * sm / (n * n);
    }
    else {
      tresh = 0;
    }
    for (ip = 1; ip <= n - 1; ip++) {
      for (iq = ip + 1; iq <= n; iq++) {
	g = Integral::abs(a[ip][iq]) * 100;
	if ((i > 4) && ((double)(Integral::abs(d[ip]) + g) == (double)Integral::abs(d[ip])) && ((double)(Integral::abs(d[iq]) + g) == (double)Integral::abs(d[iq]))) {
	  a[ip][iq] = 0;
	}
	else if (Integral::abs(a[ip][iq]) > tresh) {
	  h = d[iq] - d[ip];
	  if ((double)(Integral::abs(h) + g) == (double)Integral::abs(h)) {
	    t = a[ip][iq] / h;
	  }
	  else {
	    theta = 0.5 * h / a[ip][iq];
	    t = 1.0 / (Integral::abs(theta) + sqrt(1.0 + theta * theta));
	    if (theta < 0) {
	      t = -t;
	    }
	  }
	  c = 1.0 / sqrt(1 + t * t);
	  s = t * c;
	  tau = s / (1.0 + c);
	  h = t * a[ip][iq];
	  z[ip] -= h;
	  z[iq] += h;
	  d[ip] -= h;
	  d[iq] += h;
	  a[ip][iq] = 0;
	  for (j = 1; j <= ip - 1; j++) {
	    ROTATE(a, j, ip, j, iq);
	  }
	  for (j = ip + 1; j <= iq - 1; j++) {
	    ROTATE(a, ip, j, j, iq);
	  }
	  for (j = iq + 1; j <= n; j++) {
	    ROTATE(a, ip, j, iq, j);
	  }
	  for (j = 1; j <= n; j++) {
	    ROTATE(v, j, ip, j, iq);
	  }

	  ++nrot;
	}
      }
    }
    for (ip = 1; ip <= n; ip++) {
      b[ip] += z[ip];
      d[ip] = b[ip];
      z[ip] = 0;
    }
  }
  printf("error: too many ieterations");

 ab_exit:

  // sort eigenvalues
  //
  long k;
  double p;
  for (i = 1; i < n; i++) {
    p = d[k = i];
    for (j = i + 1; j <= n; j++) {
      if (d[j] >= p) {
	p = d[k = j];
      }
    }
    if (k != i) {
      d[k] = d[i];
      d[i] = p;
      for (j = 1; j <= n; j++) {
	p = v[j][i];
	v[j][i] = v[j][k];
	v[j][k] = p;
      }
    }
  }
  
  eigenvals_a.setLength(n);
  eigenvecs_a.setDimensions(n, n, false, Integral::FULL);

  for (i = 0; i < n; i++) {
    eigenvals_a(i) = d[i + 1];
    for (j = 0; j < n; j++) {
      eigenvecs_a.setValue(i, j, v[j + 1][i + 1]);
    }
  }
  
  delete z;
  delete b;

  for (i = 0; i <= n; i++) {
    delete v[i];
    delete a[i];
  }

  delete v;
  delete a;

  delete d;

  return true;
}
	      
// method: chooseEigenNum
//
//  const VectorDouble& eigenvals: (output) output eigenvalues
//
// return: a long value indicating the effective eigen values
//
// this method chooses effective eigen values
//
long chooseEigenNum(const VectorDouble& eigenvals) {

  // declare local variables
  //
  long len = eigenvals.length();
  Double sum = eigenvals.sum();
  Double thresh = sum * EIGEN_RATIO;
  
  for (long i = len - 1; i >= 0; i--) {

    // sub the eigen value from the sum
    //
    sum -= eigenvals(i);
    if (sum <= thresh) {

      // those eigen values before(including) the current value must
      // be important
      //
      return (i + 1);
    }
  }

  // exit gracefully
  //
  return 1;
}
int main(int argc, char** argv) {

  if (EVAL_ID == 0) {
    CLASS_NUM = 11;
  }
  else {
    CLASS_NUM = 5;
  }
  
  Vector< Vector<VectorDouble> > patterns_buf;
  patterns_buf.setLength(CLASS_NUM);

  for (long i = 0; i < CLASS_NUM; i++) {
    patterns_buf(i).setCapacity(TRAIN_NUM);
    patterns_buf(i).setLength(0);
  }

  // check the parameters
  //
  if (argc < 2) {
    Console::put(L"error on usage");
    Integral::exit();
  }

  else if (argc == 2) {
    return eval(argv[1]);
  }

  FILE* fp_ref = fopen(argv[1], "r");

  char* buf_ref;
  buf_ref = new char[999];
  if (!buf_ref) {
    printf("error in memory!");
    Integral::exit();
  }
    
  memset(buf_ref, 999, 0);
  char buf_num[8];

  char* pref;
  char* pnum;

  int num_data = 0;

  Double norm_val;
  norm_val = 0;

  while (fgets(buf_ref, 999, fp_ref) != (char*)NULL) {

    // get the chars before ':'
    //
    pref = buf_ref;
    pnum = buf_num;
    while (*pref != ':') {
      *pnum++ = *pref++;
    }
    *pnum = 0;
    pref++;

    // try to find out the PATTERN_LEN
    //
    if (PATTERN_LEN <= 0) {
      char* pCount = pref;
      PATTERN_LEN = 0;
      while(*pCount) {
	if (*pCount++ == '.') {
	  PATTERN_LEN++;
	}
      }

      // output the pattern length
      //
      printf("The pattern length = %ld\n", PATTERN_LEN);

      if (PATTERN_LEN == 10) {
	EVAL_ID = 0;
	CLASS_NUM = 11;
      }

      else if (PATTERN_LEN == 39) {
	EVAL_ID = 1;
	CLASS_NUM = 5;
      }

      else {
	delete buf_ref;
	printf("data file %s is wrong\n", argv[1]);
	return false;
      }
    }
    
    // get the patterns class index
    //
    int class_index;
    class_index = atoi(buf_num) - 1;

    // check the index
    //
    if ((class_index < 0) || (class_index >= CLASS_NUM)) {
      printf("error index %d\n", class_index);
      delete buf_ref;
      Integral::exit();
    }

    // resize the patterns_buf buffer
    //
    long len = patterns_buf(class_index).length();
    patterns_buf(class_index).setLength(len + 1, true);
    patterns_buf(class_index)(len).setLength(PATTERN_LEN);

    double temp;
    for (long i = 0; i < PATTERN_LEN; i++) {
      sscanf(pref, " %lf", &temp);
      patterns_buf(class_index)(len)(i) = temp;
      while(*pref == ' ') {
	pref++;
      }
      while(*pref != ' ') {
	pref++;
      }
    }

    num_data++;
    memset(buf_ref, 999, 0);
  }

  patterns_buf.setLength(CLASS_NUM, true);

  printf("Total %d training tokens\n", num_data);
  fclose(fp_ref);

  if (EVAL_ID == 1) {
    num_data /= 5;
    
    for (long i = 0; i < patterns_buf.length(); i++) {
      if (!merge(patterns_buf(i))) {
	printf("error in merge!");
	break;
      }
    }
  }

  // copy the patterns data from Vector< Vector<VectorDouble> > to
  // Vector<VectorFloat>
  //
  Vector<VectorDouble> patterns;
  patterns.setLength(num_data);
  long pat_index = 0;
  for (long i = 0; i < patterns_buf.length(); i++) {
    for (long j = 0; j < patterns_buf(i).length(); j++) {
      //      patterns(pat_index++).assign(patterns_buf(i)(j));
      patterns(pat_index).assign(patterns_buf(i)(j));
      double sum = patterns(pat_index++).sumSquare();
      if (sum > norm_val) {
	norm_val = sum;
      }
    }
  }

  // normalize patterns
  //
  norm_val.sqrt();
  norm_val.assign((double)1/norm_val);
  for (long i = 0; i < patterns.length(); i++) {
    patterns(i).mult(norm_val);
  }

  // output normalization factor
  //
  norm_val.debug(L"norm_out");
  
  // read test patterns
  //
  Vector<VectorDouble> test_patterns;
  fp_ref = fopen(argv[2], "r");

  memset(buf_ref, 999, 0);
  num_data = 0;
  while (fgets(buf_ref, 999, fp_ref) != (char*)NULL) {

    // get the chars before ':'
    //
    pref = buf_ref;
    pnum = buf_num;
    while (*pref != ':') {
      *pnum++ = *pref++;
    }
    *pnum = 0;
    pref++;

    // resize the test patterns buffer
    //
    long len = test_patterns.length();
    test_patterns.setLength(len + 1, true);
    test_patterns(len).setLength(PATTERN_LEN);

    double temp;
    for (long i = 0; i < PATTERN_LEN; i++) {
      sscanf(pref, " %lf", &temp);
      test_patterns(len)(i) = temp;
      while(*pref == ' ') {
	pref++;
      }
      while(*pref != ' ') {
	pref++;
      }
    }

    num_data++;
    memset(buf_ref, 999, 0);
  }

  if (EVAL_ID == 1) {
    num_data /= 5;
    
    if (!merge(test_patterns)) {
      printf("error in merge!");
    }
  }

  // normalize test patterns
  //
  for (long i = 0; i < test_patterns.length(); i++) {
    test_patterns(i).mult(norm_val);
  }
  
  delete buf_ref;
  buf_ref = 0;
  printf("Total %d testing tokens\n", num_data);
  fclose(fp_ref);

  printf("extracting all features...\n");

  // compute train features and test features
  //
  MatrixDouble features, test_features;
  if (!compute(features, test_features, patterns, test_patterns)) {
    Integral::exit();
  }

  printf("computing on pseudoinverse matrix...\n");

  // compute pseudoinverse matrix
  //
  long feature_len = features.getNumColumns();
  long feature_num = features.getNumRows();
  MatrixDouble Y(feature_num, feature_len + 1, Integral::FULL);
  VectorDouble vec_b(feature_num);

  for (long i = 0; i < feature_len; i++) {

    features.getColumn(vec_b, i);
    Y.setColumn(i + 1, vec_b);
  }

  vec_b.assign((double)1);
  Y.setColumn(0, vec_b);

  MatrixDouble Yt;
  Yt.transpose(Y);

  MatrixDouble Ytemp;
  Ytemp.mult(Yt, Y);
  Y.inverse(Ytemp);
  Y.mult(Yt);

  printf("computing on %ld classifiers...\n", (long)CLASS_NUM);

  // compute the matrix for classifiers (column-based)
  //
  MatrixDouble a(feature_len + 1, (long)CLASS_NUM, Integral::FULL);
  pat_index = 0;
  VectorDouble vec_tmp;
  for (long i = 0; i < patterns_buf.length(); i++) {
    vec_b.assign((double)(MARGIN_L));
    for (long j = 0; j < patterns_buf(i).length(); j++) {
      vec_b(pat_index++).assign((double)(MARGIN_H));
    }

    Y.multv(vec_tmp, vec_b);
    a.setColumn(i, vec_tmp);
  }

  // if there is no test patterns, exit
  //
  if (test_patterns.length() <= 0) {
    Integral::exit();
  }
  
  printf("computing on scores...\n");

  // compute the scores for test features
  //
  long test_feature_num = test_features.getNumRows();
  Yt.setDimensions(test_feature_num, feature_len + 1, false);
  vec_tmp.setLength(test_feature_num);
  vec_tmp.assign((double)1);
  Yt.setColumn(0, vec_tmp);

  // copy test features to Yt buffer
  //
  for (long i = 0; i < feature_len; i++) {
    test_features.getColumn(vec_tmp, i);
    Yt.setColumn(i + 1, vec_tmp);
  }

  // compute scores
  //
  Yt.mult(a);

  // count the highest score to determine the class
  //
  VectorLong vec_class(test_feature_num);
  long max_index;
  double max_val;
  double cur_val;
  for (long i = 0; i < test_feature_num; i++) {
    max_index = 0;
    max_val = Yt(i, 0);
    for (long j = 1; j < CLASS_NUM; j++) {
      cur_val = Yt(i, j);
      if (max_val < cur_val) {
	max_val = cur_val;
	max_index = j;
      }
    }
    vec_class(i) = max_index + 1;
  }

  printf("successful classification...\n");

  // write the result to data file
  //
  fp_ref = fopen(OUTFILE, "w");

  if (EVAL_ID == 0) {
    for (long i = 0; i < test_feature_num; i++) {
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
    }
  }
  else {
    for (long i = 0; i < test_feature_num; i++) {
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
      fprintf(fp_ref, "%ld\n", (long)vec_class(i));
    }
  }
  fclose(fp_ref);

  // write the features to data file for evaluation
  //
  a.write(k_file, (long)4);
  norm_val.write(k_file, (long)5);
  k_file.close();
  
  if (buf_ref) {
    delete buf_ref;
  }

  printf("successfully writing result to %s...\n", OUTFILE);
  
  // exit gracefully
  //
  Integral::exit();
}