BikeScore, xPower, and Relative Intensity updated to match Skiba's method

more or less exactly.  Also added (TM) to BikeScore in Ride Summary.

src/gui/BikeScore.cpp

@@ -3,60 +3,72 @@
 #include "Zones.h"
 #include <math.h>
 
+// NOTE: This code follows the description of xPower, Relative Intensity, and
+// BikeScore in "Analysis of Power Output and Training Stress in Cyclists: The
+// Development of the BikeScore(TM) Algorithm", page 5, by Phil Skiba:
+//
+// http://www.physfarm.com/Analysis%20of%20Power%20Output%20and%20Training%20Stress%20in%20Cyclists-%20BikeScore.pdf
+//
+// The weighting factors for the exponentially weighted average are taken from
+// a spreadsheet provided by Dr. Skiba.
+
 class XPower : public RideMetric {
     double xpower;
     double secs;
 
     friend class RelativeIntensity;
+    friend class BikeScore;
 
     public:
 
-    XPower() : xpower(0.0), secs(0) {}
+    XPower() : xpower(0.0), secs(0.0) {}
     QString name() const { return "skiba_xpower"; }
     QString units(bool) const { return "watts"; }
     double value(bool) const { return xpower; }
     void compute(const RawFile *raw, const Zones *, int,
                  const QHash<QString,RideMetric*> &) {
+
+        static const double EPSILON = 0.1;
+        static const double NEGLIGIBLE = 0.1;
+
+        double secsDelta = raw->rec_int_ms / 1000.0;
+        double sampsPerWindow = 25.0 / secsDelta;
+        double attenuation = sampsPerWindow / (sampsPerWindow + secsDelta);
+        double sampleWeight = secsDelta / (sampsPerWindow + secsDelta);
+        
+        double lastSecs = 0.0;
+        double weighted = 0.0;
+
         double total = 0.0;
         int count = 0;
-        // we're interested in a 25 second moving average
-        int movingAvgLen = (int) ceil(25000 / raw->rec_int_ms);
-        double *movingAvgWatts = new double[movingAvgLen];
-        int j = 0; // an index for the moving average function
+
         QListIterator<RawFilePoint*> i(raw->points);
-        double secsDelta = raw->rec_int_ms / 1000.0;
         while (i.hasNext()) {
             const RawFilePoint *point = i.next();
-            if (point->watts >= 0.0) {
-                secs += secsDelta;
-                // 25 second moving average
-                // populate the array in the first 25 seconds
-                if (secs <= 25)
-                    movingAvgWatts[j] = point->watts;
-                else {
-                    // add the most recent value, and calculate average
-                    movingAvgWatts[j] = point->watts;
-                    double avg = 0.0;
-                    for (int jj = 0; jj < movingAvgLen; jj++)
-                        avg += movingAvgWatts[jj];
-                    total += pow(avg / movingAvgLen, 4.0);
-                    count++;
-                }
-                j = (j + 1) % movingAvgLen;
+            while ((weighted > NEGLIGIBLE) 
+                   && (point->secs > lastSecs + secsDelta + EPSILON)) {
+                weighted *= attenuation;
+                lastSecs += secsDelta;
             }
+            weighted *= attenuation;
+            weighted += sampleWeight * point->watts;
+            lastSecs = point->secs;
+            total += pow(weighted, 4.0);
+            count++;
         }
         xpower = pow(total / count, 0.25);
-        delete [] movingAvgWatts;
+        secs = count * secsDelta;
     }
     RideMetric *clone() const { return new XPower(*this); }
 };
 
 class RelativeIntensity : public RideMetric {
     double reli;
+    double secs;
 
     public:
 
-    RelativeIntensity() : reli(0.0) {}
+    RelativeIntensity() : reli(0.0), secs(0.0) {}
     QString name() const { return "skiba_relative_intensity"; }
     QString units(bool) const { return ""; }
     double value(bool) const { return reli; }
@@ -64,9 +76,10 @@ class RelativeIntensity : public RideMetric {
                  const QHash<QString,RideMetric*> &deps) {
         if (zones) {
             assert(deps.contains("skiba_xpower"));
-            RideMetric *xp = deps.value("skiba_xpower");
+            XPower *xp = dynamic_cast<XPower*>(deps.value("skiba_xpower"));
             assert(xp);
-            reli = xp->value(true) / zones->getFTP(zoneRange);
+            reli = xp->xpower / zones->getFTP(zoneRange);
+            secs = xp->secs;
         }
     }
     RideMetric *clone() const { return new RelativeIntensity(*this); }
@@ -81,15 +94,22 @@ class BikeScore : public RideMetric {
     QString name() const { return "skiba_bike_score"; }
     QString units(bool) const { return ""; }
     double value(bool) const { return score; }
-    void compute(const RawFile *raw, const Zones *, int,
+    void compute(const RawFile *, const Zones *zones, int zoneRange,
                  const QHash<QString,RideMetric*> &deps) {
+        assert(deps.contains("skiba_xpower"));
         assert(deps.contains("skiba_relative_intensity"));
+        XPower *xp = dynamic_cast<XPower*>(deps.value("skiba_xpower"));
         RideMetric *ri = deps.value("skiba_relative_intensity");
         assert(ri);
-        score = ((raw->points.back()->secs / 60) / 60) 
-                * (pow(ri->value(true), 2)) * 100.0;
-     }
+        double normWork = xp->xpower * xp->secs;
+        double rawBikeScore = normWork * ri->value(true);
+        // TODO: use CP, not FTP here
+        double workInAnHourAtCP = zones->getFTP(zoneRange) * 3600; 
+        score = rawBikeScore / workInAnHourAtCP * 100.0;
+    }
     RideMetric *clone() const { return new BikeScore(*this); }
+    bool canAggregate() const { return true; }
+    void aggregateWith(RideMetric *other) { score += other->value(true); }
 };
 
 static bool addAllThree() {
@@ -97,7 +117,6 @@ static bool addAllThree() {
     QVector<QString> deps;
     deps.append("skiba_xpower");
     RideMetricFactory::instance().addMetric(RelativeIntensity(), &deps);
-    deps.clear();
     deps.append("skiba_relative_intensity");
     RideMetricFactory::instance().addMetric(BikeScore(), &deps);
     return true;

src/gui/RideItem.cpp

@@ -155,13 +155,14 @@ static const char *metricsXml =
     "    <metric name=\"average_cad\" display_name=\"Cadence\"\n"
     "            precision=\"0\"/>\n"
     "  </metric_group>\n"
-    "  <metric_group name=\"Dr. Skiba's BikeScore\">\n"
+    "  <metric_group name=\"BikeScore(TM)\" note=\"BikeScore is a trademark "
+    "      of Dr. Philip Friere Skiba, PhysFarm Training Systems LLC\">\n"
     "    <metric name=\"skiba_xpower\" display_name=\"xPower\"\n"
     "            precision=\"0\"/>\n"
     "    <metric name=\"skiba_relative_intensity\"\n"
     "            display_name=\"Relative Intensity\" precision=\"3\"/>\n"
     "    <metric name=\"skiba_bike_score\" display_name=\"BikeScore\"\n"
-    "            precision=\"2\"/>\n"
+    "            precision=\"0\"/>\n"
     "  </metric_group>\n"
     "</metrics>\n";
 
@@ -303,17 +304,27 @@ later:
             }
         }
 
+        QString noteString = "";
+        QString stars;
         QDomNodeList groups = doc.elementsByTagName("metric_group");
         for (int groupNum = 0; groupNum < groups.size(); ++groupNum) {
             QDomElement group = groups.at(groupNum).toElement();
             assert(!group.isNull());
             QString groupName = group.attribute("name");
+            QString groupNote = group.attribute("note");
             assert(groupName.length() > 0);
             if (groupNum % 2 == 0)
                 summary += "<table border=0 cellspacing=10><tr>";
             summary += "<td align=\"center\" width=\"45%\"><table>"
                 "<tr><td align=\"center\" colspan=2><h2>%1</h2></td></tr>";
-            summary = summary.arg(groupName);
+            if (groupNote.length() > 0) {
+                stars += "*";
+                summary = summary.arg(groupName + stars);
+                noteString += "<br>" + stars + " " + groupNote;
+            }
+            else {
+                summary = summary.arg(groupName);
+            }
             QDomNodeList metricsList = group.childNodes();
             for (int i = 0; i < metricsList.size(); ++i) {
                 QDomElement metric = metricsList.at(i).toElement();
@@ -399,6 +410,8 @@ later:
                 summary += " <li>" + i.next();
             summary += "</ul>";
         }
+        if (noteString.length() > 0)
+            summary += "<br><hr width=\"80%\">" + noteString;
         summary += "</center>";
     }