msb.c

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

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

   - 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. 

   msb.c is the version of ms.c designed for big sparse matrices A
   that are represented by double lists. 

   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. 

    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;
  double *m;
} 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 ) ;
/* making the data */ 

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


static int taps , tapsmin=2 , tapsmax=0 ; 
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 int **mlist, **nlist, *num_mlist, *num_nlist ,
           *l_up_to , *norder , *suspicion ; 
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 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 taps_approx = 0.0 ; 
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 */
static double betahuge = 10000 ; 
static double betap = 1.0 ; /* factor for the prior bit of the energy, */
static double betap_huge = 100.0 ; /* only enabled when betahuge is too */

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 = 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 ;
  thing *my ; 
  FILE *fp; 

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

/*  new things for ms */
  int biggest_num_m ,   biggest_num_n ;

/*  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("msb 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 {
    if( MS_DEMO || verbose >= 1 )     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 || taps > tapsmax ) ; 
  if( MS_DEMO || verbose >= 1 )   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" ) ; 
  }
  if( MS_DEMO || verbose >= 1 ) 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" ) ; 

  if( MS_DEMO || verbose >= 1 )   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" ) ; 

/* one idea for bounding the list: */
/*    biggest_num_nlist += MIN((taps+1), (M0 - two_to_p * N0)) ; */
  biggest_num_m = (taps+1) ; 
  biggest_num_n = (taps+1) * pmax ; 
  num_mlist = ivector ( 1 , M ) ; 
  mlist = imatrix ( 1 , M , 1 , biggest_num_m ) ; 
  num_nlist = ivector ( 1 , N ) ; 
  nlist = imatrix ( 1 , N , 1 , biggest_num_n ) ; 
  l_up_to = ivector ( 1 , M ) ; 


  if( MS_DEMO || verbose >= 1 )   printf ( "generating matrix in list representation\n" ) ; 
  for ( m = 1 ; m <= M ; m++ ) {
    num_mlist[m] = taps + 1 ; 
  }
  for ( n = 1 ; n <= N ; n++ ) {
    num_nlist[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 ++ ) {
      u = 1 ; 
      for ( l = n , i = 1 ; i <= N0 + 1 ; i ++ , l += two_to_p ) {
	if (  pol[i]       ) {
	  mlist[m][u] = l ;
	  u ++ ; 
	  num_nlist[l] ++ ;
	  nlist[l][ num_nlist[l] ] = m ; 
	}
      }
      if ( u != biggest_num_m + 1 ) printf ("eek %d %d\n",u,biggest_num_m+1);
      m++ ; 
    }
  }
  if ( MS_DEMO || verbose >= 2  )   {
    printf ( "nlist numbers :\n" ) ;
    printoutivector ( num_nlist , 1 , N ) ; 
    printf ( "nlist :\n" ) ;
    for ( n = 1 ; n <= N ; n++ ) {
      printf ( "%3d : " , n  ) ;
      printoutivector ( nlist[n] , 1 , num_nlist[n] ) ; 
    }
    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_index ( mlist[m] , t0 , 1 , num_mlist[m] ) ) ; 
    printf ( "\n" ) ; 
    printf ( "error checksum\n" ) ; 
  }
  for ( m = 1 ; m <= M ; m++ ) {
    t[m] = cdotprod_mod2_index ( mlist[m] , s , 1 , num_mlist[m] ) ; 
    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_index ( mlist[m] , z , 1 ,  num_mlist[m] ) ) ; 
    printf ( "\n" ) ; 

    printf ( "Error suspicion vector\n" ) ; 
  }
  suspicion = ivector ( 1 , N ) ; 
  for ( n = 1 ; n <= N ; n++ ) {
    suspicion[n] =  cdotprod_index ( nlist[n] , t , 1 , num_nlist[n] ) ;
    if ( MS_DEMO || verbose >= 2  ) printf ( "%d " , suspicion[n] ) ; 
  }
  if ( MS_DEMO || verbose >= 2  ) {
    printf ( "\nand the noise again:\n" ) ; 
    printoutcvector ( s , 1 , M0 ) ; 
    printf ( "\n" ) ; 
  }
/*  if ( MS_DEMO ) exit(0) ; */

  /* set up the norder */
  norder = ivector ( 1, N ) ; 
  if ( NORDER == 3 || NORDER == 2 || NORDER == 1 ) {
    my = (thing * ) malloc ( N * sizeof(thing)) ; 
    for ( n = 0 ; n <= N-1 ; n ++ ) {
      my[n].index = n+1 ;