gh.c

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

/*
   GH.c
   mncN.c
   mnc6.c                                    (c) DJCM 94 07 04
                                                      94 10 26
						      94 10 31
						      95 01 02
						      95 02 23 (from mnc2 to mnc3)
						      95 04 01 continued work on mnc3
						      95 05 01 to mnc4
						      95 05 27 to mnc5
						      95 08 31 to mnc6
						      95 11 04 to mnc7
						      95 11 16 to mncN
						      97 04 19 to GH

   - mncN = quick hack to do all bits as if they are noise bits.

   - free energy minimization for the MacKay-Neal Error Correcting Code - 
   - or belief net decoder                                              -
   - mnc7 includes option of gaussian channel                           -

   - space is made for the fe min or for the bnd, then the scoring 
     routine gets to look at the results via pointers in the mnc_vector
     structure. Thus scoring is divorced from fe. 

     crucial var files: both fe_var6 and RMdc_var

   ======================================================================

   This code is (c) David J.C. MacKay 1994, 1995. It is free software 
   as defined by the free software foundation. 

   What this program does: 

     reads in a matrix A in alist format.
     do {
       generates a  random vector x, 
       shoves it through A, adds y if appropriate, -> z = Ax + y
       do {
         gives b, g (z) to the minimization machine, or A,z to the bnd
         checks the answer. 
	 next strategy
       }
       until ( success OR no more strategies )
     }
     while ( not many failures )

*/

#include "./ansi/r.h"
#include "./ansi/rand2.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./ansi/macopt.h"
#include "./thing_sort.h"
#include "./fe6.h"
#include "./RM.h" /* defines trelises */
#include "./RMbnd.h"
/* this must come last: */
#include "./GH.h"

static void totalline ( FILE * , mnc_all * , data_creation_param * );
static void dc_defaults ( data_creation_param * ) ; 
static int    process_command ( int , char **  , mnc_all * ) ; 
static void   print_usage ( char ** , FILE * ,  mnc_all *  );
static void   fe_print_usage ( char ** , FILE * ,  mnc_all *  );
static int make_sense ( data_creation_param * , mnc_all * ) ; 
static int make_space ( data_creation_param * , mnc_vectors * , mnc_all * ) ;
static int make_vectors_quick ( data_creation_param * , mnc_vectors * ,  mnc_all * ) ;
static int make_vectors_quickER ( data_creation_param * , mnc_vectors * ,  mnc_all * ) ;
static int make_gaussian_noise_bits_and_fix_biases ( unsigned char * ,
						    double *,
						    double  ,
						    int , int  ) ;
static int make_gaussian_noise_bits_ONLY ( unsigned char * ,
						    double *,
						    double  ,
						    int , int  ) ;
static void set_up_priors ( double * , mnc_vectors * , data_creation_param * ) ; 
static int score ( data_creation_param * , mnc_vectors  * , mnc_all * ) ; 
static int evaluate_feasibility ( data_creation_param * ) ;

static double bern ( int  , int  , double * , double * ,  double * , double );
static void write_histo ( FILE * , mnc_all * ) ;
static void snappyline ( data_creation_param * ) ;
static void   finalline (  FILE * , mnc_all * , data_creation_param * ) ;

static void mnc_free ( mnc_all * ) ; 
static int check_alist_MN ( alist_matrix * , mnc_vectors *vec ) ; 
static int hook_mnc_vec_to_bnd ( bnd_param *p , mnc_vectors *vec ) ;
static void bnd_print_usage ( char **argv , FILE * fp , mnc_all *all );
static double h2 ( double ) ;

int   main ( int , char ** ) ;

/*
        MAIN
                     */
int main ( int argc, char *argv[] )
{
  FILE   *fp  ; 

  fe_min_param        p    ;
  fe_min_control      c    ;
  bnd_control         bndc ; 
  bnd_param           bndp ;
  data_creation_param dc   ;
  mnc_vectors         vec  ; 
  mnc_all             all  ;
  linear_trellis      lt   ;

  all.dc   = &dc ; 
  all.p    = &p ; 
  all.c    = &c ; 
  all.bndp = &bndp ; 
  all.bndc = &bndc ; 
  all.vec  = &vec ; 
  all.a    = &(p.a) ; 
  bndp.a   = &(p.a) ; 
  bndp.lt  = &lt ;

  fe_defaults_p ( &p ) ; 
  fe_defaults_c ( &c ) ; 
  dc_defaults   ( &dc ) ; 
  bnd_defaults  ( &bndp , & bndc ) ; /* this sets the trellis defaults */

  if ( process_command (argc, argv, &all ) < 0 ) exit (0) ;
  if ( make_sense ( &dc , &all ) < 0 ) exit (0) ;

  fprintf(stderr,"GH N0=%d, M0=%d, N=%d, M=%d, x=%6.3g, fn=%6.3g, fs=%6.3g, nn=%d, ns=%d\n",
	  dc.N , dc.M , vec.N , vec.M , dc.gcx , dc.fn , dc.fs , dc.nn , dc.ns) ;
  fflush(stderr);

  if ( make_space ( &dc , &vec , &all ) < 0 ) exit (0) ;
  fprintf(stderr,".");  fflush(stderr);
  make_linear_trellis ( &lt , bndc.trellisN , bndc.trellisT , bndc.trellisV ) ; /* this allocates memory */
  fprintf(stderr,".");  fflush(stderr);

  if ( read_allocate_alist ( &(p.a) , dc.afile ) < 0 ) exit (0) ; 
  fprintf(stderr,"A");  fflush(stderr);
  if ( check_alist_MN ( &(p.a) , &vec ) < 0 ) exit (0) ; 
  fprintf(stderr,".");  fflush(stderr);

  if ( dc.use_bnd >= 1 ) {
    hook_mnc_vec_to_bnd (  &bndp ,  &vec ) ;
    bnd_allocate ( &bndp , &bndc ) ; 
    set_up_priors ( bndp.bias , &vec , &dc ) ; 
  }
  fprintf(stderr,".");  fflush(stderr);

  if ( evaluate_feasibility ( &dc ) < 0 ) exit ( 0 ) ; 
  fprintf(stderr,".");  fflush(stderr);

  if ( c.writelog ) {
    fp = fopen ( c.logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , c.logfile ) ; 
      c.writelog = 0 ; 
    } else fclose (fp ) ; 
  }

  fprintf(stderr,".");  fflush(stderr);
  if ( bndc.writelog ) {
    fp = fopen ( bndc.logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , bndc.logfile ) ; 
      bndc.writelog = 0 ; 
    } else fclose (fp ) ; 
  }

  fprintf(stderr,"*\n");  fflush(stderr);
  if ( dc.error_log ) {
    fp = fopen ( dc.error_logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , dc.error_logfile ) ; 
      dc.error_log = 0 ; 
    } else fclose (fp ) ; 
  }

  /*
      MAIN LOOP
                 */
  /* I will adopt the convention that z = 0 always */

  ran_seed ( dc.vseed ) ; 
  dc.message = 1 ;
  if ( dc.skip ) {
    fprintf ( stderr , "skipping %d", dc.skip ) ; 
    for (  ; dc.message <= dc.skip ; dc.message ++ ) {
      make_vectors_quickER ( &dc , &vec , &all ) ; 
      if ( !(dc.message % 200) ) { fprintf ( stderr , "." ) ; fflush ( stderr ) ; }
    }fprintf ( stderr , "\n" ) ; 
  }
  for ( ;       ( dc.message <= dc.MESSAGE ) && 
       ( ( dc.failures==0 ) || ( dc.failcount < dc.failures ) ) ; 
       dc.message ++ ) {

    snappyline( &dc ) ; /* added august 99 from gallager.c */
    fflush(stdout) ; 
    make_vectors_quick ( &dc , &vec , &all ) ; 

    all.controller = 0 ; 
    bndecode ( &bndp , &bndc ) ;

    if ( score ( &dc , &vec , &all ) < 0 ) exit ( 0 ) ; 

    if ( dc.verbose > 0 ) finalline ( stdout , &all , &dc ) ;

    if ( c.printout ) { /* append */
      fp = fopen ( c.outfile , ((c.outappend)? "a":"w") ) ;
      if( !fp ) {
	fprintf( stderr, "No such file: %s\n", c.outfile ) ;
	finalline ( stderr , &all , &dc ) ;
      }      else 	{
	finalline ( fp , &all , &dc ) ;
	fclose ( fp ) ;
      }
    }
    if ( c.printtot ) { /* write */
      fp = fopen ( c.totoutfile , "w" ) ;
      if( !fp ) {
	fprintf( stderr, "No such file: %s\n", c.totoutfile ) ;
	totalline ( stderr , &all , &dc ) ;
      }      else 	{
	totalline ( fp , &all , &dc ) ;
	fclose ( fp ) ;
      }
    }
          if ( c.printhisto && dc.block_valid ) { /* update histogram file */
	fp = fopen ( c.histofile , "w" ) ;
	if( !fp ) {
	  fprintf( stderr, "No such file: %s\n", c.histofile ) ;
	}      else 	{
	  write_histo ( fp , &all ) ;
	  fclose ( fp ) ;
	}
      }

  }

  if ( dc.use_bnd >= 1 ) {
    bnd_free ( &bndp , &bndc ) ;
  }
  mnc_free ( &all ) ; 
  return (1) ; /* Sat 15/12/01 got warning: return type of `main' is not `int'
		  so made change of type of main to int */
}

static void write_histo ( FILE *fp , mnc_all *c ) {
  int l;
  double t , cum = 0.0 ;
  double  tot = (double) c->dc->block_valid ; 
  int *histo = c->vec->histo ; 
  int loops = c->bndc->loops ;

  fprintf ( fp , "# total valid blocks %d\n" , (int) tot ) ; 
  for ( l = 1 ; l <= loops ; l ++ ) {
    t= (double) histo[l] ;
    cum += t ;
    fprintf ( fp , "%d\t%d\t%d\t%9.4g\t%9.4g\n" , l , (int)(t) , (int)(cum) , 
	      t/tot , cum/tot ) ; 
  }

}
static void snappyline ( data_creation_param  *c ) {
  printf ( "%d:%du%dd%dl/%d\t" , c->failcount ,
	   c->block_undet, c->block_det,
	   c->block_detlw , c->message - 1 ) ;
 fflush ( stdout ) ; 
}


static void finalline ( FILE *fp , mnc_all *all , data_creation_param *dc )
{
  mnc_vectors *v = all->vec ;
  fe_min_control *c = all->c ; 
  if ( (!(c->outappend)) || ( dc->pheading_period && 
      ! ( ( ++ dc->pheading_count ) % dc->pheading_period ))  ) 
    fprintf ( fp , 
                 "#es tr_wt recn_wt viols N   M    t    fs     fn    gc nn bstl  seed  mseed  vseed rho rate capacity  gcx      tfn       beta nl<fe>iter bndloops  ms# failures\n") ; 

  fprintf ( fp , "%4d %4d %4d %4d %6d %6d %d " , 
	   v->count , v->count_s , v->count_high , v->viols , 
	   dc->N , dc->M , dc->tc ) ;
  fprintf ( fp , "%7.4g %7.4g %2d %2d %d %1ld %8ld %8ld " ,
	   dc->fs , dc->fn , dc->gc , dc->nn , c->betastyle , 
	   c->seed  , dc->mseed  , dc->vseed  );
  fprintf ( fp , "%7.4g %7.4g %7.4g " ,
	   dc->rho , dc->rate , dc->capacity ) ; 
  fprintf ( fp , "%7.4g %7.4g %7.4g %3d %3d %3d %5d %2d\n" , 
	   dc->gcx, dc->true_fn , all->p->beta ,all->c->nl , 
	   all->c->iter , all->bndc->loop , dc->message , dc->failcount ) ;
}

static void totalline ( FILE *fp , mnc_all *all , data_creation_param *dc )
/* records total number of its used, total number of failures of type I
   (there exist violations and bits wrong; maximum number of its was used)
   and of type II (no violations but some bits wrong; early stop)
   total number of bots wrong (in both categories).

   Also computes error bars on the error rates
   and reports the value of Eb/No

   number of failures
   number of trials
   num fail type I
   num fail type II
   L M Rate (K/N, K=L-M, N=L)
   x
   Capacity of channel at this x
   Eb/No
   Distance from shannon limit for this rate
   itmax
   number of its used in successes
   number of its used in failures
   number of bits wrong
   num bits type I 
   num bits type II
   block error rate: ML, point, upper, lower
   bit error rate: ML, point, upper, lower

   THERE ARE GRAND INTENTIONS HERE< BUT NOT IMPLEMENTED!
*/
{
  mnc_vectors *v = all->vec ;
  fe_min_control *c = all->c ; 
  if ( dc->pheading_period && 
      ! ( ( ++ dc->pheading_count ) % dc->pheading_period )  ) 
    fprintf ( fp , 
                 "#es tr_wt recn_wt viols N   M    t    fs     fn    gc nn bstl  seed  mseed  vseed rho rate capacity  gcx      tfn       beta nl<fe>iter bndloops  ms# failures\n") ; 

  fprintf ( fp , "%4d %4d %4d %4d %6d %6d %d " , 
	   v->count , v->count_s , v->count_high , v->viols , 
	   dc->N , dc->M , dc->tc ) ;
  fprintf ( fp , "%7.4g %7.4g %2d %2d %d %1ld %8ld %8ld " ,
	   dc->fs , dc->fn , dc->gc , dc->nn , c->betastyle , 
	   c->seed  , dc->mseed  , dc->vseed  );
  fprintf ( fp , "%7.4g %7.4g %7.4g " ,
	   dc->rho , dc->rate , dc->capacity ) ; 
  fprintf ( fp , "%7.4g %7.4g %7.4g %3d %3d %3d %5d %2d\n" , 
	   dc->gcx, dc->true_fn , all->p->beta ,all->c->nl , 
	   all->c->iter , all->bndc->loop , dc->message , dc->failcount ) ;
}

/* be nice to record the number of iterations used? */

static int make_sense ( data_creation_param *dc , mnc_all *all ) 
{ /*  correct silly control parameters */
/* 
   first the data creation  - note M,N conventions changed relative to mncN
                             */
  int status = 0 ; 

                                 /* there is potential confusion 
				    here:

				    the matrix is A is M * N
				    */				    

  all->vec->M = dc->M ; 
  all->vec->N = dc->N ; 
  all->vec->NS = 0 ; 
  all->vec->NN = dc->N ; 
  all->vec->nsfrom = 1 ; 
  all->vec->nsto   = 1 ; 
  all->vec->nnfrom = 1 ; 
  all->vec->nnto   = dc->N  ; 

  if ( dc->bM == 0 ) { /* assuming only one trellis */
    dc->bM = dc->M * all->bndp->lt->M ; 
  }

  return status ; 
}

static int make_space ( data_creation_param *dc , mnc_vectors *vec  , mnc_all *all ) {
/* see also mnc_free */
  int status = 0 ; 

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

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

  vec->bias = dvector ( 1 , vec->N ) ;
  vec->n = dvector ( 1 , vec->N ) ;
  vec->nsum = ivector ( -1 , vec->N ) ;
  set_ivector_const ( vec->nsum ,  -1 , vec->N , 0 ) ;

  vec->histo = ivector ( 0 , all->bndc->loops+1 ) ; /* histogram of number of loops taken  */
  set_ivector_const ( vec->histo ,  0 , all->bndc->loops+1 , 0 ) ; 

  return status ; 
}

static int hook_mnc_vec_to_bnd ( bnd_param *p , mnc_vectors *vec ) {
  int status = 0 ; 
  p->M = vec->M ;