default:    status -- ;    fprintf ( stderr , "invalid NORDER %d \n" , c->NORDER ) ;     break ;   }  return status ; }/* needs personalising for fe or depersonalising for bnd */int make_vectors_quick ( data_creation_param *dc , mnc_vectors *v , mnc_all *all ){  int status = 0 ;   int count_s = 0 , count_n = 0 ;   /* the two main options used in 1998 are     dc->gc = 1 (ordinary gaussian channel )     and      dc->gc = 0 (BSC)     other minor options are: dc->nsn  fixed weight noise ; dc->fgc  fixed gaussian noise power.     */  if ( dc->nsn ) { /* a fixed weight vector for all s and n is required */    fixed_wt_cvector ( v->x , dc->nsn , v->nsfrom , v->nnto ) ;    v->count_s = dc->nsn ;     dc->gcx =  (double)  dc->nsn / (double) v->N ;    dc->gcxact = dc->gcx ;  } else { /* other data generation rules */    if ( !dc->ns ) { /* in NMN codes there are source bits to make */       /* normally ns=0 , and nsfrom=1 and nsto=1 */      count_s = random_cvector ( v->x , dc->fs , v->nsfrom , v->nsto ) ;     } else {      fixed_wt_cvector (  v->x , dc->ns , v->nsfrom , v->nsto ) ;      count_s = dc->ns ;     }        if ( dc->fgc ) { /* fixed noise level first, assumed noise level second */      count_n = make_fgaussian_noise_bits_and_fix_biases 	( v->x , v->bias , v->n , dc->fgcx, dc->fgcnf , dc->gcx , v->nnfrom , v->nnto ) ;    }    else if ( dc->gc ) { /* Here is assumed that using bnd */      count_n = make_gaussian_noise_bits_and_fix_biases 	( v->x , v->bias , dc->gcx , v->nnfrom , v->nnto ) ;      dc->fgcx = dc->gcx ;       dc->gcxact = dc->gcx ;  /* ideally this should be set to the 				 actual noise level each time.... */    } else { /* ordinary BSC stuff */      if ( !dc->nn ) {	count_n = random_cvector ( v->x , dc->fn , v->nnfrom , v->nnto ) ; 	dc->gcx = dc->fn ; /* write fn into gcx */      } else { /* this makes the noise fixed weight, independent of what the {s} would be (for MN code) */	fixed_wt_cvector (  v->x , dc->nn , v->nnfrom , v->nnto ) ;	count_n = dc->nn ; /* the desired weight */	dc->gcx =  (double)  dc->nn / (double) v->N ;      }      dc->gcxact =  (double) count_n / (double) v->N ;    }    v->count_s = count_s + count_n ;  }      dc->true_fs = 0.0 ;  dc->true_fn = (double) count_n / (double) dc->M ;       alist_times_cvector_mod2 ( all->a , v->x , v->z ) ;   return status ; }int make_vectors_quickER ( data_creation_param *dc , mnc_vectors *v , mnc_all *all ) /* just makes the calls to the rand number generator */{  int status = 0 ;   int count_s = 0 , count_n = 0 ;   if ( !dc->ns ) {    count_s = random_cvector ( v->x , dc->fs , v->nsfrom , v->nsto ) ;   } else {    fixed_wt_cvector (  v->x , dc->ns , v->nsfrom , v->nsto ) ;  }  if ( dc->gc ) { /* Here is assumed that using bnd */    count_n = make_gaussian_noise_bits_ONLY       ( v->x , v->bias , dc->gcx , v->nnfrom , v->nnto ) ;  }  else {    if ( !dc->nn ) {      count_n = random_cvector ( v->x , dc->fn , v->nnfrom , v->nnto ) ;     } else {      fixed_wt_cvector (  v->x , dc->nn , v->nnfrom , v->nnto ) ;    }  }  return status ; }static void process_vector_fe ( mnc_vectors *v , fe_min_param *p ) {  int m ;   p->true_s = 1 ;   for ( m = 1 ; m <= p->M ; m++ ) {    p->g[m] = ( v->z[m] ) ? 1.0 : -1.0 ;   }}static int make_gaussian_noise_bits_and_fix_biases ( unsigned char *x ,						    double *b,						    double gcx ,						    int lo , int hi ) {  int i , c = 0 ;  double z , p ;  for ( i = lo ; i <= hi ; i ++ ) {    /* make a random normal variate with s.d. 1.0 and mean gcx */    z = gcx + rann() ;     p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcx ) ) ;     if ( z <= 0.0 ) { /* a noise bit is high */      x[i] = 1 ; c ++ ;      b[i] = p ;       if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %f\n" , b[i] ) ;    } else {      x[i] = 0 ;       b[i] = 1.0 - p ;       if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %d %f\n" , x[i] , b[i] ) ;    }  }  return c ; }/* fixed noise - by making n a unit-like vector */ static int make_fgaussian_noise_bits_and_fix_biases ( unsigned char *x ,						    double *b,double *n,						    double fgcx ,						    double fgcnf ,						    double gcx ,						    int lo , int hi ) {  int i , c = 0 ;  double z , p , power=0.0 , factor , count = 0.0 ;  for ( i = lo ; i <= hi ; i ++ ) {    /* make a random normal variate with s.d. 1.0 and mean gcx */    n[i] = rann() ;     power += n[i] * n[i] ;     count += 1.0 ;   }  factor = sqrt( count / power ) ;   fprintf ( stdout , "f=%8.6g -> %5g : " , factor , fgcnf ) ;   fflush(stdout) ;   factor *= fgcnf ; /* boost or reduce noise for importance sampling */  for ( i = lo ; i <= hi ; i ++ ) {    /* make a random normal variate with s.d. 1.0 and mean gcx */    z = fgcx + n[i] * factor ; /* fixed signal to noise level */    p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcx ) ) ; /* gcx - assumed value */     if ( z <= 0.0 ) { /* a noise bit is high */      x[i] = 1 ; c ++ ;      b[i] = p ;       if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %f\n" , b[i] ) ;    } else {      x[i] = 0 ;       b[i] = 1.0 - p ;       if ( b[i] > 0.5 ) fprintf ( stderr , "something wrong %d %f\n" , x[i] , b[i] ) ;    }  }  return c ; }static int make_gaussian_noise_bits_ONLY ( unsigned char *x ,						    double *b,						    double gcx ,						    int lo , int hi ) {  int i , c = 0 ;  double z ;  for ( i = lo ; i <= hi ; i ++ ) {    /* make a random normal variate with s.d. 1.0 and mean gcx */    z = rann() ;   }  return c ; }static void set_up_biases ( double *b , mnc_vectors *v , data_creation_param *dc ) {  int n ;  double lfn1fn = log ( dc->fn / ( 1 - dc->fn ) ) ;  double lfs1fs = log ( dc->fs / ( 1 - dc->fs ) ) ;  for ( n = v->nsfrom ; n <= v->nsto ; n++ ) {    b[n] = lfs1fs ;  }  for ( n = v->nnfrom ; n <= v->nnto ; n++ ) {    b[n] = lfn1fn ;  }}  static void set_up_priors ( double *b , mnc_vectors *v , data_creation_param *dc ) {  int n ;  for ( n = v->nsfrom ; n <= v->nsto ; n++ ) {    b[n] = dc->fs ;  }  for ( n = v->nnfrom ; n <= v->nnto ; n++ ) {    b[n] = dc->fn ;  }}  static int evaluate_feasibility ( data_creation_param *dc )    {  int status = 0 ;/*   compute Shannon capacity and compare with info content  */  dc->h2fs =  h2(dc->fs) ;  dc->h2fn =  h2(dc->fn) ;  dc->rho =  (double ) dc->N / (double ) dc->M  ;  dc->rate = dc->rho * dc->h2fs ;  dc->capacity = 1.0 - dc->h2fn ;  if ( !dc->gc ) {    if ( dc->rate >   dc->capacity ) {      printf ( "task impossible, signed Shannon\n" ) ;       if ( dc->notabovecap ) {status -- ; }    }    if ( status < 0 || dc->verbose )       printf ( "H(s) = %f, capacity of channel = %f.\n" , 	      dc->h2fs , dc->capacity ) ;   }  return status ; }static void dc_defaults ( data_creation_param *dc ) {#include "dc_var6_def.c"}static int process_command ( int argc , char **argv , mnc_all *all ) {  data_creation_param *dc = all->dc ;   fe_min_control *c = all->c ;   fe_min_control *c2 = all->c2 ;   fe_min_control *c3 = all->c3 ;   bnd_control *bndc = all->bndc ;   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 )        {      c->verbose = 1;      dc->verbose = 1;    }    else  if ( strcmp (argv[i], "-VV") == 0 )        {      c->verbose = 2;      dc->verbose = 2;    }#include "fe_var6_clr.c"#include "fed2_var_clr.c"#include "dc_var6_clr.c"#include "bnd_var_clr.c"    else if ( strcmp (argv[i], "-b") == 0 ) {      if ( i + 3 >= argc ) { ERROR2;      }      else {	cs *= sscanf(argv[++i], "%d", &(c->betastyle)); 	cs *= sscanf(argv[++i], "%lf", &(c->beta0)); 	cs *= sscanf(argv[++i], "%lf", &(c->beta1));       }    }    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 , all ) ;  return ( status ) ;}#undef ERROR1#undef ERROR2#undef ERRORREG#define DNT fprintf( fp, "\n        ")#define NLNE  fprintf( fp, "\n")static void print_usage ( char **argv , FILE * fp ,			 mnc_all *all ){  data_creation_param *dc = all->dc ;     fprintf( fp, "Usage: %s ",argv[0]);  fprintf( fp, " [optional arguments]");  DNT;  fprintf( fp, "-V | -VV                 (verbose or very verbose)");   NLNE; fprintf( fp, " Data creation:" ) ; #include "dc_var6_usg.c"  pause_for_return();  NLNE; fprintf( fp, " Inference:" ) ;   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" );   fprintf( fp, "\n");  fe_print_usage ( argv , fp , all ) ;   pause_for_return();  bnd_print_usage ( argv , fp , all ) ;   return ;}static void fe_print_usage ( char **argv , FILE * fp ,			    mnc_all *all ){  fe_min_control *c = all->c ;   fe_min_control *c2 = all->c2 ;   fe_min_control *c3 = all->c3 ;   NLNE; fprintf( fp, " Further free energy minimization stuff: <defaults>");#include "fe_var6_usg.c"#include "fed2_var_usg.c"  fprintf( fp, "\n");  return ;}static void bnd_print_usage ( char **argv , FILE * fp ,			    mnc_all *all ){  bnd_control *bndc = all->bndc ;   NLNE; fprintf( fp, " Belief Net decoder:     <defaults>");#include "bnd_var_usg.c"  fprintf( fp, "\n");  return ;}#undef DNT#undef NLNEstatic void mnc_free ( mnc_all *all ) {/*  data_creation_param *dc = all->dc ;*/  mnc_vectors *vec = all->vec ;   free_cvector ( vec->xo, 1 , vec->N ) ;   free_cvector ( vec->t , 1 , vec->M ) ;   free_cvector ( vec->x , 1 , vec->N ) ;   free_cvector ( vec->y , 1 , vec->M ) ;   free_cvector ( vec->z , 1 , vec->M ) ;   free_dvector ( vec->bias , 1 , vec->N ) ;   free_dvector ( vec->n , 1 , vec->N ) ; }double h2 ( double x ) {  double tmp ;   if ( x <= 0.0 || x>= 1.0 ) tmp = 0.0 ;  tmp = x * log ( x ) + ( 1.0 - x ) * log ( 1.0 - x ) ;    return - tmp / log ( 2.0) ; }/*<!-- hhmts start -->Last modified: Sat Aug 23 16:46:51 1997<!-- hhmts end -->*/