[ML] Handle NaNs when detrending seasonal components (#408)

While it is as yet unknown how it can occur, NaN values have been
observed to be present in seasonal components when detrending, resulting
in a chain of errors in the logs as the NaN values percolate up the call
stack.

This PR attempts to detect the NaN values as close as possible to their
source and ensure that they do not unduly affect the analysis.

The approach is to check for bad values in seasonal and calendar
components from within
CTimeSeriesDecompositionDetail::CComponents::interpolate, which is
called periodically. If invalid values are detected in a component then
it is removed from the decomposition.

Author: edsavage

docs/CHANGELOG.asciidoc

@@ -39,9 +39,11 @@
 
 == {es} version 7.1.0
 
-Remove hard limit for maximum forecast interval and limit based on the time interval of data added
+* Remove hard limit for maximum forecast interval and limit based on the time interval of data added
 to the model. (See {pull}214[#214].)
 
+* Handle NaNs when detrending seasonal components. {ml-pull}408[#408]
+
 == {es} version 7.0.0-alpha1
 
 == {es} version 6.7.0

include/maths/CAdaptiveBucketing.h

@@ -133,6 +133,9 @@ class MATHS_EXPORT CAdaptiveBucketing {
     //! Name of component
     virtual std::string name() const = 0;
 
+    //! Check that the state is valid.
+    virtual bool isBad() const = 0;
+
 protected:
     using TRestoreFunc = std::function<bool(core::CStateRestoreTraverser&)>;
     using TPersistFunc = std::function<void(core::CStatePersistInserter&)>;

include/maths/CCalendarComponent.h

@@ -140,6 +140,9 @@ class MATHS_EXPORT CCalendarComponent : private CDecompositionComponent {
     //! Get the memory used by this component.
     std::size_t memoryUsage() const;
 
+    //! Check that the state is valid.
+    bool isBad() const;
+
 private:
     //! Create by traversing a state document.
     bool acceptRestoreTraverser(double decayRate,

include/maths/CCalendarComponentAdaptiveBucketing.h

@@ -100,6 +100,9 @@ class MATHS_EXPORT CCalendarComponentAdaptiveBucketing final : public CAdaptiveB
     //! Name of component
     std::string name() const override;
 
+    //! Check that the state is valid.
+    bool isBad() const override;
+
 private:
     using TFloatMeanVarVec = std::vector<TFloatMeanVarAccumulator>;
 

include/maths/CMathsFuncsDetail.h

@@ -67,16 +67,6 @@ bool CMathsFuncs::isNan(const CSymmetricMatrixNxN<double, N>& val) {
     return anElement(static_cast<bool (*)(double)>(&isNan), val);
 }
 
-template<std::size_t N>
-bool CMathsFuncs::isNan(const core::CSmallVector<double, N>& val) {
-    for (std::size_t i = 0u; i < val.size(); ++i) {
-        if (isNan(val[i])) {
-            return true;
-        }
-    }
-    return false;
-}
-
 template<std::size_t N>
 bool CMathsFuncs::isInf(const CVectorNx1<double, N>& val) {
     return aComponent(static_cast<bool (*)(double)>(&isInf), val);
@@ -87,16 +77,6 @@ bool CMathsFuncs::isInf(const CSymmetricMatrixNxN<double, N>& val) {
     return anElement(static_cast<bool (*)(double)>(&isInf), val);
 }
 
-template<std::size_t N>
-bool CMathsFuncs::isInf(const core::CSmallVector<double, N>& val) {
-    for (std::size_t i = 0u; i < val.size(); ++i) {
-        if (isInf(val[i])) {
-            return true;
-        }
-    }
-    return false;
-}
-
 template<std::size_t N>
 bool CMathsFuncs::isFinite(const CVectorNx1<double, N>& val) {
     return everyComponent(static_cast<bool (*)(double)>(&isFinite), val);
@@ -106,16 +86,6 @@ template<std::size_t N>
 bool CMathsFuncs::isFinite(const CSymmetricMatrixNxN<double, N>& val) {
     return everyElement(static_cast<bool (*)(double)>(&isFinite), val);
 }
-
-template<std::size_t N>
-bool CMathsFuncs::isFinite(const core::CSmallVector<double, N>& val) {
-    for (std::size_t i = 0u; i < val.size(); ++i) {
-        if (!isFinite(val[i])) {
-            return false;
-        }
-    }
-    return true;
-}
 }
 }
 

include/maths/CRegression.h

@@ -22,6 +22,8 @@
 #include <cstddef>
 #include <cstdint>
 
+class CNanInjector;
+
 namespace ml {
 namespace core {
 class CStatePersistInserter;
@@ -400,6 +402,9 @@ class MATHS_EXPORT CRegression {
         //! regression. There are 3N - 1 in total, for the distinct
         //! values in the design matrix and vector.
         TVectorMeanAccumulator m_S;
+
+        //! Befriend a helper class used by the unit tests
+        friend class ::CNanInjector;
     };
 
     //! Get the predicted value of \p r at \p x.

include/maths/CSeasonalComponent.h

@@ -19,6 +19,8 @@
 #include <string>
 #include <vector>
 
+class CNanInjector;
+
 namespace ml {
 namespace core {
 class CStatePersistInserter;
@@ -186,6 +188,9 @@ class MATHS_EXPORT CSeasonalComponent : private CDecompositionComponent {
     //! Get the memory used by this component.
     std::size_t memoryUsage() const;
 
+    //! Check that the state is valid.
+    bool isBad() const { return m_Bucketing.isBad(); }
+
 private:
     //! Create by traversing a state document.
     bool acceptRestoreTraverser(double decayRate,
@@ -202,6 +207,9 @@ class MATHS_EXPORT CSeasonalComponent : private CDecompositionComponent {
 
     //! Regression models for a collection of buckets covering the period.
     CSeasonalComponentAdaptiveBucketing m_Bucketing;
+
+    //! Befriend a helper class used by the unit tests
+    friend class ::CNanInjector;
 };
 
 //! Create a free function which will be picked up in Koenig lookup.

include/maths/CSeasonalComponentAdaptiveBucketing.h

@@ -21,6 +21,8 @@
 
 #include <stdint.h>
 
+class CNanInjector;
+
 namespace ml {
 namespace core {
 class CStatePersistInserter;
@@ -133,6 +135,9 @@ class MATHS_EXPORT CSeasonalComponentAdaptiveBucketing final : public CAdaptiveB
     //! Name of component
     std::string name() const override;
 
+    //! Check that the state is valid.
+    bool isBad() const override;
+
 private:
     using TSeasonalTimePtr = std::unique_ptr<CSeasonalTime>;
 
@@ -149,6 +154,9 @@ class MATHS_EXPORT CSeasonalComponentAdaptiveBucketing final : public CAdaptiveB
 
         uint64_t checksum(uint64_t seed) const;
 
+        //! Check that the state is valid.
+        bool isBad() const;
+
         TRegression s_Regression;
         CFloatStorage s_Variance;
         core_t::TTime s_FirstUpdate;
@@ -198,6 +206,9 @@ class MATHS_EXPORT CSeasonalComponentAdaptiveBucketing final : public CAdaptiveB
 
     //! The buckets.
     TBucketVec m_Buckets;
+
+    //! Befriend a helper class used by the unit tests
+    friend class ::CNanInjector;
 };
 
 //! Create a free function which will be found by Koenig lookup.

include/maths/CTimeSeriesDecomposition.h

@@ -14,7 +14,7 @@
 
 #include <memory>
 
-class CTimeSeriesDecompositionTest;
+class CNanInjector;
 
 namespace ml {
 namespace core {
@@ -252,6 +252,9 @@ class MATHS_EXPORT CTimeSeriesDecomposition : public CTimeSeriesDecompositionInt
 
     //! The state for modeling the components of the decomposition.
     CComponents m_Components;
+
+    //! Befriend a helper class used by the unit tests
+    friend class ::CNanInjector;
 };
 }
 }

include/maths/CTimeSeriesDecompositionDetail.h

@@ -28,6 +28,8 @@
 #include <memory>
 #include <vector>
 
+class CNanInjector;
+
 namespace ml {
 namespace maths {
 class CExpandingWindow;
@@ -373,12 +375,12 @@ class MATHS_EXPORT CTimeSeriesDecompositionDetail {
         //! Set whether or not we're testing for a change.
         void testingForChange(bool value);
 
-        //! Start observing for new components.
-        void observeComponentsAdded();
+        //! Start observing for modifications to the components.
+        void observeComponentsModified();
 
-        //! Check if any components were added since observeComponentsAdded
+        //! Check if any components were added or removed since observeComponentsModified
         //! was last called.
-        bool componentsAdded();
+        bool componentsModified();
 
         //! Apply \p shift to the level at \p time and \p value.
         void shiftLevel(core_t::TTime time, double value, double shift);
@@ -653,12 +655,18 @@ class MATHS_EXPORT CTimeSeriesDecompositionDetail {
             //! Get the memory used by this object.
             std::size_t memoryUsage() const;
 
+            //! Remove any components with invalid values
+            bool removeComponentsWithBadValues(core_t::TTime);
+
         private:
             //! The components.
             maths_t::TSeasonalComponentVec m_Components;
 
             //! The components' prediction errors.
             TComponentErrorsVec m_PredictionErrors;
+
+            //! Befriend a helper class used by the unit tests
+            friend class ::CNanInjector;
         };
 
         using TSeasonalPtr = std::unique_ptr<CSeasonal>;
@@ -728,6 +736,9 @@ class MATHS_EXPORT CTimeSeriesDecompositionDetail {
             //! Get the memory used by this object.
             std::size_t memoryUsage() const;
 
+            //! Remove any components with invalid values
+            bool removeComponentsWithBadValues(core_t::TTime time);
+
         private:
             //! The calendar components.
             maths_t::TCalendarComponentVec m_Components;
@@ -846,8 +857,11 @@ class MATHS_EXPORT CTimeSeriesDecompositionDetail {
         //! Set to true if the trend model should be used for prediction.
         bool m_UsingTrendForPrediction = false;
 
-        //! Set to true when new components are added.
-        bool m_ComponentsAdded = false;
+        //! Set to true when new components are added or removed
+        bool m_ComponentsModified = false;
+
+        //! Befriend a helper class used by the unit tests
+        friend class ::CNanInjector;
     };
 };
 

lib/maths/CCalendarComponent.cc

@@ -181,5 +181,9 @@ std::size_t CCalendarComponent::memoryUsage() const {
     return core::CMemory::dynamicSize(m_Bucketing) +
            core::CMemory::dynamicSize(this->splines());
 }
+
+bool CCalendarComponent::isBad() const {
+    return m_Bucketing.isBad();
+}
 }
 }

lib/maths/CCalendarComponentAdaptiveBucketing.cc

@@ -15,6 +15,7 @@
 
 #include <maths/CBasicStatisticsPersist.h>
 #include <maths/CChecksum.h>
+#include <maths/CMathsFuncs.h>
 #include <maths/CTools.h>
 
 #include <boost/bind.hpp>
@@ -356,5 +357,13 @@ std::string CCalendarComponentAdaptiveBucketing::name() const {
     return "Calendar[" + std::to_string(this->decayRate()) + "," +
            std::to_string(this->minimumBucketLength()) + "]";
 }
+
+bool CCalendarComponentAdaptiveBucketing::isBad() const {
+    // check for bad values in both the means and the variances
+    return std::any_of(m_Values.begin(), m_Values.end(), [](const auto& value) {
+        return ((CMathsFuncs::isFinite(CBasicStatistics::mean(value)) == false) ||
+                (CMathsFuncs::isFinite(CBasicStatistics::variance(value))) == false);
+    });
+}
 }
 }

lib/maths/CDecompositionComponent.cc

@@ -107,6 +107,7 @@ void CDecompositionComponent::interpolate(const TDoubleVec& knots,
                                           const TDoubleVec& values,
                                           const TDoubleVec& variances) {
     m_Splines.interpolate(knots, values, variances, m_BoundaryCondition);
+
     m_MeanValue = this->valueSpline().mean();
     m_MeanVariance = this->varianceSpline().mean();
 }
@@ -123,9 +124,9 @@ TDoubleDoublePr CDecompositionComponent::value(double offset, double n, double c
     // asymptotically normal with mean equal to the sample mean
     // and variance equal to the sample variance divided by root
     // of the number of samples.
-
     if (this->initialized()) {
         double m{this->valueSpline().value(offset)};
+
         if (confidence == 0.0) {
             return {m, m};
         }

lib/maths/CPeriodicityHypothesisTests.cc

@@ -1981,6 +1981,13 @@ bool CPeriodicityHypothesisTests::testPartition(const TTimeTimePr2Vec& partition
         return false;
     }
 
+    // It's possible that none of the candidates are <= 1.05 times the minimum,
+    // this would be the case if a NaN were present in the values say.
+    // NaNs are detected and purged elsewhere so we simply return false here.
+    if (best.count() == 0) {
+        return false;
+    }
+
     startOfPartition = (m_StartTime + best[0].second) % repeat;
     double v1{varianceAtPercentile(correction * minimum[0].first, df1,
                                    50.0 + CONFIDENCE_INTERVAL / 2.0)};

lib/maths/CSeasonalComponentAdaptiveBucketing.cc

@@ -19,6 +19,8 @@
 #include <maths/CLinearAlgebra.h>
 #include <maths/CLinearAlgebraPersist.h>
 #include <maths/CLinearAlgebraTools.h>
+#include <maths/CMathsFuncs.h>
+#include <maths/CMathsFuncsDetail.h>
 #include <maths/CRegression.h>
 #include <maths/CRegressionDetail.h>
 #include <maths/CSeasonalTime.h>
@@ -54,6 +56,13 @@ double gradient(const TRegression& r) {
     return params[1];
 }
 
+//! Safely promote the type of the vector statistics value returned by
+//! CBasicStatistics::mean to a type accepted by CMathsFuncs::isFinite
+template<typename T>
+typename SPromoted<T>::Type promote(const T& vectorStat) {
+    return (vectorStat);
+}
+
 // Version 6.3
 const std::string VERSION_6_3_TAG("6.3");
 const std::string ADAPTIVE_BUCKETING_6_3_TAG{"a"};
@@ -578,6 +587,11 @@ std::string CSeasonalComponentAdaptiveBucketing::name() const {
            std::to_string(this->minimumBucketLength()) + "]";
 }
 
+bool CSeasonalComponentAdaptiveBucketing::isBad() const {
+    return std::any_of(m_Buckets.begin(), m_Buckets.end(),
+                       [](const auto& bucket) { return bucket.isBad(); });
+}
+
 double CSeasonalComponentAdaptiveBucketing::observedInterval(core_t::TTime time) const {
     return m_Time->regressionInterval(
         std::min_element(m_Buckets.begin(), m_Buckets.end(),
@@ -627,5 +641,11 @@ uint64_t CSeasonalComponentAdaptiveBucketing::SBucket::checksum(uint64_t seed) c
     seed = CChecksum::calculate(seed, s_FirstUpdate);
     return CChecksum::calculate(seed, s_LastUpdate);
 }
+
+bool CSeasonalComponentAdaptiveBucketing::SBucket::isBad() const {
+    return ((CMathsFuncs::isFinite(
+                 promote(CBasicStatistics::mean(s_Regression.statistic()))) == false) ||
+            (CMathsFuncs::isFinite(s_Variance) == false));
+}
 }
 }