diff --git a/bp_test.c b/bp_test.c
new file mode 100644
index 0000000..25db53f
--- /dev/null
+++ b/bp_test.c
@@ -0,0 +1,334 @@
+//Weilei Apr 5, write a test for bp deocoding, using steane 5 qubit code
+
+
+// copied from my_bp_1.c
+#include <future>
+#include <string>
+#include <iostream>
+#include <itpp/itbase.h>
+#include <itpp/itcomm.h>
+#include <sstream>
+#include <stdio.h>
+#include "my_lib.h"
+#include <math.h>
+using namespace std;
+using namespace itpp;
+
+
+
+// Read the code from files and do BP decoding
+//input:source file for stabilzier matrix; error propability p ;
+
+int decode( GF2mat G, GF2mat H, double p, mat * data, int col_index, int row_index);
+
+
+int main(int argc, char **argv){
+  Parser parser;
+  parser.init(argc,argv);
+  //p.set_silentmode(true);
+  string filename_G, filename_H, filename_result;
+  //parser.get(filename_G,"filename_G");
+  //parser.get(filename_H,"filename_H");
+  //parser.get(filename_result,"filename_result");
+  
+  vector<future<int>> pool;
+  vector<future<int>>::size_type pool_size=15;
+  std::chrono::milliseconds span (100);
+  //change parameter p, code size
+  //char * filename_result=argv[3];//prefix for the file
+  double p;
+  //parser.get(p,"p");
+  int sizes[]= {13,11,9,7,5};
+  string stabilizer_folder="data/toric/stabilizer";
+  //string error_folder="data/toric/bp_decoding4";
+  mat data(50,5*5);   //return result in a mat, 5 columns for each size. format defines in header
+  data.zeros();
+  int col_index=-5;
+  for ( int size : sizes){
+    //    cout<<"size = "<<size<<endl;
+    col_index +=5;//5 col for each size
+    filename_G = stabilizer_folder + "/toric_S_x_size_" + to_string(size) + ".mm";
+    filename_H = stabilizer_folder + "/toric_S_z_size_" + to_string(size) + ".mm";
+    //filename_result = error_folder + "/toric_S_size_" + to_string(size) + ".mm_rate";
+    int row_index=-1;
+    for ( double ip=-1;ip > -3;ip-=0.1){
+      //cout<<"ip = "<<ip<<endl;
+      row_index ++;
+      //atof(argv[4]);
+      p=pow(10,ip);
+      //p=ip/100000.0;//previous use 1000 division. Now use 100,000 division cause the thershold for toric codes seems to be around 0.1%.
+      //char filename_result_p[255];
+      //sprintf( filename_result_p,"%s%.5f",filename_result.c_str(),p);//append p in the file name
+      //string filename_result_p_string(filename_result_p);
+      //cout<<filename_result_p_string<<endl;
+      GF2mat G =MM_to_GF2mat(filename_G);      
+      GF2mat H =MM_to_GF2mat(filename_H);
+
+      //manage thread within limit of pool size
+      while ( pool.size() >= pool_size ){
+	//wait until some of them finish
+	//break;
+	 for(vector <future<int>> :: iterator it = pool.begin(); it != pool.end(); ++it){
+	  //printf("%d", *it);
+	   if ( it->wait_for(span) == future_status::ready){
+	    cout<<"."<<endl;
+	    pool.erase(it);	    
+	    cout<<"remove thread. pool_size="<<pool.size()<<endl;
+	    it--;
+	    //break;
+	  }
+	  }
+      }
+    
+      // future<int> fut = async(launch::async, decode,G, H, p, filename_result_p);
+      //pool.push_back(move(fut));
+      pool.push_back(async(launch::async, decode,G, H, p, &data,col_index, row_index));
+      
+      cout<<"add new thread. pool_size="<<pool.size()
+	  <<", size = "<<size
+	  <<", p = "<<p
+	  <<", row_index = "<<row_index
+	  <<", col_index = "<<col_index
+	  <<endl;
+    }
+  }
+  //wait all process to finish
+  for(vector <future<int>> :: iterator it = pool.begin(); it != pool.end(); ++it){
+    it->get();
+  }
+  //print result
+  //cout<<"data = "<<data<<endl;
+  string header="#header\nsizes: ";
+  for (int size : sizes){
+    header += to_string(size) + ", ";
+  }
+  header += ("\n p, P_c converge rate, zero error vectors, weight one error, weight 2 error");
+  string filename_data="gnuplot/result/iteration0-cycle1000.gnudat";
+  mat2gnudata(data,filename_data,header);
+  cout<<"save data to "<<filename_data.c_str()<<endl;
+  return 0;
+}
+
+//*********************** bp decoding function by itpp
+
+int decode( GF2mat G, GF2mat H, double p,  mat * data, int col_index, int row_index){
+  //parameter setup
+  int cycles=1000;//10000;//number of cycles: fro toric code, 10000 give reletively clear result
+  int exit_at_iteration=50;//parameter for bp decoding set_exit_condition()
+  //  int bound= (int) -4096 * log (p/(1-p));// -300; see note.pdf on how to choose bound
+  int bound = 0;
+  int max_repetition = 0;//10 for best result. supposed to be zero in our set up, just for a check
+  
+  
+  Real_Timer timer;
+  LDPC_Code C=GF2mat_to_LDPC_Code(H);  //convert GF2mat saved in .mm file to LDPC_Code
+  //    LDPC_Code C=MM_to_LDPC_Code(filename_G);  //convert GF2mat saved in .mm file to LDPC_Code
+  //LDPC_Code C = get_test_LDPC_Code();
+  
+  //bp decoding set up
+  C.set_exit_conditions(exit_at_iteration,true,true);//high perperformance. This number of iteration would not affect small errors which converges and stop very fast. Only affect those doesn't converge soon or not converge at all.
+  //  C.set_exit_conditions(50,true,true);  // 50 iterations,check syndrome always
+  //C.set_llrcalc(LLR_calc_unit(12,300,7)); //(recommended settings with "exact" computation via high resolution lookup table)
+  C.set_llrcalc(LLR_calc_unit(12,0,7)); // fast mode; see "llr.h" //not good resolution for the result, produce more errors
+  //  cout << C << endl;
+  int N = C.get_nvar(); // number of bits per codeword; N=2 x n x n is the size of the toric code.
+  
+  bvec bitsin = zeros_b(N);//original zero vector
+  GF2mat E_input_converge(bitsin, false),E_input_nonconverge(bitsin, false),E_output_converge(bitsin,false),E_output_nonconverge(bitsin, false);//false for row vector;//These two GF2mat saves the input errors and their output after bp_decoding. They are clasified into convergeds cases and non converged cases.
+  RNG_randomize();
+  timer.tic(); 
+
+  //double pp=p;//=0.05;//test error probability //=pvals(j);   //  double pp=pow(10.0,-pvals(j));
+  //pp=0.01  ans~1;pp=0.03  ans=-2500
+  BSC bsc(p); // initialize BSC, error channel with probability pp
+  int ans=0; //the value return by bp_decoding(), which is the number of iterations and with a minus sign if it doesn't converge
+      
+
+
+  int counts_converge=0;
+  //int counts_weight_zero=0;
+  //int counts_weight_one=0;
+  //int counts_weight_two=0;
+  for (int i=0;i<cycles;i++){
+    //bvec rec_bits = bitsin;   // input vector with manually-input errors, for test
+    //    C.set_exit_conditions(2500);
+    bvec rec_bits0 = bsc(bitsin);    // input vector with errors from bsc chanel
+    //bvec rec_bits = find_error(rec_bits0,H);//switch to en error with same syndrome
+    bvec rec_bits=rec_bits0;
+    
+    //    bvec rec_bits( error_transform(rec_bits0,MM_to_GF2mat(filename_G)) ); //for test
+
+    //manual input error for test
+    //    rec_bits = bitsin; //zero vector
+    //horizontally aligned double error 
+    //    rec_bits(0)=1;    rec_bits(1)=1; //not converge, loop between edges in a vertex
+    //perpendicular double error
+    //    rec_bits(1+5)=1;    rec_bits(25+1)=1; //not converge, the error does not change, r vibirate around -600 or -300
+	
+    //rec_bits(1)=1;    rec_bits(25+2)=1; //same as above
+    //vertically aligned double error
+    //rec_bits(25+1)=1;    rec_bits(25+5+1)=1; // not converge but give the right result when seting delta K to be 30 or -30, need ieration max = 2500.
+    //check logical error. weight 5 causes decoding failure, but weight 4 can be reduced to smaller weight, and can then be decoded.
+    //rec_bits(0)=1;
+    //rec_bits(5)=1;
+    //rec_bits(10)=1;
+    //rec_bits(16)=1;
+      //rec_bits(20)=1;
+      //rec_bits(21)=1;
+      //rec_bits(16+25)=1;
+      //rec_bits(21+25)=1;
+
+    //check vertex operator. good result.
+    //    rec_bits(5)=1;
+    //    rec_bits(5+1)=1;
+    //    rec_bits(25+1)=1;
+    //    rec_bits(25+5+1)=1;
+	
+
+    //check for 743 code
+    //rec_bits(0)=1;
+    //rec_bits(1)=1;
+    //rec_bits(2)=1;
+    //rec_bits(4)=1;
+    //rec_bits(5)=1;
+
+    //    rec_bits(38-1)=1;
+    //rec_bits(43-1)=1;
+    
+    //cout<<"rec_bits="<<rec_bits<<endl;
+    vec s(N), s0(N); // input LLR version 
+    QLLRvec llr_output;
+    for(int i1=0;i1<N;i1++)
+      s(i1)=s0(i1)=(rec_bits(i1)==1)?-log(1.0/p-1.0):log(1.0/p-1.0);//convert to input format for bp_decode  
+    //cout<<"pp="<<pp<<"\t log(1.0/pp-1.0)="<<log(1.0/pp-1.0)<<endl;
+    //cout<<"s="<<s<<endl;
+  
+    QLLRvec llr_input=C.get_llrcalc().to_qllr(s);
+    //cout<<"llr_input="<<llr_input<<endl;
+    ans=C.bp_decode(llr_input, llr_output);
+    //cout<<"iteration: ans="<<ans<<", ";
+    //    int bound= (int) 4096 * log (p/(1-p));// -300; see note.pdf on how to choose bound
+    //    bound=0;
+    //cout<<"bound value: "<<bound<<", below which means error"<<endl;
+    bvec bitsout = llr_output < bound; //30;//0;//output
+
+    //debug: check why bound increase the weight
+    /*if (true){
+      bvec bound0 = llr_output <0;
+      bvec bounds = llr_output <-6000;
+      int w = weight(bounds)-weight(bound0);
+      if ( w >0){
+	cout<<endl;
+	cout<<"------------------------------------------------------------------------------------Find w < 0 "<<endl;
+	cout<<endl;
+      }
+      }*/
+    
+    
+    //bvec bitsout2 = llr < -bound;// for a check
+    //it turns out that 0 is not a good choice, but 30 or -30 is a good choice. The idea is, when BP deesn't converge, it loop between some bonds in a cycle, with small K value close to zero. By choosing a critical value greater or smaller than those small value, the output bits will either include those bits or discard those bits.
+
+    //    int max_iteration = 10;
+    //int current
+
+    //int max_repetition = 0;//5;//100; usually converge to zero with 3
+    int current_iteration = ans;
+    for (int i_iteration = 0;i_iteration<max_repetition;i_iteration++){
+      if (ans >0){
+	break;
+      }
+      // usleep(100000);
+      //cout<<"iteration: "<<i_iteration<<", current_iteration = "<<current_iteration<<", ";
+      //      cout<<"bits_out, bound = "<<bound<<endl;
+      //      draw_toric_x_error(bitsout);
+      //cout<<"effective output bits, bound = "<<bound<<endl;
+      //draw_toric_x_error(rec_bits+bitsout+rec_bits0);
+
+      //check value of K, to compare value of bound
+      //cout<<"min abs in llr: "<<min(abs(llr_output))<<endl; //it is either large, or smaller than bound
+
+      //repeat decoding
+      //      C.get_llrcalc().to_qllr(s)
+      llr_input = C.get_llrcalc().to_qllr( C.get_llrcalc().to_double(llr_output));//not sure if needed, just want to resign llr_output to llr_input
+      ans = C.bp_decode(llr_input,llr_output);
+      bitsout = llr_output < bound;
+    }
+
+    bitsout=bitsout+rec_bits+rec_bits0;//cancel input vector
+    //bitsout=reduce_weight( bitsout, G);//remove trivial cycles
+    rec_bits=rec_bits0;//replace for compatibility for saving data
+
+    //count weight
+    /*switch ( weight(rec_bits) ){
+    case 0:
+      counts_weight_zero++;break;
+    case 1:
+      counts_weight_one++;break;
+    case 2:
+      counts_weight_two++;break;
+      }*/
+      
+    if(ans>=0){
+      if ( current_iteration < 0 ){
+	//cout<<endl;
+	//cout<<"-------------------------------make it converge after iteration."<<endl;
+	//cout<<endl;
+      }
+      /*
+      cout<<"llr="<<llr_output<<endl;
+      cout<<"rec_bits="<<rec_bits<<endl;
+      cout<<"bitsout ="<<bitsout<<endl;
+      //cout<<"diff   = "<<bitsin+bitsout<<endl;
+      cout<<"diff    ="<<rec_bits+bitsout<<endl;
+      */
+      //check min weight of bitsout?
+      //E_input_converge = append_vector(E_input_converge, rec_bits);
+      //E_output_converge = append_vector(E_output_converge, bitsout);//maybe all zero, but no harm to check            
+      counts_converge++;
+    }else{
+      /*
+       cout<<"not converge when i="<<i<<endl
+	  <<"rec_bits= "<<rec_bits<<endl
+	  <<"bitsout = "<<bitsout<<endl;
+      cout<<"diff    ="<<rec_bits+bitsout<<endl;
+      cout<<"llr = "<<llr_output<<endl;
+      cout<<"NC, \# of errors in rec_bits = "<<BERC::count_errors(bitsin,rec_bits);
+      cout<<", \# of errors in bitsout  = "<<BERC::count_errors(bitsin,bitsout)<<endl;
+      */
+
+      //E_input_nonconverge = append_vector(E_input_nonconverge,rec_bits);
+      //E_output_nonconverge = append_vector(E_output_nonconverge,bitsout);//this would be the input error for the random window decoder.
+      // break;//break here to just print the first non converged error
+    }
+    //    cout<<"rec_bits0"<<endl;
+    //draw_toric_x_error(rec_bits0);
+    /*
+cout<<"rec_bits"<<endl;
+    draw_toric_x_error(rec_bits);
+    cout<<"bits_out, bound = "<<bound<<endl;
+    draw_toric_x_error(bitsout);
+    //    cout<<"bits_out2, bound = "<< -bound<<endl;
+    //draw_toric_x_error(bitsout2);
+    cout<<"iteration  ans = "<<ans<<endl;
+    */
+  }
+  //cout<<"counts_converge="<<counts_converge<<endl;
+  double rate_converge=counts_converge*1.0/cycles;
+  cout<<"N = "<<N<<" =  2 x "<< N/2;
+  cout<<", p = "<<p;
+  cout<<", Converge rate ="<<rate_converge<<endl;
+  //  save_result(p,rate_converge,filename_result_p);//no need to save this. can get it by counting the size of the matrix when doing gnuplot
+
+  //save result to mat data, suing row_index and col_index
+  //data entries: p, rate
+  data->set(row_index,col_index,p);
+  data->set(row_index,col_index+1,rate_converge);
+  //data->set(row_index,col_index+2,1.0*counts_weight_zero/cycles);
+  //data->set(row_index,col_index+3,1.0*counts_weight_one/cycles);
+  //data->set(row_index,col_index+4,1.0*counts_weight_two/cycles);
+  //data->set(row_index,col_index+2,rate_converge);
+  
+  timer.toc_print();
+  return 0;
+}
+   
diff --git a/bp_test.c~ b/bp_test.c~
new file mode 100644
index 0000000..17f9df4
--- /dev/null
+++ b/bp_test.c~
@@ -0,0 +1,331 @@
+//Weilei  March 26, 2020. yse c++ async to manage multi threads
+// copied from bp_decoding4.c
+#include <future>
+#include <string>
+#include <iostream>
+#include <itpp/itbase.h>
+#include <itpp/itcomm.h>
+#include <sstream>
+#include <stdio.h>
+#include "my_lib.h"
+#include <math.h>
+using namespace std;
+using namespace itpp;
+
+
+
+// Read the code from files and do BP decoding
+//input:source file for stabilzier matrix; error propability p ;
+
+int decode( GF2mat G, GF2mat H, double p, mat * data, int col_index, int row_index);
+  
+int main(int argc, char **argv){
+  Parser parser;
+  parser.init(argc,argv);
+  //p.set_silentmode(true);
+  string filename_G, filename_H, filename_result;
+  //parser.get(filename_G,"filename_G");
+  //parser.get(filename_H,"filename_H");
+  //parser.get(filename_result,"filename_result");
+  
+  vector<future<int>> pool;
+  vector<future<int>>::size_type pool_size=15;
+  std::chrono::milliseconds span (100);
+  //change parameter p, code size
+  //char * filename_result=argv[3];//prefix for the file
+  double p;
+  //parser.get(p,"p");
+  int sizes[]= {13,11,9,7,5};
+  string stabilizer_folder="data/toric/stabilizer";
+  //string error_folder="data/toric/bp_decoding4";
+  mat data(50,5*5);   //return result in a mat, 5 columns for each size. format defines in header
+  data.zeros();
+  int col_index=-5;
+  for ( int size : sizes){
+    //    cout<<"size = "<<size<<endl;
+    col_index +=5;//5 col for each size
+    filename_G = stabilizer_folder + "/toric_S_x_size_" + to_string(size) + ".mm";
+    filename_H = stabilizer_folder + "/toric_S_z_size_" + to_string(size) + ".mm";
+    //filename_result = error_folder + "/toric_S_size_" + to_string(size) + ".mm_rate";
+    int row_index=-1;
+    for ( double ip=-1;ip > -3;ip-=0.1){
+      //cout<<"ip = "<<ip<<endl;
+      row_index ++;
+      //atof(argv[4]);
+      p=pow(10,ip);
+      //p=ip/100000.0;//previous use 1000 division. Now use 100,000 division cause the thershold for toric codes seems to be around 0.1%.
+      //char filename_result_p[255];
+      //sprintf( filename_result_p,"%s%.5f",filename_result.c_str(),p);//append p in the file name
+      //string filename_result_p_string(filename_result_p);
+      //cout<<filename_result_p_string<<endl;
+      GF2mat G =MM_to_GF2mat(filename_G);      
+      GF2mat H =MM_to_GF2mat(filename_H);
+
+      //manage thread within limit of pool size
+      while ( pool.size() >= pool_size ){
+	//wait until some of them finish
+	//break;
+	 for(vector <future<int>> :: iterator it = pool.begin(); it != pool.end(); ++it){
+	  //printf("%d", *it);
+	   if ( it->wait_for(span) == future_status::ready){
+	    cout<<"."<<endl;
+	    pool.erase(it);	    
+	    cout<<"remove thread. pool_size="<<pool.size()<<endl;
+	    it--;
+	    //break;
+	  }
+	  }
+      }
+    
+      // future<int> fut = async(launch::async, decode,G, H, p, filename_result_p);
+      //pool.push_back(move(fut));
+      pool.push_back(async(launch::async, decode,G, H, p, &data,col_index, row_index));
+      
+      cout<<"add new thread. pool_size="<<pool.size()
+	  <<", size = "<<size
+	  <<", p = "<<p
+	  <<", row_index = "<<row_index
+	  <<", col_index = "<<col_index
+	  <<endl;
+    }
+  }
+  //wait all process to finish
+  for(vector <future<int>> :: iterator it = pool.begin(); it != pool.end(); ++it){
+    it->get();
+  }
+  //print result
+  //cout<<"data = "<<data<<endl;
+  string header="#header\nsizes: ";
+  for (int size : sizes){
+    header += to_string(size) + ", ";
+  }
+  header += ("\n p, P_c converge rate, zero error vectors, weight one error, weight 2 error");
+  string filename_data="gnuplot/result/iteration0-cycle1000.gnudat";
+  mat2gnudata(data,filename_data,header);
+  cout<<"save data to "<<filename_data.c_str()<<endl;
+  return 0;
+}
+
+//*********************** bp decoding function by itpp
+
+int decode( GF2mat G, GF2mat H, double p,  mat * data, int col_index, int row_index){
+  //parameter setup
+  int cycles=1000;//10000;//number of cycles: fro toric code, 10000 give reletively clear result
+  int exit_at_iteration=50;//parameter for bp decoding set_exit_condition()
+  //  int bound= (int) -4096 * log (p/(1-p));// -300; see note.pdf on how to choose bound
+  int bound = 0;
+  int max_repetition = 0;//10 for best result. supposed to be zero in our set up, just for a check
+  
+  
+  Real_Timer timer;
+  LDPC_Code C=GF2mat_to_LDPC_Code(H);  //convert GF2mat saved in .mm file to LDPC_Code
+  //    LDPC_Code C=MM_to_LDPC_Code(filename_G);  //convert GF2mat saved in .mm file to LDPC_Code
+  //LDPC_Code C = get_test_LDPC_Code();
+  
+  //bp decoding set up
+  C.set_exit_conditions(exit_at_iteration,true,true);//high perperformance. This number of iteration would not affect small errors which converges and stop very fast. Only affect those doesn't converge soon or not converge at all.
+  //  C.set_exit_conditions(50,true,true);  // 50 iterations,check syndrome always
+  //C.set_llrcalc(LLR_calc_unit(12,300,7)); //(recommended settings with "exact" computation via high resolution lookup table)
+  C.set_llrcalc(LLR_calc_unit(12,0,7)); // fast mode; see "llr.h" //not good resolution for the result, produce more errors
+  //  cout << C << endl;
+  int N = C.get_nvar(); // number of bits per codeword; N=2 x n x n is the size of the toric code.
+  
+  bvec bitsin = zeros_b(N);//original zero vector
+  GF2mat E_input_converge(bitsin, false),E_input_nonconverge(bitsin, false),E_output_converge(bitsin,false),E_output_nonconverge(bitsin, false);//false for row vector;//These two GF2mat saves the input errors and their output after bp_decoding. They are clasified into convergeds cases and non converged cases.
+  RNG_randomize();
+  timer.tic(); 
+
+  //double pp=p;//=0.05;//test error probability //=pvals(j);   //  double pp=pow(10.0,-pvals(j));
+  //pp=0.01  ans~1;pp=0.03  ans=-2500
+  BSC bsc(p); // initialize BSC, error channel with probability pp
+  int ans=0; //the value return by bp_decoding(), which is the number of iterations and with a minus sign if it doesn't converge
+      
+
+
+  int counts_converge=0;
+  //int counts_weight_zero=0;
+  //int counts_weight_one=0;
+  //int counts_weight_two=0;
+  for (int i=0;i<cycles;i++){
+    //bvec rec_bits = bitsin;   // input vector with manually-input errors, for test
+    //    C.set_exit_conditions(2500);
+    bvec rec_bits0 = bsc(bitsin);    // input vector with errors from bsc chanel
+    //bvec rec_bits = find_error(rec_bits0,H);//switch to en error with same syndrome
+    bvec rec_bits=rec_bits0;
+    
+    //    bvec rec_bits( error_transform(rec_bits0,MM_to_GF2mat(filename_G)) ); //for test
+
+    //manual input error for test
+    //    rec_bits = bitsin; //zero vector
+    //horizontally aligned double error 
+    //    rec_bits(0)=1;    rec_bits(1)=1; //not converge, loop between edges in a vertex
+    //perpendicular double error
+    //    rec_bits(1+5)=1;    rec_bits(25+1)=1; //not converge, the error does not change, r vibirate around -600 or -300
+	
+    //rec_bits(1)=1;    rec_bits(25+2)=1; //same as above
+    //vertically aligned double error
+    //rec_bits(25+1)=1;    rec_bits(25+5+1)=1; // not converge but give the right result when seting delta K to be 30 or -30, need ieration max = 2500.
+    //check logical error. weight 5 causes decoding failure, but weight 4 can be reduced to smaller weight, and can then be decoded.
+    //rec_bits(0)=1;
+    //rec_bits(5)=1;
+    //rec_bits(10)=1;
+    //rec_bits(16)=1;
+      //rec_bits(20)=1;
+      //rec_bits(21)=1;
+      //rec_bits(16+25)=1;
+      //rec_bits(21+25)=1;
+
+    //check vertex operator. good result.
+    //    rec_bits(5)=1;
+    //    rec_bits(5+1)=1;
+    //    rec_bits(25+1)=1;
+    //    rec_bits(25+5+1)=1;
+	
+
+    //check for 743 code
+    //rec_bits(0)=1;
+    //rec_bits(1)=1;
+    //rec_bits(2)=1;
+    //rec_bits(4)=1;
+    //rec_bits(5)=1;
+
+    //    rec_bits(38-1)=1;
+    //rec_bits(43-1)=1;
+    
+    //cout<<"rec_bits="<<rec_bits<<endl;
+    vec s(N), s0(N); // input LLR version 
+    QLLRvec llr_output;
+    for(int i1=0;i1<N;i1++)
+      s(i1)=s0(i1)=(rec_bits(i1)==1)?-log(1.0/p-1.0):log(1.0/p-1.0);//convert to input format for bp_decode  
+    //cout<<"pp="<<pp<<"\t log(1.0/pp-1.0)="<<log(1.0/pp-1.0)<<endl;
+    //cout<<"s="<<s<<endl;
+  
+    QLLRvec llr_input=C.get_llrcalc().to_qllr(s);
+    //cout<<"llr_input="<<llr_input<<endl;
+    ans=C.bp_decode(llr_input, llr_output);
+    //cout<<"iteration: ans="<<ans<<", ";
+    //    int bound= (int) 4096 * log (p/(1-p));// -300; see note.pdf on how to choose bound
+    //    bound=0;
+    //cout<<"bound value: "<<bound<<", below which means error"<<endl;
+    bvec bitsout = llr_output < bound; //30;//0;//output
+
+    //debug: check why bound increase the weight
+    /*if (true){
+      bvec bound0 = llr_output <0;
+      bvec bounds = llr_output <-6000;
+      int w = weight(bounds)-weight(bound0);
+      if ( w >0){
+	cout<<endl;
+	cout<<"------------------------------------------------------------------------------------Find w < 0 "<<endl;
+	cout<<endl;
+      }
+      }*/
+    
+    
+    //bvec bitsout2 = llr < -bound;// for a check
+    //it turns out that 0 is not a good choice, but 30 or -30 is a good choice. The idea is, when BP deesn't converge, it loop between some bonds in a cycle, with small K value close to zero. By choosing a critical value greater or smaller than those small value, the output bits will either include those bits or discard those bits.
+
+    //    int max_iteration = 10;
+    //int current
+
+    //int max_repetition = 0;//5;//100; usually converge to zero with 3
+    int current_iteration = ans;
+    for (int i_iteration = 0;i_iteration<max_repetition;i_iteration++){
+      if (ans >0){
+	break;
+      }
+      // usleep(100000);
+      //cout<<"iteration: "<<i_iteration<<", current_iteration = "<<current_iteration<<", ";
+      //      cout<<"bits_out, bound = "<<bound<<endl;
+      //      draw_toric_x_error(bitsout);
+      //cout<<"effective output bits, bound = "<<bound<<endl;
+      //draw_toric_x_error(rec_bits+bitsout+rec_bits0);
+
+      //check value of K, to compare value of bound
+      //cout<<"min abs in llr: "<<min(abs(llr_output))<<endl; //it is either large, or smaller than bound
+
+      //repeat decoding
+      //      C.get_llrcalc().to_qllr(s)
+      llr_input = C.get_llrcalc().to_qllr( C.get_llrcalc().to_double(llr_output));//not sure if needed, just want to resign llr_output to llr_input
+      ans = C.bp_decode(llr_input,llr_output);
+      bitsout = llr_output < bound;
+    }
+
+    bitsout=bitsout+rec_bits+rec_bits0;//cancel input vector
+    //bitsout=reduce_weight( bitsout, G);//remove trivial cycles
+    rec_bits=rec_bits0;//replace for compatibility for saving data
+
+    //count weight
+    /*switch ( weight(rec_bits) ){
+    case 0:
+      counts_weight_zero++;break;
+    case 1:
+      counts_weight_one++;break;
+    case 2:
+      counts_weight_two++;break;
+      }*/
+      
+    if(ans>=0){
+      if ( current_iteration < 0 ){
+	//cout<<endl;
+	//cout<<"-------------------------------make it converge after iteration."<<endl;
+	//cout<<endl;
+      }
+      /*
+      cout<<"llr="<<llr_output<<endl;
+      cout<<"rec_bits="<<rec_bits<<endl;
+      cout<<"bitsout ="<<bitsout<<endl;
+      //cout<<"diff   = "<<bitsin+bitsout<<endl;
+      cout<<"diff    ="<<rec_bits+bitsout<<endl;
+      */
+      //check min weight of bitsout?
+      //E_input_converge = append_vector(E_input_converge, rec_bits);
+      //E_output_converge = append_vector(E_output_converge, bitsout);//maybe all zero, but no harm to check            
+      counts_converge++;
+    }else{
+      /*
+       cout<<"not converge when i="<<i<<endl
+	  <<"rec_bits= "<<rec_bits<<endl
+	  <<"bitsout = "<<bitsout<<endl;
+      cout<<"diff    ="<<rec_bits+bitsout<<endl;
+      cout<<"llr = "<<llr_output<<endl;
+      cout<<"NC, \# of errors in rec_bits = "<<BERC::count_errors(bitsin,rec_bits);
+      cout<<", \# of errors in bitsout  = "<<BERC::count_errors(bitsin,bitsout)<<endl;
+      */
+
+      //E_input_nonconverge = append_vector(E_input_nonconverge,rec_bits);
+      //E_output_nonconverge = append_vector(E_output_nonconverge,bitsout);//this would be the input error for the random window decoder.
+      // break;//break here to just print the first non converged error
+    }
+    //    cout<<"rec_bits0"<<endl;
+    //draw_toric_x_error(rec_bits0);
+    /*
+cout<<"rec_bits"<<endl;
+    draw_toric_x_error(rec_bits);
+    cout<<"bits_out, bound = "<<bound<<endl;
+    draw_toric_x_error(bitsout);
+    //    cout<<"bits_out2, bound = "<< -bound<<endl;
+    //draw_toric_x_error(bitsout2);
+    cout<<"iteration  ans = "<<ans<<endl;
+    */
+  }
+  //cout<<"counts_converge="<<counts_converge<<endl;
+  double rate_converge=counts_converge*1.0/cycles;
+  cout<<"N = "<<N<<" =  2 x "<< N/2;
+  cout<<", p = "<<p;
+  cout<<", Converge rate ="<<rate_converge<<endl;
+  //  save_result(p,rate_converge,filename_result_p);//no need to save this. can get it by counting the size of the matrix when doing gnuplot
+
+  //save result to mat data, suing row_index and col_index
+  //data entries: p, rate
+  data->set(row_index,col_index,p);
+  data->set(row_index,col_index+1,rate_converge);
+  //data->set(row_index,col_index+2,1.0*counts_weight_zero/cycles);
+  //data->set(row_index,col_index+3,1.0*counts_weight_one/cycles);
+  //data->set(row_index,col_index+4,1.0*counts_weight_two/cycles);
+  //data->set(row_index,col_index+2,rate_converge);
+  
+  timer.toc_print();
+  return 0;
+}
+   
diff --git a/gnuplot/result/iteration0-cycle1000.gnudat b/gnuplot/result/iteration0-cycle1000.gnudat
index 6bb902a..8f1727f 100644
--- a/gnuplot/result/iteration0-cycle1000.gnudat
+++ b/gnuplot/result/iteration0-cycle1000.gnudat
@@ -1,26 +1,26 @@
 #header
 sizes: 13, 11, 9, 7, 5, 
  p, P_c converge rate, zero error vectors, weight one error, weight 2 error
-0.100000	0.006000	0.000000	0.000000	0.000000	0.100000	0.030000	0.000000	0.000000	0.000000	0.100000	0.101000	0.000000	0.000000	0.000000	0.100000	0.267000	0.000000	0.000000	0.003000	0.100000	0.535000	0.005000	0.041000	0.059000	
-0.079433	0.020000	0.000000	0.000000	0.000000	0.079433	0.078000	0.000000	0.000000	0.000000	0.079433	0.186000	0.000000	0.000000	0.000000	0.079433	0.420000	0.000000	0.002000	0.016000	0.079433	0.660000	0.019000	0.066000	0.140000	
-0.063096	0.077000	0.000000	0.000000	0.000000	0.063096	0.172000	0.000000	0.000000	0.000000	0.063096	0.321000	0.000000	0.001000	0.003000	0.063096	0.539000	0.005000	0.006000	0.035000	0.063096	0.733000	0.030000	0.126000	0.246000	
-0.050119	0.164000	0.000000	0.000000	0.000000	0.050119	0.277000	0.000000	0.000000	0.000000	0.050119	0.441000	0.001000	0.003000	0.011000	0.050119	0.610000	0.003000	0.026000	0.078000	0.050119	0.816000	0.069000	0.207000	0.265000	
-0.039811	0.303000	0.000000	0.000000	0.000000	0.039811	0.430000	0.000000	0.000000	0.003000	0.039811	0.581000	0.001000	0.011000	0.030000	0.039811	0.735000	0.013000	0.074000	0.161000	0.039811	0.860000	0.112000	0.283000	0.267000	
-0.031623	0.446000	0.001000	0.000000	0.002000	0.031623	0.550000	0.000000	0.002000	0.014000	0.031623	0.685000	0.004000	0.028000	0.068000	0.031623	0.805000	0.046000	0.130000	0.212000	0.031623	0.920000	0.209000	0.333000	0.269000	
-0.025119	0.584000	0.000000	0.000000	0.010000	0.025119	0.706000	0.002000	0.014000	0.043000	0.025119	0.772000	0.015000	0.071000	0.147000	0.025119	0.844000	0.087000	0.193000	0.242000	0.025119	0.945000	0.260000	0.364000	0.232000	
-0.019953	0.694000	0.001000	0.007000	0.020000	0.019953	0.782000	0.009000	0.036000	0.089000	0.019953	0.857000	0.036000	0.117000	0.236000	0.019953	0.914000	0.154000	0.274000	0.264000	0.019953	0.965000	0.357000	0.387000	0.178000	
-0.015849	0.805000	0.004000	0.026000	0.066000	0.015849	0.843000	0.023000	0.084000	0.142000	0.015849	0.888000	0.067000	0.202000	0.262000	0.015849	0.950000	0.225000	0.328000	0.259000	0.015849	0.969000	0.451000	0.370000	0.129000	
-0.012589	0.857000	0.010000	0.061000	0.150000	0.012589	0.897000	0.056000	0.144000	0.224000	0.012589	0.923000	0.128000	0.277000	0.254000	0.012589	0.956000	0.291000	0.336000	0.229000	0.012589	0.975000	0.540000	0.314000	0.114000	
-0.010000	0.906000	0.033000	0.111000	0.199000	0.010000	0.940000	0.091000	0.211000	0.260000	0.010000	0.938000	0.166000	0.316000	0.273000	0.010000	0.980000	0.381000	0.364000	0.177000	0.010000	0.987000	0.633000	0.289000	0.068000	
-0.007943	0.930000	0.068000	0.171000	0.250000	0.007943	0.961000	0.146000	0.301000	0.271000	0.007943	0.979000	0.292000	0.374000	0.200000	0.007943	0.979000	0.469000	0.370000	0.123000	0.007943	0.997000	0.649000	0.294000	0.049000	
-0.006310	0.960000	0.115000	0.255000	0.255000	0.006310	0.968000	0.194000	0.331000	0.255000	0.006310	0.977000	0.355000	0.377000	0.179000	0.006310	0.988000	0.526000	0.339000	0.106000	0.006310	0.991000	0.725000	0.226000	0.046000	
-0.005012	0.976000	0.170000	0.311000	0.286000	0.005012	0.977000	0.291000	0.372000	0.228000	0.005012	0.992000	0.473000	0.341000	0.138000	0.005012	0.995000	0.623000	0.301000	0.060000	0.005012	0.995000	0.766000	0.195000	0.033000	
-0.003981	0.980000	0.263000	0.360000	0.225000	0.003981	0.988000	0.350000	0.401000	0.190000	0.003981	0.991000	0.494000	0.363000	0.109000	0.003981	0.997000	0.690000	0.246000	0.054000	0.003981	0.997000	0.828000	0.156000	0.016000	
-0.003162	0.991000	0.329000	0.357000	0.219000	0.003162	0.993000	0.447000	0.382000	0.129000	0.003162	0.993000	0.610000	0.299000	0.075000	0.003162	1.000000	0.735000	0.229000	0.033000	0.003162	1.000000	0.853000	0.133000	0.014000	
-0.002512	0.995000	0.456000	0.331000	0.156000	0.002512	0.997000	0.531000	0.357000	0.090000	0.002512	0.995000	0.671000	0.250000	0.071000	0.002512	0.999000	0.785000	0.198000	0.015000	0.002512	0.998000	0.886000	0.103000	0.010000	
-0.001995	0.996000	0.520000	0.331000	0.119000	0.001995	0.999000	0.613000	0.301000	0.069000	0.001995	0.999000	0.737000	0.222000	0.037000	0.001995	1.000000	0.817000	0.167000	0.016000	0.001995	1.000000	0.911000	0.084000	0.005000	
-0.001585	0.995000	0.583000	0.309000	0.088000	0.001585	0.999000	0.677000	0.274000	0.047000	0.001585	0.999000	0.756000	0.216000	0.028000	0.001585	0.998000	0.869000	0.121000	0.009000	0.001585	1.000000	0.926000	0.071000	0.003000	
-0.001259	0.998000	0.637000	0.284000	0.071000	0.001259	1.000000	0.740000	0.229000	0.029000	0.001259	0.999000	0.817000	0.160000	0.022000	0.001259	1.000000	0.904000	0.090000	0.006000	0.001259	0.999000	0.925000	0.072000	0.003000	
+0.100000	0.006000	0.000000	0.000000	0.000000	0.100000	0.037000	0.000000	0.000000	0.000000	0.100000	0.102000	0.000000	0.000000	0.000000	0.100000	0.259000	0.000000	0.000000	0.000000	0.100000	0.548000	0.000000	0.000000	0.000000	
+0.079433	0.018000	0.000000	0.000000	0.000000	0.079433	0.077000	0.000000	0.000000	0.000000	0.079433	0.184000	0.000000	0.000000	0.000000	0.079433	0.377000	0.000000	0.000000	0.000000	0.079433	0.627000	0.000000	0.000000	0.000000	
+0.063096	0.065000	0.000000	0.000000	0.000000	0.063096	0.150000	0.000000	0.000000	0.000000	0.063096	0.329000	0.000000	0.000000	0.000000	0.063096	0.522000	0.000000	0.000000	0.000000	0.063096	0.731000	0.000000	0.000000	0.000000	
+0.050119	0.160000	0.000000	0.000000	0.000000	0.050119	0.265000	0.000000	0.000000	0.000000	0.050119	0.430000	0.000000	0.000000	0.000000	0.050119	0.626000	0.000000	0.000000	0.000000	0.050119	0.834000	0.000000	0.000000	0.000000	
+0.039811	0.266000	0.000000	0.000000	0.000000	0.039811	0.405000	0.000000	0.000000	0.000000	0.039811	0.587000	0.000000	0.000000	0.000000	0.039811	0.703000	0.000000	0.000000	0.000000	0.039811	0.859000	0.000000	0.000000	0.000000	
+0.031623	0.449000	0.000000	0.000000	0.000000	0.031623	0.553000	0.000000	0.000000	0.000000	0.031623	0.663000	0.000000	0.000000	0.000000	0.031623	0.810000	0.000000	0.000000	0.000000	0.031623	0.900000	0.000000	0.000000	0.000000	
+0.025119	0.554000	0.000000	0.000000	0.000000	0.025119	0.659000	0.000000	0.000000	0.000000	0.025119	0.791000	0.000000	0.000000	0.000000	0.025119	0.856000	0.000000	0.000000	0.000000	0.025119	0.936000	0.000000	0.000000	0.000000	
+0.019953	0.716000	0.000000	0.000000	0.000000	0.019953	0.795000	0.000000	0.000000	0.000000	0.019953	0.835000	0.000000	0.000000	0.000000	0.019953	0.911000	0.000000	0.000000	0.000000	0.019953	0.945000	0.000000	0.000000	0.000000	
+0.015849	0.767000	0.000000	0.000000	0.000000	0.015849	0.865000	0.000000	0.000000	0.000000	0.015849	0.915000	0.000000	0.000000	0.000000	0.015849	0.939000	0.000000	0.000000	0.000000	0.015849	0.976000	0.000000	0.000000	0.000000	
+0.012589	0.888000	0.000000	0.000000	0.000000	0.012589	0.921000	0.000000	0.000000	0.000000	0.012589	0.947000	0.000000	0.000000	0.000000	0.012589	0.950000	0.000000	0.000000	0.000000	0.012589	0.986000	0.000000	0.000000	0.000000	
+0.010000	0.906000	0.000000	0.000000	0.000000	0.010000	0.943000	0.000000	0.000000	0.000000	0.010000	0.953000	0.000000	0.000000	0.000000	0.010000	0.979000	0.000000	0.000000	0.000000	0.010000	0.990000	0.000000	0.000000	0.000000	
+0.007943	0.943000	0.000000	0.000000	0.000000	0.007943	0.963000	0.000000	0.000000	0.000000	0.007943	0.969000	0.000000	0.000000	0.000000	0.007943	0.980000	0.000000	0.000000	0.000000	0.007943	0.989000	0.000000	0.000000	0.000000	
+0.006310	0.967000	0.000000	0.000000	0.000000	0.006310	0.970000	0.000000	0.000000	0.000000	0.006310	0.987000	0.000000	0.000000	0.000000	0.006310	0.992000	0.000000	0.000000	0.000000	0.006310	0.996000	0.000000	0.000000	0.000000	
+0.005012	0.966000	0.000000	0.000000	0.000000	0.005012	0.983000	0.000000	0.000000	0.000000	0.005012	0.987000	0.000000	0.000000	0.000000	0.005012	0.993000	0.000000	0.000000	0.000000	0.005012	0.995000	0.000000	0.000000	0.000000	
+0.003981	0.985000	0.000000	0.000000	0.000000	0.003981	0.990000	0.000000	0.000000	0.000000	0.003981	0.990000	0.000000	0.000000	0.000000	0.003981	0.997000	0.000000	0.000000	0.000000	0.003981	0.998000	0.000000	0.000000	0.000000	
+0.003162	0.989000	0.000000	0.000000	0.000000	0.003162	0.990000	0.000000	0.000000	0.000000	0.003162	0.992000	0.000000	0.000000	0.000000	0.003162	0.999000	0.000000	0.000000	0.000000	0.003162	1.000000	0.000000	0.000000	0.000000	
+0.002512	0.993000	0.000000	0.000000	0.000000	0.002512	0.996000	0.000000	0.000000	0.000000	0.002512	0.999000	0.000000	0.000000	0.000000	0.002512	1.000000	0.000000	0.000000	0.000000	0.002512	0.999000	0.000000	0.000000	0.000000	
+0.001995	0.996000	0.000000	0.000000	0.000000	0.001995	0.996000	0.000000	0.000000	0.000000	0.001995	0.998000	0.000000	0.000000	0.000000	0.001995	0.999000	0.000000	0.000000	0.000000	0.001995	1.000000	0.000000	0.000000	0.000000	
+0.001585	1.000000	0.000000	0.000000	0.000000	0.001585	0.999000	0.000000	0.000000	0.000000	0.001585	0.999000	0.000000	0.000000	0.000000	0.001585	1.000000	0.000000	0.000000	0.000000	0.001585	1.000000	0.000000	0.000000	0.000000	
+0.001259	0.998000	0.000000	0.000000	0.000000	0.001259	0.999000	0.000000	0.000000	0.000000	0.001259	1.000000	0.000000	0.000000	0.000000	0.001259	0.998000	0.000000	0.000000	0.000000	0.001259	1.000000	0.000000	0.000000	0.000000	
 0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	
 0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	
 0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	
diff --git a/gnuplot/toric_rate.gnu b/gnuplot/toric_rate.gnu
index 2d1001b..1e93598 100755
--- a/gnuplot/toric_rate.gnu
+++ b/gnuplot/toric_rate.gnu
@@ -8,6 +8,7 @@ data_folder='result'
 data_title='iteration0-cycle1000'
 
 plot_title='-weight'
+plot_title=''
 
 data_file=data_folder.'/'.data_title.'.gnudat'
 
@@ -27,9 +28,9 @@ set logscale y
 set format y "10^{%L}"
 set logscale x
 set format x "10^{%L}"
-set xrange [0.001:0.1]
-set yrange [0.001:1]
-set title "converge rate of BP for toric code ".data_title
+set xrange [0.01:0.1]
+set yrange [0.01:1]
+set title "converge rate of BP for toric code ".data_title.plot_title
 
 set linetype cycle 5
 
@@ -43,12 +44,10 @@ size_set="13 11 9 7 5"
 color_set="violet blue red yellow black green"
 
 # plot rate
-#plot for [i=1:5] data_file using (column((i-1)*5+1)):(column((i-1)*5+2)) title 'size '.word(size_set,i)
+plot for [i=1:5] data_file using (column((i-1)*5+1)):(column((i-1)*5+2)) title 'size '.word(size_set,i) with linespoints ps 0.5 lc i
 
 # plot weight of input error
-plot for [i=1:5] data_file using (column((i-1)*5+1)):(column((i-1)*5+2)) title 'size '.word(size_set,i).' converge' ps 0.5,\
-for [i=1:5] data_file using (column((i-1)*5+1)):(column((i-1)*5+3) + column((i-1)*5+4)+ column((i-1)*5+5)) title 'size '.word(size_set,i).' weight 0+1+2' with lines linecolor i,\
-for [i=1:5] data_file using (column((i-1)*5+1)):(column((i-1)*5+3)) title 'size '.word(size_set,i).' weight 0' linecolor i dashtype ".-" with lines,
+# plot for [i=1:5] data_file using (column((i-1)*5+1)):(column((i-1)*5+2)) title 'size '.word(size_set,i).' converge' ps 0.5,for [i=1:5] data_file using (column((i-1)*5+1)):(column((i-1)*5+3) + column((i-1)*5+4)+ column((i-1)*5+5)) title 'size '.word(size_set,i).' weight 0+1+2' with lines linecolor i,for [i=1:5] data_file using (column((i-1)*5+1)):(column((i-1)*5+3)) title 'size '.word(size_set,i).' weight 0' linecolor i dashtype ".-" with lines,
 
 
 #plot data_file using 1:2 title 'size 13' with linespoints,\
diff --git a/makefile b/makefile
index 590c7ff..be885b9 100644
--- a/makefile
+++ b/makefile
@@ -13,10 +13,10 @@ END=`pkg-config --libs itpp`
 
 #update:
 INC_DIR=~/working/weilei_lib
-files=$(INC_DIR)/mm_read.c $(INC_DIR)/mm_read.h $(INC_DIR)/mmio.c $(INC_DIR)/mmio.h $(INC_DIR)/mm_write.c $(INC_DIR)/mm_write.h $(INC_DIR)/lib.cpp $(INC_DIR)/lib.h $(INC_DIR)/dist.c $(INC_DIR)/dist.h $(INC_DIR)/concatenation_lib.c $(INC_DIR)/concatenation_lib.h $(INC_DIR)/my_lib.h makefile
-#command=$(CXX) -I../../weilei_lib $(START) -o $@ $< $(word 2,$^) $(word 4, $^) $(word 6, $^) $(word 8, $^) $(word 10, $^) $(word 12, $^) $(END)
+files=$(INC_DIR)/mm_read.c $(INC_DIR)/mm_read.h $(INC_DIR)/mmio.c $(INC_DIR)/mmio.h $(INC_DIR)/mm_write.c $(INC_DIR)/mm_write.h $(INC_DIR)/lib.cpp $(INC_DIR)/lib.h $(INC_DIR)/dist.c $(INC_DIR)/dist.h $(INC_DIR)/concatenation_lib.c $(INC_DIR)/concatenation_lib.h $(INC_DIR)/bp.c $(INC_DIR)/bp.h $(INC_DIR)/my_lib.h makefile
+#command=$(CXX) -I../../weilei_lib $(START) -o $@ $< $(word 2,$^) $(word 4, $^) $(word 6, $^) $(word 8, $^) $(word 10, $^) $(word 12, $^) $(word 14, $^) $(END)
 
-command=$(CXX) -std=c++11 -pthread -I $(INC_DIR) $(START) -o $@ $< $(word 2,$^) $(word 4, $^) $(word 6, $^) $(word 8, $^) $(word 10, $^) $(word 12, $^) $(END) ; echo finish making
+command=$(CXX) -std=c++11 -pthread -I $(INC_DIR) $(START) -o $@ $< $(word 2,$^) $(word 4, $^) $(word 6, $^) $(word 8, $^) $(word 10, $^) $(word 12, $^) $(word 14, $^) $(END) ; echo finish making
 
 
 
@@ -39,10 +39,10 @@ data_collect_perm_check.out:data_collect_perm_check.c $(files)
 	$(command)
 convert_data.out:convert_data.c $(files)
 	$(command)
-test.out:test.c
-	g++ -std=c++11 -pthread -o test.out test.c
-#test.out:test.c $(files)
-#	$(command)
+#test.out:test.c
+#	g++ -std=c++11 -pthread -o test.out test.c
+test.out:test.c $(files)
+	$(command)
 	./test.out
 #	$(CXX) $(START) -o $@ $< $(word 2,$^) $(word 4, $^) $(word 6, $^) $(END)
 error_analysis.out:error_analysis.c mm_read.c mm_read.h mmio.c mmio.h mm_write.c mm_write.h my_lib.h
diff --git a/my_bp1.c b/my_bp1.c
index 37ba9ae..17f9df4 100644
--- a/my_bp1.c
+++ b/my_bp1.c
@@ -143,9 +143,9 @@ int decode( GF2mat G, GF2mat H, double p,  mat * data, int col_index, int row_in
 
 
   int counts_converge=0;
-  int counts_weight_zero=0;
-  int counts_weight_one=0;
-  int counts_weight_two=0;
+  //int counts_weight_zero=0;
+  //int counts_weight_one=0;
+  //int counts_weight_two=0;
   for (int i=0;i<cycles;i++){
     //bvec rec_bits = bitsin;   // input vector with manually-input errors, for test
     //    C.set_exit_conditions(2500);
@@ -252,18 +252,18 @@ int decode( GF2mat G, GF2mat H, double p,  mat * data, int col_index, int row_in
     }
 
     bitsout=bitsout+rec_bits+rec_bits0;//cancel input vector
-    bitsout=reduce_weight( bitsout, G);//remove trivial cycles
+    //bitsout=reduce_weight( bitsout, G);//remove trivial cycles
     rec_bits=rec_bits0;//replace for compatibility for saving data
 
     //count weight
-    switch ( weight(rec_bits) ){
+    /*switch ( weight(rec_bits) ){
     case 0:
       counts_weight_zero++;break;
     case 1:
       counts_weight_one++;break;
     case 2:
       counts_weight_two++;break;
-    }
+      }*/
       
     if(ans>=0){
       if ( current_iteration < 0 ){
@@ -320,9 +320,9 @@ cout<<"rec_bits"<<endl;
   //data entries: p, rate
   data->set(row_index,col_index,p);
   data->set(row_index,col_index+1,rate_converge);
-  data->set(row_index,col_index+2,1.0*counts_weight_zero/cycles);
-  data->set(row_index,col_index+3,1.0*counts_weight_one/cycles);
-  data->set(row_index,col_index+4,1.0*counts_weight_two/cycles);
+  //data->set(row_index,col_index+2,1.0*counts_weight_zero/cycles);
+  //data->set(row_index,col_index+3,1.0*counts_weight_one/cycles);
+  //data->set(row_index,col_index+4,1.0*counts_weight_two/cycles);
   //data->set(row_index,col_index+2,rate_converge);
   
   timer.toc_print();
diff --git a/my_bp1.out b/my_bp1.out
index 82ae267..eec3906 100755
Binary files a/my_bp1.out and b/my_bp1.out differ
diff --git a/temp/async_test.c b/old/async_test.c
similarity index 100%
rename from temp/async_test.c
rename to old/async_test.c
diff --git a/temp/async_test.out b/old/async_test.out
similarity index 100%
rename from temp/async_test.out
rename to old/async_test.out
diff --git a/temp/makefile b/old/makefile
similarity index 100%
rename from temp/makefile
rename to old/makefile
diff --git a/temp/prime_calculator.c b/old/prime_calculator.c
similarity index 100%
rename from temp/prime_calculator.c
rename to old/prime_calculator.c
diff --git a/pattern_observer.out b/pattern_observer.out
index 141037e..0bfa72e 100755
Binary files a/pattern_observer.out and b/pattern_observer.out differ
diff --git a/readme.md b/readme.md
index e1bdaf7..dd87beb 100644
--- a/readme.md
+++ b/readme.md
@@ -10,5 +10,9 @@
 - Apr 4
   - [x] plot result.
   - [x] plotted counts for errors with weight 0, 1, 2 respectively. Those counts are smaller than the converged error. Because in large-weight errors, there will be single error as well. Hence these counts are not meaningful. For a meaningful comparison, one should plot the weight of the errors, instead of counts. This comparison has been done before, hence not shown here. For reference see "Weilei Research Note.pdf"
-  - [ ] replace bp decoding
+  - [x] check iteration, see improvement about 20-30% in converge rate. no threshold.
+- APr 5
+  - [x] replace bp decoding using syndrome based and LLR simplied. Passing test with 7 qubit code
+  - [ ] test repetition code
+  - [ ] test toric code and check result
 
diff --git a/test.c b/test.c
index f8fa815..6cb38f1 100644
--- a/test.c
+++ b/test.c
@@ -1,9 +1,17 @@
 
 // CPP program to demonstrate multithreading
 // using three different callables.
+#include <future>
 #include <iostream>
 #include <thread>
+#include <itpp/itbase.h>
+#include <itpp/itcomm.h>
+#include "my_lib.h"
+#include <math.h>
+
 using namespace std;
+using namespace itpp;
+
 
 // A dummy function
 void foo(int Z)
@@ -25,7 +33,9 @@ class thread_obj {
   }
 };
 
-int main()
+//cannot compile unless it is in main()
+/*
+int thread_test()
 {
       cout << "Threads 1 and 2 and 3 "
 	"operating independently" << endl;
@@ -64,3 +74,48 @@ int main()
 
       return 0;
 } 
+*/
+
+int bp_test(){
+  cout<<"Test the bp_syndrome_llr function"<<endl;
+  GF2mat H = get_check(3,7);
+  //  GF2mat H = get_check(2,7);
+  H = H.get_submatrix(0,0,2,6);
+  cout<<"parity check H = \n"<<H<<endl;
+  bvec error = zeros_b(H.cols());
+  error.set(0,1);
+  error.set(1,1);
+  error.set(2,1);
+  cout<<"input error ="<<error<<endl;
+  bvec syndrome = (H*error);
+  cout<<"syndrome = "<<syndrome<<endl;
+  vec LLRin(H.cols());
+  double p =0.1;
+  double LLR=log( (1-p)/p);
+  LLRin.ones();
+  LLRin = LLRin*LLR;    
+  vec LLRout(LLRin);//default, zero syndrome give zero error
+  //LLRout.zeros();
+  cout<<"LLRin  = "<< LLRin<<endl;
+  //  cout<<"LLRout = "<< LLRout<<endl;
+  int exit_iteration = 5;
+  int iteration = bp_syndrome_llr(H,syndrome,LLRin, LLRout, exit_iteration);
+  cout<<" iteration = "<<iteration <<endl;
+  cout<<" LLRout = "<<LLRout<<endl;
+  bvec bits_out = LLRout < 0;
+  cout<<"bits_out = "<<bits_out<<endl;
+  cout<<"error    = "<<error<<endl;
+  bvec error_residue = error + bits_out;
+  bvec syndrome_residue = H*error_residue;
+  cout<<"error_residue = "<< error_residue <<endl;
+  cout<<"syndrome_residue = "<< syndrome_residue <<endl;
+  return 0;
+}
+
+
+int main(){
+  //thread_test();
+  bp_test();
+  
+
+}
diff --git a/test.out b/test.out
index c052e66..14e2285 100755
Binary files a/test.out and b/test.out differ