classdef incrementalPCA < matlab.mixin.internal.Scalar ...
        & matlab.mixin.CustomDisplay & matlab.mixin.CustomCompactDisplayProvider
    %

    %   Copyright 2023-2024 The MathWorks, Inc.

    properties (GetAccess = public, SetAccess = protected)
        IsWarm = false
        NumTrainingObservations(1,1) double = 0
        WarmupPeriod(1,1) {mustBeInteger,mustBeNumeric,mustBeNonnegative, mustBeNonsparse} = 0
    end

    properties (GetAccess=public, SetAccess=protected, Dependent)
        Mu
        Sigma
        ExplainedVariance
        EstimationPeriod
        Latent
        Coefficients
        VariableWeights
        NumComponents
        NumPredictors
    end

    properties (Access = {?matlab.unittest.TestCase})
        Impl = [];
    end

    properties (Access=protected)
        PredictorClass(1,:) char = 'double'
        PrivPredictorClassFixed = false
        DOF
    end

    methods
        function obj = incrementalPCA(varargin)

            [varargin{:}] = convertStringsToChars(varargin{:});
            internal.stats.checkNotTall(upper(mfilename),0,varargin{:});

            [modelParams,implParams,parametersAreFromOffline,remainingArgs] = incrementalPCA.extractPCAParameters(mfilename,varargin{:});

            obj.WarmupPeriod = modelParams.WarmupPeriod;
            obj.WarmupPeriod = double(obj.WarmupPeriod);

            validateattributes(implParams.NumPredictors,{'numeric'},{'scalar','integer','nonnegative','nonsparse'},mfilename,'NumPredictors');
            validateattributes(implParams.NumComponents,{'numeric'},{'scalar','integer','nonnegative','nonsparse'},mfilename,'NumComponents');

            implParams.NumPredictors = double(implParams.NumPredictors);
            implParams.NumComponents = double(implParams.NumComponents);

            obj.PredictorClass = modelParams.PredictorClass;
            if parametersAreFromOffline
                obj.PrivPredictorClassFixed = true;
            end


            if ~isempty(implParams.VariableWeights)
                if 0 < implParams.NumPredictors
                    validateattributes(implParams.VariableWeights,{'double','single'},{'vector','numel',implParams.NumPredictors,'finite','real','positive'},mfilename,'VariableWeights');
                else
                    validateattributes(implParams.VariableWeights,{'double','single'},{'vector','finite','real','positive'},mfilename,'VariableWeights');
                    implParams.NumPredictors = numel(implParams.VariableWeights);
                end
                implParams.VariableWeights = implParams.VariableWeights(:)';
            end

            if 0 < implParams.NumPredictors
                validateattributes(implParams.NumComponents,{'double','single'},{'scalar','<=',implParams.NumPredictors},mfilename,'NumComponents');
            end

            if implParams.NumComponents == 0
                implParams.NumComponents = implParams.NumPredictors;
            end

            if parametersAreFromOffline
                if isempty(implParams.VariableWeights)
                    implParams.VariableWeights = ones(1,implParams.NumPredictors,'like',modelParams.Coefficients);
                end
                PhiSqrt = sqrt(implParams.VariableWeights);
                coeff = modelParams.Coefficients;
                % validate orthonormality within eps
                R = abs(coeff'*diag(implParams.VariableWeights)*coeff);
                epsV = 100*max(eps(R),[],"all");
                n = size(R,1);
                if epsV<norm(eye(n)-R,"fro")
                    error(message('stats:incremental:PCA:InvalidCoefficients',n));
                end
                implParams.U = coeff.*PhiSqrt';
                % S will be truncated to have NumComponents elements
                if implParams.OfflineDOF == 0
                    implParams.S = sqrt(modelParams.Latent);
                else
                    implParams.S = sqrt(modelParams.Latent.*implParams.OfflineDOF);
                end
            else
                implParams.U = zeros(implParams.NumPredictors,implParams.NumComponents,obj.PredictorClass);
                % S will be truncated to have NumComponents elements
                implParams.S = zeros(implParams.NumComponents,1,obj.PredictorClass);
                if isempty(implParams.VariableWeights)
                    implParams.VariableWeights = ones(1,implParams.NumPredictors,'like',implParams.U);
                end
            end

            obj.Impl = incremental.impl.PCAImpl.make(implParams,parametersAreFromOffline);
            [obj.Impl,remainingArgs] = configureForEstimation(obj.Impl,parametersAreFromOffline,obj.NumPredictors,obj.PredictorClass,remainingArgs);

            % No parameters should remain
            f = fieldnames(remainingArgs);
            if ~isempty(f)
                error(message('stats:incremental:Learner:UnrecognizedNV',f{1}));
            end

            if ~parametersAreFromOffline
                obj = updateWarmupPeriod(obj);
            end
            obj = updateIsWarmIfNeeded(obj,0);
        end

        function [varargout] = fit(obj,X,varargin)

            [varargin{:}] = convertStringsToChars(varargin{:});
            internal.stats.checkNotTall(upper(mfilename),0,X,varargin{:});
            varargin = incremental.utils.parseUnsupportedParameters('ObservationsIn', varargin{:});
            [W,observationsIn] = obj.parseCommonInputs(varargin{:});
            obsInRows = obj.validateObservationsIn(observationsIn,mfilename);
            [varargout{1:nargout}] = fitPCA(obj,X,W,obsInRows,mfilename);
        end

        function [Xtransformed] = transform(obj,X,varargin)
            if ~canTransform(obj)
                error(message('stats:incremental:PCA:ModelNotFullyDefined'));
            end
            [varargin{:}] = convertStringsToChars(varargin{:});
            internal.stats.checkNotTall(upper(mfilename),0,X,varargin{:});
            varargin = incremental.utils.parseUnsupportedParameters({'ObservationsIn','Weights'}, varargin{:});
            [W,observationsIn] = obj.parseCommonInputs(varargin{:});
            obsInRows = obj.validateObservationsIn(observationsIn,mfilename);
            if obsInRows
                nobsinitial = size(X,1);
            else
                nobsinitial = size(X,2);
            end
            [~,X,~,removedIndices] = prepareData(obj,X,W,obsInRows,false,mfilename);
            [Xtransformed] = transformWithoutValidation(obj,X,[],nobsinitial,removedIndices,obsInRows);

        end

        function varweights = get.VariableWeights(obj)
            varweights = obj.Impl.VariableWeights;
        end

        function numpredictors = get.NumPredictors(obj)
            numpredictors = obj.Impl.NumPredictors;
        end

        function obj = set.NumPredictors(obj,numpredictors)
            obj.Impl.NumPredictors = numpredictors;
        end

        function numcomponents = get.NumComponents(obj)
            numcomponents = obj.Impl.NumComponents;
        end

        function obj = set.NumComponents(obj,numcomponents)
            obj.Impl.NumComponents = numcomponents;
        end

        function latent = get.Latent(obj)

            if obj.DOF == 0
                partiallatent = obj.Impl.S.^2;
            else
                partiallatent = obj.Impl.S.^2./obj.DOF;
            end
            if all(isnan(partiallatent),"all")
                latent = nan(obj.NumComponents,1,'like',partiallatent);
                return;
            end
            n = numel(partiallatent);
            if n < obj.NumComponents
                latent = zeros(obj.NumComponents,1,'like',partiallatent);
                latent(1:n) = partiallatent;
            else
                latent = partiallatent;
            end

        end

        function dof = get.DOF(obj)
            dof = obj.NumTrainingObservations+obj.Impl.DOF;
            dof = max(0,dof);
        end

        function coefficients = get.Coefficients(obj)
            coefficients = obj.Impl.Coefficients;
        end

        function estimationperiod = get.EstimationPeriod(obj)
            estimationperiod = obj.Impl.EstimationPeriod;
        end

        function explainedvariance = get.ExplainedVariance(obj)
            if obj.DOF == 0
                fulllatent = obj.Impl.SSquared;
            else
                fulllatent = obj.Impl.SSquared./obj.DOF;
            end
            p = min(obj.NumComponents,numel(fulllatent));
            sumlatent = sum(fulllatent);
            temp = 100*fulllatent(1:p,:);
            if 0 < sumlatent
                temp = temp/sumlatent;
            end
            if numel(temp) < obj.NumComponents
                if all(isnan(temp),"all")
                    explainedvariance = nan(obj.NumComponents,1,'like',temp);
                else
                    explainedvariance = zeros(obj.NumComponents,1,'like',temp);
                    explainedvariance(1:numel(temp)) = temp;
                end
            else
                explainedvariance = temp;
            end
        end


        function mu = get.Mu(obj)
            mu = anomaly.utils.getMu(obj.Impl.Mu,[]);
        end

        function obj = set.Mu(obj,mu)
            obj.Impl.Mu = mu;
        end

        function sigma = get.Sigma(obj)
            sigma = anomaly.utils.getSigma(obj.Impl.Sigma,[]);
        end

        function obj = reset(obj,varargin)
            narginchk(1,1);

            if obj.WarmupPeriod == 0
                obj.WarmupPeriod = 1000;
            end
            if obj.WarmupPeriod < obj.NumComponents
                obj.WarmupPeriod = obj.NumComponents;
            end
            obj.NumTrainingObservations = 0;
            obj.Impl = reset(obj.Impl,varargin{:});
            obj.IsWarm = false;
            obj = updateIsWarmIfNeeded(obj,obj.NumTrainingObservations);
        end
    end

    methods (Access = protected)

        function [Xtransformed] = transformWithoutValidation(obj,X,W,N,removedIndices,obsInRows)
            scores = score(obj.Impl,X,W,obsInRows);
            if obsInRows
                Xtransformed = nan(N,obj.NumComponents,'like',X);
                Xtransformed(~removedIndices,:) = scores;
            else
                Xtransformed = nan(obj.NumComponents,N,'like',X);
                Xtransformed(:,~removedIndices) = scores;
            end
        end

        function [obj,Xtransformed] = fitPCA(obj,X,W,obsInRows,functionname,varargin)

            if obj.InputIsEmpty(X,obsInRows)
                if nargout == 2
                    Xtransformed = X;
                end
                return;
            end

            if obsInRows
                nobsinitial = size(X,1);
            else
                nobsinitial = size(X,2);
            end

            [obj,X,W,removedIndices,nobs] = prepareForFit(obj,X,W,obsInRows,functionname, varargin{:});

            if canTrain(obj)
                % if model can train, call fit method for model

                % if no incremental training has been done, it's first call
                % to fit.
                firstCalltoFit = obj.NumTrainingObservations == 0;
                obj.Impl = fit(obj.Impl,X,[],W,true,obsInRows,obj.NumTrainingObservations,firstCalltoFit);

                if obj.IsWarm
                    if nargout == 2
                        Xtransformed = transformWithoutValidation(obj,X,[],nobsinitial,removedIndices,obsInRows);
                    end
                else
                    if nargout == 2
                        if obsInRows
                            Xtransformed = nan(nobsinitial,obj.NumComponents,obj.PredictorClass);
                        else
                            Xtransformed = nan(obj.NumComponents,nobsinitial,obj.PredictorClass);
                        end
                        if obj.WarmupPeriod < nobs
                            chunksize = nobs-obj.WarmupPeriod;
                            partialtransform = true;
                        else
                            chunksize = obj.NumTrainingObservations+nobs-obj.WarmupPeriod;
                            partialtransform = 0 < chunksize;
                        end
                        if partialtransform
                            outputindices = (1:nobsinitial)';
                            outputindices = outputindices(~removedIndices);
                            n = numel(outputindices);
                            outputindices = outputindices(n-chunksize+1:n);
                            transformindices = (1:nobs)';
                            n = numel(transformindices);
                            transformindices = transformindices(n-chunksize+1:n);
                            n = numel(transformindices);
                            if obsInRows
                                Xtransformed(outputindices,:) = transformWithoutValidation(obj,X(transformindices,:),[],n,false(n,1),obsInRows);
                            else
                                Xtransformed(:,outputindices) = transformWithoutValidation(obj,X(:,transformindices),[],n,false(n,1),obsInRows);
                            end
                        end
                    end
                end
                obj.NumTrainingObservations = incremental.utils.incrementObservations(obj.NumTrainingObservations,nobs);
            else
                if nargout == 2
                    if obsInRows
                        Xtransformed = nan(nobsinitial,obj.NumComponents,obj.PredictorClass);
                    else
                        Xtransformed = nan(obj.NumComponents,nobsinitial,obj.PredictorClass);
                    end
                end
            end
            obj = updateIsWarmIfNeeded(obj,obj.NumTrainingObservations);
        end

        function [obj,X,W,removedIndices,nobs] = prepareForFit(obj,X,W,obsInRows,functionname,varargin)
            % prepare data for training and update internal states

            [obj,X,W,removedIndices,nobs] = prepareData(obj,X,W,obsInRows,true,functionname, varargin{:});
            obj.Impl = allocateEstimationBuffers(obj.Impl,obj.NumPredictors,obj.PredictorClass);
            [obj.Impl,X,~,W,removedIndices] = prepareDataAndUpdateEstimationState(obj.Impl,X,[],W,obsInRows,nobs,removedIndices);
            obj = synchronizePredictorInformation(obj,X,obsInRows);
            if obsInRows
                nobs = size(X,1);
            else
                nobs = size(X,2);
            end
        end

        function [obj,X,W,removedIndices,N] = prepareData(obj,X,W, obsInRows,prepareForFit,functionname,varargin)

            if ~obsInRows
                [D,N] = size(X);
            else
                [N,D] = size(X);
            end

            removedIndices = false(N,1);

            if (~isfloat(X) || ~ismatrix(X)) % error if input is a table or it's not a floating-point matrix
                error(message('stats:classreg:learning:FullClassificationRegressionModel:prepareDataCR:BadXType'));
            end

            % verify that X is not sparse, complex, or a GPU array.
            internal.stats.checkSupportedNumeric('X',X);

            if prepareForFit
                if obj.NumPredictors == 0
                    obj.PredictorClass = class(X);
                    obj.NumPredictors = D;
                else
                    if obj.PrivPredictorClassFixed
                        validateattributes(X,{obj.PredictorClass},{},functionname,'X');
                    else % do not set obj.PrivPredictorClassFixed here
                        obj.PredictorClass = class(X);
                    end
                    classreg.learning.internal.numPredictorsCheck(X,...
                        obj.NumPredictors,obsInRows);
                end
                if obj.NumComponents == 0
                    obj.NumComponents = obj.NumPredictors;
                else
                    if obj.NumPredictors < obj.NumComponents
                        error(message('stats:incremental:PCA:InvalidNumPredictors'));
                    end
                end
            else
                classreg.learning.internal.numPredictorsCheck(X,...
                    obj.NumPredictors,obsInRows);
            end

            % validate Weights
            W = obj.validateWeights(X,W,N,functionname);

            % remove any observations that have nan values or corresponding
            % zero or nan weight values
            anynanrows = any(isnan(X),2) | any(isnan(W)) | any(W==0);
            if any(anynanrows)
                removedIndices = anynanrows;
                if obsInRows
                    X(anynanrows,:) = [];
                    N = size(X,1);
                else
                    X(:,anynanrows) = [];
                    N = size(X,2);
                end
                W(anynanrows) = [];
            end
        end

        function statesdefined = statesDefined(obj)
            statesdefined = statesDefined(obj.Impl);
        end

        function obj = updateStateParameters(obj,numPredictors,predictorClass)
            obj.Impl = updateStateParameters(obj.Impl,numPredictors,predictorClass);
        end

        function [numPredictors,numComponents,predictorClass] = predictorAttributes(obj)
            [numPredictors,numComponents,predictorClass] = predictorAttributes(obj.Impl);
        end

        function obj = updatePredictorInformation(obj)
            [~,~,predictorClass] = predictorAttributes(obj);
            if ~isequal(obj.PredictorClass,predictorClass)
                obj.PredictorClass = predictorClass;
            end
        end

        function obj = updateWarmupPeriod(obj)
            if 0 < obj.NumComponents && obj.WarmupPeriod < obj.NumComponents
                obj.WarmupPeriod = obj.NumComponents;
            end
        end

        function cantrain = canTrain(obj)
            cantrain = canTrain(obj.Impl);
        end

        function obj = updateIsWarmIfNeeded(obj,numObservations)
            if ~obj.IsWarm
                if canTransform(obj) && canTrain(obj) && obj.WarmupPeriod <= numObservations
                    obj.IsWarm = true;
                end
            end
        end

        function obj = synchronizePredictorInformation(obj,X,obsInRows)
            % If model can not predict, information from incoming data
            % needs to extracted. After extraction, state information is
            % updated. Predictor information between derived model and
            % anomaly detector is synchronized.

            if obsInRows
                D = size(X,2);
            else
                D = size(X,1);
            end

            if (~canTransform(obj) || ~obj.PrivPredictorClassFixed) && D ~=0
                if isempty(X)
                    obj = updateStateParameters(obj,obj.NumPredictors,obj.PredictorClass);
                    obj = updatePredictorInformation(obj);
                    obj = updateWarmupPeriod(obj);
                    obj = updateIsWarmIfNeeded(obj,obj.NumTrainingObservations);
                else
                    obj = updateStateParameters(obj,obj.NumPredictors,obj.PredictorClass);
                    obj = updatePredictorInformation(obj);
                    obj.PrivPredictorClassFixed = true;
                    obj = updateWarmupPeriod(obj);
                    obj = updateIsWarmIfNeeded(obj,obj.NumTrainingObservations);
                end
            end
        end

        function cantransform = canTransform(obj)
            % can predict if states are defined, NumPredictors is not 0
            if statesDefined(obj) && obj.NumPredictors ~= 0
                cantransform = true;
            else
                cantransform = false;
            end
        end

    end

    methods (Static, Access = protected)

        function inputisempty = InputIsEmpty(X,obsInRows)
            inputisempty = false;
            if isequal(X,[])
                inputisempty = true;
                return;
            end
            if obsInRows
                if size(X,2) == 0
                    inputisempty = true;
                    return;
                end
            else
                if size(X,1) == 0
                    inputisempty = true;
                    return;
                end
            end
        end

        function [W,ObservationsIn] = parseCommonInputs(NameValueArgs)
            % Parse allowed inputs.
            % ObservationsIn is set to 'rows' as 'columns' is not
            % supported.
            arguments
                NameValueArgs.Weights = []
            end
            W = NameValueArgs.Weights;
            ObservationsIn = 'rows';
        end

        function obsin = validateObservationsIn(obs,functionname)
            obsval = validatestring(obs,{'rows','columns'},functionname,'ObservationsIn');
            if obsval(1) == 'r'
                obsin = true;
            else
                obsin = false;
            end
        end


        function W = validateWeights(X,W,N,functionname)

            if isempty(W)
                W = ones(N,1,'like',X);
            end

            internal.stats.checkSupportedNumeric('W',W);

            validateattributes(W,{'single','double'},{'vector','numel',N},functionname,'Weights');

            if ~isempty(W) && (any(W<0) || all(W==0))
                error(message('stats:classreg:learning:FullClassificationRegressionModel:prepareDataCR:NegativeWeights'));
            end
            W = W(:);
        end

        function [offlineParams,parametersAreFromOffline,remainingArgs] = extractOfflineParameters(functionname,varargin)

            p = inputParser;
            p.StructExpand = true;
            p.KeepUnmatched = true;


            addParameter(p,'Coefficients',[]);
            addParameter(p,'Latent',[]);
            addParameter(p,'NumObservations',[]);
            addParameter(p,'Means',[]);

            parse(p,varargin{:});
            args = p.Results;
            remainingArgs = p.Unmatched;

            offlineParams = struct;
            % if any positional inputs exist, first two must be nonempty
            coeffIsEmpty = isempty(args.Coefficients);
            latentIsEmpty = isempty(args.Latent);
            meansIsEmpty = isempty(args.Means);
            numObsIsEmpty = isempty(args.NumObservations);
            parametersAreFromOffline = ~coeffIsEmpty || ~latentIsEmpty || ~meansIsEmpty || ~numObsIsEmpty;

            if parametersAreFromOffline
                if coeffIsEmpty || latentIsEmpty
                    error(message('stats:incremental:PCA:CoefficientsLatentMissing'));
                end

                % inputs are from pcacov output (~coeffIsEmpty && ~latentIsEmpty must be true)
                pcaFromCov = meansIsEmpty && numObsIsEmpty;                

                validateattributes(args.Coefficients,{'double','single'},{'2d','nonempty','real','finite'}, functionname,'Coefficients');



                internal.stats.checkSupportedNumeric('Coefficients',args.Coefficients);
                numPredictors = size(args.Coefficients,1);
                numComponents = size(args.Coefficients,2);
                predictorClass = class(args.Coefficients);

                coeff = args.Coefficients;

                if pcaFromCov % if pcaFromCov, coeff must be a square matrix
                    if numPredictors~=numComponents
                        error(message('stats:incremental:PCA:NonSquareCoefficients'));
                    end
                end                

                if ~pcaFromCov
                    if numObsIsEmpty
                        error(message('stats:incremental:PCA:InputMissing','NumObservations'));
                    end
                end

                % delay checking size of Latent
                if ~meansIsEmpty
                    validateattributes(args.Means,{predictorClass},{'vector','nonempty','real','finite','numel',numPredictors},functionname,'Means');
                end

                mu = args.Means(:)';
                internal.stats.checkSupportedNumeric('Means',mu);

                if ~numObsIsEmpty
                    validateattributes(args.NumObservations,{'numeric'},{'scalar','integer','positive'},functionname,'NumObservations');
                end

                args.NumObservations = double(args.NumObservations);
                
                DOF = args.NumObservations;

                internal.stats.checkSupportedNumeric('NumObservations',DOF);

                if isempty(DOF)
                    DOF = 0;
                end

                trivialMeans = meansIsEmpty || all(mu==0);
                if ~trivialMeans
                    DOF = max(0,DOF-1);
                end

                % if Means is provided, Latent from offline PCA can either
                % have NumPredictors number of elements if Economy = false in offline PCA, or
                % it can have NumPredictors-1 number of elements if Economy = true and DOF < NumPredictors.
                % Both are valid inputs, so try both.

                % When Means is not provided, check for the following cases:
                %  * Latent can have size equal to NumPredictors
                %  * Latent can have size equal to DOF
                %  * Latent can have size equal to NumComponents
                %  If Latent has size NumPredictors-1 or DOF-1, provide an
                %  error about Means missing.
                % When Means is provided, check for the following cases:
                %  * Latent can have size equal to NumPredictors
                %  * Latent can have size equal to NumPredictors-1
                %  * Latent can have size equal to DOF
                %  * Latent can have size equal to NumComponents

                validateattributes(args.Latent,{predictorClass},{'vector','real','finite','nonincreasing'},functionname,'Latent');
                if pcaFromCov
                    if numel(args.Latent)~=numPredictors
                        error(message('stats:incremental:PCA:InvalidPCACovLatentLength'));
                    end
                    latent = args.Latent(:);
                else
                    latentsize = numel(args.Latent(:));
                    if trivialMeans
                        if ~(latentsize== numPredictors || latentsize == DOF || latentsize == numComponents)
                            if (latentsize == numPredictors-1) || (latentsize == DOF-1)
                                error(message('stats:incremental:PCA:InputMissing','Means'));
                            else
                                vals = strjoin(string(unique(sort([numPredictors,numComponents,DOF]))),', ');
                                error(message('stats:incremental:PCA:InvalidLatentLengthWithoutMeans',vals));
                            end
                        end
                    else
                        if ~(latentsize== numPredictors || latentsize == DOF || latentsize == numComponents || latentsize == numPredictors-1)
                            vals = strjoin(string(unique(sort([numPredictors,numComponents,DOF,numPredictors-1]))),', ');
                            error(message('stats:incremental:PCA:InvalidLatentLengthWithMeans',vals));
                        end
                    end
                    latent = zeros(numPredictors,1,predictorClass);
                    latent(1:latentsize) = args.Latent(:);
                end


                internal.stats.checkSupportedNumeric('Latent',latent);

                offlineParams.NumPredictors = numPredictors;
                offlineParams.NumComponents = numComponents;
                offlineParams.Coefficients = coeff;
                offlineParams.Latent = latent;
                offlineParams.Mu = mu;
                offlineParams.DOF = DOF;
                offlineParams.PredictorClass = predictorClass;
                offlineParams.CenterData = false;
                offlineParams.StandardizeData = false;
                if ~isempty(offlineParams.Mu)
                    offlineParams.CenterData = true;
                end

            end
            remainingArgs = incremental.utils.convertStructToCell(remainingArgs);
        end

        function [modelParams,implParams,parametersAreFromOffline,remainingArgs] = extractPCAParameters(functionname,varargin)

            % extract NumPredictors
            [offlineParams,parametersAreFromOffline,remainingArgs] = incrementalPCA.extractOfflineParameters(functionname,varargin{:});

            modelParams = struct;
            implParams = struct;


            if parametersAreFromOffline
                dfltWarmupPeriod = 0;
                dfltVariableWeights = [];
                p = inputParser;
                p.StructExpand = true;
                p.KeepUnmatched = true;
                addParameter(p,'VariableWeights',dfltVariableWeights);
                addParameter(p,'WarmupPeriod',dfltWarmupPeriod);
                parse(p,remainingArgs{:});
                args = p.Results;
                remainingArgs = p.Unmatched;

                implParams.NumPredictors = offlineParams.NumPredictors;
                implParams.NumComponents = offlineParams.NumComponents;
                implParams.OfflineDOF = offlineParams.DOF;
                modelParams.Coefficients = offlineParams.Coefficients;
                modelParams.Latent = offlineParams.Latent;

                implParams.Mu = offlineParams.Mu;
                implParams.CenterData = offlineParams.CenterData;
                implParams.StandardizeData = offlineParams.StandardizeData;
                modelParams.PredictorClass = offlineParams.PredictorClass;
                implParams.VariableWeights = args.VariableWeights;
                modelParams.WarmupPeriod = args.WarmupPeriod;
            else
                dfltNumPredictors = 0;
                dfltNumComponents = 0;
                dfltCenterData = false;
                dfltStandardizeData = false;
                dfltWarmupPeriod = 1000;
                dfltVariableWeights = [];

                p = inputParser;
                p.StructExpand = true;
                p.KeepUnmatched = true;
                addParameter(p,'VariableWeights',dfltVariableWeights);
                addParameter(p,'WarmupPeriod',dfltWarmupPeriod);
                addParameter(p,'NumPredictors',dfltNumPredictors);
                addParameter(p,'NumComponents',dfltNumComponents);
                addParameter(p,'StandardizeData',dfltStandardizeData);
                addParameter(p,'CenterData',dfltCenterData);
                parse(p,remainingArgs{:});
                args = p.Results;
                remainingArgs = p.Unmatched;
                modelParams.WarmupPeriod = args.WarmupPeriod;
                implParams.NumPredictors = args.NumPredictors;
                implParams.NumComponents = args.NumComponents;
                modelParams.WarmupPeriod = args.WarmupPeriod;
                modelParams.PredictorClass = 'double';
                standardize = internal.stats.parseOnOff(args.StandardizeData,'StandardizeData');
                center = internal.stats.parseOnOff(args.CenterData,'CenterData');
                implParams.StandardizeData = standardize;
                implParams.CenterData = center;
                implParams.VariableWeights = args.VariableWeights;
                implParams.OfflineDOF = 0;
            end
            checkGPUArray(implParams);
            checkGPUArray(modelParams);
        end
    end

    methods (Access=protected) % custom display methods

        function propList = propsForDisplay(obj)
            % common properties
            propList = struct;
            propList.IsWarm = obj.IsWarm;
            propList.NumComponents = obj.NumComponents;
            propList.Coefficients = obj.Coefficients;
            propList.Latent = obj.Latent;
        end

        function propGrp = getPropertyGroups(obj)
            % override used by the CustomDisplay interface to display properties
            propList =  propsForDisplay(obj);
            propGrp =  matlab.mixin.util.PropertyGroup(propList);
        end

        function header = getHeader(obj)
            header = [];
            internal.stats.displayClassName(obj);
        end

        function footerText = getFooter(obj)
            % override used by the CustomDisplay interface to display the footer
            % use the same footer defined in
            % internal.stats.displayMethodsProperties for classreg objects
            footerText = [];
            internal.stats.displayMethodsProperties(obj);
        end
    end

    methods(Hidden)
        function rep = compactRepresentationForSingleLine(obj, varargin)
            % Hidden method to dispatch to compactRepresentationForSingleLine
            %   for matlab.mixin.CustomCompactDisplayProvider.
            rep = compactRepresentationForSingleLine@matlab.mixin.CustomCompactDisplayProvider(obj, varargin{:});
        end

        function rep = compactRepresentationForColumn(obj, varargin)
            % Hidden method to dispatch to compactRepresentationForColumn
            %   for matlab.mixin.CustomCompactDisplayProvider.
            rep = compactRepresentationForColumn@matlab.mixin.CustomCompactDisplayProvider(obj, varargin{:});
        end
    end
end

% LocalWords:  orthonormality nonincreasing


function checkGPUArray(input,name)

if isstruct(input)
    fnames = fieldnames(input);
    for ii = 1:numel(fnames)
        fname = fnames{ii};
        if isa(input.(fname), 'gpuArray')
            error(message('stats:internal:utils:NoGPUArrayNamed',fname));
        end
    end
else
    narginchk(2,2);
    assert(~iscell(input));
    if isa(input, 'gpuArray')
        error(message('stats:internal:utils:NoGPUArrayNamed',name));
    end
end

end