File: AdaptiveSingularSpectrumSequenceModeler.cs
Web Access
Project: src\src\Microsoft.ML.TimeSeries\Microsoft.ML.TimeSeries.csproj (Microsoft.ML.TimeSeries)
// 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);
                }
            }
        }
    }
}