mnc4new.c

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

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

   - 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.

mnc4new:

   a hacked version to allow M not equal to N in MNC4

mnc4:

   Removed the neural net bits. 
   stripped out the code generation into the program code4.

   see mnc3.c for a more comprehensive description.

   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
       gives b, g (z) to the minimization machine, 
       checks the answer. 
     }
     while ( not many failures )

*/

#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"

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 * , fe_min_control * , mnc_all * ) ; 
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_quick ( data_creation_param * , fe_min_param * , 
			mnc_vectors * ) ;
static int score ( fe_min_param * , 
			  fe_min_control * , data_creation_param * , 
			  mnc_vectors  * , mnc_code * ) ; 
static int evaluate_feasibility ( data_creation_param * , fe_min_param * ) ;

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

static void mnc_free ( mnc_all * ) ; 
static int check_alist_p ( fe_min_param * ) ; 


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_all             all  ;

  int n ;
  double tmpd ; 

  all.dc = &dc ; 
  all.p = &p ; 
  all.c = &c ; 
  all.code = &code ;  
  all.vec = &vec ; 

  fe_defaults ( &p , &c ) ; 
  dc_defaults ( &dc ) ; 

  if ( process_command (argc, argv, &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_sense ( &dc , &c, &all ) < 0 ) exit (0) ;
  if ( make_space ( &dc , &p , &code , &vec ) < 0 ) exit (0) ;

/* read in alist */

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

  fe_allocate ( &p , &c ) ; 
  if ( make_norder ( &(p.a) , &c ) < 0 ) exit ( 0 ) ; 

  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 ) ; 
  }

/* toronto */
  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 ) ; 
  }

  ran_seed ( dc.vseed ) ; 
  if ( dc.skip ) {
    for ( dc.message = 1 ; dc.message <= dc.skip ; dc.message ++ ) {
      make_vectors_quick ( &dc , &p , &vec ) ; 
      if ( ! c.init_given ) {
	if ( c.seed != 0 ) ran_seed ( c.seed ) ;  /* is this right? */
	for  ( n = 1 ; n <= p.N ; n++ ) { /* initial condition in fe_min */
	  tmpd = rann() ; 
	}
      }

    }
  }
  for ( dc.message = 1 ; /* is this right? */
       ( dc.message <= dc.MESSAGE ) && 
       ( ( dc.failures==0 ) || ( dc.failcount < dc.failures ) ) ; 
       dc.message ++ ) {

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

    fe_min ( &p , &c ) ; 

    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 ) ; 

}

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 tru_wt recon_wt viols N    t    fs     fn    ns nn  bstl  seed  mseed  vseed tfs      tfn       beta   MNC ms# failures\n") ; 

  fprintf ( fp , "%4d %6.4g %4d %4d %4d %6d %d %7.4g %7.4g %2d %2d %d %1ld %8ld %8ld %7.4g %7.4g %7.4g %d %4d %4d\n" ,
                  p->count , p->v , p->count_s , p->count_high , p->viols , 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 , dc->MNC , dc->message , dc->failcount ) ;
}

static int make_sense ( data_creation_param *dc ,  fe_min_control *c , 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 */
    }

/* default name for afile */

  sprintf ( junk , "MNC%d/%d.%d.%ld" , 
	   dc->MNC , dc->N , dc->pC1.per_row ,
	   dc->mseed 
	   ) ; 
  if ( dc->reada==1 )  
    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->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 ; 

  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 = dc->M + dc->N ; 
    break ; 
  default:
    fprintf ( stderr , "invalid MNC = %d\n" , dc->MNC ) ; status -- ; 
    break ; 
  } 

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

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

/* space for some of the optimizer's stuff */
  p->x = dvector ( 1 , p->N ) ; 
/* this is the real parameter 
   vector of the problem 
   don't confuse with vec->x */
  p->g = dvector ( 1 , p->M ) ; 
  p->bias = dvector ( 1 , p->N ) ; 
  p->s = vec->x ; 
/* where the true vector lives */
  p->so = cvector ( 1 , p->N ) ;
/* where the answer is spat */

  return status ; 
}

int check_alist_p ( fe_min_param *p ) {
  int status = 0 ; 
  if (
      p->N != p->a.N 
      ||       p->M != p->a.M
      ) {
    fprintf ( stderr , "eek %d %d %d %d\n" , p->N , p->a.N , p->M , p->a.M ) ; 
    status -- ; 
  }
  return status ; 
}


static int score ( fe_min_param *p , fe_min_control *c , 
		  data_creation_param *dc , 
       		  mnc_vectors  *v , mnc_code *code 
) /* compares so with s and works out quality of solution if wrong. */
{
  int status = 0 ; 
  int n , nmax , count , m ; 
  FILE *fp ;

  if ( dc->verbose ) printf ( "Scoring answer --- " ) ; 
  p->count = 0 ; 
  nmax = ( dc->MNC == 4 ) ? p->N : dc->N ; 
  for ( n = 1 ; n <= nmax ; n++ ) { /* toronto: this now works out total err */
    if ( p->so[n] != v->x[n] ) {    /* not just the signal bits */
      p->count ++ ; 
    }
  }
  if ( p->count > 0 ) {
    if ( dc->verbose ) printf ( "fail\n" ) ; 
    dc->failcount ++ ; 

    /* toronto - count violations */
    alist_times_cvector_mod2 ( &(p->a) , v->x , p->t ) ; 
    for ( count = 0 , m = 1 ; m <= p->M ; m ++ ) {
      if ( p->t[m] != v->z[m] ) count ++ ;
    }
    p->viols = count ; 
    p->v = (double) (p->count_s - p->count_high) ; /* if viols = 0 and v > 0
						      then the recon is better */

/* toronto */
    if ( p->count && 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 , c , dc , p ) ;
	fprintf ( fp , "# " ) ;
	count = 0 ; 
	for ( n = 1 ; n <= p->N && ( ( count <= dc->maxcount ) 
				    || dc->maxcount < 0 ) ; n++ ) { 
	  if ( p->so[n] != v->x[n] ) {    
	    count ++ ;
	    if ( p->so[n] ) {
	      fprintf ( fp , "+%-5d " , n ) ; /* + indicates extra bit given in so */
	    } else {
	      fprintf ( fp , "-%-5d " , n ) ; /* - indicates bit lacking     */
	    }
	  }
	}
	fprintf ( fp , "\n" ) ; 
	fclose (fp ) ; 
      }
    }
  
  }  else {
    p->v = 0.0 ; 
    p->viols = 0 ; 
    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 )
{
  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 ; 
}


/* this routine should have the same effect as make_vectors followed
   by make_pvectors */

int make_vectors_quick ( data_creation_param *dc , fe_min_param *p , 
			        mnc_vectors *v )
{
  int status = 0 ; 
  int m,n , u = 1 ; 
  double lfn1fn = log ( dc->fn / ( 1 - dc->fn ) ) ;
  double lfs1fs = log ( dc->fs / ( 1 - dc->fs ) ) ;

  int count_s = 0 , count_n = 0 ; 

  switch ( dc->MNC ) {
  case ( 1) : 
    for ( n = 1 ; n <= dc->N ; n++ ) {
      /* NB, p->s[n] = v->x[n] is automatic */
      v->x[n] =     ( ranu() < dc->fs ) ? 1 : 0 ; 
    }
    for ( m = 1 ; m <= dc->M ; m++ , u ++ ) {
      v->y[m] = ( ranu() < dc->fs ) ? 1 : 0 ; 
    }
    p->true_s = 1 ; 
    break ;
  case ( 2) : 
    fprintf ( stderr, "unchecked code in use\n");
    p->true_s = 0 ; 
    break ;
  case(3): 
  case(4): 
    for ( n = 1 ; n <= dc->N ; n++ , u ++ ) {
      v->x[u] = ( ranu() < dc->fs ) ? 1 : 0 ; 
      count_s += v->x[u] ;
    }
    for ( m = 1 ; m <= dc->M ; m++ , u++) {
      v->x[u] = ( ranu() < dc->fn ) ? 1 : 0  ;