[ML] Correct ageing of windowed seasonal components (#88)

In particular, we were ageing them too fast and their decay rate should be prorated by the 
fraction of values with which they are updated.

Author: tveasey

docs/CHANGELOG.asciidoc

@@ -22,6 +22,8 @@ model state on disk ({pull}89[#89])
 
 === Bug Fixes
 
+Age seasonal components in proportion to the fraction of values with which they're updated ({pull}88[#88])
+
 === Regressions
 
 === Known Issues

lib/maths/CSeasonalComponentAdaptiveBucketing.cc

@@ -240,7 +240,8 @@ void CSeasonalComponentAdaptiveBucketing::propagateForwardsByTime(double time, b
     if (time < 0.0) {
         LOG_ERROR(<< "Can't propagate bucketing backwards in time");
     } else if (this->initialized()) {
-        double factor{std::exp(-this->CAdaptiveBucketing::decayRate() * time)};
+        double factor{std::exp(-this->CAdaptiveBucketing::decayRate() *
+                               m_Time->fractionInWindow() * time)};
         this->CAdaptiveBucketing::age(factor);
         for (auto& bucket : m_Buckets) {
             bucket.s_Regression.age(factor, meanRevert);

lib/maths/CTimeSeriesDecompositionDetail.cc

@@ -1103,9 +1103,8 @@ void CTimeSeriesDecompositionDetail::CComponents::handle(const SAddValue& messag
         for (std::size_t i = 1u; i <= m; ++i) {
             CSeasonalComponent* component{seasonalComponents[i - 1]};
             CComponentErrors* error_{seasonalErrors[i - 1]};
-            double wi{weight / component->time().fractionInWindow()};
-            component->add(time, values[i], wi);
-            error_->add(error, predictions[i - 1], wi);
+            component->add(time, values[i], weight);
+            error_->add(error, predictions[i - 1], weight);
         }
         for (std::size_t i = m + 1; i <= m + n; ++i) {
             CCalendarComponent* component{calendarComponents[i - m - 1]};

lib/maths/unittest/CForecastTest.cc

@@ -83,14 +83,14 @@ void CForecastTest::testDailyNoLongTermTrend() {
 
     test::CRandomNumbers rng;
 
-    auto trend = [&y, bucketLength](core_t::TTime time, double noise) {
-        core_t::TTime i{(time % 86400) / bucketLength};
+    TTrend trend = [&y, bucketLength](core_t::TTime time, double noise) {
+        core_t::TTime i{(time % core::constants::DAY) / bucketLength};
         double alpha{static_cast<double>(i % 6) / 6.0};
         double beta{1.0 - alpha};
         return 40.0 + alpha * y[i / 6] + beta * y[(i / 6 + 1) % y.size()] + noise;
     };
 
-    this->test(trend, bucketLength, 60, 64.0, 5.0, 0.14);
+    this->test(trend, bucketLength, 63, 64.0, 5.0, 0.14);
 }
 
 void CForecastTest::testDailyConstantLongTermTrend() {
@@ -99,18 +99,18 @@ void CForecastTest::testDailyConstantLongTermTrend() {
                  80.0, 100.0, 110.0, 120.0, 110.0, 100.0, 90.0, 80.0,
                  30.0, 15.0,  10.0,  8.0,   5.0,   3.0,   2.0,  0.0};
 
-    auto trend = [&y, bucketLength](core_t::TTime time, double noise) {
-        core_t::TTime i{(time % 86400) / bucketLength};
+    TTrend trend = [&y, bucketLength](core_t::TTime time, double noise) {
+        core_t::TTime i{(time % core::constants::DAY) / bucketLength};
         return 0.25 * static_cast<double>(time) / static_cast<double>(bucketLength) +
                y[i] + noise;
     };
 
-    this->test(trend, bucketLength, 60, 64.0, 16.0, 0.02);
+    this->test(trend, bucketLength, 63, 64.0, 14.0, 0.02);
 }
 
 void CForecastTest::testDailyVaryingLongTermTrend() {
     core_t::TTime bucketLength{3600};
-    double day{86400.0};
+    double day{static_cast<double>(core::constants::DAY)};
     TDoubleVec times{0.0,         5.0 * day,   10.0 * day,  15.0 * day,
                      20.0 * day,  25.0 * day,  30.0 * day,  35.0 * day,
                      40.0 * day,  45.0 * day,  50.0 * day,  55.0 * day,
@@ -124,13 +124,13 @@ void CForecastTest::testDailyVaryingLongTermTrend() {
     maths::CSpline<> trend_(maths::CSplineTypes::E_Cubic);
     trend_.interpolate(times, values, maths::CSplineTypes::E_Natural);
 
-    auto trend = [&trend_](core_t::TTime time, double noise) {
+    TTrend trend = [&trend_](core_t::TTime time, double noise) {
         double time_{static_cast<double>(time)};
         return trend_.value(time_) +
                8.0 * std::sin(boost::math::double_constants::two_pi * time_ / 43200.0) + noise;
     };
 
-    this->test(trend, bucketLength, 100, 9.0, 13.0, 0.04);
+    this->test(trend, bucketLength, 98, 9.0, 14.0, 0.042);
 }
 
 void CForecastTest::testComplexNoLongTermTrend() {
@@ -140,13 +140,13 @@ void CForecastTest::testComplexNoLongTermTrend() {
                  60.0, 40.0,  30.0,  20.0,  10.0,  10.0,  5.0,  0.0};
     TDoubleVec scale{1.0, 1.1, 1.05, 0.95, 0.9, 0.3, 0.2};
 
-    auto trend = [&y, &scale, bucketLength](core_t::TTime time, double noise) {
-        core_t::TTime d{(time % 604800) / 86400};
-        core_t::TTime h{(time % 86400) / bucketLength};
+    TTrend trend = [&y, &scale, bucketLength](core_t::TTime time, double noise) {
+        core_t::TTime d{(time % core::constants::WEEK) / core::constants::DAY};
+        core_t::TTime h{(time % core::constants::DAY) / bucketLength};
         return scale[d] * (20.0 + y[h] + noise);
     };
 
-    this->test(trend, bucketLength, 60, 24.0, 28.0, 0.13);
+    this->test(trend, bucketLength, 63, 24.0, 8.0, 0.13);
 }
 
 void CForecastTest::testComplexConstantLongTermTrend() {
@@ -156,19 +156,19 @@ void CForecastTest::testComplexConstantLongTermTrend() {
                  60.0, 40.0,  30.0,  20.0,  10.0,  10.0,  5.0,  0.0};
     TDoubleVec scale{1.0, 1.1, 1.05, 0.95, 0.9, 0.3, 0.2};
 
-    auto trend = [&y, &scale, bucketLength](core_t::TTime time, double noise) {
-        core_t::TTime d{(time % 604800) / 86400};
-        core_t::TTime h{(time % 86400) / bucketLength};
+    TTrend trend = [&y, &scale, bucketLength](core_t::TTime time, double noise) {
+        core_t::TTime d{(time % core::constants::WEEK) / core::constants::DAY};
+        core_t::TTime h{(time % core::constants::DAY) / bucketLength};
         return 0.25 * static_cast<double>(time) / static_cast<double>(bucketLength) +
                scale[d] * (20.0 + y[h] + noise);
     };
 
-    this->test(trend, bucketLength, 60, 24.0, 14.0, 0.01);
+    this->test(trend, bucketLength, 63, 24.0, 8.0, 0.01);
 }
 
 void CForecastTest::testComplexVaryingLongTermTrend() {
     core_t::TTime bucketLength{3600};
-    double day{86400.0};
+    double day{static_cast<double>(core::constants::DAY)};
     TDoubleVec times{0.0,         5.0 * day,   10.0 * day,  15.0 * day,
                      20.0 * day,  25.0 * day,  30.0 * day,  35.0 * day,
                      40.0 * day,  45.0 * day,  50.0 * day,  55.0 * day,
@@ -186,14 +186,14 @@ void CForecastTest::testComplexVaryingLongTermTrend() {
     maths::CSpline<> trend_(maths::CSplineTypes::E_Cubic);
     trend_.interpolate(times, values, maths::CSplineTypes::E_Natural);
 
-    auto trend = [&trend_, &y, &scale, bucketLength](core_t::TTime time, double noise) {
-        core_t::TTime d{(time % 604800) / 86400};
-        core_t::TTime h{(time % 86400) / bucketLength};
+    TTrend trend = [&trend_, &y, &scale, bucketLength](core_t::TTime time, double noise) {
+        core_t::TTime d{(time % core::constants::WEEK) / core::constants::DAY};
+        core_t::TTime h{(time % core::constants::DAY) / bucketLength};
         double time_{static_cast<double>(time)};
         return trend_.value(time_) + scale[d] * (20.0 + y[h] + noise);
     };
 
-    this->test(trend, bucketLength, 60, 4.0, 28.0, 0.053);
+    this->test(trend, bucketLength, 63, 4.0, 42.0, 0.06);
 }
 
 void CForecastTest::testNonNegative() {
@@ -400,7 +400,6 @@ void CForecastTest::test(TTrend trend,
                          double noiseVariance,
                          double maximumPercentageOutOfBounds,
                          double maximumError) {
-
     //std::ofstream file;
     //file.open("results.m");
     //TDoubleVec actual;
@@ -423,7 +422,8 @@ void CForecastTest::test(TTrend trend,
     TDouble2VecWeightsAryVec weights{maths_t::CUnitWeights::unit<TDouble2Vec>(1)};
     for (std::size_t d = 0u; d < daysToLearn; ++d) {
         TDoubleVec noise;
-        rng.generateNormalSamples(0.0, noiseVariance, 86400 / bucketLength, noise);
+        rng.generateNormalSamples(0.0, noiseVariance,
+                                  core::constants::DAY / bucketLength, noise);
 
         for (std::size_t i = 0u; i < noise.size(); ++i, time += bucketLength) {
             maths::CModelAddSamplesParams params;
@@ -450,7 +450,8 @@ void CForecastTest::test(TTrend trend,
 
     for (std::size_t i = 0u; i < prediction.size(); /**/) {
         TDoubleVec noise;
-        rng.generateNormalSamples(0.0, noiseVariance, 86400 / bucketLength, noise);
+        rng.generateNormalSamples(0.0, noiseVariance,
+                                  core::constants::DAY / bucketLength, noise);
         TDoubleVec day;
         for (std::size_t j = 0u; i < prediction.size() && j < noise.size();
              ++i, ++j, time += bucketLength) {

lib/maths/unittest/CSeasonalComponentAdaptiveBucketingTest.cc

@@ -11,6 +11,7 @@
 #include <core/CRapidXmlParser.h>
 #include <core/CRapidXmlStatePersistInserter.h>
 #include <core/CRapidXmlStateRestoreTraverser.h>
+#include <core/Constants.h>
 
 #include <maths/CSeasonalComponentAdaptiveBucketing.h>
 #include <maths/CSeasonalTime.h>
@@ -36,14 +37,14 @@ using TMaxAccumulator = maths::CBasicStatistics::SMax<double>::TAccumulator;
 }
 
 void CSeasonalComponentAdaptiveBucketingTest::testInitialize() {
-    maths::CDiurnalTime time(0, 1, 101, 100);
+    maths::CGeneralPeriodTime time(100);
     maths::CSeasonalComponentAdaptiveBucketing bucketing(time);
 
     CPPUNIT_ASSERT(!bucketing.initialize(0));
 
     const std::string expectedEndpoints("[0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100]");
-    const std::string expectedKnots("[0, 4, 14, 24, 34, 44, 54, 64, 74, 84, 94, 100]");
-    const std::string expectedValues("[41, 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 41]");
+    const std::string expectedKnots("[0, 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 100]");
+    const std::string expectedValues("[50, 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 50]");
 
     CPPUNIT_ASSERT(bucketing.initialize(10));
     const TFloatVec& endpoints = bucketing.endpoints();
@@ -62,7 +63,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testInitialize() {
 }
 
 void CSeasonalComponentAdaptiveBucketingTest::testSwap() {
-    maths::CDiurnalTime time1(0, 0, 100, 100);
+    maths::CGeneralPeriodTime time1(100);
     maths::CSeasonalComponentAdaptiveBucketing bucketing1(time1, 0.05);
 
     test::CRandomNumbers rng;
@@ -82,7 +83,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testSwap() {
         bucketing1.propagateForwardsByTime(1.0);
     }
 
-    maths::CDiurnalTime time2(10, 10, 120, 110);
+    maths::CGeneralPeriodTime time2(120);
     maths::CSeasonalComponentAdaptiveBucketing bucketing2(time2, 0.1);
 
     uint64_t checksum1 = bucketing1.checksum();
@@ -109,7 +110,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testRefine() {
         double function[] = {10, 10,  10, 10, 100, 90, 80,
                              90, 100, 20, 10, 10,  10, 10};
 
-        maths::CDiurnalTime time(0, 0, 86400, 86400);
+        maths::CDiurnalTime time(0, 0, core::constants::WEEK, core::constants::DAY);
         maths::CSeasonalComponentAdaptiveBucketing bucketing1(time);
         maths::CSeasonalComponentAdaptiveBucketing bucketing2(time);
 
@@ -188,7 +189,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testRefine() {
         // Test that the variance in each bucket is approximately equal.
         // The test function is y = (x - 50)^2 / 50.
 
-        maths::CDiurnalTime time(0, 0, 100, 100);
+        maths::CGeneralPeriodTime time(100);
         maths::CSeasonalComponentAdaptiveBucketing bucketing(time, 0.05);
 
         bucketing.initialize(10);
@@ -294,7 +295,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testPropagateForwardsByTime() {
     // the bucket values and that the rate at which the total
     // count is reduced uniformly.
 
-    maths::CDiurnalTime time(0, 0, 100, 100);
+    maths::CGeneralPeriodTime time(100);
     maths::CSeasonalComponentAdaptiveBucketing bucketing(time, 0.2);
 
     bucketing.initialize(10);
@@ -331,7 +332,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testMinimumBucketLength() {
 
     core_t::TTime period = static_cast<core_t::TTime>(n) *
                            static_cast<core_t::TTime>(bucketLength);
-    maths::CDiurnalTime time(0, 0, period, period);
+    maths::CGeneralPeriodTime time(period);
     maths::CSeasonalComponentAdaptiveBucketing bucketing1(time, 0.0, 0.0);
     maths::CSeasonalComponentAdaptiveBucketing bucketing2(time, 0.0, 3000.0);
     bucketing1.initialize(n);
@@ -389,7 +390,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testUnintialized() {
     // Check that all the functions work and return the expected
     // values on an uninitialized bucketing.
 
-    maths::CDiurnalTime time(0, 0, 10, 10);
+    maths::CGeneralPeriodTime time(10);
     maths::CSeasonalComponentAdaptiveBucketing bucketing(time, 0.1);
 
     bucketing.add(0, 1.0, 1.0);
@@ -436,7 +437,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testKnots() {
 
     LOG_DEBUG(<< "*** Values ***");
     {
-        maths::CDiurnalTime time(0, 0, 86400, 86400);
+        maths::CDiurnalTime time(0, 0, core::constants::WEEK, core::constants::DAY);
         maths::CSeasonalComponentAdaptiveBucketing bucketing(time, 0.1, 864.0);
 
         bucketing.initialize(20);
@@ -479,7 +480,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testKnots() {
     }
     LOG_DEBUG(<< "*** Variances ***");
     {
-        maths::CDiurnalTime time(0, 0, 86400, 86400);
+        maths::CDiurnalTime time(0, 0, core::constants::WEEK, core::constants::DAY);
         maths::CSeasonalComponentAdaptiveBucketing bucketing(time, 0.1, 864.0);
 
         bucketing.initialize(20);
@@ -530,7 +531,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testLongTermTrendKnots() {
 
     test::CRandomNumbers rng;
 
-    maths::CDiurnalTime time(0, 0, 86400, 86400);
+    maths::CDiurnalTime time(0, 0, core::constants::WEEK, core::constants::DAY);
     maths::CSeasonalComponentAdaptiveBucketing bucketing(time, 0.1, 864.0);
     maths::CSeasonalComponentAdaptiveBucketing::TFloatMeanAccumulatorVec empty;
 
@@ -588,7 +589,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testShiftValue() {
     // Test that applying a shift translates the predicted values
     // but doesn't alter the slope or predicted variances.
 
-    maths::CDiurnalTime time(0, 0, 86400, 86400);
+    maths::CDiurnalTime time(0, 0, core::constants::WEEK, core::constants::DAY);
     maths::CSeasonalComponentAdaptiveBucketing bucketing(time, 0.1, 600.0);
     maths::CSeasonalComponentAdaptiveBucketing::TFloatMeanAccumulatorVec empty;
 
@@ -632,7 +633,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testShiftValue() {
 void CSeasonalComponentAdaptiveBucketingTest::testSlope() {
     // Test that the slope increases by the shift.
 
-    maths::CDiurnalTime time(0, 0, 86400, 86400);
+    maths::CDiurnalTime time(0, 0, core::constants::WEEK, core::constants::DAY);
     maths::CSeasonalComponentAdaptiveBucketing bucketing(time, 0.1, 600.0);
     maths::CSeasonalComponentAdaptiveBucketing::TFloatMeanAccumulatorVec empty;
 
@@ -669,7 +670,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testPersist() {
     double decayRate = 0.1;
     double minimumBucketLength = 1.0;
 
-    maths::CDiurnalTime time(0, 0, 86400, 86400);
+    maths::CDiurnalTime time(0, 0, core::constants::WEEK, core::constants::DAY);
     maths::CSeasonalComponentAdaptiveBucketing origBucketing(time, decayRate, minimumBucketLength);
 
     origBucketing.initialize(10);
@@ -723,7 +724,7 @@ void CSeasonalComponentAdaptiveBucketingTest::testUpgrade() {
     double decayRate = 0.1;
     double minimumBucketLength = 1.0;
 
-    maths::CDiurnalTime time(0, 0, 86400, 86400);
+    maths::CDiurnalTime time(0, 0, core::constants::WEEK, core::constants::DAY);
     maths::CSeasonalComponentAdaptiveBucketing expectedBucketing(time, decayRate, minimumBucketLength);
 
     expectedBucketing.initialize(10);