ms.c

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

void beta_opt (double *x ) /* optimization by `reestimation' jumps, 
			     with the value of beta being controlled
			     by looking at magnitude of jumps and 
			     at the inner product of successive jumps*/
{
  double tmpd , gg , sum_dx , ip_dx , sdx1=0.0 , sdxp ; 
  double *gr  , *dx , *dxp , *xp ;
  int i = 0  , moratorium = 0 , all_clear ; 
#define START_M 5

  gr = dvector ( 1 , N ) ;
  dx = dvector ( 1 , N ) ;/* MONITOR */
  dxp = dvector ( 1 , N ) ;/* MONITOR */
  xp = dvector ( 1 , N ) ;/* MONITOR */
    for ( n = 1 ; n <= N ; n++ ) {
      dxp[n] = 0.0 ; /* MONITOR */
      xp[n] = x[n] ; 
    }

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

    find_q_S ( x , 1 ) ; /* zero if we don't want the entropy; 
			switch this to 1 to make a routine that makes both 
			the objective function and its derivative */
    E = 0.0 ; 

    m = 1 ;
    for ( ; m <= M ; m ++ ) { 
      forward_pass ( ) ; 
      E += p1[N] * g[m] ; 

      backward_pass ( gr ) ;  
    }

    sum_dx = 0.0 ; 
    ip_dx = 0.0 ; 
    for ( n = 1 , gg = 0.0  ; n <= N ; n++ ) {
      gr[n] *= beta ;             /* temperature effect */
      gr[n] -= bias[n] ; /* prior effect */
      tmpd = gr[n] ; 
      gr[n] += x[n] ;               /* derivative of entropy */
      dx[n] = - tmpd -  x[n] ; /* MONITOR */
      sum_dx += dx[n] * dx[n] ;/* MONITOR */ 
      ip_dx += dx[n] * dxp[n] ;/* MONITOR */
      gg += gr[n] * gr[n] ; 
    }
    tmpd = sqrt(sum_dx) ;
    if (tmpd*sdx1 > 0) ip_dx *= 1.0/(tmpd*sdx1) ; 
    all_clear =  ( tmpd < stepmax && ip_dx >= 0.0 ) ? 1 : 0 ; 
    sdx1 = tmpd ;       /* MONITOR */
    if (  all_clear ) {
    }
    else {
      if ( ip_dx < 0.0 ) { 
	tmpd = sludge ; 
      } else {
	tmpd = stepmax/(tmpd) ; 
      }
      sdx1 *= tmpd ; 

      for ( n = 1 , gg = 0.0  ; n <= N ; n++ ) {
	dx[n] *= tmpd ; 
      }
    }
    for ( n = 1 , gg = 0.0  ; n <= N ; n++ ) {
      xp[n] = x[n] ; 
      dxp[n] = dx[n] ; /* MONITOR */
      x[n] += dx[n] ; 
    }
    
    printf("%6.3g %6.3g %6.3g : ",beta,sdx1, ip_dx ); 
    pdv( x , 21 , 34 , 63 ) ;  /* MONITOR */
    if(!(i%1))printf("\n");
    fflush(stdout);/* MONITOR */

    if ( !all_clear || i==1 ) { moratorium = START_M ; } 
    if ( moratorium > 0 ) { moratorium -- ; }
    else  {
      /* If direction reversed, decreased beta a lot */
/*      if ( ip_dx < 0.0 ) { beta *= beta_dec2 ; moratorium = START_M ; }*/
      /* If step size too small, increase beta */
      if ( ( sdx1 < stepmin ) ) /*  && ( ip_dx > ip_close ) ) */
	{ beta *= beta_inc ; }
      /* If step size too big, decrease beta */
    }
  }
  while ( beta < beta1 && i < itmax ) ; /* continue until top beta reached */
  if ( MS_DEMO ) printf ("%3d Jumps, beta %6.3g. ",i,beta ) ; 

  free_dvector ( gr , 1 , N ) ; 
  free_dvector ( dx , 1 , N ) ; /* MONITOR */
  free_dvector ( dxp , 1 , N ) ; /* MONITOR */
  free_dvector ( xp , 1 , N ) ; /* MONITOR */
}

void seq_opt (double *x , param *controlp ) /* optimization by `reestimation', sequential */
{
  double tmpd , grnn , v0 , v1 ; 
  int i = 0  , nn ; 

  v1 = objective ( x , controlp ) ; /* this sets up the qs for us too */
  v0 = v1 + 2 * ftol ; /* hack to prevent premature terminatino */

  do {
    mega_forward_pass ( ) ; 
    i ++ ; 
    for ( nn = N ; nn >= 1 ; nn-- ) {
      grnn = 0.0 ; 
      for ( m = 1 ; m <= M ; m ++ ) {
	if ( A[m][nn] ) { /* accumulate gradient  */
	  tmpd = mp1[m][nn-1] * mp1[m][nn+1] +  mp0[m][nn-1] * mp0[m][nn+1] -
	      mp0[m][nn-1] * mp1[m][nn+1] -  mp1[m][nn-1] * mp0[m][nn+1]  ;
	  grnn -= g[m] * tmpd ; 
	}
      }
      grnn *= beta ;
      grnn -= bias[nn] ; /* prior effect */
      x[nn] = - grnn ; /* evaluate the new value of q */
      x_to_q ( x[nn] , &q0[nn] , &q1[nn] , 0 ) ;
      mega_backward_step ( nn ) ;
    }
    if ( !(i%seq_period) ) {
      if ( verbose ) printf ("Seq %d %8.5g " , i , v1  ) ; 
      v0 = v1 ; 
      v1 = objective ( x , controlp ) ;
      if ( verbose ) printf ("%8.5g\n" , v1  ) ; 
      if ( v0 - v1 < -ftol  ) fprintf ( stderr , "Warning, F increased significantly %6.3g\n" , v0 - v1 ) ; 
    }
  }
  while ( v0 - v1 > ftol && i < itmax ) ;
  if ( MS_DEMO ) printf ("%3d Jumps, beta %6.3g. ",i,beta ) ; 

}

void print_state ( double tmpd , double gg , double *x ) {
/*  printf ( "%6.3g = %6.3g - %6.3g: " , tmpd , S , E ) ; */
  for ( n = 1 ; n <= N ; n++ ) printf ( "%2.0f " , 99.0 * q1[n] ) ; 
  printf ( ":%6.3g\nx: " , gg ) ;
  for ( n = 1 ; n <= N ; n++ ) printf ( "%6.3g " , x[n] ) ; 
  printf ( "\n" ) ; 
}


void printall ( FILE *fp ) 
{
  printoutcmatrix ( A , 1 , M , 1 , N ) ; 
   printf ( "\n" ) ; 
  printoutcvector ( s , 1 , N ) ; 
   printf ( "\n" ) ; 
  printoutcvector ( t , 1 , M ) ; 
  printoutcvector ( to , 1 , M ) ; 
  pdv ( g , 1 , M , 4 ) ;   printf ( "\n" ) ; 
}

/*          That's the end of the interesting code          */
/*     Here follow routines to parse the input command line */

#define DNT fprintf( fp, "\n        ")
#define NLNE  fprintf( fp, "\n")

static void print_usage ( char **argv , FILE * fp )
{
  fprintf( fp, "Usage: %s ",argv[0]);
  fprintf( fp, "datafile [optional arguments]");

  NLNE; fprintf( fp, " Arguments affecting verbosity of run-time output: <defaults>");
  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)"); 
  DNT;  fprintf( fp, "-CG                      (check gradient) "); 
  DNT;  fprintf( fp, "-MS_DEMO                 (only do demo of M-S) "); 
  DNT;  fprintf( fp, "[-e epsilon]       <%g> (step size for gradient)", epsilon); 
  NLNE; fprintf( fp, " Data creation:" ) ; 
  DNT;  fprintf( fp, "-err err           <%g> (error probability) " , err ); 
  DNT;  fprintf( fp, "-fpol fpol         <%g>  (fraction of bits high in poly)", fpol ); 
  DNT;  fprintf( fp, "-fs fs             <%g>  (fraction in shift register seed)", fs); 
  NLNE; fprintf( fp, " Inference:" ) ; 
  DNT;  fprintf( fp, "-sdx sdx           <%3g> (sd of initial state)" , sdx ); 
  DNT;  fprintf( fp, "-b betastyle beta0 beta1 (what to do with beta)" ); 
  DNT;  fprintf( fp, "       betastyle 0: const; 1: linear; 2: multiply" ); 
  DNT;  fprintf( fp, "-opt optimizer     <%d>   (0=frp,1=macopt,2=jumpopt,3=seqopt,4=betaopt)",opt); 
  DNT;  fprintf( fp, "-itmax itmax     <%d>   (num of linmins or jumps per optimization)",itmax); 
  DNT;  fprintf( fp, "-nl loops          <%d>   (number of cg runs)",NL); 
  DNT;  fprintf( fp, "-ftol ftol         <%g> ", ftol ); 
  DNT;  fprintf( fp, "-seed seed "); 
  DNT;  fprintf( fp, "-seqp seq_period   <%d>  (how many loops per check of delta F)" , seq_period ); 
  NLNE; fprintf( fp, " Files:");
  DNT;  fprintf( fp, "-o outfile               "); 
  fprintf( fp, "\n");
  return ;
}
#undef DNT
#undef NLNE

static void make_sense ( void ) 
{ /* routine to correct silly control parameters */
  if ( M < N ) fprintf ( stderr , "Warning, M < N\n" ) ;
  if ( betastyle == 1 && NL > 1 ) 	
    betaf = ( beta1 - beta ) / ((double) (NL - 1) ) ; 
  else if ( betastyle == 2 && NL > 1 )
    betaf = exp ( log ( beta1 / beta ) / ((double) (NL - 1) )) ; 
/*  if ( opt > 4 ) opt = 4 ; */
}

static int process_command ( int argc , char **argv )
{
  int p_usage = 0 ;
  int status = 0 ;
  int cs , i ;

  if ( argc < 1 )     {
    p_usage = 1 ; 
    status -- ;
  }
#define ERROR1 fprintf ( stderr , "arg to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERROR2 fprintf ( stderr , "args to `%s' missing\n" , argv[i] ) ; \
               status --
#define ERRORREG fprintf ( stderr , "regtype must be defined before `%s'\n" , argv[i] ) ; \
               status --
  for (i = 1 ; i < argc; i++)    {
    cs = 1 ;
    if ( strcmp (argv[i], "-V") == 0 )        {
      verbose = 1;
    }
    else  if ( strcmp (argv[i], "-VV") == 0 )        {
      verbose = 2;
    }
    else  if ( strcmp (argv[i], "-CG") == 0 )        {
      CG = 1;
    }
    else  if ( strcmp (argv[i], "-NL") == 0 )        {
      NL = 0; 
    }
    else  if ( strcmp (argv[i], "-MS_DEMO") == 0 )        {
      MS_DEMO = 1; 
    }
    else  if ( strcmp (argv[i], "-nl") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &NL); 
    }
    else  if ( strcmp (argv[i], "-opt") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &opt); 
    }
    else  if ( strcmp (argv[i], "-itmax") == 0 ) {
      if ( i + 1 == argc ) { ERROR1; }
      else cs *= sscanf(argv[++i], "%d", &itmax ) ; 
    }
    else   if ( strcmp (argv[i], "-seed") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%ld", &seed); 
    }  
    else  if ( strcmp (argv[i], "-seqp") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &seq_period); 
    }
/*    else   if ( strcmp (argv[i], "-m") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else {
	cs *= sscanf(argv[++i], "%d", &M); 
      }
    } */ /* NB M is defined by M0 */
    else   if ( strcmp (argv[i], "-n0") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else {
	cs *= sscanf(argv[++i], "%d", &N0); 
      }
    }
    else   if ( strcmp (argv[i], "-m0") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else {
	cs *= sscanf(argv[++i], "%d", &M0); 
      }
    }
/*    else   if ( strcmp (argv[i], "-n") == 0 )        {
      if ( i + 1 == argc )             { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &N); 
    }*/ /* NB N is defined to be M0 */
    else if ( strcmp (argv[i], "-e") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &epsilon); 
    }
    else if ( strcmp (argv[i], "-err") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &err ); 
    }
    else if ( strcmp (argv[i], "-sdx") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &sdx); 
    }
    else if ( strcmp (argv[i], "-b") == 0 ) {
      if ( i + 3 >= argc ) { ERROR2;      }
      else {
	cs *= sscanf(argv[++i], "%d", &betastyle); 
	cs *= sscanf(argv[++i], "%lf", &beta); 
	cs *= sscanf(argv[++i], "%lf", &beta1); 
      }
    }
    else if ( strcmp (argv[i], "-ftol") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &ftol); 
    }
    else if ( strcmp (argv[i], "-fs") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &fs); 
    }
    else if ( strcmp (argv[i], "-fpol") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &fpol); 
    }
    else if ( strcmp (argv[i], "-pheading") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &pheading); 
    }
    else  if ( strcmp (argv[i], "-o") == 0 )	{
      if ( i + 1 == argc ) 	    { ERROR1;      }
      else {
	strcpy(outfile, argv[++i]);
	if ( !printout ) printout = 1 ;
      }
    }
    else {
      fprintf ( stderr , "arg `%s' not recognised\n" , argv[i] ) ; 
      p_usage = 1 ;
      status -- ;
    }
    if ( cs == 0 ) {
      fprintf ( stderr , "arg at or before `%s' has incorrect format\n" , 
	       argv[i] ) ;
      p_usage = 1 ;
      status -- ;
    }
  }
  if ( p_usage ) print_usage ( argv , stderr ) ;
  return ( status ) ;
}
#undef ERROR1
#undef ERROR2
#undef ERRORREG