fed6.c

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

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

   - 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

    ===========
    SUMMARY
    ===========
    Input : in p->g and p->bias
    Output : in p->so, which must be allocated elsewhere.
*/

#include "./ansi/r.h"
#include "./ansi/rand2.h" 
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./ansi/macopt.h"
#include "./ansi/strangeopt.h"
#include "./fe6.h"

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

  p->x = dvector ( 1 , N ) ;  /* state of fe algm */
/* this is the real parameter 
   vector of the problem 
   don't confuse with vec->x */

  p->qd = dvector ( 1 , N ) ; 
/*  if ( c->metric ) p->q0q1 = dvector ( 1 , N ) ;  */

  if ( c ) { /* 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 ; 

  p->g = dvector ( 1 , M ) ;  /* likelihood contributions */
  p->bias = dvector ( 1 , N ) ;  /* prior contributinos */
  p->a.norder = ivector ( 1, N ) ; 
/*   p->so must be allocated and pointed to by an outside authority - 
     where the answer is spat */

}

void fe_defaults_c ( fe_min_control *c ) {
#include "fe_var6_def.c"
  macopt_defaults ( &(c->macarg) ) ; 
}

void fe_defaults_p ( fe_min_param *p ) {

  p->h_dev = 0.0 ; 
  p->h_bar = 0.0 ; 
  p->h_var = 0.0 ; 
  p->beta = 1.0 ; 
}

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

  free_dvector ( p->x , 1 , N ) ;
  free_dvector ( p->qd , 1 , N ) ; 

  if ( c ) { 
    free_dmatrix ( p->mpd , 1 , M , 0 , p->a.biggest_num_m + 1 ) ;
  }
  free_dvector ( p->pd , 0 , p->a.biggest_num_m + 1 ) ; 

  free_dvector ( p->g , 1 , M ) ; 
  free_dvector ( p->bias   , 1 , N ) ; 
  free_ivector ( p->a.norder , 1, N ) ; 
}

int fe_min ( fe_min_param *p , fe_min_control *c ) { /* returns number of diffs */
  int N = p->N ;
  int n ; 
  double *x , v ; 
  int no_sol ; 

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

  x = p->x ;

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

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

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

/*   Main loop        loop NL times    */

  no_sol = 1 ; 
  for ( c->nl = 1 ; no_sol && c->nl <= c->NL ; c->nl ++ ) {
    if( c->NL > 1 && c->verbose > 0 ) 
      printf ( "Loop %d of %d ; tol %5g ; beta %5g\n" , 
	      c->nl , c->NL , c->ftol , p->beta ) ;

    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 ) )  no_sol = 0 ; /* ensure stop when 
					     any valid solution is found */
    if (( c->verbose >= 1 ) || c->DEMO ||  ( c->nl == c->NL ) || (!no_sol) ) {
      if ( c->opt ==5 || c->opt == 6 ) c->iter = c->macarg.its ; 
      printf( " diffs %3d  true %3d  recon %3d   its %5d nl=%d\n", p->count , p->count_s , p->count_high , c->iter , c->nl ) ; 
    }

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

  return (p->count) ; 
}

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

void find_qd_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_qd ( x[n] , &(p->qd[n]) , sflag ) ;
  }
  if ( sflag ) p->S = tmpd ; 
}

double objective ( double *x , void *f_args ) 
{
  double tmpd = 0.0 ; 
  int n , N ; 
  fe_min_param *p = ( fe_min_param * ) ( f_args ) ; 

  N = p->N ;

  find_qd_S ( x , 1 , p ) ; /* puts S in S */