markov bench
5 changed files+500−2
.gitignoremodified+4−2View file
@@ -1,2 +1,4 @@
1-benchmark_resuts
2-.benchmark_cache
1+*.out
2+benchmark_results
3+.benchmark_cache
4+__pycache__
devel/compile_markov_bench_cpp.shadded+3−0View file
@@ -0,0 +1,3 @@
1+# sudo apt-get install libeigen3-dev
2+
3+g++ -std=c++17 -O2 -I /usr/include/eigen3 markov_bench.cpp -o markov_bench.out
devel/format_code.shadded+11−0View file
@@ -0,0 +1,11 @@
1+#!/bin/bash
2+
3+# Format Python code
4+echo "Formatting Python code..."
5+black zia_benchmark/src/zia_benchmark
6+
7+# Format TypeScript/JavaScript code
8+echo "Formatting TypeScript/JavaScript code..."
9+cd web-ui && npm run format
10+
11+echo "Code formatting complete!"
devel/markov_bench.cppadded+219−0View file
@@ -0,0 +1,219 @@
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+ *************************************************/
14+
15+#include <Eigen/Dense>
16+#include <iostream>
17+#include <vector>
18+#include <random>
19+#include <chrono>
20+#include <cmath>
21+#include <cstdint>
22+
23+//---------------------------------------------------------------------------
24+// 1) Helper for timing
25+//---------------------------------------------------------------------------
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();
31+}
32+
33+//---------------------------------------------------------------------------
34+// 2) Generate random int16_t data
35+// We create N random integers in [-30000, 30000].
36+//---------------------------------------------------------------------------
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;
47+}
48+
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+//---------------------------------------------------------------------------
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+ }
73+ }
74+
75+ return A;
76+}
77+
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+//---------------------------------------------------------------------------
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;
90+}
91+
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+//---------------------------------------------------------------------------
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]);
113+ }
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;
134+}
135+
136+//---------------------------------------------------------------------------
137+// 6) Main benchmark function
138+//---------------------------------------------------------------------------
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;
219+}
devel/markov_bench.pyadded+263−0View file
@@ -0,0 +1,263 @@
1+import numpy as np
2+import time
3+from sklearn.linear_model import LinearRegression
4+
5+########################################################################
6+# 1) Functions to build the design matrix in float32
7+########################################################################
8+
9+def build_design_matrix_naive(x: np.ndarray, M: int):
10+ """
11+ Naive Python approach to build design matrix A for Markov model.
12+ x is int16, but we'll cast to float32 in A.
13+ A has shape (N - M, M + 1).
14+ Row j = [1, x_{j+M-1}, x_{j+M-2}, ..., x_{j+M-M}],
15+ for j in 0..(N-M-1).
16+ """
17+ N = len(x)
18+ A = np.ones((N - M, M + 1), dtype=np.float32) # store as float32
19+ for row_idx in range(N - M):
20+ for k in range(1, M + 1):
21+ # cast x (int16) to float32 on assignment
22+ A[row_idx, k] = x[M + row_idx - k]
23+ return A
24+
25+def build_design_matrix_vectorized(x: np.ndarray, M: int):
26+ """
27+ Vectorized approach (using slicing) to build the same design matrix in float32.
28+ x is int16; cast slices to float32.
29+ """
30+ N = len(x)
31+ rows = N - M
32+ # First column of 1's in float32
33+ A_cols = [np.ones(rows, dtype=np.float32)]
34+ # Next columns are shifted versions of x, cast to float32
35+ for shift in range(1, M + 1):
36+ A_cols.append(x[M - shift : N - shift].astype(np.float32))
37+ A = np.column_stack(A_cols)
38+ return A
39+
40+########################################################################
41+# 2) Functions to solve for coefficients (stored as float32)
42+########################################################################
43+
44+def solve_lstsq_numpy(A: np.ndarray, y: np.ndarray):
45+ """
46+ Solve for coefficients using NumPy's lstsq solver in float32.
47+ Returns a 1D array of shape (M+1,) in float32.
48+ """
49+ # lstsq may do computations in float64 internally, but we can cast back
50+ c, _, _, _ = np.linalg.lstsq(A, y, rcond=None)
51+ return c.astype(np.float32)
52+
53+def solve_lstsq_sklearn(A: np.ndarray, y: np.ndarray):
54+ """
55+ Solve for coefficients using scikit-learn's LinearRegression.
56+ (fit_intercept=False because we already have the 1-column in A).
57+ We store final coefficients as float32.
58+ """
59+ reg = LinearRegression(fit_intercept=False)
60+ reg.fit(A, y)
61+ # sklearn might do double precision internally; cast to float32
62+ return reg.coef_.astype(np.float32)
63+
64+########################################################################
65+# 3) Functions to predict a new sequence
66+########################################################################
67+
68+def predict_naive(coeffs: np.ndarray, seed: np.ndarray, N: int):
69+ """
70+ Naive Python prediction loop with float32 arithmetic.
71+ - coeffs: float32 array of shape (M+1,) => [c0, c1, ..., cM].
72+ - seed: last M known values (int16), but we'll cast to float32 internally.
73+ - N: number of output points to produce.
74+ The final predictions are rounded to int16.
75+ """
76+ M = len(coeffs) - 1
77+ # We'll store predictions in float32, then round to int16 at the end
78+ pred = np.zeros(N, dtype=np.float32)
79+ # Copy seed (cast to float32)
80+ pred[:M] = seed.astype(np.float32)
81+
82+ # Predict forward
83+ for j in range(M, N):
84+ val = coeffs[0]
85+ for k in range(1, M + 1):
86+ val += coeffs[k] * pred[j - k]
87+ pred[j] = val
88+
89+ # Now round to int16
90+ return np.round(pred).astype(np.int16)
91+
92+# Optional: Numba-accelerated version
93+try:
94+ from numba import njit
95+
96+ @njit
97+ def _predict_markov_numba(pred, coeffs, M):
98+ for j in range(M, len(pred)):
99+ val = coeffs[0]
100+ for k in range(1, M + 1):
101+ val += coeffs[k] * pred[j - k]
102+ pred[j] = val
103+
104+ def predict_numba(coeffs: np.ndarray, seed: np.ndarray, N: int):
105+ """
106+ Numba-accelerated version of predict_naive (float32 internal arithmetic).
107+ Rounds final results to int16.
108+ """
109+ M = len(coeffs) - 1
110+ pred = np.zeros(N, dtype=np.float32)
111+ pred[:M] = seed.astype(np.float32)
112+ _predict_markov_numba(pred, coeffs, M)
113+ return np.round(pred).astype(np.int16)
114+
115+except ImportError:
116+ def predict_numba(coeffs: np.ndarray, seed: np.ndarray, N: int):
117+ """
118+ Fallback if Numba is not installed: use naive approach.
119+ """
120+ print("Numba not installed. Falling back to naive prediction.")
121+ return predict_naive(coeffs, seed, N)
122+
123+########################################################################
124+# 4) Main benchmarking function
125+########################################################################
126+
127+def benchmark_markov_model(N=1_000_000, M=5, seed=0):
128+ """
129+ 1) Generate random data x of length N as int16.
130+ 2) Fit the Markov model in multiple ways:
131+ - Naive (build + solve) in float32
132+ - Vectorized+NumPy (float32)
133+ - Vectorized+sklearn (float32)
134+ Print total fitting time for each method, plus throughput in MB/s
135+ (using 2*N bytes as the data size).
136+ 3) Predict a new sequence in multiple ways (naive vs numba-accelerated),
137+ using float32 arithmetic internally, then rounding to int16.
138+ Print total prediction time + throughput in MB/s.
139+ """
140+ np.random.seed(seed)
141+ # Generate random integers in range [-32768, 32767], but let's just do typical range
142+ x = np.random.randint(-30000, 30000, size=N, dtype=np.int16)
143+
144+ # We define "bytes processed" as 2*N because each x_i is an int16 (2 bytes).
145+ bytes_processed = 2 * N
146+ # Build target array (still in float32 for the solver)
147+ # y corresponds to x[M:], but we cast to float32
148+ y = x[M:].astype(np.float32)
149+
150+ print(f"\n=== Benchmarking Markov Model (N={N}, M={M}, data=int16, coeffs=float32) ===")
151+
152+ ########################################################################
153+ # Fitting: Naive approach
154+ ########################################################################
155+ fit_start = time.perf_counter()
156+
157+ # a) build design matrix (naive)
158+ t0 = time.perf_counter()
159+ A_naive = build_design_matrix_naive(x, M)
160+ dt_build_naive = time.perf_counter() - t0
161+
162+ # b) solve
163+ t0 = time.perf_counter()
164+ c_naive = solve_lstsq_numpy(A_naive, y)
165+ dt_solve_naive = time.perf_counter() - t0
166+
167+ total_naive = time.perf_counter() - fit_start
168+
169+ throughput_build_naive = bytes_processed / (dt_build_naive * 1e6)
170+ throughput_solve_naive = bytes_processed / (dt_solve_naive * 1e6)
171+ throughput_total_naive = bytes_processed / (total_naive * 1e6)
172+
173+ print("\n--- Fitting (Naive) ---")
174+ print(f" Build matrix: {dt_build_naive:.3f}s, ~{throughput_build_naive:.2f} MB/s")
175+ print(f" Solve: {dt_solve_naive:.3f}s, ~{throughput_solve_naive:.2f} MB/s")
176+ print(f" TOTAL: {total_naive:.3f}s, ~{throughput_total_naive:.2f} MB/s")
177+
178+ ########################################################################
179+ # Fitting: Vectorized + NumPy
180+ ########################################################################
181+ fit_start = time.perf_counter()
182+
183+ # a) build design matrix (vectorized)
184+ t0 = time.perf_counter()
185+ A_vec = build_design_matrix_vectorized(x, M)
186+ dt_build_vec = time.perf_counter() - t0
187+
188+ # b) solve
189+ t0 = time.perf_counter()
190+ c_vec = solve_lstsq_numpy(A_vec, y)
191+ dt_solve_vec = time.perf_counter() - t0
192+
193+ total_vec = time.perf_counter() - fit_start
194+
195+ throughput_build_vec = bytes_processed / (dt_build_vec * 1e6)
196+ throughput_solve_vec = bytes_processed / (dt_solve_vec * 1e6)
197+ throughput_total_vec = bytes_processed / (total_vec * 1e6)
198+
199+ print("\n--- Fitting (Vectorized + NumPy) ---")
200+ print(f" Build matrix: {dt_build_vec:.3f}s, ~{throughput_build_vec:.2f} MB/s")
201+ print(f" Solve: {dt_solve_vec:.3f}s, ~{throughput_solve_vec:.2f} MB/s")
202+ print(f" TOTAL: {total_vec:.3f}s, ~{throughput_total_vec:.2f} MB/s")
203+
204+ ########################################################################
205+ # Fitting: Vectorized + scikit-learn
206+ ########################################################################
207+ fit_start = time.perf_counter()
208+
209+ # We can re-use A_vec from above, but to measure a true "total time"
210+ # for the vectorized+sklearn path, let's rebuild it.
211+ t0 = time.perf_counter()
212+ A_vec_sklearn = build_design_matrix_vectorized(x, M)
213+ dt_build_vec_sklearn = time.perf_counter() - t0
214+
215+ t0 = time.perf_counter()
216+ c_sklearn = solve_lstsq_sklearn(A_vec_sklearn, y)
217+ dt_solve_sklearn = time.perf_counter() - t0
218+
219+ total_sklearn = time.perf_counter() - fit_start
220+
221+ throughput_build_sklearn = bytes_processed / (dt_build_vec_sklearn * 1e6)
222+ throughput_solve_sklearn = bytes_processed / (dt_solve_sklearn * 1e6)
223+ throughput_total_sklearn = bytes_processed / (total_sklearn * 1e6)
224+
225+ print("\n--- Fitting (Vectorized + Sklearn) ---")
226+ print(f" Build matrix: {dt_build_vec_sklearn:.3f}s, ~{throughput_build_sklearn:.2f} MB/s")
227+ print(f" Solve: {dt_solve_sklearn:.3f}s, ~{throughput_solve_sklearn:.2f} MB/s")
228+ print(f" TOTAL: {total_sklearn:.3f}s, ~{throughput_total_sklearn:.2f} MB/s")
229+
230+ ########################################################################
231+ # Prediction
232+ ########################################################################
233+ # We can pick c_vec from above as "our" coefficients for prediction.
234+ # We'll seed from the last M points of x (int16).
235+ seed_vals = x[-M:]
236+
237+ # 1) Naive Prediction
238+ t0 = time.perf_counter()
239+ pred_naive = predict_naive(c_vec, seed_vals, N)
240+ dt_pred_naive = time.perf_counter() - t0
241+ throughput_pred_naive = bytes_processed / (dt_pred_naive * 1e6)
242+
243+ print("\n--- Prediction (Naive) ---")
244+ print(f" Predict: {dt_pred_naive:.3f}s, ~{throughput_pred_naive:.2f} MB/s")
245+
246+ # 2) Numba-Accelerated Prediction
247+ t0 = time.perf_counter()
248+ pred_numba = predict_numba(c_vec, seed_vals, N)
249+ dt_pred_numba = time.perf_counter() - t0
250+ throughput_pred_numba = bytes_processed / (dt_pred_numba * 1e6)
251+
252+ print("\n--- Prediction (Numba) ---")
253+ print(f" Predict: {dt_pred_numba:.3f}s, ~{throughput_pred_numba:.2f} MB/s")
254+
255+ print("\n=== Done ===\n")
256+
257+########################################################################
258+# Script entry point
259+########################################################################
260+if __name__ == "__main__":
261+ # Example usage
262+ benchmark_markov_model(N=5_000_000, M=5, seed=42)
263+ # Adjust N, M, and seed as needed.