10namespace py = pybind11;
12template <typename T>
13std::tuple<py::array_t<float>, py::array_t<T>, py::array_t<T>>
14lpc_predict_impl(py::array_t<T> x, size_t M, size_t num_training_samples) {
15 // Get array buffer
16 auto x_buf = x.request();
17 T *x_ptr = static_cast<T *>(x_buf.ptr);
18 size_t N = x_buf.shape[0];
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<T> initial(initial_shape);
24 py::buffer_info initial_buf = initial.request(true);
25 T *initial_ptr = static_cast<T *>(initial_buf.ptr);
27 // Copy initial values
28 for (size_t i = 0; i < initial_size; i++) {
29 initial_ptr[i] = x_ptr[i];
30 }
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 =
38 resid_size > num_training_samples ? resid_size / num_training_samples : 1;
40 Eigen::MatrixXf predictors(num_samples_for_fit, M - 1);
41 Eigen::VectorXf target(num_samples_for_fit);
43 // Fill predictors matrix and target vector with strided training samples
44 for (size_t i = 0; i < num_samples_for_fit; i++) {
45 size_t idx = i * stride;
46 for (size_t j = 0; j < M - 1; j++) {
47 predictors(i, j) = static_cast<float>(x_ptr[idx + j]);
48 }
49 target(i) = static_cast<float>(x_ptr[idx + M - 1]);
50 }
52 // Add constant term column (ones) to predictors
53 Eigen::MatrixXf X(predictors.rows(), predictors.cols() + 1);
54 X << predictors, Eigen::VectorXf::Ones(predictors.rows());
56 // Solve least squares problem: X * coeffs = target
57 Eigen::VectorXf coeffs = X.colPivHouseholderQr().solve(target);
59 // Create coefficients array
60 std::vector<ssize_t> coeffs_shape = {static_cast<ssize_t>(M)};
61 py::array_t<float> coeffs_array(coeffs_shape);
62 py::buffer_info coeffs_buf = coeffs_array.request(true);
63 float *coeffs_ptr = static_cast<float *>(coeffs_buf.ptr);
65 // Copy coefficients
66 for (size_t i = 0; i < M - 1; i++) {
67 coeffs_ptr[i] = coeffs(i);
68 }
69 coeffs_ptr[M - 1] = coeffs(M - 1); // bias term
71 // Calculate residuals
72 std::vector<ssize_t> resid_shape = {static_cast<ssize_t>(resid_size)};
73 py::array_t<T> residuals(resid_shape);
74 py::buffer_info resid_buf = residuals.request(true);
75 T *resid_ptr = static_cast<T *>(resid_buf.ptr);
77 for (size_t i = 0; i < resid_size; i++) {
78 float prediction = 0.0f;
79 for (size_t j = 0; j < M - 1; j++) {
80 float term = coeffs_ptr[j] * static_cast<float>(x_ptr[i + j]);
81 prediction += term;
82 }
83 prediction += coeffs_ptr[M - 1]; // bias term
84 float rounded_prediction = std::round(prediction);
85 resid_ptr[i] = x_ptr[i + M - 1] - static_cast<T>(rounded_prediction);
86 }
88 return std::make_tuple(coeffs_array, initial, residuals);
89}