Restructure C++ for new version of paper and implement other estimators.

cpp/CMakeLists.txt

@@ -1,17 +1,20 @@
 cmake_minimum_required(VERSION 3.18...3.31)
 project(
   ${SKBUILD_PROJECT_NAME}
-  DESCRIPTION "Deterministic Kaczmarz Algorithm for Online Paramter Estimation."
+  DESCRIPTION "Online estimation methods."
   VERSION ${SKBUILD_PROJECT_VERSION}
   LANGUAGES CXX
 )
 
 set(
   SOURCES
   src/bindings.cpp
-  src/estimators/deka.cpp
+  src/estimators/grk.cpp
+  src/estimators/kf.cpp
+  src/estimators/rk.cpp
   src/estimators/rls.cpp
-  src/estimators/ekf.cpp
+  src/estimators/tagrk.cpp
+  src/estimators/tark.cpp
 )
 
 set(PYBIND11_FINDPYTHON ON)

cpp/README.md

@@ -15,5 +15,28 @@ pip install .
 ## Usage
 The estimators can be imported like python classes:
 ```python
-from deka_cpp._core import RLS, EKF, DEKA
+from online_estimation._core import RLS, RLS_Robust, KF, KF_Robust, RK, GRK, TARK, TAGRK
+```
+
+Each estimator exposes an iterate method. For example, with the RLS estimator:
+```python
+n = 5
+
+A = ...
+b = ...
+x0 = np.zeros((n, 1))
+
+kwargs = {
+    "lambda": 0.99,
+    "p_coeff": 1000,
+    "x0": None
+}
+estimator = RLS(n, **kwargs)
+
+x_estimate = estimator.iterate(A, b, x0=x0)
+```
+
+You can also modify the internal state without passing it to the iterate function:
+```python
+estimator.state = x_new  # performs check on if dimensions match current state
 ```

cpp/pyproject.toml

@@ -3,9 +3,9 @@ requires = ["scikit-build-core", "pybind11", "numpy"]
 build-backend = "scikit_build_core.build"
 
 [project]
-name = "deka_cpp"
+name = "online_estimation"
 version = "0.0.1"
-description = "C++ bindings for DeKA."
+description = "Python bindings for online_estimation methods."
 readme = "README.md"
 requires-python = ">=3.8"
 dependencies = []

cpp/src/bindings.cpp

@@ -2,83 +2,83 @@
 #include <pybind11/eigen/matrix.h>
 #include <pybind11/stl.h>
 
+#include "estimator_types.h"
+#include "estimators/estimator.h"
+#include "estimators/grk.h"
+#include "estimators/kf.h"
+#include "estimators/rk.h"
 #include "estimators/rls.h"
-#include "estimators/deka.h"
-#include "estimators/ekf.h"
+#include "estimators/tagrk.h"
+#include "estimators/tark.h"
 
 #define STRINGIFY(x) #x
 #define MACRO_STRINGIFY(x) STRINGIFY(x)
 
 namespace py = pybind11;
 
-
 PYBIND11_MODULE(_core, m) {    
     m.doc() = R"pbdoc(
-        DEKA Estimation Methods
+        Online estimation methods.
         -----------------------
 
-        .. currentmodule:: deka
+        .. currentmodule:: online_estimation
 
         .. autosummary::
          :toctree: _generate
-
-         DEKA
     )pbdoc";
 
-    py::class_<deka::DEKA>(m, "DEKA")
-        .def(py::init<const int, const double, const double, const double,
-                const double, const double>(),
-             py::arg("n"),
-             py::arg("damping") = 0.1,
-             py::arg("reg") = 1e-6,
-             py::arg("smoothing") = 0.9,
-             py::arg("tol_max") = 1e-4,
-             py::arg("tol_min") = 1e-7
-        )
-        .def("iterate", &deka::DEKA::iterate, py::arg("A"), py::arg("b"),
-             py::arg("x0") = py::none(), py::arg("num_iters") = 1000)
-        .def_property_readonly("x", [](const deka::DEKA &self) { return self.x; })
-        .def_property_readonly("x_smooth", [](const deka::DEKA &self) { return self.x_smooth; });
-
-    py::class_<deka::RLS>(m, "RLS")
-        .def(py::init<const int, const double, const double>(),
-             py::arg("n"),
-             py::arg("lambda"),
-             py::arg("coeff")
-        )
-        .def("iterate", &deka::RLS::iterate, py::arg("A"), py::arg("b"),
-             py::arg("x0") = py::none(), py::arg("num_iters") = 1000)
-        .def_property_readonly("x", [](const deka::RLS &self) { return self.x; });
-
-    py::class_<deka::RLS_Robust>(m, "RLS_Robust")
-        .def(py::init<const int, const double, const double>(),
-             py::arg("n"),
-             py::arg("lambda"),
-             py::arg("coeff")
-        )
-        .def("iterate", &deka::RLS_Robust::iterate, py::arg("A"), py::arg("b"),
-             py::arg("x0") = py::none(), py::arg("num_iters") = 1000)
-        .def_property_readonly("x", [](const deka::RLS_Robust &self) { return self.x; });
-
-    py::class_<deka::EKF>(m, "EKF")
-        .def(py::init<int, const Eigen::MatrixXd&, const Eigen::MatrixXd&, double>(),
-             py::arg("n"), py::arg("process_noise"),
-             py::arg("meas_noise"), py::arg("coeff"))
-        .def("iterate", &deka::EKF::iterate, py::arg("A"), py::arg("b"),
-             py::arg("x0") = py::none(), py::arg("num_iters") = 1)
-        .def_property_readonly("x", [](const deka::EKF &self) { return self.x; })
-        .def_property_readonly("P", [](const deka::EKF &self) { return self.P; });
-
-    py::class_<deka::EKF_Robust>(m, "EKF_Robust")
-        .def(py::init<int, const Eigen::MatrixXd&, const Eigen::MatrixXd&, double>(),
-             py::arg("n"), py::arg("process_noise"),
-             py::arg("meas_noise"), py::arg("coeff"))
-        .def("iterate", &deka::EKF_Robust::iterate, py::arg("A"), py::arg("b"),
-             py::arg("x0") = py::none(), py::arg("num_iters") = 1)
-        .def_property_readonly("x", [](const deka::EKF_Robust &self) { return self.x; })
-        .def_property_readonly("P", [](const deka::EKF_Robust &self) { return self.P; });
-
-
+    /// Estimator interface
+    py::class_<estim::Estimator>(m, "Estimator")
+        .def_property("state", &estim::Estimator::get_state, &estim::Estimator::set_state);
+
+    /// Estimator implementations
+    py::class_<estim::RLS, estim::Estimator>(m, "RLS")
+        .def(py::init<int, double, double, const std::optional<estim::RowMatrixXd>&>(),
+             py::arg("n"), py::arg("lambda"), py::arg("p_coeff"), py::arg("x0"))
+        .def("iterate", &estim::RLS::iterate, py::arg("A"), py::arg("b"), py::arg("x0") = std::nullopt,
+             py::return_value_policy::reference_internal);
+
+    py::class_<estim::RLS_Robust, estim::Estimator>(m, "RLS_Robust")
+        .def(py::init<int, double, double, const std::optional<estim::RowMatrixXd>&>(),
+             py::arg("n"), py::arg("lambda"), py::arg("p_coeff"), py::arg("x0"))
+        .def("iterate", &estim::RLS_Robust::iterate, py::arg("A"), py::arg("b"), py::arg("x0") = std::nullopt,
+             py::return_value_policy::reference_internal);
+
+    py::class_<estim::KF, estim::Estimator>(m, "KF")
+        .def(py::init<int, const estim::RowMatrixXd&, const estim::RowMatrixXd&, double, const std::optional<estim::RowMatrixXd>&>(),
+             py::arg("n"), py::arg("process_noise"), py::arg("measurement_noise"), py::arg("p_coeff"), py::arg("x0"))
+        .def("iterate", &estim::KF::iterate, py::arg("A"), py::arg("b"), py::arg("x0") = std::nullopt,
+             py::return_value_policy::reference_internal);
+
+    py::class_<estim::KF_Robust, estim::Estimator>(m, "KF_Robust")
+        .def(py::init<int, const estim::RowMatrixXd&, const estim::RowMatrixXd&, double, const std::optional<estim::RowMatrixXd>&>(),
+             py::arg("n"), py::arg("process_noise"), py::arg("measurement_noise"), py::arg("p_coeff"), py::arg("x0"))
+        .def("iterate", &estim::KF_Robust::iterate, py::arg("A"), py::arg("b"), py::arg("x0") = std::nullopt,
+             py::return_value_policy::reference_internal);
+
+    py::class_<estim::RK, estim::Estimator>(m, "RK")
+        .def(py::init<int, const std::optional<estim::RowMatrixXd>&>(),
+             py::arg("n"), py::arg("x0"))
+        .def("iterate", &estim::RK::iterate, py::arg("A"), py::arg("b"), py::arg("x0") = std::nullopt,
+             py::return_value_policy::reference_internal);
+
+    py::class_<estim::GRK, estim::Estimator>(m, "GRK")
+        .def(py::init<int, double, const std::optional<estim::RowMatrixXd>&>(),
+             py::arg("n"), py::arg("tolerance"), py::arg("x0"))
+        .def("iterate", &estim::GRK::iterate, py::arg("A"), py::arg("b"), py::arg("x0") = std::nullopt,
+             py::return_value_policy::reference_internal);
+
+    py::class_<estim::TARK, estim::Estimator>(m, "TARK")
+        .def(py::init<int, int, const std::optional<estim::RowMatrixXd>&>(),
+             py::arg("n"), py::arg("burnin_steps"), py::arg("x0"))
+        .def("iterate", &estim::TARK::iterate, py::arg("A"), py::arg("b"), py::arg("x0") = std::nullopt,
+             py::return_value_policy::reference_internal);
+
+    py::class_<estim::TAGRK, estim::Estimator>(m, "TAGRK")
+        .def(py::init<int, int, double, const std::optional<estim::RowMatrixXd>&>(),
+             py::arg("n"), py::arg("burnin_steps"), py::arg("tolerance"), py::arg("x0"))
+        .def("iterate", &estim::TAGRK::iterate, py::arg("A"), py::arg("b"), py::arg("x0") = std::nullopt,
+             py::return_value_policy::reference_internal);
 
 #ifdef VERSION_INFO
     m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO);

cpp/src/estimator_types.h

@@ -0,0 +1,11 @@
+#pragma once
+
+#include <Eigen/Dense>
+
+namespace estim {
+
+// RowMatrixXd is the default type as it is more easily mapped into by by numpy.
+// The various Kaczmarz methods are also generally row-based, so this assists with locality.
+using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
+} // namespace estim