mnc.c

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

    for ( u = N ; u >= 1 ; u-- ) { /* reverse sequence used in mega_for */
      nn = alist.norder[u] ; 
      grnn = 0.0 ; 
      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] ;
	tmpd = mp1[m][ll-1] * mp1[m][ll+1] +  mp0[m][ll-1] * mp0[m][ll+1] -
	       mp0[m][ll-1] * mp1[m][ll+1] -  mp1[m][ll-1] * mp0[m][ll+1]  ;
	grnn -= g[m] * tmpd ; 
      }
      grnn *= beta ;
      grnn -= betap * bias[nn] ; /* prior effect */
      x[nn] = - grnn ; /* evaluate the new value of q */
      x_to_q ( x[nn] , &q0[nn] , &q1[nn] , 0 ) ;
      if ( seqverbose >= 2 ) printf("%4d %8.3g\n",nn,x[nn] ) ; 
/* print_state ( 0.0 , 0.0 , x ) ;  */
      mega_backward_step ( nn ) ;
    }
    if ( seqverbose >= 2 ) printf("\n");
    if ( !(i%seq_period) ) { /* default period = 10 */
/*      if ( seqverbose ) printf ("Seq %d %8.5g " , i , v1  ) ; */
      v0 = v1 ; 
      v1 = objective ( x , controlp ) ;
      if ( seqverbose ) printf ("%8.5g %8.5g\n" , v1 , (v0-v1) ) ; 
      if ( v0 - v1 < -ftol  ) fprintf ( stderr , "Warning, F increased significantly %6.3g\n" , v0 - v1 ) ; 
    }
  }
  while ( v0 - v1 > ftol && i < itmax ) ;
  if (seqverbose) 
    printf ("seq_opt done: %8.3g %3d Js, beta %6.3g\n",(v0-v1),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 ) 
{
  printoutcvector ( s , 1 , N ) ; 
   printf ( "\n" ) ; 
  printoutcvector ( t , 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, " Verbosity of run-time output: <defaults>");
  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)"); 
  DNT;  fprintf( fp, "-SV | -SVV               (verbose in seq_opt)"); 
  DNT;  fprintf( fp, "-CG                      (check gradient) "); 
  DNT;  fprintf( fp, "-DEMO                 (only do demo of MNC) "); 
  DNT;  fprintf( fp, "[-e epsilon]       <%g> (step size for gradient)", epsilon); 
  NLNE; fprintf( fp, " Data creation:" ) ; 
  DNT;  fprintf( fp, "-n N             <%d> (number of bits in string) " , N ); 
  DNT;  fprintf( fp, "-a astyle        <%d> (style of A creation) " , how_to_gen_A ); 
  DNT;  fprintf( fp, "-u ustyle        <%d> (if astyle=1 : style of bit-add matrix creation) " , ustyle ); 
  DNT;  fprintf( fp, "-k K             <%d> (ulist style: number of bit-adds in code [even]) " , K ); 
  DNT;  fprintf( fp, "-t t             <%d> (if astyle=3: number of 1s per row) " , num_per_row ); 
  DNT;  fprintf( fp, "-dens dens         <%g> (LU style: density of LU) " , dens ); 
  DNT;  fprintf( fp, "-err err           <%g> (error probability) " , err ); 
  DNT;  fprintf( fp, "-fs fs             <%g>  (fraction in string)", fs); 
  DNT;  fprintf( fp, "-pbm pbm_ofile     <%s>  (write A.pbm)", pbm_ofile); 
  DNT;  fprintf( fp, "-xv 0/1            <%d>  (xv it too)", pbm_xv ); 

  NLNE; fprintf( fp, " Inference:" ) ; 
  DNT;  fprintf( fp, "-RR                   (use Radford Representation) "); 
  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; 22: multiply and go wild on last loop" ); 
  DNT;  fprintf( fp, "-betap betap       <%3g> (prior's beta)" , betap ); 
  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, "-norder n        <%d>   (for seqopt - 0/1/2 = 1..N / random / suspicion / suspicion/num )",NORDER); 
  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 */
  M = 2 * N ; 
  if ( betastyle == 1 && NL > 1 ) 	
    betaf = ( beta1 - beta ) / ((double) (NL - 1) ) ; 
  else if ( ( betastyle == 2 || betastyle == 22 ) && NL > 1 )
    betaf = exp ( log ( beta1 / beta ) / ((double) (NL - 1) )) ; 
/*  if ( opt > 4 ) opt = 4 ; */
  if ( opt != 3 ) {
    fprintf ( stderr , "opt = 3 override \n" ) ; 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], "-SV") == 0 )        {
      seqverbose = 1;
    }
    else  if ( strcmp (argv[i], "-SVV") == 0 )        {
      seqverbose = 2;
    }
    else  if ( strcmp (argv[i], "-CG") == 0 )        {
      CG = 1;
    }
    else  if ( strcmp (argv[i], "-NL") == 0 )        {
      NL = 0; 
    }
    else  if ( strcmp (argv[i], "-DEMO") == 0 )        {
      DEMO = 1; 
    }
    else  if ( strcmp (argv[i], "-RR") == 0 )        {
      Radford_rep = 1; 
    }
    else  if ( strcmp (argv[i], "-LUDEMO") == 0 )        {
      LUDEMO = 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], "-n") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &N); 
    }
    else  if ( strcmp (argv[i], "-u") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &ustyle); 
    }
    else  if ( strcmp (argv[i], "-a") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &how_to_gen_A ); 
    }
    else  if ( strcmp (argv[i], "-opt") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &opt); 
    }
    else  if ( strcmp (argv[i], "-norder") == 0 ) {
      if ( i + 1 == argc ) { ERROR1; }
      else cs *= sscanf(argv[++i], "%d", &NORDER ) ; 
    }
    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], "-k") == 0 ) {
      if ( i + 1 == argc ) { ERROR1;      }
      else {
	cs *= sscanf(argv[++i], "%d", &K); 
      }
    }
    else if ( strcmp (argv[i], "-e") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &epsilon); 
    }
    else if ( strcmp (argv[i], "-dens") == 0 )      {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &dens ); 
    }
    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], "-betap") == 0 )        {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%lf", &betap); 
    }
    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], "-t") == 0 ) {
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &num_per_row ); 
    }
    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  if ( strcmp (argv[i], "-xv") == 0 )	{
      if ( i + 1 == argc ) { ERROR1;      }
      else cs *= sscanf(argv[++i], "%d", &pbm_xv ); 
    }
    else  if ( strcmp (argv[i], "-pbm") == 0 )	{
      if ( i + 1 == argc ) 	    { ERROR1;      }
      else {
	strcpy(pbm_ofile, argv[++i]);
	if ( !pbm_o ) pbm_o = 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

unsigned char cdotprod_mod2_index 
( 
 int *a , unsigned char *b , int lo , int hi ) {
  unsigned char ans = 0 ; 
  int i ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans ^= b[a[i]] ;
  }
  return ans ;
}


int cdotprod_index 
( 
 int *a , unsigned char *b , int lo , int hi ) {
  int ans = 0 ; 
  int i ; 
  for ( i = lo ; i <= hi ; i ++ ) {
    ans += b[a[i]] ;
  }
  return ans ;
}

static void initialize_alist ( alist_matrix *a , int NN , int MM , 
			      int biggestn , int biggestm ) {
  int mm , nn ; 

  a->N = NN ; a->M = MM ; 
  a->biggest_num_m =    biggestm ;
  a->biggest_num_n = biggestn ;
    
  a->num_mlist = ivector ( 1 , NN ) ; 
  a->mlist = imatrix ( 1 , NN , 1 , a->biggest_num_m ) ; 
  a->num_nlist = ivector ( 1 , NN ) ; 
  a->nlist = imatrix ( 1 , NN , 1 , a->biggest_num_n ) ; 
  a->l_up_to = ivector ( 1 , NN ) ; 
  
  for ( mm = 1 ; mm <= NN ; mm++ ) {
    a->num_mlist[mm] = 0 ; 
  }
  for ( nn = 1 ; nn <= NN ; nn++ ) {
    a->num_nlist[nn] = 0 ; 
  }
}

static void cmatrix_2_alist ( unsigned char **A , alist_matrix *a , int NN ) {
  int mm , nn , uu ; 
  
  for ( nn = 1 ; nn <= NN ; nn++ ) {
    for ( uu = 1 , mm = 1 ; mm <= NN ; mm++ ) {
      if ( A[mm][nn] ) {
	a->nlist[nn][uu] = mm ;
	uu ++ ; 
	a->num_mlist[mm] ++ ;
	a->mlist[mm][ a->num_mlist[mm] ] = nn ; 
      }
    }
    uu -- ;
    a->num_nlist[nn] = uu ;
  }
  finish_off_alist ( a ) ;
}

static void cmatrices_2_alist ( unsigned char **A ,  unsigned char **B ,
			       alist_matrix *a , int NN ) {
  int mm , nn , uu , mat , tn ; 
  unsigned char **C ;

  tn = 1 ; 
  for ( C = A , mat = 1 ; mat <= 2 ; mat++ , C = B ) {
    for ( nn = 1 ; nn <= NN ; nn++ , tn++ ) {
      for ( uu = 1 , mm = 1 ; mm <= NN ; mm++ ) {
	if ( C[mm][nn] ) {
	  a->nlist[tn][uu] = mm ;
	  uu ++ ; 
	  a->num_mlist[mm] ++ ;
	  a->mlist[mm][ a->num_mlist[mm] ] = tn ; 
	}
      }
      uu -- ;
      a->num_nlist[tn] = uu ;
    }
    finish_off_alist ( a ) ;
  }
}

static void finish_off_alist ( alist_matrix *a ) {
  int nn , mm , uu ; 
  int tbm = 0  , tbn = 0  ; /* keeps track of the true biggest_m , n */

  for ( mm  = 1 ; mm <= a->M ; mm++ ) {
    uu = a->num_mlist[mm] ; 
    if ( uu > tbm ) tbm = uu ; 
  }
  a->tot = 0 ; 
  for ( nn  = 1 ; nn <= a->N ; nn++ ) {
    uu = a->num_nlist[nn] ; 
    if ( uu > tbn ) tbn = uu ; 
    a->tot += uu ;
  }
  if ( tbm >= a->biggest_num_m ) {
    fprintf (stderr , "eek list overrun %d %d\n",tbm,a->biggest_num_m);
    pause_for_return() ;
  }
  else a->biggest_num_m = tbm ; 
  if ( tbn >= a->biggest_num_n ) {
    fprintf ( stderr , "eek list overrun %d %d\n",tbn,a->biggest_num_n);
    pause_for_return() ;
  }
  else a->biggest_num_n = tbn ; 
}


int cmp_thing ( const void *a , const void *b ) {
  const thing *ma = (const thing*)a ;
  const thing *mb = (const thing*)b ;
  return ( ( ma->val - mb->val ) ) ; 
}

double h2 ( double x ) {
  double tmp ; 
  tmp = x * log ( x ) + ( 1.0 - x ) * log ( 1.0 - x ) ;
  return - tmp / log ( 2.0) ; 
}