#include "psotracker.h"

int bdb_transcache_itr_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_transcache_slice_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_transcache_size_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_transcache_hit_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);

int bdb_mi_score_slice_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_mi_score_xy_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_mi_score_score_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);

int bdb_particle_setup_psoslice_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);

int bdb_particle_state_vx_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_particle_state_vy_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_particle_state_speed_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_particle_state_iter_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);
int bdb_particle_state_pid_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey);

int bdb_int_cmp(Db *dbp, const Dbt *a, const Dbt *b);
int bdb_double_cmp(Db *dbp, const Dbt *a, const Dbt *b);

/// custom int compare
int bdb_int_cmp(Db *dbp, const Dbt *a, const Dbt *b)
{
  dbp = NULL;
  int ai, bi;
  // Returns: 
  // < 0 if a < b 
  // = 0 if a = b 
  // > 0 if a > b 
  memcpy(&ai, a->get_data(), sizeof(int)); 
  memcpy(&bi, b->get_data(), sizeof(int)); 
  return (ai - bi); 
}

/// custom double compare
int bdb_double_cmp(Db *dbp, const Dbt *a, const Dbt *b)
{
  dbp = NULL;
  double ai, bi;
  // Returns: 
  // < 0 if a < b 
  // = 0 if a = b 
  // > 0 if a > b 
  memcpy(&ai, a->get_data(), sizeof(double)); 
  memcpy(&bi, b->get_data(), sizeof(double)); 
  double res;
  res = ai - bi;
  if(res < 0 )
    return -1;
  else if( res > 0 )
    return 1;
  else
    return 0;
}

int bdb_transcache_itr_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::trans_cache_key *tck;
  tck = (psotracker::trans_cache_key *)(pkey->get_data());
  skey->set_size( sizeof( int ) );
  skey->set_data( &(tck->iteration) );
  sdbp = NULL;
  pdata = NULL;
  return 0;
}

int bdb_transcache_slice_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::trans_cache_key *tck;
  tck = (psotracker::trans_cache_key *)(pkey->get_data());
  skey->set_size( sizeof( int ) );
  skey->set_data( &(tck->psoslice) );
  sdbp = NULL;
  pdata = NULL;  
  return 0;
}

int bdb_transcache_size_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::trans_cache_value *tcv;
  tcv = (psotracker::trans_cache_value *)(pdata->get_data());
  skey->set_size( sizeof( int ) );
  skey->set_data( &(tcv->size) );
  sdbp = NULL;
  pkey = NULL;
  return 0;
}

int bdb_transcache_hit_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::trans_cache_value *tcv;
  tcv = (psotracker::trans_cache_value *)(pdata->get_data());
  skey->set_size( sizeof( double ) );
  skey->set_data( &(tcv->hit) );
  sdbp = NULL;
  pkey = NULL;
  return 0;
}

int bdb_mi_score_slice_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::mi_score_key *msk;
  msk = (psotracker::mi_score_key *)(pkey->get_data());
  skey->set_size( sizeof( int ) );
  skey->set_data( &(msk->slice) );
  sdbp = NULL;
  pdata = NULL;
  return 0;
}

int bdb_mi_score_xy_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::mi_score_key *msk;
  msk = (psotracker::mi_score_key *)(pkey->get_data());
  skey->set_size( 2*sizeof( int ) );
  skey->set_data( &(msk->x) );
  sdbp = NULL;
  pdata = NULL;
  return 0;
}

int bdb_mi_score_score_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  double *msk;
  msk = (double *)(pdata->get_data());
  skey->set_size( sizeof( double ) );
  skey->set_data( msk );
  sdbp = NULL;
  pkey = NULL;
  return 0;
}

int bdb_particle_setup_psoslice_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::particle_setup_value *msk;
  msk = (psotracker::particle_setup_value *)(pkey->get_data());
  skey->set_size( sizeof( int ) );
  skey->set_data( &(msk->psoslice) );
  sdbp = NULL;
  pdata = NULL;
  return 0;
}

int bdb_particle_state_vx_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::particle_state_value *msk;
  msk = (psotracker::particle_state_value *)(pdata->get_data());
  skey->set_size( sizeof( double ) );
  skey->set_data( &(msk->vx) );
  sdbp = NULL;
  pkey = NULL;
  return 0;
}

int bdb_particle_state_vy_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::particle_state_value *msk;
  msk = (psotracker::particle_state_value *)(pdata->get_data());
  skey->set_size( sizeof( double ) );
  skey->set_data( &(msk->vy) );
  sdbp = NULL;
  pkey = NULL;
  return 0;
}

int bdb_particle_state_speed_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::particle_state_value *msk;
  msk = (psotracker::particle_state_value *)(pdata->get_data());
  skey->set_size( sizeof( double ) );
  skey->set_data( &(msk->speed) );
  sdbp = NULL;
  pkey = NULL;
  return 0;
}

int bdb_particle_state_iter_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::particle_state_key *msk;
  msk = (psotracker::particle_state_key *)(pkey->get_data());
  skey->set_size( sizeof( int ) );
  skey->set_data( &(msk->iteration) );
  sdbp = NULL;
  pdata = NULL;
  return 0;
}

int bdb_particle_state_pid_extract(Db *sdbp, const Dbt *pkey, const Dbt *pdata, Dbt *skey)
{
  psotracker::particle_state_key *msk;
  msk = (psotracker::particle_state_key *)(pkey->get_data());
  skey->set_size( sizeof( int ) );
  skey->set_data( &(msk->pid) );
  sdbp = NULL;
  pdata = NULL;
  return 0;
}

psotracker::psotracker()
{
  pso_paramdb = NULL;
  slice_setupdb = NULL;

  trans_cachedb = NULL;        ///< db for storing cache
  trans_cache_iteridx = NULL;
  trans_cache_sliceidx = NULL;
  trans_cache_sizeidx = NULL;
  trans_cache_hitidx = NULL;

  mi_scoredb = NULL;
  mi_score_sliceidx = NULL;
  mi_score_xyidx = NULL;
  mi_score_scoreidx = NULL;

  particle_setupdb = NULL;
  particle_setup_psosliceidx = NULL;

  particle_statedb = NULL;
  particle_state_vxidx = NULL;
  particle_state_vyidx = NULL;
  particle_state_speedidx = NULL;
  particle_state_iteridx = NULL;
  particle_state_pididx = NULL;

  iteration_statdb = NULL;

  cluster_positiondb = NULL;
}

psotracker::~psotracker()
{
  close();
}

int psotracker::close()
{
  if( pso_paramdb != NULL )
    {
      //sync
      pso_paramdb->sync(0);
      slice_setupdb->sync(0);
      
      trans_cache_iteridx->sync(0);
      trans_cache_sliceidx->sync(0);
      trans_cache_sizeidx->sync(0);
      trans_cache_hitidx->sync(0);
      trans_cachedb->sync(0);
      
      mi_score_sliceidx->sync(0);
      mi_score_xyidx->sync(0);
      mi_score_scoreidx->sync(0);
      mi_scoredb->sync(0);

      particle_setup_psosliceidx->sync(0);
      particle_setupdb->sync(0);
      
      particle_state_vxidx->sync(0);
      particle_state_vyidx->sync(0);
      particle_state_speedidx->sync(0);
      particle_state_iteridx->sync(0);
      particle_state_pididx->sync(0);
      particle_statedb->sync(0);

      iteration_statdb->sync(0);

      cluster_positiondb->sync(0);

      // close
      pso_paramdb->close(0);
      slice_setupdb->close(0);
      
      trans_cache_iteridx->close(0);
      trans_cache_sliceidx->close(0);
      trans_cache_sizeidx->close(0);
      trans_cache_hitidx->close(0);
      trans_cachedb->close(0);
      
      mi_score_sliceidx->close(0);
      mi_score_xyidx->close(0);
      mi_score_scoreidx->close(0);
      mi_scoredb->close(0);

      particle_setup_psosliceidx->close(0);
      particle_setupdb->close(0);
      
      particle_state_vxidx->close(0);
      particle_state_vyidx->close(0);
      particle_state_speedidx->close(0);
      particle_state_iteridx->close(0);
      particle_state_pididx->close(0);
      particle_statedb->close(0);

      iteration_statdb->close(0);

      cluster_positiondb->close(0);
      // delete the memory
      delete pso_paramdb;
      delete slice_setupdb;
      delete trans_cache_iteridx;
      delete trans_cache_sliceidx;
      delete trans_cache_sizeidx;
      delete trans_cache_hitidx;
      delete trans_cachedb;
      delete mi_score_sliceidx;
      delete mi_score_xyidx;
      delete mi_score_scoreidx;
      delete mi_scoredb;
      delete particle_setup_psosliceidx;
      delete particle_setupdb; 
      delete particle_state_vxidx;
      delete particle_state_vyidx;
      delete particle_state_speedidx;
      delete particle_state_iteridx;
      delete particle_state_pididx;
      delete particle_statedb;
      delete iteration_statdb;
      delete cluster_positiondb;
      bdb.close();
    }

  pso_paramdb = NULL;
  slice_setupdb = NULL;

  trans_cachedb = NULL;        ///< db for storing cache
  trans_cache_iteridx = NULL;
  trans_cache_sliceidx = NULL;
  trans_cache_sizeidx = NULL;
  trans_cache_hitidx = NULL;

  mi_scoredb = NULL;
  mi_score_sliceidx = NULL;
  mi_score_xyidx = NULL;
  mi_score_scoreidx = NULL;

  particle_setupdb = NULL;
  particle_setup_psosliceidx = NULL;

  particle_statedb = NULL;
  particle_state_vxidx = NULL;
  particle_state_vyidx = NULL;
  particle_state_speedidx = NULL;
  particle_state_iteridx = NULL;
  particle_state_pididx = NULL;

  iteration_statdb = NULL;

  cluster_positiondb = NULL;
  return 0;
}

int psotracker::upgrade( std::string dbfile )
{
  dbfilename = dbfile;
  close();
  Db *tdb;
  bdbenv mybdb;
  DbEnv *env;
  std::cout << "upgrading database\n";
  mybdb.open( dbfile );
  env = mybdb.envptr();
  tdb = new Db( env, 0 );
  tdb->upgrade( "pso_param.db", 0 );
  tdb->upgrade( "slice_setup.db", 0 );
  tdb->upgrade( "trans_cache.db", 0 );
  tdb->upgrade( "mi_score.db", 0 );
  tdb->upgrade( "particle_setup.db", 0 );
  tdb->upgrade( "particle_state.db", 0 );
  tdb->upgrade( "iteration_stat.db", 0 );
  tdb->upgrade( "cluster_state.db", 0 );
  tdb->upgrade( "trans_cache_iteridx.db", 0 );
  tdb->upgrade( "trans_cache_sliceidx.db", 0 );
  tdb->upgrade( "trans_cache_sizeidx.db", 0 );
  tdb->upgrade( "trans_cache_hitidx.db", 0 );
  tdb->upgrade( "mi_score_sliceidx.db", 0 );
  tdb->upgrade( "mi_score_xyidx.db", 0 );
  tdb->upgrade( "mi_score_scoreidx.db", 0 );
  tdb->upgrade( "particle_setup_psosliceidx.db", 0 );
  tdb->upgrade( "particle_state_vxidx.db", 0 );
  tdb->upgrade( "particle_state_vyidx.db", 0 );
  tdb->upgrade( "particle_state_speedidx.db", 0 );
  tdb->upgrade( "particle_state_iteridx.db", 0 );
  tdb->upgrade( "particle_state_pididx.db", 0 );
  delete tdb;
  mybdb.close();
  std::cout << "upgrade done\n";
  return 0;
}

int psotracker::init( std::string dbfile, int rebuildindex )
{
  dbfilename = dbfile;
  close();
  bdb.setopenflags(DB_CREATE | DB_THREAD | DB_INIT_MPOOL );
  bdb.setcache( 1, 512 );
  bdb.open( dbfile );
  DbEnv *env;
  env = bdb.envptr();
  // create db objects
  pso_paramdb = new Db( env, 0 );
  slice_setupdb = new Db( env, 0 );
  trans_cachedb = new Db( env, 0 );
  trans_cache_iteridx = new Db( env, 0 );
  trans_cache_sliceidx = new Db( env, 0 );
  trans_cache_sizeidx = new Db( env, 0 );
  trans_cache_hitidx = new Db( env, 0 );
  mi_scoredb = new Db( env, 0 );
  mi_score_sliceidx = new Db( env, 0 );
  mi_score_xyidx = new Db( env, 0 );
  mi_score_scoreidx = new Db( env, 0 );
  particle_setupdb = new Db( env, 0 );
  particle_setup_psosliceidx = new Db( env, 0 );
  particle_statedb = new Db( env, 0 );
  particle_state_vxidx = new Db( env, 0 );
  particle_state_vyidx = new Db( env, 0 );
  particle_state_speedidx = new Db( env, 0 );
  particle_state_iteridx = new Db( env, 0 );
  particle_state_pididx = new Db( env, 0 );
  iteration_statdb = new Db( env, 0 );
  cluster_positiondb = new Db( env, 0 );
  // open databases
  pso_paramdb->open( NULL, "pso_param.db", NULL, 
                     DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  slice_setupdb->set_bt_compare( bdb_int_cmp );
  slice_setupdb->open( NULL, "slice_setup.db", NULL, 
                       DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  trans_cachedb->open( NULL, "trans_cache.db", NULL, 
                       DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  mi_scoredb->open( NULL, "mi_score.db", NULL, 
                    DB_HASH, DB_CREATE | DB_THREAD, 0 );

  
  particle_setupdb->set_bt_compare( bdb_int_cmp );
  particle_setupdb->open( NULL, "particle_setup.db", NULL, 
                          DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  particle_statedb->open( NULL, "particle_state.db", NULL, 
                          DB_HASH, DB_CREATE | DB_THREAD, 0 );

  iteration_statdb->set_bt_compare( bdb_int_cmp );
  iteration_statdb->open( NULL, "iteration_stat.db", NULL, 
                          DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  cluster_positiondb->open( NULL, "cluster_state.db", NULL, 
                            DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  if( rebuildindex )
    {
      bdb.removedb( "trans_cache_iteridx.db" );
      bdb.removedb( "trans_cache_sizeidx.db"); 
      bdb.removedb( "trans_cache_hitidx.db"); 
      bdb.removedb( "mi_score_sliceidx.db"); 
      bdb.removedb( "mi_score_xyidx.db"); 
      bdb.removedb( "mi_score_scoreidx.db"); 
      bdb.removedb( "particle_setup_psosliceidx.db"); 
      bdb.removedb( "particle_state_vxidx.db"); 
      bdb.removedb( "particle_state_vyidx.db"); 
      bdb.removedb( "particle_state_speedidx.db"); 
      bdb.removedb( "particle_state_iteridx.db"); 
      bdb.removedb( "particle_state_pididx.db");
    }
  // open secondary index

  // trans_cache secondary databases
  trans_cache_iteridx->set_flags( DB_DUP );
  trans_cache_iteridx->set_bt_compare( bdb_int_cmp );
  trans_cache_iteridx->open( NULL, "trans_cache_iteridx.db", NULL, 
                             DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  trans_cache_sliceidx->set_flags( DB_DUP );
  trans_cache_sliceidx->set_bt_compare( bdb_int_cmp );
  trans_cache_sliceidx->open( NULL, "trans_cache_sliceidx.db", NULL,
                              DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  trans_cache_sizeidx->set_flags( DB_DUP );
  trans_cache_sizeidx->set_bt_compare( bdb_int_cmp );
  trans_cache_sizeidx->open( NULL, "trans_cache_sizeidx.db", NULL,
                             DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  trans_cache_hitidx->set_flags( DB_DUP );
  trans_cache_hitidx->set_bt_compare( bdb_double_cmp );
  trans_cache_hitidx->open( NULL, "trans_cache_hitidx.db", NULL,
                            DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  //////////////////////////////////////////////////////

  // mi_scoredb secondary databases
  mi_score_sliceidx->set_flags( DB_DUP );
  mi_score_sliceidx->set_bt_compare( bdb_int_cmp );
  mi_score_sliceidx->open( NULL, "mi_score_sliceidx.db", NULL,
                           DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  mi_score_xyidx->set_flags( DB_DUP );
  mi_score_xyidx->open( NULL, "mi_score_xyidx.db", NULL,
                        DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  mi_score_scoreidx->set_flags( DB_DUP );
  mi_score_scoreidx->set_bt_compare( bdb_double_cmp );
  mi_score_scoreidx->open( NULL, "mi_score_scoreidx.db", NULL,
                           DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  //////////////////////////////////////////////////////

  // particle_setupdb secondary databases
  particle_setup_psosliceidx->set_flags( DB_DUP );
  particle_setup_psosliceidx->set_bt_compare( bdb_int_cmp );
  particle_setup_psosliceidx->open( NULL, "particle_setup_psosliceidx.db", NULL,
                                    DB_BTREE, DB_CREATE | DB_THREAD, 0 );

  //////////////////////////////////////////////////////

  // particle_statdb secondary databases
  particle_state_vxidx->set_flags( DB_DUP );
  particle_state_vxidx->set_bt_compare( bdb_double_cmp );
  particle_state_vxidx->open( NULL, "particle_state_vxidx.db", NULL,
                              DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  particle_state_vyidx->set_flags( DB_DUP );
  particle_state_vyidx->set_bt_compare( bdb_double_cmp );
  particle_state_vyidx->open( NULL, "particle_state_vyidx.db", NULL,
                              DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  particle_state_speedidx->set_flags( DB_DUP );
  particle_state_speedidx->set_bt_compare( bdb_double_cmp );
  particle_state_speedidx->open( NULL, "particle_state_speedidx.db", NULL,
                              DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  particle_state_iteridx->set_flags( DB_DUP);
  particle_state_iteridx->set_bt_compare( bdb_int_cmp );
  particle_state_iteridx->open( NULL, "particle_state_iteridx.db", NULL,
                              DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  particle_state_pididx->set_flags( DB_DUP);
  particle_state_pididx->set_bt_compare( bdb_int_cmp );
  particle_state_pididx->open( NULL, "particle_state_pididx.db", NULL,
                              DB_BTREE, DB_CREATE | DB_THREAD, 0 );
  
  ////////////////////////////////////////////
#ifdef __PSOTRACKER_SECONDARYDB__
  DbTxn *dbtxn;
  dbtxn = NULL;
  //  dbtxn = bdb.begintransaction();
  std::cout << "associating secondary index trans_cache_itridx\n";
  trans_cachedb->associate( dbtxn, trans_cache_iteridx,
                            bdb_transcache_itr_extract, DB_CREATE );
  std::cout << "associating secondary index trans_cache_sliceidx\n";
  trans_cachedb->associate( dbtxn, trans_cache_sliceidx, 
                            bdb_transcache_slice_extract, DB_CREATE );
  std::cout << "associating secondary index trans_cache_sizeidx\n";
  trans_cachedb->associate( dbtxn, trans_cache_sizeidx, 
                            bdb_transcache_size_extract, DB_CREATE );
  std::cout << "associating secondary index trans_cache_hitidx\n";
  trans_cachedb->associate( dbtxn, trans_cache_hitidx, 
                            bdb_transcache_hit_extract, DB_CREATE );
  std::cout << "associating secondary index mi_score_sliceidx\n";
  mi_scoredb->associate( dbtxn, mi_score_sliceidx,
                         bdb_mi_score_slice_extract, DB_CREATE );
  std::cout << "associating secondary index mi_score_xyidx\n";
  mi_scoredb->associate( dbtxn, mi_score_xyidx,
                         bdb_mi_score_xy_extract, DB_CREATE );
  std::cout << "associating secondary index mi_score_scoreidx\n";
  mi_scoredb->associate( dbtxn, mi_score_scoreidx,
                         bdb_mi_score_score_extract, DB_CREATE );
  std::cout << "associating secondary index particle_setup_psosliceidx\n";
  particle_setupdb->associate( dbtxn, particle_setup_psosliceidx, 
                               bdb_particle_setup_psoslice_extract, DB_CREATE );
  std::cout << "associating secondary index particle_state_vxidx\n";
  particle_statedb->associate( dbtxn, particle_state_vxidx,
                               bdb_particle_state_vx_extract, DB_CREATE );
  std::cout << "associating secondary index particle_state_vyidx\n";
  particle_statedb->associate( dbtxn, particle_state_vyidx,
                               bdb_particle_state_vy_extract, DB_CREATE );
  std::cout << "associating secondary index particle_state_speedidx\n";
  particle_statedb->associate( dbtxn, particle_state_speedidx,
                               bdb_particle_state_speed_extract, DB_CREATE );
  std::cout << "associating secondary index particle_state_itridx\n";
  particle_statedb->associate( dbtxn, particle_state_iteridx,
                               bdb_particle_state_iter_extract, DB_CREATE );
  std::cout << "associating secondary index particle_state_pididx\n";
  particle_statedb->associate( dbtxn, particle_state_pididx,
                               bdb_particle_state_pid_extract, DB_CREATE );
  //  dbtxn->commit(0);
#endif
  return 0;
}

double psotracker::getparam( std::string name )
{
  double myval;
  Dbt key, value;
  key.set_data( (char*) name.c_str() );
  key.set_size( name.length() );
  value.set_data( &myval );
  value.set_flags( DB_DBT_USERMEM );
  value.set_ulen( sizeof( double ) );
  int ret;
  ret = pso_paramdb->get( NULL, &key, &value, 0 );
  if( ret != 0 )
    {
      std::cout << "ERROR can not find param\"" << name << "\"\n";
      return 0;
    }
  return myval;
}

int psotracker::setparam( std::string name, double val )
{
  double myval = val;
  Dbt key, value;
  key.set_data( (char*) name.c_str() );
  key.set_size( name.length() );
  value.set_data( &myval );
  value.set_size( sizeof( double ) );
  pso_paramdb->put( mytxn, &key, &value, 0 );
  return 0;
}

int psotracker::set_runconfig( std::vector< psotrack_slice > &psoslicedata,
                               int parperslice,
                               box searchbox )
{
  droptables();
  //  mytxn = bdb.begintransaction();
  mytxn = NULL;
  std::cout << "inserting parperslice data\n";
  setparam( "numpars", parperslice );
  setparam( "searchlx", searchbox.lowx );
  setparam( "searchly", searchbox.lowy );
  setparam( "searchhx", searchbox.highx );
  setparam( "searchhy", searchbox.highy );
  setparam( "searchspace", (searchbox.highx-searchbox.lowx)*
            (searchbox.highy-searchbox.lowy) );
  slice_setup_value ssv;
  int psoslice;
  Dbt key, value;
  key.set_data( &psoslice );
  key.set_size( sizeof(int) );
  value.set_data( &ssv );
  value.set_size( sizeof( slice_setup_value ) );
  std::cout << "inserting psoslice data\n";
  for( unsigned i = 0; i < psoslicedata.size(); i++ )
    {
      psoslice = psoslicedata[i].psoslice;
      ssv.slice = psoslicedata[i].dbslice;
      ssv.dbtx1 = psoslicedata[i].dbtilex1;
      ssv.dbty1 = psoslicedata[i].dbtiley1;
      ssv.dba1 = psoslicedata[i].dbalignedtile1;
      ssv.dbtx2 = psoslicedata[i].dbtilex2;
      ssv.dbty2 = psoslicedata[i].dbtiley2;
      ssv.dba2 = psoslicedata[i].dbalignedtile2;
      ssv.channel = psoslicedata[i].dbchannel;
      slice_setupdb->put( mytxn, &key, &value, 0 );
    }

  int counter;
  counter = 0;
  particle_setup_value psv;
  key.set_data( &counter );
  key.set_size( sizeof(int) );
  value.set_data( &psv );
  value.set_size( sizeof( particle_setup_value ) );
  std::cout << "setting up particle data\n";
  for( unsigned i = 0; i < psoslicedata.size(); i++ )
    {
      for( int j = 0; j < parperslice; j++ )
        {
          psv.psoslice = i;
          psv.parslice = j;
          particle_setupdb->put( mytxn, &key, &value, 0 );
          counter++;
        }
    }
  //  mytxn->commit(DB_TXN_NOSYNC);
  mytxn = NULL;
  std::cout << "done with setup\n";
  myslices = (int) psoslicedata.size();
  myparperslice = parperslice;
  return 0;
}

int psotracker::save_iterationdata( int iteration, std::vector< psoparticle* > &pvec,
                                    std::vector< psotrack_cluster > &kmeans, 
                                    std::vector< psotrack_cachestat > &cache,
                                    std::vector< transformation > *newtrans )
{
  timer mytimer;
  mytimer.reset();
  mytimer.start();
  //  mytxn = bdb.begintransaction();
  mytxn = NULL;
  std::cout << "updating mutual information table\n";
  // save new transformations
  for( int i = 0; i < myslices; i++ )
    for( unsigned j = 0; j < newtrans[i].size(); j++ )
      savemutualinformation( newtrans[i][j], i );
      
  // save the particles
  std::cout << "updating particles\n";
  for( unsigned i = 0; i < pvec.size(); i++ )
    save_pso_particles( iteration, pvec[i]);
  std::cout << "updating kmeans\n";
  for( unsigned i = 0; i < kmeans.size(); i++ )
    save_kmeans( iteration, i, kmeans[i] );
  std::cout << "updating cache\n";
  for( unsigned i = 0; i < cache.size(); i++ )
    save_cache( iteration, cache[i] );
  std::cout << "iteration data saved\n";
  //  mytxn->commit(DB_TXN_NOSYNC);
  mytimer.stop();
  std::cout << "updated " << pvec.size() << " particles in "
            << mytimer.getseconds() << " s\n";
  mytxn = NULL;
  return 0;
}

int psotracker::save_cache( int iteration,  psotrack_cachestat &cs )
{
  trans_cache_key tck;
  trans_cache_value tcv;
  tck.iteration = iteration;
  tck.psoslice = cs.psoslice;
  tcv.size = cs.size;
  tcv.hit = cs.hit;

#ifndef __PSO_TRACK_DEBUG__
  if( !isfinite( tcv.hit ) )
    {
      std::cout << "ERROR pso cache stat not normal\n"
                << tck.iteration << " " << tck.psoslice 
                << " : " << tcv.size << " " << tcv.hit << "\n";
      exit(0);
    }
#endif

  Dbt key, value;
  key.set_data( &tck );
  key.set_size( sizeof( trans_cache_key ) );
  value.set_data( &tcv );
  value.set_size( sizeof( trans_cache_value ) );

  trans_cachedb->put( mytxn, &key, &value, 0 );
  return 0;
}

int psotracker::save_kmeans( int iteration, int cid, psotrack_cluster &c )
{
  cluster_table_key ctk;
  cluster_table_value ctv;

  ctk.iteration = iteration;
  ctk.cid = cid;
  ctv.x = c.x;
  ctv.y = c.y;
#ifndef __PSO_TRACK_DEBUG__
  if( !isfinite( ctv.x ) || !isfinite( ctv.y ) )
    {
      std::cout << "ERROR pso kmeans not normal\n";
      std::cout << ctk.iteration << " " << ctk.cid 
                << " : " << ctv.x << " " << ctv.y << "\n";
      exit(0);
    }
#endif
  Dbt key, value;
  key.set_data( &ctk );
  key.set_size( sizeof( cluster_table_key ) );
  value.set_data( &ctv );
  value.set_size( sizeof( cluster_table_value ) );

  cluster_positiondb->put( mytxn, &key, &value, 0 );
  return 0;
}

int psotracker::save_pso_particles( int iteration, psoparticle *p )
{
  particle_state_value psv;
  particle_state_key psk;
  psk.pid = p->imparnum + (p->imnum*myparperslice);
  psk.iteration = iteration;
  psv.x = p->x;
  psv.y = p->y;
  psv.vx = p->vx;
  psv.vy = p->vy;
  psv.speed = pow(p->vx,2)+pow(p->vy,2);
#ifndef __PSO_TRACK_DEBUG__
  if( !isfinite( psv.x ) || !isfinite( psv.y ) ||
      !isfinite( psv.vx ) || !isfinite( psv.vy ) ||
      !isfinite( psv.speed ) || !isfinite( psk.pid ) )
    {
      std::cout << "ERROR pso particle state value not normal\n";
      std::cout << psk.pid << " " << p->imparnum << ":" << p->imnum << "\n";
      std::cout << psv.x << " " << psv.y << " " << psv.vx << " " 
                << psv.vy << " " << psv.speed << "\n";
      exit(0);
    }
#endif
  Dbt key, value;
  key.set_data( &psk );
  key.set_size( sizeof(particle_state_key) );
  value.set_data( &psv );
  value.set_size( sizeof( particle_state_value) );
  particle_statedb->put( mytxn, &key, &value, 0 );
  return 0;
}

int psotracker::savemutualinformation( transformation &t, int slice )
{
  Dbt key, value;
  mi_score_key msk;
  double msv;
  value.set_data( &msv );
  value.set_size( sizeof(double) );
  key.set_data( &msk );
  key.set_size( sizeof( mi_score_key ) );
  msv = t.score;
  msk.slice = slice;
  msk.x = (int) t.transx;
  msk.y = (int) t.transy;
#ifndef __PSO_TRACK_DEBUG__
  if( !isfinite( msv ) )
    {
      std::cout << "ERROR pso mutual information not normal\n";
      std::cout << msk.slice << " " << msk.x << " " << msk.y 
                << ":" << msv << "\n";
      exit(0);
    }
#endif
  mi_scoredb->put( mytxn, &key, &value, 0 );
  return 0;
}

int psotracker::droptables()
{
  unsigned int counter;
  std::cout << "removing from database ";
  pso_paramdb->truncate( NULL, &counter, 0 );
  std::cout << counter << " ";
  slice_setupdb->truncate( NULL, &counter, 0 );
  std::cout << counter << " ";
  trans_cachedb->truncate( NULL, &counter, 0 );
  std::cout << counter << " ";
  mi_scoredb->truncate( NULL, &counter, 0 );
  std::cout << counter << " ";
  particle_setupdb->truncate( NULL, &counter, 0 );
  std::cout << counter << " ";
  particle_statedb->truncate( NULL, &counter, 0 );
  std::cout << counter << " ";
  iteration_statdb->truncate( NULL, &counter, 0 );
  std::cout << counter << " ";
  cluster_positiondb->truncate( NULL, &counter, 0 );
  std::cout << counter << "\n";
  return 0;
}

double psotracker::searched()
{
  Dbt key, pkey, value;
  mi_score_key msk;
  double score;
  value.set_data( &score );
  value.set_ulen( sizeof( double ) );
  value.set_flags( DB_DBT_USERMEM );
  int xy[2];
  key.set_data( xy );
  key.set_ulen( 2*sizeof(int) );
  key.set_flags( DB_DBT_USERMEM );

  pkey.set_data( &msk );
  pkey.set_ulen( sizeof( mi_score_key ) );
  pkey.set_flags( DB_DBT_USERMEM );
  int uniquexy;
  uniquexy = 0;
  Dbc *cur;
  int ret;
  mi_score_xyidx->cursor( NULL, &cur, 0 );
  ret = cur->pget( &key, &pkey, &value, DB_FIRST );
  while( ret == 0 )
    {
      uniquexy++;
      ret = cur->pget( &key, &pkey, &value, DB_NEXT_NODUP );
    }
  cur->close();
  double tosp;
  tosp = getparam( "searchspace" );
  tosp = ( (double)uniquexy )/ tosp;
  return tosp;
}

int psotracker::printavgmiscore()
{
  Dbt pkey, key, value;
  mi_score_key msk;
  double score;
  double lowx, lowy, highx, highy;
  lowx = getparam( "searchlx" );
  highx = getparam( "searchhx" );
  lowy = getparam( "searchly" );
  highy = getparam( "searchhy" );

  int imszx, imszy;
  imszx = (int)( highx - lowx );
  imszy = (int)( highy - lowy );

  image2d ss, stdim;
  ss.create( imszx, imszy );
  ss.clear();
  stdim.create( imszx, imszy );
  stdim.clear();

  value.set_data( &score );
  value.set_ulen( sizeof( double ) );
  value.set_flags( DB_DBT_USERMEM );
  int xy[2];
  key.set_data( xy );
  key.set_size( 2*sizeof(int) );
  key.set_ulen( 2*sizeof(int) );
  key.set_flags( DB_DBT_USERMEM );

  pkey.set_data( &msk );
  pkey.set_ulen( sizeof( mi_score_key ) );
  pkey.set_flags( DB_DBT_USERMEM );
  Dbc *cur;
  mi_score_xyidx->cursor( NULL, &cur, 0 );
  int ret;
  std::cout << "generating avg and std images...\n";
  double minavg;
  minavg = 10000;
  for( int i = 0; i < imszx; i++ )
    {
      std::cout << "." << std::flush;
      for( int j = 0; j < imszy; j++ )
        {
          xy[0] = (int)( i+lowx );
          xy[1] = (int)( j+lowy );
          ret = cur->pget( &key, &pkey, &value, DB_SET );
          std::vector< double > vlist;
          vlist.clear();
          while( ret == 0 )
            {
              vlist.push_back( score );
              ret = cur->pget( &key, &pkey, &value, DB_NEXT_DUP );
            }
          double avg, std;
          if( avg > 0.0001 && minavg > avg )
            minavg = avg;
          avgstd( vlist, avg, std );
          ss.setpixel(i, j, avg );
          stdim.setpixel( i, j, std );
        }
    }
  std::cout << "\n";
  std::cout << "min score: " << minavg << "\n";
  for( int i = 0; i < imszx; i++ )
    for( int j = 0; j < imszy; j++ )
      if( ss.getpixel(i,j) < minavg/2.0 )
        ss.setpixel(i,j, minavg );
  cur->close();
#ifndef _NOMAGICK_
  Magick::Image mim;
  ss.getimage_normalized( mim );
  overlaygrid( mim );
  mim.write( "avg.png" );
  stdim.getimage_normalized( mim );
  overlaygrid( mim );
  mim.write( "std.png" );
#endif
  return 0;
}

#ifndef _NOMAGICK_
int psotracker::overlaygaussians( Magick::Image &mim,
                                  std::vector< gaussianmodel* > gms )
{
  double lowx, lowy, highx, highy;
  lowx = getparam( "searchlx" );
  highx = getparam( "searchhx" );
  lowy = getparam( "searchly" );
  highy = getparam( "searchhy" );
  Magick::Geometry geo;
  geo = mim.size();
  Magick::ColorRGB rgb("green");
  rgb.alpha( 1 ); 
  mim.strokeColor( "blue" );
  mim.fillColor( rgb );
  for( unsigned i = 0; i < gms.size(); i++ )
    {
      matrix mean, covar;
      gms[i]->getmean( mean );
      gms[i]->getvariance( covar );
      int x, y;
      x = mean.get(0, 0) - lowx;
      y = mean.get(1, 0) - lowy;
      std::cout << "overlaying cluster at " << x << " " << y
                << " "<< mean.get(0,0) << " " << mean.get(0,1 ) << "\n";
      mim.draw( Magick::DrawableCircle(y, x, y+5, x) );
    }
  return 0;
}
#endif

int psotracker::avgstd( std::vector< double > &vlist,
                        double &avg, double &std )
{
  avg = std = 0;
  if( vlist.size() == 0 )
    return 0;
  for( unsigned i = 0; i < vlist.size(); i++ )
    avg = avg + vlist[i];
  avg = avg / ( (double)(vlist.size()) );
  for( unsigned i = 0; i < vlist.size(); i++ )
    std = std + pow( vlist[i]-avg, 2 );
  if( vlist.size () > 1 )
    std = std / ( (double)(vlist.size()-1) );
    
  return 0;
}

int psotracker::printsearched( )
{
  Dbt key, value;
  mi_score_key msk;
  double score;
  double lowx, lowy, highx, highy;
  lowx = getparam( "searchlx" );
  highx = getparam( "searchhx" );
  lowy = getparam( "searchly" );
  highy = getparam( "searchhy" );

  int imszx, imszy;
  imszx = (int)( highx - lowx );
  imszy = (int)( highy - lowy );

  image2d ss;
  ss.create( imszx+1, imszy+1 );
  ss.clear();

  value.set_data( &score );
  value.set_ulen( sizeof( double ) );
  value.set_flags( DB_DBT_USERMEM );

  key.set_data( &msk );
  key.set_ulen( sizeof( mi_score_key ) );
  key.set_flags( DB_DBT_USERMEM );

  Dbc *cur;
  int ret;
  std::cout << "getting overlapping search regions\n";
  mi_scoredb->cursor( NULL, &cur, 0 );
  ret = cur->get( &key, &value, DB_FIRST );
  while( ret == 0 )
    {
      double np;
      int x, y;
      x = (int)( msk.x - lowx );
      y = (int)( msk.y - lowy );
      np = ss.getpixel( x, y ) + 1;
      ss.setpixel( x, y, np );
      ret = cur->get( &key, &value, DB_NEXT );
    }
  cur->close();
#ifndef _NOMAGICK_
  Magick::Image mim;
  ss.getimage_normalized( mim );
  overlaygrid( mim );
  mim.write( "searchregion.png" );
#endif
  return 0;
}
#ifndef _NOMAGICK_
int psotracker::overlaygrid( Magick::Image &mim )
{
  double lowx, lowy, highx, highy;
  lowx = getparam( "searchlx" );
  highx = getparam( "searchhx" );
  lowy = getparam( "searchly" );
  highy = getparam( "searchhy" );
  mim.colorSpace( Magick::RGBColorspace );
  mim.type( Magick::TrueColorType );
  mim.strokeColor( "red" );
  mim.strokeWidth( 1 );
  Magick::Geometry geo;
  geo = mim.size();
  char buffer[512];
  std::cout << lowx << " " << lowy << " : " << highx << " " << highy << "\n";
  for( double i = 0; i < 3; i++ )
    {
      double z;
      z = lowx + (highx-lowx)*( (i+1)*0.25 );
      sprintf( buffer, "%03.1f", z );
      int h;
      h = geo.height() * ( (i+1)*0.25 );
      mim.draw( Magick::DrawableLine(0, h, geo.width(), h ) );
      mim.draw( Magick::DrawableText(0, h-5, buffer ) );
    }
  mim.strokeColor( "green" );
  for( double i = 0; i < 3; i++ )
    {
      double z;
      z = lowy + (highy-lowy)*( (i+1)*0.25 );
      sprintf( buffer, "%03.1f", z );
      int w;
      w = geo.width() * ((i+1)*0.25);
      mim.draw( Magick::DrawableLine(w, 0, w, geo.height() ) );
      mim.draw( Magick::DrawableText(w, 20, buffer ) );
    }
  return 0;
}
#endif

int psotracker::printpoints( int iteration )
{
  double lowx, lowy, highx, highy;

#ifndef _NOMAGICK_
  Magick::Geometry geo;
#endif

  lowx = getparam( "searchlx" );
  highx = getparam( "searchhx" );
  lowy = getparam( "searchly" );
  highy = getparam( "searchhy" );

#ifndef _NOMAGICK_
  geo.width( highy - lowy + 1 );
  geo.height( highx - lowx + 1 );
  Magick::Image mim( geo, Magick::Color("black") );
#endif

  std::list< particle_state_value > plist;
  filterparticles( iteration, -1, plist );
  filteroutliers( 2, plist );
#ifndef _NOMAGICK_  
  std::list<Magick::Drawable> drawList;
#endif
  std::vector< gaussianmodel* > gms;
  clusterparticles( plist, gms );

  // generate a list of colors for the gaussian models
#ifndef _NOMAGICK_
  std::vector< Magick::ColorRGB > colors;
  generatecolors( colors, gms.size() );

  drawList.push_back(Magick::DrawableStrokeWidth( 2 ));
#endif

  std::vector< int > csizes;
  csizes.resize( gms.size(), 0 );

  for( std::list< particle_state_value >::iterator it = plist.begin();
       it != plist.end(); it++ )
    {
      int x, y;
      int cid;
      cid = closestgaussian( it->x, it->y, gms );
      csizes[cid] = csizes[cid] + 1;
#ifndef _NOMAGICK_
      drawList.push_back(Magick::DrawableStrokeColor( colors[cid] ));
      drawList.push_back(Magick::DrawableFillColor( colors[cid] ));
#endif
      x = (int)( it->x - lowx );
      y = (int)( it->y - lowy );
#ifndef _NOMAGICK_
      if( x < 0 || (int) geo.height() <= x  || 0 > y || (int) geo.width() <= y )
        {
          std::cout << "ERROR x,y not correct " << x << " " << y << " : " 
                    << geo.height() << " " << geo.width() << "\n";
        }

      drawList.push_back( Magick::DrawablePoint(y, x) );
#endif
    }

#ifndef _NOMAGICK_
  mim.draw(drawList);
  overlaygaussians( mim, gms );
#endif

  std::cout << "cluster sizes\n";
  for( unsigned i = 0; i < gms.size(); i++ )
    {
      std::cout << "------------------------------\n";
      std::cout << *(gms[i]) << "\n";
      std::cout << "size: " << csizes[i] << "\n";
      std::cout << "------------------------------\n";
      delete gms[i];
    }
#ifndef _NOMAGICK_
  mim.write( "points.png" );
#endif
  return 0;
}

int psotracker::filteroutliers( int radius,
                                std::list< particle_state_value > &plist )
{
  std::map< std::pair< int, int>, int > pmap;
  pmap.clear();
  for( std::list<particle_state_value>::iterator it = plist.begin();
       it != plist.end(); it++ )
    {
      std::pair< int, int > key;
      key.first = (int)round( it->x );
      key.second = (int)round( it->y );
      if( pmap.count(key) == 0 )
        pmap[key] = 1;
      else
        pmap[key] = pmap[key]+1;
    }
  std::vector< std::list<particle_state_value>::iterator > todel;
  for( std::list<particle_state_value>::iterator it = plist.begin();
       it != plist.end(); it++ )
    {
      std::pair< int, int > key;
      key.first = (int)round( it->x );
      key.second = (int)round( it->y );
      if( pmap[key] < 3 )
        {
          int c;
          c = 0;
          for( int i = (int)round( it->x ) - radius;
               i < (int)round( it->x )+1+radius && c < 3; i++ )
            {
              for( int j = (int)round( it->y ) - radius;
                   j < (int)round( it->y )+1+radius && c < 3; j++ )
                {
                  if( i != (int)round( it->x ) || j != (int)round( it->y ) )
                    {
                      key.first = i;
                      key.second = j;
                      if( pmap.count( key ) > 0 )
                        c++;
                    }
                }
            }
          if( c < 2 )
            todel.push_back( it );
        }
        
    }
  std::cout << "deleting: " << todel.size() << "\n";
  for( unsigned i = 0; i < todel.size(); i++ )
    plist.erase( todel[i] );
  std::cout << "filtered list to: " << plist.size() << "\n";
  return 0;
}

int psotracker::closestgaussian( double x, double y,
                                 std::vector< gaussianmodel* > &gsm )
{
  int res;
  res = 0;
  long double sc;
  matrix m;
  m.init( 2, 1 );
  m.set(x, 0, 0 );
  m.set(y, 1, 0 );
  sc = gsm[0]->prob( m );
  for( unsigned i = 1; i < gsm.size(); i++ )
    {
      long double tsc;
      tsc = gsm[i]->prob( m );
      if( tsc > sc )
        {
          res = (int)i;
          sc = tsc;
        }
    }
  return res;
}

int psotracker::filterparticles( int iteration, double speed,
                                 std::list< particle_state_value > &plist )
{
  particle_state_value psv;
  particle_state_key psk;
  Dbt pkey, key, value;
  Dbc *cur; 
  int sid, ret;

  double lowx, lowy, highx, highy;
  lowx = getparam( "searchlx" );
  highx = getparam( "searchhx" );
  lowy = getparam( "searchly" );
  highy = getparam( "searchhy" ); 

  sid = iteration;
  pkey.set_data( &psk );
  pkey.set_ulen( sizeof( particle_state_key ) );
  pkey.set_flags( DB_DBT_USERMEM );
  key.set_data( &sid );
  key.set_size( sizeof( int ) );
  value.set_data( &psv );
  value.set_ulen( sizeof( particle_state_value ) );
  value.set_flags( DB_DBT_USERMEM );
  particle_state_iteridx->cursor( NULL, &cur, 0 );

  ret = cur->pget( &key, &pkey, &value, DB_SET );
  double avg;
  int cc;
  avg = 0;
  cc = 0;
  while( ret == 0 )
    {
      if( (psv.speed < speed || speed < 0) && psv.speed < 1000 &&
          ( (lowx + 1) <= psv.x && (lowy + 1) <= psv.y &&
            (highx -1) >= psv.x && (highy -1) >= psv.y ) )
        {
          plist.push_back( psv );
          avg = avg + psv.speed;
          cc++;
        }
      ret = cur->pget( &key, &pkey, &value, DB_NEXT_DUP );

    }
  cur->close();

  if( speed < 0 )
    {
      avg = avg / ( (double)cc );
      std::cout << "average speed: " << avg << "\n";
      std::list< particle_state_value >::iterator it, it1;
      it = plist.begin();
      while( it != plist.end() )
        {
          if( it->speed > avg )
            {
              it1 = it;
              it++;
              plist.erase( it1 );
            }
          else
            it++;
        }
      std::cout << "filtered to: " << plist.size() << "\n";
    }
  return 0;
}

#ifndef _NOMAGICK_
int psotracker::generatecolors( std::vector< Magick::ColorRGB > &colors,
                                int ncolors )
{
  int ncomb;
  int cstep;
  cstep = 2;
  ncomb = (int)pow(cstep, 3)-cstep;
  while(ncomb < ncolors )
    {
      cstep++;
      ncomb = (int)pow(cstep, 3)-cstep;
    }
  double stepsize;
  stepsize = 1.0/( (double)(cstep-1) );
  std::cout << "colors\n";
  for( double i = 0; i < 1.0+(stepsize/2.0); i = i + stepsize )
    {
      for( double j = 0; j < 1.0+(stepsize/2.0); j = j + stepsize )
        {
          for( double k = 0; k < 1.0+(stepsize/2.0); k = k + stepsize )
            {
              if( i != j || j != k )
                {
                  std::cout << "(" << i << " " << j << " " << k << ")\n";
                  Magick::ColorRGB rgb(1,1,1);
                  rgb.red(i); rgb.green(j); rgb.blue(k);
                  colors.push_back( rgb );
                }
            }
        }
    }
  std::cout << "generated " << colors.size() << " colors\n";
  return 0;
}
#endif

int psotracker::clusterparticles( std::list< particle_state_value > &lpsv,
                                  std::vector< gaussianmodel* > &gms )
{
  double **positions;
  positions = new double*[ lpsv.size() ];
  int c;
  c = 0;
  for( std::list<particle_state_value>::iterator it = lpsv.begin();
       it != lpsv.end(); it++ )
    {
      positions[c] = new double[2];
      positions[c][0] = it->x;
      positions[c][1] = it->y;
      c++;
    }

  emcluster em;
  em.loaddata( positions, (int)lpsv.size(), 2, 1 );
  for( int i = 0; i < (int)lpsv.size(); i++ )
    {
      if( positions[i] == NULL )
        {
          std::cout << "?!?!?!?\n";
          exit(0);
        }
      delete [] positions[i];
    }
  delete [] positions;
  int nclusters;
  nclusters = 3;
  //nclusters = em.bicrun( 2, 7 );
  std::cout << "located " << nclusters << "\n";
  em.run( nclusters );
  gms = em.getmodels();
  em.printmodels();
  return 0;
}

int psotracker::printsliceparticles( double &lowx, double &lowy,
                                     double &highx, double &highy)
{
  lowx = getparam( "searchlx" );
  highx = getparam( "searchhx" );
  lowy = getparam( "searchly" );
  highy = getparam( "searchhy" );
  return 0;
}

int psotracker::printsliceparticles( int iteration )
{
  particle_state_value psv;
  particle_state_key psk;
  particle_setup_value pstv;
  Dbc *cur;
  Dbt key, value, skey;
  Dbt pskey, psvalue;
  Dbt value1;
  int sid;
  int iter;
  double lowx, lowy, highx, highy;

  lowx = getparam( "searchlx" );
  highx = getparam( "searchhx" );
  lowy = getparam( "searchly" );
  highy = getparam( "searchhy" );

  int sx, sy;
  sx = (int)( highx - lowx + 1 );
  sy = (int)( highy - lowy + 1 );

  iter = iteration;
  skey.set_data( &iter );
  skey.set_size( sizeof(int) );
  key.set_data( &psk );
  key.set_flags( DB_DBT_USERMEM );
  key.set_ulen( sizeof( particle_state_key ) );
  value.set_data( &psv );
  value.set_flags( DB_DBT_USERMEM );
  value.set_ulen( sizeof( particle_state_value ) );
  pskey.set_data( &(psk.pid) );
  pskey.set_flags( DB_DBT_USERMEM );
  pskey.set_size( sizeof(int) );
  pskey.set_ulen( sizeof(int) );
  psvalue.set_data( &pstv );
  psvalue.set_flags( DB_DBT_USERMEM );
  psvalue.set_ulen( sizeof(particle_setup_value) );
  int ret;
  particle_state_iteridx->cursor( NULL, &cur, 0 );
  ret = cur->pget( &skey, &key, &value, DB_SET );
  std::map< int, image2d* > ims;
  while( ret == 0 )
    {
      particle_setupdb->get( NULL, &pskey, &psvalue, 0 );    
      if( ims.count( pstv.psoslice ) == 0 )
        {
          std::cout << "found pso slice: " << pstv.psoslice << "\n";
          ims[pstv.psoslice] = new image2d;
          ims[pstv.psoslice]->create( sx, sy );
          ims[pstv.psoslice]->clear();
        }
      int x, y;
      x = (int)round( psv.x - lowx );
      y = (int)round( psv.y - lowy );
      ims[pstv.psoslice]->setpixel( x, y, 1 );
      ret = cur->pget( &skey, &key, &value, DB_NEXT_DUP );
    }
  cur->close();
  slice_setup_value ssv;
  value.set_data( &ssv );
  value.set_ulen( sizeof( slice_setup_value ) );
  key.set_data( &sid );
  key.set_size( sizeof(int) );
  for( std::map< int, image2d* >::iterator it = ims.begin(); 
       it != ims.end(); it++ )
    {
      sid = it->first;
      slice_setupdb->get( NULL, &key, &value, 0 );
      char buffer[512];
      std::cout << sid << " : " << ssv.slice << " " 
                << ssv.dba1 << " " << ssv.dba2 << " " << ssv.channel << "\n";
      sprintf(buffer, "%i-%i-%i-%i.png", ssv.slice, ssv.dba1, ssv.dba2, ssv.channel );
#ifndef _NOMAGICK_
      Magick::Image mim;
      it->second->getimage(mim);
      mim.write( buffer );
#endif
    }
  for( std::map< int, image2d* >::iterator it = ims.begin(); 
       it != ims.end(); it++ )
    delete it->second;
  return 0;
}

int psotracker::getparticles( int iteration, int psoslice, 
                              std::vector< psoparticle > &particles )
{
  // find the particle ids for a given slice
  Dbc *cur;
  particle_setup_psosliceidx->cursor( NULL, &cur, 0 );
  Dbt key, value, pkey;
  Dbt key1, value1;
  int slice;
  int pid;
  particle_setup_value psv;
  particle_state_key psk;
  particle_state_value pstatev;
  slice = psoslice;
  key.set_data( &slice );
  key.set_size( sizeof(int) );
  value.set_data( &psv );
  value.set_ulen( sizeof( particle_setup_value ) );
  value.set_flags( DB_DBT_USERMEM );
  pkey.set_data( &pid );
  pkey.set_ulen( sizeof(int) );
  pkey.set_flags( DB_DBT_USERMEM );
  int ret;

  key1.set_data( &psk );
  key1.set_size( sizeof( particle_state_key ) );
  value1.set_data( &pstatev );
  value1.set_ulen( sizeof( particle_state_value ) );
  value1.set_flags( DB_DBT_USERMEM );

  ret = cur->pget( &key, &pkey, &value, DB_SET );
  while( ret == 0 )
    {
      //get the particle
      psk.iteration = iteration;
      psk.pid = pid;
      particle_statedb->get( NULL, &key1, &value1, 0 );
      psoparticle ppt;
      ppt.x = pstatev.x;
      ppt.y = pstatev.y;
      ppt.vx = pstatev.vx;
      ppt.vy = pstatev.vy;
      ppt.pid = pid;
      particles.push_back( ppt );
      ret = cur->pget( &key, &pkey, &value, DB_NEXT_DUP );
    }

  cur->close();
  return 0;
}

int psotracker::getparticles( int iteration,
                              std::vector< psoparticle > &particles )
{
  particle_state_value psv;
  particle_state_key psk;
  Dbt pkey, key, value;
  Dbc *cur; 
  int sid, ret;

  sid = iteration;
  pkey.set_data( &psk );
  pkey.set_ulen( sizeof( particle_state_key ) );
  pkey.set_flags( DB_DBT_USERMEM );
  key.set_data( &sid );
  key.set_size( sizeof( int ) );
  value.set_data( &psv );
  value.set_ulen( sizeof( particle_state_value ) );
  value.set_flags( DB_DBT_USERMEM );
  particle_state_iteridx->cursor( NULL, &cur, 0 );

  ret = cur->pget( &key, &pkey, &value, DB_SET );
  while( ret == 0 )
    {
      psoparticle ppt;
      ppt.x = psv.x;
      ppt.y = psv.y;
      ppt.vx = psv.vx;
      ppt.vy = psv.vy;
      ppt.pid = psk.pid;
      particles.push_back( ppt );
      ret = cur->pget( &key, &pkey, &value, DB_NEXT_DUP );
    }

  cur->close();
  return 0;
}

