cg.c

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

/* 
   cg.c                                      (c) DJCM 94 07 04

   - free energy minimization for inferring state of shift register - 

   This program first generates fake data, i.e. a random A and s, and 
   a noisy version of t = A s mod 2. 

          You get to choose (on command line): 
	  N (state length),                                      (e.g. 100)
	  M (number of observations),                            (e.g. 1000)
	  diagno (number of rows of A that are singlet)          (e.g. 100)
	  fA (fraction of remaining elements of A that are high) (e.g. 0.3)
	  err (probability of error in observation of t = As)    (e.g. 0.3)
	  seed (for random numbers)

    An optimizer then runs to try to infer s. 

          optimizers available are:
	  frp (standard conjugate gradients algm)
	  macopt (my possibly-faster cg algm)
	  jumpopt (proceeds by re-estimation formula instead of 
	           creeping downhill) 
	  seqopt (sequential re-estimation so that decrease in F is 
	          guaranteed)
	  All of these optimizers can optionally be run at temperatures
	  other than beta = 1.  

	  The jumpopt can optionally omit the diagonal entries 
	  from the objective function calculations after one iteration. 

    A summary report is written. There are three possible outcomes:

          1  answer is correct
          0  answer is `wrong', and has higher energy than true state
	  2  `ok' answer is wrong, but has lower energy, so is (conditionally)
	       a perfectly good answer to an ill-posed problem. 

	       These are distinguished by the values of count (0 if right)
	       and v (negative if `ok')
	       
     05 07 94     corrected temperature-dependence, 
                  putting it into the F = S - beta E
		  instead of having it in the activation 
		  function 1/[1+exp(-beta x)]

		  NB, both S and E have sign opposite to S and E in the fe.tex
		  paper of 07 07 94

     19 07 94    monte carlo addition: 
                 after optimization, optionally go hopping around for a load of iterns
		 changing state and checking for change in energy. 
		 Just need to count how many relations a bit participates in and 
		 how many of those are satisfied. 

      01 08 94   Modified data creation to omit singlets.
*/

#include "./ansi/r.h"
#include "./ansi/rand.h"
#include "./ansi/mynr.h"


typedef struct {
  int n;
  double *m;
} param;

static int    process_command ( int , char ** ) ;
static void jump_opt (double * ) ; /* optimization by `reestimation' */
static void seq_opt (double * , param * ) ; /* optimization by `reestimation' */
static void   print_usage ( char ** , FILE * );
static void make_sense ( void ) ; 
static void print_state ( double , double , double * )  ;
static void printall ( FILE * )  ;
static void finalline ( int , double , FILE * )  ;
static void forward_pass ( void ) ;
static void backward_pass ( double * ) ;
void mega_forward_pass ( void ) ;
void mega_backward_step ( int );  /* nn should start at N and go backwards */
static void find_q_S ( double * , int ) ;
static double x_to_q ( double , double * , double * , int )  ;
static double objective ( double *, void *) ;
static void   d_objective ( double *, double *, void *) ;
void   main ( int , char ** ) ;

/* making the data */ 

static int n , N ; /* runs over state bits */
static int m , M ; /* runs over the data */
static int **A ;   /* binary matrix A[m][n] */
static int *s ;    /* the true binary state vector s[n] */
static int sordered = 0 ; /* makes an s of the form 111111110000000000 */
static int *t ;    /* the binary vector t[m] = A s mod 2 */
static double *er ;   /* the error probability vector er[m] */
static int *to ;   /* the observed binary vector 
		      to[m] = t[m] with prob 1-er[m] */
static double *g ;  /* log odds g[m] = +/- log er[m]/(1-er[m]) depending
		       whether to[m] = 1 or 0 . */
static int diagno = -1 ; /* number of rows of A that are of the 
			    form 0,0,0,1,0..0 
			    Set this to zero to get zero. Or to -1 to get 
			    the max number, N (gets reset automatically to N */
static int omit_diagonal = 0 ; /* whether to switch off the diagonal influences
				  after first iteration */
static double fA = 0.5 ;  /* fraction of bits in rest of A that are high */
static double true_fA ;   /* actual fraction */
static double true_err ;   /* actual fraction */
static double fs = 0.5 ;  /* fraction of bits in s that are high */ 
static double err = 0.25 ;       /* default error probability */

/* the model */

/* static double *x ;  x[n=1..N] 
		      the real log probs (also known as theta in my paper ) 
		      This one is not global because it gets created by the 
		      optimizer */
static double sdx = 0.0 ; /* initial standard deviation */
static double *q0 , *q1 ; /* q1 = 1 / 1+ exp (- x) , etc. */
static double *p0 , *p1 ; /* probabilities used in the forward and backward
			     passes, p0[0..N+1], p1[0..N+1] */
static double **mp0 , **mp1 ; /* mega-probabilities mp0[1..M][0..N+1] */
static int *so ; /* the answer so[1..N] */
double S , E ; /* entropy and energy of solution */

/* control */

static double beta = 1.0 ;  /* the temperature of the party */
static double beta1 = 1.0 ; 
static double betaf = 0.0 ;
static int betastyle = 0 ;  
static int     verbose = 0 ;
static int CG=0; /* CG = 1 causes gradient to be checked */
static int nl , NL = 3 ; /* Number of loops. (annealing happens within these) */
static int rich = 0 ; 
static int opt = 2 ; 
static int iter , itmax = 100 ; 

static  char 	outfile[50] = "_out" ;
static  int     printout = 0 ;

static  long int seed = 123 ;
static  double  ftol=0.0001;
static  double  flintol=0.001;
static  double  epsilon = 0.0001 ;

/* sequential stuff */
static int seq_period = 10 ; /* how frequently the free energy is evaluated 
				to decide whether to 
				stop the sequential optimizer */

/* monte carlo bit */
static int MC = 0 ; /* whether, and how many MC loops to do */
static int verbosemc = 0 ; 
static int mcm = 0 ; /* whether to do mc metropolis or mc gibbs (0) */
static int mcr = 0 ; /* whether to select candidate bit randomly or in sequence (0) */
static int mcstop = 0 ; /* whether to stop mc when count = 0 . Also whether
			   to enter mc even when count = 0 */
/* if mcstop == 1 then a stop can also occur if no flips take place in a period
   longer than mcboredom times the number of bits */
static int mcboredom = 2 , lastc = 0 ;
static  int     printmc = 0 ;
static  char 	mcfile[50] = "_mc" ;



void main(int argc, char *argv[])
{
  param control;  /*  ignore this, it's for future structured program */
  double *x , v ;
  int count , count_to , total_a_count , err_count ; 
  double tmpd ; 
  FILE *fp; 

  double deltae , field ;
  int nn = 1 ; 
/*  char 	junk[100] ; */

  M = 20  ;
  N = 10  ;

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

  printf("cg M = %d, N = %d, fA = %6.3g, err = %6.3g, opt = %d, bstyle = %d; MC = %d\n",
	 M , N , fA , err , opt , betastyle , MC ) ;

  /* Set up memory , etc. */

  A  = imatrix ( 1 , M , 1 , N ) ; 
  s  = ivector ( 1 , N ) ; 
  so = ivector ( 1 , N ) ; 
  t  = ivector ( 1 , M ) ; 
  to = ivector ( 1 , M ) ; 
  er = dvector ( 1 , M ) ; 
  g  = dvector ( 1 , M ) ; 
  x  = dvector ( 1 , N ) ; 
  q0 = dvector ( 1 , N ) ; 
  q1 = dvector ( 1 , N ) ; 
  if ( opt == 3 ) { /* mega state vectors needed */
    mp0 = dmatrix ( 1 , M , 0 , N + 1 ) ;
    mp1 = dmatrix ( 1 , M , 0 , N + 1 ) ;
    for ( m = 1 ; m <= M ; m ++ ) {
      mp0[m][0] = 1.0 ; mp1[m][0] = 0.0 ; 
      mp0[m][N+1] = 1.0 ; mp1[m][N+1] = 0.0 ; 
    }
  }
  p0 = dvector ( 0 , N + 1 ) ; 
  p1 = dvector ( 0 , N + 1 ) ; 
  p0[0] = 1.0 ; p1[0] = 0.0 ; 
  p0[N+1] = 1.0 ; p1[N+1] = 0.0 ; 

  /* make data */

  ran_seed ( seed ) ; 
  count = 0 ; 
  for ( n = 1 ; n <= N ; n++ ) 
    count += (     s[n] = ( ranu() < fs ) ? 1 : 0  ) ; 
  if ( sordered ) {
    for ( n = 1 ; n <= count ; n++ ) 
      s[n] = 1 ;
    for (  ; n <= N ; n++ ) 
      s[n] = 0 ;
  }
  for ( m = 1 ; m <= diagno ; m++ ) {
    for ( n = 1 ; n <= N ; n++ ) {
      A[m][n] = ( n == m ) ? 1 : 0 ;
    }
  }
  total_a_count = 0 ; 
  for (  ; m <= M ; m++ ) {
    count = 0 ; 
    for ( n = 1 ; n <= N ; n++ ) {
      A[m][n] = ( ranu() < fA ) ? 1 : 0 ; 
      count += A[m][n] ;
    }
    if ( count <= 1 ) m-- ; 
    else total_a_count += count ;
  }
  /* Here could print actual fA, or modify actual fA */
  true_fA = (double) total_a_count / (double) ( ( M - diagno ) * N ) ; 
  err_count = 0 ; 
  for ( m = 1 ; m <= M ; m++ ) {
    t[m] = 0 ; 
    for ( n = 1 ; n <= N ; n++ ) {
      if ( A[m][n] && s[n] ) 
	t[m] ++ ;
    }
    t[m] = t[m] % 2 ;
    if  ( ranu() < err ) {
      to[m] = 1 - t[m] ; 
      err_count ++ ; 
    } else { 
      to[m] = t[m] ; 
    }
    g[m] = log ( err / ( 1 - err ) ) ; 
    if ( to[m] ) g[m] = - g[m] ; /* log odds contributed by observation */
  }
  true_err = (double) err_count / (double) M ; 

/* set up model */ 

  for  ( n = 1 ; n <= N ; n++ ) { /* initial condition */
    x[n] = sdx * rann() ; 
  }
  if ( verbose > 1 ) printall ( stdout ) ; 


  if( CG )
    checkgrad ( x , N , epsilon , objective , &control , 
	       d_objective , &control);

/*   Main loop        loop NL times, doing a conjugate minimization     */

  for ( nl = 1 ; nl <= NL ; nl ++ ) {
    if( NL > 1 && verbose > 0 ) 
      printf ( "Loop %d of %d ; tol %5g ; beta %5g\n" , 
	      nl , NL , ftol , beta ) ;
    if ( opt == 0 ) 
      frprmn ( x , N , ftol , flintol , &iter , itmax , &v ,
	     objective , &control , d_objective , &control ) ;
    else if ( opt == 1 ) 
      macopt ( x , N , rich , ftol , &iter , itmax , 
	     d_objective , &control ) ;
    else if ( opt == 2 ) 
      jump_opt ( x ) ;
    else 
      seq_opt ( x , &control ) ;

    if( CG )
      checkgrad ( x , N , epsilon , objective , &control , 
		 d_objective , &control);

    if ( betastyle == 1 ) 	
      beta += betaf ;
    else if ( betastyle == 2 )
      beta *= betaf ; 

  }

/* End of main loop */

/* produce answer */

  count = 0 ; 
  count_to = 0 ; 
  for ( n = 1 ; n <= N ; n++ )  {
    so[n] = ( x[n] > 0 ) ? 1 : 0 ;
    if ( so[n] != s[n] ) {
      count ++ ;
    }
    if ( n <= diagno && ( so[n] != to[n] ) ) count_to ++ ;
  }
  if ( verbose > 0 ) {
    printf ( "\n          " ) ;
    printoutivector ( s , 1 , N ) ; 
  }
  if ( verbose > 0 ) 
    printf ( "ERRORS: %d; disagreements with immediate cues %d / %d\n" , 
	    count , count_to , diagno ) ; 
  if ( count ) {
    beta = 1.0 ; /* Return beta to standard value */
    for ( n = 1 ; n <= N ; n++ )  {
      x[n] = ( x[n] > 0 ) ?  100 : -100 ;
    }
    v = objective ( x , &control ) ;
    if ( verbose > 0 ) 
      printf ( "score: %6g " , v ) ; 
    for ( n = 1 ; n <= N ; n++ )  {
      x[n] = ( s[n] > 0 ) ?  100 : -100 ;
    }
    tmpd = objective ( x , &control ) ;
    if ( verbose > 0 ) 
      printf ( "c.f. correct state's score: %6g", tmpd ) ;
    if ( verbose > 0 )  {
      if ( v < tmpd ) {
	printf ( " THIS STATE BETTER\n" ) ; 
      } else {
	printf ( " \n" ) ;
      }
    }
    v -= tmpd ; /* v contains myscore - correctscore ; if negative or zero, good. */
  }
  else {
    if ( verbose > 0 ) printf ("CORRECT!\n" ) ; 
    v = 0.0 ;
  }
  if ( verbose > 0 ) finalline ( count , v , stdout ) ;
  if ( printout ) {
    fp = fopen ( outfile , "a" ) ;
    if( !fp )   fprintf( stderr, "No such file: %s\n", outfile ) ;
    else {
      finalline ( count , v , fp ) ;
      fclose ( fp ) ;
    }
  }

  /* That was the end of the show. Now for some monte carlo 
     --- start from initial state given by so */
  if ( MC && ( count || (mcstop==0) ) ) {  /* mcstop means stop when count is 0 */
    if ( verbosemc >= 0 ) printf ( "MC..." ) ;
    nn = 0 ; 
    for ( nl = 1 ; nl <= MC ; nl ++ ) {
      /* Each loop, randomly select a bit, or increment */
      if ( mcr ) nn = rani ( N ) + 1 ; 
      else { nn ++ ; if ( nn > N ) nn = 1 ; }
      
      if ( verbosemc == 2 ) printf ( "[%d] " , nn ) ; 
      /* See how many relations it is involved in and how many are not satisfied */
      for ( m = 1 , field = 0.0 ; m <= M ; m++ ) {
	if ( A[m][nn] ) {
	  t[m] = 0 ; 
	  for ( n = 1 ; n <= N ; n++ ) {
	    if ( A[m][n] && so[n] ) 
		t[m] ++ ;
	  }
	  t[m] = t[m] % 2 ;
	  if ( t[m] ) field += g[m] ;
	  else field -= g[m] ; 
	}
      }
      /*     field *= 0.5 ; */
      deltae =  ( field ) ; /* * ( so[nn] ? 1.0 : -1.0 ) ; */
      if ( ( mcm && ( (deltae < 0.0) || (ranu() < exp(-deltae)) ) )  /* Metropolis rule */ 
	  || ( !mcm &&  ( ranu() < 1.0 / ( 1.0 + exp(deltae) ) ) ) )  {  /* Gibbs rule */
	so[nn] = 1 - so[nn] ;                      
	for ( count = 0 , n = 1 ; n <= N ; n++ )  {
	  if ( so[n] != s[n] ) {
	    count ++ ;
	  }
	}
	if ( verbosemc )  { 
	  printf ("*") ;
	  if ( count == 0 ) printf ("YIPPEE\n") ; 
	}
	lastc = 1 ; 
      }
      else lastc ++ ; /* keeps count of when the last change occurred */
      if ( verbosemc ) printf ("(%6.4g) %d ", deltae , count ) ; 
      if ( verbosemc && ( ( !(nl%7)&&mcr ) || ( !mcr&& ( nn == N || !(nn%8) ) ) ) )
	  printf ("\n") ;
      if ( ( count == 0 && mcstop ) || ( mcstop == 1 && lastc > mcboredom )  ) break ; 
      /* end of MC loop */
    }
    
    if ( verbosemc >= 0 )    printf ("MC ended after %d its\n" , nl ) ;
    if ( printmc ) {
      fp = fopen ( mcfile , "a" ) ;
      if( !fp )   fprintf( stderr, "No such file: %s\n", mcfile ) ;
      else {
	finalline ( count , (double) (nl) , fp ) ;
	fclose ( fp ) ;
      }
    }
  }
}

void finalline ( int count , double v , FILE *fp ) 
{
  fprintf ( fp , "%2d %7.4g %3d %3d %7.4g %7.4g %d %ld %7.4g %7.4g\n" ,
                  count , v , N , M , fA , err , betastyle , seed , true_fA , true_err ) ;
}

double x_to_q ( double x , double *qq0 , double *qq1 , 
	       int sflag ) /* if sflag = 1, returns the binary entropy */
{
  double tmpd , l0 , l1 ;

/*  x *= beta ;  temperature was at first implemented here, but 
                 for the gradients it is better to put it elsewhere */

  if ( x > 0.0 ) {
    tmpd = 1.0 + exp ( - x ) ;
    if ( sflag ) {
      l1 = - log ( tmpd ) ; 
      l0 = -x + l1 ; 
    }
    *qq1 = 1.0/tmpd ; *qq0 = 1 - *qq1 ; 
  } else {
    tmpd = 1.0 + exp ( x ) ;
    if ( sflag ) {
      l0 = - log ( tmpd ) ; 
      l1 = x + l0 ; 
    }
    *qq0 = 1.0/tmpd ; *qq1 = 1 - *qq0 ; 
  }
  if ( sflag ) {