fe.c

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

double objective ( double *x , void *f_args ) /* this could be turned into param *controlp */ 
{
  double tmpd = 0.0 ; 
  int n , N ; 
  fe_min_param *p = ( fe_min_param * ) ( f_args ) ; 

  N = p->N ;

  find_q_S ( x , 1 , p ) ; /* puts S in S */
  p->EL = 0.0 ; 
  p->EP = 0.0 ; 

  if ( p->betap != 0.0 ) {
    for ( n = 1 ; n <= N ; n++ ) {
      p->EP += p->bias[n] * p->q1[n] ; 
    }
  }

  if ( p->beta != 0.0 ) {
    for ( p->m = 1 ; p->m <= p->M ; p->m ++ ) {
      p->EL += forward_pass ( p ) * p->g[p->m] ; 
    }
  }

  p->F = tmpd  = 
    ( p->zero_temperature ? - p->EL :
     p->S - p->beta * p->EL - p->betap * p->EP ) ; 
  if ( p->c->verbose > 0 ) 
    printf ( "%6.4g " , tmpd ) ; 
  return ( tmpd ) ;
}

void new_forward_pass ( fe_min_param *p ) /* m is passed inside p */
{ /* differs in that p->EL is incremented internally
     and other forms of E_L can be handled too */
  int l , n  , m = p->m ; 
  int lmax = p->a.num_mlist[m] ;
  int *mlist = p->a.mlist[m] ; 
  double *q0 = p->q0 ; 
  double *q1 = p->q1 ; 
  double *p0 = p->p0 ; 
  double *p1 = p->p1 ; 
  double ans ; 
  double dgamma = p->c->dgamma ;
#define EPSILON 1e-5 
#define RECTIFY; if ( ans < EPSILON ) { ans = EPSILON ; }
/* was: fprintf ( stdout , "" ) ;    ans < EPSILON in new_forward_pass */

  for ( l = 1 ; l <= lmax ; l ++ ) {
    n = mlist[l] ; 

    p0[l] = q0[n] * p0[l-1] + q1[n] * p1[l-1] ;
    p1[l] = q0[n] * p1[l-1] + q1[n] * p0[l-1] ;
  }

/* compute contribution to E_L; NB, EL is back to front in sign, by 
   an absurd convention of this program */
  if ( p->c->gamma == -1 ) { /* traditional answer */
    p->EL +=  p1[lmax] * p->g[m] ;
  } else {
    /* look at the vector `z' to decide what the score is. 
       And Load up g[m] with the appropriate value for when the 
       derivative is calculated. 
       */
    if ( p->z[m] ) {
      ans = p1[lmax] ; RECTIFY; 
      p->g[m] =   pow ( ans , - dgamma - 1.0 ) ; 
    } else {
      ans = p0[lmax] ; RECTIFY; 
      p->g[m] = - pow ( ans , - dgamma - 1.0 ) ; 
    }
    p->EL += ( ( p->c->gamma > 0 ) ?
	      ( - pow ( ans , - dgamma ) / dgamma ) :
	      ( log ( ans ) ) ) ; 
  }
  return ; 
}  

double forward_pass ( fe_min_param *p ) /* m is passed inside p */
{
  int l , n  , m = p->m ; 
  int lmax = p->a.num_mlist[m] ;
  int *mlist = p->a.mlist[m] ; 
  double *q0 = p->q0 ; 
  double *q1 = p->q1 ; 
  double *p0 = p->p0 ; 
  double *p1 = p->p1 ; 

  for ( l = 1 ; l <= lmax ; l ++ ) {
    n = mlist[l] ; 

    p0[l] = q0[n] * p0[l-1] + q1[n] * p1[l-1] ;
    p1[l] = q0[n] * p1[l-1] + q1[n] * p0[l-1] ;
  }
  return p1[lmax] ; 
}

void mega_forward_pass ( fe_min_param *p ) 
{
  int l , ll , u ; 
  int m , n ;
  double *p0 , *p1 , q0 , q1 , p0o , p1o ; 
  int *norder = p->a.norder ; 
  int *lupto  = p->a.l_up_to ; 
  int *nlist    ; 
  int lmax ;

  for ( m = 1 ; m <= p->M ; m++ ) 
    lupto[m] = 0 ; 

  for ( u = 1 ; u <= p->N ; u++ ) {
    n = norder[u] ; 
    lmax = p->a.num_nlist[n] ; 
    nlist = p->a.nlist[n] ; 
    q0 = p->q0[n] ;
    q1 = p->q1[n] ;
    for ( l = 1 ; l <= lmax ; l ++ ) {
      m = nlist[l] ; 
      p0 = p->mp0[m] ; 
      p1 = p->mp1[m] ; 

      ll = lupto[m] ; 
      p0o = p0[ll] ; p1o = p1[ll] ; 

      lupto[m] ++ ;
      ll ++ ; 

      p0[ll] = q0 * p0o + q1 * p1o ;
      p1[ll] = q0 * p1o + q1 * p0o ;
    }
  }
}

void backward_pass ( double *gr , fe_min_param *p ) 
{
  double tmpd ; 
  int l ; 
  int n  , m = p->m ; 

  l = p->a.num_mlist[m] ;
  p->p0[ l+1 ] = 1.0 ; p->p1[ l+1 ] = 0.0; 
  for ( ; l >= 1 ; l -- ) {
      n = p->a.mlist[m][l] ; 
    
      p->p0[l] = p->q0[n] * p->p0[l+1] + p->q1[n] * p->p1[l+1] ;
      p->p1[l] = p->q0[n] * p->p1[l+1] + p->q1[n] * p->p0[l+1] ;

      tmpd = p->p1[l-1] * p->p1[l+1] +  p->p0[l-1] * p->p0[l+1] -
	     p->p0[l-1] * p->p1[l+1] -  p->p1[l-1] * p->p0[l+1]  ;
      gr[n] -= p->g[m] * tmpd ; 
  }
}

void mega_backward_step ( int nn  , fe_min_param *p )  /* nn should go backwards 
				     through the sequence used in megaforward */
{
  int l , ll , m ; 
  double q0 = p->q0[nn] ; 
  double q1 = p->q1[nn] ; 
  double *p0 , *p1 , p0o , p1o ; 
  int lmax = p->a.num_nlist[nn] ; 
  int *nlist  = p->a.nlist[nn] ;
  int *lupto = p->a.l_up_to ; 

  for ( l = 1 ;  l <= lmax ; l ++ ) {
    m = nlist[l] ; 
    p0 = p->mp0[m] ; 
    p1 = p->mp1[m] ; 
    
    ll = lupto[m] ;
    p0o = p0[ll+1] ; 
    p1o = p1[ll+1] ; 
    
    lupto[m] -- ;
    
    p0[ll] = q0 * p0o  + q1 * p1o ;
    p1[ll] = q0 * p1o  + q1 * p0o ;
  }
}

void d_objective( double *x, double *gr, void *f_args)
{
  double tmpd ;
  tmpd = dd_objective( x, gr, f_args) ;
}

double dd_objective( double *x, double *gr, void *f_args)
{
  double tmpd , gg ; 
  int n ; 
  fe_min_param *p = ( fe_min_param * ) ( f_args ) ; 
  int metric = p->c->metric ; 

  for ( n = 1 ; n <= p->N ; n++ ) gr[n] = 0.0 ; 

  find_q_S ( x , 1 , p ) ;
  p->EL = 0.0 ; 

  for ( p->m = 1 ; p->m <= p->M ; p->m ++ ) {
    new_forward_pass ( p ) ; 

    backward_pass ( gr , p ) ; 
  }

  for ( n = 1 , gg = 0.0  ; n <= p->N ; n++ ) {
    if ( !p->zero_temperature ) {
      gr[n] *= p->beta ;			  /* temperature effect */
      gr[n] -= p->betap * p->bias[n] ;		  /* prior effect */
      gr[n] += x[n] ;                             /* derivative of entropy */
    }
    if ( metric ) gr[n] *= p->q0[n] * p->q1[n] ;  /* factor omitted from all terms */
    gg += gr[n] * gr[n] ; 
  }

  tmpd = p->F = 
    ( p->zero_temperature ? - p->EL :
     p->S - p->beta * p->EL - p->betap * p->EP ) ;
  if ( p->c->verbose ) {
    print_state ( tmpd , gg , x , p ) ; 
  }
  if ( p->c->writelog ) {
    fprint_state ( tmpd , gr , x , p ) ; 
  }
  return tmpd ; 
}

void jump_opt (double *x, fe_min_param *p, fe_min_control *c ) /* optimization by `reestimation' */
{
  double tmpd , gg ; 
  double *gr  ;
  int i = 0 , n   ; 

  gr = dvector ( 1 , p->N ) ;

  do {
    i ++ ; 
    for ( n = 1 ; n <= p->N ; n++ ) {
      gr[n] = 0.0 ; 
    }

    find_q_S ( x , 1 , p ) ;
    p->EL = 0.0 ; 

    
    for ( p->m = 1 ; p->m <= p->M ; p->m ++ ) { 
      p->EL +=     forward_pass ( p ) * p->g[p->m] ; 

      backward_pass ( gr , p ) ;  
    }

    for ( n = 1 , gg = 0.0  ; n <= p->N ; n++ ) {
      gr[n] *= p->beta ;             /* temperature effect */
      gr[n] -= p->betap * p->bias[n] ; /* prior effect */
      tmpd = gr[n] ; 
      gr[n] += x[n] ;               /* derivative of entropy */
      x[n] = - tmpd ; 
      gg += gr[n] * gr[n] ; 
    }

    tmpd = p->S - p->beta * p->EL ; 
    if ( c->verbose ) {
      print_state ( tmpd , gg , x , p ) ; 
    }
  }
  while ( gg > c->ftol && i < c->itmax ) ;
  if ( c->DEMO ) printf ("%3d Jumps, beta %6.3g. ",i,p->beta ) ; 

  free_dvector ( gr , 1 , p->N ) ; 
}


void seq_opt (double *x , fe_min_param *p , fe_min_control *c )
/* optimization by `reestimation', sequential */
/* an advantage of this approach is that there is a free energy
   improvement guarantee. Also we might expect some sort of metric
   effect-- i.e. unlike steepest descent, here, the change of x[n+1]
   will be influenced by the change just made in x[n] */
{
  double tmpd , grnn , v0 , v1 ; 
  int m ; 

  int nn , l , ll ; 
  int u ; /* runs over the norder */
/*  alist_matrix *a = controlp->a ;  */

  c->iter = 0 ; 
  v1 = objective ( x , p ) ; /* this sets up the qs for us too */
  v0 = v1 + 2 * c->ftol ; /* hack to prevent premature termination */

  do {
    mega_forward_pass ( p ) ; 
    c->iter ++ ; 
    for ( u = p->N ; u >= 1 ; u-- ) { /* reverse sequence used in mega_for */
      nn = p->a.norder[u] ; 
      grnn = 0.0 ; 
      for ( l = 1 ;  l <= p->a.num_nlist[nn] ; l ++ ) {
	m = p->a.nlist[nn][l] ; 
	ll = p->a.l_up_to[m] ;
	tmpd = p->mp1[m][ll-1] * p->mp1[m][ll+1] +  p->mp0[m][ll-1] * p->mp0[m][ll+1] -
	       p->mp0[m][ll-1] * p->mp1[m][ll+1] -  p->mp1[m][ll-1] * p->mp0[m][ll+1]  ;
	grnn -= p->g[m] * tmpd ; 
      }
      grnn *= p->beta ;
      grnn -= p->betap * p->bias[nn] ; /* prior effect */
      x[nn] = - grnn ; /* evaluate the new value of q */
      x_to_q ( x[nn] , &(p->q0[nn]) , &(p->q1[nn]) , 0 ) ;
      if ( c->seqverbose >= 2 ) printf("%4d %8.3g\n",nn,x[nn] ) ; 
      mega_backward_step ( nn , p ) ;
    }
    if ( c->seqverbose >= 2 ) printf("\n");
    if ( ! ( c->iter % c->seq_period ) ) { /* default period = 10 */
      v0 = v1 ; 
      v1 = objective ( x , p ) ;
      if ( c->seqverbose ) printf ("%8.5g %8.5g\n" , v1 , ( v0 - v1 ) ) ; 
      if ( v0 - v1 < -c->ftol  ) fprintf ( stderr , "Warning, F increased significantly %6.3g\n" , v0 - v1 ) ; 
    }
  }
  while ( v0 - v1 > c->ftol && c->iter < c->itmax ) ;
  if (c->seqverbose) 
    printf ("seq_opt done: %8.3g %3d Js, beta %6.3g\n",(v0-v1),
	    c->iter,p->beta ) ; 

}

void print_state ( double tmpd , double gg , double *x , fe_min_param *p ) {
  int n ;

  for ( n = 1 ; n <= p->N ; n++ ) printf ( "%2.0f " , 99.0 * p->q1[n] ) ; 
  printf ( ":%6.3g\nx: " , gg ) ;
  for ( n = 1 ; n <= p->N ; n++ ) printf ( "%6.3g " , x[n] ) ; 
  printf ( "\n" ) ; 
}

void fprint_state (  double tmpd , double *gr , double *x , fe_min_param *p ) {
  FILE *fp ;
  int n ;

  fp = fopen ( p->c->logfile , "a" ) ;
  for ( n = 1 ; n <= p->c->writelog ; n++ ) 
    fprintf ( fp ,  "%-2.0f " , 99.0 * p->q1[n] ) ; 
  fprintf ( fp ,  "%-9.4g " , tmpd ) ;
  for ( n = 1 ; n <= p->c->writelog ; n++ ) 
    fprintf ( fp ,  "%-9.4g " , x[n] ) ; 
  fprintf ( fp ,  "%-9.4g " , tmpd ) ;
  for ( n = 1 ; n <= p->c->writelog ; n++ ) 
    fprintf ( fp ,  "%-9.4g " , gr[n] ) ; 
  fprintf ( fp ,  "\n" ) ; 
  fclose ( fp ) ; 
}


void printall ( FILE *fp , fe_min_param *p ) 
{
  printoutcvector ( p->s , 1 , p->N ) ; 
  printf ( "\n" ) ; 
  printoutcvector ( p->t , 1 , p->M ) ; 
  pdv ( p->g , 1 , p->M , 4 ) ;   printf ( "\n" ) ; 
}

/*
<!-- hhmts start -->
Last modified: Sat May 27 19:32:46 1995
<!-- hhmts end -->
*/