bp_decoder.h

//Weilei Zeng, April 8 , 2020
//copied from bp.c

#include "dist.h"
#include <itpp/itbase.h>
#include <stdexcept> //for invalid argument


const double INF_BP=10000; ///< infinity belief

/**return weight of bvec
 * should use itpp::BERC::count_errors
 */
int summation(itpp::bvec u){
  int result = 0;
  for ( int i=0; i<u.size(); i++){
    if (u(i)) result ++;
  }
  return result;
}



/**\class BP_Decoder
 * belief propagation decoder
 */
class BP_Decoder{
 public:
  bool silent_mode = false;///<not in use
  bool debug = false;

  itpp::GF2mat H; ///<parity check matrix
  int nvar, ncheck, nedge;
  int exit_iteration=50;
  double alpha=1.0; ///< used for normalized decoder, alpha=1.25
  int decode_mode = 2;
  std::string decode_mode_str="min sum";

  //should be replaced by constant. not in use yet
  const int STANDARD=1;
  const int MIN_SUM=2;
  const int NORMALIZATION=3;

  bool is_initialized=false;
  bool is_H_valid(itpp::GF2mat H_temp);
  void init(itpp::GF2mat H_temp);

  void set_exit_iteration(int exit_iteration_temp);
  void set_decode_mode(int decode_mode_temp);
  void set_decode_mode_str(std::string decode_mode_str_temp);
  void set_silent_mode(bool silent_mode_temp);
  void set_debug_mode(bool debug_mode_temp);
  void print_info();

  
  int decode(itpp::bvec syndrome, const itpp::vec & LLRin, itpp::vec &LLRout);
  int bp_syndrome_llr(itpp::bvec syndrome,  const itpp::vec & LLRin, itpp::vec   & LLRout);
  int bp_schedule(itpp::bvec syndrome,  const itpp::vec & LLRin, itpp::vec   & LLRout);
  bool match_syndrome(itpp::vec LLR, itpp::bvec syndrome);
  int bp_flexible( itpp::bvec syndrome,  const itpp::vec & LLRin, itpp::vec   & LLRout);
  
  int schedule_mode=0;
  std::string schedule_mode_str="no schedule, parallel";
  itpp::mat schedule;
  void set_schedule_mode(int schedule_mode_temp);
  //  void set_schedule(int schedule_mode_temp);
};

void BP_Decoder::set_silent_mode(bool silent_mode_temp){
  silent_mode = silent_mode_temp;  
  return;
}
void BP_Decoder::set_debug_mode(bool debug_mode_temp){
  debug = debug_mode_temp;
  return;
}


bool BP_Decoder::is_H_valid(itpp::GF2mat H_temp){
  //check degree >=2
  for ( int i =0; i<H_temp.rows();i++){
    if ( summation(H_temp.get_row(i)) < 2 ) return false;      
  }
  //variable degree >=1
  for ( int j = 0; j< H_temp.cols(); j++){
    if ( summation(H_temp.get_col(j)) < 1 ) return false;
  }
  return true;
}

/**set up parity check matrix*/
void BP_Decoder::init(itpp::GF2mat H_temp){
  H = H_temp;
  nvar = H.cols();
  ncheck = H.rows();  
  if ( ! is_H_valid(H) ) throw std::invalid_argument( "BP_Decoder: invalid parity check matrix H" );
  is_initialized = true;
  nedge=0;
  for ( int i =0; i<ncheck; i++){
    for ( int j = 0; j<nvar;j++){
      if (H(i,j)) nedge++;
    }
  }
  return;
}

void BP_Decoder::print_info(){
  if (is_initialized){
    std::cout<<"--- BP_Decoder --- nvar = "<<nvar
	     <<", ncheck ="<<ncheck
	     <<", schedule_mode ("<<schedule_mode<<") "<<schedule_mode_str.c_str()
	     <<", decode mode ("<<decode_mode<<") "<<decode_mode_str.c_str()
	     <<", exit_iteration = "<<exit_iteration
	     <<std::endl;

  }else{
    std::cout<<"decoder is not initialized"<<std::endl;
  }
  return;

}

void BP_Decoder::set_exit_iteration(int exit_iteration_temp){
  exit_iteration = exit_iteration_temp;
  return;
}



/**set decode_mode
 * should replace those int and string by class const int, not implemented yet
 */
void BP_Decoder::set_decode_mode(int decode_mode_temp){
  decode_mode = decode_mode_temp;
  alpha=1.0;//reset first
  switch (decode_mode_temp){
  case 1:
    decode_mode_str = "standard";
    break;
  case 2:
    decode_mode_str = "min sum";
    break;
  case 3:// part of min sum, just change alpha
    decode_mode_str = "normalization";
    alpha = 1.25;
    break;
  default:
    throw std::invalid_argument( "BP_Decoder: illegal decode mode" );
  }
  std::cout<<"BP_Dcoder: set decode mode "<<decode_mode<<" - "<<decode_mode_str.c_str()<<std::endl;
  return;
}

   
void BP_Decoder::set_decode_mode_str(std::string decode_mode_str_temp){
  decode_mode_str = decode_mode_str_temp;
  if (decode_mode_str_temp =="standard")    decode_mode = 1;   
  else if (decode_mode_str_temp == "min sum") decode_mode = 2;  
  else if (decode_mode_str_temp == "normalization") decode_mode = 3;
  else throw std::invalid_argument( "BP_Decoder: illegal decode mode string" );
  
  if (debug)  std::cout<<"BP_Dcoder: set decode mode "<<decode_mode<<" - "<<decode_mode_str.c_str()<<std::endl;
  return;
}

void BP_Decoder::set_schedule_mode(int schedule_mode_temp){
  schedule_mode = schedule_mode_temp;
  int s;
  //schedule mode is used in decode();
  switch ( schedule_mode ){
  case 0:
    //no schedule, parrallel schedule
    //    set_decode_mode_str("standard");
    schedule_mode_str = "no schedule, parallel; using bp_syndrome_llr()";
    break;
  case 1: // default schedule mode 1, edge by edge
    //set_decode_mode_str("standard");
    schedule_mode_str  = "edge by edge, using bp_schedule()";
    break;
  case 2: //flexible
    {// edge by edge
      schedule_mode_str ="edge by edge, using bp_flexible()";
      schedule.set_size(2*nedge,3);
      schedule.zeros();
      s=0;
      for ( int i =0; i<ncheck; i++ ){
	for ( int j = 0; j<nvar; j ++){
	  if (H(i,j)){
	    // 0 for c->v, 1 for v->c
	    schedule.set(s,0,0);
	    schedule.set(s,1,i);
	    schedule.set(s,2,j);
	    schedule.set(s+1,0,1);
	    schedule.set(s+1,1,i);
	    schedule.set(s+1,2,j);
	    s +=2;
	  }
	}
      }
      break;
    }
  case 3: //flexible
    { // edge by edge, switch var and check
      schedule_mode_str ="edge by edge, switch var and check order. using bp_flexible()";
      schedule.set_size(2*nedge,3);
      schedule.zeros();
      s=0;
      for ( int j = 0; j<nvar; j ++){
	for ( int i =0; i<ncheck; i++ ){	
	  if (H(i,j)){
	    schedule.set(s,0,0);
	    schedule.set(s,1,i);
	    schedule.set(s,2,j);	
	    schedule.set(s+1,0,1);
	    schedule.set(s+1,1,i);
	    schedule.set(s+1,2,j);
	    s +=2;
	  }
	}
      }
      break;
    }
  case 4: //flexible
    { // random
      set_schedule_mode(2);
      schedule_mode = 4;
      schedule_mode_str ="random purterbation of of mode 2, using bp_flexible()";
      //random permutate
      int permutation=2*nedge;//number of swaps
      itpp::RNG_randomize();//get randome seed 
      itpp::ivec source=itpp::randi(permutation,0,nedge*2-1);
      itpp::ivec target=itpp::randi(permutation,0,nedge*2-1);
      for ( int i =0; i< permutation; i++){
	schedule.swap_rows(source(i),target(i));
      }
      break;
    }
  case 5: //flexible
    { // edge by edge, 
      schedule_mode_str ="edge by edge, follow the paper. using bp_flexible()";
      schedule.set_size(2*nedge,3);
      schedule.zeros();
      s=0;
      std::cout<<"start updating schedule"<<std::endl;
      for ( int j = 0; j<nvar; j ++){
	for ( int i =0; i<ncheck; i++ ){	
	  if (H(i,j)){
	    schedule.set(s,0,0);
	    schedule.set(s,1,i);
	    schedule.set(s,2,j);
	    s ++;
	  }	
	}
	for ( int i =0; i<ncheck; i++ ){	
	  if (H(i,j)){
	    schedule.set(s,0,1);
	    schedule.set(s,1,i);
	    schedule.set(s,2,j);
	    s ++;
	  }	
	}
      }
      std::cout<<"finish updating schedule"<<std::endl;
      break;
    }
    
  default:
    throw std::invalid_argument( "BP_Decoder: illegal schedule" );
  }
  return;
}


int BP_Decoder::decode( itpp::bvec syndrome,  const itpp::vec & LLRin, itpp::vec   & LLRout){
  // a wrapper to determine using which decoding function
  switch ( schedule_mode){
  case 0: // no schedule
    return bp_syndrome_llr(syndrome,  LLRin,  LLRout);
  case 1://same position for u and v, one by one
    return bp_schedule( syndrome, LLRin, LLRout);
  case 2://same position for u and v, one by one
  case 3://same position for v and u, one by one,
  case 4: //random permutation of case 2
  case 5: //follow the paper
    return bp_flexible( syndrome, LLRin, LLRout);
  default:
    throw std::invalid_argument( "BP_Decoder: illegal schedule" );
    //return bp_syndrome_llr(syndrome,  LLRin,  LLRout);
  }  
}

bool BP_Decoder::match_syndrome(itpp::vec LLR, itpp::bvec syndrome){
  itpp::bvec  error = LLR < 0;
  return itpp::GF2mat(H*error - syndrome).is_zero();
}

//int BP_Decoder::bp_syndrome_llr(const GF2mat H,  itpp::bvec syndrome,  itpp::vec & LLRin, itpp::vec   & LLRout, int exit_iteration, int decode_mode){
int BP_Decoder::bp_syndrome_llr( itpp::bvec syndrome,  const itpp::vec & LLRin, itpp::vec   & LLRout){  
  // input: syndrome vector s, loglikelihood ratio LLRin and LLRout
  // exit_iteration: max number of iteration
  // decode_mode 1: standard, 2: min sum
  //LLR(x) = log( p(x)/ (1-p(x)) )
  // initially we assume zero error, so LLR = LLR(0)=log ( (1-p)/p )>0
  // bits_out = LLRout < 0;
  //output: number of iteration, negative if not converge.

  if ( itpp::GF2mat(syndrome).is_zero() ){
    //return zero error vector, which is the default input for LLRout
    return 1;
  }
  
  // bool debug = false;// enable/disable printing
  //initialize
  //int nvar = H.cols(), ncheck = H.rows();
  if (debug) std::cout<<"nvar = "<<nvar <<", ncheck = "<<ncheck<<std::endl;
  itpp::mat llrs = itpp::zeros(ncheck, nvar), LLRs = itpp::zeros(ncheck, nvar);  
  //LLRout.set_size(nvar);should be the same size
  //  bool match_syndrome = false;// a flag to indicate whether the syndrome has been satisfied
  
  for ( int i = 0; i< ncheck ; i++){
    for ( int j=0; j<nvar; j++){
      if (H(i,j)) {
	llrs.set(i,j,LLRin(j));
	//llrs.set(i,j,log( (1-p)/p ));
      }
    }
  }
  //  if (debug) std::cout<<"finish initialize"<<std::endl;

  int update_count=0;
  double sum=0;
  double llr,LLR;
  std::string str="";
  int sign = 1;
  double prod=1.0;
  //  bool degree_one=true;
  while ( update_count < exit_iteration ){
    //check to variable update, LLR
    switch (decode_mode){
    case 1://standard
      for ( int i = 0; i< ncheck ; i++){
	for ( int j=0; j<nvar; j++){
	  if (H(i,j)) {
	    prod=1.0;
	    str = "prod list: i,j,k,prod,llr:";	  
	    for ( int k=0; k<nvar; k++){
	      if ( H(i,k) ){
		if ( k != j ) {
		  prod = prod * tanh( llrs(i,k)/2 );
		  str += std::to_string(i)+",";
		  str += std::to_string(j)+",";
		  str += std::to_string(k)+",";
		  str += std::to_string(prod)+",";
		  str += std::to_string(llrs(i,k))+",";
		  //if (debug) std::cout<<",prod = "<<prod<<" i="<<i <<std::endl;
		}
	      }
	    }
	    
	    LLR = atanh(prod)*2;
	    if ( syndrome(i) ){
	      LLR = -LLR;
	    }
	    LLRs.set(i,j,LLR);

	    if (debug) if ( std::abs(LLR) > 1000000.0) std::cout<<"LLRs: LLR = "<<LLR<<", prod = "<<prod<<"\n"<<str<<std::endl<<"H.get_row(i)="<<H.get_row(i)<<std::endl <<"llrs.get_row(i)="<<llrs.get_row(i)<<std::endl;
	  }
	}
      }
    case 2://min sum
      for ( int i = 0; i< ncheck ; i++){
	for ( int j=0; j<nvar; j++){
	  if (H(i,j)) {
	    prod=INF_BP;
	    sign = 1;
	    //	    str = "prod list: i,j,k,prod,llr:";	  
	    for ( int k=0; k<nvar; k++){
	      if ( H(i,k) ){
		if ( k != j ) {
		  if (llrs(i,k)>0){
		    llr = llrs(i,k);
		  }else{
		    llr = -llrs(i,k);
		    sign = -sign;
		  }
		  prod = (llr < prod) ? llr: prod;//min( prod, llr);
			      
		    //		  prod = prod * tanh( llrs(i,k)/2 );
		    /*str += std::to_string(i)+",";
		  str += std::to_string(j)+",";
		  str += std::to_string(k)+",";
		  str += std::to_string(prod)+",";
		  str += std::to_string(llrs(i,k))+",";*/
		  //if (debug) std::cout<<",prod = "<<prod<<" i="<<i <<std::endl;
		}
	      }
	    }
	    LLR = sign * prod;
	    //	    double LLR = atanh(prod)*2;
	    if ( syndrome(i) ){
	      LLR = -LLR;
	    }
	    LLRs.set(i,j,LLR);

	    if (debug) if ( std::abs(LLR) > INF_BP ) std::cout<<"LLRs: LLR = "<<LLR<<", prod = "<<prod<<"\n"<<str<<std::endl<<"H.get_row(i)="<<H.get_row(i)<<std::endl <<"llrs.get_row(i)="<<llrs.get_row(i)<<std::endl;
	  }
	}
      }
    }
    
    //    if (debug) std::cout<<"finish check to variable update"<<std::endl;    
    //variable to check update, llr
    


    for ( int i = 0; i< ncheck ; i++){
      for ( int j=0; j<nvar; j++){
	if (H(i,j)) {
	  sum=  LLRin(j);	
	
	  for ( int t=0; t<ncheck; t++){
	    if ( H(t,j) ){
	      if ( t != i ) {
		sum += LLRs(t,j);
		
	      }
	    }
	  }
	  llrs.set(i,j,sum);
	  //	  if ( std::abs(sum) > 1000)	  std::cout<<"llrs: sum = "<<sum<<"\n"<<LLRs.get_col(j)<<std::endl;
	}

      }
    }

    //    if (debug) std::cout<<"finish variable to checkupdate"<<std::endl;
  
    // get output LLRout and check result
    //    match_syndrome = true;
    for ( int j=0; j<nvar; j++){
        sum=LLRin(j);
	//if (debug) std::cout<<" sum = "<<sum<<std::endl;
	for ( int t=0; t<ncheck; t++){
	  //if (debug) std::cout<<"t = "<<t<<", sum = "<<sum<<std::endl;
	  if ( H(t,j) ){
	      sum += LLRs(t,j);
	  }
	}
	//if (debug) std::cout<<"LLRout = "<<LLRout<<std::endl;
	LLRout.set(j,sum);
    }
    if (debug) std::cout<<"update_count = "<<update_count<<", LLRout = "<<floor(LLRout)<<std::endl ;
    //if (debug) std::cout<<"update_count = "<<update_count<<std::endl;
    //if (debug) draw_toric_x_error(LLRout<0);
    update_count++;
    if ( match_syndrome( LLRout, syndrome) ){
      break;
    }        
  }
    
  if (debug) std::cout<<"LLRout = "<<LLRout<<std::endl;

  //  if (debug) std::cout<<"llrs = "<<llrs<<std::endl;
  //if (debug) std::cout<<"LLRs = "<<LLRs<<std::endl;

  //not converge, output negative value
  if (! match_syndrome( LLRout, syndrome) ){
    update_count = - update_count;
  }
    
  return update_count ;
}


int BP_Decoder::bp_schedule( itpp::bvec syndrome,  const itpp::vec & LLRin, itpp::vec   & LLRout){
 // input: syndrome vector s, loglikelihood ratio LLRin and LLRout
  // exit_iteration: max number of iteration
  // decode_mode 1: standard, 2: min sum
  //LLR(x) = log( p(x)/ (1-p(x)) )
  // initially we assume zero error, so LLR = LLR(0)=log ( (1-p)/p )>0
  // bits_out = LLRout < 0;
  //output: number of iteration, negative if not converge.

  if ( itpp::GF2mat(syndrome).is_zero() ){
    //return zero error vector, which is the default input for LLRout
    return 1;
  }
  
  // bool debug = false;// enable/disable printing
  //initialize
  //int nvar = H.cols(), ncheck = H.rows();
  if (debug) std::cout<<"nvar = "<<nvar <<", ncheck = "<<ncheck<<std::endl;
  itpp::mat llrs = itpp::zeros(ncheck, nvar), LLRs = itpp::zeros(ncheck, nvar);  
  //LLRout.set_size(nvar);should be the same size
  //  bool match_syndrome = false;// a flag to indicate whether the syndrome has been satisfied
  
  for ( int i = 0; i< ncheck ; i++){
    for ( int j=0; j<nvar; j++){
      if (H(i,j)) {
	llrs.set(i,j,LLRin(j));
	//llrs.set(i,j,log( (1-p)/p ));
      }
    }
  }
  //  if (debug) std::cout<<"finish initialize"<<std::endl;

  // *********************************
  int update_count=0;
  double sum=0;
  double LLR;//llr is not used yet
  std::string str="";
  //  int sign; //not used yet
  double prod=1.0;
  //  bool degree_one=true;
  while ( update_count < exit_iteration ){
    //check to variable update, LLR
    //use standard updating rule, no min sum
      for ( int i = 0; i< ncheck ; i++){
	for ( int j=0; j<nvar; j++){
	  if (H(i,j)) {
	    prod=1.0;
	    //	    str = "prod list: i,j,k,prod,llr:";	  
	    for ( int k=0; k<nvar; k++){
	      if ( H(i,k) ){
		if ( k != j ) {
		  prod = prod * tanh( llrs(i,k)/2 );
		  /*str += std::to_string(i)+",";
		  str += std::to_string(j)+",";
		  str += std::to_string(k)+",";
		  str += std::to_string(prod)+",";
		  str += std::to_string(llrs(i,k))+",";*/
		  //if (debug) std::cout<<",prod = "<<prod<<" i="<<i <<std::endl;
		}
	      }
	    }
	    
	    LLR = atanh(prod)*2;
	    if ( syndrome(i) ){
	      LLR = -LLR;
	    }
	    LLRs.set(i,j,LLR);

	    if (debug) if ( std::abs(LLR) > 1000000.0) std::cout<<"LLRs: LLR = "<<LLR<<", prod = "<<prod<<"\n"<<str<<std::endl<<"H.get_row(i)="<<H.get_row(i)<<std::endl <<"llrs.get_row(i)="<<llrs.get_row(i)<<std::endl;
	    //}
	  //}
	//}

    
      //    for ( int i = 0; i< ncheck ; i++){
      //for ( int j=0; j<nvar; j++){
      //if (H(i,j)) {
	  sum=  LLRin(j);		
	  for ( int t=0; t<ncheck; t++){
	    if ( H(t,j) ){
	      if ( t != i ) {
		sum += LLRs(t,j);
		
	      }
	    }
	  }
	  llrs.set(i,j,sum);
	  //	  if ( std::abs(sum) > 1000)	  std::cout<<"llrs: sum = "<<sum<<"\n"<<LLRs.get_col(j)<<std::endl;
	}

      }
    }

    //    if (debug) std::cout<<"finish variable to checkupdate"<<std::endl;
  
    // get output LLRout and check result
    //    match_syndrome = true;
    for ( int j=0; j<nvar; j++){
        sum=LLRin(j);
	//if (debug) std::cout<<" sum = "<<sum<<std::endl;
	for ( int t=0; t<ncheck; t++){
	  //if (debug) std::cout<<"t = "<<t<<", sum = "<<sum<<std::endl;
	  if ( H(t,j) ){
	      sum += LLRs(t,j);
	  }
	}
	//if (debug) std::cout<<"LLRout = "<<LLRout<<std::endl;
	LLRout.set(j,sum);
    }
    if (debug) std::cout<<"update_count = "<<update_count<<", LLRout = "<<floor(LLRout)<<std::endl ;
    //if (debug) std::cout<<"update_count = "<<update_count<<std::endl;
    //if (debug) draw_toric_x_error(LLRout<0);
    update_count++;
    if ( match_syndrome( LLRout, syndrome) ){
      break;
    }        
  }
    
  if (debug) std::cout<<"LLRout = "<<LLRout<<std::endl;

  //  if (debug) std::cout<<"llrs = "<<llrs<<std::endl;
  //if (debug) std::cout<<"LLRs = "<<LLRs<<std::endl;

  //not converge, output negative value
  if (! match_syndrome( LLRout, syndrome) ){
    update_count = - update_count;
  }
    
  return update_count ;

}


int BP_Decoder::bp_flexible( itpp::bvec syndrome,  const itpp::vec & LLRin, itpp::vec   & LLRout){
  // a flexible function to encorporate all variations
 // input: syndrome vector s, loglikelihood ratio LLRin and LLRout
  // exit_iteration: max number of iteration
  // decode_mode 1: standard, 2: min sum
  //LLR(x) = log( p(x)/ (1-p(x)) )
  // initially we assume zero error, so LLR = LLR(0)=log ( (1-p)/p )>0
  // bits_out = LLRout < 0;
  //output: number of iteration, negative if not converge.

  if ( itpp::GF2mat(syndrome).is_zero() )     return 1;
    //return zero error vector, which is the default input for LLRout
  
  //initialize
  if (debug) std::cout<<"nvar = "<<nvar <<", ncheck = "<<ncheck<<std::endl;
  itpp::mat llrs = itpp::zeros(ncheck, nvar), LLRs = itpp::zeros(ncheck, nvar);  
  //LLRout.set_size(nvar);should be the same size
  
  for ( int i = 0; i< ncheck ; i++){
    for ( int j=0; j<nvar; j++){
      if (H(i,j)) {
	llrs.set(i,j,LLRin(j));
	//llrs.set(i,j,log( (1-p)/p ));
      }
    }
  }

  // ********************************* start updating cycle
  int update_count=0;
  double sum=0;
  double LLR, llr;
  std::string str="";
  int sign; 
  double prod=1.0;
  while ( update_count < exit_iteration ){
    int i,j;
    for ( int is = 0; is < nedge*2; is ++){
      i = schedule(is,1);
      j = schedule(is,2);      
      //      for ( int i = 0; i< ncheck ; i++){
      //for ( int j=0; j<nvar; j++){
      int direction = schedule(is,0);
      switch ( direction ){
      case 0:     //check to variable update, LLR
	{
	  switch ( decode_mode ) {
	  case 1: //standard
	    {
	      prod=1.0;
	      for ( int k=0; k<nvar; k++){
		if ( H(i,k) ){
		  if ( k != j ) prod = prod * tanh( llrs(i,k)/2 );		  
		}
	      }
	      
	      LLR = atanh(prod)*2;
	      if ( syndrome(i) ) LLR = -LLR;	      
	      LLRs.set(i,j,LLR);
	    
	      if (debug) if ( std::abs(LLR) > 1000000.0) std::cout<<"LLRs: LLR = "<<LLR<<", prod = "<<prod<<"\n"<<str<<std::endl<<"H.get_row(i)="<<H.get_row(i)<<std::endl <<"llrs.get_row(i)="<<llrs.get_row(i)<<std::endl;
	      break;
	    }

	  case 2://min sum
	    //      for ( int i = 0; i< ncheck ; i++){
	    //for ( int j=0; j<nvar; j++){
	    //  if (H(i,j)) {
	    {
	      prod=INF_BP;
	      sign = 1;
	      //	    str = "prod list: i,j,k,prod,llr:";	
	      for ( int k=0; k<nvar; k++){
		if ( H(i,k) ){
		  if ( k != j ) {
		    if (llrs(i,k)>0){
		      llr = llrs(i,k);
		    }else{
		      llr = -llrs(i,k);
		      sign = -sign;
		    }
		    prod = (llr<prod) ? llr:prod;//min( prod, llr);
		  }
		}
	      }
	      LLR = sign * prod;
	      if ( syndrome(i) ) LLR = -LLR;
	      LLRs.set(i,j,LLR);

	      if (debug) if ( std::abs(LLR) > INF_BP ) std::cout<<"LLRs: LLR = "<<LLR<<", prod = "<<prod<<"\n"<<str<<std::endl<<"H.get_row(i)="<<H.get_row(i)<<std::endl <<"llrs.get_row(i)="<<llrs.get_row(i)<<std::endl;

	      break;
	    }
	    break;
	  }
	}
      case 1:
	{
	    // variable to check, llr updating
	  sum=  LLRin(j);		
	  for ( int t=0; t<ncheck; t++){
	    if ( H(t,j) ){
	      if ( t != i ) {
		sum += LLRs(t,j);		
	      }
	    }
	  }
	  llrs.set(i,j,sum);
	  break;
	}
      }
    }
  
    // get output LLRout and check result
    for ( int j=0; j<nvar; j++){
        sum=LLRin(j);
	//if (debug) std::cout<<" sum = "<<sum<<std::endl;
	for ( int t=0; t<ncheck; t++){
	  //if (debug) std::cout<<"t = "<<t<<", sum = "<<sum<<std::endl;
	  if ( H(t,j) ){
	      sum += LLRs(t,j);
	  }
	}
	//if (debug) cout<<"LLRout = "<<LLRout<<std::endl;
	LLRout.set(j,sum);
    }
    if (debug) std::cout<<"update_count = "<<update_count<<", LLRout = "<<floor(LLRout)<<std::endl ;
    //if (debug) std::cout<<"update_count = "<<update_count<<std::endl;
    //if (debug) draw_toric_x_error(LLRout<0);
    update_count++;
    if ( match_syndrome( LLRout, syndrome) ){
      break;
    }        
  }
    
  if (debug) std::cout<<"LLRout = "<<LLRout<<std::endl;

  //  if (debug) std::cout<<"llrs = "<<llrs<<std::endl;
  //if (debug) std::cout<<"LLRs = "<<LLRs<<std::std::endl;

  //not converge, output negative value
  if (! match_syndrome( LLRout, syndrome) ){
    update_count = - update_count;
  }
    
  return update_count ;

}