diff --git a/.gitignore b/.gitignore index fec6447..987729f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,7 @@ latex/build/ latex/tmp/ +sw/cpp_modules +sw/cpp/build sw/sim_saves/ .idea __pycache__ diff --git a/sw/cpp/.clang-format b/sw/cpp/.clang-format new file mode 100644 index 0000000..4b9f225 --- /dev/null +++ b/sw/cpp/.clang-format @@ -0,0 +1,28 @@ +BasedOnStyle: LLVM +Language: Cpp + +IndentWidth: 4 +UseTab: Never +#NamespaceIndentation: All + +PointerAlignment: Left +AccessModifierOffset: -4 +AlwaysBreakTemplateDeclarations: true +LambdaBodyIndentation: Signature + +MaxEmptyLinesToKeep: 3 +# ColumnLimit: 128 + +CompactNamespaces: true +FixNamespaceComments: true + +AllowShortFunctionsOnASingleLine: false +AllowShortIfStatementsOnASingleLine: true + +AlignConsecutiveAssignments: true +AlignConsecutiveBitFields: true +AlignConsecutiveDeclarations: true +AlignConsecutiveMacros: true + +#BraceWrapping: +# BeforeElse: true diff --git a/sw/cpp/CMakeLists.txt b/sw/cpp/CMakeLists.txt new file mode 100644 index 0000000..95869fc --- /dev/null +++ b/sw/cpp/CMakeLists.txt @@ -0,0 +1,18 @@ +cmake_minimum_required (VERSION 3.0) +project(cpp_decoders) + +find_package(Eigen3 3.3 REQUIRED NO_MODULE) +find_package(pybind11 CONFIG REQUIRED) + +include_directories(${pybind11_INCLUDE_DIRS}) + + +pybind11_add_module(cpp_decoders src/cpp_decoders.cpp) +target_link_libraries(cpp_decoders PRIVATE Eigen3::Eigen) + +set(INSTALL_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../cpp_modules) + +install(TARGETS cpp_decoders ARCHIVE DESTINATION ${INSTALL_DIR} + LIBRARY DESTINATION ${INSTALL_DIR} + RUNTIME DESTINATION ${INSTALL_DIR}) + diff --git a/sw/cpp/src/cpp_decoders.cpp b/sw/cpp/src/cpp_decoders.cpp new file mode 100644 index 0000000..4a50b2d --- /dev/null +++ b/sw/cpp/src/cpp_decoders.cpp @@ -0,0 +1,156 @@ +#include +#include +#include +#include + +#include + + +/* + + 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; + +using MatrixXdR = + Eigen::Matrix; + + +class ProximalDecoder { +public: + ProximalDecoder(const Eigen::Ref& 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& 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().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 mH_zero_indices; + + + static Eigen::RowVectorXd L_awgn(const Eigen::RowVectorXd& s, + const Eigen::RowVectorXd& y) { + return s.array() - y.array(); + } + + static std::vector find_zero(MatrixXiR mat) { + std::vector 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(); + + 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_(proximal, "ProximalDecoder") + .def(pybind11::init(), + "H"_a.noconvert(), "K"_a, "omega"_a, "gamma"_a, "eta"_a) + .def("decode", &ProximalDecoder::decode, "x"_a.noconvert()); +} \ No newline at end of file