mnc7.c

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

/*
   mnc7.c
   mnc6.c                                    (c) DJCM 94 07 04
                                                      94 10 26
						      94 10 31
						      95 01 02
						      95 02 23 (from mnc2 to mnc3)
						      95 04 01 continued work on mnc3
						      95 05 01 to mnc4
						      95 05 27 to mnc5
						      95 08 31 to mnc6
						      95 11 04 to mnc7

   - free energy minimization for the MacKay-Neal Error Correcting Code - 
   - or belief net decoder                                              -
   - mnc7 includes option of gaussian channel                           -

CHANGES relative to mnc4new, mnc4:

   - specialised for MNC4 only                                          -

   - all structures defined in _var files:
       fe_var6
       dc_var6
       fed2_var (there was no need to change this to fed6_var)

   - space is made for the fe min or for the bnd, then the scoring 
     routine gets to look at the results via pointers in the mnc_vector
     structure. Thus scoring is divorced from fe. 

   ======================================================================

   This code is (c) David J.C. MacKay 1994, 1995. It is free software 
   as defined by the free software foundation. 

   What this program does: 

     reads in a matrix A in alist format.
     do {
       generates a  random vector x, 
       shoves it through A, adds y if appropriate, -> z = Ax + y
       do {
         gives b, g (z) to the minimization machine, or A,z to the bnd
         checks the answer. 
	 next strategy
       }
       until ( success OR no more strategies )
     }
     while ( not many failures )

*/

#include "./ansi/r.h"
#include "./ansi/rand2.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./ansi/macopt.h"
#include "./thing_sort.h"
#include "./fe6.h"
#include "./bnd.h"
/* this must come last: */
#include "./mnc6.h"

static void dc_defaults ( data_creation_param * ) ; 
static int    process_command ( int , char **  , mnc_all * ) ; 
static void   print_usage ( char ** , FILE * ,  mnc_all *  );
static void   fe_print_usage ( char ** , FILE * ,  mnc_all *  );
static int make_sense ( data_creation_param * , mnc_all * ) ; 
static int make_sense_c ( data_creation_param * , fe_min_control * ) ; 
static int make_space ( data_creation_param * , mnc_vectors * ) ;
static int make_norder ( alist_matrix * , fe_min_control * ) ;
static int make_vectors_quick ( data_creation_param * , mnc_vectors * ,  mnc_all * ) ;
static void process_vector_fe ( mnc_vectors *v , fe_min_param *p ) ;
static int make_gaussian_noise_bits_and_fix_biases ( unsigned char * ,
						    double *,
						    double  ,
						    int , int  ) ;
static void set_up_biases ( double * , mnc_vectors * , data_creation_param * ) ; 
static void set_up_priors ( double * , mnc_vectors * , data_creation_param * ) ; 
static int score ( data_creation_param * , mnc_vectors  * , mnc_all * ) ; 
static int evaluate_feasibility ( data_creation_param * ) ;

static void   finalline (  FILE * , mnc_all * , data_creation_param * ) ;

static void mnc_free ( mnc_all * ) ; 
static int check_alist_MN ( alist_matrix * , mnc_vectors *vec ) ; 
static int hook_mnc_vec_to_bnd ( bnd_param *p , mnc_vectors *vec ) ;
static int hook_mnc_vec_to_fe ( fe_min_param *p , mnc_vectors *vec ) ;
static void bnd_print_usage ( char **argv , FILE * fp , mnc_all *all );
static double h2 ( double ) ;

void   main ( int , char ** ) ;

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

  fe_min_param        p    ;
  fe_min_control      c    ;
  fe_min_control      c2   ; /* backup strategy */
  fe_min_control      c3   ; /* backup strategy */
  bnd_control         bndc ; 
  bnd_param           bndp ;
  data_creation_param dc   ;
  mnc_vectors         vec  ; 
  mnc_all             all  ;

  all.dc   = &dc ; 
  all.p    = &p ; 
  all.c    = &c ; 
  all.bndp = &bndp ; 
  all.bndc = &bndc ; 
  all.c2   = &c2 ; 
  all.c3   = &c3 ; 
  all.vec  = &vec ; 
  all.a    = &(p.a) ; 
  bndp.a   = &(p.a) ; 

  fe_defaults_p ( &p ) ; 
  fe_defaults_c ( &c ) ; 
  fe_defaults_c ( &c2 ) ; 
  fe_defaults_c ( &c3 ) ; 
  dc_defaults   ( &dc ) ; 
  bnd_defaults  ( &bndp , & bndc ) ; 

  if ( process_command (argc, argv, &all ) < 0 ) exit (0) ;
  if ( make_sense ( &dc , &all ) < 0 ) exit (0) ;

  fprintf(stderr,"mnc6 N0=%d, M0=%d, N=%d, M=%d, fn=%6.3g, fs=%6.3g, nn=%d, ns=%d, opt=%d:%d:%d, bstyle=%d\n",
	  dc.N , dc.M , vec.N , vec.M , dc.fn , dc.fs , dc.nn , dc.ns , c.opt , c.gamma , c.metric , c.betastyle ) ;
  fflush(stderr);
  if ( dc.use_fe && c.second ) {
    fprintf(stdout,"opt=%2d:%2d:%2d ",
	     c.opt , c.gamma , c.metric ) ;
    fprintf(stdout,"opt2=%2d:%2d:%2d ",
	     c2.opt , c2.gamma , c2.metric ) ;
    fprintf(stdout,"opt3=%2d:%2d:%2d ",
	     c3.opt , c3.gamma , c3.metric ) ;
  }

  if ( make_space ( &dc , &vec ) < 0 ) exit (0) ;

  if ( read_allocate_alist ( &(p.a) , dc.afile ) < 0 ) exit (0) ; 
  if ( check_alist_MN ( &(p.a) , &vec ) < 0 ) exit (0) ; 

  if ( dc.use_fe >= 1 ) {
    c.mega_needed = ( c.opt == 3 || 
		   ( c.second && c2.opt == 3 ) || 
		   ( c.third && c3.opt == 3 ) ) ;
    hook_mnc_vec_to_fe (  &p ,  &vec ) ;
    fe_allocate2 ( &p , c.mega_needed ) ; 
    set_up_biases ( p.bias , &vec , &dc ) ; 
    if ( make_norder ( &(p.a) , &c ) < 0 ) exit ( 0 ) ; 
  }
  if ( dc.use_bnd >= 1 ) {
    hook_mnc_vec_to_bnd (  &bndp ,  &vec ) ;
    bnd_allocate ( &bndp , &bndc ) ; 
    set_up_priors ( bndp.bias , &vec , &dc ) ; 
  }

  if ( evaluate_feasibility ( &dc ) < 0 ) exit ( 0 ) ; 

  if ( c.writelog ) {
    fp = fopen ( c.logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , c.logfile ) ; 
      c.writelog = 0 ; 
    } else fclose (fp ) ; 
  }

  if ( bndc.writelog ) {
    fp = fopen ( bndc.logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , bndc.logfile ) ; 
      bndc.writelog = 0 ; 
    } else fclose (fp ) ; 
  }

  if ( dc.error_log ) {
    fp = fopen ( dc.error_logfile , "w" ) ;
    if ( !fp ) {
      fprintf ( stderr , " couldn't open logfile %s\n" , dc.error_logfile ) ; 
      dc.error_log = 0 ; 
    } else fclose (fp ) ; 
  }

  /*
      MAIN LOOP
                 */
  ran_seed ( dc.vseed ) ; 
  for ( dc.message = 1 ; 
       ( dc.message <= dc.MESSAGE ) && 
       ( ( dc.failures==0 ) || ( dc.failcount < dc.failures ) ) ; 
       dc.message ++ ) {

    make_vectors_quick ( &dc , &vec , &all ) ; 

    if ( dc.use_fe >=1 ) {
      process_vector_fe ( &vec , &p ) ; 
      all.controller = 1 ; 
      if ( ( fe_min ( &p , &c ) > 0 ) && ( c.second ) ) {
	all.controller ++ ;
	if ( ( fe_min ( &p , &c2 ) > 0 ) && ( c.third ) ) { 
	  all.controller ++ ;
	  fe_min (  &p , &c3 ) ; 
	}
      }
    } else {
      all.controller = 0 ; 
      bndecode ( &bndp , &bndc ) ;
    }
    if ( score ( &dc , &vec , &all ) < 0 ) exit ( 0 ) ; 

    if ( dc.verbose > 0 ) finalline ( stdout , &all , &dc ) ;

    if ( c.printout ) {
      fp = fopen ( c.outfile , "a" ) ;
      if( !fp ) {
	fprintf( stderr, "No such file: %s\n", c.outfile ) ;
	finalline ( stderr , &all , &dc ) ;
      }      else 	{
	finalline ( fp , &all , &dc ) ;
	fclose ( fp ) ;
      }
    }
  }

  if ( dc.use_fe >= 1 ) {
    fe_free2 ( &p , c.mega_needed ) ; 
  }
  if ( dc.use_bnd >= 1 ) {
    bnd_free ( &bndp , &bndc ) ;
  }
  mnc_free ( &all ) ; 

}

static void finalline ( FILE *fp , mnc_all *all , data_creation_param *dc )
{
  mnc_vectors *v = all->vec ;
  fe_min_control *c = all->c ; 
  if ( dc->pheading_period && 
      ! ( ( ++ dc->pheading_count ) % dc->pheading_period )  ) 
    fprintf ( fp , 
                 "#es tr_wt recn_wt viols N   M    t    fs     fn    ns nn bstl  seed  mseed  vseed rho rate capacity  tfs      tfn       beta nl<fe>iter bndloops  ms# failures\n") ; 

  fprintf ( fp , "%4d %4d %4d %4d %6d %6d %d " , 
	   v->count , v->count_s , v->count_high , v->viols , 
	   dc->N , dc->M , dc->tc ) ;
  fprintf ( fp , "%7.4g %7.4g %2d %2d %d %1ld %8ld %8ld " ,
	   dc->fs , dc->fn , dc->ns , dc->nn , c->betastyle , 
	   c->seed  , dc->mseed  , dc->vseed  );
  fprintf ( fp , "%7.4g %7.4g %7.4g " ,
	   dc->rho , dc->rate , dc->capacity ) ; 
  fprintf ( fp , "%7.4g %7.4g %7.4g %3d %3d %3d %5d %2d\n" , 
	   dc->true_fs, dc->true_fn , all->p->beta ,all->c->nl , 
	   all->c->iter , all->bndc->loop , dc->message , dc->failcount ) ;
}

/* be nice to record the number of iterations used? */

static int make_sense ( data_creation_param *dc , mnc_all *all ) 
{ /*  correct silly control parameters */
/* 
   first the data creation 
                             */
  int status = 0 ; 
  char junk[200] ; 

  if ( dc->M == 0 ) {
    dc->M = dc->N ; /* default is a square MNC */
  }
                                 /* there is potential confusion 
				    here:
				    the "n,k" of the error correcting code
				    are called "M,N" in the dc-> structure.
				    But when it comes to decoding, 
				    the matrix is A is M * N
				    where M = dc->M (traditionally known as n)
				    and   N = dc->N + dc->M = "k + n" 
				    
				    this comment needs to go to check_alist also */
  all->vec->M = dc->M ; 
  all->vec->N = dc->M + dc->N ; 
  all->vec->NS = dc->N ; 
  all->vec->NN = dc->M ; 
  all->vec->nsfrom = 1 ; 
  all->vec->nsto   = dc->N ; 
  all->vec->nnfrom = dc->N + 1 ; 
  all->vec->nnto   = dc->M + dc->N  ; 

/* default name for afile */

  sprintf ( junk , "MNC%d/%d.%d.%ld" , 
	   dc->MNC , dc->N , dc->tc , dc->mseed ) ; 
  if ( dc->reada==1 )  /* if this is 2 then filename
			  's been put there already */
    sprintf ( dc->afile , "%s" ,  junk ) ; 
  if ( dc->reada==0 )  
    fprintf ( stderr , "oi, how can reada be zero?\n" ) ; 

/* 
   and now the free energy control
                                       */
  if ( dc->use_fe ) {
    status += make_sense_c ( dc , all->c ) ; 
    if ( all->c->second ) 
      status += make_sense_c ( dc , all->c2 ) ; 
    if ( all->c->third ) 
      status += make_sense_c ( dc , all->c3 ) ; 
  }
  return status ; 
}


static int make_sense_c (  data_creation_param *dc , fe_min_control *c ) {
  int status = 0 ; 

  if ( ( dc->MNC ==3 || dc->MNC == 4 || c->NL > 1 )  
      && ( !(c->opt==5 || c->opt==6) )
      && ( c->betastyle == 0 ) )
    fprintf ( stderr , "Warning: no annealing schedule given, MNC = %d, NL = %d\n" , 
	     dc->MNC , c->NL ) ; 
  if ( ( c->NL > 1 || c->betastyle != 0 ) && ( (c->opt==5) ) )
    fprintf ( stderr , "Warning: this optimizer doesn't need annealing or multiple loops, opt = %d, NL = %d\n" , 
	     c->opt , c->NL ) ; 

  if ( c->betastyle == 1 && c->NL > 1 ) 	
    c->betaf = ( c->beta1 - c->beta0 ) / ((double) ( c->NL - 1) ) ; 
  else if ( ( c->betastyle == 2 || c->betastyle == 22 ) && c->NL > 1 )
    c->betaf = exp ( log ( c->beta1 / c->beta0 ) / ((double) ( c->NL - 1) )) ; 

  if ( c->opt > 6 || c->opt < 0 ) {
      fprintf ( stderr , " opt %d not existo\n" , c->opt ) ;   status-- ; }

  if ( c->opt == 5 || c->opt == 6 ) {
    c->zero_temperature = 1 ; /* changed from mnc5 on */
  } 

  c->dgamma = (double) (c->gamma) ;

/* any optimizer-specific stuff? */
  switch ( c->opt ) {
  case ( 0 ) : case ( 1 ) : case(2): case(3): case (4): default: 
    break ; 
      
  case(5):    case(6):
    c->macarg.itmax = c->itmax ; 
    break ;
  }

  return status ; 
}

static int make_space ( data_creation_param *dc , mnc_vectors *vec ) {
/* see also mnc_free */
  int status = 0 ; 

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

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

  vec->bias = dvector ( 1 , vec->N ) ;

  return status ; 
}

static int hook_mnc_vec_to_fe ( fe_min_param *p , mnc_vectors *vec ) {
  int status = 0 ; 
  p->M = vec->M ;
  p->N = vec->N ;
  p->z = vec->z ; 
  p->s = vec->x ; 
  /* where the true vector lives */
  p->so = vec->xo ;