/*   gallager - software to simulate decoding of low-density parity-check codes     Copyright   (c) 2002   David J.C. MacKay     This program is free software; you can redistribute it and/or modify    it under the terms of the GNU General Public License as published by    the Free Software Foundation; either version 2 of the License, or    (at your option) any later version.    This program is distributed in the hope that it will be useful,    but WITHOUT ANY WARRANTY; without even the implied warranty of    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    GNU General Public License for more details.    You should have received a copy of the GNU General Public License    along with this program; if not, write to the Free Software    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA    GNU licenses are here :    http://www.gnu.org/licenses/licenses.html    Author contact details are here :    http://www.inference.phy.cam.ac.uk/mackay/    If you find this software useful, please feel free to make a donation to    support David MacKay's research group.*//*   gallager.c   mncN.c   mnc6.c                                    (c) DJCM 94 07 04                                                      94 10 26						      94 10 31						      95 01 02						      95 02 23 (from mnc2 to mnc3)						      95 04 01 continued work on mnc3						      95 05 01 to mnc4						      95 05 27 to mnc5						      95 08 31 to mnc6						      95 11 04 to mnc7						      95 11 16 to mncN						      97 07 23 added noise-writing                                                      98 10 12 to gallager						      98 11 19 added bsc option again   - This simulator can decode LDPC codes and other closely related codes,     for example MN codes, and LDPC codes with some bits punctured.     It does binary codes only.     It can simulate Gaussian channels and binary symmetric channels,     and BSCs with two different known noise levels.						         - gallager writes files rationally like RA.c      added features: spotting of low weight errors      and logging of det/undet and histograms.      So that the output file becomes unnecessary.   - mncN = quick hack to do all bits as if they are noise bits.   - free energy minimization for the MacKay-Neal Error Correcting Code -    - or belief net decoder                                              -   - mnc7 includes option of gaussian channel                           -   - space is made for the fe min or for the bnd, then the scoring      routine gets to look at the results via pointers in the mnc_vector     structure. Thus scoring is divorced from fe.      crucial var files: both fe_var6 and dc_var6   ======================================================================   This code is (c) David J.C. MacKay 1994, 1995. It is free software    as defined by the free software foundation.    What this program does:      reads in a matrix A in alist format.     do {       generates a  random vector x,        shoves it through A, adds y if appropriate, -> z = Ax + y       do {         gives b, g (z) to the minimization machine, or A,z to the bnd         checks the answer. 	 next strategy       }       until ( success OR no more strategies )     }     while ( not many failures )*/#include "./ansi/r.h"#include "./ansi/rand2.h"#include "./ansi/mynr.h"#include "./ansi/cmatrix.h"#include "./ansi/macopt.h"#include "./thing_sort.h"#include "./fe6.h"#include "./bnd.h"#include "./GAd.h" /* added purely for structure compatibility *//* this must come last: */#include "./mnc6.h"/* static void totalline ( FILE * , mnc_all * , data_creation_param * ); */static void dc_defaults ( data_creation_param * ) ; static int    process_command ( int , char **  , mnc_all * ) ; static void   print_usage ( char ** , FILE * ,  mnc_all *  );static void   fe_print_usage ( char ** , FILE * ,  mnc_all *  );static int make_sense ( data_creation_param * , mnc_all * ) ; static int make_sense_c ( data_creation_param * , fe_min_control * ) ; static int make_space ( mnc_all * , mnc_vectors * ) ;static int make_norder ( alist_matrix * , fe_min_control * ) ;static int make_vectors_quick ( data_creation_param * , mnc_vectors * ,  mnc_all * ) ;static int make_vectors_quickER ( data_creation_param * , mnc_vectors * ,  mnc_all * ) ;static void process_vector_fe ( mnc_vectors *v , fe_min_param *p ) ;static int make_fgaussian_noise_bits_and_fix_biases ( unsigned char * ,						    double *,double *,						    double  , double , 						    double  , int , int  ) ;static int make_gaussian_noise_bits_and_fix_biases ( unsigned char * ,						    double *,						    double  ,						    int , int  ) ;static int make_gaussian_noise_bits_ONLY ( unsigned char * ,						    double *,						    double  ,						    int , int  ) ;static void set_up_biases ( double * , mnc_vectors * , data_creation_param * ) ; static void set_up_priors ( double * , mnc_vectors * , data_creation_param * ) ; static int score ( data_creation_param * , mnc_vectors  * , mnc_all * ) ; static int evaluate_feasibility ( data_creation_param * ) ;static double bern ( int  , int  , double * , double * ,  double * , double );static void write_histo ( FILE * , mnc_all * ) ;static void snappyline ( data_creation_param * ) ;static void   finalline (  FILE * , mnc_all * , data_creation_param * , int ) ;static void mnc_free ( mnc_all * ) ; static int check_alist_MN ( alist_matrix * , mnc_vectors *vec , data_creation_param * ) ; static int hook_mnc_vec_to_bnd ( bnd_param *p , mnc_vectors *vec ) ;static int hook_mnc_vec_to_fe ( fe_min_param *p , mnc_vectors *vec ) ;static void bnd_print_usage ( char **argv , FILE * fp , mnc_all *all );static double h2 ( double ) ;void   main ( int , char ** ) ;/*        MAIN                     */void main ( int argc, char *argv[] ){  FILE   *fp  ;   fe_min_param        p    ;  fe_min_control      c    ;  fe_min_control      c2   ; /* backup strategy */  fe_min_control      c3   ; /* backup strategy */  bnd_control         bndc ;   bnd_param           bndp ;  data_creation_param dc   ;  mnc_vectors         vec  ;   mnc_all             all  ;  all.dc   = &dc ;   all.p    = &p ;   all.c    = &c ;   all.bndp = &bndp ;   all.bndc = &bndc ;   all.c2   = &c2 ;   all.c3   = &c3 ;   all.vec  = &vec ;   all.a    = &(p.a) ;   bndp.a   = &(p.a) ;   fe_defaults_p ( &p ) ;   fe_defaults_c ( &c ) ;   fe_defaults_c ( &c2 ) ;   fe_defaults_c ( &c3 ) ;   dc_defaults   ( &dc ) ;   bnd_defaults  ( &bndp , & bndc ) ;   if ( process_command (argc, argv, &all ) < 0 ) exit (0) ;  if ( make_sense ( &dc , &all ) < 0 ) exit (0) ;  fprintf(stderr,"mncN N0=%d, M0=%d, N=%d, M=%d, x=%6.3g, fn=%6.3g, fs=%6.3g, nn=%d, ns=%d\n",	  dc.N , dc.M , vec.N , vec.M , dc.gcx , dc.fn , dc.fs , dc.nn , dc.ns) ;  fflush(stderr);  if ( dc.use_fe && c.second ) {    fprintf(stdout,"opt=%2d:%2d:%2d ",	     c.opt , c.gamma , c.metric ) ;    fprintf(stdout,"opt2=%2d:%2d:%2d ",	     c2.opt , c2.gamma , c2.metric ) ;    fprintf(stdout,"opt3=%2d:%2d:%2d ",	     c3.opt , c3.gamma , c3.metric ) ;  }  if ( make_space ( &all , &vec ) < 0 ) exit (0) ;  if ( read_allocate_alist ( &(p.a) , dc.afile ) < 0 ) exit (0) ;   if ( check_alist_MN ( &(p.a) , &vec , &dc ) < 0 ) exit (0) ;   if ( dc.use_fe >= 1 ) {    c.mega_needed = ( c.opt == 3 || 		   ( c.second && c2.opt == 3 ) || 		   ( c.third && c3.opt == 3 ) ) ;    hook_mnc_vec_to_fe (  &p ,  &vec ) ;    fe_allocate2 ( &p , c.mega_needed ) ;     set_up_biases ( p.bias , &vec , &dc ) ;     if ( make_norder ( &(p.a) , &c ) < 0 ) exit ( 0 ) ;   }  if ( dc.use_bnd >= 1 ) {    hook_mnc_vec_to_bnd (  &bndp ,  &vec ) ;    bnd_allocate ( &bndp , &bndc ) ;     set_up_priors ( bndp.bias , &vec , &dc ) ;   }  if ( evaluate_feasibility ( &dc ) < 0 ) exit ( 0 ) ;   if ( c.writelog ) {    fp = fopen ( c.logfile , "w" ) ;    if ( !fp ) {      fprintf ( stderr , " couldn't open logfile %s\n" , c.logfile ) ;       c.writelog = 0 ;     } else fclose (fp ) ;   }  if ( bndc.writelog ) {    fp = fopen ( bndc.logfile , "w" ) ;    if ( !fp ) {      fprintf ( stderr , " couldn't open logfile %s\n" , bndc.logfile ) ;       bndc.writelog = 0 ;     } else fclose (fp ) ;   }  if ( dc.error_log ) {    fp = fopen ( dc.error_logfile , "w" ) ;    if ( !fp ) {      fprintf ( stderr , " couldn't open logfile %s\n" , dc.error_logfile ) ;       dc.error_log = 0 ;     } else fclose (fp ) ;   }  if ( dc.noise_sum ) {    fp = fopen ( dc.noise_sumfile , "w" ) ;    if ( !fp ) {      fprintf ( stderr , " couldn't open logfile %s\n" , dc.noise_sumfile ) ;       dc.noise_sum = 0 ;     } else fclose (fp ) ;   }  /*      MAIN LOOP                 */  ran_seed ( dc.vseed ) ;   dc.message = 1 ;  if ( dc.skip ) {    fprintf ( stderr , "skipping %d", dc.skip ) ;     for (  ; dc.message <= dc.skip ; dc.message ++ ) {      make_vectors_quickER ( &dc , &vec , &all ) ;       if ( !(dc.message % 200) ) { fprintf ( stderr , "." ) ; fflush ( stderr ) ; }    }fprintf ( stderr , "\n" ) ;   }  for ( ;       ( dc.message <= dc.MESSAGE ) &&        ( ( dc.failures==0 ) || ( dc.failcount < dc.failures ) ) ;        dc.message ++ ) {    snappyline( &dc ) ;    fflush(stdout) ;     make_vectors_quick ( &dc , &vec , &all ) ; /* see this file */    if ( dc.use_fe >=1 ) {      process_vector_fe ( &vec , &p ) ;       all.controller = 1 ;       if ( ( fe_min ( &p , &c ) > 0 ) && ( c.second ) ) {	all.controller ++ ;	if ( ( fe_min ( &p , &c2 ) > 0 ) && ( c.third ) ) { 	  all.controller ++ ;	  fe_min (  &p , &c3 ) ; 	}      }    } else {      all.controller = 0 ;       bndecode ( &bndp , &bndc ) ;    }    if ( score ( &dc , &vec , &all ) < 0 ) exit ( 0 ) ; /* score also writes error logs and keeps track of nasty noise patterns */    if ( dc.verbose > 0 ) finalline ( stdout , &all , &dc , 0 ) ;    if ( c.printout ) { /* append */      fp = fopen ( c.outfile , ((c.outappend)? "a":"w") ) ;      if( !fp ) {	fprintf( stderr, "No such file: %s\n", c.outfile ) ;	finalline ( stderr , &all , &dc , 0 ) ;      }      else 	{	finalline ( fp , &all , &dc , 0 ) ;	fclose ( fp ) ;      }    }    /* if this is going to be the last loop, or the number of messages is 30,60,... */    if ( ( !( ( dc.message+1 <= dc.MESSAGE ) && 	    ( ( dc.failures==0 ) || ( dc.failcount < dc.failures ) ) ) )  	 || ( !(dc.message % dc.big_write_period ) ) ) {      if ( c.printtot ) { /* write */	fp = fopen ( c.totoutfile , "w" ) ;	if( !fp ) {	  fprintf( stderr, "No such file: %s\n", c.totoutfile ) ;	  finalline ( stderr , &all , &dc ,  1 ) ;	}      else 	{	  finalline ( fp , &all , &dc , 1 ) ;	  fclose ( fp ) ;	}      }      if ( c.printhisto && dc.block_valid ) { /* update histogram file */	fp = fopen ( c.histofile , "w" ) ;	if( !fp ) {	  fprintf( stderr, "No such file: %s\n", c.histofile ) ;	}      else 	{	  write_histo ( fp , &all ) ;	  fclose ( fp ) ;	}      }    }  }  if ( dc.use_fe >= 1 ) {    fe_free2 ( &p , c.mega_needed ) ;   }  if ( dc.use_bnd >= 1 ) {    bnd_free ( &bndp , &bndc ) ;  }  mnc_free ( &all ) ; }static void write_histo ( FILE *fp , mnc_all *c ) {  int l;  double t , cum = 0.0 ;  double  tot = (double) c->dc->block_valid ;   int *histo = c->vec->histo ;   int loops = c->bndc->loops ;  fprintf ( fp , "# total valid blocks %d\n" , (int) tot ) ;   for ( l = 1 ; l <= loops ; l ++ ) {    t= (double) histo[l] ;    cum += t ;    fprintf ( fp , "%d\t%d\t%d\t%9.4g\t%9.4g\n" , l , (int)(t) , (int)(cum) , 	      t/tot , cum/tot ) ;   }}static void snappyline ( data_creation_param  *c ) {  printf ( "%d:%du%dd%dl/%d\t" , c->failcount ,	   c->block_undet, c->block_det,	   c->block_detlw , c->message - 1 ) ; fflush ( stdout ) ; }static double bern ( int errors , int trials , double *p_point , double *p_upper		   ,  double *p_lower , double zscore ) {  double factor ; double lowerbound = 1e-6 ;  if ( errors == 0  ) {    *p_point = lowerbound ;    *p_lower = lowerbound ;    *p_upper = 1.0 - exp( - zscore / (double)(trials) ) ;     factor = 100.0 ;   } else {    *p_point =  (double)(errors)/(double)(trials) ;    /* factor corresponding z sd's in log p */    factor = exp ( zscore * sqrt( (double)(trials - errors ) / ((double)(errors * trials ) ) )) ;    *p_upper = *p_point * factor ;     *p_lower = *p_point / factor ;   }  return factor ;}static void finalline ( FILE *fp , mnc_all *all , data_creation_param *dc , int totals ){  mnc_vectors *v = all->vec ;  fe_min_control *c = all->c ;   double factor ;