nn.c

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

	   nc->EDHtr   , safe_divide ( nc->EDHtr , nc->tot_tr ) ,
	   nc->EHwb_tr , safe_divide ( nc->EHwb_tr , nc->train_n ) ) ; 
  fprintf ( fp , "%-11.6g %-9d %-11.6g %-6d %-9.4g " ,
	   nc->ET ,  
	   nc->EDHte   , safe_divide ( nc->EDHte , nc->tot_te )  ,
	   nc->EHwb_te , safe_divide ( nc->EHwb_te , nc->test_n ) ) ;
  fprintf ( fp , "%-9d %-11.6g %-6d %-9.4g " ,
	   nc->itEH ,   safe_divide ( nc->itEH , nc->tot_dec  )  ,
	   nc->itEHwb , safe_divide ( nc->itEHwb , nc->decodn ) ) ; 
  if ( extras ) {
    fprintf ( fp , "%-11.6g " , nc->M ) ; 
    fprintf ( fp , "%-11.6g " , nc->QW ) ;  
    fprintf ( fp , "%-9.4g " , nc->tol ) ; 
    fprintf ( fp , "%-2d " , nc->loop ) ; 
    fprintf ( fp , "%-3d " , nc->macarg.its ) ; 
    fprintf ( fp , "%-5d " , nc->num_derivatives ) ; 
    for ( c = 1 ; c <= nc->RC ; c++ ) {
      fprintf ( fp , "%-9.4g " ,  (double)(nc->ninrc[c]) / nc->EW[c] ) ; 
    }
    for ( c = 1 ; c <= nc->RC ; c++ ) {
      fprintf ( fp , "%-9.4g " , nc->alpha[c] ) ; 
    }
  }
  fprintf ( fp , "\n" ) ; 
}

static double safe_divide ( int a , int n ) {
  return ( ( n > 0 ) ? ( (double)(a) / (double)(n) ) : -1.0 ) ;
}

void    nn_say_headings ( FILE *fp ,  mnc_net_control *nc , int extras ) 
{
  int c , d ; 
  fprintf ( fp , "#-Tr:-logP  bit_ers   3::frac_bit blk_es 5::fr_blk "  ) ; 
  fprintf ( fp , "6-Te:-logP  bit_ers   8::frac_bit blk_es 10:fr_blk "  ) ; 
  fprintf ( fp ,             "11itv_ers 12:frac_bit blk_es 14:fr_blk "  ) ; 
  if ( extras ) {
    c = 15 ; 
    fprintf ( fp , "%2d-M        %2d-QW       " , c , c + 1 ) ; c+=2 ; 
    fprintf ( fp , "%2d-tol    lp its der's " , c ) ; c+=4 ;
    for ( d = 1 ; d <= nc->RC ; d++ ) {
      fprintf ( fp , "%2d-^alph%d " , c , d ) ; c++ ; 
    }
    for ( d = 1 ; d <= nc->RC ; d++ ) {
      fprintf ( fp , "%2d-alpha%d " , c , d ) ; c++ ;
    }
  }
  fprintf ( fp , "\n" ) ; 
}

double nn_objective ( double *w , void *arg ) 
{
  mnc_all *all = ( mnc_all * ) arg ; 
  mnc_net *net = all->net ; 
  mnc_net_control *nc = all->nc ; 

  net->w = w ; 

  ran_seed ( nc->trseed ) ; 

  all->dc->message = 0 ; /* ensure seed left alone */
  nc->ED = nn_error_on ( net , all , nc->train_n ) ; 
  nc->QW = nn_weight_energy ( net , nc ) ; 
  return ( nc->M = nc->ED + nc->QW ) ; 
}

/*
   this bit is split off from the rest so as to make it easy to evaluate 
   a test error
*/
double nn_error_on ( mnc_net *net , mnc_all *all , int n ) 
{
  double ED = 0.0 ; 
  int evalH = all->nc->evalH , EDH=0 , EHwb = 0 , tmpi ; 

  net->x = all->vec->z ; /* input  = output of A x + y = z */
  net->t = all->vec->x ; /* target = original x vector */

  for (  ; n >= 1 ; n-- ) {
    make_vectors_quick ( all->dc , all->p , all->vec) ; 
    nn_forward_pass ( net ) ;
    ED += nn_error_signal ( net ) ; 
    if ( evalH ) {
      EDH  += ( tmpi = nn_hard_errors ( net ) ) ; 
      EHwb += ( tmpi > 0 ) ? 1 : 0 ; 
    }
  }
  if ( evalH ) {
    all->nc->EDH = EDH ; 
    all->nc->EHwb = EHwb ;
  }
  return ( ED ) ; 
}

  /* 
     How to compute the gradient.

     1: make the net's current w vector point where demanded, and g too. 
     2: initialize the signal generating machine
     3: loop-
       - generate a signal, and corresponding encoded string
       - these are traditionally in p->s[] and in p->g[m] (implicitly), or
           in v->r[]. these become the targets and the inputs respectively
	   of the net. Now they are in p->x amd p->z. 
       - do a forward pass
       - score reconstruction, make errors
       - backward pass
     4: add regularizer
     5: done
                                                       */
void nn_derivative ( double *w , double *g , void *arg ) 
{
  mnc_all *all = ( mnc_all * ) arg ; 
  mnc_net *net = all->net ; 
  mnc_net_control *nc = all->nc ; 
  int n ; 
  nc->ED = 0.0 ; 

  net->w = w ; 
  net->g = g ; 

  net->x = all->vec->z ; /* input  = output of A x + y = z */
  net->t = all->vec->x ; /* target = original x vector */

  ran_seed ( nc->trseed ) ; 

  set_dvector_const ( g , 1 , net->K , 0.0 ) ;
  all->dc->message = 0 ; /* ensure seed left alone */
  for ( n = nc->train_n ; n >= 1 ; n-- ) {
    if ( make_vectors_quick ( all->dc , all->p , all->vec) < 0 ) exit (0) ; 
    nn_forward_pass ( net ) ;
    nc->ED += nn_error_signal ( net ) ; 
    nn_backward_pass ( net ) ; 
  }
  nn_add_regularizer ( net , all->nc ) ; 
  /* g now contains the gradient */
  nc->num_derivatives ++ ; 
}

void nn_forward_pass ( mnc_net *net ) {
  double *a = net->a ; 
  double *ca = net->ca ; 
  double *y = net->y ; 
  double *w = net->w ;
  unsigned char *x = net->x ; 

  int o , i ; 

  w++ ; 
  for ( o = 1 ; o <= net->O ; o ++ ) {
    a[o] = *w ; w ++ ; 
    for ( i = 1 ; i <= net->I ; i++ , w++ ) {
      if ( x[i] ) a[o] += *w ;
    }
    y[o] = nn_act_fun ( 10 , &a[o] , &ca[o] ) ; 
  }
}

/* 
   This routine puts in the hit list all bits with 
   probability greater than nc->hitlist_thresh, 
   or else the top nc->hitlist_n bits.
   Or rather, at least one bit, and up to that many, stopping 
   at the threshold nc->hitlist_low
*/
int nn_hitlist ( mnc_net *net , mnc_all *all ) {
  double *y = net->y ; 
  unsigned char *xpartial = all->vec->xpartial ; 
  int *touches = all->vec->touches ; 
  mnc_net_control *nc = all->nc ; 
  int count = 0 , u , N = net->O , o , n , min_hit=1 , max_hit=N ; 
  thing *my ; 

  for ( o = 1 ; o <= net->O ; o ++ ) {
    xpartial[o] = 0 ; 
  }

  if ( nc->decverbose > 1 ) 
    printf ( "making hitlist\n" ) ; 

  if ( nc->hitlist_policy == 1 ) { /* hope to get away without sorting */
    for ( o = 1 ; o <= N ; o ++ ) {
      if ( y[o] > nc->hitlist_thresh ) {
	xpartial[o] = 1 ; 
	count ++ ; 
      }
    }
    if ( nc->decverbose > 1 ) { printf ( "%d hits\n" , count ) ; }
    if ( count > 0 ) { 
      return count ;
    }
    /* else no candidates, time to sort */
  }

  /* need an ordering machine */
  my = (thing * ) malloc ( N * sizeof(thing)) ; 
  for ( n = 0 ; n <= N-1 ; n ++ ) {
    my[n].index = n+1 ; 
    my[n].val = (int) ( 100000000.0 * y[n+1] )  ;
  }
  qsort ( my , (unsigned int) (N) , sizeof(thing) , cmp_thing ) ; 
  if ( nc->hitlist_policy == 1 ) max_hit = 1 ;
  else max_hit = nc->hitlist_n ;

  if ( nc->decverbose > 1 ) 
    printf ( "sorted buggers:\n" ) ;
  /* 
     run down the list until the max number of required hits is reached
     or until we hit the bottom or until
     the min number of hits is reached AND the score is below _low.

     only add a dude to the list if he's not been touched twice already 
     */
  for ( u = 0 , n = N-1 ; n >= 0 ;  n-- )  {
    o = my[n].index ; 
    if ( nc->decverbose > 1 ) 
      printf ( "%d: %15d ( = %9.5g )\n" , n+1 , my[n].val , y[o] ) ; 
    if ( touches[o] < 2 ) {
      xpartial[o] = 1 ; 
      u++ ; 
      touches[o] ++ ; 
    } else { touches[o] += 10 ; } /* this leaves a record of a strong desire */
    if ( ( u >= max_hit )
	|| ( ( u >= min_hit ) && ( y[o] < nc->hitlist_low ) )
	) break ; 
  } 
  if ( u < min_hit ) { /* we're not making enough changes ! */
    if ( nc->decverbose > 0 ) 
      printf ( " stuck (%d), resetting touches to zero\n" , u ) ; 
    set_ivector_const ( touches , 1 , N , 0 ) ;
  }
  free((char * ) (my ) ) ; 
  return 0 ;
}

int nn_flip_and_score ( mnc_net *net , mnc_all *all ) {
  FILE *fp ; 
  unsigned char *x = all->vec->x ; 
  unsigned char *xpartial = all->vec->xpartial ; 
  mnc_net_control *nc = all->nc ; 

  int o , N = net->O , flipped = 0  , corrected = 0  , spoilt = 0 , remain = 0 ;
  for ( o = 1 ; o <= N ; o ++ ) {
    if ( xpartial[o] == 1 ) {
      flipped ++ ;
      if ( x[o] == 1 ) { x[o] = 0 ; corrected ++ ; }
      else { x[o] = 1 ; spoilt ++ ; remain ++ ; }
    } else if ( x[o] == 1 ) { remain ++ ; }
  }
  
  if ( nc->decverbose > 0 ) {
    fprintf ( stdout , "%-4d %-4d %-4d %-4d\n" ,
	     flipped   , corrected   , spoilt  , remain ) ; 
  }
  if ( nc->decwrite ) {
    fp = fopen ( nc->decfile , "a" ) ; 
    if ( fp ) {
      fprintf ( fp , "%-4d %-4d %-4d %-4d\n" ,
	       flipped   , corrected   , spoilt  , remain ) ; 
      fclose ( fp ) ; 
    } else fprintf ( stderr , "couldn't open %s\n" , nc->decfile ) ; 
  }
  return ( remain ) ; 
}

int nn_hard_errors ( mnc_net *net ) {
  double *y = net->y ; 
  unsigned char *t = net->t ; 
  unsigned char *h = net->h ; 
  int E = 0 ; 
  int o ; 
  double thresh = net->thresh ; 

  for ( o = 1 ; o <= net->O ; o ++ ) {
    h[o] = ( y[o] >= thresh ) ? 1 : 0 ; 
    E += h[o]^t[o] ; 
  }
  return E ; 
}

double nn_error_signal ( mnc_net *net ) {
  double *y = net->y ; 
  double *e = net->e ; 
  double *a = net->a ; 
  double *ca = net->ca ; 
  unsigned char *t = net->t ; 
  double E = 0.0 ; 
  int o ; 

  for ( o = 1 ; o <= net->O ; o ++ ) {
    if ( t[o] ) { e[o] = 1.0 - y[o] ; E -= a[o] ; }
    else { e[o] = - y[o] ; E -= ca[o] ; }
  }
  return E ; 
}              /* by definition E = - log prob (data ) */
               /* dE/da = -(t-y)  */

void nn_backward_pass ( mnc_net *net ) {
  double *e = net->e ; 
  double *g = net->g ;
  unsigned char *x = net->x ; 

  int o , i ; 

  g++ ; 
  for ( o = 1 ; o <= net->O ; o ++ ) {
    *g -= e[o] ; g ++ ; 
    for ( i = 1 ; i <= net->I ; i++ , g++ ) {
      if ( x[i] ) *g -= e[o] ; 
    }
  }
}

double nn_weight_energy ( mnc_net *net , mnc_net_control *nc ) {
  double *alpha = nc->alpha ; 
  double *w = net->w ;
  double *EW = nc->EW ; 
  double Q = 0.0 ; 
  int  i , c ; 
  int *rc = net->rc ; 

  for ( c = 1 ; c <= nc->RC ; c++ ) {
    EW[c] = 0.0 ; 
  }
  w++ ;
  for ( i = 1 ; i <= net->K ; i++ , w++ ) {
    c = rc[i] ;
    EW[c] += (*w) * (*w) ;
  }
  for ( c = 1 ; c <= nc->RC ; c++ ) {
    Q += 0.5 * alpha[c] *  EW[c] ; 
  }
  return Q ; 
}

void nn_add_regularizer ( mnc_net *net , mnc_net_control *nc ) {
  double *alpha = nc->alpha ; 
  double *g = net->g ;
  double *w = net->w ;
  int *rc = net->rc ; 
  int  i ; 

  g++ ; w++ ;
  for ( i = 1 ; i <= net->K ; i++ , g++ , w++ ) {
    *g += alpha[rc[i]] * (*w) ; /* gradient points up hill */
  }
}

/* this routine computes the activity and records the values of log y
   and log (1-y) in the two locations 
   (handy for evaluating the objective function) 
   [from glo4.c]
*/
double  nn_act_fun( int number , double *act , double *cact ) 
{
  double tmpd , tmpd2 ;

  switch(number){
  case(1):
    return tanh(*act);
    break;
  case(2):
    return 1.0/(1.0+exp(-2.0 * (*act)));
    break;
  case(3):
  case(10):/* 10 is the choice of af for classification problems'
	      output units. */
    /* anti[0]=1.0/(1.0+exp(*act)); */
    if ( *act > 0.0 ) {
      tmpd = 1.0 + exp ( - *act ) ;
      tmpd2 = - log ( tmpd ) ;
      *cact = - (*act) + tmpd2 ;
      *act = tmpd2 ; 
      return 1.0/tmpd ;
    } else {
      tmpd = 1.0 + exp ( *act ) ;
      *cact = - log ( tmpd ) ;
      *act += *cact ; 
      return 1.0 - 1.0/tmpd ;
    }
    break;
  default:
  case(0):
    return  *act;
    break;
  }
}

double  nn_d_act_fun(int number, double act , double xact , double antixact)
{
  switch(number){
  case(1):
    return  1.0-xact*xact;
    break;
  case(2):
    return 2.0*(xact)*(antixact); /* check this */
    break;
  case(3):
    return (xact)*(antixact); /* check this */
    break;
  case(4):
  case(5):
    return antixact ;
    break;
  case(0):
  case(10):       /* 10 is the choice of af for classification problems'
		     output units. It is a 0/1 sigmoid going forward
		     but going backward it is linear because 
		     of the G-info objective function */
  default:
    return  1.0;
    break;
  }
}


/*
<!-- hhmts start -->
Last modified: Sat Apr 15 14:00:18 1995
<!-- hhmts end -->
*/