fed6.c

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

  p->EL = 0.0 ; 
  p->EP = 0.0 ; 

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

  if ( p->beta != 0.0 ) {
    for ( p->m = 1 ; p->m <= p->M ; p->m ++ ) {
      new_forward_pass ( p ) ; 
    }
  }

  p->F = tmpd  = 
    ( p->c->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 *qd = p->qd ; 
  double *ppd = p->pd ; 
  double p0 , 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] ; 

    ppd[l] = qd[n] * ppd[l-1] ;
  }
  p0 = ( ppd[lmax] + 1.0 ) * 0.5 ; 
  p1 = ( 1.0 - ppd[lmax] ) * 0.5 ; 

/* 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 * 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 ; RECTIFY; 
      p->g[m] =   pow ( ans , - dgamma - 1.0 ) ; 
    } else {
      ans = p0 ; RECTIFY; 
      p->g[m] = - pow ( ans , - dgamma - 1.0 ) ; 
    }
    p->EL += ( ( p->c->gamma > 0 ) ?
	      ( - pow ( ans , - dgamma ) / dgamma ) :
	      ( log ( ans ) ) ) ; 
  }
  return ; 
}  

void mega_forward_pass ( fe_min_param *p ) 
{
  int l , ll , u ; 
  int m , n ;
  double *ppd , qd ; 
  double **mpd = p->mpd ;
  int *norder = p->a.norder ; 
  int *lupto  = p->a.l_up_to ; 
  int *nlist    ; 
  int lmax ;
  alist_matrix *a = &(p->a) ;

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

  for ( u = 1 ; u <= p->N ; u++ ) {
    n = norder[u] ; 
    lmax = a->num_nlist[n] ; 
    nlist = a->nlist[n] ; 
    qd = p->qd[n] ;

    for ( l = 1 ; l <= lmax ; l ++ ) {
      m = nlist[l] ; 
      ppd = mpd[m] ; 

      ll = lupto[m] ; 
      lupto[m] ++ ;

      ppd[ll+1] = qd * ppd[ll] ;
    }
  }
}

void backward_pass ( double *gr , fe_min_param *p ) 
{
  int l ; 
  int n  , m = p->m ; 
  double *ppd = p->pd , *qd = p->qd , *g = p->g ;
  int  *mlistm = p->a.mlist[m] ;

  l = p->a.num_mlist[m] ;
  ppd[ l+1 ] = 1.0 ;
  for ( ; l >= 1 ; l -- ) {
      n = mlistm[l] ; 
    
      ppd[l] = qd[n] * ppd[l+1] ;

      gr[n] -= g[m] * ppd[l-1] * ppd[l+1]  ; 
  }
}

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 qd = p->qd[nn] ; 
  double *ppd ; 
  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] ; 
    ppd = p->mpd[m] ; 
    
    ll = lupto[m] ;
    lupto[m] -- ;
    
    ppd[ll] = qd * ppd[ll+1] ;
  }
}

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 ; 
  double *qd = p->qd ; 

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

  find_qd_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->c->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] *= (1.0 - qd[n]) * (qd[n] + 1.0) * 0.25 ;
    gg += gr[n] * gr[n] ; 
  }

  tmpd = p->F = 
    ( p->c->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_qd_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++ ) {
      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 grnn , v0 , v1 ; 
  int m ; 

  int nn , l , ll ; 
  int u ; /* runs over the norder */
/*  alist_matrix *a = controlp->a ;  */
  int *num_nlist = p->a.num_nlist ;
  double **mpd = p->mpd ; 
  double *g = p->g ; 
  int *l_up_to = p->a.l_up_to ;
  int **nlist = p->a.nlist ; 
  int *nlistnn ; 

  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] ; 
      nlistnn = nlist[nn] ;
      grnn = 0.0 ; 
      for ( l = 1 ;  l <= num_nlist[nn] ; l ++ ) {
	m = nlistnn[l] ; 
	ll = l_up_to[m] ;
	grnn -= g[m] * mpd[m][ll-1] * mpd[m][ll+1] ;
      }
      grnn *= p->beta ;
      grnn -= p->betap * p->bias[nn] ; /* prior effect */
      x[nn] = - grnn ; /* evaluate the new value of q */
      x_to_qd ( x[nn] , &(p->qd[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 " , 9.0 * p->qd[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 ,  "%-3.0f " , 99.0 * p->qd[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: Thu Sep 21 14:46:36 1995
<!-- hhmts end -->
*/