156 lines
4.3 KiB
C++
156 lines
4.3 KiB
C++
#include <Eigen/Cholesky>
|
|
#include <iostream>
|
|
#include <pybind11/numpy.h>
|
|
#include <pybind11/pybind11.h>
|
|
|
|
#include <pybind11/eigen.h>
|
|
|
|
|
|
/*
|
|
|
|
import numpy as np
|
|
|
|
|
|
class ProximalDecoder:
|
|
"""Class implementing the Proximal Decoding algorithm. See "Proximal
|
|
Decoding for LDPC Codes"
|
|
by Tadashi Wadayama, and Satoshi Takabe.
|
|
"""
|
|
|
|
def decode(self, y: np.array) -> np.array:
|
|
"""Decode a received signal. The algorithm is detailed in 3.2, p.3.
|
|
|
|
This function assumes a BPSK modulated signal and an AWGN channel.
|
|
|
|
:param y: Vector of received values. (y = x + w, where 'x' is
|
|
element of [-1, 1]^n and 'w' is noise)
|
|
:return: Most probably sent codeword (element of [0, 1]^n). If
|
|
decoding fails, the returned value is 'None'
|
|
"""
|
|
s = np.zeros(self._n)
|
|
x_hat = np.zeros(self._n)
|
|
for k in range(self._K):
|
|
r = s - self._step_size * self._L_awgn(s, y)
|
|
|
|
s = r - self._gamma * self._grad_h(r)
|
|
s = self._projection(s) # Equation (15)
|
|
|
|
x_hat = np.sign(s)
|
|
|
|
# Map the codeword from [ -1, 1]^n to [0, 1]^n
|
|
x_hat = (x_hat == -1) * 1
|
|
|
|
if self._check_parity(x_hat):
|
|
return x_hat
|
|
|
|
return None
|
|
|
|
|
|
* */
|
|
|
|
namespace py11 = pybind11;
|
|
using namespace pybind11::literals;
|
|
|
|
|
|
using MatrixXiR =
|
|
Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
|
|
|
|
using MatrixXdR =
|
|
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
|
|
|
|
|
|
class ProximalDecoder {
|
|
public:
|
|
ProximalDecoder(const Eigen::Ref<const MatrixXiR>& H, int K, double omega,
|
|
double gamma, double eta)
|
|
: mN(H.cols()), mK(K), mOmega(omega), mGamma(gamma), mEta(eta), mH(H),
|
|
mH_zero_indices(find_zero(H)) {
|
|
}
|
|
|
|
const Eigen::RowVectorXi
|
|
decode(const Eigen::Ref<const Eigen::VectorXd>& y) {
|
|
Eigen::RowVectorXd s = Eigen::RowVectorXd::Zero(mH.cols());
|
|
Eigen::RowVectorXi x_hat;
|
|
Eigen::RowVectorXd r;
|
|
|
|
for (std::size_t i = 0; i < mK; ++i) {
|
|
r = s - mOmega * L_awgn(s, y);
|
|
s = projection(r - mGamma * grad_H(r));
|
|
|
|
x_hat = s.cast<int>().cwiseSign();
|
|
x_hat = (x_hat.array() - 1).matrix() / (-2);
|
|
|
|
if (check_parity(x_hat)) {
|
|
return x_hat;
|
|
}
|
|
}
|
|
|
|
return x_hat;
|
|
// TODO: Return 'None'
|
|
}
|
|
|
|
private:
|
|
const int mN;
|
|
const int mK;
|
|
const double mOmega;
|
|
const double mGamma;
|
|
const double mEta;
|
|
|
|
const MatrixXiR mH;
|
|
|
|
const std::vector<Eigen::Index> mH_zero_indices;
|
|
|
|
|
|
static Eigen::RowVectorXd L_awgn(const Eigen::RowVectorXd& s,
|
|
const Eigen::RowVectorXd& y) {
|
|
return s.array() - y.array();
|
|
}
|
|
|
|
static std::vector<Eigen::Index> find_zero(MatrixXiR mat) {
|
|
std::vector<Eigen::Index> indices;
|
|
|
|
for (Eigen::Index i = 0; i < mat.size(); ++i)
|
|
if (mat(i) == 0) indices.push_back(i);
|
|
|
|
return indices;
|
|
}
|
|
|
|
Eigen::RowVectorXd grad_H(const Eigen::RowVectorXd& x) {
|
|
MatrixXdR A_prod_matrix = x.replicate(mH.rows(), 1);
|
|
|
|
for (const auto& index : mH_zero_indices)
|
|
A_prod_matrix(index) = 1;
|
|
MatrixXdR A_prods = A_prod_matrix.rowwise().prod();
|
|
|
|
|
|
Eigen::RowVectorXd B_sums =
|
|
(A_prods.array().pow(2) - A_prods.array()).matrix().transpose();
|
|
B_sums = B_sums * mH.cast<double>();
|
|
|
|
Eigen::RowVectorXd result = 4 * (x.array().pow(2) - 1) * x.array() +
|
|
(2 * x.array().inverse()) * B_sums.array();
|
|
|
|
return result;
|
|
}
|
|
|
|
bool check_parity(const Eigen::RowVectorXi& x_hat) {
|
|
Eigen::RowVectorXi syndrome =
|
|
(mH * x_hat.transpose()).unaryExpr([](int i) { return i % 2; });
|
|
|
|
return !(syndrome.count() > 0);
|
|
}
|
|
|
|
Eigen::RowVectorXd projection(const Eigen::RowVectorXd& v) {
|
|
return v.cwiseMin(mEta).cwiseMax(-mEta);
|
|
}
|
|
};
|
|
|
|
|
|
PYBIND11_MODULE(cpp_decoders, proximal) {
|
|
proximal.doc() = "Proximal decoder";
|
|
|
|
pybind11::class_<ProximalDecoder>(proximal, "ProximalDecoder")
|
|
.def(pybind11::init<MatrixXiR, int, double, double, double>(),
|
|
"H"_a.noconvert(), "K"_a, "omega"_a, "gamma"_a, "eta"_a)
|
|
.def("decode", &ProximalDecoder::decode, "x"_a.noconvert());
|
|
} |