#include "GLCP.h"
#include "dbimage.h"
#include "bdbenv.h"
#include "dbalnimage.h"
#include "dbtmpfeature.h"
#include <tclap/CmdLine.h>
#include "emcluster.h"
#include "kldivergence.h"
#include <string>

int calcmean( std::vector< double *> &vecs, 
              std::vector< double > &mean );

int calcvar( std::vector< double *> &vecs, 
             std::vector< double > &mean,
             std::vector< double > &var );

int calcjsd( std::vector< double *> &vecsa, 
             std::vector< double *> &vecsb, 
             std::vector< double > &var,
             int param);

int calcmean( std::vector< double *> &vecs, 
              std::vector< double > &mean )
{
  for( int j = 0; j < mean.size(); j++ )
    mean[j] = 0;

  for( int i = 0; i < vecs.size(); i++ )
    for( int j = 0; j < mean.size(); j++ )
      mean[j] = mean[j] + vecs[i][j];
  
  for( int j = 0; j < mean.size(); j++ )
    mean[j] = mean[j]/((double) vecs.size());
  
  return 0;
}

int calcvar( std::vector< double *> &vecs,
             std::vector< double > &mean,
             std::vector< double > &var )
{
  for( int j = 0; j < var.size(); j++ )
    var[j] = 0;

  for( int i = 0; i < vecs.size(); i++ )
    for( int j = 0; j < var.size(); j++ )
      var[j] = var[j] + pow( mean[j] - vecs[i][j], 2);
  
  for( int j = 0; j < var.size(); j++ )
    var[j] = var[j]/((double) vecs.size() - 1.0);
  
  return 0;
}

int calcjsd( std::vector< double *> &vecsa, 
             std::vector< double *> &vecsb, 
             std::vector< double > &var,
             int param)
{
  std::vector< double > a, b;
  for( unsigned int i = 0; i < vecsa.size(); i++ )
    a.push_back( vecsa[i][param] );
  for( unsigned int i = 0; i < vecsb.size(); i++ )
    b.push_back( vecsb[i][param] );
  kldivergence jsd;
#ifdef __DEBUG__
  std::cout << "param " << param << " size: " << a.size() << " " << b.size() << "\n";
#endif
  double res;
  if( var.size() <= param )
    {
      std::cout << "var size error\n"; 
      exit(0);
    }
  res = jsd.calc( a, b, sqrt(var[param])/5.0 );
  std::cout << res << "\n";
  return 0;
}

int main( int argc, char *argv[] )
{

#ifndef _NOMAGICK_
  try
    {
      ///set up the command line parsing
      TCLAP::CmdLine cmd("Image alignment program", ' ', "0.001");

      TCLAP::ValueArg<std::string> dbdir("d", "db", "database dirctory", true, "", "databaase");
      TCLAP::ValueArg<std::string> tdbg("g", "good", "good feature db name", true,
                                        "", "file");
      TCLAP::ValueArg<std::string> tdbb("b", "bad", "bad feature db name", true,
                                        "", "file");


      cmd.add(dbdir);
      cmd.add( tdbg );
      cmd.add( tdbb );

      bdbenv bdb;
      dbtmpfeature dbtf_g, dbtf_b;

      cmd.parse(argc, argv);

      bdb.open( dbdir.getValue() );

      dbtf_g.setname( tdbg.getValue() );
      dbtf_b.setname( tdbb.getValue() );
      dbtf_g.init( bdb.envptr(), 0 );
      dbtf_b.init( bdb.envptr(), 0 );

      bdb.adddb( &dbtf_g );
      bdb.adddb( &dbtf_b );
      bdb.cleanlog();

      std::vector< double * > vecsg, vecsb, vecst;
      unsigned int maxid;

      maxid = dbtf_g.maxid();
      int dbrecsz;
      dbtf_g.getrecsize( dbrecsz );
      std::cout << "record size: " << dbrecsz << "\n";
      char *buffer;
      buffer = new char[dbrecsz];
      GLCP myglcp;
      int numfeatures;
      for( unsigned i = 1; i < maxid; i++ )
        {
          dbtf_g.get( i, buffer );
          myglcp.deserialize( buffer, dbrecsz );
          std::vector< double * > ttt;
          myglcp.getfeaturets( ttt, numfeatures );
          for( unsigned j = 0; j < ttt.size(); j++ )
            vecsg.push_back( ttt[j] );
          ttt.clear();
        }

      maxid = dbtf_b.maxid();
      for( unsigned i = 1; i < maxid; i++ )       
        {
          dbtf_b.get( i, buffer );
          myglcp.deserialize( buffer, dbrecsz );
          std::vector< double * > ttt;
          myglcp.getfeaturets( ttt, numfeatures );
          for( unsigned j = 0; j < ttt.size(); j++ )
            vecsb.push_back( ttt[j] );
          ttt.clear();
        }
      delete [] buffer;

      dbtf_g.close();
      dbtf_b.close();
      bdb.close();

      for( unsigned i = 0; i < vecsg.size(); i++ )
        vecst.push_back( vecsg[i] );
      for( unsigned i = 0; i < vecsb.size(); i++ )
        vecst.push_back( vecsb[i] );

      std::vector< double > mean, var;
      std::cout << "number of features " << numfeatures << "\n";
      mean.resize( numfeatures );
      var.resize( numfeatures );

      calcmean( vecst, mean );
      calcvar( vecst, mean, var );

      for( int i = 0; i < numfeatures; i++ )
        std::cout << mean[i] << " \t";
      std::cout << "\n";
      for( int i = 0; i < numfeatures; i++ )
        std::cout << var[i] << " \t";
      std::cout << "\n";

      for( int i = 0; i < numfeatures; i++ )
        calcjsd( vecsg, vecsb, var, i );

      for( unsigned i = 0; i < vecst.size(); i++ )
        delete [] vecst[i];
    }
  catch (TCLAP::ArgException &e)
    {
      std::cout << "ERROR: " << e.error() 
                << " for arg " << e.argId() << "\n";
      return 0;
    }
  
#else
  std::cout << "compiled with no image magick\n";
#endif
  return 0;
}
