markov cpp
8 changed files+138−88
zia_benchmark/setup.pymodified+2−0View file
@@ -45,6 +45,8 @@ setup(
4545 ],
4646 ext_modules=[
4747 CMakeExtension("zia_benchmark.algorithms.simple_ans.markov_reconstruct_cpp_ext",
48+ sourcedir="src/zia_benchmark/algorithms/simple_ans"),
49+ CMakeExtension("zia_benchmark.algorithms.simple_ans.markov_predict_cpp_ext",
4850 sourcedir="src/zia_benchmark/algorithms/simple_ans")
4951 ],
5052 cmdclass={
zia_benchmark/src/zia_benchmark/algorithms/simple_ans/CMakeLists.txtmodified+18−3View file
@@ -1,5 +1,5 @@
11 cmake_minimum_required(VERSION 3.15)
2-project(markov_reconstruct_cpp)
2+project(markov_cpp)
33
44 set(CMAKE_CXX_STANDARD 14)
55 set(CMAKE_CXX_STANDARD_REQUIRED ON)
@@ -14,8 +14,23 @@ FetchContent_Declare(
1414 GIT_REPOSITORY https://github.com/pybind/pybind11.git
1515 GIT_TAG v2.11.1
1616 )
17-FetchContent_MakeAvailable(pybind11)
1817
18+# Fetch and include Eigen
19+FetchContent_Declare(
20+ eigen
21+ GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
22+ GIT_TAG 3.4.0
23+)
24+
25+FetchContent_MakeAvailable(pybind11 eigen)
26+
27+# Build markov_reconstruct module
1928 pybind11_add_module(markov_reconstruct_cpp_ext markov_reconstruct.cpp)
2029
21-install(TARGETS markov_reconstruct_cpp_ext DESTINATION zia_benchmark/algorithms/simple_ans)
30+# Build markov_predict module
31+pybind11_add_module(markov_predict_cpp_ext markov_predict.cpp)
32+target_link_libraries(markov_predict_cpp_ext PRIVATE Eigen3::Eigen)
33+
34+# Install both modules
35+install(TARGETS markov_reconstruct_cpp_ext markov_predict_cpp_ext
36+ DESTINATION zia_benchmark/algorithms/simple_ans)
zia_benchmark/src/zia_benchmark/algorithms/simple_ans/__init__.pymodified+5−3View file
@@ -1,8 +1,10 @@
11 import numpy as np
2-from zia_benchmark.algorithms.simple_ans.markov_reconstruct_wrapper import (
2+from zia_benchmark.algorithms.simple_ans.markov_reconstruct import (
33 markov_reconstruct as markov_reconstruct_cpp,
44 )
5-from .markov import markov_predict
5+from zia_benchmark.algorithms.simple_ans.markov_predict import (
6+ markov_predict as markov_predict_cpp,
7+)
68
79 SOURCE_FILE = "simple_ans/__init__.py"
810
@@ -152,7 +154,7 @@ def simple_ans_markov_encode(x: np.ndarray) -> bytes:
152154 from simple_ans import ans_encode
153155
154156 assert x.ndim == 1
155- coeffs, initial, resid = markov_predict(x, M=10)
157+ coeffs, initial, resid = markov_predict_cpp(x, M=6, num_training_samples=10000)
156158 # Encode just the differences
157159 encoded = ans_encode(resid)
158160 if x.dtype == np.uint8:
zia_benchmark/src/zia_benchmark/algorithms/simple_ans/markov.pydeleted+0−66View file
@@ -1,66 +0,0 @@
1-import numpy as np
2-from zia_benchmark._analysis import linear_fit
3-
4-
5-def markov_predict(x: np.ndarray, M: int) -> tuple:
6- """Predict signal using Markov model and return coefficients, initial values and residuals.
7-
8- Args:
9- x: Input signal
10- M: Number of previous samples to use for prediction (default: 20)
11-
12- Returns:
13- tuple: (coefficients, initial_values, residuals)
14- """
15- # Keep initial values for reconstruction
16- initial = x[: M - 1]
17-
18- # Create sequences of M consecutive samples
19- sequences = np.array([x[i : i + M] for i in range(len(x) - M + 1)])
20- predictors = sequences[:, : M - 1] # Use M-1 previous samples to predict
21- target = sequences[:, M - 1] # The value to predict
22-
23- # Get coefficients and prediction function using linear regression
24- coeffs, predict = linear_fit(predictors, target)
25-
26- # Make predictions using the linear model
27- predictions = predict(predictors)
28- predictions = np.round(predictions)
29-
30- # Calculate residuals (difference between actual and predicted values)
31- residuals = target - predictions
32- residuals = residuals.astype(x.dtype)
33-
34- return coeffs, initial, residuals
35-
36-
37-# C++/pybind11 implementation is a lot faster than the Numba implementation
38-# @numba.jit(nopython=True)
39-# def markov_reconstruct(
40-# coeffs: np.ndarray, initial: np.ndarray, resid: np.ndarray
41-# ) -> np.ndarray:
42-# """Reconstruct signal from Markov model parameters and residuals.
43-
44-# Args:
45-# coeffs: Model coefficients from linear regression
46-# initial: Initial values needed for prediction
47-# resid: Prediction residuals
48-
49-# Returns:
50-# np.ndarray: Reconstructed signal
51-# """
52-# M = len(initial) + 1 # Number of samples used in prediction
53-# output = np.zeros(len(resid) + len(initial), dtype=resid.dtype)
54-# output[: len(initial)] = initial # Set initial values
55-
56-# # Reconstruct signal iteratively
57-# for i in range(len(resid)):
58-# # Get previous M-1 values to make prediction
59-# prev_values = output[i : i + M - 1]
60-# # Make prediction using coefficients
61-# prediction = np.sum(coeffs[:-1] * prev_values) + coeffs[-1]
62-# prediction = np.round(prediction)
63-# # Add residual to get actual value
64-# output[i + M - 1] = prediction + resid[i]
65-
66-# return output
zia_benchmark/src/zia_benchmark/algorithms/simple_ans/markov_predict.cppadded+95−0View file
@@ -0,0 +1,95 @@
1+#include <pybind11/pybind11.h>
2+#include <pybind11/numpy.h>
3+#include <pybind11/eigen.h>
4+#include <Eigen/Dense>
5+#include <cmath>
6+#include <iostream>
7+
8+namespace py = pybind11;
9+
10+std::tuple<py::array_t<float>, py::array_t<int16_t>, py::array_t<int16_t>> markov_predict_cpp(
11+ py::array_t<int16_t> x,
12+ size_t M,
13+ size_t num_training_samples
14+) {
15+ // Get array buffer
16+ auto x_buf = x.request();
17+ int16_t* x_ptr = static_cast<int16_t*>(x_buf.ptr);
18+ size_t N = x_buf.shape[0];
19+
20+ // Keep initial values for reconstruction
21+ size_t initial_size = M - 1;
22+ std::vector<ssize_t> initial_shape = {static_cast<ssize_t>(initial_size)};
23+ py::array_t<int16_t> initial(initial_shape);
24+ py::buffer_info initial_buf = initial.request(true);
25+ int16_t* initial_ptr = static_cast<int16_t*>(initial_buf.ptr);
26+
27+ // Copy initial values
28+ for (size_t i = 0; i < initial_size; i++) {
29+ initial_ptr[i] = x_ptr[i];
30+ }
31+
32+ // Create sequences matrix for linear regression
33+ size_t resid_size = N - M + 1;
34+ // Use only num_training_samples sequences for model fitting
35+ size_t num_samples_for_fit = std::min(resid_size, num_training_samples);
36+ // Take evenly spaced samples for training
37+ size_t stride = resid_size > num_training_samples ? resid_size / num_training_samples : 1;
38+
39+ Eigen::MatrixXf predictors(num_samples_for_fit, M - 1);
40+ Eigen::VectorXf target(num_samples_for_fit);
41+
42+ // Fill predictors matrix and target vector with strided training samples
43+ for (size_t i = 0; i < num_samples_for_fit; i++) {
44+ size_t idx = i * stride;
45+ for (size_t j = 0; j < M - 1; j++) {
46+ predictors(i, j) = static_cast<float>(x_ptr[idx + j]);
47+ }
48+ target(i) = static_cast<float>(x_ptr[idx + M - 1]);
49+ }
50+
51+ // Add constant term column (ones) to predictors
52+ Eigen::MatrixXf X(predictors.rows(), predictors.cols() + 1);
53+ X << predictors, Eigen::VectorXf::Ones(predictors.rows());
54+
55+ // Solve least squares problem: X * coeffs = target
56+ Eigen::VectorXf coeffs = X.colPivHouseholderQr().solve(target);
57+
58+ // Create coefficients array
59+ std::vector<ssize_t> coeffs_shape = {static_cast<ssize_t>(M)};
60+ py::array_t<float> coeffs_array(coeffs_shape);
61+ py::buffer_info coeffs_buf = coeffs_array.request(true);
62+ float* coeffs_ptr = static_cast<float*>(coeffs_buf.ptr);
63+
64+ // Copy coefficients
65+ for (size_t i = 0; i < M - 1; i++) {
66+ coeffs_ptr[i] = coeffs(i);
67+ }
68+ coeffs_ptr[M - 1] = coeffs(M - 1); // bias term
69+
70+ // Calculate residuals
71+ std::vector<ssize_t> resid_shape = {static_cast<ssize_t>(resid_size)};
72+ py::array_t<int16_t> residuals(resid_shape);
73+ py::buffer_info resid_buf = residuals.request(true);
74+ int16_t* resid_ptr = static_cast<int16_t*>(resid_buf.ptr);
75+
76+ for (size_t i = 0; i < resid_size; i++) {
77+ float prediction = 0.0f;
78+ for (size_t j = 0; j < M - 1; j++) {
79+ float term = coeffs_ptr[j] * static_cast<float>(x_ptr[i + j]);
80+ prediction += term;
81+ }
82+ prediction += coeffs_ptr[M - 1]; // bias term
83+ float rounded_prediction = std::round(prediction);
84+ resid_ptr[i] = x_ptr[i + M - 1] - static_cast<int16_t>(rounded_prediction);
85+ }
86+
87+ return std::make_tuple(coeffs_array, initial, residuals);
88+}
89+
90+PYBIND11_MODULE(markov_predict_cpp_ext, m) {
91+ m.doc() = "C++ implementation of markov_predict using pybind11";
92+ m.def("markov_predict_cpp", &markov_predict_cpp,
93+ "Predict signal using Markov model and return coefficients, initial values and residuals",
94+ py::arg("x"), py::arg("M"), py::arg("num_training_samples") = 10000);
95+}
zia_benchmark/src/zia_benchmark/algorithms/simple_ans/markov_predict.pyadded+18−0View file
@@ -0,0 +1,18 @@
1+import numpy as np
2+from .markov_predict_cpp_ext import markov_predict_cpp
3+
4+
5+def markov_predict(x: np.ndarray, M: int, num_training_samples: int = 10000) -> tuple:
6+ """Predict signal using Markov model and return coefficients, initial values and residuals using C++ implementation.
7+
8+ Args:
9+ x: Input signal (will be converted to int16)
10+ M: Number of previous samples to use for prediction
11+ num_training_samples: Maximum number of samples to use for fitting the model coefficients (default: 10000).
12+ Using fewer samples speeds up model fitting on large inputs while maintaining accuracy.
13+
14+ Returns:
15+ tuple: (coefficients (float32), initial_values (int16), residuals (int16))
16+ """
17+ # Call C++ implementation with proper type conversion
18+ return markov_predict_cpp(x.astype(np.int16), M, num_training_samples)
zia_benchmark/src/zia_benchmark/algorithms/simple_ans/markov_reconstruct_wrapper.py →zia_benchmark/src/zia_benchmark/algorithms/simple_ans/markov_reconstruct.pyrenamed+0−0View file
No changes to the file's content.
zia_benchmark/test_markov_reconstruct.pydeleted+0−16View file
@@ -1,16 +0,0 @@
1-import numpy as np
2-from zia_benchmark.algorithms.simple_ans.markov import markov_predict, markov_reconstruct
3-from zia_benchmark.algorithms.simple_ans.markov_reconstruct_wrapper import markov_reconstruct as markov_reconstruct_cpp
4-
5-coeffs = np.array([1, 2, 3], dtype=np.float32)
6-initial = np.array([7, 5], dtype=np.int16)
7-resid = np.array([6, 7], dtype=np.int16)
8-a = markov_reconstruct(coeffs, initial, resid)
9-
10-b = markov_reconstruct_cpp(coeffs, initial, resid)
11-
12-print(a)
13-print(b)
14-print(type(b))
15-print(b.shape)
16-print(b.dtype)