Skip to content

Commit a1a4ec3

Browse files
HanSur94HanSur94
authored
Merge pull request #38 from HanSur94/claude/optimize-sensor-resolve-1Xpjc
Claude/optimize sensor resolve 1 xpjc
2 parents 1a29d6f + d3356a5 commit a1a4ec3

9 files changed

Lines changed: 1277 additions & 101 deletions

File tree

‎benchmarks/benchmark_resolve_stress.m‎

Lines changed: 284 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,284 @@
1+
function benchmark_resolve_stress()
2+
%BENCHMARK_RESOLVE_STRESS Stress test for Sensor.resolve() at 500M points.
3+
% Creates a realistic sensor with 500M datapoints and 4 thresholds,
4+
% each governed by different condition rules across 2 state channels
5+
% that change frequently (~5000 transitions each).
6+
%
7+
% This exercises the full resolve pipeline:
8+
% - Segment boundary collection with ~10K unique boundaries
9+
% - Composite state evaluation at each boundary
10+
% - Condition matching & segment index mapping
11+
% - Batch violation detection (MEX or vectorized)
12+
% - Threshold merge & step-function conversion
13+
%
14+
% Run:
15+
% >> benchmark_resolve_stress
16+
17+
repo_root = fileparts(fileparts(mfilename('fullpath')));
18+
addpath(repo_root);
19+
install();
20+
21+
N = 500e6;
22+
23+
fprintf('\n');
24+
fprintf('============================================================\n');
25+
fprintf(' Sensor.resolve() Stress Test — %.0fM datapoints\n', N / 1e6);
26+
fprintf('============================================================\n');
27+
fprintf(' Building sensor data... ');
28+
29+
tBuild = tic;
30+
31+
% ------------------------------------------------------------------
32+
% 1. Sensor data: 500M points of realistic process signal
33+
% Base signal with multiple frequency components + noise + drift
34+
% ------------------------------------------------------------------
35+
x = linspace(0, 10000, N);
36+
37+
% Layered signal: slow drift + medium oscillation + fast noise
38+
y = 50 ...
39+
+ 10 * sin(2*pi*x / 3000) ... % slow thermal drift
40+
+ 8 * sin(2*pi*x / 200) ... % process oscillation
41+
+ 5 * randn(1, N); % measurement noise
42+
43+
% Inject deliberate excursions to guarantee violations in every regime
44+
% High spikes in running states
45+
y(round(N*0.05):round(N*0.05)+2000) = 95 + 5*randn(1, 2001);
46+
y(round(N*0.35):round(N*0.35)+3000) = 90 + 5*randn(1, 3001);
47+
y(round(N*0.65):round(N*0.65)+1500) = 85 + 3*randn(1, 1501);
48+
% Low dips in evacuated states
49+
y(round(N*0.15):round(N*0.15)+2500) = 10 + 3*randn(1, 2501);
50+
y(round(N*0.55):round(N*0.55)+1000) = 5 + 2*randn(1, 1001);
51+
y(round(N*0.85):round(N*0.85)+2000) = 15 + 4*randn(1, 2001);
52+
53+
fprintf('%.1f s\n', toc(tBuild));
54+
55+
% ------------------------------------------------------------------
56+
% 2. State channels with many transitions (~5000 each)
57+
% This creates ~10K segment boundaries — the core stress factor
58+
% for the resolve pipeline beyond raw data size.
59+
% ------------------------------------------------------------------
60+
fprintf(' Building state channels... ');
61+
tState = tic;
62+
63+
% Machine state: cycles through 0 (idle) → 1 (running) → 2 (evacuated)
64+
% with irregular interval lengths to create realistic patterns
65+
nMachineTransitions = 5000;
66+
rng(42);
67+
machineIntervals = 1.5 + 1.5 * rand(1, nMachineTransitions);
68+
machineX = [0, cumsum(machineIntervals)];
69+
% Scale to fit data range
70+
machineX = machineX / machineX(end) * x(end);
71+
% Cycle through states 0→1→2→1→0→1→2→...
72+
statePattern = [0 1 2 1];
73+
machineY = zeros(1, numel(machineX));
74+
for i = 1:numel(machineX)
75+
machineY(i) = statePattern(mod(i-1, numel(statePattern)) + 1);
76+
end
77+
78+
scMachine = StateChannel('machine');
79+
scMachine.X = machineX;
80+
scMachine.Y = machineY;
81+
82+
% Operating zone: string-valued, toggles between A/B/C
83+
% Different transition frequency than machine state → many unique combos
84+
nZoneTransitions = 4000;
85+
zoneIntervals = 2.0 + 2.0 * rand(1, nZoneTransitions);
86+
zoneX = [0, cumsum(zoneIntervals)];
87+
zoneX = zoneX / zoneX(end) * x(end);
88+
zonePattern = {'A', 'B', 'C'};
89+
zoneY = cell(1, numel(zoneX));
90+
for i = 1:numel(zoneX)
91+
zoneY{i} = zonePattern{mod(i-1, numel(zonePattern)) + 1};
92+
end
93+
94+
scZone = StateChannel('zone');
95+
scZone.X = zoneX;
96+
scZone.Y = zoneY;
97+
98+
fprintf('%.1f s\n', toc(tState));
99+
fprintf(' Machine transitions: %d\n', numel(machineX));
100+
fprintf(' Zone transitions: %d\n', numel(zoneX));
101+
102+
% ------------------------------------------------------------------
103+
% 3. Build sensor with 4 threshold rules (different conditions each)
104+
% ------------------------------------------------------------------
105+
fprintf(' Configuring sensor + thresholds...\n');
106+
107+
s = Sensor('stress_test', 'Name', 'Process Stress Test', 'Units', 'bar');
108+
s.X = x;
109+
s.Y = y;
110+
s.addStateChannel(scMachine);
111+
s.addStateChannel(scZone);
112+
113+
% Rule 1: Upper alarm when running (machine==1), any zone
114+
% → Active in ~25% of timeline, single-condition
115+
s.addThresholdRule(struct('machine', 1), 75, ...
116+
'Direction', 'upper', 'Label', 'HH (running)', ...
117+
'Color', [0.9 0.1 0.1], 'LineStyle', '--');
118+
119+
% Rule 2: Lower alarm when evacuated (machine==2), any zone
120+
% → Active in ~25% of timeline, single-condition
121+
s.addThresholdRule(struct('machine', 2), 25, ...
122+
'Direction', 'lower', 'Label', 'LL (evacuated)', ...
123+
'Color', [0.1 0.1 0.9], 'LineStyle', '--');
124+
125+
% Rule 3: Upper alarm when running AND zone B (two conditions)
126+
% → Active in ~8% of timeline (intersection), stricter limit
127+
s.addThresholdRule(struct('machine', 1, 'zone', 'B'), 65, ...
128+
'Direction', 'upper', 'Label', 'HH (running+B)', ...
129+
'Color', [0.8 0.0 0.8], 'LineStyle', ':');
130+
131+
% Rule 4: Lower alarm when idle AND zone C (two conditions)
132+
% → Active in ~8% of timeline (intersection)
133+
s.addThresholdRule(struct('machine', 0, 'zone', 'C'), 30, ...
134+
'Direction', 'lower', 'Label', 'LL (idle+C)', ...
135+
'Color', [0.0 0.5 0.8], 'LineStyle', ':');
136+
137+
fprintf(' Rules: %d (%d upper, %d lower)\n', ...
138+
numel(s.ThresholdRules), 2, 2);
139+
140+
% ------------------------------------------------------------------
141+
% 4. MEX availability check
142+
% MEX files live in private/ dirs, so exist() can't see them from
143+
% here. Probe the actual file paths instead.
144+
% ------------------------------------------------------------------
145+
ext = mexext();
146+
fsPriv = fullfile(repo_root, 'libs', 'FastSense', 'private');
147+
stPriv = fullfile(repo_root, 'libs', 'SensorThreshold', 'private');
148+
mexProbes = {
149+
'binary_search_mex', fsPriv
150+
'compute_violations_mex', stPriv
151+
'to_step_function_mex', stPriv
152+
};
153+
fprintf('\n MEX status:\n');
154+
for i = 1:size(mexProbes, 1)
155+
mexFile = fullfile(mexProbes{i,2}, [mexProbes{i,1} '.' ext]);
156+
if exist(mexFile, 'file')
157+
fprintf(' %-30s compiled\n', mexProbes{i,1});
158+
else
159+
fprintf(' %-30s NOT compiled (MATLAB fallback)\n', mexProbes{i,1});
160+
end
161+
end
162+
163+
% ------------------------------------------------------------------
164+
% 5. JIT warmup — resolve a tiny sensor to force-compile all code paths
165+
% ------------------------------------------------------------------
166+
fprintf('\n JIT warmup... ');
167+
tWarm = tic;
168+
sw = Sensor('warmup');
169+
sw.X = [0 1 2 3]; sw.Y = [50 60 40 70];
170+
scw1 = StateChannel('machine'); scw1.X = [0 1 2]; scw1.Y = [0 1 2];
171+
scw2 = StateChannel('zone'); scw2.X = [0 1]; scw2.Y = {'A','B'};
172+
sw.addStateChannel(scw1); sw.addStateChannel(scw2);
173+
sw.addThresholdRule(struct('machine', 1), 55, 'Direction', 'upper');
174+
sw.addThresholdRule(struct('machine', 2), 35, 'Direction', 'lower');
175+
sw.addThresholdRule(struct('machine', 1, 'zone', 'B'), 50, 'Direction', 'upper');
176+
sw.addThresholdRule(struct('machine', 0, 'zone', 'A'), 40, 'Direction', 'lower');
177+
sw.resolve();
178+
clear sw scw1 scw2;
179+
fprintf('%.3f s\n', toc(tWarm));
180+
181+
% ------------------------------------------------------------------
182+
% 6. Run resolve() with timing
183+
% ------------------------------------------------------------------
184+
nRuns = 3;
185+
fprintf('\n Resolving (%d runs)...\n', nRuns);
186+
187+
times = zeros(1, nRuns);
188+
for r = 1:nRuns
189+
% Rebuild sensor each run to avoid cache effects
190+
sr = Sensor('stress_test');
191+
sr.X = x; sr.Y = y;
192+
sr.addStateChannel(scMachine);
193+
sr.addStateChannel(scZone);
194+
sr.addThresholdRule(struct('machine', 1), 75, ...
195+
'Direction', 'upper', 'Label', 'HH (running)');
196+
sr.addThresholdRule(struct('machine', 2), 25, ...
197+
'Direction', 'lower', 'Label', 'LL (evacuated)');
198+
sr.addThresholdRule(struct('machine', 1, 'zone', 'B'), 65, ...
199+
'Direction', 'upper', 'Label', 'HH (running+B)');
200+
sr.addThresholdRule(struct('machine', 0, 'zone', 'C'), 30, ...
201+
'Direction', 'lower', 'Label', 'LL (idle+C)');
202+
203+
tic;
204+
sr.resolve();
205+
times(r) = toc;
206+
207+
fprintf(' Run %d: %.3f s\n', r, times(r));
208+
end
209+
210+
% Use the last run's sensor for result inspection
211+
s = sr;
212+
213+
% ------------------------------------------------------------------
214+
% 7. Results summary
215+
% ------------------------------------------------------------------
216+
fprintf('\n============================================================\n');
217+
fprintf(' RESULTS\n');
218+
fprintf('============================================================\n');
219+
fprintf(' Data points: %s\n', format_size(N));
220+
fprintf(' State transitions: %d (machine) + %d (zone)\n', ...
221+
numel(machineX), numel(zoneX));
222+
fprintf(' Threshold rules: %d\n', numel(s.ThresholdRules));
223+
fprintf(' Resolved lines: %d\n', numel(s.ResolvedThresholds));
224+
225+
totalViol = 0;
226+
for v = 1:numel(s.ResolvedViolations)
227+
nv = numel(s.ResolvedViolations(v).X);
228+
totalViol = totalViol + nv;
229+
fprintf(' [%s] %s: %s violations\n', ...
230+
s.ResolvedViolations(v).Direction, ...
231+
s.ResolvedViolations(v).Label, ...
232+
format_size(nv));
233+
end
234+
fprintf(' Total violations: %s\n', format_size(totalViol));
235+
fprintf('\n');
236+
fprintf(' resolve() time:\n');
237+
fprintf(' median: %.3f s\n', median(times));
238+
fprintf(' min: %.3f s\n', min(times));
239+
fprintf(' max: %.3f s\n', max(times));
240+
fprintf(' throughput: %.1f M pts/s\n', (N / 1e6) / median(times));
241+
242+
% Memory estimate (approximate)
243+
memBytes = N * 8 * 2; % X + Y arrays
244+
for v = 1:numel(s.ResolvedViolations)
245+
memBytes = memBytes + numel(s.ResolvedViolations(v).X) * 8 * 2;
246+
end
247+
for t = 1:numel(s.ResolvedThresholds)
248+
memBytes = memBytes + numel(s.ResolvedThresholds(t).X) * 8 * 2;
249+
end
250+
fprintf(' approx memory: %.1f GB (sensor) + %.1f MB (results)\n', ...
251+
N * 8 * 2 / 1e9, (memBytes - N*8*2) / 1e6);
252+
253+
% ------------------------------------------------------------------
254+
% 8. Plot with FastSense
255+
% ------------------------------------------------------------------
256+
fprintf('\n Rendering with FastSense...\n');
257+
tPlot = tic;
258+
259+
fp = FastSense();
260+
fp.addSensor(s, 'ShowThresholds', true);
261+
fp.render();
262+
title('Sensor.resolve() Stress Test — 500M pts, 4 thresholds, ~9K state transitions');
263+
xlabel('Time [s]');
264+
ylabel('Process Value [bar]');
265+
266+
fprintf(' Render time: %.3f s\n', toc(tPlot));
267+
268+
fprintf('\n============================================================\n');
269+
fprintf(' Done.\n');
270+
fprintf('============================================================\n');
271+
end
272+
273+
274+
function s = format_size(N)
275+
if N >= 1e9
276+
s = sprintf('%.1fB', N / 1e9);
277+
elseif N >= 1e6
278+
s = sprintf('%.1fM', N / 1e6);
279+
elseif N >= 1e3
280+
s = sprintf('%.1fK', N / 1e3);
281+
else
282+
s = sprintf('%d', N);
283+
end
284+
end

‎install.m‎

Lines changed: 62 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,8 @@ function install()
1212
% 2. Adds examples, benchmarks, and tests to the path
1313
% 3. Compiles MEX accelerators if not yet built (first run only)
1414
% 4. Verifies the installation (first run only)
15+
% 5. JIT warmup — runs a tiny end-to-end workflow to force-compile
16+
% all hot code paths (once per MATLAB session)
1517
%
1618
% Directories added:
1719
% libs/FastSense — core plotting engine
@@ -60,6 +62,9 @@ function install()
6062
if needs_build(root)
6163
first_run(root);
6264
end
65+
66+
% --- Once per session: JIT warmup ---
67+
jit_warmup();
6368
end
6469

6570
function yes = needs_build(root)
@@ -69,9 +74,19 @@ function install()
6974
return;
7075
end
7176
mex_dir = fullfile(root, 'libs', 'FastSense', 'private');
72-
probe = fullfile(mex_dir, ['binary_search_mex.' mexext()]);
73-
probe_oct = fullfile(mex_dir, 'binary_search_mex.mex');
74-
yes = exist(probe, 'file') ~= 3 && exist(probe_oct, 'file') ~= 3;
77+
sensor_dir = fullfile(root, 'libs', 'SensorThreshold', 'private');
78+
% Probe a representative MEX from each library — if any are missing,
79+
% trigger build_mex() which will compile only the missing ones.
80+
probes = {
81+
fullfile(mex_dir, ['binary_search_mex.' mexext()])
82+
fullfile(mex_dir, 'binary_search_mex.mex')
83+
fullfile(sensor_dir, ['to_step_function_mex.' mexext()])
84+
fullfile(sensor_dir, 'to_step_function_mex.mex')
85+
};
86+
% Need build if either probe set is missing
87+
core_ok = exist(probes{1}, 'file') == 3 || exist(probes{2}, 'file') == 3;
88+
step_ok = exist(probes{3}, 'file') == 3 || exist(probes{4}, 'file') == 3;
89+
yes = ~core_ok || ~step_ok;
7590
end
7691

7792
function first_run(root)
@@ -159,3 +174,47 @@ function verify_installation(root)
159174
fprintf(' %d/%d checks passed, %d warnings\n', n_ok, total, n_warn);
160175
end
161176
end
177+
178+
179+
function jit_warmup()
180+
%JIT_WARMUP Force MATLAB's JIT to compile all hot code paths once per session.
181+
% Runs a tiny end-to-end workflow (sensor creation, state channels,
182+
% threshold resolve, downsampling, rendering setup) on trivial data.
183+
% Subsequent calls are no-ops. Adds ~0.1-0.3 s on first install() call
184+
% but eliminates multi-second JIT overhead on real workloads.
185+
persistent warmedUp;
186+
if ~isempty(warmedUp)
187+
return;
188+
end
189+
warmedUp = true;
190+
191+
try
192+
% --- Sensor + StateChannel + resolve pipeline ---
193+
sw = Sensor('__jit_warmup__');
194+
sw.X = [0 1 2 3 4 5];
195+
sw.Y = [50 60 40 70 30 80];
196+
197+
sc1 = StateChannel('s1');
198+
sc1.X = [0 2 4]; sc1.Y = [0 1 2];
199+
sc2 = StateChannel('s2');
200+
sc2.X = [0 3]; sc2.Y = {'A', 'B'};
201+
sw.addStateChannel(sc1);
202+
sw.addStateChannel(sc2);
203+
204+
sw.addThresholdRule(struct('s1', 1), 65, 'Direction', 'upper');
205+
sw.addThresholdRule(struct('s1', 2), 35, 'Direction', 'lower');
206+
sw.addThresholdRule(struct('s1', 1, 's2', 'B'), 55, 'Direction', 'upper');
207+
sw.addThresholdRule(struct('s1', 0, 's2', 'A'), 45, 'Direction', 'lower');
208+
sw.resolve();
209+
210+
% --- FastSense downsample + render setup (hidden figure) ---
211+
fig = figure('Visible', 'off');
212+
ax = axes('Parent', fig);
213+
fp = FastSense('Parent', ax);
214+
fp.addSensor(sw, 'ShowThresholds', true);
215+
fp.render();
216+
close(fig);
217+
catch
218+
% Warmup is best-effort — never block install on failure
219+
end
220+
end

‎libs/FastSense/build_mex.m‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@ function build_mex()
3535
% lttb_core_mex.c — Largest-Triangle-Three-Buckets kernel
3636
% violation_cull_mex.c — threshold violation culling
3737
% compute_violations_mex.c — batch violation detection for resolve()
38+
% to_step_function_mex.c — SIMD step-function conversion for thresholds
3839
% build_store_mex.c — bulk SQLite writer for DataStore init
3940
% mksqlite.c — SQLite3 MEX interface (bundled sqlite3.c)
4041
%
@@ -137,6 +138,7 @@ function build_mex()
137138
'compute_violations_mex.c', 'compute_violations_mex', {{}}, {{}}
138139
'resolve_disk_mex.c', 'resolve_disk_mex', {{sqlite3_src}}, {sqlite3_flags}
139140
'build_store_mex.c', 'build_store_mex', {{sqlite3_src}}, {sqlite3_flags}
141+
'to_step_function_mex.c', 'to_step_function_mex', {{}}, {{}}
140142
};
141143

142144
mksqlite_src = fullfile(rootDir, 'mksqlite.c');
@@ -228,6 +230,7 @@ function build_mex()
228230
copy_mex_to(outDir, sensorPrivDir, 'violation_cull_mex');
229231
copy_mex_to(outDir, sensorPrivDir, 'compute_violations_mex');
230232
copy_mex_to(outDir, sensorPrivDir, 'resolve_disk_mex');
233+
copy_mex_to(outDir, sensorPrivDir, 'to_step_function_mex');
231234
end
232235

233236
function compile_mex(src_file, out_name, outDir, include_flag, opt_flags, compiler, extra_srcs)

0 commit comments

Comments
 (0)