Fix of an array processing error in the Theil-Sen regressor

app/engine/Flywheel.test.js

@@ -1,6 +1,6 @@
 'use strict'
 /*
-  Open Rowing Monitor, https://github.com/jaapvanekris/openrowingmonitor
+  Open Rowing Monitor, https://github.com/JaapvanEkris/openrowingmonitor
 */
 import { test } from 'uvu'
 import * as assert from 'uvu/assert'
@@ -90,9 +90,9 @@ test('Correct Flywheel behaviour for a noisefree stroke', () => {
   testDeltaTime(flywheel, 0.011062297)
   testSpinningTime(flywheel, 0.077918634)
   testAngularPosition(flywheel, 8.377580409572781)
-  testAngularVelocity(flywheel, 94.76231358849583)
-  testAngularAcceleration(flywheel, 28.980404808837132)
-  testTorque(flywheel, 3.975668304221995)
+  testAngularVelocity(flywheel, 94.77498684553687)
+  testAngularAcceleration(flywheel, 28.980405331480235)
+  testTorque(flywheel, 3.975932584148498)
   testDragFactor(flywheel, 0.00011)
   testIsDwelling(flywheel, false)
   testIsUnpowered(flywheel, false)
@@ -115,9 +115,9 @@ test('Correct Flywheel behaviour for a noisefree stroke', () => {
   testDeltaTime(flywheel, 0.010722165)
   testSpinningTime(flywheel, 0.23894732900000007)
   testAngularPosition(flywheel, 24.085543677521745)
-  testAngularVelocity(flywheel, 97.13471664858164)
-  testAngularAcceleration(flywheel, -29.657593800236377)
-  testTorque(flywheel, -2.0198310711803433)
+  testAngularVelocity(flywheel, 97.12541571421204)
+  testAngularAcceleration(flywheel, -29.657604177526746)
+  testTorque(flywheel, -2.0200308891605716)
   testDragFactor(flywheel, 0.00011)
   testIsDwelling(flywheel, false)
   testIsUnpowered(flywheel, true)
@@ -140,9 +140,9 @@ test('Correct Flywheel behaviour for a noisefree stroke', () => {
   testDeltaTime(flywheel, 0.020722165)
   testSpinningTime(flywheel, 0.43343548300000007)
   testAngularPosition(flywheel, 39.79350694547071)
-  testAngularVelocity(flywheel, 50.71501160141977)
-  testAngularAcceleration(flywheel, -159.90034506799844)
-  testTorque(flywheel, -16.202804212320103)
+  testAngularVelocity(flywheel, 50.85265548983507)
+  testAngularAcceleration(flywheel, -159.89027501034317)
+  testTorque(flywheel, -16.20022817082592)
   testDragFactor(flywheel, 0.00011)
   testIsDwelling(flywheel, true)
   testIsUnpowered(flywheel, true)

app/engine/Rower.test.js

@@ -1,5 +1,6 @@
 'use strict'
 /*
+  Open Rowing Monitor, https://github.com/JaapvanEkris/openrowingmonitor
 
   This test is a test of the Rower object, that tests wether this object fills all fields correctly, given one validated rower, (the
   Concept2 RowErg) using a validated cycle of strokes. This thoroughly tests the raw physics of the translation of Angular physics
@@ -110,11 +111,11 @@ test('Test behaviour for three perfect identical strokes, including settingling
   testDriveDuration(rower, 0) // Shouldn't this one be filled after the first drive?
   testDriveLinearDistance(rower, 0.2491943602992768)
   testDriveLength(rower, 0) // Shouldn't this one be filled after the first drive?
-  testDriveAverageHandleForce(rower, 1691.793078056684)
-  testDrivePeakHandleForce(rower, 10246.062011594136)
+  testDriveAverageHandleForce(rower, 249.91096328436572)
+  testDrivePeakHandleForce(rower, 280.43473478416803)
   testRecoveryDuration(rower, 0)
   testRecoveryDragFactor(rower, 110)
-  testInstantHandlePower(rower, 372.0199762100516)
+  testInstantHandlePower(rower, 372.09477620281604)
   // Recovery initial stroke starts here
   rower.handleRotationImpulse(0.010769)
   rower.handleRotationImpulse(0.010707554)
@@ -142,8 +143,8 @@ test('Test behaviour for three perfect identical strokes, including settingling
   testDriveDuration(rower, 0.19636192600000005)
   testDriveLinearDistance(rower, 0.6407854979124261)
   testDriveLength(rower, 0.2638937829015426)
-  testDriveAverageHandleForce(rower, 851.8820525641245) // This is the first stroke, which always leads to insane data like this
-  testDrivePeakHandleForce(rower, 10246.062011594136)
+  testDriveAverageHandleForce(rower, 247.35502383653122)
+  testDrivePeakHandleForce(rower, 325.1619554833936)
   testRecoveryDuration(rower, 0)
   testRecoveryDragFactor(rower, 110)
   testInstantHandlePower(rower, 0)
@@ -178,11 +179,11 @@ test('Test behaviour for three perfect identical strokes, including settingling
   testDriveDuration(rower, 0.19636192600000005)
   testDriveLinearDistance(rower, 0.4520822644211139)
   testDriveLength(rower, 0.2638937829015426)
-  testDriveAverageHandleForce(rower, 251.04336322997108)
-  testDrivePeakHandleForce(rower, 396.7011215867992)
+  testDriveAverageHandleForce(rower, 251.12896067596512)
+  testDrivePeakHandleForce(rower, 396.7733761783577)
   testRecoveryDuration(rower, 0.152533057)
   testRecoveryDragFactor(rower, 309.02744980039836)
-  testInstantHandlePower(rower, 526.5255378434941)
+  testInstantHandlePower(rower, 526.7173432408988)
   // Recovery second stroke starts here
   rower.handleRotationImpulse(0.010769)
   rower.handleRotationImpulse(0.010707554)
@@ -210,8 +211,8 @@ test('Test behaviour for three perfect identical strokes, including settingling
   testDriveDuration(rower, 0.25056694500000004)
   testDriveLinearDistance(rower, 1.1553213424095137)
   testDriveLength(rower, 0.3371976114853044)
-  testDriveAverageHandleForce(rower, 290.98159585708896)
-  testDrivePeakHandleForce(rower, 456.9929898648157)
+  testDriveAverageHandleForce(rower, 290.9778542238004)
+  testDrivePeakHandleForce(rower, 456.96817421319486)
   testRecoveryDuration(rower, 0.152533057)
   testRecoveryDragFactor(rower, 309.02744980039836) // As we decelerate the flywheel quite fast, this is expected
   testInstantHandlePower(rower, 0)
@@ -246,11 +247,11 @@ test('Test behaviour for three perfect identical strokes, including settingling
   testDriveDuration(rower, 0.25056694500000004)
   testDriveLinearDistance(rower, 0.552544989848028)
   testDriveLength(rower, 0.3371976114853044)
-  testDriveAverageHandleForce(rower, 223.750606354492)
-  testDrivePeakHandleForce(rower, 396.7011215854034)
+  testDriveAverageHandleForce(rower, 223.8446015637509)
+  testDrivePeakHandleForce(rower, 396.7733761769528)
   testRecoveryDuration(rower, 0.09847952300000018)
   testRecoveryDragFactor(rower, 309.02744980039836)
-  testInstantHandlePower(rower, 526.5255378417136)
+  testInstantHandlePower(rower, 526.7173432390936)
   // Recovery third stroke starts here
   rower.handleRotationImpulse(0.010769)
   rower.handleRotationImpulse(0.010707554)
@@ -278,8 +279,8 @@ test('Test behaviour for three perfect identical strokes, including settingling
   testDriveDuration(rower, 0.2727572410000002)
   testDriveLinearDistance(rower, 1.2557840678364274)
   testDriveLength(rower, 0.36651914291880905)
-  testDriveAverageHandleForce(rower, 272.7765993429924)
-  testDrivePeakHandleForce(rower, 456.99298986363897)
+  testDriveAverageHandleForce(rower, 272.78784054454604)
+  testDrivePeakHandleForce(rower, 456.96817421200865)
   testRecoveryDuration(rower, 0.09847952300000018)
   testRecoveryDragFactor(rower, 309.02744980039836)
   testInstantHandlePower(rower, 0)
@@ -310,8 +311,8 @@ test('Test behaviour for three perfect identical strokes, including settingling
   testDriveDuration(rower, 0.2727572410000002)
   testDriveLinearDistance(rower, 1.2557840678364274)
   testDriveLength(rower, 0.36651914291880905)
-  testDriveAverageHandleForce(rower, 272.7765993429924)
-  testDrivePeakHandleForce(rower, 456.99298986363897)
+  testDriveAverageHandleForce(rower, 272.78784054454604)
+  testDrivePeakHandleForce(rower, 456.96817421200865)
   testRecoveryDuration(rower, 0.1430115999999999)
   testRecoveryDragFactor(rower, 309.02744980039836)
   testInstantHandlePower(rower, 0)
@@ -379,10 +380,10 @@ test('sample data for NordicTrack RX800 should produce plausible results', async
   await replayRowingSession(rower.handleRotationImpulse, { filename: 'recordings/RX800.csv', realtime: false, loop: false })
 
   testTotalMovingTimeSinceStart(rower, 17.389910236000024)
-  testTotalLinearDistanceSinceStart(rower, 62.49982252262572)
+  testTotalLinearDistanceSinceStart(rower, 62.499750609934196)
   testTotalNumberOfStrokes(rower, 8)
   // As dragFactor is dynamic, it should have changed
-  testRecoveryDragFactor(rower, 493.1277530352103)
+  testRecoveryDragFactor(rower, 493.127960064474)
 })
 
 test('A full session for SportsTech WRX700 should produce plausible results', async () => {
@@ -411,10 +412,10 @@ test('A full session for a Concept2 RowErg should produce plausible results', as
   await replayRowingSession(rower.handleRotationImpulse, { filename: 'recordings/Concept2_RowErg_Session_2000meters.csv', realtime: false, loop: false })
 
   testTotalMovingTimeSinceStart(rower, 590.111937)
-  testTotalLinearDistanceSinceStart(rower, 2029.6932502534587)
+  testTotalLinearDistanceSinceStart(rower, 2029.6932305734617)
   testTotalNumberOfStrokes(rower, 206)
   // As dragFactor isn't static, it should have changed
-  testRecoveryDragFactor(rower, 80.79039510767821)
+  testRecoveryDragFactor(rower, 80.7904433692072)
 })
 
 function testStrokeState (rower, expectedValue) {

app/engine/utils/BinarySearchTree.js

@@ -1,6 +1,6 @@
 'use strict'
 /*
-  Open Rowing Monitor, https://github.com/jaapvanekris/openrowingmonitor
+  Open Rowing Monitor, https://github.com/JaapvanEkris/openrowingmonitor
 
   This creates an ordered series with labels
   It allows for efficient determining the Median, Number of Above and Below

app/engine/utils/BinarySearchTree.test.js

@@ -1,6 +1,6 @@
 'use strict'
 /*
-  Open Rowing Monitor, https://github.com/laberning/openrowingmonitor
+  Open Rowing Monitor, https://github.com/JaapvanEkris/openrowingmonitor
 
   As this object is fundamental for most other utility objects, we must test its behaviour quite thoroughly
 */

app/engine/utils/FullTSLinearSeries.js

@@ -1,6 +1,6 @@
 'use strict'
 /*
-  Open Rowing Monitor, https://github.com/laberning/openrowingmonitor
+  Open Rowing Monitor, https://github.com/JaapvanEkris/openrowingmonitor
 
   The TSLinearSeries is a datatype that represents a Linear Series. It allows
   values to be retrieved (like a FiFo buffer, or Queue) but it also includes
@@ -33,11 +33,11 @@ function createTSLinearSeries (maxSeriesLength = 0) {
   let _B = 0
 
   function push (x, y) {
-    // Invariant: A contains all a's (as in the general formula y = a * x^2 + b * x + c)
+    // Invariant: A contains all a's (as in the general formula y = a * x + b)
     // Where the a's are labeled in the Binary Search Tree with their xi when they BEGIN in the point (xi, yi)
     if (maxSeriesLength > 0 && X.length() >= maxSeriesLength) {
-      // The maximum of the array has been reached, so when pushing the x,y the array gets shifted,
-      // thus we have to remove the a's belonging to the current position X0 as well before this value is trashed
+      // The maximum of the array has been reached, so when pushing the x,y into the arrays, they get shifted automatically,
+      // So, we have to remove the a's belonging to the current position X0 as well before this value is trashed
       A.remove(X.get(0))
     }
 
@@ -49,6 +49,8 @@ function createTSLinearSeries (maxSeriesLength = 0) {
       // There are at least two points in the X and Y arrays, so let's add the new datapoint
       let i = 0
       while (i < X.length() - 1) {
+        // Calculate the slope with all preceeding datapoints with the X.length() - 1'th datapoint (as the array starts at zero)
+        // And store it at its beginpoint (i.e. X.get(i)) to allow remove when that point gets removed from the flank
         A.push(X.get(i), calculateSlope(i, X.length() - 1))
         i++
       }
@@ -66,7 +68,7 @@ function createTSLinearSeries (maxSeriesLength = 0) {
     if (X.length() > 1) {
       // There are at least two points in the X and Y arrays, so let's calculate the intercept
       let i = 0
-      while (i < X.length() - 1) {
+      while (i < X.length()) {
         // Please note , as we need to recreate the B-tree for each newly added datapoint anyway, the label i isn't relevant
         B.push(i, (Y.get(i) - (_A * X.get(i))))
         i++
@@ -101,23 +103,34 @@ function createTSLinearSeries (maxSeriesLength = 0) {
   function goodnessOfFit () {
     // This function returns the R^2 as a goodness of fit indicator
     let i = 0
-    let ssr = 0
+    let sse = 0
     let sst = 0
+    let _goodnessOfFit = 0
     if (X.length() >= 2) {
-      while (i < X.length() - 1) {
-        ssr += Math.pow((Y.get(i) - projectX(X.get(i))), 2)
+      while (i < X.length()) {
+        sse += Math.pow((Y.get(i) - projectX(X.get(i))), 2)
         sst += Math.pow((Y.get(i) - Y.average()), 2)
         i++
       }
-      if (sst !== 0) {
-        const _goodnessOfFit = 1 - (ssr / sst)
-        return _goodnessOfFit
-      } else {
-        return 0
+      switch (true) {
+        case (sse === 0):
+          _goodnessOfFit = 1
+          break
+        case (sse > sst):
+          // This is a pretty bad fit as the error is bigger than just using the line for the average y as intercept
+          _goodnessOfFit = 0
+          break
+        case (sst !== 0):
+          _goodnessOfFit = 1 - (sse / sst)
+          break
+        default:
+          // When SST = 0, R2 isn't defined
+          _goodnessOfFit = 0
       }
     } else {
-      return 0
+      _goodnessOfFit = 0
     }
+    return _goodnessOfFit
   }
 
   function projectX (x) {

app/engine/utils/FullTSLinearSeries.test.js

@@ -1,6 +1,6 @@
 'use strict'
 /*
-  Open Rowing Monitor, https://github.com/laberning/openrowingmonitor
+  Open Rowing Monitor, https://github.com/JaapvanEkris/openrowingmonitor
 */
 import { test } from 'uvu'
 import * as assert from 'uvu/assert'

app/engine/utils/FullTSQuadraticSeries.js

@@ -1,6 +1,6 @@
 'use strict'
 /*
-  Open Rowing Monitor, https://github.com/laberning/openrowingmonitor
+  Open Rowing Monitor, https://github.com/JaapvanEkris/openrowingmonitor
 
   The FullTSQuadraticSeries is a datatype that represents a Quadratic Series. It allows
   values to be retrieved (like a FiFo buffer, or Queue) but it also includes
@@ -69,21 +69,13 @@ function createTSQuadraticSeries (maxSeriesLength = 0) {
 
         // Next, we calculate the B and C via Linear regression over the residu
         i = 0
-        while (i < X.length() - 1) {
+        while (i < X.length()) {
           linearResidu.push(X.get(i), Y.get(i) - (_A * Math.pow(X.get(i), 2)))
           i++
         }
         _B = linearResidu.coefficientA()
         _C = linearResidu.coefficientB()
         break
-      case (X.length() === 2 && X.get(1) - X.get(0) !== 0):
-        // There are only two datapoints, so we need to be creative to get to a quadratic solution
-        // As we know this is part of a 'linear' acceleration, we know that the second derivative should obey 2 * _A = angular acceleration = 2 * angular distance / (delta t)^2
-        _A = (Y.get(1) - Y.get(0)) / Math.pow(X.get(1) - X.get(0), 2)
-        // As the first derivative should match angular velocity (= angular acceleration * (delta t))
-        _B = -2 * _A * X.get(0)
-        _C = 0
-        break
       default:
         _A = 0
         _B = 0
@@ -141,23 +133,34 @@ function createTSQuadraticSeries (maxSeriesLength = 0) {
   function goodnessOfFit () {
     // This function returns the R^2 as a goodness of fit indicator
     let i = 0
-    let ssr = 0
+    let sse = 0
     let sst = 0
-    if (X.length() >= 2) {
-      while (i < X.length() - 1) {
-        ssr += Math.pow((Y.get(i) - projectX(X.get(i))), 2)
+    let _goodnessOfFit = 0
+    if (X.length() > 2) {
+      while (i < X.length()) {
+        sse += Math.pow((Y.get(i) - projectX(X.get(i))), 2)
         sst += Math.pow((Y.get(i) - Y.average()), 2)
         i++
       }
-      if (sst !== 0) {
-        const _goodnessOfFit = 1 - (ssr / sst)
-        return _goodnessOfFit
-      } else {
-        return 0
+      switch (true) {
+        case (sse === 0):
+          _goodnessOfFit = 1
+          break
+        case (sse > sst):
+          // This is a pretty bad fit as the error is bigger than just using the line for the average y as intercept
+          _goodnessOfFit = 0
+          break
+        case (sst !== 0):
+          _goodnessOfFit = 1 - (sse / sst)
+          break
+        default:
+          // When SST = 0, R2 isn't defined
+          _goodnessOfFit = 0
       }
     } else {
-      return 0
+      _goodnessOfFit = 0
     }
+    return _goodnessOfFit
   }
 
   function projectX (x) {