fed.c

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

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

   - free energy minimization algorithm  - 

   - based on differences qd = q0 - q1, pd = p0 - p1, instead of q0 and q1 &c

   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->qd = dvector ( 1 , N ) ; 
/*  if ( c->metric ) p->q0q1 = dvector ( 1 , N ) ;  */

  if ( c->opt == 3 ) { /* mega state vectors needed */
    p->mpd = dmatrix ( 1 , M , 0 , p->a.biggest_num_m + 1 ) ;
    for ( m = 1 ; m <= M ; m ++ ) {
      p->mpd[m][0] = 1.0 ; 
      p->mpd[m][p->a.num_mlist[m]+1 ] = 1.0 ; 
    }
  }

  p->pd = dvector ( 0 , p->a.biggest_num_m + 1 ) ; 

  p->pd[0] = 1.0 ; 
}

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 ; 
  c->metric = 1 ; 

  macopt_defaults ( &(c->macarg) ) ; 

  p->h_dev = 0.0 ; 
  p->h_bar = 0.0 ; 
  p->h_var = 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->qd , 1 , N ) ; 

  if ( c->opt == 3 ) { 
    free_dmatrix ( p->mpd , 1 , M , 0 , p->a.biggest_num_m + 1 ) ;
  }
  free_dvector ( p->pd , 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 ; 

/*
#ifdef _DEBUG_MALLOC_INC
  unsigned long histid1, histid2 , orig_size , current_size ; 
#endif
*/

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

    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_qd_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_qd ( double x , double *qqd , 
	       int sflag ) /* if sflag = 1, returns the binary entropy */
{
  double tmpd , l0 , l1 , q0 , q1 ;

  if ( x > 0.0 ) {
    tmpd = 1.0 + exp ( - x ) ;
    if ( sflag ) {
      l1 = - log ( tmpd ) ; 
      l0 = -x + l1 ; 
    }
    q1 = 1.0/tmpd ; q0 = 1 - q1 ; 
  } else {
    tmpd = 1.0 + exp ( x ) ;
    if ( sflag ) {
      l0 = - log ( tmpd ) ; 
      l1 = x + l0 ; 
    }
    q0 = 1.0/tmpd ; q1 = 1 - q0 ; 
  }
  *qqd = q0 - q1 ; /* this could be improved on... */
  if ( sflag ) {
    tmpd =  ( q0 * l0 + q1 * l1 ) ; 
    return tmpd ;
  }
  else return 0.0 ; 
}