macopt.c

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

      }
    }
  }
  if ( actual_itmax > 1 )
    fprintf(stderr,"Reached iteration limit (%d) in macopt; continuing.\n",a->itmax); 
  macopt_free ( a ) ;	
  return;  /* (0/-1) ; abnormal end caused by time limit */
} /* NB this leaves the best value of p in the p vector, but
     the function has not been evaluated there if rich=0     */

/* maclinmin.
   Method: 
       evaluate gradient at a sequence of points and calculate the inner 
       product with the line search direction. Continue until a 
       bracketing is achieved ( i.e a change in sign ).
Note, for a 1-d optimization, this may be disatisfactory: perhaps
a simple bracketing and linear interpolation is not accurate enough.
I wrote it this way assuming we are doing a lengthy high-d optimization
and this would be accurate enough.
*/
double maclinminII 
(
 double *p , 
 void   (*dfunc)(double *,double *, void *), /* evaluates the gradient */
 void   *arg ,
 macopt_args *a )
{
  int n = a->n ; 

  double x , y ;
  double s , t , m ;
  int    its = 1 , i ;
  double step , tmpd ; 
  double  *gx = a->gx , *gy = a->gy , *met = a->m  , *xi = a->xi ;

  if ( a->verbose >= 2 ) {
    fprintf (stderr, "doing line search in direction \n" ) ;
    for ( i = 1 ; i <= n ; i ++ ) fprintf (stderr, "%8.3g " , xi [i] ) ;
    fprintf (stderr, "\n" ) ;
  }
  /* at x=0, the gradient (uphill) satisfies s < 0 */
  if ( a->verbose > 2 ) { /* check this is true: (no need to do this really
			   as it is already checked at the end of the main
			   loop of macopt) */
/*
#define TESTS 5
    x = a->lastx / a->gtyp ;
    fprintf (stderr, "inner product at:\n" ) ; 
    for ( i = -TESTS ; i <= TESTS ; i += 2 ) {
      step = x * 2.0 * (double) i / (double) TESTS ; 
      fprintf (stderr, "%9.5g %9.5g\n" , step ,
	       tmpd = macprodII ( p , gy , step , dfunc , arg , a ) ) ; 
    }
*/
    fprintf (stderr, "inner product at 0 = %9.4g\n" ,
	     tmpd = macprodII ( p , gy , 0.0 , dfunc , arg , a ) ) ; 
    if ( tmpd > 0.0 ) { 
      a->restart = 1 ; 
      return 0.0 ; 
    }
  }

  x = a->lastx / a->gtyp ;
  s = macprodII ( p , gx , x , dfunc , arg , a ) ; 
  
  if ( s < 0 )  {  /* we need to go further */
    do {
      y = x * a->linmin_g1 ;
      t = macprodII ( p , gy , y , dfunc , arg , a ) ; 
      if ( a->verbose > 1 ) 
	printf ("s = %6.3g: t = %6.3g; x = %6.3g y = %6.3g\n",s, t , x , y );
      if ( t >= 0.0 ) break ;
      x = y ; s = t ; a->gunused = gx ; gx = gy ; gy = a->gunused ; 
      its++ ;
    }
    while ( its <= a->linmin_maxits ) ;
  } else if ( s > 0 ) { /* need to step back inside interval */
    do {
      y = x * a->linmin_g3 ;
      t = macprodII ( p , gy , y , dfunc , arg , a ) ; 
      if ( a->verbose > 1 ) 
	printf ("s = %6.3g: t = %6.3g; x = %6.3g y = %6.3g\n",s, t , x , y );
      if ( t <= 0.0 ) break ;
      x = y ; s = t ; a->gunused = gx ; gx = gy ; gy = a->gunused ; 
      its ++ ;
    } while ( its <= a->linmin_maxits ) ;
  } else { /* hole in one s = 0.0 */
    t = 1.0 ; y = x;
  }

  if ( its > a->linmin_maxits )  {
    fprintf (stderr, "Warning! maclinmin overran" );
/* this can happen where the function goes \_ and doesn't buck up
 again; it also happens if the initial `gradient' does not satisfy
 gradient.`gradient' > 0, so that there is no minimum in the supposed
 downhill direction.  I don't know if this actually happens... If it
 does then I guess a->rich should be 1.

 If the overrun is because too big a step was taken then
 the interpolation should be made between zero and the most 
 recent measurement. 

 If the overrun is because too small a step was taken then 
 the best place to go is the most distant point. 
 I will assume that this doesn't happen for the moment.

 Also need to check up what happens to t and s in the case of overrun.
 And gx and gy. 

 Maybe sort this out when writing a macopt that makes use of the gradient 
 at zero? 
*/
    fprintf (stderr, "- inner product at 0 = %9.4g\n" ,
	     tmpd = macprodII ( p , gy , 0.0 , dfunc , arg , a ) ) ; 
    if ( tmpd > 0 && a->rich == 0 ) {
      fprintf (stderr, "setting rich to 1\n" ) ;       a->rich = 1 ; 
    }
    if ( tmpd > 0 ) a->restart = 1 ; 
  }

 /*  Linear interpolate between the last two. 
     This assumes that x and y do bracket. */
  if ( s < 0.0 ) s = - s ;
  if ( t < 0.0 ) t = - t ;
  m = ( s + t ) ;
  s /= m ; t /= m ;
  
  m =  s * y + t * x ; 
  /* evaluate the step length, not that it necessarily means anything */
  for ( step = 0.0 , i = 1 ; i <= n ; i ++ ) {
    tmpd = m * a->xi[i] ;
    p[i] += tmpd ; /* this is the point where the parameter vector steps */
    a->xi[i] = s * gy[i] + t * gx[i] ;
    if ( a->metric ) { /* macoptIIc : covariant measure of step length */
      step +=  tmpd * tmpd / met[i] ; 
      a->mg[i] = - met[i] * a->xi[i] ;
    }    else {
      step += fabs ( tmpd ) ;
    }
/* send back the estimated gradient in xi (NB not like linmin) */
  }
  a->lastx = m * a->linmin_g2 *  a->gtyp ;
  step =  step / (double) ( n )  ;
  if ( a->metric ) { /* macoptIIc step length */
    step = sqrt(step) ; 
  }
  return ( step ) ; 
}

double macprodII 
( 
 double *p , double *gy , double y , 
 void   (*dfunc)(double *,double *, void *), 
 void   *arg , 
 macopt_args *a
) {
  double *pt = a->pt ; 
  double *xi = a->xi ; 
  /* finds pt = p + y xi and gets gy there, 
				       returning gy . xi */
  int n = a->n ; 

  int i;
  double s = 0.0 ;

  for ( i = 1 ; i <= n ; i ++ ) 
    pt[i] = p[i] + y * xi[i] ;
  
  dfunc( pt , gy , arg ) ;

  for ( i = 1 ; i <= n ; i ++ ) 
    s += gy[i] * xi[i] ;

  return s ;
}
  
void macopt_defaults ( macopt_args *a ) {
  a->verbose = 2 ; /* if verbose = 1 then there is one report for each
		      line minimization.
		      if verbose = 2 then there is an additional report for
		      each step of the line minimization.
		      if verbose = 3 then extra debugging routines kick in.
		   */

  a->tol = 0.001 ;           /* Do fiddle with this */
  a->end_if_small_step = 1 ; /* Change this to 0/1 if you prefer */
  a->itmax = 100 ;           /* You may wish to change this */
  a->rich = 1 ;              /* if this is 1, then the program runs a bit slower */
  a->stepmax = 0.5 ; 

  a->grad_tol_tiny = 1e-16 ; /* Probably not worth fiddling with */
  a->step_tol_tiny = 0.0 ;   /* Probably not worth fiddling with */
  a->linmin_maxits = 20 ;    /* Probably not worth fiddling with */
  a->lastx = 0.01 ;          /* only has a transient effect      */
  a->lastx_default = 0.01 ;  /* -- defines typical distance in parameter
				space at which the line minimum is expected;
				both these should be set. the default is 
				consulted if something goes badly wrong and 
				a reset is demanded. */

  a->do_newitfunc = 0 ;     /* this says whether newitfunc is set to something */

/* don't fiddle with the following, unless you really mean it */
  a->linmin_g1 = 2.0 ; 
  a->linmin_g2 = 1.25 ; 
  a->linmin_g3 = 0.5 ; 
  a->restart = 0 ;
  a->metric = 0 ; /* whether we are doing things the macoptIIc way */
}

void macopt_allocate_metric (  macopt_args *a , int n ) {
  /* this routine illustrates how the metric should be allocated
     and set, except the 1.0s should be more interesting */
  a->m = dvector ( 1 , n ) ; /* metric for macoptIIc */
  for ( ; n >= 1 ; n -- ) a->m[n] = 1.0 ; 
  /*    free_dvector ( a->m , 1 , n ) ;   */
}
void macopt_allocate (  macopt_args *a , int n ) {
  a->n = n ; 
  a->g = dvector ( 1 , n ) ; /* vectors as in NR code */
  a->h = dvector ( 1 , n ) ; /*                       */
  a->xi = dvector ( 1 , n ) ;/*                       */
  a->pt = dvector ( 1 , n ) ; /* scratch vector for sole use of macprod */
  a->gx = dvector ( 1 , n ) ; /* scratch gradients             */
  a->gy = dvector ( 1 , n ) ; /* used by maclinmin and macprod */
  if ( a->metric ) {  /* macoptIIc */
    a->mg = dvector ( 1 , n ) ; /* natural gradient (contravariant) */
  }
}
void macopt_free ( macopt_args *a ) 
{
  int n = a->n ; 
  free_dvector ( a->xi , 1 , n ) ;
  free_dvector ( a->h  , 1 , n ) ;
  free_dvector ( a->g  , 1 , n ) ;  
  free_dvector ( a->pt , 1 , n ) ;   
  free_dvector ( a->gx , 1 , n ) ;  
  free_dvector ( a->gy , 1 , n ) ;
  if ( a->metric ) { /* macoptIIc */
    free_dvector ( a->mg , 1 , n ) ;
  }
}

void macopt_restart ( macopt_args *a , int start ) 
/* if start == 1 then this is the start of a fresh macopt, not a restart */
{
  int j , n=a->n ; 
  double *g, *h, *xi , *mg ;
  g = a->g ; mg = a->mg ;  h = a->h ;  xi = a->xi ; 

  if ( start == 0 ) a->lastx = a->lastx_default ; 
  /* it is assumed that (*dfunc)( p , xi , dfunc_arg ) ;  has happened, and setting mg = -m*xi */
  for ( j = 1 ; j <= n ; j ++ ) {
    if ( a->restart != 2 ) g[j] = -xi[j] ;
    if (  a->metric )     xi[j] = h[j] = mg[j] ;
    else     xi[j] = h[j] = g[j] ;
  }
  a->restart = 0 ; 
}

void maccheckgrad 
/* Examines objective function and d_objective function to see if 
   they agree for a step of size epsilon */
  (double *p,
   int    n,
   double epsilon,
   double (*func)(double *, void *),
   void   *func_arg,
   void   (*dfunc)(double *,double *, void *),
   void   *dfunc_arg ,
   int    stopat          /* stop at this component. If 0, do the lot. */
)
{
  int j;
  double f1;
  double *g,*h;
  double tmpp ; 
  
  h=dvector(1,n);
  g=dvector(1,n);
  f1=(*func)(p,func_arg);
  (*dfunc)(p,g,dfunc_arg);
  if ( stopat <= 0 || stopat > n ) stopat = n ; 

  printf("Testing gradient evaluation  (epsilon = %9.5g)\n", epsilon);
  printf("      analytic   1st_diffs     difference\n");
  for ( j = 1 ; j <= stopat ; j ++ ) {
    tmpp = p[j] ; 
    p[j] += epsilon ;
    h[j] = (*func)(p,func_arg) - f1 ;
    p[j] =  tmpp ;

    printf("%2d %12.5g %12.5g %12.5g\n" , j , g[j] , h[j]/epsilon , g[j] - h[j]/epsilon );
    fflush(stdout) ; 
  }
  free_dvector(h,1,n);
  free_dvector(g,1,n);
  printf("      --------     ---------\n");
}

/*****************************************************
  find second derivative of the log likelihood using
  the first-derivative function 
  ****************************************************/
void    evaluate_hessian 
( double **H , /* put the hessian here */
  double *p ,  /* point for evaluation */
  int n ,               /* number of dimensions                           */
  double epsilon ,
  void   (*dfunc)(double *,double *, void *), 
                      /* evaluates the gradient of the optimized function */
  void   *dfunc_arg,    /* arguments that get passed to dfunc             */
  int verbose )
{
  int i,j;
  double *g,*h;
  double tmpp ; 
  
  h=dvector(1,n); /* this will store the original gradient */
  g=dvector(1,n); /* and this the new one */
  (*dfunc)(p,h,dfunc_arg);

  if ( verbose >= 1  ) {
    printf("Evaluating Hessian (epsilon = %9.5g)\n", epsilon);
  }
  for ( j = 1 ; j <= n ; j ++ ) {
    tmpp = p[j] ; 
    p[j] += epsilon ;
    (*dfunc)(p,g,dfunc_arg);
    for ( i = 1 ; i <= n ; i ++ ) {
      H[i][j] = ( g[i] - h[i] )/ epsilon ; 
      if ( verbose >= 1  ) {
	printf ("%10.3g\t", H[i][j] ) ; 
      }
    }
    p[j] =  tmpp ;

    if ( verbose >= 1  ) {
      printf ("\n") ; 
      fflush(stdout) ; 
    }
  }
  free_dvector(h,1,n);
  free_dvector(g,1,n);
}

/*
<!-- hhmts start -->
Last modified: Mon Nov 24 21:18:02 1997
<!-- hhmts end -->
*/