msb.c

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

      if ( (NORDER==2) ) 
	my[n].val = suspicion[n+1] ;
      else if ( NORDER==1 )
	my[n].val =   rani(10*N) ; 
      else /* NORDER==3 */
	my[n].val = ( suspicion[n+1] * 1000000) / num_nlist[n+1] ; 
      /* 	if ==1,  replace suspicions by random numbers   */
    }
    /* sort norder by suspicion; biggest suspicion should go to norder[N], 
     so that it is the first to be touched in the backward pass */
    qsort ( my , (unsigned int) (N) , sizeof(thing) , cmp_thing ) ; 
    if ( MS_DEMO || verbose >= 2  ) printf ( " Order for bits 1..N\n" ) ; 
    for ( n = 0 ; n <= N-1 ; n ++ ) {
      norder[n+1] = my[n].index ; 
      if ( MS_DEMO || verbose >= 2  ) printf ( "%2d " , norder[n+1] ) ; 
    }
    free((char * ) (my ) ) ; 
    if ( MS_DEMO || verbose >= 2  ) {
      printf ( "\n" ) ; 
      for ( n = 1 ; n <= N ; n ++ ) {
	printf ( "%2d " , suspicion[ norder[n]] ) ; 
      }
      printf ( "\n" ) ; 
    }
  }
  else  /* ( NORDER == 0 ) */ { 
    for ( n = 1 ; n <= N ; n++ ) {
      norder[n] = n ; 
    }
  }


/* old program from here */

/*  er = dvector ( 1 , M ) ;  */
  g  = dvector ( 1 , M ) ; 
  x  = dvector ( 1 , N ) ; 
  q0 = dvector ( 1 , N ) ; 
  q1 = dvector ( 1 , N ) ; 
  if ( opt == 3 ) { /* mega state vectors needed */
    mp0 = dmatrix ( 1 , M , 0 , taps + 1 + 1 ) ;
    mp1 = dmatrix ( 1 , M , 0 , taps + 1 + 1 ) ;
    for ( m = 1 ; m <= M ; m ++ ) {
      mp0[m][0] = 1.0 ; mp1[m][0] = 0.0 ; 
      mp0[m][taps+1+1] = 1.0 ; mp1[m][taps+1+1] = 0.0 ; 
/* a more general program would here use num_mlist[m] instead */
    }
  }
  p0 = dvector ( 0 , taps+1 + 1 ) ; 
  p1 = dvector ( 0 , taps+1 + 1 ) ; 
  p0[0] = 1.0 ; p1[0] = 0.0 ; 
  p0[taps+1+1] = 1.0 ; p1[taps+1+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 
      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_index ( mlist[m] , so , 1 , num_mlist[m] ) ; 
    /* cdotprod_mod2 ( A[m] , so , 1 , N ) ;  */ 
      if ( to[m] != t[m] ) count_viol ++ ; 
    }
    if ( count == 0 || count_viol == 0 )  nl = NL ; /* ensure stop when 
					     any valid solution is found */
    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 || betastyle == 22 )
	beta *= betaf ; 
    }
    if ( betastyle == 22 && nl == NL - 1 ) {
      beta = betahuge ; 
      betap = betap_huge ; 
      for  ( n = 1 ; n <= N ; n++ ) { /* back to initial condition */
	x[n] = bias[n] + sdx * rann() ; 
	norder[n] = n ; /* return to sensible order */
      }
    }
  }
  v = (double) ( count_high - count_err +
		( ( count_viol > 0 ) ? 10000 : 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 ) 
{
  ran_seed( seed ) ; 
  if ( ranu() < pheading ) fprintf ( fp , 
                 "#es rl.scre  N0   M0     M    fpol     err bstl seed taps  tr_err    beta\n") ; 

  fprintf ( fp , "%4d %6.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 ( ) ; 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 <= num_mlist[m] ; l ++ ) {
      n = 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++ ) 
    l_up_to[m] = 0 ; 

  for ( u = 1 ; u <= N ; u++ ) {
    n = norder[u] ; 
    for ( l = 1 ; l <= num_nlist[n] ; l ++ ) {
      m = nlist[n][l] ; 
      l_up_to[m] ++ ;
      ll = l_up_to[m] ; /* final value should be taps+1 */

      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 mlist 
void mega_forward_pass ( void ) 
{
  for ( m = 1 ; m <= M ; m++ ) {
    for ( l = 1 ; l <= num_mlist[m] ; l ++ ) {
      n = 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 ; 

  for ( l = num_mlist[m] ; l >= 1 ; l -- ) {
      n = 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 <= num_nlist[nn] ; l ++ ) {
      m = nlist[nn][l] ; 
/*  for ( m = 1 ; m <= M ; m++ ) { */
      ll = l_up_to[m] ;
      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 ; 
  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 {