ms.c

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

/* 
   ms.c                                      (c) DJCM 94 07 04
                                                      94 10 26

   - free energy minimization for inferring state of shift register - 

   This program is based on cg.c but uses the Meier-Staffelbach 
   approach to infer the noise rather than the sequence. 

   Data generation:

   (1) generate a random polynomial with number of taps: taps

   (2) generate random s vector (this is the old `s', not the 
       new s of this algm.      this s is called s0 and runs :1-N0 )
       and make the sequence t0= A0 s0 ( 1 - M0, here `N` = M0 ). 
       A0 is a virtual matrix, not evaluated. 

   (3) generate random e vector: 1 - M0, and make z: 1 - M0;
       z = t0 + e

   (4) generate matrix A (a.k.a. Omega) of low weight parity checks. 
       -- each based on (N0+1) bits. 
       There are (M0-N0) straight checks, 
       M0-2N0 of the next sort, 
       M0-4N0 of the next, 
       etc. until M0 < 2^p N0

   (5) Thus A has N = M0 and M = sum M0 - 2^p N0.

   Inference:

   (1) The task is now to solve A (e+z) = 0, 
       i.e. A e = A z
       This right hand side is a long list of parity checks. 
       The approach is to pretend that this problem has the form 
            A e = t + n, 
       and obtain P(e|t), then take the n-noise level to zero. 
       (i.e. the annealing goes from high temperature to temperature 
       zero, rather than 1.)

       Note it should be the case that A z is independent of s0, right?
       It's just a signature of the e. 

          You get to choose (on command line): 
	  N (state length),                                      (e.g. 100)
	  M (number of observations),                            (e.g. 1000)
	  err (probability of error in observation of t = As)    (e.g. 0.3)
	  seed (for random numbers)

    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')
	       
    28 oct 94  Modifying ms so it tries to detect oscillations in 
               jumpopt, and uses a sludge parameter to slow them away. 

*/

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


typedef struct {
  int n;
  double *m;
} param;

static int    process_command ( int , char ** ) ;
static void jump_opt (double * ) ; /* optimization by `reestimation' */
static void beta_opt (double * ) ; /* optimization by beta juggling */
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 void 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 ** ) ;

/* making the data */ 

static int MS_DEMO = 0 ; 
static int M0=30, N0=7; 


static int taps , tapsmin=2 ; 
static int n , N ; /* runs over state bits */
static int m , M ; /* runs over the data */
static unsigned char **A ;   /* binary matrix A[m][n] */
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 unsigned char *to ;   /* the observed binary vector 
		      to[m] = t[m] with prob 1-er[m] */
static double *g ;  /* log odds g[m] = +/- log er[m]/(1-er[m]) depending
		       whether to[m] = 1 or 0 . */
static double fpol = 0.1 ;
static double true_err ;   /* actual fraction */
static double fs = 0.5 ;  /* fraction of bits in s that are high */ 
static double err = 0.2 ;       /* default error probability */

/* 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 = 0.0 ; /* if annealing is used then 
			       beta goes from current value of 
			       beta to beta1 by uniform steps 
			       either arithmetic or geometric */

/* params for beta_opt juggling */
/* static double beta_dec1 = 0.71 ; */
static double beta_dec2 = 0.5 ;
static double beta_inc = 1.0712 ;
static double stepmax = 40.0 ; 
static double stepmin = 5.0 ;
static double sludge = 0.5 ;
/* static double ip_close = 0.5 ; */

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

static  char 	outfile[50] = "_out" ;
static  int     printout = 0 ;

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

/* 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; 

void main(int argc, char *argv[])
{
  param control;  /*  ignore this, it's for future structured program */
  double *x , v ;
  int count , count_high , count_viol , count_err ;
  FILE *fp; 

/* new things for ms */
  unsigned char *pol , *t0 , *z , **At ;
  int start , sum , i , j , l  ;
  int two_to_p , pmax, p ;

/*  char 	junk[100] ; */

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

  /* Given the settings of M0 and N0 plan how big things will be. */

  for ( M = 0 , pmax = 0 , two_to_p = 1 ; 
       two_to_p * N0 < M0 ; 
       two_to_p *= 2 ) {
    pmax ++ ;
    M += M0 - two_to_p * N0 ; 
  }
  N = M0 ; 

  printf("ms N = M0 = %d, N0 = %d, M = %d, pmax = %d, fpol = %6.3g, err = %6.3g, opt = %d, bstyle = %d\n",
	 M0 , N0 , M , pmax , fpol , err , opt , betastyle ) ;

  /* Set up memory , etc. */

  ran_seed ( seed ) ; 
/*  s0  = ivector ( 1 , N0 ) ;  not needed because the first N0 states
  are copied */
  t0  = cvector ( 1 , M0 ) ; /* THis is the noise-free emission */
  pol = cvector ( 1 , N0 + 1 ) ; 

  do {
    printf ( "gen'g taps " ) ; 
    pol[1] = 1 ; /* the  first bit must be 1 */
    taps = 1 ; 
    for ( n = 2 ; n <= N0 ; n ++ ) {
      taps += ( pol[n] = ( ranu() < fpol ) ? 1 : 0 ) ; 
    }
  } while ( taps < tapsmin ) ; 
  printf ( " - %d\n" , taps );
  pol[n] = 1 ; /* the N0 + 1 entry is 1 for convenience in checksum defn */
  if ( MS_DEMO || verbose >= 2 ) printoutcvector ( pol , 1 , N0 ) ; 
  printf ( "\n" ) ; 

  printf ( "generating sequence\n" ) ;
  for ( n = 1 ; n <= N0 ; n++ ) 
    t0[n] = ( ranu() < fs ) ? 1 : 0 ; 
  for ( start = 1 ; n <= M0 ; n++ , start ++ ) {
    for ( sum = 0 , i = start , j = 1 ; i < n ; i ++ , j ++ ) 
      sum ^= ( t0[i] & pol[j] ) ; 
    t0[n] = sum ; 
  }
  if ( MS_DEMO  || verbose >= 2 ) printoutcvector ( t0 , 1 , M0 ) ; 
  if ( MS_DEMO || verbose >= 2  ) printf ( "\n" ) ; 

  printf ( "generating noise\n" ) ;
  s  = cvector ( 1 , N ) ; /* this is going to be the error vector */
  so = cvector ( 1 , N ) ; 
  z = cvector ( 1 , N ) ; /* this'll be the observation */
  count_err = 0 ; 
  for ( n = 1 ; n <= N ; n++ ) {
    count_err += (     s[n] = ( ranu() < err ) ? 1 : 0  ) ; 
    so[n] = s[n] ;
    z[n] = s[n] ^ t0[n] ;
  }
  if ( MS_DEMO || verbose >= 2  )printoutcvector ( s , 1 , M0 ) ; 
  if ( MS_DEMO || verbose >= 2  )printf ( "\n" ) ; 

  printf ( "generating matrix " ) ; 
  printf ( "\n" ) ; 
  A  = cmatrix ( 1 , M , 1 , N ) ; 
  At  = cmatrix ( 1 , N , 1 , M ) ; 
  for ( m = 1 ; m <= M ; m++ ) {
    for ( n = 1 ; n <= N ; n++ ) {
      A[m][n] = 0 ;
    }
  }
  for ( m = 1 , p = 0 , two_to_p = 1 ; 
       p <= pmax ; 
       p ++ , two_to_p *= 2 ) {
    for ( n = 1 ; n <= M0 - two_to_p * N0 ; n ++ ) {
      for ( l = n , i = 1 ; i <= N0 + 1 ; i ++ , l += two_to_p ) {
	A[m][l] = pol[i] ;
      }
      m++ ; 
    }
  }
  for ( m = 1 ; m <= M ; m++ ) {
    for ( n = 1 ; n <= N ; n++ ) {
      At[n][m] = A[m][n] ;
    }
  }
  if ( MS_DEMO || verbose >= 2 ) {
    if ( verbose >= 1 )
      printoutcmatrix ( A , 1 , M , 1 , N ) ;  
    if ( verbose >= 2 ) 
      printoutcmatrix ( At , 1 , N , 1 , M ) ;  
  }
  if ( MS_DEMO || verbose >= 2  ) printf ( "\n" ) ; 

/*  Checksum time */

  t = cvector ( 1 , M ) ; 
  to = cvector ( 1 , M ) ; 

  if ( MS_DEMO || verbose >= 2  ){
    printf ( "t0 checksum\n" ) ; 
    for ( m = 1 ; m <= M ; m++ ) 
      printf ( "%d " , cdotprod_mod2 ( A[m] , t0 , 1 , N ) ) ; 
    printf ( "\n" ) ; 
    printf ( "error checksum\n" ) ; 
  }
  for ( m = 1 ; m <= M ; m++ ) {
    t[m] = cdotprod_mod2 ( A[m] , s , 1 , N ) ; 
    if ( MS_DEMO || verbose >= 2  ) printf ( "%d " , t[m] ) ; 
/* because this is in fact noise free, the observed to is the same as t */
    to[m] = t[m] ; /* this is the last time t is used in this role */
  }
  if ( MS_DEMO || verbose >= 2  ) {
    printf ( "\n" ) ; 
    printf ( "z checksum\n" ) ; 
    for ( m = 1 ; m <= M ; m++ ) 
      printf ( "%d " , cdotprod_mod2 ( A[m] , z , 1 , N ) ) ; 
    printf ( "\n" ) ; 

    printf ( "Error suspicion vector\n" ) ; 
    for ( n = 1 ; n <= N ; n++ ) 
      printf ( "%d " , cdotprod ( At[n] , t , 1 , M ) ) ; 
    printf ( "\nand the noise again:\n" ) ; 
    printoutcvector ( s , 1 , M0 ) ; 
    printf ( "\n" ) ; 
  }
/*  if ( MS_DEMO ) exit(0) ; */

    


/* old program from here */

/*  er = dvector ( 1 , M ) ;  */
  g  = dvector ( 1 , M ) ; 
  x  = dvector ( 1 , N ) ; 
  q0 = dvector ( 1 , N ) ; 
  q1 = dvector ( 1 , N ) ; 
  if ( opt == 3 ) { /* mega state vectors needed */
    mp0 = dmatrix ( 1 , M , 0 , N + 1 ) ;
    mp1 = dmatrix ( 1 , M , 0 , N + 1 ) ;
    for ( m = 1 ; m <= M ; m ++ ) {
      mp0[m][0] = 1.0 ; mp1[m][0] = 0.0 ; 
      mp0[m][N+1] = 1.0 ; mp1[m][N+1] = 0.0 ; 
    }
  }
  p0 = dvector ( 0 , N + 1 ) ;