mnc3.c

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

  if ( dc->verbose > 1 ) printf ( "s.%6ld." , dc->mseed ) ; 
  fflush (stdout ) ; 
  if ( dc->verbose ) printf ( "Making code\n" ) ; 
  ran_seed ( dc->mseed ) ; 

  switch ( dc->MNC ) {
  case ( 2) :
    code->C1 = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    code->C1I = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
  
    code->C2 = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    code->C2I = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 

    sparse_invertible_cmatrix ( code->C1 , code->C1I , dc->pC1.per_row , dc->N ) ;
    sparse_invertible_cmatrix ( code->C2I , code->C2 , dc->pC2.per_row , dc->N ) ;
  
    status += cmatrices_2_alist ( code->C1 , code->C2I , &(p->a) , dc->N , 1 ) ;
    break ; 
  case(4):
    code->C1 = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    code->C1I = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
  
    code->C2 = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    code->C2I = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 

    sparse_invertible_cmatrix ( code->C1I , code->C1 , dc->pC1.per_row , dc->N ) ;
    sparse_invertible_cmatrix ( code->C2 , code->C2I , dc->pC2.per_row , dc->N ) ;
  
    status += cmatrices_2_alist ( code->C2 , code->C1I , &(p->a) , dc->N , 0 ) ;
    break; 
  case(1): 
  
    code->C1 = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    code->C1I = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    sparse_invertible_cmatrix ( code->C1 , code->C1I , dc->pC1.per_row , dc->N ) ;
    
    status += cmatrix_2_alist ( code->C1 , &(p->a) , dc->N ) ; 
    break ;
  case(3):
  
    code->C1 = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    code->C1I = cmatrix ( 1 , dc->N , 1 , dc->N ) ; 
    sparse_invertible_cmatrix ( code->C1 , code->C1I , dc->pC1.per_row , dc->N ) ;
    
    status += cmatrix_andI_2_alist ( code->C1 , &(p->a) , dc->N , 0 ) ;
    break ;
  case ( 0):
/* make MNC=0 style of A matrix -- see cg.c, but go directly to list 
 representation */
    status += code0_matrix_alist ( &(p->a) , dc->N , dc->M , dc->pC1.per_row ) ;
    break ; 
  default : 
    status --;
    break ; 
  }

/* more work below here */
  if ( dc->pbm_o ) {
    fp = fopen ( dc->pbm_ofile , "w" ) ;
    cmatrix2pbm ( code->C1 , 1 , dc->N , 1 , dc->N , fp ) ; 
    fclose ( fp ) ; 
    if ( dc->pbm_xv ) {
      sprintf ( junk , "xv %s -geometry +10+10 &" , dc->pbm_ofile ) ;
      system ( junk ) ;
    }
  }
  return status ; 
}

static int report_code ( data_creation_param *dc , mnc_code *code , fe_min_param *p ) {
  int status = 0 ; 

  if ( dc->verbose ) printf ( "The code\n" ) ; 

  if ( dc->verbose ) {
    switch ( dc->MNC ) {
    case ( 2) :
      printoutcmatrix ( code->C1 , 1 , dc->N , 1 , dc->N ) ; 
      printoutcmatrix ( code->C2I , 1 , dc->N , 1 , dc->N ) ; 
      break ; 
    case(4):
      printoutcmatrix ( code->C1I , 1 , dc->N , 1 , dc->N ) ; 
      printoutcmatrix ( code->C2 , 1 , dc->N , 1 , dc->N ) ; 
      break; 
    case(1): 
      printoutcmatrix ( code->C1 , 1 , dc->N , 1 , dc->N ) ; 
      break ;
    case(3):
      printoutcmatrix ( code->C1 , 1 , dc->N , 1 , dc->N ) ; 
      break ;
    case ( 0):
      /* make MNC=0 style of A matrix -- see cg.c, but go directly to list 
	 representation */
      printf ( "Code is only defined by alist\n" ) ; 
      break ; 
    default : 
      status --;
      break ; 
    }
    report_alist ( &(p->a) , dc->verbose + 1 ) ; 
  }

  return status ; 
}

static int check_alist ( alist_matrix *alist , fe_min_param *p , 
			fe_min_control *c , data_creation_param *dc )
{
  int status = 0 ; 
  int tmpi = 0 , n , m , MNC = dc->MNC ; 

  if ( dc->verbose ) printf ( "Checking alist\n" ) ; 

  if ( MNC == 0 || MNC == 1 || MNC == 2 || MNC == 4 ) {
    /* then every col should have more than 1 1 in it */
    for ( n = 1 ; n <= alist->N ; n ++ ) {
      if ( alist->num_nlist[n] <= 1 ) {
	tmpi ++ ; 
	fprintf ( stderr , "%d:" , n ) ; 
      }
    }
    if ( tmpi > 0 ) {
      fprintf ( stderr , "Warning unchecked bits n : %d\n", tmpi ) ;
      status -- ; 
    }
  }

  if ( MNC == 1 || MNC == 2 || MNC == 3 || MNC == 4 ) {
    /* then every row should have more than 1 1 in it */

    tmpi = 0 ; 
    for ( m = 1 ; m <= alist->M ; m ++ )  {
      if ( alist->num_mlist[m] <= 1 ) {
	tmpi ++ ; 
	fprintf ( stderr , "%d:" , m ) ; 
      }
    }
    if ( tmpi > 0 ) {
      fprintf ( stderr , "Warning unchecked bits n : %d\n", tmpi ) ;
      status -= 10 ; 
    }
  }
  if ( status < 0 && dc->verbose ) {
    fprintf ( stderr , "Failing ... here is the offending alist\n" ) ;
    report_alist ( alist , 2 ) ; 
  }
  else if ( dc->DEMO || dc->verbose >= 2 )   {
    report_alist ( alist , dc->verbose ) ; 
  }
  return status ; 
}

static int make_pvectors ( fe_min_param *p , 
			  fe_min_control *c , data_creation_param *dc , 
			  mnc_vectors  *v )
{
  int status = 0 ; 
  int m , n , nn ; 
  double lfn1fn = log ( dc->fn / ( 1 - dc->fn ) ) ;
  double lfs1fs = log ( dc->fs / ( 1 - dc->fs ) ) ;

  if ( dc->verbose ) printf ( "Making vectors for optimizer\n" ) ; 
  /* 
     log likelihood 
                     */
  switch ( dc->MNC ) {
  case(0): case(1): 
    for ( m = 1 ; m <= p->M ; m++ ) {
      p->g[m] = lfn1fn ;
      if ( v->r[m] ) p->g[m] *= -1.0 ; 
    }
    break ; 
  case(3):
    for ( m = 1 ; m <= p->M ; m++ ) {
      p->g[m] = -1.0 ; /* annealing will be needed here */
      if ( v->r[m] ) p->g[m] *= -1.0 ; 
      v->z[m] = v->r[m] ; 
    }
    break ; 
  case(2):
    for ( m = 1 ; m <= dc->M ; m++ ) {
      p->g[m] = lfn1fn ;
      if ( v->r[m] ) p->g[m] *= -1.0 ; 
    }
    for (  ; m <= p->M ; m++ ) {
      p->g[m] = lfs1fs ;
    }
    break ; 
  case(4):
    /* use the vector C_1^-1 r */
    for ( m = 1 ; m <= p->M ; m++ ) {
      p->g[m] = -1.0 ; /* annealing will be needed here */
      if ( v->cr[m] ) p->g[m] *= -1.0 ; 
      v->z[m] = v->cr[m] ; 
    }
    break ; 
  default:
    fprintf ( stderr , "MNC error \n" ) ; 
    break ; 
  } 
  /* 
     prior / bias
                     */
  switch ( dc->MNC ) {
  case(0): case(1):
    for ( n = 1 ; n <= p->N ; n ++ ) 
      p->bias[n] = lfs1fs ;
    break ; 
  case(2):
    for ( n = 1 ; n <= p->N ; n ++ ) 
      p->bias[n] = 0.0 ;
    break ; 
  case (3) : case ( 4) : 
    for ( n = 1 ; n <= dc->N ; n ++ ) 
      p->bias[n] = lfs1fs ;
    for (       ; n <= p->N ; n ++ ) 
      p->bias[n] = lfn1fn ;
    break ; 
  default:
    fprintf ( stderr , "MNC error \n" ) ; 
    status -- ; 
    break ; 
  } 
  /* 
     sneak the answer in 
                           */
  switch ( dc->MNC ) {
  case(0): case(1):
    p->true_s = 1 ; 
    for ( n = 1 ; n <= p->N ; n ++ ) 
      p->s[n] = v->s[n] ;
    break ; 
  case(2):
    p->true_s = 0 ; 
    break ; 
  case (3) : case ( 4) : 
    p->true_s = 1 ; 
    for ( n = 1 ; n <= dc->N ; n ++ ) 
      p->s[n] = v->s[n] ;
    for ( nn = 1 ; n <= p->N ; n ++ , nn ++ ) 
      p->s[n] = v->noise[nn] ; 
    break ; 
  default:
    fprintf ( stderr , "MNC error \n" ) ; 
    status -- ; 
    break ; 
  } 

  return status ; 
}

static int score ( fe_min_param *p , 
			  fe_min_control *c , data_creation_param *dc , 
       		  mnc_vectors  *v , mnc_code *code 
) /* compares so with s and works out quality of solution if wrong. */
{
  int status = 0 ; 
  int n ; 

  if ( dc->verbose ) printf ( "Scoring answer --- " ) ; 
  p->count = 0 ; 
  for ( n = 1 ; n <= dc->N ; n++ ) {
    if ( v->so[n] != v->s[n] ) {
      p->count ++ ; 
    }
  }
  if ( p->count > 0 ) {
    /* more work here */
    p->v = (double) (p->count) ; /* number of bits difference (used to be 
				    difference in free nergy */
    if ( dc->verbose ) printf ( "fail\n" ) ; 
    dc->failcount ++ ; 
  }  else {
    p->v = 0.0 ; 
    if ( dc->verbose ) printf ( "SUCCESS\n" ) ; 
  }
  /* could also work out the typicality of the solution */
  return status ; 
}


static int decode ( fe_min_param *p , 
			  fe_min_control *c , data_creation_param *dc , 
			  mnc_vectors  *v , mnc_code *code 
) /* 
     decode takes the x from the fe algorithm and puts the decoded 
     string in so 
                                                                   */
{
  int status = 0 ; 
  int n ; 

  if ( dc->verbose ) printf ( "Decoding output of optimizer\n" ) ; 

  switch ( dc->MNC ) {
  case(0):case(1):case(3):case(4):
    for ( n = 1 ; n <= dc->N ; n++ ) {
      v->so[n] = ( p->x[n] > 0.0 ) ? 1 : 0 ; 
    }
    break ; 
  case(2):
    for ( n = 1 ; n <= p->N ; n++ ) {
      v->uo[n] = ( p->x[n] > 0.0 ) ? 1 : 0 ; 
    }
    mult_cm_cv ( code->C2I , v->uo , v->so , 1 , dc->N ) ;
    break ; 
  default : 
    fprintf ( stderr , "MNC error \n" ) ; 
    status -- ; 
    break ; 
  } 
  return status ; 
}

static int make_norder ( alist_matrix *alist , fe_min_control *c )
{
  int status = 0 ; 
  int n , N=alist->N ; 

    /* set up the norder */
  switch ( c->NORDER ) {
  case(0):
    
    for ( n = 1 ; n <= N ; n++ ) {
      alist->norder[n] = n ; 
    }
    break;
  default:
    status -- ;
    fprintf ( stderr , "invalid NORDER %d \n" , c->NORDER ) ; 
    break ; 
  }

  return status ; 
}

static int make_vectors ( data_creation_param *dc , fe_min_param *p , 
		       mnc_code *code , mnc_vectors *v )
{
  int status = 0 ; 
  int m,n ; 
  int count_s , count_err , count_t=0 ; 

  if ( dc->message == 1 ) ran_seed ( dc->vseed ) ; 

  for ( n = 1, count_s = 0 ; n <= dc->N ; n++ ) {
    v->s[n] = ( ranu() < dc->fs ) ? 1 : 0 ; 
    count_s += v->s[n] ; 
  }
  dc->true_fs = (double) count_s / (double) dc->N ; 
  if( dc->DEMO || dc->verbose >= 1 ) {
    printf ( "generating string\n" ) ;
    printoutcvector ( v->s , 1 , dc->N ) ; 
  }

  /* Use A to generate t = A s */
  if ( dc->MNC==2 || dc->MNC==4 ) {
    mult_cm_cv ( code->C2 , v->s , v->u , 1 , dc->N ) ;
    mult_cm_cv ( code->C1 , v->u , v->t , 1 , dc->N ) ;
  } else if ( dc->MNC == 1 || dc->MNC == 3 ) {
    mult_cm_cv ( code->C1 , v->s , v->t , 1 , dc->N ) ;
  } else {
    alist_times_cvector_mod2 ( &(p->a) , v->s , v->t ) ; 
  }
  for ( m = 1 ; m <= dc->M ; m++ ) {
    count_t += v->t[m] ; 
  }
  dc->fte = (double) count_t / (double) dc->M ; 

  for ( count_err = 0 , m = 1 ; m <= dc->M ; m++ ) {
    count_err += ( v->noise[m] = ( ranu() < dc->fn ) ? 1 : 0 ) ; 
    v->r[m] = v->noise[m] ^ v->t[m] ;
  }
  dc->true_fn = (double) count_err / (double) dc->M ; 

  if ( dc->MNC==4 ) {
    mult_cm_cv ( code->C1I , v->r , v->cr , 1 , dc->N ) ;
  }

  if( dc->DEMO || dc->verbose >= 1 ) {
    printf ( "generating noise -> received:\n" ) ;
    printoutcvector ( v->noise , 1 , dc->N ) ; 
    printoutcvector ( v->r , 1 , dc->N ) ; 
  }

  return status ; 
}

int make_vectors_quick ( data_creation_param *dc , fe_min_param *p , 
			        mnc_vectors *v )
{
  int status = 0 ; 
  int m,n , u = 1 ; 

  switch ( dc->MNC ) {
  case(3): 
  case(4): 
/* Make the signal and noise vectors */  
/* Make the `x' vector from the signal and noise vectors */  
    for ( n = 1 ; n <= dc->N ; n++ , u ++ ) {
      v->x[u] = v->s[n] = ( ranu() < dc->fs ) ? 1 : 0 ; 
    }
    for ( m = 1 ; m <= dc->M ; m++ , u++) {
      v->x[u] = v->noise[m] = ( ranu() < dc->fn ) ? 1 : 0  ;
    }
    break;
  default :
    pause_for_return();exit(0);
    break;
  }

  alist_times_cvector_mod2 ( &(p->a) , v->x , v->z ) ; 

/* is there a y vector ? */
/* for the moment I assume not, this is just MNC = 3 or 4 */

  return status ; 
}

/* 
   This routine takes the vector in v->xpartial, which represents 
   a hit list produced by a neural net, say, of bits to flip.
   It attacks this vector with a and subtracts the result from 
   v->z. 
   It notices if the resulting z is blank.
   It adds xpartial to xreconstructed.
                                                                     */
int update_z ( data_creation_param *dc , fe_min_param *p , 
			        mnc_vectors *v )
{
  int status = 0 ; 
  return status ;