From b11a63a8fdcd3315440b5765948ddd79dd3d6024 Mon Sep 17 00:00:00 2001
From: Jaap van Ekris <82339657+JaapvanEkris@users.noreply.github.com>
Date: Wed, 24 Jan 2024 15:50:09 +0100
Subject: [PATCH] Change to use Lin. TS for Drag calculation

Changing the drag and recovery slope calculation from OLS to Linear Theil-Sen
---
 app/engine/Flywheel.js | 40 +++++++++++++++++++++++-----------------
 1 file changed, 23 insertions(+), 17 deletions(-)

diff --git a/app/engine/Flywheel.js b/app/engine/Flywheel.js
index 5eaab18..e8230d0 100644
--- a/app/engine/Flywheel.js
+++ b/app/engine/Flywheel.js
@@ -23,6 +23,7 @@
 import loglevel from 'loglevel'
 import { createStreamFilter } from './utils/StreamFilter.js'
 import { createOLSLinearSeries } from './utils/OLSLinearSeries.js'
+import { createTSLinearSeries } from './utils/FullTSLinearSeries.js'
 import { createTSQuadraticSeries } from './utils/FullTSQuadraticSeries.js'
 import { createWeighedSeries } from './utils/WeighedSeries.js'
 
@@ -35,13 +36,14 @@ function createFlywheel (rowerSettings) {
   const minumumTorqueBeforeStroke = rowerSettings.minumumForceBeforeStroke * (rowerSettings.sprocketRadius / 100)
   const currentDt = createStreamFilter(rowerSettings.smoothing, rowerSettings.maximumTimeBetweenImpulses)
   const _deltaTime = createOLSLinearSeries(flankLength)
+//  const _deltaTime = createTSLinearSeries(flankLength)
   const _angularDistance = createTSQuadraticSeries(flankLength)
   const _angularVelocityMatrix = []
   const _angularAccelerationMatrix = []
-  const drag = createStreamFilter(rowerSettings.dragFactorSmoothing, (rowerSettings.dragFactor / 1000000))
-  const recoveryDeltaTime = createOLSLinearSeries()
+  const drag = createWeighedSeries(rowerSettings.dragFactorSmoothing, (rowerSettings.dragFactor / 1000000))
+  const recoveryDeltaTime = createTSLinearSeries()
   const strokedetectionMinimalGoodnessOfFit = rowerSettings.minimumStrokeQuality
-  const minumumRecoverySlope = createStreamFilter(rowerSettings.dragFactorSmoothing, rowerSettings.minumumRecoverySlope)
+  const minumumRecoverySlope = createWeighedSeries(rowerSettings.dragFactorSmoothing, rowerSettings.minumumRecoverySlope)
   let _deltaTimeBeforeFlank
   let _angularVelocityAtBeginFlank
   let _angularVelocityBeforeFlank
@@ -117,8 +119,8 @@ function createFlywheel (rowerSettings) {
     }
 
     // Let's make room for a new set of values for angular velocity and acceleration
-    _angularVelocityMatrix[_angularVelocityMatrix.length] = createWeighedSeries()
-    _angularAccelerationMatrix[_angularAccelerationMatrix.length] = createWeighedSeries()
+    _angularVelocityMatrix[_angularVelocityMatrix.length] = createWeighedSeries(flankLength, 0)
+    _angularAccelerationMatrix[_angularAccelerationMatrix.length] = createWeighedSeries(flankLength, 0)
 
     let i = 0
     const goodnessOfFit = _angularDistance.goodnessOfFit()
@@ -133,7 +135,7 @@ function createFlywheel (rowerSettings) {
     _angularAccelerationAtBeginFlank = _angularAccelerationMatrix[0].weighedAverage()
 
     // And finally calculate the torque
-    _torqueAtBeginFlank = (rowerSettings.flywheelInertia * _angularAccelerationAtBeginFlank + drag.clean() * Math.pow(_angularVelocityAtBeginFlank, 2))
+    _torqueAtBeginFlank = (rowerSettings.flywheelInertia * _angularAccelerationAtBeginFlank + drag.weighedAverage() * Math.pow(_angularVelocityAtBeginFlank, 2))
   }
 
   function maintainStateOnly () {
@@ -153,24 +155,28 @@ function createFlywheel (rowerSettings) {
     // Completion of the recovery phase
     inRecoveryPhase = false
 
+    // As goodnessOfFit is calculated in-situ (for 220 datapoints on a C2) and is CPU intensive, we only calculate it only once and reuse the cached value
+    const goodnessOfFit = recoveryDeltaTime.goodnessOfFit()
+
     // Calculation of the drag-factor
-    if (rowerSettings.autoAdjustDragFactor && recoveryDeltaTime.length() > minimumDragFactorSamples && recoveryDeltaTime.slope() > 0 && (!drag.reliable() || recoveryDeltaTime.goodnessOfFit() >= rowerSettings.minimumDragQuality)) {
-      drag.push(slopeToDrag(recoveryDeltaTime.slope()))
-      log.debug(`*** Calculated drag factor: ${(slopeToDrag(recoveryDeltaTime.slope()) * 1000000).toFixed(4)}, no. samples: ${recoveryDeltaTime.length()}, Goodness of Fit: ${recoveryDeltaTime.goodnessOfFit().toFixed(4)}`)
+    if (rowerSettings.autoAdjustDragFactor && recoveryDeltaTime.length() > minimumDragFactorSamples && recoveryDeltaTime.slope() > 0 && (!drag.reliable() || goodnessOfFit >= rowerSettings.minimumDragQuality)) {
+      drag.push(slopeToDrag(recoveryDeltaTime.slope()), goodnessOfFit)
+
+      log.debug(`*** Calculated drag factor: ${(slopeToDrag(recoveryDeltaTime.slope()) * 1000000).toFixed(4)}, no. samples: ${recoveryDeltaTime.length()}, Goodness of Fit: ${goodnessOfFit.toFixed(4)}`)
       if (rowerSettings.autoAdjustRecoverySlope) {
         // We are allowed to autoadjust stroke detection slope as well, so let's do that
-        minumumRecoverySlope.push((1 - rowerSettings.autoAdjustRecoverySlopeMargin) * recoveryDeltaTime.slope())
-        log.debug(`*** Calculated recovery slope: ${recoveryDeltaTime.slope().toFixed(6)}, Goodness of Fit: ${recoveryDeltaTime.goodnessOfFit().toFixed(4)}`)
+        minumumRecoverySlope.push((1 - rowerSettings.autoAdjustRecoverySlopeMargin) * recoveryDeltaTime.slope(), goodnessOfFit)
+        log.debug(`*** Calculated recovery slope: ${recoveryDeltaTime.slope().toFixed(6)}, Goodness of Fit: ${goodnessOfFit.toFixed(4)}`)
       } else {
         // We aren't allowed to adjust the slope, let's report the slope to help help the user configure it
-        log.debug(`*** Calculated recovery slope: ${recoveryDeltaTime.slope().toFixed(6)}, Goodness of Fit: ${recoveryDeltaTime.goodnessOfFit().toFixed(4)}, not used as autoAdjustRecoverySlope isn't set to true`)
+        log.debug(`*** Calculated recovery slope: ${recoveryDeltaTime.slope().toFixed(6)}, Goodness of Fit: ${goodnessOfFit.toFixed(4)}, not used as autoAdjustRecoverySlope isn't set to true`)
       }
     } else {
       if (!rowerSettings.autoAdjustDragFactor) {
         // autoAdjustDampingConstant = false, thus the update is skipped, but let's log the dragfactor anyway
         log.debug(`*** Calculated drag factor: ${(slopeToDrag(recoveryDeltaTime.slope()) * 1000000).toFixed(4)}, slope: ${recoveryDeltaTime.slope().toFixed(8)}, not used because autoAdjustDragFactor is not true`)
       } else {
-        log.debug(`*** Calculated drag factor: ${(slopeToDrag(recoveryDeltaTime.slope()) * 1000000).toFixed(4)}, not used because reliability was too low. no. samples: ${recoveryDeltaTime.length()}, fit: ${recoveryDeltaTime.goodnessOfFit().toFixed(4)}`)
+        log.debug(`*** Calculated drag factor: ${(slopeToDrag(recoveryDeltaTime.slope()) * 1000000).toFixed(4)}, not used because reliability was too low. no. samples: ${recoveryDeltaTime.length()}, fit: ${goodnessOfFit.toFixed(4)}`)
       }
     }
   }
@@ -217,7 +223,7 @@ function createFlywheel (rowerSettings) {
 
   function dragFactor () {
     // Ths function returns the current dragfactor of the flywheel
-    return drag.clean()
+    return drag.weighedAverage()
   }
 
   function isDwelling () {
@@ -243,9 +249,9 @@ function createFlywheel (rowerSettings) {
   }
 
   function isUnpowered () {
-    if (deltaTimeSlopeAbove(minumumRecoverySlope.clean()) && torqueAbsent() && _deltaTime.length() >= flankLength) {
+    if (deltaTimeSlopeAbove(minumumRecoverySlope.weighedAverage()) && torqueAbsent() && _deltaTime.length() >= flankLength) {
       // We reached the minimum number of increasing currentDt values
-      // log.info(`*** INFO: recovery detected based on due to slope exceeding recoveryslope = ${deltaTimeSlopeAbove(minumumRecoverySlope.clean())}, exceeding minumumForceBeforeStroke = ${torqueAbsent()}`)
+      // log.info(`*** INFO: recovery detected based on due to slope exceeding recoveryslope = ${deltaTimeSlopeAbove(minumumRecoverySlope.weighedAverage())}, exceeding minumumForceBeforeStroke = ${torqueAbsent()}`)
       return true
     } else {
       return false
@@ -253,7 +259,7 @@ function createFlywheel (rowerSettings) {
   }
 
   function isPowered () {
-    if ((deltaTimeSlopeBelow(minumumRecoverySlope.clean()) && torquePresent()) || _deltaTime.length() < flankLength) {
+    if ((deltaTimeSlopeBelow(minumumRecoverySlope.weighedAverage()) && torquePresent()) || _deltaTime.length() < flankLength) {
       return true
     } else {
       return false