mnc3.c

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

/* 
   mnc3.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

   - free energy minimization for the MacKay-Neal Error Correcting Code - 

   - general purpose program handling many different codes              -

     main's central bit sends many messages using a single code 

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

   Notes for mnc, mnc2 are in the 1995 red book, 95 01 02.

mnc3: 

   New addition = neural net deconvolution options. 
   Strategy is based on M-S and Gallager ideas. Have an estimator 
   for the worst bit, and `correct' it (a decision); then recompute
   the parity checks, and correct another. 

   This then needs an initial training phase / weight loading.

mnc2:

   Data generation:

          For MNC 1-4, The code C is defined by one or two
	  N*N matrices represented explicitly. 

          A second representation is used for the binary matrix A[m][n]
	  involved in the inference problem As + n = z. It is represented by 
	  a pair of arrays of lists. In the `mlist' mlist[m][l], 
	  l=1...num_mlist[m], each element is an `n' such that A[m][n] = 1.
	  Thus the list lists for each m where the 1's are in row m.
	  There is an equivalent `nlist' listing where the 1's are in each 
	  column. 

	  There are four MN codes, each with a different construction
	  for C and / or a different definition of A,s,n,z. 

	  The code is applied to a sparse string s, giving t = Cs. 
	  Noise is added to this, so that Cs + n = z. 
	  But these s, n and z's are not necessarily the same as the 
	  s, n, z of the fe algm.

	  Code 1:

	  to define A = C and plug the problem directly into the 
	  free energy machine, with the approximating distribution 
	  defined on s: Q(s) = prod q_i(s_i). In this case A is an 
	  N * N matrix. 

	  Code 2:

	  The `Radford representation' has separate C into 
	  a product C = C1 C2, where C1 is sparse and C2 has a sparse
	  inverse, maybe of  similar density 
	  to C1.

	  Then we can rewrite the equations

	  Cs + n = z
	  Is + s = 0

	  as 

	  [A1]u + n = z
	  [A2]    s   0

	  where A1 = C1, A2 = C2^{-1}, and u = C2 s. 
	  We can then apply the approximation in the intermediate
	  representation u. The matrix A is an M*N matrix with M=2N. 

	  Each code calls the fe_min routine when things are set up.


          You get to choose (on command line): 
	  N (string length),                                      (e.g. 100)
	  fn  (probability of error in observation of t = As)    (e.g. 0.1)
	  fs  (probability of high bit in s itself)              (e.g. 0.1)
	  seed (for random numbers)



*/

#include "./ansi/r.h"
#include "./ansi/rand.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./ansi/macopt.h"
#include "./thing_sort.h"
#include "./fe.h"
#include "./mnc.h"
#include "./mnc3.h"
#include "./nn.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 void   nn_print_usage ( char ** , FILE * ,  mnc_all *  );
static int make_sense ( data_creation_param * , fe_min_control * , mnc_all * ) ; 
static int make_code ( data_creation_param * , mnc_code * , fe_min_param * ) ;
static int report_code ( data_creation_param * , mnc_code * , fe_min_param * ) ;
static int make_space ( data_creation_param * , fe_min_param * , 
		       mnc_code * , mnc_vectors * ) ;
static int make_norder ( alist_matrix * , fe_min_control * ) ;
static int make_vectors ( data_creation_param * , fe_min_param * , 
		       mnc_code * , mnc_vectors * ) ;
static int decode ( fe_min_param * , 
			  fe_min_control * , data_creation_param * , 
			  mnc_vectors  * , mnc_code * ) ; 
static int score ( fe_min_param * , 
			  fe_min_control * , data_creation_param * , 
			  mnc_vectors  * , mnc_code * ) ; 
static int make_pvectors ( fe_min_param * , 
			  fe_min_control * , data_creation_param * , 
			  mnc_vectors  * );
static int evaluate_feasibility ( data_creation_param * , fe_min_param * ) ;
static int check_alist ( alist_matrix * , fe_min_param * , 
			fe_min_control * , data_creation_param * ) ;

static void   finalline (  FILE * , fe_min_control * , data_creation_param * , fe_min_param * ) ;

static void mnc_free ( mnc_all * ) ; 

void   main ( int , char ** ) ;

static double h2 ( double ) ;

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

  fe_min_param        p    ;
  fe_min_control      c    ;
  data_creation_param dc   ;
  mnc_code            code ; 
  mnc_vectors         vec  ; 
  mnc_net             net ; 
  mnc_net_control     netc ;
  mnc_all             all  ;

  unsigned long histid1, histid2 , orig_size , current_size ; 
#ifdef _DEBUG_MALLOC_INC

  union dbmalloptarg  m;
  m.i = M_HANDLE_CORE | M_HANDLE_DUMP;
  dbmallopt(MALLOC_WARN,&m);

  m.i = M_HANDLE_ABORT;
  dbmallopt(MALLOC_FATAL,&m);

/*  m.str = "log";
  dbmallopt(MALLOC_ERRFILE,&m); */
#endif

  all.net = &net ; 
  all.nc = &netc ; 
  all.dc = &dc ; 
  all.p = &p ; 
  all.c = &c ; 
  all.code = &code ;  
  all.vec = &vec ; 

  fe_defaults ( &p , &c ) ; 
  dc_defaults ( &dc ) ; 
  nn_defaults ( &net , &netc ) ;

  if ( process_command (argc, argv, &all ) < 0 ) exit (0) ;
  if ( make_sense ( &dc , &c, &all ) < 0 ) exit (0) ;
  if ( make_space ( &dc , &p , &code , &vec ) < 0 ) exit (0) ;
  orig_size = malloc_inuse ( &histid1 ) ; 
  if ( dc.usenet == 1 ) { 
    if ( nn_make_sense ( &netc , &all ) < 0 ) exit (0) ;
  }

  if ( dc.verbose ) {
    fprintf(stderr,"mnc MNC=%d, 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.MNC , dc.N , dc.M , p.N , p.M ,  dc.fn , dc.fs , dc.nn , dc.ns , c.opt , c.gamma , c.metric , c.betastyle ) ;
    fflush(stderr);
  } else {
    fprintf(stdout,"MNC=%d opt=%2d:%2d:%2d ",
	    dc.MNC , c.opt , c.gamma , c.metric ) ;
  }

  if ( make_code ( &dc , &code , &p ) < 0 ) { 
    report_code ( &dc , &code , &p ) ; exit ( 0 ) ;     
  }
  if ( check_alist ( &(p.a) , &p , &c , &dc ) < 0 ) { 
    report_code ( &dc , &code , &p ) ; exit ( 0 ) ;
  }
  if ( dc.fe == 1 ) {
    fe_allocate ( &p , &c ) ; 
    if ( make_norder ( &(p.a) , &c ) < 0 ) exit ( 0 ) ; 
  } else if ( dc.usenet == 1 ) {
    nn_allocate ( &p , &net , &netc , &all ) ; 
    nn_initialize ( &net , &netc , &all ) ; 
                           /* this includes training the thing if necessary,
			      and saving its weights */
  }

  if ( evaluate_feasibility ( &dc , &p ) < 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 ) ; 
  }

  for ( dc.message = 1 ; 
       ( dc.message <= dc.MESSAGE ) && 
       ( ( dc.failures==0 ) || ( dc.failcount < dc.failures ) ) ; 
       dc.message ++ ) {
    if ( dc.usenet == 1 ) {
      exit(0);
    }
    else {
      if ( make_vectors ( &dc  , &p, &code, &vec) < 0 ) exit (0) ; 
      if ( make_pvectors ( &p , &c , &dc , &vec ) < 0 ) exit ( 0 ) ; 
      
      fe_min ( &p , &c ) ; 

    }

    if ( decode ( &p , &c , &dc , &vec , &code ) < 0 ) exit ( 0 ) ; 
    if ( score ( &p , &c , &dc , &vec , &code ) < 0 ) exit ( 0 ) ; 
    if ( dc.verbose > 0 ) finalline ( stdout , &c , &dc , &p ) ;

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

  fe_free ( &p , &c ) ; 
  mnc_free ( &all ) ; 
#ifdef _DEBUG_MALLOC_INC
  current_size = malloc_inuse ( &histid2 ) ; 
  if ( current_size != orig_size && dc.verbose ) {
/*    malloc_list ( 2 , histid1 , histid2 ) ; */
    fprintf ( stderr , "memory cockup\n" ) ; 
    printf ( "mZ:%d\n" , malloc_chain_check(1) ) ; fflush(stdout) ;
  }
#endif
}

void finalline ( FILE *fp , fe_min_control *c , data_creation_param *dc , fe_min_param *p ) 
{
  if ( c->pheading_period && ! ( ( c->pheading_count ++ ) % c->pheading_period )  ) 
    fprintf ( fp , 
                 "#es rl.scre  N    t    fs     fn    ns nn  bstl  seed  mseed  vseed tfs      tfn       beta   h     h_dev  MNC ms#\n") ; 

  fprintf ( fp , "%4d %6.4g %3d %3d %7.4g %7.4g %2d %2d %d %8ld %8ld %8ld %7.4g %7.4g %7.4g %7.4g %7.4g %d %4d %4d\n" ,
                  p->count , p->v , dc->N , dc->pC1.per_row , dc->fs , dc->fn , dc->ns , dc->nn , c->betastyle , c->seed  , dc->mseed  , dc->vseed  , dc->true_fs, dc->true_fn , p->beta , p->h_soln , p->h_dev , dc->MNC , dc->message , dc->failcount ) ;
}


static int make_sense ( data_creation_param *dc ,  fe_min_control *c , mnc_all *all ) 
{ /* routine to correct silly control parameters */
/* 
   first the data creation 
                             */
  int status = 0 ; 

  switch ( dc->MNC ) {
  case(0) :
    if ( dc->M == 0 ) {
/* if MNC == 0, M and N better had got set on command line */
      fprintf ( stderr , " No M specified \n" ) ;  status-- ; }
    break ; 
  case ( 1) : 
  case ( 2) : 
  case ( 3) : 
  case ( 4) : 
    dc->M = dc->N ; 
    break ; 
  default : 
    fprintf ( stderr , "invalid MNC = %d\n" , dc->MNC ) ; 
    break ; 
  }

/* 
   and now the free energy control
                                       */
  if ( ( dc->MNC ==3 || dc->MNC == 4 || c->NL > 1 )  
      && ( !(c->opt==5 || c->opt==6) )
      && ( !(dc->usenet) ) && 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 ) {
    all->p->zero_temperature = 1 ; 
  }

  c->dgamma = (double) (c->gamma) ;
  return status ; 
}

static int make_space ( data_creation_param *dc , fe_min_param *p , 
		       mnc_code *code , mnc_vectors *vec ) {
/* see also mnc_free */
  int status = 0 ; 
  int tr , tc ; 

  malloc_enter ("make_space") ;
  if ( dc->verbose ) printf ( "Making space\n" ) ; 
  /*
     set up size of the A matrix
                                   */
  switch ( dc->MNC ) {
  case ( 0 ) : 
    p->M = dc->M ;     p->N = dc->N ; 
    break ; 
  case ( 1) :
    p->M = dc->N ;     p->N = dc->N ; 
    break ; 
  case( 2) :
    p->M = 2 * dc->N ; p->N = dc->N ; 
    break ; 
  case(3) : case ( 4 ) : 
    p->M = dc->N ;     p->N = 2 * dc->N ; 
    break ; 
  default:
    fprintf ( stderr , "invalid MNC = %d\n" , dc->MNC ) ; status -- ; 
    break ; 
  } 

/* 
   Now estimate size of alist matrix 

     MNC1: Number per row and column is almost fixed to dc.pC1.per_row
     MNC2: Per row is the max of C1 and C2. Per col is the sum plus a bit.
     MNC3: As MNC1, plus 1. 
     MNC4: converse of MNC2.
     MNC0: identity matrix followed by a load of sparse junk 
*/
  switch ( dc->MNC ) {
  case ( 1 ) : 
    tr = dc->pC1.per_row + 2 ; 
    tc = tr + 10 ; break ;
  case ( 2 ) : 
    tr = MAX ( dc->pC1.per_row , dc->pC2.per_row ) + 2 ; 
    tc = 2 * tr + 10 ; break ;
  case ( 3 ) : 
    tr = dc->pC1.per_row + 3 ; 
    tc = tr + 10 ; break ;
  case ( 4 ) : 
    tc = MAX ( dc->pC1.per_row , dc->pC2.per_row ) + 10 ; 
    tr = 2 * tc ; break ;
  default : 
  case ( 0 ) :
    tr = dc->tr + 1 ;
    tc = 1 + ( p->M - p->N ) * tr * 3 / ( 2 * p->N )  ; /* 3/2 for safety */
    break ;
  }
  dc->tr = tr ; dc->tc = tc ; 

  /* Set up memory , etc. */

  vec->s = cvector ( 1 , dc->N ) ;
  if ( dc->MNC == 2 || dc->MNC == 4 ) 
    vec->u = cvector ( 1 , dc->N ) ;
  if ( dc->MNC == 2 ) 
    vec->uo = cvector ( 1 , dc->N ) ; /* u-out */
  else if ( dc->MNC == 4 || dc->MNC == 3 )
    vec->sno = cvector ( 1 , p->N ) ;

  vec->so = cvector ( 1 , dc->N ) ; /* s-out */
  vec->t = cvector ( 1 , dc->M ) ; 
  vec->noise = cvector ( 1 , dc->M ) ; 
  vec->r = cvector ( 1 , dc->M ) ; 
  if ( dc->MNC == 4 ) vec->cr = cvector ( 1 , dc->M ) ; 

  initialize_alist ( &(p->a) , p->N , p->M , tc , tr ) ;
  vec->x = cvector ( 1 , p->N ) ; 
  vec->xpartial = cvector ( 1 , p->N ) ; 
  vec->xrecon = cvector ( 1 , p->N ) ; 
/*  vec->y = cvector ( 1 , p->M ) ;  */
  vec->z = cvector ( 1 , p->M ) ; p->z = vec->z ; 
  vec->touches = ivector ( 1 , p->M ) ; 

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

/* space for some of the optimizer's stuff */
  p->g = dvector ( 1 , p->M ) ; 
  p->bias = dvector ( 1 , p->N ) ; 
  p->s = cvector ( 1 , p->N ) ;
  p->so = cvector ( 1 , p->N ) ;
  p->t = cvector ( 1 , p->M ) ;


  malloc_leave ("make_space");
  return status ; 

}

static int make_code ( data_creation_param *dc , mnc_code *code , fe_min_param *p ) {
  int status = 0 ; 
  FILE   *fp  ; 
  char 	junk[100] ;