fe.c

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

/* fe.c                                      (c) DJCM 94 07 04
                                                      94 10 26
						      94 10 31
						      95 01 02
						      95 01 09

   - free energy minimization algorithm  - 

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

   -----
   fe.c:
   -----

   This is a set of routines to do a free energy optimization for 
 an arbitrary problem 

 A s + n = z 

 defined by 

 A, represented in an alist_matrix. 
 a prior on s, given in a vector b
 a prior on n, given in a vector g which implicitly contains z too, 
               in the sign of g. (z is passed explicitly in some cases)
 an initial condition
 
 Also if you tell fe_min the right answer it can tell you how it did. 

   The principal routine, fe_min, is also given control parameters 
   specifying which optimizer to use and what annealing procedures to use. 

   Inference:

    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)
	  macoptII (my possibly-faster cg algm)
	  strangeopt ( a gradient descent algm intended for the infinite
	               temperature case )

	  All of these optimizers can optionally be run at temperatures
	  other than beta = 1.  

    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')
	       
    Output of this program can be turned into convenient gnuplot ps files
    by running codeperfms.p outputfile | gnuplot
    The ps file is named by munging the name of outputfile


*/

#include "./ansi/r.h"
#include "./ansi/rand.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./ansi/macopt.h"
#include "./ansi/strangeopt.h"
#include "./fe.h"

void fe_allocate ( fe_min_param *p , fe_min_control *c ) {
  int M = p->M ; 
  int N = p->N ;
  int m ; 

/* toronto */
  p->t = cvector ( 1 , M ) ; 
  p->x = dvector ( 1 , N ) ; 

  p->q0 = dvector ( 1 , N ) ; 
  p->q1 = dvector ( 1 , N ) ; 

  if ( c->opt == 3 ) { /* mega state vectors needed */
    p->mp0 = dmatrix ( 1 , M , 0 , p->a.biggest_num_m + 1 ) ;
    p->mp1 = dmatrix ( 1 , M , 0 , p->a.biggest_num_m + 1 ) ;
    for ( m = 1 ; m <= M ; m ++ ) {
      p->mp0[m][0] = 1.0 ; p->mp1[m][0] = 0.0 ; 
      p->mp0[m][p->a.num_mlist[m]+1 ] = 1.0 ; p->mp1[m][p->a.num_mlist[m]+1 ] = 0.0; 
    }
  }
/* NB can't use p0 and p1 with fixed length for backward pass */
  p->p0 = dvector ( 0 , p->a.biggest_num_m + 1 ) ; 
  p->p1 = dvector ( 0 , p->a.biggest_num_m + 1 ) ; 

  p->p0[0] = 1.0 ; p->p1[0] = 0.0 ; 
/*
  p->p0[p->a.biggest_num_m + 1] = 1.0 ; p->p1[p->a.biggest_num_m+1] = 0.0 ; 
  This is handled at the beginning of backpass.
*/
}

void fe_defaults ( fe_min_param *p , fe_min_control *c ) {
#include "fe_var_def.c"
  p->beta = 1.0 ; 
  c->betap = 1.0 ; 
  p->c = c ; 

  c->init_given = 0 ; /* whether an initial condition has been 
			  given in the accompanying fe_min_param */
  c->sdx = 0.0 ;   /* initial standard deviation */
  c->beta0 = 1.0 ; 
  c->beta1 = 1.0 ; 
  c->betastyle = 0 ;  
  c->verbose = 0 ;
  c->seqverbose = 0 ;
  c->CG=0;            /* CG = 1 causes gradient to be checked */
  c->NL = 1 ;         /* Number of loops. (annealing happens within these) */
  c->rich = 0 ; 
  c->opt = 3 ;        /* NB opt must be 3 */
  c->itmax = 100 ; 
  
  c->printout = 0 ;
  c->seed = 0 ;
  c->ftol=0.0001;
  c->flintol=0.001;
  c->epsilon = 0.0001 ;
  c->seq_period = 10 ; /* how frequently the free energy is evaluated 
			   to decide whether to 
			   stop the sequential optimizer */
  
  c->pheading_count=0; 
  c->pheading_period=20; 
  c->NORDER = 0 ; 
  c->DEMO = 0 ; 

  macopt_defaults ( &(c->macarg) ) ; 

  p->h_dev = 0.0 ; 
  p->zero_temperature = 0 ; 
}

void fe_free ( fe_min_param *p , fe_min_control *c ) {
  int M = p->M ; 
  int N = p->N ;

  free_dvector ( p->q0 , 1 , N ) ; 
  free_dvector ( p->q1 , 1 , N ) ; 

  if ( c->opt == 3 ) { 
    free_dmatrix ( p->mp0 , 1 , M , 0 , p->a.biggest_num_m + 1 ) ;
    free_dmatrix ( p->mp1 , 1 , M , 0 , p->a.biggest_num_m + 1 ) ;
  }
  free_dvector ( p->p0 , 0 , p->a.biggest_num_m + 1 ) ; 
  free_dvector ( p->p1 , 0 , p->a.biggest_num_m + 1 ) ; 
/*  should free p->x , 1 , N   also? */
/* toronto */
  free_cvector (   p->t , 1 , M ) ;
}

int fe_min ( fe_min_param *p , fe_min_control *c ) {
  int N = p->N ;
  int n ; 
  double *x , v ; 
  int nl ; 

  x = p->x ;

  p->beta = c->beta0 ;
  p->betap = c->betap ;

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

  p->count_s = 0 ; 
  if ( p->true_s ) {
    for  ( n = 1 ; n <= N ; n++ ) { 
      p->count_s += p->s[n] ;
    }
  }
      
  if ( ! c->init_given ) {
    if ( c->seed != 0 ) ran_seed ( c->seed ) ; 
    for  ( n = 1 ; n <= N ; n++ ) { /* initial condition */
      x[n] = p->bias[n] + c->sdx * rann() ; 
    }
  }

  if ( c->verbose > 1 ) printall ( stdout , p ) ; 

/*
  orig_size = malloc_inuse ( &histid1 ) ; 
*/

  if( c->CG ) {
    switch ( c->opt ) {
    case ( 0 ) : case ( 1 ) : case(2): case(3): case (4): default: 
      checkgrad ( x , N , c->epsilon , objective , p , 
		 d_objective , p);
      break ; 
      
    case(5):
    case(6):
      ddcheckgrad ( x , N , c->epsilon ,  dd_objective , p  , 0 ) ; 
      break ;
    }
  }

/*
  current_size = malloc_inuse ( &histid2 ) ; 
  if ( current_size != orig_size ) {
    fprintf ( stderr , "checkgrad memory cockup\n" ) ; 
    malloc_list ( 2 , histid1 , histid2 ) ; 
  }
  else  fprintf ( stderr , "checkgrad memory OK\n" ) ; 
*/

/*   Main loop        loop NL times    */

  for ( nl = 1 ; nl <= c->NL ; nl ++ ) {
    if( c->NL > 1 && c->verbose > 0 ) 
      printf ( "Loop %d of %d ; tol %5g ; beta %5g\n" , 
	      nl , c->NL , c->ftol , p->beta ) ;
/*
    orig_size = malloc_inuse ( &histid1 ) ; 
*/
    switch ( c->opt ) {
    case ( 0 ) :
      frprmn ( x , N , c->ftol , c->flintol , &(c->iter) , c->itmax , &v ,
	     objective , p , d_objective , p ) ;
      break ; 
    case ( 1 ) :
      macopt ( x , N , c->rich , c->ftol , &(c->iter) , c->itmax , 
	     d_objective , p ) ;
      break; 
    case(2):
      jump_opt ( x , p , c ) ;
      break ; 
    case(3):
    default :
      seq_opt ( x ,  p , c) ;
      break ; 
    case (4):
      /* will put macoptII here (what is betaopt?) */
      break;
    case(5):
      strangeopt ( x , N , dd_objective , p , &(c->macarg) ) ; 
      break ; 
    case(6):
      creepopt ( x , N , dd_objective , p , &(c->macarg) ) ; 
      break ; 
    }

/*    current_size = malloc_inuse ( &histid2 ) ; 
    if ( current_size != orig_size ) {
      fprintf ( stderr , "opt memory cockup\n" ) ; 
      malloc_list ( 2 , histid1 , histid2 ) ; 
    }
    else  fprintf ( stderr , "opt memory OK\n" ) ; 
*/

    if( c->CG ) {
      switch ( c->opt ) {
      case ( 0 ) : case ( 1 ) : case(2): case(3): case (4): default: 
	checkgrad ( x , N , c->epsilon , objective , p , 
		   d_objective , p);
	break ; 
	
      case(5):
      case(6):
	ddcheckgrad ( x , N , c->epsilon ,  dd_objective , p  , 0 ) ; 
	break ;
      }
    }

/* End of main loop */

/* produce answer */

    find_q_S ( x , 0 , p ) ; 
    if ( c->verbose >= 1 ) print_state ( 0.0 , 0.0, x , p  ) ; 

    for ( p->count = 0 , p->count_high = 0 , n = 1 ; n <= N ; n++ )  {
      p->so[n] = ( x[n] > 0 ) ? 1 : 0 ;
      if ( p->true_s && (p->so[n] != p->s[n] ) ) {
	p->count ++ ;
      }
      p->count_high += p->so[n] ;
    }
    if ( c->verbose > 0 ) {
      if(p->true_s){
	printf ( "\n          " ) ;
	printoutcvector ( p->s , 1 , N ) ; 
      }
      printf ( "\n          " ) ;
      printoutcvector ( p->so , 1 , N ) ; 
    }
    if ( p->true_s && (  p->count == 0 ) )  nl = c->NL ; /* ensure stop when 
					     any valid solution is found */
    if (( c->verbose >= 1 ) || c->DEMO ||  ( nl == c->NL )) {
      if ( c->opt ==5 || c->opt == 6 ) c->iter = c->macarg.its ; 
      printf( " diffs %3d  true %3d  recon %3d   its %5d\n", p->count , p->count_s , p->count_high , c->iter ) ; 
    }

/* this belongs at end of main loop */
    if ( nl < c->NL ) {
      if ( c->betastyle == 1 ) 	
	p->beta += c->betaf ;
      else if ( c->betastyle%10 == 2 )
	p->beta *= c->betaf ; 
    }
  }

  return (1) ; 
}

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

  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 ) {
    tmpd =  ( *qq0 * l0 + *qq1 * l1 ) ; 
    return tmpd ;
  }
  else return 0.0 ; 
}

void find_q_S ( double *x , int sflag , fe_min_param *p )
{
  double tmpd = 0.0 ; 
  int n , N = p->N ;

  for ( n = 1 ; n <= N ; n++ ) {
    tmpd += x_to_q ( x[n] , &(p->q0[n]) , &(p->q1[n]) , sflag ) ;
  }
  if ( sflag ) p->S = tmpd ; 
}