gh.c

快速傅立叶变换程序代码,学信号的同学,可要注意了

  p->N = vec->N ;
  p->z = vec->z ; 
  /* where the true vector lives */
  p->x = vec->x ; 
  p->bias = vec->bias ;
  /* where the output lives */
  p->xo = vec->xo ; 
  return status ; 
}

int check_alist_MN ( alist_matrix *a , mnc_vectors *v ) {
  int status = 0 ; 
  if (     v->N != a->N 
      ||   v->M != a->M      ) {
    fprintf ( stderr , "eek %d %d %d %d\n" , v->N , a->N , v->M , a->M ) ; 
    status -- ; 
  }
  return status ; 
}

/* I adopt the convention that z = 0 always */
static int score ( data_creation_param *dc , 
       		   mnc_vectors  *v , mnc_all *all 
) /* compares so with s and works out quality of solution if wrong. */
{
  int status = 0 ; 
  int n , count = 0 ; 
  FILE *fp ;
 int *nsum = v->nsum ; 

  if ( dc->verbose ) printf ( "Scoring answer --- " ) ; 
  v->count = 0 ; v->count_high = 0 ; v->viols = all->bndp->count_viol ;
  for ( n = 1 ; n <= v->N ; n++ ) { /* this compares all bits */
    if ( v->xo[n] != v->x[n] ) {    
      v->count ++ ;  /* count is the errors */
    }
    v->count_high += v->xo[n] ; /* count_high is the number of highs returned */
  }
  if ( v->count > 0 ) { /* this doesn't look right to me?  aug 97 */
    alist_times_cvector_mod2 ( (all->bndp->a) , v->xo , v->t ) ; 

    if ( dc->verbose ) printf ( "fail\n" ) ; 
    dc->failcount ++ ;
    /* added 99 08: */

    dc->bit_errs += v->count ;
    if ( v->viols ) {
      dc->block_det ++ ;
      dc->bit_det += v->count ;
      /* check for near-codeword */
      if ( v->count < dc->lowweight ) {
	dc->block_detlw ++ ;
      }
    } else {
      printf ( "UNDETECTED ERROR\n" ) ; 
      dc->block_undet ++ ;
      dc->bit_undet += v->count ;
    } 

    
    if ( dc->maxcount != 0 && dc->error_log ) {
      fp = fopen ( dc->error_logfile , "a" ) ;
      if ( !fp ) {
	fprintf ( stderr , " couldn't open logfile %s\n" , dc->error_logfile ) ; 
	dc->error_log = 0 ; 
      } else {
	/* tell me more about this failure */ 	    /* write to a file */
	finalline ( fp , all , dc ) ;
	fprintf ( fp , "# " ) ;
	for ( n = 1 ; n <= v->N && ( ( count <= dc->maxcount ) 
				    || dc->maxcount < 0 ) ; n++ ) { 
	  if ( v->xo[n] != v->x[n] ) {    
	    count ++ ;
	    if ( v->xo[n] ) {
	      fprintf ( fp , "+%4d " , n ) ; /* + indicates extra bit given in so */
	    } else {
	      fprintf ( fp , "-%4d " , n ) ; /* - indicates bit lacking     */
	    }
	  }
	}
	fprintf ( fp , "\n" ) ;
	fclose (fp ) ; 
      }
      /* added 99 08 */
            if ( dc->noise_sum ) { 
/* Keep track of the noise vectors that caused the errors (mean only) */
	fp = fopen ( dc->noise_sumfile , "w" ) ;
	if ( !fp ) {
	  fprintf ( stderr , " couldn't open logfile %s\n" , 
		   dc->noise_sumfile ) ; 
	  dc->noise_sum = 0 ; 
	} else { /* pop the total number of failures and number of bits wrong
		     at the top in nsum[0], nsum[-1] */
	  nsum[0] ++ ;
	  for ( n = 1 ; n <= v->N ; n++ ) { 
	    if ( v->x[n] ) {    
	      nsum[n] ++ ; nsum[-1] ++ ;
	    }
	  }
	  fprintf ( fp , "#noise %d blocks %d bits %d block_error_rate %d / %d \n" , nsum[0] , nsum[-1] , v->N  , dc->failcount, dc->message ) ;
	  for ( n = 1 ; n <= v->N ; n++ ) { 
	    fprintf ( fp , "%d" , nsum[n] ) ;
	    fprintf ( fp , "\n" ) ;
	  }
	  fclose (fp ) ; 
	}
      }
    }
  }  else {
    if ( dc->verbose ) printf ( "SUCCESS\n" ) ; 
  }
  /* added 99 08  */
  if ( !(v->viols) ) {
    dc->block_valid ++ ;
    dc->totloops += all->bndc->loop ; 
    v->histo[all->bndc->loop] ++ ;
  }
  
  /* could also work out the typicality of the solution */
  return status ; 
}

/* needs personalising for fe or depersonalising for bnd */
int make_vectors_quick ( data_creation_param *dc , mnc_vectors *v , mnc_all *all )
{
  int status = 0 ; 
  int count_s = 0 , count_n = 0 , i ; 
  int verbose = 0 ;

  if ( verbose  ) { printf("making a bias vector:\n" ) ; }
  if ( dc->gc ) { /* Here is assumed that using bnd */
    count_n = make_gaussian_noise_bits_and_fix_biases 
      ( v->x , v->bias , dc->gcx , v->nnfrom , v->nnto ) ;
    /* this fixes x to zero */
  } else { /* BSC */
    if ( dc->nsn ) { /* a fixed weight vector is required */
      count_n = fixed_wt_cvector ( v->x , dc->nsn , v->nsfrom , v->nnto ) ;
    } else { /* use fn other data generation rules */
      count_n = random_cvector ( v->x , dc->fn , v->nnfrom , v->nnto ) ; 
    }
    for ( i =  v->nnfrom ; i <=  v->nnto ; i ++) {
      if ( v->x[i] ) {
	v->bias[i] = 1.0 - dc->fn ;
	v->x[i] = 0 ; /* NB I am fixing the true x to be all zeroes */
      } else {
	v->bias[i] =  dc->fn ;
      }
      if ( verbose ) {
	printf ( "%9.4g " , v->bias[i] ) ;
	if ( !(i%8) ) printf ("\n");
      }
    }
  }
  v->count_s = count_s + count_n ;
    
  dc->true_fs = 0.0 ;
  dc->true_fn = (double) v->count_s / (double) dc->M ; 
    
  return status ; 
}
int make_vectors_quickER ( data_creation_param *dc , mnc_vectors *v , mnc_all *all ) /* just makes the calls to the rand number generator */
{
  int status = 0 ; 
  int count_s = 0 , count_n = 0 ; 

  if ( !dc->ns ) {
    count_s = random_cvector ( v->x , dc->fs , v->nsfrom , v->nsto ) ; 
  } else {
    fixed_wt_cvector (  v->x , dc->ns , v->nsfrom , v->nsto ) ;
  }

  if ( dc->gc ) { /* Here is assumed that using bnd */
    count_n = make_gaussian_noise_bits_ONLY 
      ( v->x , v->bias , dc->gcx , v->nnfrom , v->nnto ) ;
  }
  else {
    if ( !dc->nn ) {
      count_n = random_cvector ( v->x , dc->fn , v->nnfrom , v->nnto ) ; 
    } else {
      fixed_wt_cvector (  v->x , dc->nn , v->nnfrom , v->nnto ) ;
    }
  }

  return status ; 
}

static int make_gaussian_noise_bits_and_fix_biases ( unsigned char *x ,
						    double *b,
						    double gcx ,
						    int lo , int hi ) {
  int i , c = 0 ;
  double z , p ;

  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    z = gcx + rann() ; 
    p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcx ) ) ;  /* likely to be large */
    x[i] = 0 ; 
    b[i] = 1.0 - p ; 
    if ( p < 0.5 ) c++ ;
  }
  return c ; 
}

static int make_gaussian_noise_bits_ONLY ( unsigned char *x ,
						    double *b,
						    double gcx ,
						    int lo , int hi ) {
  int i , c = 0 ;
  double z ;

  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    z = rann() ; 
  }
  return c ; 
}

static void set_up_priors ( double *b , mnc_vectors *v , data_creation_param *dc ) {
  int n ;

  for ( n = v->nsfrom ; n <= v->nsto ; n++ ) {
    b[n] = dc->fs ;
  }
  for ( n = v->nnfrom ; n <= v->nnto ; n++ ) {
    b[n] = dc->fn ;
  }
}
  
static int evaluate_feasibility ( data_creation_param *dc )    
{
  int status = 0 ;
/*
   compute Shannon capacity and compare with info content  
*/

  dc->h2fs = h2(dc->fs) ;
  dc->h2fn = h2(dc->fn) ;
  dc->rho =  (double ) dc->bM / (double ) dc->N  ;
  dc->rate = dc->rho * dc->h2fs ;
  dc->capacity = 1.0 - dc->h2fn ;
  if ( !dc->gc ) {
    if ( dc->rate >   dc->capacity ) {
      printf ( "task impossible, signed Shannon\n" ) ; 
      if ( dc->notabovecap ) {status -- ; }
    }
    if ( status < 0 || dc->verbose ) 
      printf ( "H(s) = %f, capacity of channel = %f.\n" , 
	      dc->h2fs , dc->capacity ) ; 
  }
  return status ; 
}


static void dc_defaults ( data_creation_param *dc ) 
{
#include "RMdc_var_def.c"
}

static int process_command ( int argc , char **argv , mnc_all *all ) {
  data_creation_param *dc = all->dc ; 
  fe_min_control *c = all->c ; 
  bnd_control *bndc = all->bndc ; 

  int p_usage = 0 ;
  int status = 0 ;
  int cs , i ;

  if ( argc < 1 )     {
    p_usage = 1 ; 
    status -- ;
  }

#define ERROR1 fprintf ( stderr , "arg to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERROR2 fprintf ( stderr , "args to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERRORREG fprintf ( stderr , "regtype must be defined before `%s'\n" , argv[i] ) ; \
               status --
  for (i = 1 ; i < argc; i++)    {
    cs = 1 ;
    if ( strcmp (argv[i], "-V") == 0 )        {
      c->verbose = 1;
      dc->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-VV") == 0 )        {
      c->verbose = 2;
      dc->verbose = 2;
    }
#include "fe_var6_clr.c"
#include "RMdc_var_clr.c"
#include "RMbnd_var_clr.c"
    else if ( strcmp (argv[i], "-b") == 0 ) {
      if ( i + 3 >= argc ) { ERROR2;      }
      else {
	cs *= sscanf(argv[++i], "%d", &(c->betastyle)); 
	cs *= sscanf(argv[++i], "%lf", &(c->beta0)); 
	cs *= sscanf(argv[++i], "%lf", &(c->beta1)); 
      }
    }
    else {
      fprintf ( stderr , "arg `%s' not recognised\n" , argv[i] ) ; 
      p_usage = 1 ;
      status -- ;
    }
    if ( cs == 0 ) {
      fprintf ( stderr , "arg at or before `%s' has incorrect format\n" , 
	       argv[i] ) ;
      p_usage = 1 ;
      status -- ;
    }
  }
  if ( p_usage ) print_usage ( argv , stderr , all ) ;
  return ( status ) ;
}
#undef ERROR1
#undef ERROR2
#undef ERRORREG

#define DNT fprintf( fp, "\n        ")
#define NLNE  fprintf( fp, "\n")

static void print_usage ( char **argv , FILE * fp ,
			 mnc_all *all 
)
{
  data_creation_param *dc = all->dc ; 
  
  fprintf( fp, "Usage: %s ",argv[0]);
  fprintf( fp, " [optional arguments]");

  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)"); 
  NLNE; fprintf( fp, " Data creation:" ) ; 
#include "RMdc_var_usg.c"
  pause_for_return();
  NLNE; fprintf( fp, " Inference:" ) ; 
  DNT;  fprintf( fp, "-b betastyle beta0 beta1 (what to do with beta)" ); 
  DNT;  fprintf( fp, "       betastyle 0: const; 1: linear; 2: multiply; 22: multiply and go wild on last loop" ); 
  fprintf( fp, "\n");
  fe_print_usage ( argv , fp , all ) ; 
  pause_for_return();
  bnd_print_usage ( argv , fp , all ) ; 
  return ;
}

static void fe_print_usage ( char **argv , FILE * fp ,
			    mnc_all *all 
)
{
  fe_min_control *c = all->c ; 

  NLNE; fprintf( fp, " Further free energy minimization stuff: <defaults>");
#include "fe_var6_usg.c"
  fprintf( fp, "\n");
  return ;
}

static void bnd_print_usage ( char **argv , FILE * fp ,
			    mnc_all *all 
)
{
  bnd_control *bndc = all->bndc ; 

  NLNE; fprintf( fp, " RM Belief Net decoder:     <defaults>");
#include "RMbnd_var_usg.c"
  fprintf( fp, "\n");
  return ;
}

#undef DNT
#undef NLNE

static void mnc_free ( mnc_all *all ) {
/*  data_creation_param *dc = all->dc ;*/
  mnc_vectors *vec = all->vec ; 

  free_cvector ( vec->xo, 1 , vec->N ) ; 
  free_cvector ( vec->t , 1 , vec->M ) ; 

  free_cvector ( vec->x , 1 , vec->N ) ; 
  free_cvector ( vec->y , 1 , vec->M ) ; 
  free_cvector ( vec->z , 1 , vec->M ) ; 

  free_dvector ( vec->bias , 1 , vec->N ) ; 

}

double h2 ( double x ) {
  double tmp ; 
  if ( x <= 0.0 || x>= 1.0 ) tmp = 0.0 ;
  tmp = x * log ( x ) + ( 1.0 - x ) * log ( 1.0 - x ) ;
  
  return - tmp / log ( 2.0) ; 
}

/*
<!-- hhmts start -->
Last modified: Sat Jun  7 17:36:27 1997
<!-- hhmts end -->
*/