1/*************************************************
2 * markov_bench.cpp
3 *
4 * A C++ program demonstrating:
5 * 1. Generating int16_t data (size N).
6 * 2. Building a design matrix in float32.
7 * 3. Solving for coefficients in float32 with Eigen.
8 * 4. Predicting a new sequence (float32 -> round to int16_t).
9 * 5. Measuring times & throughput (MB/s).
10 *
11 * Requires Eigen for linear algebra:
12 * https://eigen.tuxfamily.org/
13 *************************************************/
23//---------------------------------------------------------------------------
24// 1) Helper for timing
25//---------------------------------------------------------------------------
26inline 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();
31}
33//---------------------------------------------------------------------------
34// 2) Generate random int16_t data
35// We create N random integers in [-30000, 30000].
36//---------------------------------------------------------------------------
37std::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);
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;
47}
49//---------------------------------------------------------------------------
50// 3) Build the design matrix A (float32) in a naive manner.
51//
52// For a Markov model of order M, we want rows for j in [M..N-1]:
53// [1, x_{j-1}, x_{j-2}, ..., x_{j-M} ]
54// We'll store these in A, shape: (N - M) x (M + 1).
55//---------------------------------------------------------------------------
56Eigen::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;
62 Eigen::MatrixXf A(rows, cols); // float32
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 }
73 }
75 return A;
76}
78//---------------------------------------------------------------------------
79// 4) Solve for coefficients c in float32 using Eigen's SVD.
80//
81// c = argmin_c ||A c - y||^2
82// We store and return c as Eigen::VectorXf (float32).
83//---------------------------------------------------------------------------
84Eigen::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;
90}
92//---------------------------------------------------------------------------
93// 5) Predict a new sequence using naive float32 arithmetic and round to int16_t.
94//
95// - coeffs: c0..cM (length M+1)
96// - seed: last M real data points (int16_t)
97// - Output: N predicted points, as int16_t.
98//
99// predicted[j] = c0 + c1*predicted[j-1] + ... + cM*predicted[j-M]
100//---------------------------------------------------------------------------
101std::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
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]);
113 }
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 }
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 }
133 return predInt;
134}
136//---------------------------------------------------------------------------
137// 6) Main benchmark function
138//---------------------------------------------------------------------------
139int main()
140{
141 // Parameters
142 std::size_t N = 5*1000*1000;
143 std::size_t M = 5;
144 unsigned seed = 42;
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);
152 // We'll define "bytes processed" as 2*N for throughput (since int16_t=2 bytes).
153 double bytesProcessed = double(2 * N);
155 std::cout << "\n=== C++ Markov Model Benchmark ===\n"
156 << "N=" << N << ", M=" << M << ", data=int16_t, coeffs=float32\n\n";
158 std::cout << "Data generation: " << dtGen << " s\n";
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);
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 }
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);
184 double dtFittingTotal = dtBuild + dtSolve;
186 double buildThroughput = bytesProcessed / (dtBuild * 1.0e6);
187 double solveThroughput = bytesProcessed / (dtSolve * 1.0e6);
188 double totalThroughput = bytesProcessed / (dtFittingTotal * 1.0e6);
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";
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 }
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);
210 std::cout << "\n--- Prediction (Naive float32 -> round int16) ---\n";
211 std::cout << " Predict: " << dtPredict << " s, ~" << predThroughput << " MB/s\n";
213 //----------------------------------------------------------------------
214 // Done
215 //----------------------------------------------------------------------
216 std::cout << "\n=== Done ===\n\n";
218 return 0;
219}