cg.c

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

    tmpd =  ( *qq0 * l0 + *qq1 * l1 ) ; 
    return tmpd ;
  }
  else return 0.0 ; 
}

void find_q_S ( double *x , int sflag )
{
  double tmpd = 0.0 ; 

  for ( n = 1 ; n <= N ; n++ ) {
    tmpd += x_to_q ( x[n] , &q0[n] , &q1[n] , sflag ) ;
  }
  if ( sflag ) S = tmpd ; 
}

double objective ( double *x , void *f_args ) 
{
  double tmpd = 0.0 ; 

  find_q_S ( x , 1 ) ; /* puts S in S */
  E = 0.0 ; 

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

  tmpd = S - beta * E ; 
  if ( verbose > 0 ) 
    printf ( "%6.4g " , tmpd ) ; 
/*    print_state ( tmpd , 0.0 , x ) ; */
  return ( tmpd ) ;
}

void forward_pass ( void ) 
{
  for ( n = 1 ; n <= N ; n++ ) {
    if ( A[m][n] == 0  ) {
      p0[n] = p0[n-1] ;
      p1[n] = p1[n-1] ;
    }
    else {
      p0[n] = q0[n] * p0[n-1] + q1[n] * p1[n-1] ;
      p1[n] = q0[n] * p1[n-1] + q1[n] * p0[n-1] ;
    }
  }
}

void mega_forward_pass ( void ) 
{
  for ( m = 1 ; m <= M ; m++ ) {
    for ( n = 1 ; n <= N ; n++ ) {
      if ( A[m][n] == 0  ) {
	mp0[m][n] = mp0[m][n-1] ;
	mp1[m][n] = mp1[m][n-1] ;
      }
      else {
	mp0[m][n] = q0[n] * mp0[m][n-1] + q1[n] * mp1[m][n-1] ;
	mp1[m][n] = q0[n] * mp1[m][n-1] + q1[n] * mp0[m][n-1] ;
      }
    }
  }
}

void backward_pass ( double *gr ) 
{
  double tmpd ; 
  for ( n = N ; n >= 1 ; n-- ) {
    if ( A[m][n] == 0 ) {
      p0[n] = p0[n+1] ;
      p1[n] = p1[n+1] ;
    }
    else {
      p0[n] = q0[n] * p0[n+1] + q1[n] * p1[n+1] ;
      p1[n] = q0[n] * p1[n+1] + q1[n] * p0[n+1] ;
      /* gradient accumulation only if A_mn = 1 */
      tmpd = p1[n-1] * p1[n+1] +  p0[n-1] * p0[n+1] -
	     p0[n-1] * p1[n+1] -  p1[n-1] * p0[n+1]  ;
/*      if ( tmpd != 0.0 )   will I regret this time saver? */
      gr[n] -= g[m] * tmpd ; 
    }
  }
}

void mega_backward_step ( int nn )  /* nn should start at N and go backwards */
{
  for ( m = 1 ; m <= M ; m++ ) {
      if ( A[m][nn] == 0 ) {
	mp0[m][nn] = mp0[m][nn+1] ;
	mp1[m][nn] = mp1[m][nn+1] ;
      }
      else {
	mp0[m][nn] = q0[nn] * mp0[m][nn+1] + q1[nn] * mp1[m][nn+1] ;
	mp1[m][nn] = q0[nn] * mp1[m][nn+1] + q1[nn] * mp0[m][nn+1] ;
    }
  }
}

void d_objective( double *x, double *gr, void *f_args)
{
  double tmpd , gg ; 
  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 ; 

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

    backward_pass ( gr ) ; 
  }

  for ( n = 1 , gg = 0.0  ; n <= N ; n++ ) {
    gr[n] *= beta ;             /* temperature effect */
    gr[n] += x[n] ;               /* derivative of entropy */
    gr[n] *= q0[n] * q1[n] ;      /* factor omitted from all terms */
    gg += gr[n] * gr[n] ; 
  }

  tmpd = S - E ; 
  if ( verbose ) {
    print_state ( tmpd , gg , x ) ; 
  }
}

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

  gr = dvector ( 1 , 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 ; 

    if ( i > 1 && omit_diagonal )  /* omit diagonal elements - Note this 
				      is not correct in the annealing case, 
				      where on subsequent loops omission 
				      should also occur. */
      m = diagno + 1 ;
    else 
      m = 1 ;
    for ( ; m <= M ; m ++ ) { 
      forward_pass ( ) ; 
      E += p1[N] * g[m] ; 

      backward_pass ( gr ) ;  
    }

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

    tmpd = S - beta * E ; 
    if ( verbose ) {
      print_state ( tmpd , gg , x ) ; 
    }
  }
  while ( gg > ftol && i < itmax ) ;

  free_dvector ( gr , 1 , N ) ; 
}

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

}

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 ) 
{
  printoutimatrix ( A , 1 , M , 1 , N ) ; 
   printf ( "\n" ) ; 
  printoutivector ( s , 1 , N ) ; 
   printf ( "\n" ) ; 
  printoutivector ( t , 1 , M ) ; 
  printoutivector ( 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, "[-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, "-fA fA             <%g>  (fraction of bits high in A)", fA ); 
  DNT;  fprintf( fp, "-d diagno          <%d>   (number of rows of A that are singlet)", diagno ); 
  DNT;  fprintf( fp, "-OD                     (whether to omit diagonal entries after first jumpopt)" ) ; 
  DNT;  fprintf( fp, "-fs fs             <%g> ", fs); 
  DNT;  fprintf( fp, "-SO                <%d> (whether to order s thus: 11110000)", sordered ); 
  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)",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               "); 
  NLNE; fprintf( fp, " Monte Carlo options (after optimization):");
  DNT;  fprintf( fp, "-mc mc_loops               "); 
  DNT;  fprintf( fp, "-mco mcfile               "); 
  DNT;  fprintf( fp, "-mcv verbosemc               "); 
  DNT;  fprintf( fp, "-mcm metropolis(1)or_gibbs(0)  "); 
  DNT;  fprintf( fp, "-mcr random_selection(1/0) "); 
  fprintf( fp, "\n");
  return ;
}
#undef DNT
#undef NLNE

static void make_sense ( void ) 
{ /* routine to correct silly control parameters */
  if ( ( diagno < 0 ) || ( diagno > N )  ) diagno = N ; 
  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 ( MC ) mcboredom *= N ; 
  if ( opt > 3 ) opt = 3 ;
}

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], "-SO") == 0 )        {
      sordered = 1;
    }
    else  if ( strcmp (argv[i], "-OD") == 0 )        {
      omit_diagonal = 1;
    }
    else  if ( strcmp (argv[i], "-NL") == 0 )        {
      NL = 0; 
    }
    else  if ( strcmp (argv[i], "-nl") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &NL); 
    }
    else  if ( strcmp (argv[i], "-mc") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &MC); 
    }
    else  if ( strcmp (argv[i], "-d") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &diagno); 
    }
    else  if ( strcmp (argv[i], "-mcv") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &verbosemc); 
    }
    else  if ( strcmp (argv[i], "-mcm") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &mcm ); 
    }
    else  if ( strcmp (argv[i], "-mcstop") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &mcstop ); 
    }
    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); 
      }
    }
    else   if ( strcmp (argv[i], "-n") == 0 )        {
      if ( i + 1 == argc )             { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &N); 
    }
    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], "-fA") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &fA); 
    }
    else  if ( strcmp (argv[i], "-o") == 0 )	{
      if ( i + 1 == argc ) 	    { ERROR1;      }
      else {
	strcpy(outfile, argv[++i]);
	if ( !printout ) printout = 1 ;
      }
    }
    else  if ( strcmp (argv[i], "-mco") == 0 )	{
      if ( i + 1 == argc ) 	    { ERROR1;      }
      else {
	strcpy(mcfile, argv[++i]);
	if ( !printmc ) printmc = 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