// ----------------------------------------------------------------------------
//
// Copyright (C) Columbia University, 1994. All Rights Reserved.
// Sameer A. Nene, Shree K. Nayar, Hiroshi Murase
//
// See file LICENSE for details of software license agreement.
//
// ----------------------------------------------------------------------------
//
// eigen.C
//
// Author:               Sameer Nene
// Date:                 06/07/94
// Version:              1.0.0
// Modification History:
//
// Classes:
//
// Notes:
//   This command is used for computing eigen vectors from vector sets. See the
//  man page for more details.
//   
// ----------------------------------------------------------------------------

#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include "errorscope.h"
#include "vector.h"
#include "matrix.h"
#include "parser.h"
#include "fileiter.h"
#include "ioutil.h"
#include "reporter.h"
#include "stringary.h"

int main(char argc, char **argv)
{
  int half, i, j, k, l, n, size;
  double t;
  Vector v1, v2;
  CommandParser cp(argc, argv, "a:e:l:m:n:AD:LV");
  ErrorScope::Error error;

  if((error = cp.error()) != ErrorScope::OK) {
    cerr << error << endl << "Exiting" << endl;
    exit(1);
  }

  if(cp.getVerboseFlag() != 0)
    cerr << "Copyright \xA9 Columbia University, 1994. All Rights Reserved." << endl << "Sameer Nene, Shree Nayar, Hiroshi Murase" << endl << endl;
  
  FileListAry fa(cp.getObjectDirectory(), cp.getObjectList(), "vec");
  if((error = fa.error()) != ErrorScope::OK) {
    cerr << error << endl << "Exiting" << endl;
    exit(1);
  }

  j = fa.getSize();
  
  if(cp.getAverageFlag() != 0) {
    if(cp.getVerboseFlag() != 0)
      cerr << "Computing Average" << endl << endl;
    for(i = 0, k = 0; i < j; ++i)
      for(FileListIter it(fa[i]); it; ++it, ++k) {
	if(cp.getVerboseFlag() != 0)
	  Reporter::display_status();
	if((error = IOUtil::get(&v2, it())) != ErrorScope::OK) {
	  cerr << '\b' << error << ": " << it() << endl << "Exiting" << endl;
	  exit(1);
	}
	if(!k)
	  v1 = v2;
	else
	  v1 += v2;
      }
  }
  
  if(cp.getAverageFlag() != 0) {
    v1 /= k;
    if((error = IOUtil::put(v1, cp.getAverageVectorFileName())) != ErrorScope::OK) {
      cerr << '\b' << error << ": " << cp.getAverageVectorFileName() << endl << "Exiting" << endl;
      exit(1);
    }
  }

  for(i = 0, k = 0; i < j; ++i)
    k += fa[i].getSize();
  
  StringAry sa(k);
  for(i = 0, k = 0; i < j; ++i)
    for(FileListIter it(fa[i]); it; ++it, ++k)
      sa.set(k, it());
  
  if(cp.getVerboseFlag() != 0)
    cerr << '\b' << "Computing Implicit Covariance Matrix" << endl << endl;

  size = k;
  Matrix m(size, size);

  {
    if((error = IOUtil::get(&v2, sa[0])) != ErrorScope::OK) {
      cerr << '\b' << error << ": " << sa[0] << endl << "Exiting" << endl;
      exit(1);
    }
    half = ((cp.getMemRequired() * 1048576) >> 1) / (v2.length() * sizeof(double));
    VectAry va(half << 1);
    VectAryIter it1(va), it2(va);
    
    for(i = 0; i < size; i += half) {
      for(l = 0, k = 0, it1.set(0); l < half && i + l < size; ++l, ++k, ++it1) {
	if((error = IOUtil::get(&it1(), sa[i + l])) != ErrorScope::OK) {
	  cerr << '\b' << error << ": " << sa[i + l] << endl << "Exiting" << endl;
	  exit(1);
	}
	if(cp.getAverageFlag() != 0) 
	  it1() -= v1;
      }
      
      for(l = 0, it1.set(0); l < k; ++l, ++it1) {
	if(cp.getVerboseFlag() != 0)
	  Reporter::display_status();
	for(n = l, it2.set(l); n < k; ++n, ++it2) {
	  t = it1() * it2();
	  m.element(i + l, i + n, t);
	  m.element(i + n, i + l, t);
	}
      }
      
      for(j = i + half; j < size; j += half) {
	for(l = 0, k = 0, it1.set(half); l < half && j + l < size; ++l, ++k, ++it1) {
	  if((error = IOUtil::get(&it1(), sa[j + l])) != ErrorScope::OK) {
	    cerr << '\b' << error << ": " << sa[j + l] << endl << "Exiting" << endl;
	    exit(1);
	  }
	  if(cp.getAverageFlag() != 0) 
	    it1() -= v1;
	}
	
	for(l = 0, it1.set(0); l < half; ++l, ++it1) {
	  if(cp.getVerboseFlag() != 0)
	    Reporter::display_status();
	  for(n = 0, it2.set(half); n < k; ++n, ++it2) {
	    t = it1() * it2();
	    m.element(i + l, j + n, t);
	    m.element(j + n, i + l, t);
	  }
	}
      }
    }
  }

  if(cp.getVerboseFlag() != 0)
    cerr << '\b' << "Computing Implicit Eigen Vectors" << endl << endl;

  Vector ev(cp.getNumEigenVectors());
  VectAry va(cp.getNumEigenVectors()), vb(cp.getNumEigenVectors()),
  vc(n = half << 1);
  VectAryIter ita(va), itb(vb), itc(vc);
  
  MatrixUtil::getMaxEigenVectors(&m, &va, &ev);

  if(cp.getEigenValueFlag() != 0)
    if((error = IOUtil::put(ev, cp.getEigenValueFileName())) != ErrorScope::OK) {
      cerr << '\b' << error << ": " << cp.getEigenValueFileName() << endl << "Exiting" << endl;
      exit(1);
    }    
  
  if(cp.getVerboseFlag() != 0)
    cerr << '\b' << "Computing Eigen Vectors" << endl << endl;

  for(i = 0; i < size; i += n) {
    for(j = 0, itc.set(0); j < n && i + j < size; ++j, ++itc) {
      if((error = IOUtil::get(&itc(), sa[i + j])) != ErrorScope::OK) {
	cerr << '\b' << error << ": " << sa[i + j] << endl << "Exiting" << endl;
	exit(1);
      }
      if(cp.getAverageFlag() != 0) 
	itc() -= v1;
    }
    
    for(ita.set(0), itb.set(0); ita; ++ita, ++itb) {
      if(cp.getVerboseFlag() != 0)
	Reporter::display_status();
      for(j = 0, itc.set(0); j < n && i + j < size; ++j, ++itc)
	if(i == 0 && j == 0)
	  itb() = itc() * ita().element(i);
	else
	  itb() += itc() * ita().element(i + j);
    }
  }
  
  for(itb.set(0); itb; ++itb)
    itb() /= itb()();
  
  if((error = IOUtil::put(vb, cp.getEigenVectorFileName())) != ErrorScope::OK) {
    cerr << '\b' << error << ": " << cp.getEigenVectorFileName() << endl << "Exiting" << endl;
    exit(1);
  }

  if(cp.getVerboseFlag() != 0)
    cerr << '\b' << "Done" << endl;

  return 0;
}
