mnc2.c

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

/* 
   mnc2.c                                    (c) DJCM 94 07 04
                                                      94 10 26
						      94 10 31
						      95 01 02

   - free energy minimization for the MacKay-Neal Error Correcting Code - 

   - general purpose program handling many different codes              -

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

   Notes are in the 1995 red book, 95 01 02.

   Data generation:

          For MNC 1-4, The code C is defined by one or two
	  N*N matrices represented explicitly. 

          A second representation is used for the binary matrix A[m][n]
	  involved in the inference problem As + n = z. It is represented by 
	  a pair of arrays of lists. In the `mlist' mlist[m][l], 
	  l=1...num_mlist[m], each element is an `n' such that A[m][n] = 1.
	  Thus the list lists for each m where the 1's are in row m.
	  There is an equivalent `nlist' listing where the 1's are in each 
	  column. 

	  There are four MN codes, each with a different construction
	  for C and / or a different definition of A,s,n,z. 

	  The code is applied to a sparse string s, giving t = Cs. 
	  Noise is added to this, so that Cs + n = z. 
	  But these s, n and z's are not necessarily the same as the 
	  s, n, z of the fe algm.

	  Code 1:

	  to define A = C and plug the problem directly into the 
	  free energy machine, with the approximating distribution 
	  defined on s: Q(s) = prod q_i(s_i). In this case A is an 
	  N * N matrix. 

	  Code 2:

	  The `Radford representation' has separate C into 
	  a product C = C1 C2, where C1 is sparse and C2 has a sparse
	  inverse, maybe of  similar density 
	  to C1.

	  Then we can rewrite the equations

	  Cs + n = z
	  Is + s = 0

	  as 

	  [A1]u + n = z
	  [A2]    s   0

	  where A1 = C1, A2 = C2^{-1}, and u = C2 s. 
	  We can then apply the approximation in the intermediate
	  representation u. The matrix A is an M*N matrix with M=2N. 

	  Each code calls the fe_min routine when things are set up.


          You get to choose (on command line): 
	  N (string length),                                      (e.g. 100)
	  fn  (probability of error in observation of t = As)    (e.g. 0.1)
	  fs  (probability of high bit in s itself)              (e.g. 0.1)
	  seed (for random numbers)



*/

#include "./ansi/r.h"
#include "./ansi/rand.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./fe.h"

typedef struct {
  int index ;
  int val ; 
} thing ; 
static int cmp_thing ( const void * , const void * ) ; 
/* the above for sorting */

typedef struct {
  alist_matrix *a ;
} param;

typedef struct {
  int per_row ;     /* how many 1s to have per row eg 6 */
  int changed ;     /* how many rows got changed during the matrix creation */
} sparse_cmatrix_param ; 

typedef struct {
  int MNC  ;  /* which code to use: 1,2,3,4 = MNC 1-4
		 MNC 1b is equivalent mathematically to MNC 1, 
		 so not available. 
		 if MNC = 0 then an M*N code as in MacKay 1995 is 
		 defined. Except an exact number per row can be invoked.*/
  sparse_cmatrix_param pC1 , pC2 ; 
  int N ;      /* length of the unknown original vector s */
  int M ;      /* length of the transmitted vector t if not N */
  long int mseed ; /* seed for matrices in the code */
  long int vseed ; /* seed for vectors */
  int pbm_o ;  /* whether to write the code matrix as a bit map */
  int pbm_xv ; /* whether to xv it also */
  char pbm_ofile[50] ; /* the pbm file */
  double fpol ; /* fraction of bits high in polynomial (not used) */
  int npol ; 
  double fs ; /* fraction in s */
  double fn ; /* error probability on channel */
  int ns ;          /* if !0, number of signal bits */
  int nn ;          /* if !0, number of noise bits */
  double true_fn ;  /* actual fraction */
  double true_fs ;  /* actual fraction */
  int verbose ;    /* whether the data creation process should be verbose */
  int DEMO ; 
  int tr ; /* typical number in row and column of code */
  int tc ; 
  double fte ; 
  int message ; /* which message we are on */
  int MESSAGE ; /* how many messages to send */
} data_creation_param ;

typedef struct {
  unsigned char **A , **B , **C1 , **C2 , **C1I , **C2I ;
} mnc_code ; 

typedef struct {
  unsigned char *r , *s , *so , *noise , *t , *u , *sno , *uo , *cr ;
  double ft , fte ;  /* empirical fraction of bits in t high */
} mnc_vectors ; 

static void dc_defaults ( data_creation_param * ) ; 
static int    process_command ( int , char **  , data_creation_param * 
			       , fe_min_control * ) ;
static void   print_usage ( char ** , FILE * , data_creation_param * 
			       , fe_min_control * );
static int make_sense ( data_creation_param * , fe_min_control * ) ; 
static int make_code ( data_creation_param * , mnc_code * , fe_min_param * ) ;
static int report_code ( data_creation_param * , mnc_code * , fe_min_param * ) ;
static int make_space ( data_creation_param * , fe_min_param * , 
		       mnc_code * , mnc_vectors * ) ;
static int make_norder ( alist_matrix * , fe_min_control * ) ;
static int make_vectors ( data_creation_param * , fe_min_param * , 
		       mnc_code * , mnc_vectors * ) ;
static int decode ( fe_min_param * , 
			  fe_min_control * , data_creation_param * , 
			  mnc_vectors  * , mnc_code * ) ; 
static int score ( fe_min_param * , 
			  fe_min_control * , data_creation_param * , 
			  mnc_vectors  * , mnc_code * ) ; 
static int make_pvectors ( fe_min_param * , 
			  fe_min_control * , data_creation_param * , 
			  mnc_vectors  * );
static int evaluate_feasibility ( data_creation_param * , fe_min_param * ) ;
static int check_alist ( alist_matrix * , fe_min_param * , 
			fe_min_control * , data_creation_param * ) ;

static void   finalline (  FILE * , fe_min_control * , data_creation_param * , fe_min_param * ) ;

void   main ( int , char ** ) ;

static double h2 ( double ) ;

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

  fe_min_param        p    ;
  fe_min_control      c    ;
  data_creation_param dc   ;
  mnc_code            code ; 
  mnc_vectors         vec  ; 

  fe_defaults ( &p , &c ) ; 
  dc_defaults ( &dc ) ; 

  if ( process_command (argc, argv, &dc , &c) < 0 ) exit (0) ;
  if ( make_sense ( &dc , &c ) < 0 ) exit (0) ;
  if ( make_space ( &dc , &p , &code , &vec ) < 0 ) exit (0) ;

  fprintf(stderr,"mnc MNC=%d, N0=%d, M0=%d, N=%d, M=%d, fn=%6.3g, fs=%6.3g, nn=%d, ns=%d, opt=%d, bstyle=%d\n",
	 dc.MNC , dc.N , dc.M , p.N , p.M ,  dc.fn , dc.fs , dc.nn , dc.ns , c.opt , c.betastyle ) ;
  fflush(stderr);

  if ( make_code ( &dc , &code , &p ) < 0 ) { report_code ( &dc , &code , &p ) ; 
					      exit ( 0 ) ;     }
  if ( check_alist ( &(p.a) , &p , &c , &dc ) < 0 ) { report_code ( &dc , &code , &p ) ; 
						      exit ( 0 ) ;     }
  fe_allocate ( &p , &c ) ; 
  if ( make_norder ( &(p.a) , &c ) < 0 ) exit ( 0 ) ; 
  if ( evaluate_feasibility ( &dc , &p ) < 0 ) exit ( 0 ) ; 

  /* send many messages using a single code */
  for ( dc.message = 1 ; dc.message <= dc.MESSAGE ; dc.message ++ ) {
    if ( make_vectors ( &dc  , &p, &code, &vec) < 0 ) exit (0) ; 
    if ( make_pvectors ( &p , &c , &dc , &vec ) < 0 ) exit ( 0 ) ; 

    fe_min ( &p , &c ) ; 

    if ( decode ( &p , &c , &dc , &vec , &code ) < 0 ) exit ( 0 ) ; 
    if ( score ( &p , &c , &dc , &vec , &code ) < 0 ) exit ( 0 ) ; 
    if ( dc.verbose > 0 ) finalline ( stdout , &c , &dc , &p ) ;

    if ( c.printout ) {
      fp = fopen ( c.outfile , "a" ) ;
      if( !fp ) {
	fprintf( stderr, "No such file: %s\n", c.outfile ) ;
	fp = stderr ;
      }
      else {
	finalline ( fp , &c , &dc , &p ) ;
	fclose ( fp ) ;
      }
    }
  }

  fe_free ( &p , &c ) ; 
}

void finalline ( FILE *fp , fe_min_control *c , data_creation_param *dc , fe_min_param *p ) 
{
  if ( c->pheading_period && ! ( ( c->pheading_count ++ ) % c->pheading_period )  ) 
    fprintf ( fp , 
                 "#es rl.scre  N    t    fs     fn    ns nn  bstl  seed  mseed  vseed tfs      tfn       beta   h     h_dev  MNC ms#\n") ; 

  fprintf ( fp , "%4d %6.4g %3d %3d %7.4g %7.4g %2d %2d %d %8ld %8ld %8ld %7.4g %7.4g %7.4g %7.4g %7.4g %d %4d\n" ,
                  p->count , p->v , dc->N , dc->pC1.per_row , dc->fs , dc->fn , dc->ns , dc->nn , c->betastyle , c->seed  , dc->mseed  , dc->vseed  , dc->true_fs, dc->true_fn , p->beta , p->h_soln , p->h_dev , dc->MNC , dc->message ) ;
}


static int make_sense ( data_creation_param *dc ,  fe_min_control *c ) 
{ /* routine to correct silly control parameters */
/* 
   first the data creation 
                             */
  int status = 0 ; 

  switch ( dc->MNC ) {
  case(0) :
    if ( dc->M == 0 ) {
/* if MNC == 0, M and N better had got set on command line */
      fprintf ( stderr , " No M specified \n" ) ;  status-- ; }
    break ; 
  case ( 1) : 
  case ( 2) : 
  case ( 3) : 
  case ( 4) : 
    dc->M = dc->N ; 
    break ; 
  default : 
    fprintf ( stderr , "invalid MNC = %d\n" , dc->MNC ) ; 
    break ; 
  }

/* 
   and now the free energy control
                                       */
  if ( ( dc->MNC ==3 || dc->MNC == 4 || c->NL > 1 )  &&  c->betastyle == 0 )
    fprintf ( stderr , "Warning: no annealing schedule given, MNC = %d, NL = %d\n" , 
	     dc->MNC , c->NL ) ; 
  if ( c->betastyle == 1 && c->NL > 1 ) 	
    c->betaf = ( c->beta1 - c->beta0 ) / ((double) ( c->NL - 1) ) ; 
  else if ( ( c->betastyle == 2 || c->betastyle == 22 ) && c->NL > 1 )
    c->betaf = exp ( log ( c->beta1 / c->beta0 ) / ((double) ( c->NL - 1) )) ; 
  if ( c->opt > 4 || c->opt < 0 ) {
      fprintf ( stderr , " opt %d not existo\n" , c->opt ) ;   status-- ; }
  if ( c->opt != 3 ) {
    fprintf ( stderr , "opt = 3 override \n" ) ; c->opt = 3 ; 
  }

  return status ; 
}

static int make_space ( data_creation_param *dc , fe_min_param *p , 
		       mnc_code *code , mnc_vectors *vec ) {
  int status = 0 ; 
  int tr , tc ; 

  if ( dc->verbose ) printf ( "Making space\n" ) ; 
  /*
     set up size of the A matrix
                                   */
  switch ( dc->MNC ) {
  case ( 0 ) : 
    p->M = dc->M ;     p->N = dc->N ; 
    break ; 
  case ( 1) :
    p->M = dc->N ;     p->N = dc->N ; 
    break ; 
  case( 2) :
    p->M = 2 * dc->N ; p->N = dc->N ; 
    break ; 
  case(3) : case ( 4 ) : 
    p->M = dc->N ;     p->N = 2 * dc->N ; 
    break ; 
  default:
    fprintf ( stderr , "invalid MNC = %d\n" , dc->MNC ) ; status -- ; 
    break ; 
  }

/* 
   Now estimate size of alist matrix 

     MNC1: Number per row and column is almost fixed to dc.pC1.per_row
     MNC2: Per row is the max of C1 and C2. Per col is the sum plus a bit.
     MNC3: As MNC1, plus 1. 
     MNC4: converse of MNC2.
     MNC0: identity matrix followed by a load of sparse junk 
*/
  switch ( dc->MNC ) {
  case ( 1 ) : 
    tr = dc->pC1.per_row + 2 ; 
    tc = tr + 10 ; break ;
  case ( 2 ) : 
    tr = MAX ( dc->pC1.per_row , dc->pC2.per_row ) + 2 ; 
    tc = 2 * tr + 10 ; break ;
  case ( 3 ) : 
    tr = dc->pC1.per_row + 3 ; 
    tc = tr + 10 ; break ;
  case ( 4 ) : 
    tc = MAX ( dc->pC1.per_row , dc->pC2.per_row ) + 10 ; 
    tr = 2 * tc ; break ;
  default : 
  case ( 0 ) :
    tr = dc->tr + 1 ;
    tc = 1 + ( p->M - p->N ) * tr * 3 / ( 2 * p->N )  ; /* 3/2 for safety */
    break ;
  }
  dc->tr = tr ; dc->tc = tc ; 

  /* Set up memory , etc. */

  vec->s = cvector ( 1 , dc->N ) ;
  if ( dc->MNC == 2 || dc->MNC == 4 ) 
    vec->u = cvector ( 1 , dc->N ) ;
  if ( dc->MNC == 2 ) 
    vec->uo = cvector ( 1 , dc->N ) ; /* u-out */
  else if ( dc->MNC == 4 || dc->MNC == 3 )
    vec->sno = cvector ( 1 , p->N ) ;

  vec->so = cvector ( 1 , dc->N ) ; /* s-out */
  vec->t = cvector ( 1 , dc->M ) ; 
  vec->noise = cvector ( 1 , dc->M ) ; 
  vec->r = cvector ( 1 , dc->M ) ; 
  if ( dc->MNC == 4 ) vec->cr = cvector ( 1 , dc->M ) ; 

  initialize_alist ( &(p->a) , p->N , p->M , tc , tr ) ;