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 *************************************************/
67290a2markov benchJeremy Magland 22
23//---------------------------------------------------------------------------
24// 1) Helper for timing
25//---------------------------------------------------------------------------
27secondsBetween(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();
32//---------------------------------------------------------------------------
33// 2) Generate random int16_t data
34// We create N random integers in [-30000, 30000].
35//---------------------------------------------------------------------------
03f201aformat C++Jeremy Magland 36std::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);
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;
47//---------------------------------------------------------------------------
48// 3) Build the design matrix A (float32) in a naive manner.
49//
50// For a Markov model of order M, we want rows for j in [M..N-1]:
51// [1, x_{j-1}, x_{j-2}, ..., x_{j-M} ]
52// We'll store these in A, shape: (N - M) x (M + 1).
53//---------------------------------------------------------------------------
03f201aformat C++Jeremy Magland 54Eigen::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;
60 Eigen::MatrixXf A(rows, cols); // float32
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
67290a2markov benchJeremy Magland 72
76//---------------------------------------------------------------------------
77// 4) Solve for coefficients c in float32 using Eigen's SVD.
78//
79// c = argmin_c ||A c - y||^2
80// We store and return c as Eigen::VectorXf (float32).
81//---------------------------------------------------------------------------
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;
91//---------------------------------------------------------------------------
03f201aformat C++Jeremy Magland 92// 5) Predict a new sequence using naive float32 arithmetic and round to
93// int16_t.
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//---------------------------------------------------------------------------
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
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 }
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];
120 }
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 }
134 return predInt;
137//---------------------------------------------------------------------------
138// 6) Main benchmark function
139//---------------------------------------------------------------------------
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
153 // bytes).
154 double bytesProcessed = double(2 * N);
156 std::cout << "\n=== C++ Markov Model Benchmark ===\n"
157 << "N=" << N << ", M=" << M << ", data=int16_t, coeffs=float32\n\n";
159 std::cout << "Data generation: " << dtGen << " s\n";
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);
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 }
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);
185 double dtFittingTotal = dtBuild + dtSolve;
187 double buildThroughput = bytesProcessed / (dtBuild * 1.0e6);
188 double solveThroughput = bytesProcessed / (dtSolve * 1.0e6);
189 double totalThroughput = bytesProcessed / (dtFittingTotal * 1.0e6);
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";
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 }
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);
214 std::cout << "\n--- Prediction (Naive float32 -> round int16) ---\n";
215 std::cout << " Predict: " << dtPredict << " s, ~" << predThroughput
216 << " MB/s\n";
218 //----------------------------------------------------------------------
219 // Done
220 //----------------------------------------------------------------------
221 std::cout << "\n=== Done ===\n\n";
223 return 0;