format C++
5 changed files+294−275
devel/check_format.shmodified+11−0View file
@@ -21,5 +21,16 @@ if ! npm run format:check; then
2121 exit 1
2222 fi
2323
24+# Check C++ code formatting
25+echo "Checking C++ code formatting..."
26+cpp_files=$(find "$PROJECT_ROOT" -type f \( -name "*.cpp" -o -name "*.h" -o -name "*.hpp" \))
27+for file in $cpp_files; do
28+ if ! clang-format --dry-run -Werror "$file" > /dev/null 2>&1; then
29+ echo "C++ code is not properly formatted!"
30+ echo "Please run './devel/format_code.sh' to format the code"
31+ exit 1
32+ fi
33+done
34+
2435 echo "All code is properly formatted!"
2536 exit 0
devel/format_code.shmodified+4−0View file
@@ -12,4 +12,8 @@ black "$PROJECT_ROOT/zia_benchmark/src/zia_benchmark"
1212 echo "Formatting TypeScript/JavaScript code..."
1313 cd "$PROJECT_ROOT/web-ui" && npm run format
1414
15+# Format C++ code
16+echo "Formatting C++ code..."
17+find "$PROJECT_ROOT" -type f \( -name "*.cpp" -o -name "*.h" -o -name "*.hpp" \) -exec clang-format -i {} \;
18+
1519 echo "Code formatting complete!"
devel/markov_bench.cppmodified+160−155View file
@@ -13,37 +13,35 @@
1313 *************************************************/
1414
1515 #include <Eigen/Dense>
16-#include <iostream>
17-#include <vector>
18-#include <random>
1916 #include <chrono>
2017 #include <cmath>
2118 #include <cstdint>
19+#include <iostream>
20+#include <random>
21+#include <vector>
2222
2323 //---------------------------------------------------------------------------
2424 // 1) Helper for timing
2525 //---------------------------------------------------------------------------
26-inline double secondsBetween(
27- const std::chrono::high_resolution_clock::time_point& start,
28- const std::chrono::high_resolution_clock::time_point& end)
29-{
30- return std::chrono::duration<double>(end - start).count();
26+inline double
27+secondsBetween(const std::chrono::high_resolution_clock::time_point &start,
28+ const std::chrono::high_resolution_clock::time_point &end) {
29+ return std::chrono::duration<double>(end - start).count();
3130 }
3231
3332 //---------------------------------------------------------------------------
3433 // 2) Generate random int16_t data
3534 // We create N random integers in [-30000, 30000].
3635 //---------------------------------------------------------------------------
37-std::vector<int16_t> generateData(std::size_t N, unsigned seed = 0)
38-{
39- std::mt19937_64 rng(seed);
40- std::uniform_int_distribution<int> dist(-30000, 30000);
41-
42- std::vector<int16_t> x(N);
43- for (std::size_t i = 0; i < N; ++i) {
44- x[i] = static_cast<int16_t>(dist(rng));
45- }
46- return x;
36+std::vector<int16_t> generateData(std::size_t N, unsigned seed = 0) {
37+ std::mt19937_64 rng(seed);
38+ std::uniform_int_distribution<int> dist(-30000, 30000);
39+
40+ std::vector<int16_t> x(N);
41+ for (std::size_t i = 0; i < N; ++i) {
42+ x[i] = static_cast<int16_t>(dist(rng));
43+ }
44+ return x;
4745 }
4846
4947 //---------------------------------------------------------------------------
@@ -53,26 +51,26 @@ std::vector<int16_t> generateData(std::size_t N, unsigned seed = 0)
5351 // [1, x_{j-1}, x_{j-2}, ..., x_{j-M} ]
5452 // We'll store these in A, shape: (N - M) x (M + 1).
5553 //---------------------------------------------------------------------------
56-Eigen::MatrixXf buildDesignMatrixNaive(const std::vector<int16_t>& x, std::size_t M)
57-{
58- const std::size_t N = x.size();
59- const std::size_t rows = N - M;
60- const std::size_t cols = M + 1;
61-
62- Eigen::MatrixXf A(rows, cols); // float32
63-
64- // Fill A
65- for (std::size_t row = 0; row < rows; ++row) {
66- // First column = 1.0
67- A(row, 0) = 1.0f;
68- // Next columns: x_{(row + M) - k}, for k=1..M
69- for (std::size_t k = 1; k <= M; ++k) {
70- int16_t val = x[(row + M) - k];
71- A(row, k) = static_cast<float>(val); // int16 -> float32
72- }
54+Eigen::MatrixXf buildDesignMatrixNaive(const std::vector<int16_t> &x,
55+ std::size_t M) {
56+ const std::size_t N = x.size();
57+ const std::size_t rows = N - M;
58+ const std::size_t cols = M + 1;
59+
60+ Eigen::MatrixXf A(rows, cols); // float32
61+
62+ // Fill A
63+ for (std::size_t row = 0; row < rows; ++row) {
64+ // First column = 1.0
65+ A(row, 0) = 1.0f;
66+ // Next columns: x_{(row + M) - k}, for k=1..M
67+ for (std::size_t k = 1; k <= M; ++k) {
68+ int16_t val = x[(row + M) - k];
69+ A(row, k) = static_cast<float>(val); // int16 -> float32
7370 }
71+ }
7472
75- return A;
73+ return A;
7674 }
7775
7876 //---------------------------------------------------------------------------
@@ -81,16 +79,18 @@ Eigen::MatrixXf buildDesignMatrixNaive(const std::vector<int16_t>& x, std::size_
8179 // c = argmin_c ||A c - y||^2
8280 // We store and return c as Eigen::VectorXf (float32).
8381 //---------------------------------------------------------------------------
84-Eigen::VectorXf solveCoeffsEigen(const Eigen::MatrixXf& A, const Eigen::VectorXf& y)
85-{
86- // SVD-based solve (ComputeThinU|V for full solution in least-squares sense)
87- // Alternatively: A.colPivHouseholderQr().solve(y), etc.
88- Eigen::VectorXf c = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(y);
89- return c;
82+Eigen::VectorXf solveCoeffsEigen(const Eigen::MatrixXf &A,
83+ const Eigen::VectorXf &y) {
84+ // SVD-based solve (ComputeThinU|V for full solution in least-squares sense)
85+ // Alternatively: A.colPivHouseholderQr().solve(y), etc.
86+ Eigen::VectorXf c =
87+ A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(y);
88+ return c;
9089 }
9190
9291 //---------------------------------------------------------------------------
93-// 5) Predict a new sequence using naive float32 arithmetic and round to int16_t.
92+// 5) Predict a new sequence using naive float32 arithmetic and round to
93+// int16_t.
9494 //
9595 // - coeffs: c0..cM (length M+1)
9696 // - seed: last M real data points (int16_t)
@@ -98,122 +98,127 @@ Eigen::VectorXf solveCoeffsEigen(const Eigen::MatrixXf& A, const Eigen::VectorXf
9898 //
9999 // predicted[j] = c0 + c1*predicted[j-1] + ... + cM*predicted[j-M]
100100 //---------------------------------------------------------------------------
101-std::vector<int16_t> predictNaive(
102- const Eigen::VectorXf& coeffs,
103- const std::vector<int16_t>& seed,
104- std::size_t N)
105-{
106- std::size_t M = coeffs.size() - 1; // since c has M+1 entries
107- std::vector<float> predFloat(N, 0.0f); // store intermediate in float32
108- std::vector<int16_t> predInt(N, 0); // final integer output
109-
110- // Initialize first M from seed (cast to float)
111- for (std::size_t i = 0; i < M; ++i) {
112- predFloat[i] = static_cast<float>(seed[i]);
101+std::vector<int16_t> predictNaive(const Eigen::VectorXf &coeffs,
102+ const std::vector<int16_t> &seed,
103+ std::size_t N) {
104+ std::size_t M = coeffs.size() - 1; // since c has M+1 entries
105+ std::vector<float> predFloat(N, 0.0f); // store intermediate in float32
106+ std::vector<int16_t> predInt(N, 0); // final integer output
107+
108+ // Initialize first M from seed (cast to float)
109+ for (std::size_t i = 0; i < M; ++i) {
110+ predFloat[i] = static_cast<float>(seed[i]);
111+ }
112+
113+ // Predict forward
114+ for (std::size_t j = M; j < N; ++j) {
115+ float val = coeffs(0); // c0
116+ for (std::size_t k = 1; k <= M; ++k) {
117+ val += coeffs(k) * predFloat[j - k];
113118 }
114-
115- // Predict forward
116- for (std::size_t j = M; j < N; ++j) {
117- float val = coeffs(0); // c0
118- for (std::size_t k = 1; k <= M; ++k) {
119- val += coeffs(k) * predFloat[j - k];
120- }
121- predFloat[j] = val;
122- }
123-
124- // Round to int16_t
125- for (std::size_t i = 0; i < N; ++i) {
126- float r = std::round(predFloat[i]);
127- // clamp into int16 range if you want to be safe, but ignoring extremes here:
128- if (r > 32767.f) r = 32767.f;
129- if (r < -32768.f) r = -32768.f;
130- predInt[i] = static_cast<int16_t>(r);
131- }
132-
133- return predInt;
119+ predFloat[j] = val;
120+ }
121+
122+ // Round to int16_t
123+ for (std::size_t i = 0; i < N; ++i) {
124+ float r = std::round(predFloat[i]);
125+ // clamp into int16 range if you want to be safe, but ignoring extremes
126+ // here:
127+ if (r > 32767.f)
128+ r = 32767.f;
129+ if (r < -32768.f)
130+ r = -32768.f;
131+ predInt[i] = static_cast<int16_t>(r);
132+ }
133+
134+ return predInt;
134135 }
135136
136137 //---------------------------------------------------------------------------
137138 // 6) Main benchmark function
138139 //---------------------------------------------------------------------------
139-int main()
140-{
141- // Parameters
142- std::size_t N = 5*1000*1000;
143- std::size_t M = 5;
144- unsigned seed = 42;
145-
146- // Generate data
147- auto startAll = std::chrono::high_resolution_clock::now();
148- auto x = generateData(N, seed);
149- auto endAll = std::chrono::high_resolution_clock::now();
150- double dtGen = secondsBetween(startAll, endAll);
151-
152- // We'll define "bytes processed" as 2*N for throughput (since int16_t=2 bytes).
153- double bytesProcessed = double(2 * N);
154-
155- std::cout << "\n=== C++ Markov Model Benchmark ===\n"
156- << "N=" << N << ", M=" << M << ", data=int16_t, coeffs=float32\n\n";
157-
158- std::cout << "Data generation: " << dtGen << " s\n";
159-
160- //----------------------------------------------------------------------
161- // Build design matrix (naive)
162- //----------------------------------------------------------------------
163- auto t0 = std::chrono::high_resolution_clock::now();
164- Eigen::MatrixXf A = buildDesignMatrixNaive(x, M);
165- auto t1 = std::chrono::high_resolution_clock::now();
166- double dtBuild = secondsBetween(t0, t1);
167-
168- // Prepare y (float32) = x[M..N-1]
169- // We'll store it in an Eigen vector
170- std::size_t rows = N - M;
171- Eigen::VectorXf y(rows);
172- for (std::size_t i = 0; i < rows; ++i) {
173- y(i) = static_cast<float>(x[i + M]);
174- }
175-
176- //----------------------------------------------------------------------
177- // Solve for coefficients
178- //----------------------------------------------------------------------
179- auto t2 = std::chrono::high_resolution_clock::now();
180- Eigen::VectorXf coeffs = solveCoeffsEigen(A, y);
181- auto t3 = std::chrono::high_resolution_clock::now();
182- double dtSolve = secondsBetween(t2, t3);
183-
184- double dtFittingTotal = dtBuild + dtSolve;
185-
186- double buildThroughput = bytesProcessed / (dtBuild * 1.0e6);
187- double solveThroughput = bytesProcessed / (dtSolve * 1.0e6);
188- double totalThroughput = bytesProcessed / (dtFittingTotal * 1.0e6);
189-
190- std::cout << "--- Fitting (Naive build + Eigen solve) ---\n";
191- std::cout << " Build matrix: " << dtBuild << " s, ~" << buildThroughput << " MB/s\n";
192- std::cout << " Solve: " << dtSolve << " s, ~" << solveThroughput << " MB/s\n";
193- std::cout << " TOTAL: " << dtFittingTotal << " s, ~" << totalThroughput << " MB/s\n";
194-
195- //----------------------------------------------------------------------
196- // Prediction
197- //----------------------------------------------------------------------
198- // We'll seed from the last M points of x
199- std::vector<int16_t> seedVals(M);
200- for (std::size_t i = 0; i < M; ++i) {
201- seedVals[i] = x[N - M + i];
202- }
203-
204- auto t4 = std::chrono::high_resolution_clock::now();
205- auto predictions = predictNaive(coeffs, seedVals, N);
206- auto t5 = std::chrono::high_resolution_clock::now();
207- double dtPredict = secondsBetween(t4, t5);
208- double predThroughput = bytesProcessed / (dtPredict * 1.0e6);
209-
210- std::cout << "\n--- Prediction (Naive float32 -> round int16) ---\n";
211- std::cout << " Predict: " << dtPredict << " s, ~" << predThroughput << " MB/s\n";
212-
213- //----------------------------------------------------------------------
214- // Done
215- //----------------------------------------------------------------------
216- std::cout << "\n=== Done ===\n\n";
217-
218- return 0;
140+int main() {
141+ // Parameters
142+ std::size_t N = 5 * 1000 * 1000;
143+ std::size_t M = 5;
144+ unsigned seed = 42;
145+
146+ // Generate data
147+ auto startAll = std::chrono::high_resolution_clock::now();
148+ auto x = generateData(N, seed);
149+ auto endAll = std::chrono::high_resolution_clock::now();
150+ double dtGen = secondsBetween(startAll, endAll);
151+
152+ // We'll define "bytes processed" as 2*N for throughput (since int16_t=2
153+ // bytes).
154+ double bytesProcessed = double(2 * N);
155+
156+ std::cout << "\n=== C++ Markov Model Benchmark ===\n"
157+ << "N=" << N << ", M=" << M << ", data=int16_t, coeffs=float32\n\n";
158+
159+ std::cout << "Data generation: " << dtGen << " s\n";
160+
161+ //----------------------------------------------------------------------
162+ // Build design matrix (naive)
163+ //----------------------------------------------------------------------
164+ auto t0 = std::chrono::high_resolution_clock::now();
165+ Eigen::MatrixXf A = buildDesignMatrixNaive(x, M);
166+ auto t1 = std::chrono::high_resolution_clock::now();
167+ double dtBuild = secondsBetween(t0, t1);
168+
169+ // Prepare y (float32) = x[M..N-1]
170+ // We'll store it in an Eigen vector
171+ std::size_t rows = N - M;
172+ Eigen::VectorXf y(rows);
173+ for (std::size_t i = 0; i < rows; ++i) {
174+ y(i) = static_cast<float>(x[i + M]);
175+ }
176+
177+ //----------------------------------------------------------------------
178+ // Solve for coefficients
179+ //----------------------------------------------------------------------
180+ auto t2 = std::chrono::high_resolution_clock::now();
181+ Eigen::VectorXf coeffs = solveCoeffsEigen(A, y);
182+ auto t3 = std::chrono::high_resolution_clock::now();
183+ double dtSolve = secondsBetween(t2, t3);
184+
185+ double dtFittingTotal = dtBuild + dtSolve;
186+
187+ double buildThroughput = bytesProcessed / (dtBuild * 1.0e6);
188+ double solveThroughput = bytesProcessed / (dtSolve * 1.0e6);
189+ double totalThroughput = bytesProcessed / (dtFittingTotal * 1.0e6);
190+
191+ std::cout << "--- Fitting (Naive build + Eigen solve) ---\n";
192+ std::cout << " Build matrix: " << dtBuild << " s, ~" << buildThroughput
193+ << " MB/s\n";
194+ std::cout << " Solve: " << dtSolve << " s, ~" << solveThroughput
195+ << " MB/s\n";
196+ std::cout << " TOTAL: " << dtFittingTotal << " s, ~"
197+ << totalThroughput << " MB/s\n";
198+
199+ //----------------------------------------------------------------------
200+ // Prediction
201+ //----------------------------------------------------------------------
202+ // We'll seed from the last M points of x
203+ std::vector<int16_t> seedVals(M);
204+ for (std::size_t i = 0; i < M; ++i) {
205+ seedVals[i] = x[N - M + i];
206+ }
207+
208+ auto t4 = std::chrono::high_resolution_clock::now();
209+ auto predictions = predictNaive(coeffs, seedVals, N);
210+ auto t5 = std::chrono::high_resolution_clock::now();
211+ double dtPredict = secondsBetween(t4, t5);
212+ double predThroughput = bytesProcessed / (dtPredict * 1.0e6);
213+
214+ std::cout << "\n--- Prediction (Naive float32 -> round int16) ---\n";
215+ std::cout << " Predict: " << dtPredict << " s, ~" << predThroughput
216+ << " MB/s\n";
217+
218+ //----------------------------------------------------------------------
219+ // Done
220+ //----------------------------------------------------------------------
221+ std::cout << "\n=== Done ===\n\n";
222+
223+ return 0;
219224 }
zia_benchmark/src/zia_benchmark/algorithms/simple_ans/markov_predict.cppmodified+71−71View file
@@ -1,95 +1,95 @@
1-#include <pybind11/pybind11.h>
2-#include <pybind11/numpy.h>
3-#include <pybind11/eigen.h>
41 #include <Eigen/Dense>
52 #include <cmath>
63 #include <iostream>
4+#include <pybind11/eigen.h>
5+#include <pybind11/numpy.h>
6+#include <pybind11/pybind11.h>
77
88 namespace py = pybind11;
99
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];
10+std::tuple<py::array_t<float>, py::array_t<int16_t>, py::array_t<int16_t>>
11+markov_predict_cpp(py::array_t<int16_t> x, size_t M,
12+ size_t num_training_samples) {
13+ // Get array buffer
14+ auto x_buf = x.request();
15+ int16_t *x_ptr = static_cast<int16_t *>(x_buf.ptr);
16+ size_t N = x_buf.shape[0];
1917
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);
18+ // Keep initial values for reconstruction
19+ size_t initial_size = M - 1;
20+ std::vector<ssize_t> initial_shape = {static_cast<ssize_t>(initial_size)};
21+ py::array_t<int16_t> initial(initial_shape);
22+ py::buffer_info initial_buf = initial.request(true);
23+ int16_t *initial_ptr = static_cast<int16_t *>(initial_buf.ptr);
2624
27- // Copy initial values
28- for (size_t i = 0; i < initial_size; i++) {
29- initial_ptr[i] = x_ptr[i];
30- }
25+ // Copy initial values
26+ for (size_t i = 0; i < initial_size; i++) {
27+ initial_ptr[i] = x_ptr[i];
28+ }
3129
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;
30+ // Create sequences matrix for linear regression
31+ size_t resid_size = N - M + 1;
32+ // Use only num_training_samples sequences for model fitting
33+ size_t num_samples_for_fit = std::min(resid_size, num_training_samples);
34+ // Take evenly spaced samples for training
35+ size_t stride =
36+ resid_size > num_training_samples ? resid_size / num_training_samples : 1;
3837
39- Eigen::MatrixXf predictors(num_samples_for_fit, M - 1);
40- Eigen::VectorXf target(num_samples_for_fit);
38+ Eigen::MatrixXf predictors(num_samples_for_fit, M - 1);
39+ Eigen::VectorXf target(num_samples_for_fit);
4140
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]);
41+ // Fill predictors matrix and target vector with strided training samples
42+ for (size_t i = 0; i < num_samples_for_fit; i++) {
43+ size_t idx = i * stride;
44+ for (size_t j = 0; j < M - 1; j++) {
45+ predictors(i, j) = static_cast<float>(x_ptr[idx + j]);
4946 }
47+ target(i) = static_cast<float>(x_ptr[idx + M - 1]);
48+ }
5049
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());
50+ // Add constant term column (ones) to predictors
51+ Eigen::MatrixXf X(predictors.rows(), predictors.cols() + 1);
52+ X << predictors, Eigen::VectorXf::Ones(predictors.rows());
5453
55- // Solve least squares problem: X * coeffs = target
56- Eigen::VectorXf coeffs = X.colPivHouseholderQr().solve(target);
54+ // Solve least squares problem: X * coeffs = target
55+ Eigen::VectorXf coeffs = X.colPivHouseholderQr().solve(target);
5756
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);
57+ // Create coefficients array
58+ std::vector<ssize_t> coeffs_shape = {static_cast<ssize_t>(M)};
59+ py::array_t<float> coeffs_array(coeffs_shape);
60+ py::buffer_info coeffs_buf = coeffs_array.request(true);
61+ float *coeffs_ptr = static_cast<float *>(coeffs_buf.ptr);
6362
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
63+ // Copy coefficients
64+ for (size_t i = 0; i < M - 1; i++) {
65+ coeffs_ptr[i] = coeffs(i);
66+ }
67+ coeffs_ptr[M - 1] = coeffs(M - 1); // bias term
6968
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);
69+ // Calculate residuals
70+ std::vector<ssize_t> resid_shape = {static_cast<ssize_t>(resid_size)};
71+ py::array_t<int16_t> residuals(resid_shape);
72+ py::buffer_info resid_buf = residuals.request(true);
73+ int16_t *resid_ptr = static_cast<int16_t *>(resid_buf.ptr);
7574
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);
75+ for (size_t i = 0; i < resid_size; i++) {
76+ float prediction = 0.0f;
77+ for (size_t j = 0; j < M - 1; j++) {
78+ float term = coeffs_ptr[j] * static_cast<float>(x_ptr[i + j]);
79+ prediction += term;
8580 }
81+ prediction += coeffs_ptr[M - 1]; // bias term
82+ float rounded_prediction = std::round(prediction);
83+ resid_ptr[i] = x_ptr[i + M - 1] - static_cast<int16_t>(rounded_prediction);
84+ }
8685
87- return std::make_tuple(coeffs_array, initial, residuals);
86+ return std::make_tuple(coeffs_array, initial, residuals);
8887 }
8988
9089 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);
90+ m.doc() = "C++ implementation of markov_predict using pybind11";
91+ m.def("markov_predict_cpp", &markov_predict_cpp,
92+ "Predict signal using Markov model and return coefficients, initial "
93+ "values and residuals",
94+ py::arg("x"), py::arg("M"), py::arg("num_training_samples") = 10000);
9595 }
zia_benchmark/src/zia_benchmark/algorithms/simple_ans/markov_reconstruct.cppmodified+48−49View file
@@ -1,69 +1,68 @@
1-#include <pybind11/pybind11.h>
2-#include <pybind11/numpy.h>
31 #include <cmath>
42 #include <iostream>
3+#include <pybind11/numpy.h>
4+#include <pybind11/pybind11.h>
55
66 namespace py = pybind11;
77
8-py::array_t<int16_t> markov_reconstruct_cpp(
9- py::array_t<float> coeffs,
10- py::array_t<int16_t> initial,
11- py::array_t<int16_t> resid
12-) {
13- // Get array buffers
14- auto coeffs_buf = coeffs.request();
15- auto initial_buf = initial.request();
16- auto resid_buf = resid.request();
8+py::array_t<int16_t> markov_reconstruct_cpp(py::array_t<float> coeffs,
9+ py::array_t<int16_t> initial,
10+ py::array_t<int16_t> resid) {
11+ // Get array buffers
12+ auto coeffs_buf = coeffs.request();
13+ auto initial_buf = initial.request();
14+ auto resid_buf = resid.request();
1715
18- // Get raw pointers to data
19- float* coeffs_ptr = static_cast<float*>(coeffs_buf.ptr);
20- int16_t* initial_ptr = static_cast<int16_t*>(initial_buf.ptr);
21- int16_t* resid_ptr = static_cast<int16_t*>(resid_buf.ptr);
16+ // Get raw pointers to data
17+ float *coeffs_ptr = static_cast<float *>(coeffs_buf.ptr);
18+ int16_t *initial_ptr = static_cast<int16_t *>(initial_buf.ptr);
19+ int16_t *resid_ptr = static_cast<int16_t *>(resid_buf.ptr);
2220
23- // Calculate dimensions
24- size_t M = initial_buf.shape[0] + 1; // Number of samples used in prediction
25- size_t output_size = resid_buf.shape[0] + initial_buf.shape[0];
21+ // Calculate dimensions
22+ size_t M = initial_buf.shape[0] + 1; // Number of samples used in prediction
23+ size_t output_size = resid_buf.shape[0] + initial_buf.shape[0];
2624
27- // Create output array with explicit shape and memory ownership
28- std::vector<ssize_t> shape = {static_cast<ssize_t>(output_size)};
29- py::array_t<int16_t> output(shape);
30- py::buffer_info output_buf = output.request(true); // Request writable buffer
31- int16_t* output_ptr = static_cast<int16_t*>(output_buf.ptr);
25+ // Create output array with explicit shape and memory ownership
26+ std::vector<ssize_t> shape = {static_cast<ssize_t>(output_size)};
27+ py::array_t<int16_t> output(shape);
28+ py::buffer_info output_buf = output.request(true); // Request writable buffer
29+ int16_t *output_ptr = static_cast<int16_t *>(output_buf.ptr);
3230
33- // Copy initial values with bounds check
34- for (size_t i = 0; i < initial_buf.shape[0] && i < output_size; i++) {
35- output_ptr[i] = initial_ptr[i];
36- }
31+ // Copy initial values with bounds check
32+ for (size_t i = 0; i < initial_buf.shape[0] && i < output_size; i++) {
33+ output_ptr[i] = initial_ptr[i];
34+ }
3735
38- size_t resid_size = resid_buf.shape[0];
36+ size_t resid_size = resid_buf.shape[0];
3937
40- // Reconstruct signal iteratively
41- for (size_t i = 0; i < resid_size; i++) {
42- float prediction = 0.0f;
38+ // Reconstruct signal iteratively
39+ for (size_t i = 0; i < resid_size; i++) {
40+ float prediction = 0.0f;
4341
44- // Calculate prediction using coefficients (excluding bias term)
45- for (size_t j = 0; j < M - 1; j++) {
46- float term = coeffs_ptr[j] * static_cast<float>(output_ptr[i + j]);
47- prediction += term;
48- }
42+ // Calculate prediction using coefficients (excluding bias term)
43+ for (size_t j = 0; j < M - 1; j++) {
44+ float term = coeffs_ptr[j] * static_cast<float>(output_ptr[i + j]);
45+ prediction += term;
46+ }
4947
50- // Add bias term separately
51- prediction += coeffs_ptr[coeffs_buf.shape[0] - 1];
48+ // Add bias term separately
49+ prediction += coeffs_ptr[coeffs_buf.shape[0] - 1];
5250
53- // Round prediction to nearest integer
54- float rounded_prediction = std::round(prediction);
51+ // Round prediction to nearest integer
52+ float rounded_prediction = std::round(prediction);
5553
56- // Add residual and store result
57- int16_t final_value = static_cast<int16_t>(rounded_prediction + static_cast<float>(resid_ptr[i]));
58- output_ptr[i + M - 1] = final_value;
59- }
54+ // Add residual and store result
55+ int16_t final_value = static_cast<int16_t>(
56+ rounded_prediction + static_cast<float>(resid_ptr[i]));
57+ output_ptr[i + M - 1] = final_value;
58+ }
6059
61- return output;
60+ return output;
6261 }
6362
6463 PYBIND11_MODULE(markov_reconstruct_cpp_ext, m) {
65- m.doc() = "C++ implementation of markov_reconstruct using pybind11";
66- m.def("markov_reconstruct_cpp", &markov_reconstruct_cpp,
67- "Reconstruct signal from Markov model parameters and residuals",
68- py::arg("coeffs"), py::arg("initial"), py::arg("resid"));
64+ m.doc() = "C++ implementation of markov_reconstruct using pybind11";
65+ m.def("markov_reconstruct_cpp", &markov_reconstruct_cpp,
66+ "Reconstruct signal from Markov model parameters and residuals",
67+ py::arg("coeffs"), py::arg("initial"), py::arg("resid"));
6968 }