bnd.c

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

/* bnd.c                                      (c) DJCM 95 08 31

   - belief net decoder  -

   Fri 16/2/01 - added gallager A decoding options to this and
   saved as GAd.c

   This code is (c) David J.C. MacKay 1995. It is free software 
   as defined by the free software foundation. 

   ------
   bnd.c:
   ------

   This is a set of routines to do Gallager / MS like  iterative 
 solution of 

 A x = z 

 defined by 

 A, represented in an alist_matrix. 
 a prior on x, given in a vector bias
 an initial condition

 No use is made of the true answer.
 
*/
/*
   How a bit in the A matrix is labelled.

   A one bit is labelled in two ways. 

   1) Its (m,n) pair, i.e. its row (m) and column (n).

   2) Its l and u values, which counts how many 1 bits it is along the row (l)
       and how many 1 bits it is down the column (u).

   A bit can be identified by its (m,l) value, its (n,u) value, or by its (m,n)
   value.

   e.g. the last bit in the bottom row of the (6*12) matrix below 
        has m=6 and n=10, and u=2, and l=4. 

   A = [ 0 1 0 0 0 1 0 1 0 0 0 1 ]
       [ 1 0 0 0 1 0 1 0 0 0 1 0 ]
       [ 0 0 1 1 0 0 0 1 0 0 0 1 ]
       [ 0 1 0 0 1 0 0 0 0 1 1 0 ]
       [ 1 0 0 1 0 1 0 0 1 0 0 0 ]
       [ 0 0 1 0 0 0 1 0 1 1 0 0 ]

   Notice that all rows in this matrix have equal weight lmax=4 and all 
   columns have equal weight umax=2 .

   A = [   1       1   1       1 ]
       [ 1       1   1       1   ]
       [     1 1       1       1 ]
       [   1     1         1 1   ]
       [ 1     1   1     1       ]
       [     1       1   1 1     ]

*/
/* CONTENTS

void bnd_allocate ( bnd_param *p , bnd_control *c ) {
void bnd_defaults ( bnd_param *p , bnd_control *bndc ) {
void bnd_free ( bnd_param *p , bnd_control *c ) {
void   bnd_load_dqc ( bnd_param *p , bnd_control *c ) {
void bnd_fprint_state ( bnd_param *p , bnd_control *c ) {
int bndecode ( bnd_param *p , bnd_control *c ) {
void bnd_score_state ( bnd_param *p ) {
void horizontal_pass ( bnd_param *p  , bnd_control *c  ) 
void vertical_pass ( bnd_param *p , bnd_control *c  ) 
void vertical_pass2 ( bnd_param *p , bnd_control *c  ) 
void bnd_tweak_prior ( bnd_param *p , bnd_control *c  ) 
void fudges ( double *d , double f )   NOT USED

*/


#include "./ansi/r.h"
#include "./ansi/mynr.h"
#include "./ansi/cmatrix.h"
#include "./bnd.h"

/* malloc memory storage space */
void bnd_allocate ( bnd_param *p , bnd_control *c ) {
  int M = p->M ; 
  int N = p->N ;
  int m ; 

/* get z, xo to point to a vector allocated elsewhere; x is not used except in the printstate diagnostic*/
  p->t = cvector ( 1 , M ) ; 

/*  p->bias = dvector ( 1 , N ) ; */ /* cut out in mnc7 onwards; this is
                                         allocated outside -- NB this makes
					 this code incompatible with mnc6.c */
  p->lbias = dvector ( 1 , N ) ; /* logit of bias vector */
  p->q1 = dvector ( 1 , N ) ; 

  p->dpf = dmatrix ( 1 , M , 0 , p->a->biggest_num_m + 1 ) ;
  p->dpr = dmatrix ( 1 , M , 0 , p->a->biggest_num_m + 1 ) ;

  p->pc0 = dmatrix ( 1 , N , 1 , p->a->biggest_num_n ) ;
  p->pc1 = dmatrix ( 1 , N , 1 , p->a->biggest_num_n ) ;

  p->dqc = dmatrix ( 1 , N , 0 , p->a->biggest_num_n + 1 ) ;

  p->qc0 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qc1 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qt0 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qt1 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qb0 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;
  p->qb1 = dvector ( 0 , p->a->biggest_num_n + 1 ) ;

  for ( m = 1 ; m <= M ; m ++ ) {
    p->dpf[m][0] = 1.0 ; 
    p->dpr[m][p->a->num_mlist[m]+1 ] = 1.0 ; 
  }
  /* vertical passes must be initialized too */
}

void bnd_defaults ( bnd_param *p , bnd_control *bndc ) {
#include "bnd_var_def.c"
}

void bnd_free ( bnd_param *p , bnd_control *c ) {
  int M = p->M ; 
  int N = p->N ;

  free_cvector ( p->t , 1, M ) ; 

  free_dvector ( p->lbias , 1 , N ) ;   /* changed in mnc7 */
  free_dvector (   p->q1 , 1 , N ) ; 

  free_dmatrix ( p->dpf , 1 , M , 0 , p->a->biggest_num_m + 1 ) ;
  free_dmatrix ( p->dpr , 1 , M , 0 , p->a->biggest_num_m + 1 ) ;

  free_dmatrix ( p->pc0 , 1 , N , 1 , p->a->biggest_num_n ) ;
  free_dmatrix ( p->pc1 , 1 , N , 1 , p->a->biggest_num_n ) ;

  free_dmatrix ( p->dqc , 1 , N , 0 , p->a->biggest_num_n + 1 ) ;

  free_dvector ( p->qc0 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qc1 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qt0 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qt1 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qb0 , 0 , p->a->biggest_num_n + 1 ) ;
  free_dvector ( p->qb1 , 0 , p->a->biggest_num_n + 1 ) ;

}

void   bnd_load_dqc ( bnd_param *p , bnd_control *c ) {
  int u , n , N=p->N ; 
  int junk=-1 ;
  
  for ( n = 1 ; n <= N ; n++ )  {
    for ( u = 1 ; u <= p->a->biggest_num_n ; u++ )  {
      p->dqc[n][u] = 1.0 - 2.0 * p->bias[n] ;
    }
  }
  if ( c->vpl ) { /* using vertical pass log version   Mon 10/6/02 */
    for ( n = 1 ; n <= N ; n++ )  { /* logit */
      p->lbias[n] = safelog( ( p->bias[n] / ( 1.0 - p->bias[n] ) ) , &junk ) ;
    }
  }
  
}
double safelog ( double x , int *c ) {
  if (x>0.0) { return log(x) ; }
  else {
    if ( *c< 0 ) {
      fprintf (stderr , "SL%d: can't take log of %g\n" , - *c , x ) ;
    } else if ( *c == 0 ) {
      fprintf (stderr , "L" ) ;
      *c = *c + 1 ;
    } /* else if *c>0 do nothing */
    return -2.0 ; /* reasonable replacement for log(0) :-) */
  }
}

int bndecode ( bnd_param *p , bnd_control *c ) {

  bnd_load_dqc ( p , c ) ; /* sets all dqc entries to the prior to start off */

  for ( c->loop = 1 , p->not_decoded = 1 ; 
       ( c->loop <= c->loops ) && p->not_decoded ; 
       c->loop ++ ) {

    horizontal_pass ( p , c ) ;  /* turns dqc into dpr and dpf 
			     then reads out dpc -> pc0 and pc1 */
    if ( c->vpl )  /* use log version of vertical pass */
      vertical_pass2 ( p , c ) ;
    else
      vertical_pass ( p , c ) ;    /* turns pc0 and pc1 into qt0, qb0, etc, 
				      and works out qc0, qc1 -> dqc.
				      Also works out the latest marginal belief */
    bnd_score_state ( p ) ;  /* takes belief vector, turns it into an
			     x, and runs it through Ax =? z .
			     If equal then decoding is halted.
			     */
    if ( c->writelog ) 
      bnd_fprint_state ( p , c ) ; 
    /* added 1997 02 */
    if ( c->wis && !(c->loop%c->wis) ) 
      bnd_write_state ( p , c ) ;  /* write the whole intermediate state */


    if ( c->dotweak && ( c->loop >= c->dotweak ) ) bnd_tweak_prior (p,c) ;

  }
  c->loop -- ; 
  printf( "\tBND:\tviols %3d  recon %3d  its %2d\n", p->count_viol , p->count_high  , c->loop ) ; 

  return ( p->count_viol ) ; /* zero if success */
}

void bnd_write_state ( bnd_param *p , bnd_control *c ) {
  FILE *fp ;
  int n ;
  double *q1 = p->q1 ; 
  int l = c->loop ; 
  char junk[1000] ; 
 
  sprintf ( junk , "%s.%02d" , c->isfile , l ) ;
  fp = fopen ( junk , "w" ) ;
  fprintf ( stdout , "writing intermediate state %d to %s, grey = %d\n" , l , junk , c->wisgrey ) ; 
  if ( fp ) {
    if (c->wisgrey) 
      for ( n = 1 ; n <= p->N ; n++ )  {
	fprintf ( fp ,  "%-3d " , (int)(255.0 * q1[n]) ) ; 
	if (!(n%c->wisline)) fprintf ( fp ,  "\n" ) ; 
      }
    else 
      for ( n = 1 ; n <= p->N ; n++ ) {
	fprintf ( fp ,  "%-1d " , p->xo[n] ) ; 
	if (!(n%c->wisline)) fprintf ( fp ,  "\n" ) ; 
      }
    fprintf ( fp ,  "\n" ) ; 
    fclose ( fp ) ; 
  }
  else 
    fprintf ( stderr , "can't open %s \n" , junk ) ;
}

void bnd_fprint_state ( bnd_param *p , bnd_control *c ) {
  FILE *fp ;
  int n ;
  double *q1 = p->q1 ; 

  fp = fopen ( c->logfile , "a" ) ;
  if ( fp ) {
    for ( n = 1 ; n <= c->writelog ; n++ ) 
      fprintf ( fp ,  "%-2.0f " , 99.0 * q1[n] ) ; 
    fprintf ( fp ,  "\n" ) ; 
    for ( n = 1 ; n <= c->writelog ; n++ ) 
      fprintf ( fp ,  "%-2d " , p->x[n] ) ; 
    fprintf ( fp ,  "\n" ) ; 
    fclose ( fp ) ; 
  }
  else 
    fprintf ( stderr , "can't open %s \n" , c->logfile ) ;
}

void bnd_score_state ( bnd_param *p ) {
  double *q1 = p->q1 ;
  double *bias = p->bias ;
  int N = p->N ;
  int M = p->M ;
  int m , n ;

  
  for ( p->count_high = 0 , n = 1 ; n <= N ; n++ )  {
    p->xo[n] = ( q1[n] >= 0.5 ) ? 1 : 0 ;
    p->count_high += p->xo[n] ; /* seems that this is not being counted right */
  }
	if ( 0 ) {
	  printf ( "Proposed answer has %d high / %d bits\n" , p->count_high , N) ;
	}

  alist_times_cvector_mod2 ( (p->a) , p->xo , p->t ) ; 

  for ( p->count_viol = 0 , m = 1 ; m <= M ; m ++ ) {
    if ( p->t[m] != p->z[m] ) p->count_viol ++ ; 
  }

  /* Tue 4/6/02 addition */
  p->loglike = 0.0 ; /* failures are reported as loglike of zero */
  if ( p->count_viol == 0 ) {
    p->not_decoded = 0 ;
    /* found a valid decoding, let's get its score. Typical bias is 0.1. -> if(0),log(1-b) */
    for ( n = 1 ; n <= N ; n++ )  {
      p->loglike += p->xo[n] ?   log ( bias[n] ) :  log ( 1.0 - bias[n] );
    }
  }

}


void horizontal_pass ( bnd_param *p  , bnd_control *c  ) 
{
  int l , u , umax ; 
  int m , n ;
  double *dpf ;
  double *dqc ;
  double dpc ;