|
// Licensed to the .NET Foundation under one or more agreements.
// The .NET Foundation licenses this file to you under the MIT license.
// See the LICENSE file in the project root for more information.
using System;
using System.Collections.Generic;
using System.Numerics;
using Microsoft.ML;
using Microsoft.ML.CommandLine;
using Microsoft.ML.Data;
using Microsoft.ML.Internal.CpuMath;
using Microsoft.ML.Internal.Utilities;
using Microsoft.ML.Runtime;
using Microsoft.ML.Transforms.TimeSeries;
[assembly: LoadableClass(typeof(AdaptiveSingularSpectrumSequenceModelerInternal),
typeof(AdaptiveSingularSpectrumSequenceModelerInternal), null, typeof(SignatureLoadModel),
"SSA Sequence Modeler",
AdaptiveSingularSpectrumSequenceModelerInternal.LoaderSignature)]
namespace Microsoft.ML.Transforms.TimeSeries
{
/// <summary>
/// Ranking selection method for the signal.
/// </summary>
public enum RankSelectionMethod
{
Fixed,
Exact,
Fast
}
/// <summary>
/// Growth ratio. Defined as Growth^(1/TimeSpan).
/// </summary>
public struct GrowthRatio
{
[Argument(ArgumentType.AtMostOnce, HelpText = "Time span of growth ratio. Must be strictly positive.", SortOrder = 1)]
public int TimeSpan;
[Argument(ArgumentType.AtMostOnce, HelpText = "Growth. Must be non-negative.", SortOrder = 2)]
public Double Growth;
public Double Ratio { get { return Math.Pow(Growth, 1d / TimeSpan); } }
}
/// <summary>
/// This class implements basic Singular Spectrum Analysis (SSA) model for modeling univariate time-series.
/// For the details of the model, refer to http://arxiv.org/pdf/1206.6910.pdf.
/// </summary>
internal sealed class AdaptiveSingularSpectrumSequenceModelerInternal : SequenceModelerBase<Single, Single>
{
internal const string LoaderSignature = "SSAModel";
internal sealed class SsaForecastResult : ForecastResultBase<Single>
{
public VBuffer<Single> ForecastStandardDeviation;
public VBuffer<Single> UpperBound;
public VBuffer<Single> LowerBound;
public Single ConfidenceLevel;
internal bool CanComputeForecastIntervals;
internal Single BoundOffset;
public bool IsVarianceValid { get { return CanComputeForecastIntervals; } }
}
public sealed class ModelInfo
{
/// <summary>
/// The singular values of the SSA of the input time-series
/// </summary>
public Single[] Spectrum;
/// <summary>
/// The roots of the characteristic polynomial after stabilization (meaningful only if the model is stabilized.)
/// </summary>
public Complex[] RootsAfterStabilization;
/// <summary>
/// The roots of the characteristic polynomial before stabilization (meaningful only if the model is stabilized.)
/// </summary>
public Complex[] RootsBeforeStabilization;
/// <summary>
/// The rank of the model
/// </summary>
public int Rank;
/// <summary>
/// The window size used to compute the SSA model
/// </summary>
public int WindowSize;
/// <summary>
/// The auto-regressive coefficients learned by the model
/// </summary>
public Single[] AutoRegressiveCoefficients;
/// <summary>
/// The flag indicating whether the model has been trained
/// </summary>
public bool IsTrained;
/// <summary>
/// The flag indicating a naive model is trained instead of SSA
/// </summary>
public bool IsNaiveModelTrained;
/// <summary>
/// The flag indicating whether the learned model has an exponential trend (meaningful only if the model is stabilized.)
/// </summary>
public bool IsExponentialTrendPresent;
/// <summary>
/// The flag indicating whether the learned model has a polynomial trend (meaningful only if the model is stabilized.)
/// </summary>
public bool IsPolynomialTrendPresent;
/// <summary>
/// The flag indicating whether the learned model has been stabilized
/// </summary>
public bool IsStabilized;
/// <summary>
/// The flag indicating whether any artificial seasonality (a seasonality with period greater than the window size) is removed
/// (meaningful only if the model is stabilized.)
/// </summary>
public bool IsArtificialSeasonalityRemoved;
/// <summary>
/// The exponential trend magnitude (meaningful only if the model is stabilized.)
/// </summary>
public Double ExponentialTrendFactor;
}
private ModelInfo _info;
/// <summary>
/// Returns the meta information about the learned model.
/// </summary>
public ModelInfo Info
{
get { return _info; }
}
private Single[] _alpha;
private Single[] _state;
private readonly FixedSizeQueue<Single> _buffer;
private CpuAlignedVector _x;
private CpuAlignedVector _xSmooth;
private int _windowSize;
private readonly int _seriesLength;
private readonly RankSelectionMethod _rankSelectionMethod;
private readonly Single _discountFactor;
private readonly int _trainSize;
private int _maxRank;
private readonly Double _maxTrendRatio;
private readonly bool _shouldStablize;
private readonly bool _shouldMaintainInfo;
private readonly IHost _host;
private CpuAlignedMatrixRow _wTrans;
private Single _observationNoiseVariance;
private Single _observationNoiseMean;
private Single _autoregressionNoiseVariance;
private Single _autoregressionNoiseMean;
private int _rank;
private CpuAlignedVector _y;
private Single _nextPrediction;
/// <summary>
/// Determines whether the confidence interval required for forecasting.
/// </summary>
public bool ShouldComputeForecastIntervals { get; set; }
/// <summary>
/// Returns the rank of the subspace used for SSA projection (parameter r).
/// </summary>
public int Rank { get { return _rank; } }
/// <summary>
/// Returns the smoothed (via SSA projection) version of the last series observation fed to the model.
/// </summary>
public Single LastSmoothedValue { get { return _state[_windowSize - 2]; } }
/// <summary>
/// Returns the last series observation fed to the model.
/// </summary>
public Single LastValue { get { return _buffer.Count > 0 ? _buffer[_buffer.Count - 1] : Single.NaN; } }
private static VersionInfo GetVersionInfo()
{
return new VersionInfo(
modelSignature: "SSAMODLR",
//verWrittenCur: 0x00010001, // Initial
verWrittenCur: 0x00010002, // Added saving _state and _nextPrediction
verReadableCur: 0x00010002,
verWeCanReadBack: 0x00010001,
loaderSignature: LoaderSignature,
loaderAssemblyName: typeof(AdaptiveSingularSpectrumSequenceModelerInternal).Assembly.FullName);
}
private const int VersionSavingStateAndPrediction = 0x00010002;
/// <summary>
/// The constructor for Adaptive SSA model.
/// </summary>
/// <param name="env">The exception context.</param>
/// <param name="trainSize">The length of series from the beginning used for training (parameter N).
/// Must be at least twice the windowSize.</param>
/// <param name="seriesLength">This parameter must be greater than windowSize.</param>
/// <param name="windowSize">The length of the window on the series for building the trajectory matrix (parameter L).
/// We recommend you set this to be more than twice the maximum seasonality in the data. For example,
/// if your data can exhibit both monthly and yearly seasonality, and you have data points from each day,
/// set this to be twice the number of days in a year.</param>
/// <param name="discountFactor">The discount factor in [0,1] used for online updates (default = 1).</param>
/// <param name="rankSelectionMethod">The rank selection method (default = Exact).</param>
/// <param name="rank">The desired rank of the subspace used for SSA projection (parameter r). This parameter should be in the range in [1, windowSize].
/// If set to null, the rank is automatically determined based on prediction error minimization. (default = null)</param>
/// <param name="maxRank">The maximum rank considered during the rank selection process. If not provided (i.e. set to null), it is set to windowSize - 1.</param>
/// <param name="shouldComputeForecastIntervals">The flag determining whether the confidence bounds for the point forecasts should be computed. (default = true)</param>
/// <param name="shouldstablize">The flag determining whether the model should be stabilized.</param>
/// <param name="shouldMaintainInfo">The flag determining whether the meta information for the model needs to be maintained.</param>
/// <param name="maxGrowth">The maximum growth on the exponential trend.</param>
public AdaptiveSingularSpectrumSequenceModelerInternal(IHostEnvironment env, int trainSize, int seriesLength, int windowSize, Single discountFactor = 1,
RankSelectionMethod rankSelectionMethod = RankSelectionMethod.Exact, int? rank = null, int? maxRank = null,
bool shouldComputeForecastIntervals = true, bool shouldstablize = true, bool shouldMaintainInfo = false, GrowthRatio? maxGrowth = null)
: base()
{
Contracts.CheckValue(env, nameof(env));
_host = env.Register(LoaderSignature);
_host.CheckParam(windowSize >= 2, nameof(windowSize), "The window size should be at least 2."); // ...because otherwise we have nothing to autoregress on
_host.CheckParam(seriesLength > windowSize, nameof(seriesLength), "The series length should be greater than the window size.");
_host.Check(trainSize > 2 * windowSize, "The input size for training should be greater than twice the window size.");
_host.CheckParam(0 <= discountFactor && discountFactor <= 1, nameof(discountFactor), "The discount factor should be in [0,1].");
_host.CheckParam(maxGrowth == null || maxGrowth.Value.TimeSpan > 0, nameof(GrowthRatio.TimeSpan), "The time span must be strictly positive.");
_host.CheckParam(maxGrowth == null || maxGrowth.Value.Growth >= 0, nameof(GrowthRatio.Growth), "The growth must be non-negative.");
if (maxRank != null)
{
_maxRank = (int)maxRank;
_host.CheckParam(1 <= _maxRank && _maxRank < windowSize, nameof(maxRank),
"The max rank should be in [1, windowSize).");
}
else
_maxRank = windowSize - 1;
_rankSelectionMethod = rankSelectionMethod;
if (_rankSelectionMethod == RankSelectionMethod.Fixed)
{
if (rank != null)
{
_rank = (int)rank;
_host.CheckParam(1 <= _rank && _rank < windowSize, nameof(rank), "The rank should be in [1, windowSize).");
}
else
_rank = _maxRank;
}
_seriesLength = seriesLength;
_windowSize = windowSize;
_trainSize = trainSize;
_discountFactor = discountFactor;
_buffer = new FixedSizeQueue<Single>(seriesLength);
_alpha = new Single[windowSize - 1];
_state = new Single[windowSize - 1];
_x = new CpuAlignedVector(windowSize, CpuMathUtils.GetVectorAlignment());
_xSmooth = new CpuAlignedVector(windowSize, CpuMathUtils.GetVectorAlignment());
ShouldComputeForecastIntervals = shouldComputeForecastIntervals;
_observationNoiseVariance = 0;
_autoregressionNoiseVariance = 0;
_observationNoiseMean = 0;
_autoregressionNoiseMean = 0;
_shouldStablize = shouldstablize;
_maxTrendRatio = maxGrowth == null ? Double.PositiveInfinity : ((GrowthRatio)maxGrowth).Ratio;
_shouldMaintainInfo = shouldMaintainInfo;
if (_shouldMaintainInfo)
{
_info = new ModelInfo();
_info.WindowSize = _windowSize;
}
}
/// <summary>
/// The copy constructor.
/// </summary>
/// <param name="model">An object whose contents are copied to the current object.</param>
private AdaptiveSingularSpectrumSequenceModelerInternal(AdaptiveSingularSpectrumSequenceModelerInternal model)
{
_host = model._host.Register(LoaderSignature);
_host.Assert(model._windowSize >= 2);
_host.Assert(model._seriesLength > model._windowSize);
_host.Assert(model._trainSize > 2 * model._windowSize);
_host.Assert(0 <= model._discountFactor && model._discountFactor <= 1);
_host.Assert(1 <= model._rank && model._rank < model._windowSize);
_rank = model._rank;
_maxRank = model._maxRank;
_rankSelectionMethod = model._rankSelectionMethod;
_seriesLength = model._seriesLength;
_windowSize = model._windowSize;
_trainSize = model._trainSize;
_discountFactor = model._discountFactor;
ShouldComputeForecastIntervals = model.ShouldComputeForecastIntervals;
_observationNoiseVariance = model._observationNoiseVariance;
_autoregressionNoiseVariance = model._autoregressionNoiseVariance;
_observationNoiseMean = model._observationNoiseMean;
_autoregressionNoiseMean = model._autoregressionNoiseMean;
_nextPrediction = model._nextPrediction;
_maxTrendRatio = model._maxTrendRatio;
_shouldStablize = model._shouldStablize;
_shouldMaintainInfo = model._shouldMaintainInfo;
_info = model._info;
_buffer = model._buffer.Clone();
_alpha = new Single[_windowSize - 1];
Array.Copy(model._alpha, _alpha, _windowSize - 1);
_state = new Single[_windowSize - 1];
Array.Copy(model._state, _state, _windowSize - 1);
_x = new CpuAlignedVector(_windowSize, CpuMathUtils.GetVectorAlignment());
_xSmooth = new CpuAlignedVector(_windowSize, CpuMathUtils.GetVectorAlignment());
if (model._wTrans != null)
{
_y = new CpuAlignedVector(_rank, CpuMathUtils.GetVectorAlignment());
_wTrans = new CpuAlignedMatrixRow(_rank, _windowSize, CpuMathUtils.GetVectorAlignment());
_wTrans.CopyFrom(model._wTrans);
}
}
public AdaptiveSingularSpectrumSequenceModelerInternal(IHostEnvironment env, ModelLoadContext ctx)
{
Contracts.CheckValue(env, nameof(env));
_host = env.Register(LoaderSignature);
// *** Binary format ***
// int: _seriesLength
// int: _windowSize
// int: _trainSize
// int: _rank
// float: _discountFactor
// RankSelectionMethod: _rankSelectionMethod
// bool: isWeightSet
// float[]: _alpha
// float[]: _state
// bool: ShouldComputeForecastIntervals
// float: _observationNoiseVariance
// float: _autoregressionNoiseVariance
// float: _observationNoiseMean
// float: _autoregressionNoiseMean
// float: _nextPrediction
// int: _maxRank
// bool: _shouldStablize
// bool: _shouldMaintainInfo
// double: _maxTrendRatio
// float[]: _wTrans (only if _isWeightSet == true)
_seriesLength = ctx.Reader.ReadInt32();
// Do an early check. We'll have the stricter check later.
_host.CheckDecode(_seriesLength > 2);
_windowSize = ctx.Reader.ReadInt32();
_host.CheckDecode(_windowSize >= 2);
_host.CheckDecode(_seriesLength > _windowSize);
_trainSize = ctx.Reader.ReadInt32();
_host.CheckDecode(_trainSize > 2 * _windowSize);
_rank = ctx.Reader.ReadInt32();
_host.CheckDecode(1 <= _rank && _rank < _windowSize);
_discountFactor = ctx.Reader.ReadSingle();
_host.CheckDecode(0 <= _discountFactor && _discountFactor <= 1);
byte temp;
temp = ctx.Reader.ReadByte();
_rankSelectionMethod = (RankSelectionMethod)temp;
bool isWeightSet = ctx.Reader.ReadBoolByte();
_alpha = ctx.Reader.ReadFloatArray();
_host.CheckDecode(Utils.Size(_alpha) == _windowSize - 1);
if (ctx.Header.ModelVerReadable >= VersionSavingStateAndPrediction)
{
_state = ctx.Reader.ReadFloatArray();
_host.CheckDecode(Utils.Size(_state) == _windowSize - 1);
}
else
_state = new Single[_windowSize - 1];
ShouldComputeForecastIntervals = ctx.Reader.ReadBoolByte();
_observationNoiseVariance = ctx.Reader.ReadSingle();
_host.CheckDecode(_observationNoiseVariance >= 0);
_autoregressionNoiseVariance = ctx.Reader.ReadSingle();
_host.CheckDecode(_autoregressionNoiseVariance >= 0);
_observationNoiseMean = ctx.Reader.ReadSingle();
_autoregressionNoiseMean = ctx.Reader.ReadSingle();
if (ctx.Header.ModelVerReadable >= VersionSavingStateAndPrediction)
_nextPrediction = ctx.Reader.ReadSingle();
_maxRank = ctx.Reader.ReadInt32();
_host.CheckDecode(1 <= _maxRank && _maxRank <= _windowSize - 1);
_shouldStablize = ctx.Reader.ReadBoolByte();
_shouldMaintainInfo = ctx.Reader.ReadBoolByte();
if (_shouldMaintainInfo)
{
_info = new ModelInfo();
_info.AutoRegressiveCoefficients = new Single[_windowSize - 1];
Array.Copy(_alpha, _info.AutoRegressiveCoefficients, _windowSize - 1);
_info.IsStabilized = _shouldStablize;
_info.Rank = _rank;
_info.WindowSize = _windowSize;
}
_maxTrendRatio = ctx.Reader.ReadDouble();
_host.CheckDecode(_maxTrendRatio >= 0);
if (isWeightSet)
{
var tempArray = ctx.Reader.ReadFloatArray();
_host.CheckDecode(Utils.Size(tempArray) == _rank * _windowSize);
_wTrans = new CpuAlignedMatrixRow(_rank, _windowSize, CpuMathUtils.GetVectorAlignment());
int i = 0;
_wTrans.CopyFrom(tempArray, ref i);
tempArray = ctx.Reader.ReadFloatArray();
i = 0;
_y = new CpuAlignedVector(_rank, CpuMathUtils.GetVectorAlignment());
_y.CopyFrom(tempArray, ref i);
}
_buffer = TimeSeriesUtils.DeserializeFixedSizeQueueSingle(ctx.Reader, _host);
_x = new CpuAlignedVector(_windowSize, CpuMathUtils.GetVectorAlignment());
_xSmooth = new CpuAlignedVector(_windowSize, CpuMathUtils.GetVectorAlignment());
}
private protected override void SaveModel(ModelSaveContext ctx)
{
_host.CheckValue(ctx, nameof(ctx));
ctx.CheckAtModel();
ctx.SetVersionInfo(GetVersionInfo());
_host.Assert(_windowSize >= 2);
_host.Assert(_seriesLength > _windowSize);
_host.Assert(_trainSize > 2 * _windowSize);
_host.Assert(0 <= _discountFactor && _discountFactor <= 1);
_host.Assert(1 <= _rank && _rank < _windowSize);
_host.Assert(Utils.Size(_alpha) == _windowSize - 1);
_host.Assert(_observationNoiseVariance >= 0);
_host.Assert(_autoregressionNoiseVariance >= 0);
_host.Assert(1 <= _maxRank && _maxRank <= _windowSize - 1);
_host.Assert(_maxTrendRatio >= 0);
// *** Binary format ***
// int: _seriesLength
// int: _windowSize
// int: _trainSize
// int: _rank
// float: _discountFactor
// RankSelectionMethod: _rankSelectionMethod
// bool: _isWeightSet
// float[]: _alpha
// float[]: _state
// bool: ShouldComputeForecastIntervals
// float: _observationNoiseVariance
// float: _autoregressionNoiseVariance
// float: _observationNoiseMean
// float: _autoregressionNoiseMean
// float: _nextPrediction
// int: _maxRank
// bool: _shouldStablize
// bool: _shouldMaintainInfo
// double: _maxTrendRatio
// float[]: _wTrans (only if _isWeightSet == true)
ctx.Writer.Write(_seriesLength);
ctx.Writer.Write(_windowSize);
ctx.Writer.Write(_trainSize);
ctx.Writer.Write(_rank);
ctx.Writer.Write(_discountFactor);
ctx.Writer.Write((byte)_rankSelectionMethod);
ctx.Writer.WriteBoolByte(_wTrans != null);
ctx.Writer.WriteSingleArray(_alpha);
ctx.Writer.WriteSingleArray(_state);
ctx.Writer.WriteBoolByte(ShouldComputeForecastIntervals);
ctx.Writer.Write(_observationNoiseVariance);
ctx.Writer.Write(_autoregressionNoiseVariance);
ctx.Writer.Write(_observationNoiseMean);
ctx.Writer.Write(_autoregressionNoiseMean);
ctx.Writer.Write(_nextPrediction);
ctx.Writer.Write(_maxRank);
ctx.Writer.WriteBoolByte(_shouldStablize);
ctx.Writer.WriteBoolByte(_shouldMaintainInfo);
ctx.Writer.Write(_maxTrendRatio);
if (_wTrans != null)
{
// REVIEW: this may not be the most efficient way for serializing an aligned matrix.
var tempArray = new Single[_rank * _windowSize];
int iv = 0;
_wTrans.CopyTo(tempArray, ref iv);
ctx.Writer.WriteSingleArray(tempArray);
tempArray = new float[_rank];
iv = 0;
_y.CopyTo(tempArray, ref iv);
ctx.Writer.WriteSingleArray(tempArray);
}
TimeSeriesUtils.SerializeFixedSizeQueue(_buffer, ctx.Writer);
}
private static void ReconstructSignal(TrajectoryMatrix tMat, Single[] singularVectors, int rank, Single[] output)
{
Contracts.Assert(tMat != null);
Contracts.Assert(1 <= rank && rank <= tMat.WindowSize);
Contracts.Assert(Utils.Size(singularVectors) >= tMat.WindowSize * rank);
Contracts.Assert(Utils.Size(output) >= tMat.SeriesLength);
var k = tMat.SeriesLength - tMat.WindowSize + 1;
Contracts.Assert(k > 0);
var v = new Single[k];
int i;
for (i = 0; i < tMat.SeriesLength; ++i)
output[i] = 0;
for (i = 0; i < rank; ++i)
{
// Reconstructing the i-th eigen triple component and adding it to output
tMat.MultiplyTranspose(singularVectors, v, false, tMat.WindowSize * i, 0);
tMat.RankOneHankelization(singularVectors, v, 1, output, true, tMat.WindowSize * i, 0, 0);
}
}
private static void ReconstructSignalTailFast(Single[] series, TrajectoryMatrix tMat, Single[] singularVectors, int rank, Single[] output)
{
Contracts.Assert(tMat != null);
Contracts.Assert(1 <= rank && rank <= tMat.WindowSize);
Contracts.Assert(Utils.Size(singularVectors) >= tMat.WindowSize * rank);
int len = 2 * tMat.WindowSize - 1;
Contracts.Assert(Utils.Size(output) >= len);
Single v;
/*var k = tMat.SeriesLength - tMat.WindowSize + 1;
Contracts.Assert(k > 0);
var v = new Single[k];*/
Single temp;
int i;
int j;
int offset1 = tMat.SeriesLength - 2 * tMat.WindowSize + 1;
int offset2 = tMat.SeriesLength - tMat.WindowSize;
for (i = 0; i < len; ++i)
output[i] = 0;
for (i = 0; i < rank; ++i)
{
// Reconstructing the i-th eigen triple component and adding it to output
v = 0;
for (j = offset1; j < offset1 + tMat.WindowSize; ++j)
v += series[j] * singularVectors[tMat.WindowSize * i - offset1 + j];
for (j = 0; j < tMat.WindowSize - 1; ++j)
output[j] += v * singularVectors[tMat.WindowSize * i + j];
temp = v * singularVectors[tMat.WindowSize * (i + 1) - 1];
v = 0;
for (j = offset2; j < offset2 + tMat.WindowSize; ++j)
v += series[j] * singularVectors[tMat.WindowSize * i - offset2 + j];
for (j = tMat.WindowSize; j < 2 * tMat.WindowSize - 1; ++j)
output[j] += v * singularVectors[tMat.WindowSize * (i - 1) + j + 1];
temp += v * singularVectors[tMat.WindowSize * i];
output[tMat.WindowSize - 1] += (temp / 2);
}
}
private static void ComputeNoiseMoments(Single[] series, Single[] signal, Single[] alpha, out Single observationNoiseVariance, out Single autoregressionNoiseVariance,
out Single observationNoiseMean, out Single autoregressionNoiseMean, int startIndex = 0)
{
Contracts.Assert(Utils.Size(alpha) > 0);
Contracts.Assert(Utils.Size(signal) > 2 * Utils.Size(alpha)); // To assure that the autoregression noise variance is unbiased.
Contracts.Assert(Utils.Size(series) >= Utils.Size(signal) + startIndex);
var signalLength = Utils.Size(signal);
var windowSize = Utils.Size(alpha) + 1;
var k = signalLength - windowSize + 1;
Contracts.Assert(k > 0);
var y = new Single[k];
int i;
observationNoiseMean = 0;
observationNoiseVariance = 0;
autoregressionNoiseMean = 0;
autoregressionNoiseVariance = 0;
// Computing the observation noise moments
for (i = 0; i < signalLength; ++i)
observationNoiseMean += (series[i + startIndex] - signal[i]);
observationNoiseMean /= signalLength;
for (i = 0; i < signalLength; ++i)
observationNoiseVariance += (series[i + startIndex] - signal[i]) * (series[i + startIndex] - signal[i]);
observationNoiseVariance /= signalLength;
observationNoiseVariance -= (observationNoiseMean * observationNoiseMean);
// Computing the auto-regression noise moments
TrajectoryMatrix xTM = new TrajectoryMatrix(null, signal, windowSize - 1, signalLength - 1);
xTM.MultiplyTranspose(alpha, y);
for (i = 0; i < k; ++i)
autoregressionNoiseMean += (signal[windowSize - 1 + i] - y[i]);
autoregressionNoiseMean /= k;
for (i = 0; i < k; ++i)
{
autoregressionNoiseVariance += (signal[windowSize - 1 + i] - y[i] - autoregressionNoiseMean) *
(signal[windowSize - 1 + i] - y[i] - autoregressionNoiseMean);
}
autoregressionNoiseVariance /= (k - windowSize + 1);
Contracts.Assert(autoregressionNoiseVariance >= 0);
}
private static int DetermineSignalRank(Single[] series, TrajectoryMatrix tMat, Single[] singularVectors, Single[] singularValues,
Single[] outputSignal, int maxRank)
{
Contracts.Assert(tMat != null);
Contracts.Assert(Utils.Size(series) >= tMat.SeriesLength);
Contracts.Assert(Utils.Size(outputSignal) >= tMat.SeriesLength);
Contracts.Assert(Utils.Size(singularVectors) >= tMat.WindowSize * tMat.WindowSize);
Contracts.Assert(Utils.Size(singularValues) >= tMat.WindowSize);
Contracts.Assert(1 <= maxRank && maxRank <= tMat.WindowSize - 1);
var inputSeriesLength = tMat.SeriesLength;
var k = inputSeriesLength - tMat.WindowSize + 1;
Contracts.Assert(k > 0);
var x = new Single[inputSeriesLength];
var y = new Single[k];
var alpha = new Single[tMat.WindowSize - 1];
var v = new Single[k];
Single nu = 0;
Double minErr = Double.MaxValue;
int minIndex = maxRank - 1;
int evaluationLength = Math.Min(Math.Max(tMat.WindowSize, 200), k);
TrajectoryMatrix xTM = new TrajectoryMatrix(null, x, tMat.WindowSize - 1, inputSeriesLength - 1);
int i;
int j;
int n;
Single temp;
Double error;
Double sumSingularVals = 0;
Single lambda;
Single observationNoiseMean;
FixedSizeQueue<Single> window = new FixedSizeQueue<float>(tMat.WindowSize - 1);
for (i = 0; i < tMat.WindowSize; ++i)
sumSingularVals += singularValues[i];
for (i = 0; i < maxRank; ++i)
{
// Updating the auto-regressive coefficients
lambda = singularVectors[tMat.WindowSize * i + tMat.WindowSize - 1];
for (j = 0; j < tMat.WindowSize - 1; ++j)
alpha[j] += lambda * singularVectors[tMat.WindowSize * i + j];
// Updating nu
nu += lambda * lambda;
// Reconstructing the i-th eigen triple component and adding it to x
tMat.MultiplyTranspose(singularVectors, v, false, tMat.WindowSize * i, 0);
tMat.RankOneHankelization(singularVectors, v, 1, x, true, tMat.WindowSize * i, 0, 0);
observationNoiseMean = 0;
for (j = 0; j < inputSeriesLength; ++j)
observationNoiseMean += (series[j] - x[j]);
observationNoiseMean /= inputSeriesLength;
for (j = inputSeriesLength - evaluationLength - tMat.WindowSize + 1; j < inputSeriesLength - evaluationLength; ++j)
window.AddLast(x[j]);
error = 0;
for (j = inputSeriesLength - evaluationLength; j < inputSeriesLength; ++j)
{
temp = 0;
for (n = 0; n < tMat.WindowSize - 1; ++n)
temp += alpha[n] * window[n];
temp /= (1 - nu);
temp += observationNoiseMean;
window.AddLast(temp);
error += Math.Abs(series[j] - temp);
if (error > minErr)
break;
}
if (error < minErr)
{
minErr = error;
minIndex = i;
Array.Copy(x, outputSignal, inputSeriesLength);
}
}
return minIndex + 1;
}
internal override void InitState()
{
for (int i = 0; i < _windowSize - 2; ++i)
_state[i] = 0;
_buffer.Clear();
}
private static int DetermineSignalRankFast(Single[] series, TrajectoryMatrix tMat, Single[] singularVectors, Single[] singularValues, int maxRank)
{
Contracts.Assert(tMat != null);
Contracts.Assert(Utils.Size(series) >= tMat.SeriesLength);
Contracts.Assert(Utils.Size(singularVectors) >= tMat.WindowSize * tMat.WindowSize);
Contracts.Assert(Utils.Size(singularValues) >= tMat.WindowSize);
Contracts.Assert(1 <= maxRank && maxRank <= tMat.WindowSize - 1);
var inputSeriesLength = tMat.SeriesLength;
var k = inputSeriesLength - tMat.WindowSize + 1;
Contracts.Assert(k > 0);
var x = new Single[tMat.WindowSize - 1];
var alpha = new Single[tMat.WindowSize - 1];
Single v;
Single nu = 0;
Double minErr = Double.MaxValue;
int minIndex = maxRank - 1;
int evaluationLength = Math.Min(Math.Max(tMat.WindowSize, 200), k);
int i;
int j;
int n;
int offset;
Single temp;
Double error;
Single lambda;
Single observationNoiseMean;
FixedSizeQueue<Single> window = new FixedSizeQueue<float>(tMat.WindowSize - 1);
for (i = 0; i < maxRank; ++i)
{
// Updating the auto-regressive coefficients
lambda = singularVectors[tMat.WindowSize * i + tMat.WindowSize - 1];
for (j = 0; j < tMat.WindowSize - 1; ++j)
alpha[j] += lambda * singularVectors[tMat.WindowSize * i + j];
// Updating nu
nu += lambda * lambda;
// Reconstructing the i-th eigen triple component and adding it to x
v = 0;
offset = inputSeriesLength - evaluationLength - tMat.WindowSize + 1;
for (j = offset; j <= inputSeriesLength - evaluationLength; ++j)
v += series[j] * singularVectors[tMat.WindowSize * i - offset + j];
for (j = 0; j < tMat.WindowSize - 1; ++j)
x[j] += v * singularVectors[tMat.WindowSize * i + j];
// Computing the empirical observation noise mean
observationNoiseMean = 0;
for (j = offset; j < inputSeriesLength - evaluationLength; ++j)
observationNoiseMean += (series[j] - x[j - offset]);
observationNoiseMean /= (tMat.WindowSize - 1);
for (j = 0; j < tMat.WindowSize - 1; ++j)
window.AddLast(x[j]);
error = 0;
for (j = inputSeriesLength - evaluationLength; j < inputSeriesLength; ++j)
{
temp = 0;
for (n = 0; n < tMat.WindowSize - 1; ++n)
temp += alpha[n] * window[n];
temp /= (1 - nu);
temp += observationNoiseMean;
window.AddLast(temp);
error += Math.Abs(series[j] - temp);
if (error > minErr)
break;
}
if (error < minErr)
{
minErr = error;
minIndex = i;
}
}
return minIndex + 1;
}
private class SignalComponent
{
public Double Phase;
public int Index;
public int Cluster;
public SignalComponent(Double phase, int index)
{
Phase = phase;
Index = index;
}
}
private bool Stabilize()
{
if (Utils.Size(_alpha) == 1)
{
if (_shouldMaintainInfo)
_info.RootsBeforeStabilization = new[] { new Complex(_alpha[0], 0) };
if (_alpha[0] > 1)
_alpha[0] = 1;
else if (_alpha[0] < -1)
_alpha[0] = -1;
if (_shouldMaintainInfo)
{
_info.IsStabilized = true;
_info.RootsAfterStabilization = new[] { new Complex(_alpha[0], 0) };
_info.IsExponentialTrendPresent = false;
_info.IsPolynomialTrendPresent = false;
_info.ExponentialTrendFactor = Math.Abs(_alpha[0]);
}
return true;
}
var coeff = new Double[_windowSize - 1];
Complex[] roots = null;
bool trendFound = false;
bool polynomialTrendFound = false;
Double maxTrendMagnitude = Double.NegativeInfinity;
Double maxNonTrendMagnitude = Double.NegativeInfinity;
Double eps = 1e-9;
Double highFrequenceyBoundry = Math.PI / 2;
var sortedComponents = new List<SignalComponent>();
var trendComponents = new List<int>();
int i;
// Computing the roots of the characteristic polynomial
for (i = 0; i < _windowSize - 1; ++i)
coeff[i] = -_alpha[i];
if (!PolynomialUtils.FindPolynomialRoots(coeff, ref roots))
return false;
if (_shouldMaintainInfo)
{
_info.RootsBeforeStabilization = new Complex[_windowSize - 1];
Array.Copy(roots, _info.RootsBeforeStabilization, _windowSize - 1);
}
// Detecting trend components
for (i = 0; i < _windowSize - 1; ++i)
{
if (roots[i].Magnitude > 1 && (Math.Abs(roots[i].Phase) <= eps || Math.Abs(Math.Abs(roots[i].Phase) - Math.PI) <= eps))
trendComponents.Add(i);
}
// Removing unobserved seasonalities and consequently introducing polynomial trend
for (i = 0; i < _windowSize - 1; ++i)
{
if (roots[i].Phase != 0)
{
if (roots[i].Magnitude >= 1 && 2 * Math.PI / Math.Abs(roots[i].Phase) > _windowSize)
{
/*if (roots[i].Real > 1)
{
polynomialTrendFound = true;
roots[i] = new Complex(Math.Max(1, roots[i].Magnitude), 0);
maxPolynomialTrendMagnitude = Math.Max(maxPolynomialTrendMagnitude, roots[i].Magnitude);
}
else
roots[i] = Complex.FromPolarCoordinates(1, 0);
//roots[i] = Complex.FromPolarCoordinates(0.99, roots[i].Phase);*/
/* if (_maxTrendRatio > 1)
{
roots[i] = new Complex(roots[i].Real, 0);
polynomialTrendFound = true;
}
else
roots[i] = roots[i].Imaginary > 0 ? new Complex(roots[i].Real, 0) : new Complex(1, 0);*/
roots[i] = new Complex(roots[i].Real, 0);
polynomialTrendFound = true;
if (_shouldMaintainInfo)
_info.IsArtificialSeasonalityRemoved = true;
}
else if (roots[i].Magnitude > 1)
sortedComponents.Add(new SignalComponent(roots[i].Phase, i));
}
}
if (_maxTrendRatio > 1)
{
// Combining the close exponential-seasonal components
if (sortedComponents.Count > 1 && polynomialTrendFound)
{
sortedComponents.Sort((a, b) => a.Phase.CompareTo(b.Phase));
var clusterNum = 0;
for (i = 0; i < sortedComponents.Count - 1; ++i)
{
if ((sortedComponents[i].Phase < 0 && sortedComponents[i + 1].Phase < 0) ||
(sortedComponents[i].Phase > 0 && sortedComponents[i + 1].Phase > 0))
{
sortedComponents[i].Cluster = clusterNum;
if (Math.Abs(sortedComponents[i + 1].Phase - sortedComponents[i].Phase) > 0.05)
clusterNum++;
sortedComponents[i + 1].Cluster = clusterNum;
}
else
clusterNum++;
}
int start = 0;
bool polynomialSeasonalityFound = false;
Double largestSeasonalityPhase = 0;
for (i = 1; i < sortedComponents.Count; ++i)
{
if (sortedComponents[i].Cluster != sortedComponents[i - 1].Cluster)
{
if (i - start > 1) // There are more than one point in the cluster
{
Double avgPhase = 0;
Double avgMagnitude = 0;
for (var j = start; j < i; ++j)
{
avgPhase += sortedComponents[j].Phase;
avgMagnitude += roots[sortedComponents[j].Index].Magnitude;
}
avgPhase /= (i - start);
avgMagnitude /= (i - start);
for (var j = start; j < i; ++j)
roots[sortedComponents[j].Index] = Complex.FromPolarCoordinates(avgMagnitude,
avgPhase);
if (!polynomialSeasonalityFound && avgPhase > 0)
{
largestSeasonalityPhase = avgPhase;
polynomialSeasonalityFound = true;
}
}
start = i;
}
}
}
// Combining multiple exponential trends into polynomial ones
if (!polynomialTrendFound)
{
var ind1 = -1;
var ind2 = -1;
foreach (var ind in trendComponents)
{
if (Math.Abs(roots[ind].Phase) <= eps)
{
ind1 = ind;
break;
}
}
for (i = 0; i < _windowSize - 1; ++i)
{
if (Math.Abs(roots[i].Phase) <= eps && 0.9 <= roots[i].Magnitude && i != ind1)
{
ind2 = i;
break;
}
}
if (ind1 >= 0 && ind2 >= 0 && ind1 != ind2)
{
roots[ind1] = Complex.FromPolarCoordinates(1, 0);
roots[ind2] = Complex.FromPolarCoordinates(1, 0);
polynomialTrendFound = true;
}
}
}
if (polynomialTrendFound) // Suppress the exponential trend
{
maxTrendMagnitude = Math.Min(1, _maxTrendRatio);
foreach (var ind in trendComponents)
roots[ind] = Complex.FromPolarCoordinates(0.99, roots[ind].Phase);
}
else
{
// Spotting the exponential trend components and finding the max trend magnitude
for (i = 0; i < _windowSize - 1; ++i)
{
if (roots[i].Magnitude > 1 && Math.Abs(roots[i].Phase) <= eps)
{
trendFound = true;
maxTrendMagnitude = Math.Max(maxTrendMagnitude, roots[i].Magnitude);
}
else
maxNonTrendMagnitude = Math.Max(maxNonTrendMagnitude, roots[i].Magnitude);
}
if (!trendFound)
maxTrendMagnitude = 1;
maxTrendMagnitude = Math.Min(maxTrendMagnitude, _maxTrendRatio);
}
// Squeezing all components below the maximum trend magnitude
var smallTrendMagnitude = Math.Min(maxTrendMagnitude, (maxTrendMagnitude + 1) / 2);
for (i = 0; i < _windowSize - 1; ++i)
{
if (roots[i].Magnitude >= maxTrendMagnitude)
{
if ((highFrequenceyBoundry < roots[i].Phase && roots[i].Phase < Math.PI - eps) ||
(-Math.PI + eps < roots[i].Phase && roots[i].Phase < -highFrequenceyBoundry))
roots[i] = Complex.FromPolarCoordinates(smallTrendMagnitude, roots[i].Phase);
else
roots[i] = Complex.FromPolarCoordinates(maxTrendMagnitude, roots[i].Phase);
}
}
// Correcting all the other trend components
for (i = 0; i < _windowSize - 1; ++i)
{
var phase = roots[i].Phase;
if (Math.Abs(phase) <= eps)
roots[i] = new Complex(roots[i].Magnitude, 0);
else if (Math.Abs(phase - Math.PI) <= eps)
roots[i] = new Complex(-roots[i].Magnitude, 0);
else if (Math.Abs(phase + Math.PI) <= eps)
roots[i] = new Complex(-roots[i].Magnitude, 0);
}
// Computing the characteristic polynomial from the modified roots
try
{
if (!PolynomialUtils.FindPolynomialCoefficients(roots, ref coeff))
return false;
}
catch (OverflowException)
{
return false;
}
// Updating alpha
for (i = 0; i < _windowSize - 1; ++i)
_alpha[i] = (Single)(-coeff[i]);
if (_shouldMaintainInfo)
{
_info.RootsAfterStabilization = roots;
_info.IsStabilized = true;
_info.IsPolynomialTrendPresent = polynomialTrendFound;
_info.IsExponentialTrendPresent = maxTrendMagnitude > 1;
_info.ExponentialTrendFactor = maxTrendMagnitude;
}
return true;
}
/// <summary>
/// Feeds the next observation on the series to the model and as a result changes the state of the model.
/// </summary>
/// <param name="input">The next observation on the series.</param>
/// <param name="updateModel">Determines whether the model parameters also need to be updated upon consuming the new observation (default = false).</param>
internal override void Consume(ref Single input, bool updateModel = false)
{
if (Single.IsNaN(input))
return;
int i;
if (_wTrans == null)
{
_y = new CpuAlignedVector(_rank, CpuMathUtils.GetVectorAlignment());
_wTrans = new CpuAlignedMatrixRow(_rank, _windowSize, CpuMathUtils.GetVectorAlignment());
Single[] vecs = new Single[_rank * _windowSize];
for (i = 0; i < _rank; ++i)
vecs[(_windowSize + 1) * i] = 1;
i = 0;
_wTrans.CopyFrom(vecs, ref i);
}
// Forming vector x
if (_buffer.Count == 0)
{
for (i = 0; i < _windowSize - 1; ++i)
_buffer.AddLast(_state[i]);
}
int len = _buffer.Count;
for (i = 0; i < _windowSize - len - 1; ++i)
_x[i] = 0;
for (i = Math.Max(0, len - _windowSize + 1); i < len; ++i)
_x[i - len + _windowSize - 1] = _buffer[i];
_x[_windowSize - 1] = input;
// Computing y: Eq. (11) in https://hal-institut-mines-telecom.archives-ouvertes.fr/hal-00479772/file/twocolumns.pdf
CpuAligenedMathUtils<CpuAlignedMatrixRow>.MatTimesSrc(_wTrans, _x, _y);
// Updating the state vector
CpuAligenedMathUtils<CpuAlignedMatrixRow>.MatTranTimesSrc(_wTrans, _y, _xSmooth);
_nextPrediction = _autoregressionNoiseMean + _observationNoiseMean;
for (i = 0; i < _windowSize - 2; ++i)
{
_state[i] = ((_windowSize - 2 - i) * _state[i + 1] + _xSmooth[i + 1]) / (_windowSize - 1 - i);
_nextPrediction += _state[i] * _alpha[i];
}
_state[_windowSize - 2] = _xSmooth[_windowSize - 1];
_nextPrediction += _state[_windowSize - 2] * _alpha[_windowSize - 2];
if (updateModel)
{
// REVIEW: to be implemented in the next version based on the FAPI algorithm
// in https://hal-institut-mines-telecom.archives-ouvertes.fr/hal-00479772/file/twocolumns.pdf.
}
_buffer.AddLast(input);
}
/// <summary>
/// Train the model parameters based on a training series.
/// </summary>
/// <param name="data">The training time-series.</param>
internal override void Train(FixedSizeQueue<Single> data)
{
_host.CheckParam(data != null, nameof(data), "The input series for training cannot be null.");
_host.CheckParam(data.Count >= _trainSize, nameof(data), "The input series for training does not have enough points for training.");
Single[] dataArray = new Single[_trainSize];
int i;
int count;
for (i = 0, count = 0; count < _trainSize && i < data.Count; ++i)
if (!Single.IsNaN(data[i]))
dataArray[count++] = data[i];
if (_shouldMaintainInfo)
{
_info = new ModelInfo();
_info.WindowSize = _windowSize;
}
if (count <= 2 * _windowSize)
{
#if !TLCSSA
using (var ch = _host.Start("Train"))
ch.Warning(
"Training cannot be completed because the input series for training does not have enough points.");
#endif
}
else
{
if (count != _trainSize)
Array.Resize(ref dataArray, count);
TrainCore(dataArray, count);
}
}
#if !TLCSSA
/// <summary>
/// Train the model parameters based on a training series.
/// </summary>
/// <param name="data">The training time-series.</param>
internal override void Train(RoleMappedData data)
{
_host.CheckValue(data, nameof(data));
_host.CheckParam(data.Schema.Feature.HasValue, nameof(data), "Must have features column.");
var featureCol = data.Schema.Feature.Value;
if (featureCol.Type != NumberDataViewType.Single)
throw _host.ExceptSchemaMismatch(nameof(data), "feature", featureCol.Name, "Single", featureCol.Type.ToString());
Single[] dataArray = new Single[_trainSize];
int count = 0;
using (var cursor = data.Data.GetRowCursor(featureCol))
{
var getVal = cursor.GetGetter<Single>(featureCol);
Single val = default;
while (cursor.MoveNext() && count < _trainSize)
{
getVal(ref val);
if (!Single.IsNaN(val))
dataArray[count++] = val;
}
}
if (_shouldMaintainInfo)
{
_info = new ModelInfo();
_info.WindowSize = _windowSize;
}
if (count <= 2 * _windowSize)
{
using (var ch = _host.Start("Train"))
ch.Warning("Training cannot be completed because the input series for training does not have enough points.");
}
else
{
if (count != _trainSize)
Array.Resize(ref dataArray, count);
TrainCore(dataArray, count);
}
}
#endif
private void TrainCore(Single[] dataArray, int originalSeriesLength)
{
_host.Assert(Utils.Size(dataArray) > 0);
Single[] singularVals;
Single[] leftSingularVecs;
var learnNaiveModel = false;
var signalLength = _rankSelectionMethod == RankSelectionMethod.Exact ? originalSeriesLength : 2 * _windowSize - 1;//originalSeriesLength;
var signal = new Single[signalLength];
int i;
// Creating the trajectory matrix for the series
TrajectoryMatrix tMat = new TrajectoryMatrix(_host, dataArray, _windowSize, originalSeriesLength);
// Computing the SVD of the trajectory matrix
if (!tMat.ComputeSvd(out singularVals, out leftSingularVecs))
learnNaiveModel = true;
else
{
for (i = 0; i < _windowSize * _maxRank; ++i)
{
if (Single.IsNaN(leftSingularVecs[i]))
{
learnNaiveModel = true;
break;
}
}
}
// Checking for standard eigenvectors, if found reduce the window size and reset training.
if (!learnNaiveModel)
{
for (i = 0; i < _windowSize; ++i)
{
var v = leftSingularVecs[(i + 1) * _windowSize - 1];
if (v * v == 1)
{
if (_windowSize > 2)
{
_windowSize--;
_maxRank = _windowSize / 2;
_alpha = new Single[_windowSize - 1];
_state = new Single[_windowSize - 1];
_x = new CpuAlignedVector(_windowSize, CpuMathUtils.GetVectorAlignment());
_xSmooth = new CpuAlignedVector(_windowSize, CpuMathUtils.GetVectorAlignment());
TrainCore(dataArray, originalSeriesLength);
return;
}
else
{
learnNaiveModel = true;
break;
}
}
}
}
// Learn the naive (averaging) model in case the eigen decomposition is not possible
if (learnNaiveModel)
{
#if !TLCSSA
using (var ch = _host.Start("Train"))
ch.Warning("The precise SSA model cannot be trained.");
#endif
_rank = 1;
var temp = (Single)(1f / Math.Sqrt(_windowSize));
for (i = 0; i < _windowSize; ++i)
leftSingularVecs[i] = temp;
}
else
{
// Computing the signal rank
if (_rankSelectionMethod == RankSelectionMethod.Exact)
_rank = DetermineSignalRank(dataArray, tMat, leftSingularVecs, singularVals, signal, _maxRank);
else if (_rankSelectionMethod == RankSelectionMethod.Fast)
_rank = DetermineSignalRankFast(dataArray, tMat, leftSingularVecs, singularVals, _maxRank);
}
// Setting the the y vector
_y = new CpuAlignedVector(_rank, CpuMathUtils.GetVectorAlignment());
// Setting the weight matrix
_wTrans = new CpuAlignedMatrixRow(_rank, _windowSize, CpuMathUtils.GetVectorAlignment());
i = 0;
_wTrans.CopyFrom(leftSingularVecs, ref i);
// Setting alpha
Single nu = 0;
for (i = 0; i < _rank; ++i)
{
_y[i] = leftSingularVecs[_windowSize * (i + 1) - 1];
nu += _y[i] * _y[i];
}
CpuAligenedMathUtils<CpuAlignedMatrixRow>.MatTranTimesSrc(_wTrans, _y, _xSmooth);
for (i = 0; i < _windowSize - 1; ++i)
_alpha[i] = _xSmooth[i] / (1 - nu);
// Stabilizing the model
if (_shouldStablize && !learnNaiveModel)
{
if (!Stabilize())
{
#if !TLCSSA
using (var ch = _host.Start("Train"))
ch.Warning("The trained model cannot be stablized.");
#endif
}
}
// Computing the noise moments
if (ShouldComputeForecastIntervals)
{
if (_rankSelectionMethod != RankSelectionMethod.Exact)
ReconstructSignalTailFast(dataArray, tMat, leftSingularVecs, _rank, signal);
ComputeNoiseMoments(dataArray, signal, _alpha, out _observationNoiseVariance, out _autoregressionNoiseVariance,
out _observationNoiseMean, out _autoregressionNoiseMean, originalSeriesLength - signalLength);
_observationNoiseMean = 0;
_autoregressionNoiseMean = 0;
}
// Setting the state
_nextPrediction = _autoregressionNoiseMean + _observationNoiseMean;
if (_buffer.Count > 0) // Use the buffer to set the state when there are data points pushed into the buffer using the Consume() method
{
int len = _buffer.Count;
for (i = 0; i < _windowSize - len; ++i)
_x[i] = 0;
for (i = Math.Max(0, len - _windowSize); i < len; ++i)
_x[i - len + _windowSize] = _buffer[i];
}
else // use the training data points otherwise
{
for (i = originalSeriesLength - _windowSize; i < originalSeriesLength; ++i)
_x[i - originalSeriesLength + _windowSize] = dataArray[i];
}
CpuAligenedMathUtils<CpuAlignedMatrixRow>.MatTimesSrc(_wTrans, _x, _y);
CpuAligenedMathUtils<CpuAlignedMatrixRow>.MatTranTimesSrc(_wTrans, _y, _xSmooth);
for (i = 1; i < _windowSize; ++i)
{
_state[i - 1] = _xSmooth[i];
_nextPrediction += _state[i - 1] * _alpha[i - 1];
}
if (_shouldMaintainInfo)
{
_info.IsTrained = true;
_info.WindowSize = _windowSize;
_info.AutoRegressiveCoefficients = new Single[_windowSize - 1];
Array.Copy(_alpha, _info.AutoRegressiveCoefficients, _windowSize - 1);
_info.Rank = _rank;
_info.IsNaiveModelTrained = learnNaiveModel;
_info.Spectrum = singularVals;
}
}
/// <summary>
/// Forecasts the future values of the series up to the given horizon.
/// </summary>
/// <param name="result">The forecast result.</param>
/// <param name="horizon">The forecast horizon.</param>
internal override void Forecast(ref ForecastResultBase<Single> result, int horizon = 1)
{
_host.CheckParam(horizon >= 1, nameof(horizon), "The horizon parameter should be greater than 0.");
if (result == null)
result = new SsaForecastResult();
var str = "The result argument must be of type " + typeof(SsaForecastResult).ToString();
_host.CheckParam(result is SsaForecastResult, nameof(result), str);
var output = result as SsaForecastResult;
var resEditor = VBufferEditor.Create(ref result.PointForecast, horizon);
int i;
int j;
int k;
// Computing the point forecasts
resEditor.Values[0] = _nextPrediction;
for (i = 1; i < horizon; ++i)
{
k = 0;
resEditor.Values[i] = _autoregressionNoiseMean + _observationNoiseMean;
for (j = i; j < _windowSize - 1; ++j, ++k)
resEditor.Values[i] += _state[j] * _alpha[k];
for (j = Math.Max(0, i - _windowSize + 1); j < i; ++j, ++k)
resEditor.Values[i] += resEditor.Values[j] * _alpha[k];
}
// Computing the forecast variances
if (ShouldComputeForecastIntervals)
{
var sdEditor = VBufferEditor.Create(ref output.ForecastStandardDeviation, horizon);
var lastCol = new FixedSizeQueue<Single>(_windowSize - 1);
for (i = 0; i < _windowSize - 3; ++i)
lastCol.AddLast(0);
lastCol.AddLast(1);
lastCol.AddLast(_alpha[_windowSize - 2]);
sdEditor.Values[0] = _autoregressionNoiseVariance + _observationNoiseVariance;
for (i = 1; i < horizon; ++i)
{
Single temp = 0;
for (j = 0; j < _windowSize - 1; ++j)
temp += _alpha[j] * lastCol[j];
lastCol.AddLast(temp);
sdEditor.Values[i] = sdEditor.Values[i - 1] + _autoregressionNoiseVariance * temp * temp;
}
for (i = 0; i < horizon; ++i)
sdEditor.Values[i] = (float)Math.Sqrt(sdEditor.Values[i]);
output.ForecastStandardDeviation = sdEditor.Commit();
}
result.PointForecast = resEditor.Commit();
output.CanComputeForecastIntervals = ShouldComputeForecastIntervals;
output.BoundOffset = 0;
}
/// <summary>
/// Predicts the next value on the series.
/// </summary>
/// <param name="output">The prediction result.</param>
internal override void PredictNext(ref Single output)
{
output = _nextPrediction;
}
internal override SequenceModelerBase<Single, Single> Clone()
{
return new AdaptiveSingularSpectrumSequenceModelerInternal(this);
}
/// <summary>
/// Computes the forecast intervals for the input forecast object at the given confidence level. The results are stored in the forecast object.
/// </summary>
/// <param name="forecast">The input forecast object</param>
/// <param name="confidenceLevel">The confidence level in [0, 1)</param>
internal static void ComputeForecastIntervals(ref SsaForecastResult forecast, Single confidenceLevel = 0.95f)
{
Contracts.CheckParam(0 <= confidenceLevel && confidenceLevel < 1, nameof(confidenceLevel), "The confidence level must be in [0, 1).");
Contracts.CheckValue(forecast, nameof(forecast));
Contracts.Check(forecast.CanComputeForecastIntervals, "The forecast intervals cannot be computed for this forecast object.");
var meanForecast = forecast.PointForecast.GetValues();
var horizon = meanForecast.Length;
var sdForecast = forecast.ForecastStandardDeviation.GetValues();
Contracts.Check(sdForecast.Length >= horizon, "The forecast standard deviation values are not available.");
forecast.ConfidenceLevel = confidenceLevel;
if (horizon == 0)
return;
var upper = VBufferEditor.Create(ref forecast.UpperBound, horizon);
var lower = VBufferEditor.Create(ref forecast.LowerBound, horizon);
var z = ProbabilityFunctions.Probit(0.5 + confidenceLevel / 2.0);
double temp;
for (int i = 0; i < horizon; ++i)
{
temp = z * sdForecast[i];
upper.Values[i] = (Single)(meanForecast[i] + forecast.BoundOffset + temp);
lower.Values[i] = (Single)(meanForecast[i] + forecast.BoundOffset - temp);
}
forecast.UpperBound = upper.Commit();
forecast.LowerBound = lower.Commit();
}
public void Train(IDataView dataView, string inputColumnName) => Train(new RoleMappedData(dataView, null, inputColumnName));
public float[] Forecast(int horizon)
{
ForecastResultBase<float> result = null;
Forecast(ref result, horizon);
return result.PointForecast.GetValues().ToArray();
}
public void ForecastWithConfidenceIntervals(int horizon, out float[] forecast, out float[] confidenceIntervalLowerBounds, out float[] confidenceIntervalUpperBounds, float confidenceLevel = 0.95f)
{
ForecastResultBase<float> result = null;
Forecast(ref result, horizon);
SsaForecastResult ssaResult = (SsaForecastResult)result;
ComputeForecastIntervals(ref ssaResult, confidenceLevel);
forecast = result.PointForecast.GetValues().ToArray();
confidenceIntervalLowerBounds = ssaResult.LowerBound.GetValues().ToArray();
confidenceIntervalUpperBounds = ssaResult.UpperBound.GetValues().ToArray();
}
public void Update(IDataView dataView, string inputColumnName)
{
_host.CheckParam(dataView != null, nameof(dataView), "The input series for updating cannot be null.");
var data = new RoleMappedData(dataView, null, inputColumnName);
if (data.Schema.Feature.Value.Type != NumberDataViewType.Single)
throw _host.ExceptUserArg(nameof(data.Schema.Feature.Value.Name), "The time series input column has " +
"type '{0}', but must be a float.", data.Schema.Feature.Value.Type);
var col = data.Schema.Feature.Value;
using (var cursor = data.Data.GetRowCursor(data.Data.Schema))
{
var getVal = cursor.GetGetter<Single>(col);
Single val = default(Single);
while (cursor.MoveNext())
{
getVal(ref val);
if (!Single.IsNaN(val))
Consume(ref val);
}
}
}
}
}
|