/********************************************************************/
/**                                                                **/
/**   debug.c                                                      **/
/**                                                                **/
/**   A universal generator for multivariate T-concave             **/
/**   distributions                                                **/
/**                                                                **/
/**   print debugging information                                  **/
/**                                                                **/
/**   author: by Josef Leydold                                     **/
/**   last modified: Mon Mar  9 13:38:00 MET 1998                  **/
/**                                                                **/
/********************************************************************/              

/* the header file */
#include "tdrmv.h"

/*------------------------------------------------------------------*/

void fatal(char *msg, char *filename, int line)
/*
  fatal error, exit 
*/
{
  fprintf(stderr,"fatal error in %s, line %d: %s\n",filename,line,msg);
#if DEBUG > 0
  fprintf(LOG,"fatal error in %s, line %d: %s\n",filename,line,msg);
#endif
  exit(1);
} /* end if fatal() */

/*------------------------------------------------------------------*/

void warning(char *msg, char *filename, int line)
/*
  warning message 
*/
{
  fprintf(stderr,"warning in %s, line %d: %s\n",filename,line,msg);
#if DEBUG > 0
  fprintf(LOG,"warning in %s, line %d: %s\n",filename,line,msg);
#endif
} /* end if warning() */

/**----------------------------------------------------------------**/

#if DEBUG > 0

/** local definitions                                              **/

/* timer                                                            */
/* estimate time                                                    */
#define TICKS          1000000.0
#define F_TO_TV(f, tv) (((tv).tv_sec = (long) floor(f)), ((tv).tv_usec = (long) (((f) - floor(f)) * TICKS)))
#define TV_TO_F(tv)    ((tv).tv_sec + (tv).tv_usec / TICKS)

struct timeval timer;        /* structure needed for measuring time */
double time_start;           /* time of start of program            */
double time_start_run;       /* time of start of sampling           */

/*------------------------------------------------------------------*/

/* internal functions */
static void db_print_vertices(C_ROOT *);             /* print list of vertices */
static void db_print_cones(C_ROOT *);                /* print list of cones */
static void db_print_cone_params(CONE *,int);        /* print parameters of cone */
static void db_print_edges(C_ROOT *);                /* print edges for each cone */
static void db_print_guide(C_ROOT *);                /* print guide table statistics */
#if DEBUG > 1000
static void db_print_cones_stat_volume(C_ROOT *);    /* print statistics of volumes of cones */
#endif

void db_print_vector(char *, double *, int, char *); /* print double vector */
void db_print_string(char *);                        /* print character string into log file */

static double vec_dist(double *, double *, int);     /* distance of two vectors */
static int *int_zero_alloc(int n);                   /* allocate pointer to integer array and initialize */

/* external functions */
void db_tdrg(C_ROOT *cr);                            /* print parameters of TDRG gamma generator into log files */


/*-----------------------------------------------------------------*/

/*******************************************************************
  
  main

******************************************************************/

/*-----------------------------------------------------------------*/

void db_initialize()
/* 
   initialize debuging variables 
*/
/* (all debuging variables are global!) */
{

  /* open log file */
  if( (LOG = fopen(LOG_file_name,"w")) == NULL )
      FATAL("cannot open log file");

  /* write header into log file */
  fprintf(LOG,"TDRMV random number generator\n\n");
  fprintf(LOG,"==================================================\n\n");
  
  /* print parameters */
  fprintf(LOG,"dimension = %d\n",N);
  fprintf(LOG,"\n==================================================\n\n");

  /* start timer */
  gettimeofday(&timer, NULL);
  time_start = TV_TO_F(timer);

} /* end of db_initialize() */

/*-----------------------------------------------------------------*/

void db_hat(C_ROOT *cr)
/* write statistics about hat generation into log file */
{
  double time_stop;

  /* stop timer */
  gettimeofday(&timer, NULL);
  time_stop = TV_TO_F(timer);

#if MODE == 1
  db_print_vector("mode = ",cr->mode,N,"\n");
#else
  fprintf(LOG,"mode = origin\n");
#endif  
  fprintf(LOG,"\n==================================================\n\n");

#if RECTANGLE == 1
  fprintf(LOG,"domain = rectangle\n");
  {
    int i;
    for( i=0; i<N; i++ )
      fprintf(LOG,"\t[ %f , %f ]\n",(cr->hc)->rl[i],(cr->hc)->ru[i]);
  }
#else
  fprintf(LOG,"domain = R^n\n");
#endif  
  fprintf(LOG,"\n==================================================\n\n");

  /* print data */
  fprintf(LOG,"* Setup *\n\n");
  fprintf(LOG,"time = %g sec\n",time_stop - time_start);
  fprintf(LOG,"Volume: hat =     %f\n",cr->Htot);
  if( cr->volume_density > 0. ) {
    fprintf(LOG,"\tdensity = %f\n",cr->volume_density);
    fprintf(LOG,"\tratio =   %f %%\n",( cr->volume_density / cr->Htot )*100);
  }

  /* triangulation steps */
  fprintf(LOG,"minimum triangulation level = %d\n",cr->steps_min);
  fprintf(LOG,"maximum triangulation level = %d\n",cr->steps_max);
  fprintf(LOG,"highest level for finding touching points = %d\n",cr->step_tp);

  /* vertices */
  db_print_vertices(cr);

  /* cones */
  db_print_cones(cr);

  /* edges */
  db_print_edges(cr);

  /* guide table */
  db_print_guide(cr);

  fprintf(LOG,"\n==================================================\n\n");

#if GAMMA == 2
  /* print parameters of TDRG gamma generator into log files */
  if( DEBUG & DB_GAMMA ) db_tdrg(cr);
#endif

  /* start timer for sampling */
  gettimeofday(&timer, NULL);
  time_start_run = TV_TO_F(timer);

} /* end of db_hat */

/*-----------------------------------------------------------------*/

void db_close(C_ROOT *cr)
/*
  write statistics into log file and close log file 
*/
{
  double time_stop, time_run;
  int db_n_rpoint_total;

  /* stop timer */
  gettimeofday(&timer, NULL);
  time_stop = TV_TO_F(timer);
  time_run = time_stop - time_start_run;

  /* print data */
  fprintf(LOG,"\n==================================================\n\n");

  fprintf(LOG,"* Run *\n\n");
  fprintf(LOG,"time = %g sec\n\n",time_run);

  fprintf(LOG,"find cones: use search table\n\n");

  db_n_rpoint_total = cr->db_n_rpoint_accept + cr->db_n_rpoint_reject;
  fprintf(LOG,"points: generated = %d\n",db_n_rpoint_total);
  fprintf(LOG,"\taccepted = %d \t(%f %%)\n",cr->db_n_rpoint_accept,
	  ((double)cr->db_n_rpoint_accept)/((double)db_n_rpoint_total)*100.);
  fprintf(LOG,"\trejected = %d\n",cr->db_n_rpoint_reject);
#if RECTANGLE == 1
  fprintf(LOG,"\t(outside domain = %d)\n",cr->db_n_rpoint_reject_outside);
#endif
  fprintf(LOG,"\nrate of acceptance = %f %%",
 	  ((double)cr->db_n_rpoint_accept)/((double)db_n_rpoint_total)*100.);

  if( cr->volume_density > 0. )
    fprintf(LOG,"  \t(expected = %f %%)\n",
	    ( cr->volume_density / cr->Htot )*100.);
  fprintf(LOG,"\n");
  fprintf(LOG,"rejection constant = %f",
 	  ((double)db_n_rpoint_total)/((double)cr->db_n_rpoint_accept));

  if( cr->volume_density > 0. )
    fprintf(LOG,"  \t(expected = %f)",
	    ( cr->Htot / cr->volume_density ));
  fprintf(LOG,"\n");
  fprintf(LOG,"\ngenerated points per second = %f\n",
	  ((double)db_n_rpoint_total)/time_run);
  fprintf(LOG," accepted points per second = %f\n",
	  ((double)cr->db_n_rpoint_accept)/time_run);

#if GAMMA == 2
  if( DEBUG & DB_GAMMA ) {
    fprintf(LOG,"\n==================================================\n\n");
    fprintf(LOG,"TDRG generator for gamma variates:\n");
    fprintf(LOG,"\tgenerated = %d\n",cr->db_n_tdrg_accept+cr->db_n_tdrg_reject);
    fprintf(LOG,"\taccepted = %d \t(%f %%)\n",cr->db_n_tdrg_accept,
	    ((double)cr->db_n_tdrg_accept)/((double)cr->db_n_tdrg_accept+cr->db_n_tdrg_reject)*100.);
    fprintf(LOG,"\trejected = %d\n",cr->db_n_tdrg_reject);
  }
#endif

  fprintf(LOG,"\n==================================================\n\n");
  fprintf(LOG,"* Total *\n\n");
  fprintf(LOG,"time = %g sec\n",time_stop - time_start);

  fprintf(LOG,"\n==================================================\n\n");

  /* close LOG file */
  fclose( LOG );

} /* end of db_close() */

/*-----------------------------------------------------------------*/

/*******************************************************************
  
  internal functions

******************************************************************/

/*-----------------------------------------------------------------*/

static void db_print_vertices(C_ROOT *cr)
/*
  print list of edges 
*/
{
  VERTEX *vt;

  fprintf(LOG,"\n--------------------------------------------------\n\n");
  fprintf(LOG,"number of vertices = %d\n",cr->max_v);

  if( DEBUG & DB_VERTICES ) {
    fprintf(LOG,"\nList of vertices:\n\n");

    for( vt=cr->v; vt!=NULL; vt=vt->next ) {
      fprintf(LOG,"[%4d] = ",vt->index);
      db_print_vector("",vt->coord,N,"");
      fprintf(LOG,"\tnorm = %g\n",vt->norm);
      if( vt==cr->last_v ) break;  /* the last entry in the list */
    }
  }

} /* end of db_print_vertices() */

/*-----------------------------------------------------------------*/

static void db_print_cones(C_ROOT *cr)
/*
  print list of cones
*/
{
  int k;       /* aux counter */
  int *levels;
  CONE *c;

  fprintf(LOG,"\n--------------------------------------------------\n\n");
  fprintf(LOG,"number of cones = %d\n",cr->max_c);

  /* parameter statistics */
#if DEBUG > 1000
  if( DEBUG & DB_VOLUME )
    db_print_cones_stat_volume(cr);
#endif

  /* print statistics about triangulation levels */
  fprintf(LOG,"\n\tlevel \t cones\n");
  if( cr->steps_max >= 0 ) {
    levels = int_zero_alloc(cr->steps_max+1);  
    for( c=cr->c; c!=NULL; c=c->next ) {
      ++levels[c->level];
      if( c == cr->last_c ) break;
    }
    for( k=cr->steps_min; k<=cr->steps_max; k++ )
      fprintf(LOG,"\t%5d \t %5d  (%5.2f%%)\n",k,levels[k], 100. * ((double) levels[k]) / cr->max_c);
    free( levels );
  }
  else
    fprintf(LOG,"\n\t?????????????\n");

  /* list of cones */
  if( DEBUG & DB_CONES ) {
    fprintf(LOG,"\nList of cones:\n\n");
    for( c=cr->c,k=0; c!=NULL; c=c->next,k++ ) {
      db_print_cone_params(c,k);
      if( c == cr->last_c ) break;
    }
  }

  /* store number of cone */
  for( c=cr->c,k=0; c!=NULL; c=c->next,k++ ) {
    c->index = k;
    if( c == cr->last_c ) break;
  }

} /* end of db_print_cones() */

/*-----------------------------------------------------------------*/

static void db_print_cone_params(CONE *c,int nr)
/*
  print parameters of cone
*/
{
  int i;

  fprintf(LOG,"%4d (%2d):",nr,c->level);
  for( i=0; i<c->level; i++ )
    fprintf(LOG," ");

  /* list of vertices */
  fprintf(LOG,"( %d",(c->v[0])->index);
  for( i=1; i<N; i++ )
    fprintf(LOG,", %d",(c->v[i])->index);
  fprintf(LOG," )\n");

  if( DEBUG & DB_CPARAMS ) {
#if DEBUG > 1000
    if( DEBUG & DB_VOLUME )
      fprintf(LOG,"\tVolume = %g\n",c->Vi);
#endif
    fprintf(LOG,"\tdet/(n-1)! = %g\n",c->detf);
    db_print_vector("\tcenter = ",c->center,N,"\n");
    fprintf(LOG,"\ttp = %g\n",c->tp);
#if DEBUG > 1000
    fprintf(LOG,"\tf(tp) = %g\n",c->fp);
    if( DEBUG & DB_G )
      db_print_vector("\tg = ",c->g,N,"\n");
#endif
    db_print_vector("\tgv = ",c->gv,N,"\n");
    fprintf(LOG,"\talpha = %g\n",c->alpha);
    fprintf(LOG,"\tbeta = %g\n",c->beta);
    fprintf(LOG,"\tai = %g\n",c->ai);
    fprintf(LOG,"\tHi = %g\t(%g)\n",c->Hi,c->Hsum);
#if RECTANGLE == 1
    fprintf(LOG,"\theight = %g\n",c->height);
#endif
  }
} /* end of db_print_cone_params() */

/*-----------------------------------------------------------------*/

static void db_print_edges(C_ROOT *cr)
/*
  print edges for each cone
*/
{
  int i,j,k,n_e;
  double length;
  CONE *c;

  if( DEBUG & DB_EDGES ) {
    fprintf(LOG,"\n--------------------------------------------------\n\n");
    fprintf(LOG,"List of edges for each cone:\n\n");
    fprintf(LOG,"[ cone ] (  edge   ) length\n");

    n_e = 0; /* total number of edges */
 
    for( k=0,c=cr->c; c!=NULL; k++,c=c->next ) { /* all cones */
      fprintf(LOG,"\n");
      for( i=0; i<N; i++ )
	for( j=i+1; j<N; j++ ) {
	  /* all edges of cone */
	  ++n_e;
	  length = 1.;
	  length = vec_dist( (c->v[i])->coord,(c->v[j])->coord,N );
	  fprintf(LOG,"[%6d] (%4d,%4d) %10f\n",k,(c->v[i])->index,(c->v[j])->index,length);
	}
      if( c == cr->last_c ) break;
    }

    fprintf(LOG,"\ntotal number of edges = %d\n\n",n_e);
  }

} /* end of db_print_edges() */

/*-----------------------------------------------------------------*/

static void db_print_guide(C_ROOT *cr)
/*
  print guide table statistics
*/
{
  int guide_diff[N+2];
  int i,j;

  /* print differences */
  fprintf(LOG,"\n--------------------------------------------------\n\n");
  fprintf(LOG,"guide table: size = %d\n\n",cr->size_guide);

  /* clear array */
  for( i=0; i<N+2; i++ )
    guide_diff[i] = 0;
  
  /* compute differences between two entries */
  for( j=0; j<cr->size_guide-1; j++ ) {
    if( DEBUG & DB_GUIDE ) fprintf(LOG,"%d ",(cr->guide)[j]->index);
    i =  (cr->guide)[j+1]->index - (cr->guide)[j]->index;
    if( i > N )
	++guide_diff[N+1];
    else
	++guide_diff[i];
  }
  if( DEBUG & DB_GUIDE ) fprintf(LOG,"\n\n");

  fprintf(LOG,"   diff\n");
  for( i=0; i<N+1; i++ )
    fprintf(LOG,"     %2d = %d\n",i,guide_diff[i]);
  fprintf(LOG,"    >%2d = %d\n",N,guide_diff[N+1]);
  
} /* end of db_print_guide() */

/*-----------------------------------------------------------------*/
#if DEBUG > 1000
/*-----------------------------------------------------------------*/

static void db_print_cones_stat_volume(C_ROOT *cr)
/*
  print statistics of volumes of cones
*/
{
  CONE *c;
  double max_Vi, min_Vi, sum_Vi;
  double max_Hi, min_Hi, sum_Hi;

  max_Vi = 0.;
  min_Vi = (cr->c)->Vi;
  sum_Vi = 0.;

  for( c=cr->c; c!=NULL; c=c->next ) {
    max_Vi = MAX(max_Vi,c->Vi);
    min_Vi = MIN(min_Vi,c->Vi);
    sum_Vi += c->Vi;
    if( c == cr->last_c ) break;
  }

  fprintf(LOG,"Volume (min/max/mean) = ( %g / %g / %g )\n",min_Vi,max_Vi,sum_Vi/cr->max_c);
  fprintf(LOG,"\t max / min = %g\n\n",max_Vi/min_Vi);

  max_Hi = 0.;
  min_Hi = (cr->c)->Hi;
  sum_Hi = 0.;

  for( c=cr->c; c!=NULL; c=c->next ) {
    max_Hi = MAX(max_Hi,c->Hi);
    min_Hi = MIN(min_Hi,c->Hi);
    if( c == cr->last_c ) break;
  }

  fprintf(LOG,"Hi (min/max/mean) = ( %g / %g / %g )\n",min_Hi,max_Hi,cr->Htot/cr->max_c);
  fprintf(LOG,"\t max / min = %g\n",max_Hi/min_Hi);
  if( (cr->hc)->reject_tp_rel > 0. ) {
    fprintf(LOG,"\t median = %g (before searching for bad touching points)\n",cr->Hi_median);
    fprintf(LOG,"\t %4.2f-percentile = %g\n",(cr->hc)->reject_tp_percentile,cr->Hi_percentile);
  }

  if( (cr->hc)->reject_tp_abs > 0. )
    fprintf(LOG,"touching point rejected if Hi > %f (runtime option)\n",(cr->hc)->reject_tp_abs);
  if( (cr->hc)->reject_tp_rel > 0. )
    fprintf(LOG,"touching point rejected if Hi > median + %f * (%4.3f-percentile - median) = %f (runtime option)\n",
	    (cr->hc)->reject_tp_rel,(cr->hc)->reject_tp_percentile,
	    cr->Hi_median + (cr->hc)->reject_tp_rel * (cr->Hi_percentile - cr->Hi_median) );

} /* end of db_print_cones_stat_volume() */

/*-----------------------------------------------------------------*/
#endif
/*-----------------------------------------------------------------*/

void db_print_vector(char *begin, double *x, int n, char *end)
/*
  print double vector
*/
{
  int i;

  fprintf(LOG,"%s( %g",begin,x[0]);
  for( i=1; i<n; i++ )
    fprintf(LOG,", %g",x[i]);
  fprintf(LOG," )%s",end);

} /* end of db_print_vector() */

/*-----------------------------------------------------------------*/

void db_print_string(char *c)
/*
  print character string into log file 
*/
{
  fprintf(LOG,"%s",c);
} /* end of db_print_string() */

/*-----------------------------------------------------------------*/

static double vec_dist(double *x, double *y, int n)
/*
  distance of two vectors 
*/
{
  double dist = 0.;
  double diff;

  for( --n; n>=0; n-- ) {
    diff = x[n] - y[n];
    dist += diff * diff;
  }

  return sqrt(dist);
} /* end of vec_dist() */

/*-----------------------------------------------------------------*/
static int *int_zero_alloc(int n)
/*
  allocate pointer to integer array and initialize
*/
{
  int *ia;

  ia = (int *)malloc((size_t) (n*sizeof(int)));
  if (!ia)
    FATAL("allocation failure in int_alloc_zero()");

  /* initialize table */
  for( n-- ; n>=0; n-- )
    ia[n] = 0;

  return ia;
} /* end of int_zero_alloc() */

/*-----------------------------------------------------------------*/
#endif
/*-----------------------------------------------------------------*/