mnc4new.c

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

      count_n += v->x[u] ;
    }
    p->true_s = 1 ; 
    break;
  default :
    printf ("Error in make_vectors_quick\n" ) ; 
    pause_for_return();exit(0);
    break;
  }
  dc->true_fs = (double) count_s / (double) dc->N ; 
  dc->true_fn = (double) count_n / (double) dc->M ; 

  alist_times_cvector_mod2 ( &(p->a) , v->x , v->z ) ; 

  if ( dc->noisy ) {
    fprintf ( stderr, "unchecked code in use\n");
    for ( m = 1 ; m <= dc->M ; m++ , u++) {
      v->z[m] ^= v->y[m] ;
    }
  }

  switch ( dc->MNC ) {
  case(0): case(1): 
    for ( m = 1 ; m <= p->M ; m++ ) {
      p->g[m] = lfn1fn ;
      if ( v->z[m] ) p->g[m] *= -1.0 ; 
    }
    break ; 
  case(2):
    for ( m = 1 ; m <= dc->M ; m++ ) {
      p->g[m] = lfn1fn ;
      if ( v->z[m] ) p->g[m] *= -1.0 ; 
    }
    for (  ; m <= p->M ; m++ ) {
      p->g[m] = lfs1fs ;
    }
    break ; 
  case(3):
  case(4):
    for ( m = 1 ; m <= p->M ; m++ ) {
      p->g[m] = -1.0 ; /* annealing will be needed here */
      if ( v->z[m] ) p->g[m] *= -1.0 ; 
    }
    break ; 
  default:
    fprintf ( stderr , "MNC error \n" ) ; 
    break ; 
  } 
  /* 
     prior / bias
                     */
  switch ( dc->MNC ) {
  case(0): case(1):
    for ( n = 1 ; n <= p->N ; n ++ ) 
      p->bias[n] = lfs1fs ;
    break ; 
  case(2):
    for ( n = 1 ; n <= p->N ; n ++ ) 
      p->bias[n] = 0.0 ;
    break ; 
  case (3) : case ( 4) : 
    for ( n = 1 ; n <= dc->N ; n ++ ) 
      p->bias[n] = lfs1fs ;
    for (       ; n <= p->N ; n ++ ) 
      p->bias[n] = lfn1fn ;
    break ; 
  default:
    fprintf ( stderr , "MNC error \n" ) ; 
    status -- ; 
    break ; 
  } 

  return status ; 
}

  
static int evaluate_feasibility ( data_creation_param *dc , fe_min_param *p )    
{
  double fs =   dc->fs ;
  double fn =   dc->fn ;
  int status = 0 ;
  int MNC = dc->MNC ; 
  int tmpi ; 
  double ft , tmpd ; 

/*
   compute Shannon capacity and compare with info content  
*/

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

/* theoretical calculation of ft, using average number of 1's ; 
  only appropriate to C1 and C3 */
  if ( MNC == 1 || MNC == 3 ) {
    tmpi = p->a.tot / p->a.M ;
    for ( ft = 0.0 ; tmpi > 0 ; tmpi -- )
      ft = fs * ( 1- ft ) + (1 - fs ) * ft ; 

    tmpd = fn * (1.0 -ft) + (1-fn) * ft ; 
    printf ( "fs = %9g , ft = %9g (theo), ft = %9g (emp), fn = %9g , fz = %9g (theo)\n" , fs, ft , dc->fte , fn, tmpd) ; 
    printf ( "H(s) = %6f, H(z;t) = %6f = %6f - %6f ; capacity of channel 1 - H(n) = %6f. " , 
	    h2(fs) , h2(tmpd) - h2(fn) , h2(tmpd) , h2(fn) , 1 - h2(fn) ) ; 
    if ( (double ) dc->N * h2(fs) > (double ) dc->M * ( h2(ft) - h2(fn) ) ) {
      printf ( "task impossible because output ensemble is suboptimal\n" ) ; 
      status -= 10 ; 
    }
    else printf ( "\n" ) ; 
  }

  return status ; 
}


static void dc_defaults ( data_creation_param *dc ) 
{
#include "dc_var_def.c"
  dc->MESSAGE = 1 ; 
  dc->failures = 10 ; 
  dc->failcount = 0 ; 
  dc->mseed = 1234 ; 
  dc->vseed = 9876 ; 
  dc->N = 20 ; 
  dc->M = 0 ; 
  dc->MNC = 0 ; 
  dc->fs = 0.1 ;     /* fraction of signal bits AND/OR assumed prior */
  dc->fpol = 0.1 ; 
  dc->npol = 0 ; 
  dc->fn = 0.15 ;    /* fraction of noise bits AND/OR assumed likelihood */
  dc->nn = 0 ;                        /* if !0, number of noise bits */
  dc->ns = 0 ;                        /* if !0, number of signal bits */
  dc->DEMO = 0 ; 
  dc->verbose = 0 ;
  dc->pC1.per_row = 5 ; 
  dc->pC2.per_row = 5 ; 
  dc->fe = 1 ; 
}

static int process_command ( int argc , char **argv , mnc_all *all ) {
  data_creation_param *dc = all->dc ; 
  fe_min_control *c = all->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;
      dc->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-DV") == 0 )        {
      dc->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-FV") == 0 )        {
      c->verbose = 1;
    }
    else  if ( strcmp (argv[i], "-FVV") == 0 )        {
      c->verbose = 2 ;
    }
    else  if ( strcmp (argv[i], "-VV") == 0 )        {
      c->verbose = 2;
      dc->verbose = 2;
    }
    else  if ( strcmp (argv[i], "-SV") == 0 )        {
      c->seqverbose = 1;
    }
    else  if ( strcmp (argv[i], "-SVV") == 0 )        {
      c->seqverbose = 2;
    }
    else  if ( strcmp (argv[i], "-CG") == 0 )        {
      c->CG = 1;
    }
/* 
   fe stuff
*/
#include "fe_var_clr.c"
/*
   control stuff
*/
#include "dc_var_clr.c"
    else  if ( strcmp (argv[i], "-NL") == 0 )        {
      c->NL = 0; 
    }
    else  if ( strcmp (argv[i], "-DEMO") == 0 )        {
      c->DEMO = 1; 
    }
    else  if ( strcmp (argv[i], "-nl") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(c->NL)); 
    }
    else  if ( strcmp (argv[i], "-mess") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->MESSAGE)); 
    }
    else  if ( strcmp (argv[i], "-failures") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->failures)); 
    }
    else  if ( strcmp (argv[i], "-n") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->N)); 
    }
    else  if ( strcmp (argv[i], "-m") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->M)); 
    }
    else  if ( strcmp (argv[i], "-mnc") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->MNC) ); 
    }
    else  if ( strcmp (argv[i], "-opt") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(c->opt)); 
    }
    else  if ( strcmp (argv[i], "-norder") == 0 ) {
      if ( i + 1 == argc ) { ERROR1; }
      else cs *= sscanf(argv[++i], "%d", &(c->NORDER) ) ; 
    }
    else  if ( strcmp (argv[i], "-itmax") == 0 ) {
      if ( i + 1 == argc ) { ERROR1; }
      else cs *= sscanf(argv[++i], "%d", &(c->itmax) ) ; 
    }
    else   if ( strcmp (argv[i], "-seed") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%ld", &(c->seed)); 
    }  
    else   if ( strcmp (argv[i], "-mseed") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%ld", &(dc->mseed)); 
    }  
    else   if ( strcmp (argv[i], "-vseed") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%ld", &(dc->vseed)); 
    }  
    else  if ( strcmp (argv[i], "-seqp") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(c->seq_period)); 
    }
    else if ( strcmp (argv[i], "-e") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(c->epsilon)); 
    }
    else if ( strcmp (argv[i], "-fn") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(dc->fn) ); 
    }
    else if ( strcmp (argv[i], "-sdx") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(c->sdx)); 
    }
    else if ( strcmp (argv[i], "-betap") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(c->betap)); 
    }
    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 if ( strcmp (argv[i], "-ftol") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(c->ftol)); 
    }
    else if ( strcmp (argv[i], "-fs") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(dc->fs)); 
    }
    else if ( strcmp (argv[i], "-fpol") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &(dc->fpol)); 
    }
    else if ( strcmp (argv[i], "-t1") == 0 ) {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->pC1.per_row ) ); 
    }
    else if ( strcmp (argv[i], "-t2") == 0 ) {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(dc->pC2.per_row ) ); 
    }
    else if ( strcmp (argv[i], "-pheading") == 0 ) {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &(c->pheading_period)); 
    }
    else  if ( strcmp (argv[i], "-o") == 0 )	{
      if ( i + 1 == argc ) 	    { ERROR1;      }
      else {
	strcpy(c->outfile, argv[++i]);
	if ( !(c->printout) ) c->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 , 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 ; 
  fe_min_control      *c  = all->c  ;
  
  fprintf( fp, "Usage: %s ",argv[0]);
  fprintf( fp, " [optional arguments]");

  NLNE; fprintf( fp, " Verbosity of run-time output: <defaults>");
  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)"); 
  DNT;  fprintf( fp, "-FV | -FVV               (verbose in fe_min)"); 
  DNT;  fprintf( fp, "-DV | -DVV               (verbose in data creation)"); 
  DNT;  fprintf( fp, "-SV | -SVV               (verbose in seq_opt)"); 
  DNT;  fprintf( fp, "-DEMO                    (only do demo of MNC) "); 
  DNT;  fprintf( fp, "-CG                      (check gradient) "); 
  DNT;  fprintf( fp, "[-e epsilon]     <%g> (step size for gradient)", c->epsilon); 
  DNT;  fprintf( fp, "-pheading period <%d>" , c->pheading_period ); 

  NLNE; fprintf( fp, " Data creation:" ) ; 
  DNT;  fprintf( fp, "-mnc MNC         <%d> (which code to use [0;1,2,3,4]) " , dc->MNC ); 
  DNT;  fprintf( fp, "-n N             <%d> (number of bits in string) " , dc->N ); 
  DNT;  fprintf( fp, "-m M             <%d> (number of bits in transmission) " , dc->M ); 
  DNT;  fprintf( fp, "-t1 t1           <%d> (C1 number of 1s per row) " , dc->pC1.per_row ); 
  DNT;  fprintf( fp, "-t2 t2           <%d> (C2 [if appropriate] number of 1s per row) " , dc->pC2.per_row ); 
  DNT;  fprintf( fp, "-fn fn           <%g> (error probability) " , dc->fn ); 
  DNT;  fprintf( fp, "-fs fs           <%g>  (fraction in string)", dc->fs); 

  NLNE; fprintf( fp, " Inference:" ) ; 

  DNT;  fprintf( fp, "-sdx sdx         <%3g> (sd of initial state)" , c->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)" , c->betap ); 
  DNT;  fprintf( fp, "-opt optimizer   <%d>   (0=frp,1=macopt,2=jumpopt,3=seqopt,4=?,5=strangeopt)", c->opt); 
  DNT;  fprintf( fp, "-itmax itmax     <%d>   (num of linmins or jumps per optimization)", c->itmax); 
  DNT;  fprintf( fp, "-norder n        <%d>   (for seqopt - 0/1/2 = 1..N / random / suspicion / suspicion/num )", c->NORDER); 
  DNT;  fprintf( fp, "-mess MESSAGES   <%d>   (number of messages to send)", dc->MESSAGE); 
  DNT;  fprintf( fp, "-failures failures <%d>   (number of failed messages to stop after)", dc->failures); 
  DNT;  fprintf( fp, "-nl loops        <%d>   (number of cg runs)", c->NL); 
  DNT;  fprintf( fp, "-ftol ftol       <%g> ", c->ftol ); 
  DNT;  fprintf( fp, "-seed seed "); 
  DNT;  fprintf( fp, "-mseed mseed "); 
  DNT;  fprintf( fp, "-vseed vseed "); 
  DNT;  fprintf( fp, "-seqp seq_period   <%d>  (how many loops per check of delta F)" , c->seq_period ); 
  NLNE; fprintf( fp, " Files:");
  DNT;  fprintf( fp, "-o outfile               "); 
  fprintf( fp, "\n");
#include "dc_var_usg.c"
  fe_print_usage ( argv , fp , all ) ; 
  return ;
}

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

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

#undef DNT
#undef NLNE

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

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

  free_ivector ( p->a.norder , 1, p->N ) ; 

  free_dvector (   p->x , 1 , p->N ) ;
  free_dvector (   p->g , 1 , p->M ) ; 
  free_dvector (   p->bias    , 1 , p->N ) ; 
  free_cvector (   p->so , 1 , p->N ) ;

}

double h2 ( double x ) {
  double tmp ; 
  tmp = x * log ( x ) + ( 1.0 - x ) * log ( 1.0 - x ) ;
  return - tmp / log ( 2.0) ; 
}

/*
<!-- hhmts start -->
Last modified: Thu Jun 13 17:46:59 1996
<!-- hhmts end -->
*/