8namespace py = pybind11;
10template <typename T>
11py::array_t<T> markov_reconstruct_impl(py::array_t<float> coeffs,
12 py::array_t<T> initial,
13 py::array_t<T> resid) {
14 // Get array buffers
15 auto coeffs_buf = coeffs.request();
16 auto initial_buf = initial.request();
17 auto resid_buf = resid.request();
19 // Get raw pointers to data
20 float *coeffs_ptr = static_cast<float *>(coeffs_buf.ptr);
21 T *initial_ptr = static_cast<T *>(initial_buf.ptr);
22 T *resid_ptr = static_cast<T *>(resid_buf.ptr);
24 // Calculate dimensions
25 size_t M = initial_buf.shape[0] + 1; // Number of samples used in prediction
26 size_t output_size = resid_buf.shape[0] + initial_buf.shape[0];
28 // Create output array with explicit shape and memory ownership
29 std::vector<ssize_t> shape = {static_cast<ssize_t>(output_size)};
30 py::array_t<T> output(shape);
31 py::buffer_info output_buf = output.request(true); // Request writable buffer
32 T *output_ptr = static_cast<T *>(output_buf.ptr);
34 // Copy initial values with bounds check
35 for (size_t i = 0; i < initial_buf.shape[0] && i < output_size; i++) {
36 output_ptr[i] = initial_ptr[i];
37 }
39 size_t resid_size = resid_buf.shape[0];
41 // Reconstruct signal iteratively
42 for (size_t i = 0; i < resid_size; i++) {
43 float prediction = 0.0f;
45 // Calculate prediction using coefficients (excluding bias term)
46 for (size_t j = 0; j < M - 1; j++) {
47 float term = coeffs_ptr[j] * static_cast<float>(output_ptr[i + j]);
48 prediction += term;
49 }
51 // Add bias term separately
52 prediction += coeffs_ptr[coeffs_buf.shape[0] - 1];
54 // Round prediction to nearest integer
55 float rounded_prediction = std::round(prediction);
57 // Add residual and store result
58 T final_value =
59 static_cast<T>(rounded_prediction + static_cast<float>(resid_ptr[i]));
60 output_ptr[i + M - 1] = final_value;
61 }
63 return output;
64}