56 const std::string& format): init_flag(false)
58 if (format ==
"alist") {
62 it_error(
"LDPC_Parity::LDPC_Parity(): Only 'alist' format is supported");
87 "LDPC_Parity::set(): Wrong index i");
89 "LDPC_Parity::set(): Wrong index j");
92 int diff =
static_cast<int>(x) - static_cast<int>(
H(i, j));
112 "LDPC_Parity::display_stats(): Object not initialized");
115 vec vdeg =
zeros(cmax + 1);
116 vec cdeg =
zeros(vmax + 1);
117 for (
int col = 0; col <
nvar; col++) {
120 "LDPC_Parity::display_stats(): Internal error");
122 for (
int row = 0; row <
ncheck; row++) {
125 "LDPC_Parity::display_stats(): Internal error");
140 it_info(
"--- LDPC parity check matrix ---");
141 it_info(
"Dimension [ncheck x nvar]: " << ncheck <<
" x " << nvar);
142 it_info(
"Variable node degree distribution from node perspective:\n"
144 it_info(
"Check node degree distribution from node perspective:\n"
146 it_info(
"Variable node degree distribution from edge perspective:\n"
147 << vdegedge / edges);
148 it_info(
"Check node degree distribution from edge perspective:\n"
149 << cdegedge / edges);
151 it_info(
"--------------------------------");
163 alist.
write(alist_file);
173 for (
int i = 0; i <
ncheck; i++) {
174 for (
int j = 0; j <
nvar; j++) {
185 "LDPC_Parity::export_alist(): Object not initialized");
193 int to_j,
int godir,
int L)
const
196 "LDPC_Parity::check_connectivity(): Object not initialized");
204 if ((from_i == to_i) && (from_j == to_j) && (godir != 0)) {
208 if (
get(from_i, from_j) == 0) {
214 if (
get(from_i, to_j) == 1) {
return 0; }
217 if (
get(to_i, from_j) == 1) {
return 0; }
222 if ((godir == 1) || (godir == 0)) {
224 for (i = 0; i <
length(cj); i++) {
225 if (cj(i) != from_i) {
236 for (j = 0; j <
length(ri); j++) {
237 if (ri(j) != from_j) {
252 "LDPC_Parity::check_for_cycles(): Object not initialized");
254 if ((L&1) == 1) {
return (-1); }
255 if (L == 0) {
return (-4); }
258 for (
int i = 0; i <
nvar; i++) {
260 for (
int j = 0; j <
length(ri); j++) {
295 "LDPC_Parity::cycle_removal_MGW(): Object not initialized");
303 for (
int j = 0; j <
nvar; j++) {
305 for (
int i = 0; i < col.
nnz(); i++) {
317 Gpow(0).set(i, i, 1);
326 int cycles_found = 0;
328 for (r = 4; r <= Maxcyc; r += 2) {
330 Gpow(r / 2) = Gpow(r / 2 - 1) * G;
333 traverse_again =
false;
335 it_info_debug(
"Starting new pass of cycle elimination, target girth "
336 << (r + 2) <<
"...");
342 if (ptemp > pdone + 10) {
347 if (((Gpow(r / 2))(i, j) >= 2) && ((Gpow(r / 2 - 2))(i, j) == 0)) {
353 G.get_col(j))).get_nz_indices();
357 "LDPC_Parity_Unstructured::cycle_removal_MGW(): "
361 Ssvec rowjk = Gpow(r / 2) * (Gpow(r / 2 - 1).get_col(j)
362 + Gpow(r / 2 - 1).get_col(k));
366 for (
int s = 0; s < nvar +
ncheck; s++) {
368 if (rowjk(l) != 0) {
continue; }
369 ivec colcandi = G.get_col(l).get_nz_indices();
370 if (
length(colcandi) > 0) {
372 for (
int u = 0; u <
length(colcandi); u++) {
376 goto found_candidate_vector;
384 found_candidate_vector:
387 if (p >= ncheck) {
int z = l; l = p; p = z; }
388 if (j >= ncheck) {
int z = k; k = j; j = z; }
397 it_assert_debug((
get(j, k - ncheck) == 1) && (
get(p, l - ncheck) == 1),
398 "LDPC_Parity_Unstructured::cycle_removal_MGW(): "
402 it_assert_debug((
get(j, l - ncheck) == 0) && (
get(p, k - ncheck) == 0),
403 "LDPC_Parity_Unstructured::cycle_removal_MGW(): "
410 && G(k, j) == 1,
"G");
412 && G(k, p) == 0,
"G");
415 Ssmat Delta(ncheck + nvar, ncheck + nvar, 2);
428 && G(k, j) == 0,
"G");
430 && G(k, p) == 1,
"G");
434 for (
int z = 3; z <= r / 2; z++) {
435 Gpow(z) = Gpow(z - 1) * G;
438 traverse_again =
true;
442 if ((!traverse_again) && (cycles_found > 0)) {
447 while (cycles_found != 0);
450 <<
". Proceeding to next level.");
455 it_info_debug(
"Cycle removal (MGW algoritm) finished. Graph girth: "
456 << girth <<
". Cycles remaining on next girth level: "
475 for (
int i = 0;i < C.length();i++) {
476 for (
int j = 0; j < C(i); j++) {
477 for (
int m = 0; m < i; m++) Ne++;
488 for (
int i = 0;i < C.length();i++) {
489 for (
int j = 0; j < C(i); j++) {
490 for (
int m = 0; m < i; m++) {
500 for (
int i = 0;i < R.length();i++) {
501 for (
int j = 0; j < R(i); j++) {
502 for (
int m = 0; m < i; m++) {
513 ivec ind = sort_index(
randu(Ne));
519 for (
int i = 0; i <
length(cycopt); i++) {
520 Laim(i + 2) = cycopt(i);
522 for (
int i =
length(cycopt); i <
Nmax - 2; i++) {
523 Laim(i + 2) = cycopt(
length(cycopt) - 1);
528 const int Max_attempts = 100;
530 for (
int k = 0; k < Ne; k++) {
531 const int el = Ne - k - 2;
534 <<
". Variable node degree: " << vd(vcon(k))
535 <<
". Girth target: " << Laim(vd(vcon(k)))
536 <<
". Accumulated failures: " << failures);
538 const int c = cp(vcon(k));
539 int L = Laim(vd(vcon(k)));
542 if (attempt > 0 && attempt % apcl == 0 && L >= 6) { L -= 2; };
543 int r = rp(ccon(ind(k)));
548 int t = k + 1 +
randi(0, el - 1);
553 if (attempt == Max_attempts) {
570 int t = k + 1 +
randi(0, el - 1);
575 if (attempt == Max_attempts) {
617 if (
sum(C) != Nvar) {
619 C(ind(0)) = C(ind(0)) - (
sum(C) - Nvar);
631 R.set(ind(0), old - 1);
633 R.set(ind(0) - 1, old + 1);
637 if (ind(0) == R.length() - 1) {
642 R.set(ind(0), old - 1);
644 R.set(ind(0) + 1, old + 1);
661 const std::string& method,
664 generate(Nvar, k, l, method, options);
668 const std::string& method,
671 vec var_deg =
zeros(k);
672 vec chk_deg =
zeros(l);
681 if (method ==
"rand") {
697 const std::string& method,
700 generate(Nvar, var_deg, chk_deg, method, options);
705 const std::string& method,
711 if (method ==
"rand") {
741 for (
int r = 0; r < H_b.rows(); r++) {
742 for (
int c = 0; c < H_b.cols(); c++) {
749 for (
int i = 0; i < Z; ++i)
750 set(rz + i, cz + i, 1);
753 for (
int i = 0; i < Z; ++i)
754 set(rz + i, cz + (i + H_b(r, c)) % Z, 1);
779 calculate_base_matrix();
786 std::ifstream bm_file(filename.c_str());
787 it_assert(bm_file.is_open(),
"BLDPC_Parity::load_base_matrix(): Could not "
788 "open file \"" << filename <<
"\" for reading");
793 int line_counter = 0;
794 getline(bm_file, line);
795 while (!bm_file.eof()) {
797 std::stringstream ss(line);
804 if ((H_b.rows() == 0) || (row.size() == H_b.cols()))
807 it_warning(
"BLDPC_Parity::load_base_matrix(): Wrong size of "
808 "a parsed row number " << line_counter);
809 getline(bm_file, line);
814 if (H_b.rows() > H_b.cols())
815 H_b = H_b.transpose();
823 it_assert(H_b_valid,
"BLDPC_Parity::save_base_matrix(): Base matrix is "
825 std::ofstream bm_file(filename.c_str());
826 it_assert(bm_file.is_open(),
"BLDPC_Parity::save_base_matrix(): Could not "
827 "open file \"" << filename <<
"\" for writing");
829 for (
int r = 0; r < H_b.rows(); r++) {
830 for (
int c = 0; c < H_b.cols(); c++) {
831 bm_file << std::setw(3) << H_b(r, c);
840 void BLDPC_Parity::calculate_base_matrix()
842 std::string error_str =
"BLDPC_Parity::calculate_base_matrix(): "
843 "Invalid BLDPC matrix. Cannot calculate base matrix from it.";
844 int rows =
H.
rows() / Z;
845 int cols =
H.
cols() / Z;
847 H_b.set_size(rows, cols);
849 for (
int r = 0; r < rows; ++r) {
851 for (
int c = 0; c < cols; ++c) {
855 if (H_Z.
nnz() == 0) {
858 else if (H_Z.
nnz() == Z) {
861 while ((shift < Z) && (H_Z(0, shift) != 1))
864 for (
int i = 1; i < Z; ++i)
865 it_assert(H_Z(0, shift) == H_Z(i, (i + shift) % Z), error_str);
874 it_info(
"Base matrix calculated");
885 bool natural_ordering,
890 tmp =
construct(H, natural_ordering, ind);
895 bool natural_ordering,
896 const ivec& avoid_cols)
905 ivec col_order(nvar);
906 if (natural_ordering) {
907 for (
int i = 0; i < nvar; i++) {
913 vec col_importance =
randu(nvar);
914 for (
int i = 0; i <
length(avoid_cols); i++) {
915 col_importance(avoid_cols(i)) = (-col_importance(avoid_cols(i)));
917 col_order = sort_index(-col_importance);
920 ivec actual_ordering(nvar);
934 for (
int k = 0; k < nvar; k++) {
935 it_error_if(j1 >= nvar - ncheck,
"LDPC_Generator_Systematic::construct(): "
936 "Unable to obtain enough independent columns.");
938 bvec c = Hd.get_col(col_order(k));
941 rank = P2.
T_fact(T, U, perm);
942 actual_ordering(k) = nvar - ncheck;
949 actual_ordering(k) = nvar - ncheck + j2;
956 actual_ordering(k) = j1;
960 actual_ordering(k) = j1;
975 for (
int i = 0; i < ncheck; i++) {
976 for (
int j = 0; j < nvar; j++) {
986 "LDPC_Generator_Systematic::construct(): Incorrect generator matrix G");
993 return actual_ordering;
1001 f >>
Name(
"Fileversion") >> ver;
1003 "LDPC_Generator_Systematic::save(): Unsupported file format");
1005 f <<
Name(
"G") << G;
1014 f >>
Name(
"Fileversion") >> ver;
1016 "LDPC_Generator_Systematic::load(): Unsupported file format");
1017 std::string gen_type;
1018 f >>
Name(
"G_type") >> gen_type;
1020 "LDPC_Generator_Systematic::load(): Wrong generator type");
1021 f >>
Name(
"G") >> G;
1031 "generator not set up");
1032 it_assert(input.size() == G.
cols(),
"LDPC_Generator_Systematic::encode(): "
1033 "Improper input vector size (" << input.size() <<
" != "
1034 << G.
cols() <<
")");
1036 output =
concat(input, G * input);
1045 const std::string type):
1055 "initialized generator");
1056 it_assert(input.size() ==
K,
"BLDPC_Generator::encode(): Input vector "
1057 "length is not equal to K");
1061 output.set_size(
N,
true);
1064 for (
int k = 0; k <
Z; k++) {
1065 for (
int j = 0; j <
K; j++) {
1066 output(K + k) +=
H_enc(
M - 1 - k, j) * input(j);
1068 for (
int j = 0; j < k; j++) {
1069 output(K + k) +=
H_enc(
M - 1 - k, K + j) * output(K + j);
1074 for (
int k = 0; k <
M -
Z; k++) {
1075 for (
int j = 0; j <
K; j++) {
1076 output(K + Z + k) +=
H_enc(k, j) * input(j);
1078 for (
int j = K; j < K +
Z; j++) {
1079 output(K + Z + k) +=
H_enc(k, j) * output(j);
1081 for (
int j = K + Z; j < K + Z + k; j++) {
1082 output(K + Z + k) +=
H_enc(k, j) * output(j);
1101 for (
int i = 0; i < M -
Z; i +=
Z) {
1103 for (
int j = 0; j <
Z; j++) {
1115 for (
int c1 =
K + Z - 1; c1 >=
K; c1--) {
1118 while (
H_enc(r2, c1) == 0 && r2 < M - 1)
1125 for (r2 = r1 + 1; r2 <
M; r2++) {
1126 if (
H_enc(r2, c1) == 1) {
1144 "BLDPC_Generator::save(): Can not save not initialized generator");
1147 for (
int i = 0; i <
M / Z - 1; i++) {
1155 f >>
Name(
"Fileversion") >> ver;
1157 "Unsupported file format");
1159 f <<
Name(
"H_T") << H_T;
1160 f <<
Name(
"H_Z") << H_Z;
1161 f <<
Name(
"Z") <<
Z;
1171 f >>
Name(
"Fileversion") >> ver;
1173 "Unsupported file format");
1174 std::string gen_type;
1175 f >>
Name(
"G_type") >> gen_type;
1177 "BLDPC_Generator::load(): Wrong generator type");
1178 f >>
Name(
"H_T") >> H_T;
1179 f >>
Name(
"H_Z") >> H_Z;
1180 f >>
Name(
"Z") >>
Z;
1184 M = (H_T.rows() + 1) * Z;
1188 for (
int i = 0; i < H_T.rows(); i++) {
1189 for (
int j = 0; j <
Z; j++) {
1190 for (
int k = 0; k <
N; k++) {
1191 if (H_T(i, (k / Z)*Z + (k + Z - j) % Z)) {
1208 max_iters(50), psc(true), pisc(false),
1213 H_defined(false), G_defined(false), dec_method(
"BP"), max_iters(50),
1221 H_defined(false), G_defined(false), dec_method(
"BP"), max_iters(50),
1243 it_info_debug(
"LDPC_Code::load_code(): Loading LDPC codec from "
1248 f >>
Name(
"Fileversion") >> ver;
1250 "Unsupported file format");
1255 f >>
Name(
"C") >> C;
1256 f >>
Name(
"V") >> V;
1257 f >>
Name(
"sumX1") >> sumX1;
1258 f >>
Name(
"sumX2") >> sumX2;
1259 f >>
Name(
"iind") >> iind;
1260 f >>
Name(
"jind") >> jind;
1265 it_assert(G_in != 0,
"LDPC_Code::load_code(): Generator object is "
1266 "missing. Can not load the generator data from a file.");
1272 it_info_debug(
"LDPC_Code::load_code(): Generator data not loaded. "
1273 "Generator object will not be used.");
1276 it_info_debug(
"LDPC_Code::load_code(): Successfully loaded LDPC codec "
1277 "from " << filename);
1286 it_info_debug(
"LDPC_Code::save_to_file(): Saving LDPC codec to "
1290 f.
open(filename,
true);
1296 f <<
Name(
"C") << C;
1297 f <<
Name(
"V") << V;
1298 f <<
Name(
"sumX1") << sumX1;
1299 f <<
Name(
"sumX2") << sumX2;
1300 f <<
Name(
"iind") << iind;
1301 f <<
Name(
"jind") << jind;
1308 it_info_debug(
"LDPC_Code::save_code(): Missing generator object - "
1309 "generator data not saved");
1311 it_info_debug(
"LDPC_Code::save_code(): Successfully saved LDPC codec to "
1318 it_assert((method_in ==
"bp") || (method_in ==
"BP"),
1319 "LDPC_Code::set_decoding_method(): Not implemented decoding method");
1324 bool syndr_check_each_iter,
1325 bool syndr_check_at_start)
1328 "number of iterations can not be negative");
1330 psc = syndr_check_each_iter;
1331 pisc = syndr_check_at_start;
1361 syst_bits = (qllrout.left(
nvar -
ncheck) < 0);
1367 decode(llr_in, syst_bits);
1393 && (sumX2.size() ==
ncheck),
"LDPC_Code::bp_decode(): Wrong "
1394 "input dimensions");
1401 LLRout.set_size(LLRin.size());
1406 QLLRvec ml(max_cnd);
1407 QLLRvec mr(max_cnd);
1410 for (
int i = 0; i <
nvar; i++) {
1412 for (
int j = 0; j < sumX1(i); j++) {
1413 mvc[index] = LLRin(i);
1418 bool is_valid_codeword =
false;
1424 for (
int j = 0; j <
ncheck; j++) {
1429 it_error(
"LDPC_Code::bp_decode(): sumX2(j)=0");
1431 it_error(
"LDPC_Code::bp_decode(): sumX2(j)=1");
1433 mcv[j+
ncheck] = mvc[jind[j]];
1434 mcv[j] = mvc[jind[j+
ncheck]];
1439 QLLR m0 = mvc[jind[j0]];
1441 QLLR m1 = mvc[jind[j1]];
1443 QLLR m2 = mvc[jind[j2]];
1451 QLLR m0 = mvc[jind[j0]];
1453 QLLR m1 = mvc[jind[j1]];
1455 QLLR m2 = mvc[jind[j2]];
1457 QLLR m3 = mvc[jind[j3]];
1468 QLLR m0 = mvc[jind[j0]];
1470 QLLR m1 = mvc[jind[j1]];
1472 QLLR m2 = mvc[jind[j2]];
1474 QLLR m3 = mvc[jind[j3]];
1476 QLLR m4 = mvc[jind[j4]];
1490 QLLR m0 = mvc[jind[j0]];
1492 QLLR m1 = mvc[jind[j1]];
1494 QLLR m2 = mvc[jind[j2]];
1496 QLLR m3 = mvc[jind[j3]];
1498 QLLR m4 = mvc[jind[j4]];
1500 QLLR m5 = mvc[jind[j5]];
1516 int nodes = sumX2(j);
1517 if( nodes > max_cnd ) {
1518 std::ostringstream m_sout;
1519 m_sout <<
"check node degrees >" << max_cnd <<
" not supported in this version";
1525 m[0] = mvc[jind[jj[0]]];
1526 for(
int i = 1; i <= nodes; i++ ) {
1527 jj[i] = jj[i-1] +
ncheck;
1528 m[i] = mvc[jind[jj[i]]];
1534 for(
int i = 1; i < nodes; i++ ) {
1540 mcv[jj[0]] = mr[nodes-1];
1541 mcv[jj[nodes]] = ml[nodes-1];
1542 for(
int i = 1; i < nodes; i++ )
1549 for (
int i = 0; i <
nvar; i++) {
1552 it_error(
"LDPC_Code::bp_decode(): sumX1(i)=0");
1557 LLRout(i) = LLRin(i);
1561 QLLR m0 = mcv[iind[i]];
1563 QLLR m1 = mcv[iind[i1]];
1564 mvc[i] = LLRin(i) + m1;
1565 mvc[i1] = LLRin(i) + m0;
1566 LLRout(i) = mvc[i1] + m1;
1571 QLLR m0 = mcv[iind[i0]];
1573 QLLR m1 = mcv[iind[i1]];
1575 QLLR m2 = mcv[iind[i2]];
1576 LLRout(i) = LLRin(i) + m0 + m1 + m2;
1577 mvc[i0] = LLRout(i) - m0;
1578 mvc[i1] = LLRout(i) - m1;
1579 mvc[i2] = LLRout(i) - m2;
1584 QLLR m0 = mcv[iind[i0]];
1586 QLLR m1 = mcv[iind[i1]];
1588 QLLR m2 = mcv[iind[i2]];
1590 QLLR m3 = mcv[iind[i3]];
1591 LLRout(i) = LLRin(i) + m0 + m1 + m2 + m3;
1592 mvc[i0] = LLRout(i) - m0;
1593 mvc[i1] = LLRout(i) - m1;
1594 mvc[i2] = LLRout(i) - m2;
1595 mvc[i3] = LLRout(i) - m3;
1599 QLLR mvc_temp = LLRin(i);
1601 for (
int jp = 0; jp < sumX1(i); jp++) {
1602 mvc_temp += mcv[iind[index_iind]];
1605 LLRout(i) = mvc_temp;
1607 for (
int j = 0; j < sumX1[i]; j++) {
1608 mvc[index_iind] = mvc_temp - mcv[iind[index_iind]];
1616 is_valid_codeword =
true;
1623 return (is_valid_codeword ? iter : -iter);
1629 QLLRvec llr = 1 - 2 *
to_ivec(x);
1639 for (j = 0; j <
ncheck; j++) {
1642 for (i = 0; i < sumX2(j); i++) {
1649 if ((synd&1) == 1) {
1664 sumX1 = Hmat->
sumX1;
1665 sumX2 = Hmat->
sumX2;
1668 int cmax =
max(sumX1);
1669 int vmax =
max(sumX2);
1677 it_info_debug(
"LDPC_Code::decoder_parameterization(): Computations "
1679 for (
int i = 0; i <
nvar; i++) {
1681 for (
int j0 = 0; j0 <
length(coli); j0++) {
1682 C(j0 + cmax*i) = coli(j0);
1686 it_info_debug(
"LDPC_Code::decoder_parameterization(): Computations "
1688 it_info_debug(
"Computing decoder parameterization. Phase 2");
1689 for (
int j = 0; j <
ncheck; j++) {
1691 for (
int i0 = 0; i0 <
length(rowj); i0++) {
1692 V(j + ncheck*i0) = rowj(i0);
1696 it_info_debug(
"LDPC_Code::decoder_parameterization(): Computations "
1698 it_info_debug(
"Computing decoder parameterization. Phase 3.");
1699 for (
int j = 0; j <
ncheck; j++) {
1700 for (
int ip = 0; ip < sumX2(j); ip++) {
1701 int vip = V(j + ip * ncheck);
1704 if (C(k + vip*cmax) == j) {
1709 jind(j + ip*ncheck) = vip + k *
nvar;
1713 it_info_debug(
"LDPC_Code::decoder_parameterization(): Computations "
1715 for (
int i = 0; i <
nvar; i++) {
1716 for (
int jp = 0; jp < sumX1(i); jp++) {
1717 int cjp = C(jp + i * cmax);
1720 if (V(cjp + k*ncheck) == i) {
break; }
1723 iind(i + jp*nvar) = cjp + k *
ncheck;
1735 mvc.set_size(
max(sumX1) *
nvar);
1743 it_info_debug(
"LDPC_Code::integrity_check(): Checking integrity of "
1744 "the LDPC_Parity and LDPC_Generator data");
1751 "LDPC_Code::integrity_check(): Syndrome check failed");
1752 bv.shift_right(bv(
nvar - ncheck - 1));
1756 it_info_debug(
"LDPC_Code::integrity_check(): No generator defined "
1757 "- no check performed");
1768 for (
int i = 0; i < C.
ncheck; i++) {
1773 for (
int j = 0; j < C.
nvar; j++) {
1777 os <<
"--- LDPC codec ----------------------------------\n"
1778 <<
"Nvar : " << C.
get_nvar() <<
"\n"
1780 <<
"Rate : " << C.
get_rate() <<
"\n"
1781 <<
"Column degrees (node perspective): " << cdeg <<
"\n"
1782 <<
"Row degrees (node perspective): " << rdeg <<
"\n"
1783 <<
"-------------------------------------------------\n"
1784 <<
"Decoder parameters:\n"
1786 <<
" - max. iterations : " << C.
max_iters <<
"\n"
1787 <<
" - syndrome check at each iteration : " << C.
psc <<
"\n"
1788 <<
" - syndrome check at start : " << C.
pisc <<
"\n"
1789 <<
"-------------------------------------------------\n"