msb.c

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

    i ++ ; 
    for ( n = 1 ; n <= N ; n++ ) {
      gr[n] = 0.0 ; 
    }

    find_q_S ( x , 1 ) ; /* zero if we don't want the entropy; 
			switch this to 1 to make a routine that makes both 
			the objective function and its derivative */
    E = 0.0 ; 

    m = 1 ;
    for ( ; m <= M ; m ++ ) { 
      E +=     forward_pass ( ) * g[m] ; 

      backward_pass ( gr ) ;  
    }

/*    sum_dx = 0.0 ; 
    ip_dx = 0.0 ;  */
    for ( n = 1 , gg = 0.0  ; n <= N ; n++ ) {
      gr[n] *= beta ;             /* temperature effect */
      gr[n] -= betap * bias[n] ; /* prior effect */
      tmpd = gr[n] ; 
      gr[n] += x[n] ;               /* derivative of entropy */
/*      dx[n] = - tmpd -  x[n] ;    MONITOR */
/*      sum_dx += dx[n] * dx[n] ;   MONITOR */ 
/*      ip_dx += dx[n] * dxp[n] ;   MONITOR */
/*      dxp[n] = dx[n] ;  MONITOR */
      x[n] = - tmpd ; 
      gg += gr[n] * gr[n] ; 
    }
/*    sdxp = sdx1 ;                 MONITOR */
/*    sdx1 = sqrt(sum_dx);          MONITOR */
/*    printf("%6.3g %6.3g : ",sdx1, (sdx1*sdxp > 0) ? ip_dx/(sdx1*sdxp) : 0.0 ); */
/*    pdv( x , 1 , 10 , 63 ) ;     MONITOR */
/*     if(!(i%1))printf("\n"); */
/*    fflush(stdout);   MONITOR */

    tmpd = S - beta * E ; 
    if ( verbose ) {
      print_state ( tmpd , gg , x ) ; 
    }
  }
  while ( gg > ftol && i < itmax ) ;
  if ( MS_DEMO ) printf ("%3d Jumps, beta %6.3g. ",i,beta ) ; 

  free_dvector ( gr , 1 , N ) ; 
/*  free_dvector ( dx , 1 , N ) ;  MONITOR */
/*  free_dvector ( dxp , 1 , N ) ;  MONITOR */
}


void seq_opt (double *x , param *controlp )
/* optimization by `reestimation', sequential */
/* an advantage of this approach is that there is a free energy
   improvement guarantee. Also we can expect some sort of metric
   effect-- i.e. unlike steepest descent, here, the change of x[n+1]
   will be influenced by the change just made in x[n] */
{
  double tmpd , grnn , v0 , v1 ; 
  int i = 0  , nn , l , ll ; 
  int u ; /* runs over the norder */

  v1 = objective ( x , controlp ) ; /* this sets up the qs for us too */
  v0 = v1 + 2 * ftol ; /* hack to prevent premature termination */

  do {
    mega_forward_pass ( ) ; 
    i ++ ; 
    for ( u = N ; u >= 1 ; u-- ) { /* reverse sequence used in mega_for */
      nn = norder[u] ; 
      grnn = 0.0 ; 
      for ( l = 1 ;  l <= num_nlist[nn] ; l ++ ) {
	m = nlist[nn][l] ; 
/*      for ( m = 1 ; m <= M ; m ++ ) { */
	ll = l_up_to[m] ;
	tmpd = mp1[m][ll-1] * mp1[m][ll+1] +  mp0[m][ll-1] * mp0[m][ll+1] -
	       mp0[m][ll-1] * mp1[m][ll+1] -  mp1[m][ll-1] * mp0[m][ll+1]  ;
	grnn -= g[m] * tmpd ; 
      }
      grnn *= beta ;
      grnn -= betap * bias[nn] ; /* prior effect */
      x[nn] = - grnn ; /* evaluate the new value of q */
      x_to_q ( x[nn] , &q0[nn] , &q1[nn] , 0 ) ;
      if ( seqverbose >= 2 ) printf("%4d %8.3g\n",nn,x[nn] ) ; 
/* print_state ( 0.0 , 0.0 , x ) ;  */
      mega_backward_step ( nn ) ;
    }
    if ( seqverbose >= 2 ) printf("\n");
    if ( !(i%seq_period) ) { /* default period = 10 */
/*      if ( seqverbose ) printf ("Seq %d %8.5g " , i , v1  ) ; */
      v0 = v1 ; 
      v1 = objective ( x , controlp ) ;
      if ( seqverbose ) printf ("%8.5g %8.5g\n" , v1 , (v0-v1) ) ; 
      if ( v0 - v1 < -ftol  ) fprintf ( stderr , "Warning, F increased significantly %6.3g\n" , v0 - v1 ) ; 
    }
  }
  while ( v0 - v1 > ftol && i < itmax ) ;
  if (seqverbose) 
    printf ("seq_opt done: %8.3g %3d Js, beta %6.3g\n",(v0-v1),i,beta ) ; 

}

void print_state ( double tmpd , double gg , double *x ) {
/*  printf ( "%6.3g = %6.3g - %6.3g: " , tmpd , S , E ) ; */
  for ( n = 1 ; n <= N ; n++ ) printf ( "%2.0f " , 99.0 * q1[n] ) ; 
  printf ( ":%6.3g\nx: " , gg ) ;
  for ( n = 1 ; n <= N ; n++ ) printf ( "%6.3g " , x[n] ) ; 
  printf ( "\n" ) ; 
}


void printall ( FILE *fp ) 
{
  printoutcvector ( s , 1 , N ) ; 
   printf ( "\n" ) ; 
  printoutcvector ( t , 1 , M ) ; 
  printoutcvector ( to , 1 , M ) ; 
  pdv ( g , 1 , M , 4 ) ;   printf ( "\n" ) ; 
}

/*          That's the end of the interesting code          */
/*     Here follow routines to parse the input command line */

#define DNT fprintf( fp, "\n        ")
#define NLNE  fprintf( fp, "\n")

static void print_usage ( char **argv , FILE * fp )
{
  fprintf( fp, "Usage: %s ",argv[0]);
  fprintf( fp, "datafile [optional arguments]");

  NLNE; fprintf( fp, " Arguments affecting verbosity of run-time output: <defaults>");
  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)"); 
  DNT;  fprintf( fp, "-SV | -SVV               (verbose in seq_opt)"); 
  DNT;  fprintf( fp, "-CG                      (check gradient) "); 
  DNT;  fprintf( fp, "-MS_DEMO                 (only do demo of M-S) "); 
  DNT;  fprintf( fp, "[-e epsilon]       <%g> (step size for gradient)", epsilon); 
  NLNE; fprintf( fp, " Data creation:" ) ; 
  DNT;  fprintf( fp, "-n0 N0           <%d> (number of bits in polynomial [=`k']) " , N0 ); 
  DNT;  fprintf( fp, "-m0 M0           <%d> (number of bits of data [=`N']) " , N0 ); 
  DNT;  fprintf( fp, "-err err           <%g> (error probability) " , err ); 
  DNT;  fprintf( fp, "-fpol fpol         <%g>  (fraction of bits high in poly)", fpol ); 
  DNT;  fprintf( fp, "-taps taps (approx)        (overrides fpol)" ); 
  DNT;  fprintf( fp, "-tapsmax tapsmax           " ); 
  DNT;  fprintf( fp, "-tapsmin tapsmin   <%d>   " , tapsmin ); 
  DNT;  fprintf( fp, "-fs fs             <%g>  (fraction in shift register seed)", fs); 
  NLNE; fprintf( fp, " Inference:" ) ; 
  DNT;  fprintf( fp, "-sdx sdx           <%3g> (sd of initial state)" , sdx ); 
  DNT;  fprintf( fp, "-b betastyle beta0 beta1 (what to do with beta)" ); 
  DNT;  fprintf( fp, "       betastyle 0: const; 1: linear; 2: multiply; 22: multiply and go wild on last loop" ); 
  DNT;  fprintf( fp, "-betap betap       <%3g> (prior's beta)" , betap ); 
  DNT;  fprintf( fp, "-opt optimizer     <%d>   (0=frp,1=macopt,2=jumpopt,3=seqopt,4=betaopt)",opt); 
  DNT;  fprintf( fp, "-itmax itmax     <%d>   (num of linmins or jumps per optimization)",itmax); 
  DNT;  fprintf( fp, "-norder n        <%d>   (for seqopt - 0/1/2 = 1..N / random / suspicion / suspicion/num )",NORDER); 
  DNT;  fprintf( fp, "-nl loops          <%d>   (number of cg runs)",NL); 
  DNT;  fprintf( fp, "-ftol ftol         <%g> ", ftol ); 
  DNT;  fprintf( fp, "-seed seed "); 
  DNT;  fprintf( fp, "-seqp seq_period   <%d>  (how many loops per check of delta F)" , seq_period ); 
  NLNE; fprintf( fp, " Files:");
  DNT;  fprintf( fp, "-o outfile               "); 
  fprintf( fp, "\n");
  return ;
}
#undef DNT
#undef NLNE

static void make_sense ( void ) 
{ /* routine to correct silly control parameters */
  if ( M < N ) fprintf ( stderr , "Warning, M < N\n" ) ;
  if ( betastyle == 1 && NL > 1 ) 	
    betaf = ( beta1 - beta ) / ((double) (NL - 1) ) ; 
  else if ( ( betastyle == 2 || betastyle == 22 ) && NL > 1 )
    betaf = exp ( log ( beta1 / beta ) / ((double) (NL - 1) )) ; 
/*  if ( opt > 4 ) opt = 4 ; */
  if ( taps_approx > 0.0 ) {
    fpol = ( taps_approx - 1.0 ) / ( ( double ) ( N0 - 1 ) ) ; 
  }
  if ( tapsmax == 0 ) { tapsmax = N0 ; }
}

static int process_command ( int argc , char **argv )
{
  int p_usage = 0 ;
  int status = 0 ;
  int cs , i ;

  if ( argc < 1 )     {
    p_usage = 1 ; 
    status -- ;
  }
#define ERROR1 fprintf ( stderr , "arg to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERROR2 fprintf ( stderr , "args to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERRORREG fprintf ( stderr , "regtype must be defined before `%s'\n" , argv[i] ) ; \
               status --
  for (i = 1 ; i < argc; i++)    {
    cs = 1 ;
    if ( strcmp (argv[i], "-V") == 0 )        {
      verbose = 1;
    }
    else  if ( strcmp (argv[i], "-VV") == 0 )        {
      verbose = 2;
    }
    else  if ( strcmp (argv[i], "-SV") == 0 )        {
      seqverbose = 1;
    }
    else  if ( strcmp (argv[i], "-SVV") == 0 )        {
      seqverbose = 2;
    }
    else  if ( strcmp (argv[i], "-CG") == 0 )        {
      CG = 1;
    }
    else  if ( strcmp (argv[i], "-NL") == 0 )        {
      NL = 0; 
    }
    else  if ( strcmp (argv[i], "-MS_DEMO") == 0 )        {
      MS_DEMO = 1; 
    }
    else  if ( strcmp (argv[i], "-nl") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &NL); 
    }
    else  if ( strcmp (argv[i], "-opt") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &opt); 
    }
    else  if ( strcmp (argv[i], "-norder") == 0 ) {
      if ( i + 1 == argc ) { ERROR1; }
      else cs *= sscanf(argv[++i], "%d", &NORDER ) ; 
    }
    else  if ( strcmp (argv[i], "-itmax") == 0 ) {
      if ( i + 1 == argc ) { ERROR1; }
      else cs *= sscanf(argv[++i], "%d", &itmax ) ; 
    }
    else   if ( strcmp (argv[i], "-seed") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%ld", &seed); 
    }  
    else  if ( strcmp (argv[i], "-seqp") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &seq_period); 
    }
/*    else   if ( strcmp (argv[i], "-m") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else {
	cs *= sscanf(argv[++i], "%d", &M); 
      }
    } */ /* NB M is defined by M0 */
    else   if ( strcmp (argv[i], "-n0") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else {
	cs *= sscanf(argv[++i], "%d", &N0); 
      }
    }
    else   if ( strcmp (argv[i], "-m0") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else {
	cs *= sscanf(argv[++i], "%d", &M0); 
      }
    }
/*    else   if ( strcmp (argv[i], "-n") == 0 )        {
      if ( i + 1 == argc )             { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &N); 
    }*/ /* NB N is defined to be M0 */
    else if ( strcmp (argv[i], "-e") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &epsilon); 
    }
    else if ( strcmp (argv[i], "-err") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &err ); 
    }
    else if ( strcmp (argv[i], "-sdx") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &sdx); 
    }
    else if ( strcmp (argv[i], "-betap") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &betap); 
    }
    else if ( strcmp (argv[i], "-b") == 0 ) {
      if ( i + 3 >= argc ) { ERROR2;      }
      else {
	cs *= sscanf(argv[++i], "%d", &betastyle); 
	cs *= sscanf(argv[++i], "%lf", &beta); 
	cs *= sscanf(argv[++i], "%lf", &beta1); 
      }
    }
    else if ( strcmp (argv[i], "-ftol") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &ftol); 
    }
    else if ( strcmp (argv[i], "-fs") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &fs); 
    }
    else if ( strcmp (argv[i], "-fpol") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &fpol); 
    }
    else if ( strcmp (argv[i], "-taps") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &taps_approx); 
    }
    else if ( strcmp (argv[i], "-tapsmax") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &tapsmax); 
    }
    else if ( strcmp (argv[i], "-tapsmin") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &tapsmin ) ; 
    }
    else if ( strcmp (argv[i], "-pheading") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &pheading); 
    }
    else  if ( strcmp (argv[i], "-o") == 0 )	{
      if ( i + 1 == argc ) 	    { ERROR1;      }
      else {
	strcpy(outfile, argv[++i]);
	if ( !printout ) printout = 1 ;
      }
    }
    else {
      fprintf ( stderr , "arg `%s' not recognised\n" , argv[i] ) ; 
      p_usage = 1 ;
      status -- ;
    }
    if ( cs == 0 ) {
      fprintf ( stderr , "arg at or before `%s' has incorrect format\n" , 
	       argv[i] ) ;
      p_usage = 1 ;
      status -- ;
    }
  }
  if ( p_usage ) print_usage ( argv , stderr ) ;
  return ( status ) ;
}
#undef ERROR1
#undef ERROR2
#undef ERRORREG

unsigned char cdotprod_mod2_index 
( 
 int *a , unsigned char *b , int lo , int hi ) {
  unsigned char ans = 0 ; 
  int i ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans ^= b[a[i]] ;
  }
  return ans ;
}

int cdotprod_index 
( 
 int *a , unsigned char *b , int lo , int hi ) {
  int ans = 0 ; 
  int i ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans += b[a[i]] ;
  }
  return ans ;
}


int cmp_thing ( const void *a , const void *b ) {
  const thing *ma = (const thing*)a ;
  const thing *mb = (const thing*)b ;
  return ( ( ma->val - mb->val ) ) ; 
}