mncn.c

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

  p->xo = vec->xo ; 
  return status ; 
}

int check_alist_MN ( alist_matrix *a , mnc_vectors *v ) {
  int status = 0 ; 
  if (     v->N != a->N 
      ||   v->M != a->M      ) {
    fprintf ( stderr , "eek %d %d %d %d\n" , v->N , a->N , v->M , a->M ) ; 
    status -- ; 
  }
  return status ; 
}


static int score ( data_creation_param *dc , 
       		   mnc_vectors  *v , mnc_all *all 
) /* compares so with s and works out quality of solution if wrong. */
{
  int status = 0 ; 
  int n , count = 0 , m , M = v->M , *nsum = v->nsum ; 
  FILE *fp ; 

  if ( dc->verbose ) printf ( "Scoring answer --- " ) ; 
  v->count = 0 ; v->count_high = 0 ; v->viols = 0  ;
  for ( n = 1 ; n <= v->N ; n++ ) { /* this compares all bits */
    if ( v->xo[n] != v->x[n] ) {    
      v->count ++ ; 
    }
    v->count_high += v->xo[n] ;
  }
  if ( v->count > 0 ) {
    alist_times_cvector_mod2 ( (all->bndp->a) , v->xo , v->t ) ; 

    for ( m = 1 ; m <= M ; m ++ ) {
      if ( v->t[m] != v->z[m] ) v->viols ++ ; 
    }

    if ( dc->verbose ) printf ( "fail\n" ) ; 
/*    else fprintf ( stderr , " X" ) ;  */
    dc->failcount ++ ; 
    if ( v->count ) {
      if ( dc->maxcount != 0 && dc->error_log ) {
	fp = fopen ( dc->error_logfile , "a" ) ;
	if ( !fp ) {
	  fprintf ( stderr , " couldn't open logfile %s\n" , dc->error_logfile ) ; 
	  dc->error_log = 0 ; 
	} else {
	  /* tell me more about this failure */ 	    /* write to a file */
	  finalline ( fp , all , dc ) ;
	  fprintf ( fp , "# " ) ;
	  for ( n = 1 ; n <= v->N && ( ( count <= dc->maxcount ) 
				      || dc->maxcount < 0 ) ; n++ ) { 
	    if ( v->xo[n] != v->x[n] ) {    
	      count ++ ;
	      if ( v->xo[n] ) {
		fprintf ( fp , "+%4d " , n ) ; /* + indicates extra bit given in so */
	      } else {
		fprintf ( fp , "-%4d " , n ) ; /* - indicates bit lacking     */
	      }
	    }
	  }
	  fprintf ( fp , "\n" ) ;
	  fclose (fp ) ; 
	}
      }
      if ( dc->noise_sum ) { 
/* Keep track of the noise vectors that caused the errors (mean only) */
	fp = fopen ( dc->noise_sumfile , "w" ) ;
	if ( !fp ) {
	  fprintf ( stderr , " couldn't open logfile %s\n" , 
		   dc->noise_sumfile ) ; 
	  dc->noise_sum = 0 ; 
	} else { /* pop the total number of failures and number of bits wrong
		     at the top in nsum[0], nsum[-1] */
	  nsum[0] ++ ;
	  for ( n = 1 ; n <= v->N ; n++ ) { 
	    if ( v->x[n] ) {    
	      nsum[n] ++ ; nsum[-1] ++ ;
	    }
	  }
	  fprintf ( fp , "#noise %d blocks %d bits %d block_error_rate %d / %d \n" , nsum[0] , nsum[-1] , v->N  , dc->failcount, dc->message ) ;
	  for ( n = 1 ; n <= v->N ; n++ ) { 
	    fprintf ( fp , "%d" , nsum[n] ) ;
	    fprintf ( fp , "\n" ) ;
	  }
	  fclose (fp ) ; 
	}
      }
    }
  }  else {
    if ( dc->verbose ) printf ( "SUCCESS\n" ) ; 
  }
  /* could also work out the typicality of the solution */
  return status ; 
}



static int make_norder ( alist_matrix *alist , fe_min_control *c )
{ /*
     Once this is changed to include other orders, 
     then the command must be called every time there is a 
     switch of c (or some other solution to this problem found)
*/
  int status = 0 ; 
  int n , N=alist->N ; 

    /* set up the norder */
  switch ( c->NORDER ) {
  case(0):
    
    for ( n = 1 ; n <= N ; n++ ) {
      alist->norder[n] = n ; 
    }
    break;
  default:
    status -- ;
    fprintf ( stderr , "invalid NORDER %d \n" , c->NORDER ) ; 
    break ; 
  }

  return status ; 
}

/* needs personalising for fe or depersonalising for bnd */
int make_vectors_quick ( data_creation_param *dc , mnc_vectors *v , mnc_all *all )
{
  int status = 0 ; 
  int count_s = 0 , count_n = 0 ; 

  if ( dc->nsn ) { /* a fixed weight vector for all s and n is required */
    fixed_wt_cvector ( v->x , dc->nsn , v->nsfrom , v->nnto ) ;
    v->count_s = dc->nsn ; 
  } else { /* other data generation rules */
    if ( !dc->ns ) { /* in NMN codes there are source bits to make */
      count_s = random_cvector ( v->x , dc->fs , v->nsfrom , v->nsto ) ; 
    } else {
      fixed_wt_cvector (  v->x , dc->ns , v->nsfrom , v->nsto ) ;
      count_s = dc->ns ; 
    }
    
    if ( dc->fgc ) { /* fixed noise level first, assumed noise level second */
      count_n = make_fgaussian_noise_bits_and_fix_biases 
	( v->x , v->bias , v->n , dc->fgcx, dc->fgcnf , dc->gcx , v->nnfrom , v->nnto ) ;
    }
    else if ( dc->gc ) { /* Here is assumed that using bnd */
      count_n = make_gaussian_noise_bits_and_fix_biases 
	( v->x , v->bias , dc->gcx , v->nnfrom , v->nnto ) ;
    }
    else {
      if ( !dc->nn ) {
	count_n = random_cvector ( v->x , dc->fn , v->nnfrom , v->nnto ) ; 
      } else {
	fixed_wt_cvector (  v->x , dc->nn , v->nnfrom , v->nnto ) ;
	count_n = dc->nn ; 
      }
    }
    v->count_s = count_s + count_n ;
  }
    
  dc->true_fs = 0.0 ;
  dc->true_fn = (double) count_n / (double) dc->M ; 
    
  alist_times_cvector_mod2 ( all->a , v->x , v->z ) ; 

  return status ; 
}
int make_vectors_quickER ( data_creation_param *dc , mnc_vectors *v , mnc_all *all ) /* just makes the calls to the rand number generator */
{
  int status = 0 ; 
  int count_s = 0 , count_n = 0 ; 

  if ( !dc->ns ) {
    count_s = random_cvector ( v->x , dc->fs , v->nsfrom , v->nsto ) ; 
  } else {
    fixed_wt_cvector (  v->x , dc->ns , v->nsfrom , v->nsto ) ;
  }

  if ( dc->gc ) { /* Here is assumed that using bnd */
    count_n = make_gaussian_noise_bits_ONLY 
      ( v->x , v->bias , dc->gcx , v->nnfrom , v->nnto ) ;
  }
  else {
    if ( !dc->nn ) {
      count_n = random_cvector ( v->x , dc->fn , v->nnfrom , v->nnto ) ; 
    } else {
      fixed_wt_cvector (  v->x , dc->nn , v->nnfrom , v->nnto ) ;
    }
  }

  return status ; 
}

static void process_vector_fe ( mnc_vectors *v , fe_min_param *p ) {
  int m ; 
  p->true_s = 1 ; 
  for ( m = 1 ; m <= p->M ; m++ ) {
    p->g[m] = ( v->z[m] ) ? 1.0 : -1.0 ; 
  }
}

static int make_gaussian_noise_bits_and_fix_biases ( unsigned char *x ,
						    double *b,
						    double gcx ,
						    int lo , int hi ) {
  int i , c = 0 ;
  double z , p ;

  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    z = gcx + rann() ; 
    p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcx ) ) ; 
    if ( z <= 0.0 ) { /* a noise bit is high */
      x[i] = 1 ; c ++ ;
      b[i] = p ; 
      if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %f\n" , b[i] ) ;
    } else {
      x[i] = 0 ; 
      b[i] = 1.0 - p ; 
      if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %d %f\n" , x[i] , b[i] ) ;
    }
  }
  return c ; 
}

/* fixed noise - by making n a unit-like vector */ 
static int make_fgaussian_noise_bits_and_fix_biases ( unsigned char *x ,
						    double *b,double *n,
						    double fgcx ,
						    double fgcnf ,
						    double gcx ,
						    int lo , int hi ) {
  int i , c = 0 ;
  double z , p , power=0.0 , factor , count = 0.0 ;

  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    n[i] = rann() ; 
    power += n[i] * n[i] ; 
    count += 1.0 ; 
  }
  factor = sqrt( count / power ) ; 
  fprintf ( stdout , "f=%8.6g -> %5g : " , factor , fgcnf ) ; 
  fflush(stdout) ; 
  factor *= fgcnf ; /* boost or reduce noise for importance sampling */
  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    z = fgcx + n[i] * factor ; /* fixed signal to noise level */
    p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcx ) ) ; /* gcx - assumed value */ 
    if ( z <= 0.0 ) { /* a noise bit is high */
      x[i] = 1 ; c ++ ;
      b[i] = p ; 
      if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %f\n" , b[i] ) ;
    } else {
      x[i] = 0 ; 
      b[i] = 1.0 - p ; 
      if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %d %f\n" , x[i] , b[i] ) ;
    }
  }
  return c ; 
}

static int make_gaussian_noise_bits_ONLY ( unsigned char *x ,
						    double *b,
						    double gcx ,
						    int lo , int hi ) {
  int i , c = 0 ;
  double z ;

  for ( i = lo ; i <= hi ; i ++ ) {
    /* make a random normal variate with s.d. 1.0 and mean gcx */
    z = rann() ; 
  }
  return c ; 
}

static void set_up_biases ( double *b , mnc_vectors *v , data_creation_param *dc ) {
  int n ;
  double lfn1fn = log ( dc->fn / ( 1 - dc->fn ) ) ;
  double lfs1fs = log ( dc->fs / ( 1 - dc->fs ) ) ;

  for ( n = v->nsfrom ; n <= v->nsto ; n++ ) {
    b[n] = lfs1fs ;
  }
  for ( n = v->nnfrom ; n <= v->nnto ; n++ ) {
    b[n] = lfn1fn ;
  }
}
  
static void set_up_priors ( double *b , mnc_vectors *v , data_creation_param *dc ) {
  int n ;

  for ( n = v->nsfrom ; n <= v->nsto ; n++ ) {
    b[n] = dc->fs ;
  }
  for ( n = v->nnfrom ; n <= v->nnto ; n++ ) {
    b[n] = dc->fn ;
  }
}
  
static int evaluate_feasibility ( data_creation_param *dc )    
{
  int status = 0 ;
/*
   compute Shannon capacity and compare with info content  
*/

  dc->h2fs =  h2(dc->fs) ;
  dc->h2fn =  h2(dc->fn) ;
  dc->rho =  (double ) dc->N / (double ) dc->M  ;
  dc->rate = dc->rho * dc->h2fs ;
  dc->capacity = 1.0 - dc->h2fn ;
  if ( !dc->gc ) {
    if ( dc->rate >   dc->capacity ) {
      printf ( "task impossible, signed Shannon\n" ) ; 
      if ( dc->notabovecap ) {status -- ; }
    }
    if ( status < 0 || dc->verbose ) 
      printf ( "H(s) = %f, capacity of channel = %f.\n" , 
	      dc->h2fs , dc->capacity ) ; 
  }
  return status ; 
}


static void dc_defaults ( data_creation_param *dc ) 
{
#include "dc_var6_def.c"
}

static int process_command ( int argc , char **argv , mnc_all *all ) {
  data_creation_param *dc = all->dc ; 
  fe_min_control *c = all->c ; 
  fe_min_control *c2 = all->c2 ; 
  fe_min_control *c3 = all->c3 ; 
  bnd_control *bndc = all->bndc ; 

  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;
      dc->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-VV") == 0 )        {
      c->verbose = 2;
      dc->verbose = 2;
    }
#include "fe_var6_clr.c"
#include "fed2_var_clr.c"
#include "dc_var6_clr.c"
#include "bnd_var_clr.c"
    else if ( strcmp (argv[i], "-b") == 0 ) {
      if ( i + 3 >= argc ) { ERROR2;      }
      else {
	cs *= sscanf(argv[++i], "%d", &(c->betastyle)); 
	cs *= sscanf(argv[++i], "%lf", &(c->beta0)); 
	cs *= sscanf(argv[++i], "%lf", &(c->beta1)); 
      }
    }
    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 , all ) ;
  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 ,
			 mnc_all *all 
)
{
  data_creation_param *dc = all->dc ; 
  
  fprintf( fp, "Usage: %s ",argv[0]);
  fprintf( fp, " [optional arguments]");

  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)"); 
  NLNE; fprintf( fp, " Data creation:" ) ; 
#include "dc_var6_usg.c"
  pause_for_return();
  NLNE; fprintf( fp, " Inference:" ) ; 
  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" ); 
  fprintf( fp, "\n");
  fe_print_usage ( argv , fp , all ) ; 
  pause_for_return();
  bnd_print_usage ( argv , fp , all ) ; 
  return ;
}

static void fe_print_usage ( char **argv , FILE * fp ,
			    mnc_all *all 
)
{
  fe_min_control *c = all->c ; 
  fe_min_control *c2 = all->c2 ; 
  fe_min_control *c3 = all->c3 ; 

  NLNE; fprintf( fp, " Further free energy minimization stuff: <defaults>");
#include "fe_var6_usg.c"
#include "fed2_var_usg.c"
  fprintf( fp, "\n");
  return ;
}

static void bnd_print_usage ( char **argv , FILE * fp ,
			    mnc_all *all 
)
{
  bnd_control *bndc = all->bndc ; 

  NLNE; fprintf( fp, " Belief Net decoder:     <defaults>");
#include "bnd_var_usg.c"
  fprintf( fp, "\n");
  return ;
}

#undef DNT
#undef NLNE

static void mnc_free ( mnc_all *all ) {
/*  data_creation_param *dc = all->dc ;*/
  mnc_vectors *vec = all->vec ; 

  free_cvector ( vec->xo, 1 , vec->N ) ; 
  free_cvector ( vec->t , 1 , vec->M ) ; 

  free_cvector ( vec->x , 1 , vec->N ) ; 
  free_cvector ( vec->y , 1 , vec->M ) ; 
  free_cvector ( vec->z , 1 , vec->M ) ; 

  free_dvector ( vec->bias , 1 , vec->N ) ; 
  free_dvector ( vec->n , 1 , vec->N ) ; 

}

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

/*
<!-- hhmts start -->
Last modified: Sat Aug 23 16:46:51 1997
<!-- hhmts end -->
*/