mnc.c

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

  /* set up the norder */
  alist.norder = ivector ( 1, N ) ; 
  for ( n = 1 ; n <= N ; n++ ) {
    alist.norder[n] = n ; 
  }

   if ( DEMO ) {
     count_err = 0 ; 
     for ( n = 1 ; n <= N ; n ++ ) 
       if ( alist.num_nlist[n] == 1 ) count_err ++ ; 
     printf ("unchecked bits n : %d ", count_err ) ;
     count_err = 0 ; 
     for ( m = 1 ; m <= M ; m ++ ) 
       if ( alist.num_mlist[m] == 1 ) count_err ++ ; 
     printf (" m : %d\n", count_err ) ;
     
     exit(0) ;  
   }

/* old program from here */

/*  er = dvector ( 1 , M ) ;  */
  g  = dvector ( 1 , M ) ; 
  x  = dvector ( 1 , N ) ; 
  q0 = dvector ( 1 , N ) ; 
  q1 = dvector ( 1 , N ) ; 
/* opt = 3 is the only choice */
  if ( opt == 3 ) { /* mega state vectors needed */
    mp0 = dmatrix ( 1 , M , 0 , alist.biggest_num_m + 1 ) ;
    mp1 = dmatrix ( 1 , M , 0 , alist.biggest_num_m + 1 ) ;
    for ( m = 1 ; m <= M ; m ++ ) {
      mp0[m][0] = 1.0 ; mp1[m][0] = 0.0 ; 
      mp0[m][alist.num_mlist[m]+1 ] = 1.0 ; mp1[m][alist.num_mlist[m]+1 ] = 0.0; 
    }
  }
/* NB can't use p0 and p1 with fixed length for backward pass */
  p0 = dvector ( 0 , alist.biggest_num_m + 1 ) ; 
  p1 = dvector ( 0 , alist.biggest_num_m + 1 ) ; 

  p0[0] = 1.0 ; p1[0] = 0.0 ; 
/*
  p0[alist.biggest_num_m + 1] = 1.0 ; p1[alist.biggest_num_m+1] = 0.0 ; 
  This is handled at the beginning of backpass.
*/

/* data already made */

  for ( m = 1 ; m <= M ; m++ ) {
    g[m] = log ( err / ( 1 - err ) ) ; 
    if ( z[m] ) g[m] = - g[m] ; /* log odds contributed by observation */
  }

  bias = dvector ( 1 , N ) ; 
  for ( n = 1 ; n <= N ; n ++ ) 
    bias[n] = log ( fs / ( 1 - fs ) ) ; 
/* 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 
      seq_opt ( x , &control ) ;

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

/* 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 ) ; 
      printf ( "\n          " ) ;
      printoutcvector ( so , 1 , N ) ; 
    }
    if ( count == 0 )  nl = NL ; /* ensure stop when 
					     any valid solution is found */
    if (( verbose >= 1 ) || DEMO ||  ( nl == NL )) {
      printf( " diffs %d   true highs %d   reconstructed high bits %d\n", count , count_s , count_high ) ; 
    }

/* this belongs at end of main loop */
    if ( nl < NL ) {
      if ( betastyle == 1 ) 	
	beta += betaf ;
      else if ( betastyle%10 == 2 )
	beta *= betaf ; 
    }
  }
  v = (double) ( count_high - count_err ) ; 
  /* 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 ) 
{
  ran_seed( seed ) ; 
  if ( ranu() < pheading ) fprintf ( fp , 
                 "#es rl.scre  N  t   fs     err bstl seed taps  tr_err    beta\n") ; 

  fprintf ( fp , "%4d %6.4g %3d %4d %7.4g %7.4g %d %8ld %7.4g %7.4g %7.4g\n" ,
                  count , v , N , t , fs , err , betastyle , seed  , true_fs, 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 ) /* this could be turned into param *controlp */ 
{
  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 ( ) ; returns value of ` p1[N] ' */
    E += forward_pass ( ) * g[m] ; 
  }

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

double forward_pass ( void ) /* m is passed globally */
{
  double tmpd ;
  int l ; 

 /*  for ( n = 1 ; n <= N ; n++ ) { */
  for ( l = 1 ; l <= alist.num_mlist[m] ; l ++ ) {
      n = alist.mlist[m][l] ; 

      p0[l]        = q0[n] * p0[l-1] + q1[n] * p1[l-1] ;
      tmpd = p1[l] = q0[n] * p1[l-1] + q1[n] * p0[l-1] ;
  }
  return tmpd ; 
}

void mega_forward_pass ( void ) /* this routine allows the n order to be
				   jumbled up */
{
  int l , ll , u ; 

  for ( m = 1 ; m <= M ; m++ ) 
    alist.l_up_to[m] = 0 ; 

  for ( u = 1 ; u <= N ; u++ ) {
    n = alist.norder[u] ; 
    for ( l = 1 ; l <= alist.num_nlist[n] ; l ++ ) {
      m = alist.nlist[n][l] ; 
      alist.l_up_to[m] ++ ;
      ll = alist.l_up_to[m] ; 

      mp0[m][ll] = q0[n] * mp0[m][ll-1] + q1[n] * mp1[m][ll-1] ;
      mp1[m][ll] = q0[n] * mp1[m][ll-1] + q1[n] * mp0[m][ll-1] ;
    }
  }
}

/* this one runs over m first, then along l. But it'll be nice to jumble 
   up the n order without having to rewrite the alist.mlist 
void mega_forward_pass ( void ) 
{
  for ( m = 1 ; m <= M ; m++ ) {
    for ( l = 1 ; l <= alist.num_mlist[m] ; l ++ ) {
      n = alist.mlist[m][l] ; 
      mp0[m][l] = q0[n] * mp0[m][l-1] + q1[n] * mp1[m][l-1] ;
      mp1[m][l] = q0[n] * mp1[m][l-1] + q1[n] * mp0[m][l-1] ;
    }
  }
}
*/

void backward_pass ( double *gr ) 
{
  double tmpd ; 
  int l ; 

  l = alist.num_mlist[m] ;
  p0[ l+1 ] = 1.0 ; p1[ l+1 ] = 0.0; 
  for ( ; l >= 1 ; l -- ) {
      n = alist.mlist[m][l] ; 
    
/*  for ( n = N ; n >= 1 ; n-- ) { */

      p0[l] = q0[n] * p0[l+1] + q1[n] * p1[l+1] ;
      p1[l] = q0[n] * p1[l+1] + q1[n] * p0[l+1] ;

      tmpd = p1[l-1] * p1[l+1] +  p0[l-1] * p0[l+1] -
	     p0[l-1] * p1[l+1] -  p1[l-1] * p0[l+1]  ;
      gr[n] -= g[m] * tmpd ; 
  }
}

void mega_backward_step ( int nn )  /* nn should go backwards 
				     through the sequence used in megaforward */
{
  int l , ll ; 
  for ( l = 1 ;  l <= alist.num_nlist[nn] ; l ++ ) {
      m = alist.nlist[nn][l] ; 
/*  for ( m = 1 ; m <= M ; m++ ) { */
      ll = alist.l_up_to[m] ;
      alist.l_up_to[m] -- ;

	mp0[m][ll] = q0[nn] * mp0[m][ll+1] + q1[nn] * mp1[m][ll+1] ;
	mp1[m][ll] = q0[nn] * mp1[m][ll+1] + q1[nn] * mp0[m][ll+1] ;
  }
}

void d_objective( double *x, double *gr, void *f_args)
{
  double tmpd , gg ; 
  fprintf ( stderr , 
	   " gradient routine called but carefully tested yet \n") ;
  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 ++ ) {
    E += forward_pass ( ) * g[m] ; 

    backward_pass ( gr ) ; 
  }

  for ( n = 1 , gg = 0.0  ; n <= N ; n++ ) {
    gr[n] *= beta ;             /* temperature effect */
    gr[n] -= betap * 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 ++ ) { 
      E +=     forward_pass ( ) * 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] -= betap * 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 ( 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 */
}


void seq_opt (double *x , param *controlp )
/* optimization by `reestimation', sequential */
/* an advantage of this approach is that there is a free energy
   improvement guarantee. Also we can 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 tmpd , grnn , v0 , v1 ; 
  int i = 0  , nn , l , ll ; 
  int u ; /* runs over the norder */
/*  alist_matrix *a = controlp->a ;  */

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

  do {
    mega_forward_pass ( ) ; 
    i ++ ;