ms.c

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

  p1 = dvector ( 0 , N + 1 ) ; 
  p0[0] = 1.0 ; p1[0] = 0.0 ; 
  p0[N+1] = 1.0 ; p1[N+1] = 0.0 ; 

/* data already made */

  for ( m = 1 ; m <= M ; m++ ) {
    g[m] = - 1.0 ; 
    if ( to[m] ) g[m] = - g[m] ; /* log odds contributed by `observation' */
  }
  true_err = (double) count_err / (double) M0 ; 

  bias = dvector ( 1 , N ) ; 
  for ( n = 1 ; n <= N ; n ++ ) 
    bias[n] = log ( err / ( 1 - err ) ) ; 
/* set up model */ 

  for  ( n = 1 ; n <= N ; n++ ) { /* initial condition */
    x[n] = bias[n] + sdx * rann() ; 
  }
  if ( verbose > 1 ) printall ( stdout ) ; 

  if( CG )
    checkgrad ( x , N , epsilon , objective , &control , 
	       d_objective , &control);

/*   Main loop        loop NL times    */

  for ( nl = 1 ; nl <= NL ; nl ++ ) {
    if( NL > 1 && verbose > 0 ) 
      printf ( "Loop %d of %d ; tol %5g ; beta %5g\n" , 
	      nl , NL , ftol , beta ) ;
    if ( opt == 0 ) 
      frprmn ( x , N , ftol , flintol , &iter , itmax , &v ,
	     objective , &control , d_objective , &control ) ;
    else if ( opt == 1 ) 
      macopt ( x , N , rich , ftol , &iter , itmax , 
	     d_objective , &control ) ;
    else if ( opt == 2 ) 
      jump_opt ( x ) ;
    else if ( opt == 4 ) 
      beta_opt ( x ) ;
    else 
      seq_opt ( x , &control ) ;

    if( CG )
      checkgrad ( x , N , epsilon , objective , &control , 
		 d_objective , &control);

/*   } */  /* report every loop */

/* End of main loop */

/* produce answer */

    find_q_S ( x , 0 ) ; 
    if ( verbose >= 1 ) print_state ( 0.0 , 0.0, x ) ; 

    count = 0 ; 
    count_high = 0 ; 
    for ( n = 1 ; n <= N ; n++ )  {
      so[n] = ( x[n] > 0 ) ? 1 : 0 ;
      if ( so[n] != s[n] ) {
	count ++ ;
      }
      count_high += so[n] ;
    }
    if ( verbose > 0 ) {
      printf ( "\n          " ) ;
      printoutcvector ( s , 1 , N ) ; 
    }
    count_viol = 0 ; 
    for ( m = 1 ; m <= M ; m++ ) {
      /* because this is in fact noise free, the observed to is the same as t, 
       so I now use t to check out the answer's t */
      t[m] = cdotprod_mod2 ( A[m] , so , 1 , N ) ; 
      if ( to[m] != t[m] ) count_viol ++ ; 
    }
    if ( count == 0 )  nl = NL ; /* ensure stop now */
    if (( verbose >= 1 ) || MS_DEMO ||  ( nl == NL )) {
      if ( count_viol > 0 ) printf("INVALID SOLN - %d violations.  ",count_viol);
      else if ( count > 0 ) printf("VALID SOLN.   ");
      else { 
	printf("CORRECT!    ");
      }
      printf( " diffs %d   true errors %d   reconstructed errors %d\n", count , count_err , count_high ) ; 
    }

/* this belongs at end of main loop */
    if ( nl < NL ) {
      if ( betastyle == 1 ) 	
	beta += betaf ;
      else if ( betastyle == 2 )
	beta *= betaf ; 
    }
  }
  v = (double) ( count_high - count_err +
		( ( count_viol > 0 ) ? 1000 : 0 ) ) ; /* this marks out
							 the total failures */
  /* if v is negative then a better solution "*" found */
  if ( verbose > 0 ) finalline ( count , v , stdout ) ;
  if ( printout ) {
    fp = fopen ( outfile , "a" ) ;
    if( !fp )   fprintf( stderr, "No such file: %s\n", outfile ) ;
    else {
      finalline ( count , v , fp ) ;
      fclose ( fp ) ;
    }
  }
}

void finalline ( int count , double v , FILE *fp ) 
{
  if ( ranu() < pheading ) fprintf ( fp , 
                 "#es rl.scre N0   M0     M    fpol     err bstl seed taps tr_err    beta\n") ; 

  fprintf ( fp , "%3d %7.4g %3d %4d %5d %7.4g %7.4g %d %8ld %3d %7.4g %7.4g\n" ,
                  count , v , N0 , M0 , M , fpol , err , betastyle , seed , taps , true_err , beta ) ;
}

double x_to_q ( double x , double *qq0 , double *qq1 , 
	       int sflag ) /* if sflag = 1, returns the binary entropy */
{
  double tmpd , l0 , l1 ;

/*  x *= beta ;  temperature was at first implemented here, but 
                 for the gradients it is better to put it elsewhere */

  if ( x > 0.0 ) {
    tmpd = 1.0 + exp ( - x ) ;
    if ( sflag ) {
      l1 = - log ( tmpd ) ; 
      l0 = -x + l1 ; 
    }
    *qq1 = 1.0/tmpd ; *qq0 = 1 - *qq1 ; 
  } else {
    tmpd = 1.0 + exp ( x ) ;
    if ( sflag ) {
      l0 = - log ( tmpd ) ; 
      l1 = x + l0 ; 
    }
    *qq0 = 1.0/tmpd ; *qq1 = 1 - *qq0 ; 
  }
  if ( sflag ) {
    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 ; 
  P = 0.0 ; 

  for ( n = 1 ; n <= N ; n++ ) {
    P += bias[n] * q1[n] ; /* THIS IS RIGHT ? */
  }

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

  tmpd = S - beta * E - P ;  /* THIS IS RIGHT ? */
  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] -= bias[n] ; /* prior 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 - beta * E ;  /* THIS IS RIGHT ? (no) */
  if ( verbose ) {
    print_state ( tmpd , gg , x ) ; 
  }
}

void jump_opt (double *x ) /* optimization by `reestimation' */
{
  double tmpd , gg ; 
/*, sum_dx ; , ip_dx , sdx1=0.0 , sdxp ; */
  double *gr  ;
/* , *dx , *dxp ; */
  int i = 0  ; 

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

  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 */
/*      dxp[n] = dx[n] ;  MONITOR */
      x[n] = - tmpd ; 
      gg += gr[n] * gr[n] ; 
    }
/*    sdxp = sdx1 ;                 MONITOR */
/*    sdx1 = sqrt(sum_dx);          MONITOR */
/*    printf("%6.3g %6.3g : ",sdx1, (sdx1*sdxp > 0) ? ip_dx/(sdx1*sdxp) : 0.0 ); */
/*    pdv( x , 1 , 10 , 63 ) ;     MONITOR */
/*     if(!(i%1))printf("\n"); */
/*    fflush(stdout);   MONITOR */

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