bnd_iterate.c

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

/* 
   bnd_iterate.c

   iterates the bn decoding rule, inventing random bits as it goes 
   along. 

   at each step 100 selections of 21 states
   are made from the last list of states, and are propagated through 
   100 random x bits and the corresponding z bits to give 
   100 new states. 



*/

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

typedef struct {
  unsigned char *x ; /* x[1..N] */
  double *d ;        /*  probability diff contributed */
} bndi_state ;

typedef struct {
/*  unsigned char **x ;  x[rel][friend] */
  double *d ;         /* d[rel] product of probability diff contributed */
  unsigned char *z ;  /* z[rel] */
  unsigned char x0  ; /* state at root */
} bndi_tree ;

typedef struct { 
#include "bndi_var_str.h"
} bndi_control ;

static void print_usage ( char ** , FILE *  , bndi_control * ) ;
static void c_defaults ( bndi_control * ) ; 
static int process_command ( int  , char ** , bndi_control * ) ;
static int make_sense ( bndi_control * ) ;
static int allocate ( bndi_tree * , bndi_state * , bndi_state * , bndi_control * ) ;
static void initialize_state ( bndi_state * , bndi_control * ) ;
static void write_state ( FILE * , bndi_state * , bndi_control * ) ;
static void write_e ( FILE * , bndi_state * , bndi_control * ) ;
static void initialize_tree (  bndi_tree * , bndi_state * , bndi_control * ) ;
static void vertical_pass ( bndi_tree * , bndi_state * , bndi_control * , int ) ;
static double h2 ( double ) ;

void   main ( int , char ** ) ;

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

  bndi_state sa , sb , *s1 = &sa , *s2 = &sb , *s0 ;
  bndi_tree tree ;
  bndi_control c ;

  c_defaults ( &c ) ; 

  if ( process_command ( argc, argv, &c ) < 0 ) exit (0) ;
  if ( make_sense ( &c ) < 0 ) exit (0) ;
  
  if ( allocate ( &tree , &sb , &sa , &c ) < 0 ) exit (0) ; 

  initialize_state ( s1 , &c ) ;

  if ( c.write_entropy	) {
    sprintf ( c.efile , "bndi/H.%d.%d.%5.3f.%d" , 
	     c.t , c.tr , c.f , c.N ) ; 
    fp = fopen ( c.efile , "w" ) ; 
    if ( !fp ) {
      fprintf ( stderr , "can't open %s\n" , c.efile ) ;
      c.write_entropy = 0 ; 
    } else {
      fclose ( fp ) ; 
    }
  }

  for ( c.loop = 1 ; c.loop <= c.loops && ( c.h_non_zero ) ; c.loop ++ ) {

    for ( n = c.N ; n >=1 ; n -- ) {
      initialize_tree ( &tree , s1 , &c ) ; 
      vertical_pass ( &tree , s2 , &c , n ) ;    /* answer is dropped in s2[n] */
    }
    s0 = s1 ; s1 = s2 ; s2 = s0 ; 

    /* printout */

    if ( c.write_state && 
	( ( ( c.loop <= c.write_enda ) && ( !(c.loop% c.write_period_a )) ) ||
	 ( ( c.loop > c.write_enda ) && !(c.loop % c.write_period_b ) ) )
	) {
      sprintf ( c.statefile , "bndi/%d.%d.%5.3f.%d" , c.t , c.tr , c.f , c.loop ) ; 
      fp = fopen ( c.statefile , "w" ) ; 
      if ( !fp ) {
	fprintf ( stderr , "can't open %s\n" , c.statefile ) ;
      } else {
	write_state ( fp , s1 , &c ) ; 
	fclose ( fp ) ; 
      }
    }
    if ( c.write_entropy	) {
      fp = fopen ( c.efile , "a" ) ; 
      if ( !fp ) {
	fprintf ( stderr , "can't open %s\n" , c.efile ) ;
      } else {
	write_e ( fp , s1 , &c ) ; 
	fclose ( fp ) ; 
      }
    }
    
  }
}

static void initialize_state ( bndi_state *s , bndi_control *c ) {
  int n ; 
  double d = 1.0 - 2.0 * c->f ; /* initial difference */

  for ( n = c->N ; n >= 1 ; n -- ) {
    s->x[n] = ( ranf() < c->f ) ? 1 : 0 ;
    s->d[n] = d ;
  }
}

static void write_state ( FILE *fp , bndi_state *s , bndi_control *c ) {
  int n ; 

  for ( n = c->N ; n >= 1 ; n -- ) {
    fprintf ( fp , "%d %10.8f\n" ,     s->x[n] ,    s->d[n] ) ;
  }
}

static void write_e ( FILE *fp , bndi_state *s , bndi_control *c ) {
  int n ; 

  double h = 0.0 ; 

  for ( n = c->N ; n >= 1 ; n -- ) {
    h += h2 ( 0.5 * ( 1.0 + s->d[n] ) ) ; 
  }
  h /= (double)(c->N) ;
  fprintf ( fp , "%d %10.8f %d %d %10.8f %d\n" ,     c->loop ,    h , c->t , c->tr , c->f , c->N ) ;
  c->h = h ; 
  if ( h <= c->hlim )  c->h_non_zero = 0 ; 
}

static double h2 ( double x ) {
  double tmp ; 
  if ( x<= 0.0 || x>= 1.0 ) { tmp = 0.0 ; }
  else {
    tmp = x * log ( x ) + ( 1.0 - x ) * log ( 1.0 - x ) ;
  }
  return - tmp / log ( 2.0) ; 
}

static void initialize_tree (  bndi_tree *t , bndi_state *s , bndi_control *c ) {
  int n ; 
  int r , f ; 
  unsigned char *z = t->z ; 
  double *d = t->d ; 

  t->x0 = ( ranf() < c->f ) ? 1 : 0 ;

  for ( r = c->up_rels ; r >= 1 ; r-- ) {
    z[r] = t->x0 ; 
    d[r] = 0.5 ; 
    for ( f = c->friends ; f >= 1 ; f-- ) {
      n = rani( c->N ) + 1 ;  /* pick a state from the last level */
      z[r] ^= s->x[n] ;       /* don't need t->x[r][f] = */
      d[r] *= s->d[n] ; 
    }
    if ( z[r] ) d[r] = - d[r] ; 
  }
}

static void vertical_pass ( bndi_tree *t , bndi_state *s , bndi_control *c , int n ) {
  double q1 = c->f , q0 = 1.0 - c->f ; 
  double pc0 , pc1 , dpc , sum , dif ; 
  int r ;

  for ( r = c->up_rels ; r >= 1 ; r-- ) {
    dpc = t->d[r]  ; 
    pc0 = 0.5 + dpc ;
    pc1 = 0.5 - dpc ;
    
    q0 *= pc0 ;
    q1 *= pc1 ;
  }
  sum = q0 + q1 ;
  dif = q0 - q1 ;
  s->d[n] = dif / sum ;
  s->x[n] = t->x0 ; 
}

static int make_sense (  bndi_control *c ) {
  int status = 0 ; 
  c->up_rels = c->t - 1 ;
  c->friends = c->tr - 1 ;
  c->tot_friends = c->up_rels * c->friends ; 
  ran_seed ( c->seed ) ; 
  /* check for blatant attempt to go above shannon limit */
  if ( c->shannon_check &&
      ( ( (double)(c->tr) * h2(c->f) / (double)(c->t) ) > 1.0 ) ) {
    fprintf (stderr , "above shannon limit; exiting..\n" ) ; 
    status -- ; 
  }
  return status ; 
}

static int allocate ( bndi_tree *t , bndi_state *sb , bndi_state *sa , bndi_control *c ) {
  int status = 0 ; 
  int N = c->N ; 
  sa->x = cvector ( 1 , N ) ; 
  sb->x = cvector ( 1 , N ) ; 
  sa->d = dvector ( 1 , N ) ; 
  sb->d = dvector ( 1 , N ) ; 
/*  t->x = cmatrix ( 1 , c->up_rels , 1 , c->friends ) ;  */
  t->d = dvector ( 1 , c->up_rels ) ; 
  t->z = cvector ( 1 , c->up_rels ) ; 
  return status ;
}

static void c_defaults ( bndi_control *c ) 
{
#include "bndi_var_def.c"
}

static int process_command ( int argc , char **argv , bndi_control *c ) {
  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 )        {
      c->verbose = 1;
    }
#include "bndi_var_clr.c"
    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 , c ) ;
  return ( status ) ;
}
#undef ERROR1
#undef ERROR2
#undef ERRORREG

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

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

#include "bndi_var_usg.c"
  NLNE ;
  return ;
}

/*
<!-- hhmts start -->
Last modified: Mon Sep 11 14:28:32 1995
<!-- hhmts end -->
*/