[ML] Linear scaling change detection (#25)

This implements detection of linear scaling events. It also finishes up the unit testing of change 
detection and fixes some issues these turned up: specifically, 1) the behaviour when a change is 
detected but the trend model has no components, 2) the handling of time shifts in the trend 
model and 3) the handling of data types in the trend component change model. Finally, we are 
now more careful with the weights we apply to samples added to both the standard and change 
models. This has meant I've been able to revert scaling the changes, since the trend is less 
influenced by values during the change detection period if we're likely to detect a change.

Author: tveasey

include/maths/CCalendarComponent.h

@@ -86,6 +86,9 @@ class MATHS_EXPORT CCalendarComponent : private CDecompositionComponent
         //! Clear all data.
         void clear(void);
 
+        //! Linearly scale the component's by \p scale.
+        void linearScale(core_t::TTime time, double scale);
+
         //! Adds a value \f$(t, f(t))\f$ to this component.
         //!
         //! \param[in] time The time of the point.

include/maths/CCalendarComponentAdaptiveBucketing.h

@@ -71,6 +71,9 @@ class MATHS_EXPORT CCalendarComponentAdaptiveBucketing : private CAdaptiveBucket
         //! allocated memory.
         void clear(void);
 
+        //! Linearly scale the bucket values by \p scale.
+        void linearScale(double scale);
+
         //! Add the function value at \p time.
         //!
         //! \param[in] time The time of \p value.

include/maths/CModel.h

@@ -74,7 +74,7 @@ class MATHS_EXPORT CModelParams
         bool testForChange(core_t::TTime changeInterval) const;
 
         //! Get the minimum time to detect a change point in the model.
-        core_t::TTime minimumTimeToDetectChange(core_t::TTime timeSinceLastChangePoint) const;
+        core_t::TTime minimumTimeToDetectChange(void) const;
 
         //! Get the maximum time to test for a change point in the model.
         core_t::TTime maximumTimeToTestForChange(void) const;

include/maths/CRegression.h

@@ -235,6 +235,23 @@ class MATHS_EXPORT CRegression
                     }
                 }
 
+                //! Linearly scale the regression model.
+                //!
+                //! i.e. apply a transform such that each regression parameter maps
+                //! to \p scale times its current value.
+                //!
+                //! \param[in] scale The scale to apply to the regression parameters.
+                void linearScale(double scale)
+                {
+                    if (CBasicStatistics::count(m_S) > 0.0)
+                    {
+                        for (std::size_t i = 0u; i < N; ++i)
+                        {
+                            CBasicStatistics::moment<0>(m_S)(i+2*N-1) *= scale;
+                        }
+                    }
+                }
+
                 //! Multiply the statistics' count by \p scale.
                 CLeastSquaresOnline scaled(double scale) const
                 {
@@ -272,14 +289,13 @@ class MATHS_EXPORT CRegression
                     if (this->parameters(params, maxCondition))
                     {
                         std::ptrdiff_t n = static_cast<std::ptrdiff_t>(params.size());
-                        double xi = x;
                         for (std::ptrdiff_t i = n - 1; i >= 0; --i)
                         {
                             result[i] = params[i];
                             for (std::ptrdiff_t j = i + 1; j < n; ++j)
                             {
                                 params[j] *=  static_cast<double>(i + 1)
-                                            / static_cast<double>(j - i) * xi;
+                                            / static_cast<double>(j - i) * x;
                                 result[i] += params[j];
                             }
                         }

include/maths/CSeasonalComponent.h

@@ -104,6 +104,9 @@ class MATHS_EXPORT CSeasonalComponent : private CDecompositionComponent
         //! Shift the component's slope by \p shift.
         void shiftSlope(double shift);
 
+        //! Linearly scale the component's by \p scale.
+        void linearScale(core_t::TTime time, double scale);
+
         //! Adds a value \f$(t, f(t))\f$ to this component.
         //!
         //! \param[in] time The time of the point.

include/maths/CSeasonalComponentAdaptiveBucketing.h

@@ -97,6 +97,9 @@ class MATHS_EXPORT CSeasonalComponentAdaptiveBucketing : private CAdaptiveBucket
         //! Shift the regressions' gradients by \p shift.
         void shiftSlope(double shift);
 
+        //! Linearly scale the regressions by \p scale.
+        void linearScale(double scale);
+
         //! Add the function value at \p time.
         //!
         //! \param[in] time The time of \p value.

include/maths/CTimeSeriesChangeDetector.h

@@ -50,6 +50,7 @@ struct MATHS_EXPORT SChangeDescription
     enum EDescription
     {
         E_LevelShift,
+        E_LinearScale,
         E_TimeShift
     };
 
@@ -88,7 +89,7 @@ class MATHS_EXPORT CUnivariateTimeSeriesChangeDetector
     public:
         CUnivariateTimeSeriesChangeDetector(const TDecompositionPtr &trendModel,
                                             const TPriorPtr &residualModel,
-                                            core_t::TTime minimumTimeToDetect = 6 * core::constants::HOUR,
+                                            core_t::TTime minimumTimeToDetect = 12 * core::constants::HOUR,
                                             core_t::TTime maximumTimeToDetect = core::constants::DAY,
                                             double minimumDeltaBicToDetect = 14.0);
 
@@ -103,6 +104,13 @@ class MATHS_EXPORT CUnivariateTimeSeriesChangeDetector
         //! if there has been.
         TOptionalChangeDescription change();
 
+        //! The function used to decide whether to accept a change.
+        //! A change is accepted at a value of 1.0 for this function.
+        //!
+        //! \param[out] change Filled in with the index of the change
+        //! the most likely change.
+        double decisionFunction(std::size_t &change) const;
+
         //! Add \p samples to the change detector.
         void addSamples(const TWeightStyleVec &weightStyles,
                         const TTimeDoublePr1Vec &samples,
@@ -123,7 +131,7 @@ class MATHS_EXPORT CUnivariateTimeSeriesChangeDetector
     private:
         using TChangeModel = time_series_change_detector_detail::CUnivariateChangeModel;
         using TChangeModelPtr = boost::shared_ptr<TChangeModel>;
-        using TChangeModelPtr4Vec = core::CSmallVector<TChangeModelPtr, 4>;
+        using TChangeModelPtr5Vec = core::CSmallVector<TChangeModelPtr, 5>;
         using TMinMaxAccumulator = CBasicStatistics::CMinMax<core_t::TTime>;
 
     private:
@@ -147,7 +155,7 @@ class MATHS_EXPORT CUnivariateTimeSeriesChangeDetector
         double m_CurrentEvidenceOfChange;
 
         //! The change models.
-        TChangeModelPtr4Vec m_ChangeModels;
+        TChangeModelPtr5Vec m_ChangeModels;
 };
 
 namespace time_series_change_detector_detail
@@ -158,6 +166,7 @@ namespace time_series_change_detector_detail
 class MATHS_EXPORT CUnivariateChangeModel : private core::CNonCopyable
 {
     public:
+        using TDouble1Vec = core::CSmallVector<double, 1>;
         using TDouble4Vec = core::CSmallVector<double, 4>;
         using TDouble4Vec1Vec = core::CSmallVector<TDouble4Vec, 1>;
         using TTimeDoublePr = std::pair<core_t::TTime, double>;
@@ -189,10 +198,10 @@ class MATHS_EXPORT CUnivariateChangeModel : private core::CNonCopyable
         virtual TOptionalChangeDescription change() const = 0;
 
         //! Update the change model with \p samples.
-        virtual void addSamples(std::size_t count,
-                                const TWeightStyleVec &weightStyles,
+        virtual void addSamples(const std::size_t count,
+                                TWeightStyleVec weightStyles,
                                 const TTimeDoublePr1Vec &samples,
-                                const TDouble4Vec1Vec &weights) = 0;
+                                TDouble4Vec1Vec weights) = 0;
 
         //! Debug the memory used by this object.
         void debugMemoryUsage(core::CMemoryUsage::TMemoryUsagePtr mem) const;
@@ -206,29 +215,28 @@ class MATHS_EXPORT CUnivariateChangeModel : private core::CNonCopyable
         //! Get a checksum for this object.
         virtual uint64_t checksum(uint64_t seed) const = 0;
 
-    protected:
-        //! The sample count to initialize a change model.
-        static const std::size_t COUNT_TO_INITIALIZE{5u};
-
     protected:
         //! Restore the residual model reading state from \p traverser.
         bool restoreResidualModel(const SDistributionRestoreParams &params,
                                   core::CStateRestoreTraverser &traverser);
 
         //! Get the log-likelihood.
         double logLikelihood() const;
-        //! Update the data log-likelihood with \p logLikelihood.
-        void addLogLikelihood(double logLikelihood);
 
         //! Get the expected log-likelihood.
         double expectedLogLikelihood() const;
-        //! Update the expected data log-likelihood with \p logLikelihood.
-        void addExpectedLogLikelihood(double logLikelihood);
+
+        //! Update the log-likelihood with \p samples.
+        void updateLogLikelihood(const TWeightStyleVec &weightStyles,
+                                 const TDouble1Vec &samples,
+                                 const TDouble4Vec1Vec &weights);
+
+        //! Update the expected log-likelihoods.
+        void updateExpectedLogLikelihood(const TWeightStyleVec &weightStyles,
+                                         const TDouble4Vec1Vec &weights);
 
         //! Get the time series trend model.
         const CTimeSeriesDecompositionInterface &trendModel() const;
-        //! Get the time series trend model.
-        CTimeSeriesDecompositionInterface &trendModel();
 
         //! Get the time series residual model.
         const CPrior &residualModel() const;
@@ -275,10 +283,10 @@ class MATHS_EXPORT CUnivariateNoChangeModel final : public CUnivariateChangeMode
         virtual TOptionalChangeDescription change() const;
 
         //! Get the log likelihood of \p samples.
-        virtual void addSamples(std::size_t count,
-                                const TWeightStyleVec &weightStyles,
+        virtual void addSamples(const std::size_t count,
+                                TWeightStyleVec weightStyles,
                                 const TTimeDoublePr1Vec &samples,
-                                const TDouble4Vec1Vec &weights);
+                                TDouble4Vec1Vec weights);
 
         //! Get the static size of this object.
         virtual std::size_t staticSize() const;
@@ -312,10 +320,10 @@ class MATHS_EXPORT CUnivariateLevelShiftModel final : public CUnivariateChangeMo
         virtual TOptionalChangeDescription change() const;
 
         //! Update with \p samples.
-        virtual void addSamples(std::size_t count,
-                                const TWeightStyleVec &weightStyles,
+        virtual void addSamples(const std::size_t count,
+                                TWeightStyleVec weightStyles,
                                 const TTimeDoublePr1Vec &samples,
-                                const TDouble4Vec1Vec &weights);
+                                TDouble4Vec1Vec weights);
 
         //! Get the static size of this object.
         virtual std::size_t staticSize() const;
@@ -324,7 +332,6 @@ class MATHS_EXPORT CUnivariateLevelShiftModel final : public CUnivariateChangeMo
         virtual uint64_t checksum(uint64_t seed) const;
 
     private:
-        using TDoubleVec = std::vector<double>;
         using TMeanAccumulator = CBasicStatistics::SSampleMean<double>::TAccumulator;
 
     private:
@@ -338,6 +345,56 @@ class MATHS_EXPORT CUnivariateLevelShiftModel final : public CUnivariateChangeMo
         double m_SampleCount;
 };
 
+//! \brief Captures the likelihood of the data given an arbitrary
+//! linear scaling.
+class MATHS_EXPORT CUnivariateLinearScaleModel final : public CUnivariateChangeModel
+{
+    public:
+        CUnivariateLinearScaleModel(const TDecompositionPtr &trendModel,
+                                    const TPriorPtr &residualModel);
+
+        //! Initialize by reading state from \p traverser.
+        virtual bool acceptRestoreTraverser(const SModelRestoreParams &params,
+                                            core::CStateRestoreTraverser &traverser);
+
+        //! Persist state by passing information to \p inserter.
+        virtual void acceptPersistInserter(core::CStatePersistInserter &inserter) const;
+
+        //! The BIC of applying the level shift.
+        virtual double bic() const;
+
+        //! The expected BIC of applying the change.
+        virtual double expectedBic() const;
+
+        //! Get a description of the level shift.
+        virtual TOptionalChangeDescription change() const;
+
+        //! Update with \p samples.
+        virtual void addSamples(const std::size_t count,
+                                TWeightStyleVec weightStyles,
+                                const TTimeDoublePr1Vec &samples,
+                                TDouble4Vec1Vec weights);
+
+        //! Get the static size of this object.
+        virtual std::size_t staticSize() const;
+
+        //! Get a checksum for this object.
+        virtual uint64_t checksum(uint64_t seed) const;
+
+    private:
+        using TMeanAccumulator = CBasicStatistics::SSampleMean<double>::TAccumulator;
+
+    private:
+        //! The optimal shift.
+        TMeanAccumulator m_Scale;
+
+        //! The mode of the initial residual distribution model.
+        double m_ResidualModelMode;
+
+        //! The number of samples added so far.
+        double m_SampleCount;
+};
+
 //! \brief Captures the likelihood of the data given a specified
 //! time shift.
 class MATHS_EXPORT CUnivariateTimeShiftModel final : public CUnivariateChangeModel
@@ -364,10 +421,10 @@ class MATHS_EXPORT CUnivariateTimeShiftModel final : public CUnivariateChangeMod
         virtual TOptionalChangeDescription change() const;
 
         //! Update with \p samples.
-        virtual void addSamples(std::size_t count,
-                                const TWeightStyleVec &weightStyles,
+        virtual void addSamples(const std::size_t count,
+                                TWeightStyleVec weightStyles,
                                 const TTimeDoublePr1Vec &samples,
-                                const TDouble4Vec1Vec &weights);
+                                TDouble4Vec1Vec weights);
 
         //! Get the static size of this object.
         virtual std::size_t staticSize() const;

include/maths/CTimeSeriesDecomposition.h

@@ -120,11 +120,12 @@ class MATHS_EXPORT CTimeSeriesDecomposition : public CTimeSeriesDecompositionInt
         //! \param[in] value The value immediately before the change
         //! point.
         //! \param[in] change A description of the change to apply.
-        virtual void applyChange(core_t::TTime time, double value,
+        //! \return True if a new component was detected.
+        virtual bool applyChange(core_t::TTime time, double value,
                                  const SChangeDescription &change);
 
         //! Propagate the decomposition forwards to \p time.
-        void propagateForwardsTo(core_t::TTime time);
+        virtual void propagateForwardsTo(core_t::TTime time);
 
         //! Get the mean value of the time series in the vicinity of \p time.
         virtual double meanValue(core_t::TTime time) const;
@@ -192,10 +193,14 @@ class MATHS_EXPORT CTimeSeriesDecomposition : public CTimeSeriesDecompositionInt
         //! Get the static size of this object.
         virtual std::size_t staticSize(void) const;
 
+        //! Get the time shift which is being applied.
+        virtual core_t::TTime timeShift(void) const;
+
         //! Get the seasonal components.
         virtual const maths_t::TSeasonalComponentVec &seasonalComponents(void) const;
 
-        //! This is the latest time of any point added to this object or the time skipped to.
+        //! This is the latest time of any point added to this object or
+        //! the time skipped to.
         virtual core_t::TTime lastValueTime(void) const;
 
     private:

include/maths/CTimeSeriesDecompositionDetail.h

@@ -389,9 +389,15 @@ class MATHS_EXPORT CTimeSeriesDecompositionDetail
                 //! Create a new calendar component.
                 virtual void handle(const SDetectedCalendar &message);
 
-                //! Apply \p change at \p time.
+                //! Start using the trend for prediction.
+                void useTrendForPrediction(void);
+
+                //! Apply \p shift to the level at \p time and \p value.
                 void shiftLevel(core_t::TTime time, double value, double shift);
 
+                //! Apply a linear scale of \p scale.
+                void linearScale(core_t::TTime time, double scale);
+
                 //! Maybe re-interpolate the components.
                 void interpolate(const SMessage &message);
 
@@ -549,6 +555,9 @@ class MATHS_EXPORT CTimeSeriesDecompositionDetail
                     //! Shift the components' time origin to \p time.
                     void shiftOrigin(core_t::TTime time);
 
+                    //! Linearly scale the components' by \p scale.
+                    void linearScale(core_t::TTime time, double scale);
+
                     //! Get a checksum for this object.
                     uint64_t checksum(uint64_t seed = 0) const;
 
@@ -608,6 +617,9 @@ class MATHS_EXPORT CTimeSeriesDecompositionDetail
                     //! Remove low value components.
                     bool prune(core_t::TTime time, core_t::TTime bucketLength);
 
+                    //! Linearly scale the components' by \p scale.
+                    void linearScale(core_t::TTime time, double scale);
+
                     //! Get a checksum for this object.
                     uint64_t checksum(uint64_t seed = 0) const;
 

include/maths/CTimeSeriesDecompositionInterface.h

@@ -98,7 +98,8 @@ class MATHS_EXPORT CTimeSeriesDecompositionInterface
         //! \param[in] value The value immediately before the change
         //! point.
         //! \param[in] change A description of the change to apply.
-        virtual void applyChange(core_t::TTime time, double value,
+        //! \return True if a new component was detected.
+        virtual bool applyChange(core_t::TTime time, double value,
                                  const SChangeDescription &change) = 0;
 
         //! Propagate the decomposition forwards to \p time.
@@ -171,10 +172,14 @@ class MATHS_EXPORT CTimeSeriesDecompositionInterface
         //! Get the static size of this object.
         virtual std::size_t staticSize(void) const = 0;
 
+        //! Get the time shift which is being applied.
+        virtual core_t::TTime timeShift(void) const = 0;
+
         //! Get the seasonal components.
         virtual const maths_t::TSeasonalComponentVec &seasonalComponents(void) const = 0;
 
-        //! This is the latest time of any point added to this object or the time skipped to.
+        //! This is the latest time of any point added to this object or
+        //! the time skipped to.
         virtual core_t::TTime lastValueTime(void) const = 0;
 };