mnc.c

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

/* 
   mnc.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 - 

   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.

   This is the first really working version using the bloj construction,
   and I am running a load of tests with it. 950107

   Data generation:

          The code C is defined by an N*N matrix. This matrix is 
	  not computed explicitly but is represented in two ways;
	  first, as a `ulist', i.e., a product of a sequence 
	  of `U-matrices' such as
	            100000
	            010000
	  U(5,3) =  001000
	            000100
	            001010
	            000001  (such matrices are self-inverse, so  
		             the inverse of C = U1 U2 .. UK is just
			     UK U_K-1 .. U2 U1.)

          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. 

	  The code is applied to a sparse string s, giving t = Cs. 
	  Noise is added to this, so that Cs + n = z. 

	  We can attempt to solve for s in two ways. The first is 
	  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. 

	  The alternative `Radford representation' is to separate C into 
	  a product C = C1 C2, where C1 and C2 might have similar density 
	  to each other, and smaller density than C. 

	  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. 


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


   Inference:

    An optimizer then runs to try to infer s. 

          optimizers available are:
	  frp (standard conjugate gradients algm)
	  macopt (my possibly-faster cg algm)
	  jumpopt (proceeds by re-estimation formula instead of 
	           creeping downhill) 
	  seqopt (sequential re-estimation so that decrease in F is 
	          guaranteed)
	  All of these optimizers can optionally be run at temperatures
	  other than beta = 1.  

    A summary report is written. There are three possible outcomes:

          1  answer is correct
          0  answer is `wrong', and has higher energy than true state
	  2  `ok' answer is wrong, but has lower energy, so is (conditionally)
	       a perfectly good answer to an ill-posed problem. 

	       These are distinguished by the values of count (0 if right)
	       and v (negative if `ok')
	       
    Output of this program can be turned into convenient gnuplot ps files
    by running codeperfms.p outputfile | gnuplot
    The ps file is named by munging the name of outputfile

*/

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

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

typedef struct {
  int N ;
  int K ; 
  int **C ;
} ulist_matrix ;

typedef struct {
  int N ;
  int u ; /* whether this is an upper triangular matrix, and its inverse, 
	     or vice versa */
  int i ; /* whether the inverse has been computed */
  unsigned char **m ;
  unsigned char **mi ; /* for the inverse */
} triangular_matrix ;

typedef struct {
  int style ;
  double density ; 
  double *f ;
} te_params ; 

typedef struct {
  int N , M ;
  int **mlist;
  int **nlist;
  int *num_mlist;
  int *num_nlist;
  int *l_up_to ;
  int *norder ;
  int biggest_num_m ;
  int biggest_num_n ; 
  int tot ; 
} alist_matrix ;

typedef struct {
  alist_matrix *a ;
} param;

static int    process_command ( int , char ** ) ;
static void jump_opt (double * ) ; /* optimization by `reestimation' */
static void seq_opt (double * , param * ) ; /* optimization by `reestimation' */
static void   print_usage ( char ** , FILE * );
static void make_sense ( void ) ; 
static void print_state ( double , double , double * )  ;
static void   printall ( FILE * )  ;
static void   finalline ( int , double , FILE * )  ;
static double forward_pass ( void ) ;
static void   backward_pass ( double * ) ;
void mega_forward_pass ( void ) ;
void mega_backward_step ( int );  /* nn should start at N and go backwards */
static void find_q_S ( double * , int ) ;
static double x_to_q ( double , double * , double * , int )  ;
static double objective ( double *, void * ) ;
static void   d_objective ( double *, double *, void * ) ;
void   main ( int , char ** ) ;
unsigned char cdotprod_mod2_index (  int * , unsigned char * , int  , int  ) ;
int cdotprod_index (  int * , unsigned char * , int , int ) ;
static double h2 ( double ) ;
static void finish_off_alist ( alist_matrix * ) ;
static void cmatrix_2_alist ( unsigned char ** , alist_matrix * , int ) ;
static void cmatrices_2_alist ( unsigned char ** , unsigned char ** ,
			       alist_matrix * , int ) ;
static void initialize_alist ( alist_matrix * , int , int , int , int ) ;

/* making the data */ 

static int  num_per_row = 6 ; 
static int ustyle = -1 ; /* how to make the pseudo-ranomd U matrices */
static int DEMO = 0 , LUDEMO = 0 ; 
static int n , N = 30; /* runs over state bits */
static int m , M ; /* runs over the data */
static alist_matrix alist ;
static int NORDER = 0 ; 
static unsigned char *s ;    /* the true binary state vector s[n] */
static unsigned char *t ;    /* the binary vector t[m] = A s mod 2 */
/* static double *er ;    the error probability vector er[m] */
/* not used in this program */
static double *g ;  /* log odds g[m] = +/- log er[m]/(1-er[m]) depending
		       whether z[m] = 1 or 0 . */
static double fpol = 0.1 ;
static double true_err ;   /* actual fraction */
static double true_fs ;   /* actual fraction */
static double fs = 0.2 ;  /* fraction of bits in s that are high */ 
static double err = 0.2 ;       /* default error probability */
static double dens = 0.05 ; 
static double ft , fte ; /* empirical fraction of bits in t high */

static int K=30 ; /* number of [[1,1],[0,1]] matrices that are multiplied 
		      to give A */

/* the model */

/* static double *x ;  x[n=1..N] 
		      the real log probs (also known as theta in my paper ) 
		      This one is not global because it gets created by the 
		      optimizer */
static double sdx = 0.0 ; /* initial standard deviation */
static double *bias ; /* the prior on the s vector */
static double *q0 , *q1 ; /* q1 = 1 / 1+ exp (- x) , etc. */
static double *p0 , *p1 ; /* probabilities used in the forward and backward
			     passes, p0[0..N+1], p1[0..N+1] */
static double **mp0 , **mp1 ; /* mega-probabilities mp0[1..M][0..N+1] */
static unsigned char *so ; /* the answer so[1..N] */
double S , E , P ; /* entropy and energy of solution ; P is the prior energy bit */

/* control */

static double beta = 1.0 ;  /* the temperature of the party */
static double beta1 = 1.0 ; 
static double betaf ;       /* if annealing is used then 
			       beta goes from current value of 
			       beta to beta1 by uniform steps 
			       either arithmetic or geometric */
static double betap = 1.0 ; /* factor for the prior bit of the energy, */

static int betastyle = 0 ;  
static int     verbose = 0 ;
static int     seqverbose = 0 ;
static int CG=0; /* CG = 1 causes gradient to be checked */
static int nl , NL = 1 ; /* Number of loops. (annealing happens within these) */
static int rich = 0 ; 
static int opt = 3 ;  /* NB opt must be 3 */
static int iter , itmax = 100 ; 

static  char 	outfile[50] = "_out" ;
static  char 	pbm_ofile[50] = "A.pbm" ;
static  int     printout = 0 ;
static  int          pbm_o = 1 , pbm_xv = 0 ; 

static  long int seed = 234 ;
static  double  ftol=0.0001;
static  double  flintol=0.001;
static  double  epsilon = 0.0001 ;

static  int how_to_gen_A = 3 ; /* 1 for ulist, 2 for LU , 3 for bloj construction */

/* sequential stuff */
static int seq_period = 10 ; /* how frequently the free energy is evaluated 
				to decide whether to 
				stop the sequential optimizer */

/* monte carlo bit - cut */
/* junk */
static double pheading=0.1; 

static int Radford_rep = 0 ; /* Whether to work in the intermediate rep between
				noise space and string space */

void main ( int argc, char *argv[] )
{
  param control; 
  double *x , v ;
  int count , count_high  , count_err , count_s , count_t ;
  FILE *fp ; 
  int confirm_inverse = 0 ; 
  int *u ;

/* new things for mnc */

  double tmpd ; 
  int tmpi ; 

  unsigned char *t2 ;   /* the binary vector t[m] = A s mod 2 */

/* new things for ms */
  unsigned char *z , *noise , **A , **B , **C1 , **C2 , **C1I , **C2I ;
/*  int i , l , u  ; */

  char 	junk[100] ; 

  if ( process_command (argc, argv) < 0 ) exit (0) ;
  make_sense ( ) ;
  if ( Radford_rep ) M = 2 * N ; 
  else M = N ; 

  fprintf(stderr,"mnc N=%d, M=%d, err=%6.3g,  fs=%6.3g, opt=%d, bstyle=%d\n",
	 N , M ,  err , fs , opt , betastyle ) ;
  fflush(stderr);

  /* Set up memory , etc. */
  s = cvector ( 1 , N ) ;
  u = cvector ( 1 , N ) ;
  so = cvector ( 1 , N ) ;
  t = cvector ( 1 , M ) ; /* If Radford_rep, then M = 2 * N, 
			     else M = N. 
			     In the MNC_DEMO,  we just use the 
			     first half of it */
  noise = cvector ( 1 , N ) ; 
  z = cvector ( 1 , M ) ; /* this'll be the observation */
  if ( Radford_rep ) {
    initialize_alist ( &alist , N , M , 3 + num_per_row , 8 * num_per_row  ) ;
  } else {
    initialize_alist ( &alist , N , M , 3 + num_per_row , 4 * num_per_row  ) ;
  }

  ran_seed ( seed ) ; 

  if ( Radford_rep ) {
    C1 = cmatrix ( 1 , N , 1 , N ) ; 
    C1I = cmatrix ( 1 , N , 1 , N ) ; 
    sparse_invertible_cmatrix ( C1 , C1I , num_per_row , N ) ;
  
    C2 = cmatrix ( 1 , N , 1 , N ) ; 
    C2I = cmatrix ( 1 , N , 1 , N ) ; 
    sparse_invertible_cmatrix ( C2I , C2 , num_per_row , N ) ;
  
    cmatrices_2_alist ( C1 , C2I , &alist , N ) ;
  } else {
  
    A = cmatrix ( 1 , N , 1 , N ) ; 
    B = cmatrix ( 1 , N , 1 , N ) ; 
    sparse_invertible_cmatrix ( A , B , num_per_row , N ) ;
    
    cmatrix_2_alist ( A , &alist , N ) ;
    if ( pbm_o ) {
      fp = fopen ( pbm_ofile , "w" ) ;
      cmatrix2pbm ( A , 1 , N , 1 , N , fp ) ; 
      fclose ( fp ) ; 
      if ( pbm_xv ) {
	sprintf ( junk , "xv %s -geometry +10+10 &" , pbm_ofile ) ;
	system ( junk ) ;
      }
    }
  }  /* end if radford_rep == 0 */

  /* Use A to generate t = A s */
    
  if( DEMO || verbose >= 1 ) printf ( "generating string\n" ) ;
  count_s = 0 ; 
  for ( n = 1 ; n <= N ; n++ ) {
    t[n] = s[n] = ( ranu() < fs ) ? 1 : 0 ; 
    count_s += s[n] ; 
  }
  if ( Radford_rep ) {
    mult_cm_cv ( C2 , s , u , 1 , N ) ;
    mult_cm_cv ( C1 , u , t , 1 , N ) ;
  } else {
    mult_cm_cv ( A , s , t , 1 , N ) ;
  }
  for ( n = 1 ; n <= N ; n++ ) {
    count_t += t[n] ; 
  }
  fte = (double) count_t / (double) N ; 
  if( DEMO || verbose >= 1 )   printf ( "generating noise - received:\n" ) ;
  count_err = 0 ; 
  for ( n = 1 ; n <= N ; n++ ) {
    count_err += ( noise[n] = ( ranu() < err ) ? 1 : 0 ) ; 
    z[n] = noise[n] ^ t[n] ;
  }
  true_err = (double) count_err / (double) N ; 
  true_fs = (double) count_s / (double) N ; 
  
  if ( DEMO || verbose >= 2  )   {
    printf ( "nlist numbers :\n" ) ;
    printoutivector ( alist.num_nlist , 1 , alist.N ) ; 
    printf ( "mlist numbers :\n" ) ;
    printoutivector ( alist.num_mlist , 1 , alist.M ) ; 
    if ( verbose >= 2  )   {
      printf ( "nlist :\n" ) ;
      for ( n = 1 ; n <= N ; n++ ) {
	printf ( "%3d : " , n  ) ;
	printoutivector ( alist.nlist[n] , 1 , alist.num_nlist[n] ) ; 
      }
      printf ( "\n" ) ; 
      printf ( "mlist :\n" ) ;
      for ( m = 1 ; m <= M ; m++ ) {
	printf ( "%3d : " , m  ) ;
	printoutivector ( alist.mlist[m] , 1 , alist.num_mlist[m] ) ; 
      }
      printf ( "\n" ) ; 
    }
  }

/* theoretical calculation of ft, using average number of 1's */
  tmpi = alist.tot / alist.M ;
  for ( ft = 0.0 ; tmpi > 0 ; tmpi -- )
    ft = fs * ( 1- ft ) + (1 - fs ) * ft ; 

  tmpd = err * (1.0 -ft) + (1-err) * ft ; 
  printf ( "fs = %f , ft = %f , (fte = %f ,) err = %f , fz = %f\n" , fs, ft, fte, err, tmpd) ; 
  printf ( "H(s) = %f, H(z;t) = %f = %f - %f ; capacity of channel 1 - H(n) = %f. " , 
	  h2(fs) , h2(tmpd) - h2(err) , h2(tmpd) , h2(err) , 1 - h2(err) ) ; 
  if ( h2(fs) > h2(tmpd) - h2(err) ) {
    printf ( "task impossible\n" ) ; 
    exit ( 0 ) ; 
  }
  else printf ( "\n" ) ; 


  if ( DEMO || verbose >= 2  ){
    t2 = cvector ( 1 , M ) ; 
    printf ( "applying list representation to string\n" ) ; 
    for ( m = 1 ; m <= M ; m++ ) {
      t2[m] = cdotprod_mod2_index ( alist.mlist[m] , s , 1 , alist.num_mlist[m] ) ; 
      printf ( "%d " , t2[m] ) ; 
    }
    printf ( "\nand the t from before again:\n" ) ; 
    printoutcvector ( t , 1 , N ) ; 
    printf ( "\n" ) ; 
  }