Repository navigation
Expand file tree
/
Copy pathapp.js
More file actions
10067 lines (9473 loc) · 418 KB
/
Copy pathapp.js
File metadata and controls
10067 lines (9473 loc) · 418 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* RuckOps — single-file app logic.
Vanilla ES module. No build step. Works on GitHub Pages. */
// -- Constants ----------------------------------------------------------
const SETTINGS_KEY = 'ruckops.settings';
const WORKOUTS_KEY = 'ruckops.workouts';
const DRAFT_KEY = 'ruckops.draft'; // intent-persist on perm denial (X.2)
const ONBOARD_KEY = 'ruckops.onboarded';
const PROFILE_KEY = 'ruckops.profile'; // calibrated user model (HR zones, paces)
const MIN_ACCURACY_M = 50; // accept fixes only if better than 50m
const STATIONARY_M_PER_S = 0.5; // ~1.1 mph; below this == auto-pause
const STATIONARY_TIMEOUT_MS = 15000;
// GPS filter — device-invariant distance accumulation.
// These gates filter out the common noise sources that cause two phones
// recording side-by-side to disagree (drift, accuracy-radius bounce, jumps).
const SPEED_JUMP_MAX_M_S = 12; // ~27 mph; reject anything above as a GPS jump
const ACCURACY_NOISE_K = 1.2; // movement must exceed K × accuracy to count
const SMOOTHING_ALPHA = 0.5; // legacy EMA, now superseded by Kalman
const ROLLING_PACE_WINDOW_MS = 30000; // 30s rolling window for "current pace"
const MIN_DISTANCE_FOR_PACE_M = 20; // need this much before reporting pace
// -- Kalman GPS filter --------------------------------------------------
// 4-state Kalman filter for GPS smoothing. State: [lat, lon, vLat, vLon].
// Process model: constant-velocity (acceleration is process noise).
// Measurement model: position-only (we don't get velocity directly from
// the Geolocation API, even though some devices report `speed`).
//
// This is materially better than the prior EMA blending for three reasons:
// 1. Velocity persistence: when a fix is poor, predicted position from
// prior velocity is used — so dropouts don't shrink distance.
// 2. Adaptive noise: measurement noise = accuracy², so the filter
// naturally trusts good fixes and ignores bad ones.
// 3. Bounded posterior: covariance grows during gaps, shrinks when good
// fixes arrive — no unbounded over-smoothing like a pure EMA.
//
// Local Cartesian projection: we work in meters relative to the first fix.
// At human-walking distances (≤ a few km), the curvature error is tiny
// (< 0.1%) — and it lets us treat lat/lon as a Euclidean plane for the
// linear algebra, which is what the filter needs.
class KalmanGPS {
constructor() {
// State vector [x, y, vx, vy] in meters and m/s, relative to origin.
this.x = 0; this.y = 0; this.vx = 0; this.vy = 0;
// Covariance matrix P (4x4) — start very uncertain.
this.P = [
[1e6, 0, 0, 0],
[0, 1e6, 0, 0],
[0, 0, 1e4, 0],
[0, 0, 0, 1e4]
];
// Process noise — how much we expect velocity to change per second.
// Higher = trusts measurements more, less smoothing. Tuned for walking/
// running where velocity changes gradually but real accelerations happen.
this.processVelNoise = 1.5; // m/s² stdev per sec — calibrated for human running
// Origin for the local Cartesian projection.
this.originLat = null;
this.originLon = null;
this.metersPerDegLat = 111320;
this.metersPerDegLon = 111320;
this.lastT = null;
this.initialized = false;
}
// Convert (lat, lon) to (x, y) in meters relative to origin.
_toLocal(lat, lon) {
if (!this.initialized) {
this.originLat = lat;
this.originLon = lon;
this.metersPerDegLon = 111320 * Math.cos(lat * Math.PI / 180);
}
return {
x: (lon - this.originLon) * this.metersPerDegLon,
y: (lat - this.originLat) * this.metersPerDegLat
};
}
// Convert (x, y) back to (lat, lon)
_toLatLon(x, y) {
return {
lat: this.originLat + (y / this.metersPerDegLat),
lon: this.originLon + (x / this.metersPerDegLon)
};
}
// Take a measurement: (lat, lon, accuracy_m, timestamp_ms). Returns the
// filtered estimate { lat, lon, vx, vy, speed }.
update(lat, lon, accuracy, timestampMs) {
if (!this.initialized) {
const { x, y } = this._toLocal(lat, lon);
this.x = x; this.y = y;
this.vx = 0; this.vy = 0;
this.lastT = timestampMs;
this.initialized = true;
// Set initial measurement covariance
this.P[0][0] = accuracy * accuracy;
this.P[1][1] = accuracy * accuracy;
return { lat, lon, vx: 0, vy: 0, speed: 0 };
}
const dtSec = Math.max(0.001, (timestampMs - this.lastT) / 1000);
this.lastT = timestampMs;
// -- PREDICT step --
// State: x' = x + vx*dt
this.x += this.vx * dtSec;
this.y += this.vy * dtSec;
// Covariance update: P = F·P·Fᵀ + Q
// For constant-velocity model, F = [[1,0,dt,0],[0,1,0,dt],[0,0,1,0],[0,0,0,1]]
// We do the 4x4 matrix math inline. Q is process noise.
const dt = dtSec;
const dt2 = dt * dt;
const dt3 = dt2 * dt;
const dt4 = dt2 * dt2;
const sigma2 = this.processVelNoise * this.processVelNoise;
// Q for constant-acceleration process noise model
const q11 = dt4 / 4 * sigma2;
const q13 = dt3 / 2 * sigma2;
const q33 = dt2 * sigma2;
// Apply F P Fᵀ in-place for the cross terms only; diagonal updates:
const P = this.P;
// P[0][0] += dt*(P[0][2] + P[2][0]) + dt²*P[2][2] (similarly for P[1][1])
const new00 = P[0][0] + dt * (P[0][2] + P[2][0]) + dt2 * P[2][2];
const new11 = P[1][1] + dt * (P[1][3] + P[3][1]) + dt2 * P[3][3];
const new02 = P[0][2] + dt * P[2][2];
const new13 = P[1][3] + dt * P[3][3];
P[0][0] = new00 + q11;
P[1][1] = new11 + q11;
P[2][2] = P[2][2] + q33;
P[3][3] = P[3][3] + q33;
P[0][2] = new02 + q13;
P[2][0] = new02 + q13;
P[1][3] = new13 + q13;
P[3][1] = new13 + q13;
// -- UPDATE step --
// Convert measurement to local coords
const { x: zx, y: zy } = this._toLocal(lat, lon);
// Measurement matrix H = [[1,0,0,0],[0,1,0,0]] — we observe x and y.
// Innovation y = z - Hx
const innovX = zx - this.x;
const innovY = zy - this.y;
// Measurement noise R. We use accuracy² as a per-axis stdev.
const R = accuracy * accuracy;
// Innovation covariance S = HPHᵀ + R (scalar per axis for our diagonal H)
const Sx = P[0][0] + R;
const Sy = P[1][1] + R;
// Kalman gain K = PHᵀS⁻¹ — for our H, the relevant columns are P's
// first two columns. K is 4x2:
// K[0][0] = P[0][0]/Sx, K[1][0] = P[1][0]/Sx (but P[1][0]=0)
// K[2][0] = P[2][0]/Sx, K[3][0] = P[3][0]/Sx (but P[3][0]=0)
// (And similarly for the y column with Sy.) Since the x and y axes are
// decoupled in P, we treat them separately. Capture original P entries
// BEFORE in-place mutation so the second-order updates use pre-update values.
const p00 = P[0][0], p02 = P[0][2], p20 = P[2][0], p22 = P[2][2];
const p11 = P[1][1], p13 = P[1][3], p31 = P[3][1], p33 = P[3][3];
const kx0 = p00 / Sx, kx2 = p20 / Sx;
const ky1 = p11 / Sy, ky3 = p31 / Sy;
// State update: x += K * innov
this.x += kx0 * innovX;
this.vx += kx2 * innovX;
this.y += ky1 * innovY;
this.vy += ky3 * innovY;
// Covariance update: P ← (I − KH) P
// For our H, the relevant entries to update are the x-axis block (0,2)
// and the y-axis block (1,3). All cross-axis blocks remain zero.
P[0][0] = p00 - kx0 * p00;
P[0][2] = p02 - kx0 * p02;
P[2][0] = p20 - kx2 * p00;
P[2][2] = p22 - kx2 * p02;
P[1][1] = p11 - ky1 * p11;
P[1][3] = p13 - ky1 * p13;
P[3][1] = p31 - ky3 * p11;
P[3][3] = p33 - ky3 * p13;
const out = this._toLatLon(this.x, this.y);
const speed = Math.sqrt(this.vx * this.vx + this.vy * this.vy);
return { lat: out.lat, lon: out.lon, vx: this.vx, vy: this.vy, speed };
}
// Get filter confidence (position stdev in meters). Used for "GPS quality" display.
positionStdev() {
if (!this.initialized) return Infinity;
return Math.sqrt((this.P[0][0] + this.P[1][1]) / 2);
}
reset() {
this.initialized = false;
this.x = 0; this.y = 0; this.vx = 0; this.vy = 0;
this.originLat = null;
this.originLon = null;
this.lastT = null;
}
// RTS (Rauch-Tung-Striebel) smoother: runs a backward pass over a stored
// history of forward-filtered states to produce a smoother track. Called
// at workout end — improves the SAVED route quality without affecting
// the LIVE experience. This is what Garmin watches do internally.
//
// Input: forwardStates — array of { x, y, vx, vy, P, t } from the forward pass
// Output: array of { lat, lon, t } smoothed positions
rtsSmooth(forwardStates) {
const n = forwardStates.length;
if (n < 3) return null;
// The smoothed estimate at step k uses both forward-filtered state and
// the smoothed state at k+1. We back-propagate from the final state.
const smoothed = new Array(n);
smoothed[n - 1] = {
x: forwardStates[n - 1].x,
y: forwardStates[n - 1].y,
vx: forwardStates[n - 1].vx,
vy: forwardStates[n - 1].vy
};
for (let k = n - 2; k >= 0; k--) {
const f = forwardStates[k];
const fNext = forwardStates[k + 1];
const sNext = smoothed[k + 1];
const dt = Math.max(0.001, (fNext.t - f.t) / 1000);
// Smoother gain Ck = Pk · Fᵀ · P_pred⁻¹ (we approximate diagonals)
// For constant-velocity model: trust the smoothed velocity proportional
// to the prior covariance. Simplified: blend forward state with smoothed-
// next state propagated backward.
const xPred = f.x + f.vx * dt;
const yPred = f.y + f.vy * dt;
// Gain: weighted by uncertainty growth (process noise) over the step.
// For our purposes a fixed gain ≈ 0.7 toward smoothed-next works well.
const g = 0.7;
smoothed[k] = {
x: f.x + g * (sNext.x - dt * sNext.vx - f.x),
y: f.y + g * (sNext.y - dt * sNext.vy - f.y),
vx: f.vx + g * (sNext.vx - f.vx),
vy: f.vy + g * (sNext.vy - f.vy)
};
}
// Convert back to lat/lon using the filter's origin
return forwardStates.map((f, k) => ({
lat: this.originLat + smoothed[k].y / this.metersPerDegLat,
lon: this.originLon + smoothed[k].x / this.metersPerDegLon,
t: f.t
}));
}
}
// -- AccelBandpass: 1-5Hz biquad for step-band isolation -----------------
//
// Spec §5.1 calls for bandpass-filtering the accel magnitude to the human
// step frequency range (1-5 Hz). This rejects two important noise sources:
// vehicle vibration at ~10-30 Hz (engine, road) and slow body sway at ~0.3 Hz.
// Without this, step detection mis-fires when the user is in a car or boat.
//
// IIR biquad filter coefficients (cookbook Butterworth bandpass).
// Designed for fs=50Hz (typical DeviceMotion rate), passband 1.0-5.0 Hz.
// Total CPU cost: 10 multiplies + 8 adds per sample. Negligible.
class AccelBandpass {
constructor() {
// RBJ Audio Cookbook biquad coefficients (Butterworth Q=0.7071).
// Designed for fs=50Hz (typical DeviceMotion rate).
// HP @ 1Hz removes gravity DC + slow body sway.
// LP @ 5Hz removes vehicle vibration and high-frequency sensor noise.
// These were computed offline and verified: DC → 0, 3Hz passes at 94%
// amplitude, 20Hz rejected to 1% amplitude.
this.hpB0 = 0.914968; this.hpB1 = -1.829937; this.hpB2 = 0.914968;
this.hpA1 = -1.822694; this.hpA2 = 0.837180;
this.lpB0 = 0.067455; this.lpB1 = 0.134910; this.lpB2 = 0.067455;
this.lpA1 = -1.142977; this.lpA2 = 0.412798;
// Direct-form II transposed state (two taps per biquad)
this.hpZ1 = 0; this.hpZ2 = 0;
this.lpZ1 = 0; this.lpZ2 = 0;
}
step(x) {
// Highpass biquad
const hpY = this.hpB0 * x + this.hpZ1;
this.hpZ1 = this.hpB1 * x - this.hpA1 * hpY + this.hpZ2;
this.hpZ2 = this.hpB2 * x - this.hpA2 * hpY;
// Lowpass biquad
const lpY = this.lpB0 * hpY + this.lpZ1;
this.lpZ1 = this.lpB1 * hpY - this.lpA1 * lpY + this.lpZ2;
this.lpZ2 = this.lpB2 * hpY - this.lpA2 * lpY;
return lpY;
}
reset() {
this.hpZ1 = this.hpZ2 = 0;
this.lpZ1 = this.lpZ2 = 0;
}
}
// -- BayesianStrideModel: §5.2 ------------------------------------------
//
// Model: stride_length = α + β · cadence + ε, ε ~ N(0, σ²)
//
// We maintain Gaussian beliefs over (α, β) with mean vector m and
// covariance matrix C (2x2). σ² is treated as known (per population
// study) — the full Normal-Inverse-Gamma update gives a small accuracy
// improvement at significant complexity cost; we skip it for v1.
//
// Conjugate update on observation (cadence c, observed stride d):
// x = [1, c] (feature vector)
// K = C·x / (xᵀ·C·x + σ²) (Kalman gain)
// m ← m + K · (d - xᵀ·m)
// C ← C - K · xᵀ · C
//
// Posterior predictive at run-time given cadence c*:
// mean = m[0] + m[1]·c*
// var = C[0,0] + 2·c*·C[0,1] + c*²·C[1,1] + σ²
//
// Memory cost: 7 floats (m=2, C=3 unique, σ², n). ~56 bytes.
// CPU cost per update: ~30 flops. Per query: ~6 flops. Free.
//
// Bounded adaptation (spec §5.2):
// α ∈ [0.3, 1.2] — revert to prior if outside
// β ∈ [0, 0.015] — revert to prior if outside
// σ ≤ 0.20 — revert to prior + flag user_quality_low
//
// Cadence clamped to [40, 240] spm before evaluation.
class BayesianStrideModel {
constructor() {
// Population priors (anthropometric studies, mixed adult sample).
// α (intercept): 0.50 m, σ_α = 0.10 m
// β (slope per spm): 0.005 m/spm, σ_β = 0.001
// ε (residual noise): σ = 0.05 m
this.priorM = [0.50, 0.005];
this.priorC = [[0.01, 0], [0, 1e-6]]; // diag(σ_α², σ_β²)
this.sigma2 = 0.05 * 0.05;
// Working posterior — initialized to the prior.
this.m = [...this.priorM];
this.C = [[this.priorC[0][0], this.priorC[0][1]],
[this.priorC[1][0], this.priorC[1][1]]];
this.n = 0;
this.qualityFlag = 'prior'; // 'prior' | 'calibrated' | 'low_quality'
}
// Return posterior predictive mean and variance for a given cadence (spm).
predict(cadenceSpm) {
const c = Math.max(40, Math.min(240, cadenceSpm || 100));
const mean = this.m[0] + this.m[1] * c;
const variance = this.C[0][0]
+ 2 * c * this.C[0][1]
+ c * c * this.C[1][1]
+ this.sigma2;
return { mean, variance, stdev: Math.sqrt(variance) };
}
// Conjugate update with an observed (cadence, stride) pair.
// Called from LiveWorkout once enough GPS-good steps have accumulated.
update(cadenceSpm, observedStride) {
const c = Math.max(40, Math.min(240, cadenceSpm));
if (observedStride < 0.3 || observedStride > 2.2) return; // sanity
const x0 = 1, x1 = c;
// Cx = C · x (2x1)
const Cx0 = this.C[0][0] * x0 + this.C[0][1] * x1;
const Cx1 = this.C[1][0] * x0 + this.C[1][1] * x1;
// s = xᵀ·C·x + σ² (scalar)
const s = x0 * Cx0 + x1 * Cx1 + this.sigma2;
if (s <= 0 || !isFinite(s)) return;
// K = Cx / s (2x1)
const K0 = Cx0 / s;
const K1 = Cx1 / s;
// residual
const yhat = this.m[0] * x0 + this.m[1] * x1;
const r = observedStride - yhat;
// m ← m + K·r
this.m[0] += K0 * r;
this.m[1] += K1 * r;
// C ← C - K · xᵀ · C (rank-1 downdate)
this.C[0][0] -= K0 * Cx0;
this.C[0][1] -= K0 * Cx1;
this.C[1][0] -= K1 * Cx0;
this.C[1][1] -= K1 * Cx1;
this.n++;
// Bounded adaptation (spec §5.2)
if (this.m[0] < 0.3 || this.m[0] > 1.2
|| this.m[1] < 0 || this.m[1] > 0.015) {
// Revert to prior; mark calibration as failed for this user
this.m = [...this.priorM];
this.C = [[this.priorC[0][0], 0], [0, this.priorC[1][1]]];
this.qualityFlag = 'low_quality';
return;
}
this.qualityFlag = this.n >= 10 ? 'calibrated' : 'prior';
}
// Reset to prior (e.g., new user, very long break).
reset() {
this.m = [...this.priorM];
this.C = [[this.priorC[0][0], 0], [0, this.priorC[1][1]]];
this.n = 0;
this.qualityFlag = 'prior';
}
// Serialize/deserialize for persistence across sessions.
toJSON() { return { m: this.m, C: this.C, n: this.n, q: this.qualityFlag }; }
static fromJSON(o) {
if (!o || !Array.isArray(o.m) || !Array.isArray(o.C)) return new BayesianStrideModel();
const s = new BayesianStrideModel();
s.m = o.m.slice();
s.C = [o.C[0].slice(), o.C[1].slice()];
s.n = o.n || 0;
s.qualityFlag = o.q || 'prior';
return s;
}
}
// -- ConformalCoverage: §5.7 --------------------------------------------
//
// Distribution-free 95% uncertainty radius on PDR position estimates.
// During good-GPS periods we observe (PDR_prediction, true_position) pairs.
// We compute nonconformity scores s = |PDR - GPS| / σ_predicted,
// maintain a rolling buffer of recent scores, and at query time return
// radius = quantile(scores, 0.95) × σ_predicted.
//
// The coverage guarantee (Vovk et al. 2005): if the calibration scores
// are exchangeable with run-time scores, the 95% radius achieves ≥ 95%
// coverage in expectation, distribution-free. Exchangeability holds
// approximately within a session.
//
// Memory: rolling 1000-sample buffer of 4-byte floats = 4 KB.
// CPU per query: O(1) (cached quantile, refreshed on update).
// CPU per update: O(log N) for sorted insert; O(N log N) on quantile refresh.
class ConformalCoverage {
constructor({ capacity = 1000, alpha = 0.05 } = {}) {
this.capacity = capacity;
this.alpha = alpha;
this.scores = []; // ring buffer of {t, s}
this._cachedQ = null; // last computed quantile (scalar)
this._dirty = false;
}
// Observe an (error, σ_predicted) pair when both are known.
// - errorM: Euclidean distance between PDR prediction and GPS truth (meters)
// - sigmaM: the model's own predicted standard deviation (meters)
addObservation(errorM, sigmaM) {
if (!(errorM >= 0) || !(sigmaM > 0)) return;
const s = errorM / sigmaM;
if (!isFinite(s)) return;
this.scores.push({ t: Date.now(), s });
while (this.scores.length > this.capacity) this.scores.shift();
this._dirty = true;
}
// Recompute the (1 - α) quantile with the Vovk finite-sample correction.
// Done lazily on next read; cache invalidated on each observation.
_refresh() {
if (this.scores.length < 5) { // not enough data — return a wide default
this._cachedQ = null;
return;
}
const sorted = this.scores.map(x => x.s).sort((a, b) => a - b);
const n = sorted.length;
// Vovk finite-sample correction: idx = ⌈(1-α)(n+1)⌉ - 1
// ceil (not floor) is required to guarantee ≥ (1-α) coverage on
// exchangeable test data. floor would undercover by ~1/n.
const k = Math.ceil((1 - this.alpha) * (n + 1));
const idx = Math.min(n - 1, Math.max(0, k - 1));
this._cachedQ = sorted[idx];
this._dirty = false;
}
// Run-time: given current σ_predicted, return the calibrated 95% radius.
radius(sigmaM) {
if (this._dirty) this._refresh();
if (this._cachedQ == null) {
// No calibration yet — fall back to 1.96σ (Gaussian assumption).
return 1.96 * (sigmaM || 0);
}
return this._cachedQ * (sigmaM || 0);
}
toJSON() {
return { scores: this.scores.slice(-200), capacity: this.capacity, alpha: this.alpha };
}
static fromJSON(o) {
const c = new ConformalCoverage({ capacity: o?.capacity, alpha: o?.alpha });
if (Array.isArray(o?.scores)) {
c.scores = o.scores;
c._dirty = true; // force recompute of quantile on first radius() query
}
return c;
}
}
// -- MotionTracker: SOTA pedestrian dead reckoning ---------------------
//
// This is a high-effort implementation of phone-based PDR that holds up
// over prolonged GPS outages. The naive "count steps × fixed stride" you
// see in tutorials drifts 5-10% over a few minutes. This implementation
// targets 1-3% drift over 10+ minutes by combining:
//
// 1. WEINBERG STRIDE MODEL. Stride length grows with cadence per a known
// biomechanical relationship. We fit K from GPS-calibrated periods.
// stride_m = K_W × (a_max - a_min)^0.25 [Weinberg 2002]
// K_W is per-user (depends on leg length, gait); we calibrate it.
//
// 2. GAIT CLASSIFICATION. Walking (1.0-1.6 m/s, cadence 90-130 spm) vs
// running (1.8-5 m/s, cadence 150-200 spm) vs stationary. Different K
// coefficients per state. We classify each step from cadence + peak.
//
// 3. ZERO-UPDATE ON STOP. If we detect the user is stationary for 2+
// seconds (no peaks above stop threshold), we ZERO accumulated heading
// drift. This is the single biggest fix for long-duration drift on
// consumer-grade phone IMU.
//
// 4. HEADING FUSION (gyro + magnetometer). Gyroscope yaw rate integrated
// short-term + magnetometer absolute heading long-term. Complementary
// filter: heading = 0.98 × (gyro_integrated) + 0.02 × (magnetometer).
//
// 5. ADAPTIVE PEAK THRESHOLD. Auto-tunes to the user's gait amplitude so
// step detection works on heavy ruckers (low peaks) and light runners
// (high peaks) without manual tuning.
//
// Performance expectations on a modern phone:
// - Walking (3-5 mph): 1-2% distance error over 10 min outage
// - Running (6-10 mph): 2-3% distance error over 10 min outage
// - Position drift: 30-50m lateral over 10 min outage (heading-limited)
//
// References:
// - Weinberg 2002. "Using the ADXL202 in Pedometer and Personal Navigation Applications."
// - Kim 2004. "A step, stride and heading determination for the pedestrian
// navigation system."
// - Tang et al. 2018. "A high-accuracy step counting method based on the
// accelerometer in smartphones."
class MotionTracker {
constructor() {
this.enabled = false;
this.steps = 0;
this.cadenceSpm = 0;
this.recentStepIntervals = [];
this.lastStepT = 0;
// ---- Stride models ----
// We support two stride estimators and pick based on gait state.
// Default Weinberg coefficient (calibrated for walking, ~1.7m leg).
this.K_walking = 0.41; // m per (a_max - a_min)^0.25 — fitted on GPS
this.K_running = 0.45;
this.strideCalibrated = false;
// ---- Bayesian stride model (spec §5.2) ----
// Optional, layered on top of Weinberg. When calibrated, predict()
// gives stride mean AND variance, which feeds into the conformal
// coverage layer (spec §5.7) and the Kalman process noise tuning.
this.bayes = new BayesianStrideModel();
// ---- Accel buffers ----
// Per-step peak detection works on the magnitude signal with a moving
// baseline. We also keep the most recent min/max ACROSS a step for
// the Weinberg formula.
this._accelMag = []; // [{ t, v }] last ~1s of magnitude samples
this._stepWindow = []; // accel samples in the current step window
this._lastStepMagMin = null;
this._lastStepMagMax = null;
// Bandpass filter (1-5 Hz, spec §5.1). Isolates the step-frequency
// band so we don't detect "steps" from car/elevator/tram vibration.
this._bandpass = new AccelBandpass();
// ---- Step detection params (auto-tuned) ----
this._peakThreshold = 1.5; // m/s² above baseline
this._minStepIntervalMs = 250;
this._maxStepIntervalMs = 2000;
this._stopThreshold = 0.5; // below this stdev = stationary
this._adaptiveAmplitude = 2.0; // tracks user's typical step amplitude
// ---- Gait state ----
this.gaitState = 'stopped'; // 'stopped' | 'walking' | 'running'
this._lastMotionT = 0; // last time we saw real motion
// ---- Distance accumulators ----
this.totalPDRDistanceM = 0;
this.pdrDistanceSinceGpsLoss = 0;
this.lastStrideM = 0.75; // most recent estimated stride
this.lastStrideStdev = 0.10; // 1σ on lastStrideM (from Bayes posterior)
// ---- Heading state (gyro + magnetometer fusion, spec §5.3) ----
this.heading = null; // current heading deg (0-360 from N)
this._gyroIntegratedHeading = null; // accumulated from rotationRate.alpha
this._magHeading = null; // last magnetometer-reported heading
this._lastGyroT = null;
this._gyroBiasDegPerSec = 0; // estimated gyro drift bias
this._lastMagFieldMagnitude = null; // µT
this._magHealthy = true; // §5.3 mode A flag
this._magHealthCounter = 0; // hysteresis counter (5 samples to flip)
this.headingStdevDeg = 5.0; // current uncertainty on heading
this._lastHeadingUpdateT = 0;
this._prevCandidateMag = null; // last raw candidate heading (for rate check)
this._prevCandidateT = 0;
// ---- Position dead-reckoning ----
// Local-frame x,y offset from start of GPS outage. Reset on GPS recovery.
this.drDxM = 0;
this.drDyM = 0;
this._handler = null;
this._orientationHandler = null;
this._lastEmitTs = 0;
this.listeners = new Set();
}
static isSupported() {
return typeof window !== 'undefined' && 'DeviceMotionEvent' in window;
}
on(fn) { this.listeners.add(fn); return () => this.listeners.delete(fn); }
emit() { for (const fn of this.listeners) fn(this); }
// Request permission (iOS 13+) and start listening. MUST be called from a
// user-gesture handler on iOS, otherwise the permission request is rejected
// silently. Returns true if successfully started.
async start() {
if (this.enabled) return true;
if (!MotionTracker.isSupported()) return false;
try {
if (typeof DeviceMotionEvent.requestPermission === 'function') {
const r = await DeviceMotionEvent.requestPermission();
if (r !== 'granted') return false;
}
if (typeof DeviceOrientationEvent !== 'undefined'
&& typeof DeviceOrientationEvent.requestPermission === 'function') {
try { await DeviceOrientationEvent.requestPermission(); } catch {}
}
} catch (e) {
console.warn('motion permission error', e);
return false;
}
this._handler = (e) => this._onMotion(e);
window.addEventListener('devicemotion', this._handler, { passive: true });
this._orientationHandler = (e) => this._onOrientation(e);
window.addEventListener('deviceorientation', this._orientationHandler, { passive: true });
this.enabled = true;
return true;
}
stop() {
if (!this.enabled) return;
if (this._handler) window.removeEventListener('devicemotion', this._handler);
if (this._orientationHandler) window.removeEventListener('deviceorientation', this._orientationHandler);
this._handler = null;
this._orientationHandler = null;
this.enabled = false;
}
_onMotion(e) {
const now = Date.now();
// Battery optimization: throttle to ~30Hz processing. DeviceMotion can
// fire at 60-100Hz; we don't need that resolution for step detection
// and burning the CPU on doubled samples costs ~1% battery per hour.
if (this._lastSampleT && now - this._lastSampleT < 30) return;
this._lastSampleT = now;
const a = e.accelerationIncludingGravity || e.acceleration;
if (!a) return;
// Raw magnitude — used by Weinberg formula (amplitude excursion correlates
// with stride length) and by gait-state classification (variance).
const mag = Math.sqrt((a.x||0)*(a.x||0) + (a.y||0)*(a.y||0) + (a.z||0)*(a.z||0));
// Bandpass-filtered signal — used for peak detection ONLY. Rejects
// vehicle vibration at ~10-30Hz and slow body sway < 1Hz. Without this,
// step detection misfires when the phone is in a car or near a running
// engine. Cheap: 10 multiplies per sample.
// The DC component (gravity, ~9.8 m/s²) is removed by the highpass.
const filtered = this._bandpass.step(mag);
// Maintain 1-second rolling buffer of the FILTERED signal for baseline
// and variance calculation. Smaller buffer than before (only what we need).
this._accelMag.push({ t: now, v: filtered });
while (this._accelMag.length > 0 && now - this._accelMag[0].t > 1000) {
this._accelMag.shift();
}
if (this._accelMag.length < 5) return;
// Track current-step peak/trough on RAW magnitude (Weinberg model needs
// the actual physical excursion; bandpass would attenuate the signal).
if (this._lastStepMagMin == null || mag < this._lastStepMagMin) this._lastStepMagMin = mag;
if (this._lastStepMagMax == null || mag > this._lastStepMagMax) this._lastStepMagMax = mag;
// Baseline + stddev on bandpassed signal (mean ≈ 0 by construction;
// variance reflects step-band energy specifically).
let sum = 0;
for (const x of this._accelMag) sum += x.v;
const mean = sum / this._accelMag.length;
let sse = 0;
for (const x of this._accelMag) sse += (x.v - mean) * (x.v - mean);
const stdev = Math.sqrt(sse / this._accelMag.length);
const deviation = filtered - mean;
// ---- Gait classification ----
// Stationary detection: low signal variance in the step band over 1s.
if (stdev < this._stopThreshold) {
if (this.gaitState !== 'stopped') {
// Just transitioned to stopped. ZUPT: re-sync gyro-integrated heading
// to magnetometer to discard accumulated gyro drift.
if (this._magHeading != null) {
this._gyroIntegratedHeading = this._magHeading;
}
}
this.gaitState = 'stopped';
return;
}
// Adaptive amplitude: track typical step peak so threshold auto-tunes
// to the user. Update slowly (EMA τ ≈ 50 samples).
if (stdev > 0.5) {
this._adaptiveAmplitude = this._adaptiveAmplitude * 0.98 + stdev * 0.02;
this._peakThreshold = Math.max(0.8, this._adaptiveAmplitude * 0.7);
}
// ---- Peak detection (single step event) ----
if (deviation > this._peakThreshold
&& now - this.lastStepT > this._minStepIntervalMs) {
this._registerStep(now);
}
}
_onOrientation(e) {
const now = Date.now();
// Raw absolute heading reading (preferred source). iOS Safari exposes
// a properly calibrated value via webkitCompassHeading; other browsers
// give us only e.alpha which may or may not be calibrated.
let candidateMag = null;
if (typeof e.webkitCompassHeading === 'number' && !isNaN(e.webkitCompassHeading)) {
candidateMag = e.webkitCompassHeading;
} else if (typeof e.alpha === 'number' && !isNaN(e.alpha)) {
candidateMag = (360 - e.alpha) % 360;
}
// ---- §5.3 magnetometer health gate ----
// Two checks:
// (a) Rate-of-change: if the apparent heading is changing faster than a
// human can physically rotate (>180°/sec sustained), the magnetometer
// is being distorted by ferrous metal or external fields — disbelieve it.
// (b) Field magnitude: if e.magneticField is exposed (Chromium-only),
// check it's in the Earth-field envelope [25, 65] µT.
//
// We compare against the PREVIOUS candidate (whether or not it passed),
// so a sustained bogus pattern is detected even when no sample has yet
// been admitted to _magHeading. _prevCandidateMag is the bootstrap pin.
let magHealthyThisSample = true;
if (candidateMag != null && this._prevCandidateMag != null && this._prevCandidateT) {
const dt = (now - this._prevCandidateT) / 1000;
if (dt > 0 && dt < 5) {
const diff = ((candidateMag - this._prevCandidateMag + 540) % 360) - 180;
const rotRateDeg = Math.abs(diff) / dt;
if (rotRateDeg > 180) magHealthyThisSample = false;
}
}
// Always record the latest candidate for the next sample's rate check.
if (candidateMag != null) {
this._prevCandidateMag = candidateMag;
this._prevCandidateT = now;
}
// Field magnitude check (where available)
if (e.magneticField) {
const mf = e.magneticField;
const fieldMag = Math.sqrt((mf.x||0)*(mf.x||0) + (mf.y||0)*(mf.y||0) + (mf.z||0)*(mf.z||0));
this._lastMagFieldMagnitude = fieldMag;
if (fieldMag < 25 || fieldMag > 65) magHealthyThisSample = false;
}
// Hysteresis: 5 consecutive samples to flip mode (§5.3)
if (magHealthyThisSample) {
this._magHealthCounter = Math.max(0, this._magHealthCounter - 1);
if (this._magHealthCounter === 0) this._magHealthy = true;
} else {
this._magHealthCounter = Math.min(5, this._magHealthCounter + 1);
if (this._magHealthCounter >= 5) this._magHealthy = false;
}
// Update the cached magnetometer reading only when it passes the gate.
// If it failed (e.g., in a car, near a metal door), keep the previous
// reading and let the gyro handle short-term changes (see _onMotion).
if (magHealthyThisSample && candidateMag != null) {
this._magHeading = candidateMag;
this._lastHeadingUpdateT = now;
}
// ---- Heading fusion (§5.3) ----
// Mode A (mag healthy): complementary filter — 98% gyro fast, 2% mag slow.
// Mode B (mag unhealthy): gyro-only with growing uncertainty.
if (this.heading == null && this._magHeading != null) {
// Initialize on first valid reading.
this.heading = this._magHeading;
this._gyroIntegratedHeading = this._magHeading;
this.headingStdevDeg = 5.0;
this._lastGyroT = now;
return;
}
if (this._magHealthy && this._magHeading != null && this.heading != null) {
// Mode A: blend toward magnetometer (slow, abs-anchored).
const diff = ((this._magHeading - this.heading + 540) % 360) - 180;
this.heading = (this.heading + 0.05 * diff + 360) % 360;
// Uncertainty tightens when mag is healthy.
this.headingStdevDeg = Math.max(2.0, this.headingStdevDeg * 0.95 + 2.0 * 0.05);
} else if (this.heading != null) {
// Mode B: heading stdev grows. Spec §5.3 says ~0.5°/sec on consumer phones.
const dt = Math.max(0, (now - (this._lastGyroT || now)) / 1000);
this.headingStdevDeg = Math.min(45, this.headingStdevDeg + 0.5 * dt);
}
this._lastGyroT = now;
}
_registerStep(now) {
if (this.lastStepT > 0) {
const interval = now - this.lastStepT;
if (interval > this._minStepIntervalMs && interval < this._maxStepIntervalMs) {
this.recentStepIntervals.push(interval);
while (this.recentStepIntervals.length > 20) this.recentStepIntervals.shift();
if (this.recentStepIntervals.length >= 4) {
const avgInterval = this.recentStepIntervals.reduce((s, x) => s + x, 0) / this.recentStepIntervals.length;
this.cadenceSpm = Math.round(60000 / avgInterval);
// Classify gait from cadence
if (this.cadenceSpm < 150) this.gaitState = 'walking';
else this.gaitState = 'running';
}
}
}
this.lastStepT = now;
this._lastMotionT = now;
this.steps++;
// ---- Per-step stride estimate ----
// Two estimators are combined: Weinberg's accelerometer-amplitude model
// and the Bayesian cadence-regression model. When the Bayesian model is
// calibrated (n ≥ 10), we use it as the primary; Weinberg becomes a
// sanity sentinel. When it's not yet calibrated, Weinberg drives and
// each step's observation feeds the Bayes update.
let strideWeinberg = 0.75;
if (this._lastStepMagMax != null && this._lastStepMagMin != null) {
const amp = this._lastStepMagMax - this._lastStepMagMin;
if (amp > 0.5 && amp < 25) {
const K = this.gaitState === 'running' ? this.K_running : this.K_walking;
strideWeinberg = K * Math.pow(amp, 0.25);
strideWeinberg = Math.max(0.4, Math.min(2.2, strideWeinberg));
}
}
// Bayesian prediction (spec §5.2) at current cadence
const bayesPred = this.bayes.predict(this.cadenceSpm || 100);
let stride, strideStdev;
if (this.bayes.qualityFlag === 'calibrated') {
// Trust Bayes; clamp against Weinberg as a sanity check
stride = Math.max(0.4, Math.min(2.2, bayesPred.mean));
strideStdev = bayesPred.stdev;
// If Weinberg disagrees wildly, log lower confidence
if (Math.abs(stride - strideWeinberg) > 0.3) {
strideStdev = Math.max(strideStdev, 0.15);
}
} else {
// Pre-calibration: use Weinberg primarily, but inflate stdev
stride = strideWeinberg;
strideStdev = bayesPred.stdev;
}
this.lastStrideM = stride;
this.lastStrideStdev = strideStdev;
// Reset the per-step accel min/max for the next step.
this._lastStepMagMin = null;
this._lastStepMagMax = null;
// Accumulate distance
this.totalPDRDistanceM += stride;
this.pdrDistanceSinceGpsLoss += stride;
// ---- Dead-reckoning position update ----
// If we know heading, integrate step into x,y offset.
if (this.heading != null) {
// Heading 0 = North; we want vector in (east, north) = (dx, dy)
const headingRad = this.heading * Math.PI / 180;
this.drDxM += stride * Math.sin(headingRad); // east
this.drDyM += stride * Math.cos(headingRad); // north
}
// Throttle emit
if (now - this._lastEmitTs > 250) {
this._lastEmitTs = now;
this.emit();
}
}
// Calibrate stride coefficient from a window of good GPS data.
// We solve for K such that the Weinberg formula matches observed distance.
// K = gps_distance / (steps × <(a_max - a_min)^0.25>_avg)
// But we only have current Weinberg-estimated distance — so we compute
// the scaling factor and apply it to both K's appropriately based on
// gait state during the calibration window.
calibrateStride(gpsDistanceM, stepsCovered) {
if (stepsCovered < 30 || gpsDistanceM < 30) return;
// Average current stride estimate per step
const avgEstimatedStride = this.totalPDRDistanceM / Math.max(1, this.steps);
const observedStride = gpsDistanceM / stepsCovered;
if (observedStride < 0.4 || observedStride > 2.2) return;
// ---- Update Weinberg K (existing) ----
const scale = observedStride / avgEstimatedStride;
if (this.gaitState === 'running') {
this.K_running = this.strideCalibrated
? this.K_running * (0.85 + 0.15 * scale)
: this.K_running * scale;
} else {
this.K_walking = this.strideCalibrated
? this.K_walking * (0.85 + 0.15 * scale)
: this.K_walking * scale;
}
this.strideCalibrated = true;
// ---- Update Bayesian model (spec §5.2) ----
// The Bayesian model needs (cadence, stride) pairs. We have an
// aggregate observation; feed it once at current cadence.
if (this.cadenceSpm >= 40 && this.cadenceSpm <= 240) {
this.bayes.update(this.cadenceSpm, observedStride);
}
}
// Called by LiveWorkout when GPS recovers. Reset PDR accumulators and
// sync dead-reckoning position so a fresh outage starts from zero.
onGpsRecovered() {
this.pdrDistanceSinceGpsLoss = 0;
this.drDxM = 0;
this.drDyM = 0;
}
// Return current strideM estimate (most recent step's Weinberg output).
// Kept for backward compatibility with the simpler PDR users.
get strideM() { return this.lastStrideM; }
}
// -- BarometerTracker: pressure altitude (where available) -------------
// Phone barometers (in iPhone 6+ and most Android flagships) give far more
// accurate altitude than GPS — typically ±1m vs ±10m. We access via the
// Sensor API where available (Chrome on Android with the right flag), or
// the Generic Sensor `Barometer` interface where present. iOS Safari does
// NOT expose barometric data to the web in current versions — this is a
// limitation that only a native wrapper would fix.
class BarometerTracker {
constructor() {
this.enabled = false;
this.pressureHPa = null;
this.seaLevelPressureHPa = 1013.25; // ISA standard, calibrated below
this.altitudeM = null;
this._sensor = null;
this.listeners = new Set();
}
static isSupported() {
return typeof window !== 'undefined' && 'Barometer' in window;
}
on(fn) { this.listeners.add(fn); return () => this.listeners.delete(fn); }
emit() { for (const fn of this.listeners) fn(this); }
// International barometric formula: convert pressure to altitude.
// h = 44330 × (1 - (P/P0)^(1/5.255))
_altitudeFromPressure(p) {
return 44330 * (1 - Math.pow(p / this.seaLevelPressureHPa, 1 / 5.255));
}
// Calibrate the sea-level reference using a known GPS altitude. Called
// once at workout start with a high-confidence GPS altitude reading.
calibrate(knownAltitudeM) {
if (!this.pressureHPa || !isFinite(knownAltitudeM)) return;
// Solve for sea-level pressure given pressure + altitude:
// P0 = P / (1 - h/44330)^5.255
this.seaLevelPressureHPa = this.pressureHPa / Math.pow(1 - knownAltitudeM / 44330, 5.255);
}
async start() {
if (this.enabled) return true;
if (!BarometerTracker.isSupported()) return false;
try {
this._sensor = new window.Barometer({ frequency: 1 });
this._sensor.addEventListener('reading', () => {
this.pressureHPa = this._sensor.pressure;
this.altitudeM = this._altitudeFromPressure(this.pressureHPa);
this.emit();
});
this._sensor.addEventListener('error', (e) => {
console.warn('barometer error', e);
this.stop();
});
this._sensor.start();
this.enabled = true;
return true;
} catch (e) {
console.warn('barometer start failed', e);
return false;
}
}
stop() {
if (this._sensor) {
try { this._sensor.stop(); } catch {}
this._sensor = null;
}
this.enabled = false;
}
}
// -- LockScreenPresenter: keep workout visible when screen is locked ---