classdef LocalOutlierFactor < anomaly.AnomalyDetector
    %LocalOutlierFactor local outlier factor.
    %   LocalOutlierFactor is a model for anomaly detection.
    %
    %   An object of this class cannot be created by calling the constructor.
    %   Use lof to create a LocalOutlierFactor object by fitting the model
    %   to the data.
    %
    %   LocalOutlierFactor properties:
    %       X                     - Matrix or table of predictors used to train this model.
    %       BucketSize            - Bucket size of each leaf node in the Kd-tree.
    %       CategoricalPredictors - Indices of categorical predictors.
    %       ContaminationFraction - Fraction of anomalies in the training data.
    %       Distance              - Distance metric.
    %       DistanceParameter     - Additional distance parameters.
    %       IncludeTies           - Flag to include the tie neighbors.
    %       SearchMethod          - Method for k-nearest neighbors search.
    %       NumNeighbors          - Number of nearest neighbors.
    %       PredictorNames        - Names of predictors used for this model.
    %       ScoreThreshold        - Threshold on anomaly score used to identify anomalies.
    %
    %   LocalOutlierFactor methods:
    %       isanomaly             - Identify anomalies in data.

    %   Copyright 2021-2023 The MathWorks, Inc.

    properties(GetAccess=public,SetAccess=private)
        %Distance Distance metric
        %   A character vector specifying the built-in distance metric.
        Distance

        %IncludeTies Whether to include the tie neighbors
        %   A logical value specifying whether to include all neighbors whose
        %   distance equal the kth smallest distance.
        IncludeTies

        %SearchMethod Method for k-nearest neighbors search
        %    'kdtree' for the Kd-tree algorithm or 'exhaustive' for the
        %    exhaustive search algorithm.
        SearchMethod

        %NumNeighbors Number of nearest neighbors.
        %   A positive integer indicating the number of nearest neighbors for
        %   computing local outlier factor values.
        NumNeighbors
    end

    properties(GetAccess=public,SetAccess=private,Dependent=true)
        % other visible properties

        %BucketSize Bucket size of each leaf node in the Kd-tree.
        %   The BucketSize property is a positive integer, indicating the
        %   maximum number of data points in each leaf node of the Kd-tree.
        BucketSize

        %DistanceParameter Additional distance parameter.
        %   Value of Distance property Value
        %
        %   'minkowski'                  A positive scalar indicating
        %                                the exponent of the minkowski
        %                                distance.
        %   'mahalanobis'                A positive definite matrix
        %                                representing the covariance
        %                                matrix used for computing the
        %                                mahalanobis distance.
        %   otherwise                    Empty.
        DistanceParameter

        %X Data used to train the model.
        %   The X property contains the predictor values. It is a table for a model
        %   trained on a table, or a numeric matrix for a model trained on a
        %   matrix. It has size N-by-P, where N is the number of rows and P is the
        %   number of predictor variables or columns in the training data.
        X
    end

    properties(GetAccess=private,SetAccess=private)
        %KDistance kth Distance
        % N-by-1 numerical array for N observations in the training data.
        % This array stores the kth-Distance values for each observation in
        % the training data.
        KDistance

        %LocalReachabilityDensity Local Reachability Density
        % N-by-1 numerical array for N observations in the training data.
        % This array stores the local reachability density values for each
        % observation in the training data.
        LocalReachabilityDensity

        %NSObject Nearest Search Object
        % Either KDTreeSearcher or ExhaustiveSearcher object used to compute
        % the k-nearest neighbors.
        NSObject

        %UniqueToOriginal Unique to Original conversion of training data
        % N-by-1 numerical array for N observations in the training data.
        % This array stores group values for duplicate observations. It is
        % used to obtain the original dataset from a dataset without
        % duplicate observations.
        UniqueToOriginal

        %Weight Weight of observations in training data
        % N-by-1 numerical array for N observations in the training data.
        % This array stores the Weight (Number of Duplicates) for each
        % unique observation in the training data.
        Weight
    end

    methods
        function [isAnomaly,scores] = isanomaly(obj,X,varargin)
            %isanomaly Identify anomalies in data.
            %   ISANOM = isanomaly(LOFOBJ,TBL) returns logical array ISANOM with N
            %   elements for N observations in table TBL. TBL contains the predictor
            %   variables. LOFOBJ is an object of class LocalOutlierFactor trained on a
            %   table of predictors.
            %
            %   ISANOM = isanomaly(LOFOBJ,X) accepts X as a matrix of predictors with
            %   one row per observation and one column per predictor for model
            %   LOFOBJ trained on a predictor matrix. You must use this syntax if you
            %   passed a matrix to lof function.
            %
            %   [ISANOM,SCORES] = isanomaly(...) also returns numeric matrix SCORES of
            %   anomaly scores with N elements using any of the input arguments in the
            %   previous syntaxes. Value close to 1 indicates a normal observation,
            %   and value greater than 1 indicates an anomaly.
            %
            %   ISANOM = isanomaly(...,Name=Value) specifies optional parameter
            %   name-value arguments:
            %       ScoreThreshold      - Numeric scalar. Observations with
            %                             scores above this threshold are
            %                             identified as anomalies.
            %                             Default: value of the
            %                             ScoreThreshold property
            %       CacheSize           - A positive scalar or 'maximal'. The default is 
            %                             1e3. Applies only when 'Distance' is 
            %                             'fasteuclidean'. If numeric, this argument 
            %                             specifies the cache size in megabytes (MB) to 
            %                             allocate for an intermediate matrix. If 'maximal', 
            %                             isanomaly attempts to allocate enough memory for an 
            %                             entire intermediate matrix whose size is MT-by-MX 
            %                             (MT is the number of rows of the training data in obj,
            %                             and MX is the number of rows of the input data X). 
            %                             If the specified cache size exceeds the available 
            %                             memory, MATLAB issues an out-of-memory error. 
            %                             'CacheSize' does not have to be large enough for 
            %                             an entire intermediate matrix, but must be at 
            %                             least large enough to hold an MX-by-1 vector. 
            %                             Otherwise, isanomaly uses 'euclidean' distance.
            %
            %   See also LocalOutlierFactor, lof, ScoreThreshold.

            [isAnomaly,scores] = obj.isanomalyHelper(X,varargin{:});
        end
    end

    methods(Access=protected) % Constructor and score methods
        function obj = LocalOutlierFactor(dataSummary,modelParams,...
                nsobject,W,unq_to_orig)

            obj = obj@anomaly.AnomalyDetector(dataSummary.CategoricalPredictors,...
                dataSummary.PredictorNames,dataSummary.WasInputTable,...
                dataSummary.VariableRange,modelParams.ContaminationFraction);

            obj.NSObject                  = nsobject;
            obj.Weight                    = W;
            obj.UniqueToOriginal          = unq_to_orig;
            obj.Distance                  = modelParams.Distance;
            obj.SearchMethod              = modelParams.SearchMethod;
            obj.IncludeTies               = modelParams.IncludeTies;
            obj.NumNeighbors              = modelParams.NumNeighbors;
        end

        function [s,obj] = score(obj,X,varargin)
            %score Compute anomaly scores in data.
            %   s = score(LOFOBJ,X) returns anomaly scores S with N elements
            %   for N observations in data X.
            %
            %   [s,obj] = score(...) also returns the updated model OBJ
            %   with values for private properties KDistance and
            %   LocalReachabilityDensity when X is training data passed
            %   while using the lof function. The model OBJ is not updated
            %   when X is test data passed using the isanomaly function.
            %
            %   [s,obj] = score(...) specifies optional parameter name-value
            %   arguments:
            %       isTraining      - logical flag. Value true indicates
            %                         that X is training data and false
            %                         indicates that X is test data.

            [varargin{:}] = convertStringsToChars(varargin{:});
            args = {'isTraining','cachesize'};
            defs = {false,1e+3};
            [isTraining,cachesize] = internal.stats.parseArgs(args,defs,varargin{:});

            N          = size(X,1);
            k          = obj.NumNeighbors;
            incties    = obj.IncludeTies;
            nsobj      = obj.NSObject;
            Wdup       = obj.Weight;
            distmetric = obj.Distance;
            srchmethod = obj.SearchMethod;

            % While training, knnsearch includes the queried observation as
            % the first neighbor - so we need to find an additional neighbor
            % and remove the index and distance corresponding to the queried
            % observation under certain conditions.
            if isTraining
                kval = k + 1;
            else
                kval = k;
                if ismember(obj.Distance, {'fasteuclidean'})
        	        statslib.internal.checkArgsForFastEuc(obj.Distance,issparse(X),cachesize,size(obj.X,1),class(obj.X));
                end
            end

            % Step1: Find the k distinct neighbors
            if strcmpi(srchmethod, 'exhaustive') 
                [idxNN,distNN] = knnsearch(nsobj,X,'k',kval,'includeties',incties,'CacheSize',cachesize);
                
            else
                [idxNN,distNN] = knnsearch(nsobj,X,'k',kval,'includeties',incties);
            end

            if isTraining
                % Choose the correct k Neighbors from the k+1 Neighbors found
                [idxNN,distNN] = LocalOutlierFactor.findKNeighbors(idxNN,distNN,Wdup,incties,k);
            end

            % Remove/ignore neighbors with NaN distance.
            if ~incties
                % Wdist is used to identify neighbors with nan
                % distance to the observations.
                Wdist = true(size(distNN));
                Wdist(isnan(distNN)) = false;
            else
                for i = 1:N
                    numNeighbors = sum(~isnan(distNN{i,:}));
                    idxNN{i}  = idxNN{i}(1:numNeighbors);
                    distNN{i} = distNN{i}(1:numNeighbors);
                end
            end

            % Edge Case:
            % Throw error for observations with zero distance from all neighbors
            % This may occur when using certain distance metrics with specific data
            % properties.
            if isTraining
                if ~incties
                    % Don't error when observations has nan distance from all
                    % neighbors.
                    if any(sum(distNN,2,'omitnan') == 0 & any(Wdist,2))
                        error(message('stats:anomaly:lof:ZeroDistanceFromAllNeighbors',distmetric))
                    end
                else
                    for ij = 1:N
                        % Don't error when observations has nan distance from all
                        % neighbors.
                        if ~isempty(distNN{ij,:}) && sum(distNN{ij,:}) == 0
                            error(message('stats:anomaly:lof:ZeroDistanceFromAllNeighbors',distmetric))
                        end
                    end
                end
            end

            if isTraining
                % Step2: Find the kth-distance of all observations in the
                % training data. kth-distance is the maximum distance
                % between an observation and all its neighbors.
                % NOTE: We ignore remaining nan values in distance matrix
                % when computing KDistance for training observations.
                % k-Distance is nan when distance to all neighbors is nan.
                if ~incties
                    kdist  = max(distNN,[],2);
                else
                    % class type should be same as distance type
                    kdist  = zeros(N,1,'like',distNN{1});
                    for idx = 1:N
                        tempdist = distNN{idx,:};
                        if isempty(tempdist)
                            tempdist = nan;
                        end
                        kdist(idx)  = max(tempdist);
                    end
                end
                obj.KDistance = kdist;
            else
                % When passing a test observation, we need the kth-distance
                % values of its neighbors (points in the training data).
                kdist = obj.KDistance;
            end

            % While testing, if N=1, index, distance and wdist for that
            % observation need to be column vectors to avoid transposing
            % other variables.
            if ~isTraining && ~incties && N==1
                idxNN  = idxNN(:);
                distNN = distNN(:);
                Wdist  = Wdist(:);
                dim    = 1;
            else
                dim    = 2;
            end

            % Step 3: Compute Reachability distance and Local reachability
            % density - reachdist(obsA,obsB) = max(d(obsA,obsB),kdist(obsB))
            if ~incties
                % NOTE: We DO NOT ignore NaN values in distance matrix when
                % computing reachability distance. If we omitnan instead,
                % we incorrectly assume that the queried observation lies
                % within the neighborhood of neighboring training data.
                reachdist = max(distNN,kdist(idxNN),'includenan');
                Wnew      = Wdup(idxNN);

                % While training, if the queried observation has a
                % duplicate point, then that observation has been included
                % in the neighborhood. We decrement the weight by 1
                % to not include the weight for the queried observation.
                if isTraining
                    colidx = idxNN == (1:N)';
                    Wnew(colidx) = Wnew(colidx) -1;
                end

                % Ignore neighbors with nan distance when computing LRD and
                % LOF values.
                Wnew = Wnew.*Wdist;

                % LRD is nan when distance to all neighbors is nan.
                lrd = 1 ./(sum(reachdist.*Wnew,dim,'omitnan')./sum(Wnew,dim));
            else
                lrd = zeros(N,1,'like',distNN{1});
                for idy =1:N
                    idxnbrs  = idxNN{idy};
                    dist     = distNN{idy,:};
                    dist     = dist(:);
                    w        = Wdup(idxnbrs);
                    % While training, decrement the weights by 1 for
                    % observations with duplicates
                    if isTraining
                        ij = idxnbrs == idy;
                        w(ij) = w(ij) - 1;
                    end
                    reachdist = max(dist,kdist(idxnbrs),'includenan');
                    % LRD is nan when distance to all neighbors is nan.
                    lrd(idy)  = 1 ./ (sum(reachdist.*w,'omitnan')/sum(w));
                end
            end

            if isTraining
                obj.LocalReachabilityDensity = lrd;
                lrdneigh = lrd;
            else
                % When passing a test observation, we need the LRD values
                % of the test observation (lrd) and its neighbors (lrdneigh)
                lrdneigh = obj.LocalReachabilityDensity;
            end

            % Step 4: Compute Local Outlier factor
            if ~incties
                lrdratio = lrdneigh(idxNN)./lrd;
                % LOF is nan when distance to all neighbors is nan.
                lof = sum(lrdratio,dim,'omitnan')./sum(Wnew,dim);
            else
                lof = zeros(N,1,'like',distNN{1});
                for idz = 1:N
                    idxnbrs  = idxNN{idz};
                    w        = Wdup(idxnbrs);

                    % While training, decrement the weights by 1 for
                    % observations with duplicates
                    if isTraining
                        ij = idxnbrs == idz;
                        w(ij) = w(ij) - 1;
                    end

                    lrdratio = lrdneigh(idxnbrs)./lrd(idz);
                    % LOF is nan when distance to all neighbors is nan.
                    lof(idz) = sum(lrdratio,'omitnan')/sum(w);
                end
            end
            s = lof;
        end
    end

    methods (Access=protected) % Display and validation methods
        function propList = propsForDisplay(obj)
            propList = propsForDisplay@anomaly.AnomalyDetector(obj); % get common properties

            % obj.X is not displayed
            propList.('NumNeighbors')          = obj.NumNeighbors;
            propList.('Distance')              = obj.Distance;

            if ~isempty(obj.DistanceParameter)
                propList.('DistanceParameter') = obj.DistanceParameter;
            end

            propList.('IncludeTies')           = obj.IncludeTies;
            propList.('SearchMethod')          = obj.SearchMethod;

            if ~isempty(obj.BucketSize)
                propList.('BucketSize')        = obj.BucketSize;
            end
        end

        function X = validateXTest(obj,X)
            % Allow sparse data
            X = validateXTest@anomaly.AnomalyDetector(obj,X,'okSparse',true);
            % Map unknown categories in test data to 0 when using jaccard
            % or hamming distance
            if any(obj.CategoricalPredictors) && ...
                    any(strcmp(obj.Distance,{'hamming' 'jaccard'}))
                % NOTE: Filter correctly when mixed data is supported
                X(isnan(X)) = 0;
            end
        end
    end

    methods (Static=true,Hidden=true)
        function [mdl,isAnomaly,scores] = fitandscore(X, varargin)

            [Xuniq,W,unqtoorig,badrows,dataSummary,extraArgs] = ...
                LocalOutlierFactor.prepareData(X, varargin{:});

            % Validate lof specific parameters and hold them in a structure
            modelParams = LocalOutlierFactor.validateModelParameters(Xuniq,dataSummary,...
                extraArgs{:});

            numGoodRows = sum(~badrows);
            if size(Xuniq,1) < numGoodRows % Duplicates exist
                % If duplicates are present, number of observations should
                % be greater than or equal to 3 because NumNeighbors
                % need to be greater than 1.
                if size(Xuniq,1) < 3
                    error(message('stats:anomaly:lof:NumRowsWithDuplicates'));
                end

                % NumNeighbors should be greater than 1 because the duplicate
                % observation is treated as the first neighbor and has "zero"
                % reachability distance.
                if modelParams.NumNeighbors < 2
                    error(message('stats:anomaly:lof:NumNeighborsWithDuplicates'));
                end

                modelParams.NumNeighbors = ...
                    max(2,min(modelParams.NumNeighbors,size(Xuniq,1)-1));
            else

                unqtoorig = [];

                % If no duplicates are present, number of observations should
                % be greater than or equal to 2. Xmat is same as Xuniq here.
                if size(Xuniq,1) < 2
                    error(message('stats:anomaly:lof:NumRowsWithoutDuplicates'));
                end

                modelParams.NumNeighbors =...
                    max(1,min(modelParams.NumNeighbors,size(Xuniq,1)-1));
            end

            % Create Search Object
            if strcmpi(modelParams.SearchMethod,'kdtree')
                nsobject = KDTreeSearcher(Xuniq,'distance',modelParams.Distance,...
                    'p',modelParams.Exponent,'bucketSize',modelParams.BucketSize);
            else
                nsobject = ExhaustiveSearcher(Xuniq,'distance',modelParams.Distance,...
                    'p',modelParams.Exponent,'cov',modelParams.Cov);
            end

            % 'BucketSize','Cov','Exponent' are captured in nsobject.
            %  Distance is passed via modelParams instead of nsobject
            %  as this would be needed when supporting heteropdist2 distance
            %  measures.
            modelParams = rmfield(modelParams,{'BucketSize','Cov','Exponent'});

            % Model
            mdl = LocalOutlierFactor(dataSummary,modelParams,...
                nsobject,W,unqtoorig);

            % When computing score on training data we store the K-distance values and
            % LRD values in the object as private properties.
            [scores,mdl] = score(mdl,Xuniq,'isTraining',true,'CacheSize',modelParams.CacheSize);

            % Fill scores for duplicates.
            % Scores for duplicate observations is the same as those for
            % unique observations.
            if ~isempty(unqtoorig)
                scores = scores(unqtoorig);
            end
            % Fill scores for observations with missing in row
            % Scores are filled as NaN values
            scores = internal.stats.insertnan(badrows,scores);

            % Computation of the score threshold and isAnomaly
            [mdl,isAnomaly] = anomaliesInXTrain(mdl,scores);
        end
    end

    methods (Access=protected,Static=true)
        function validateScoreThreshold(thresh)
            validateattributes(thresh,{'single','double'},{'scalar','real','>=',0},'isanomaly','ScoreThreshold')
        end

        function [Xuniq,W,unqtoorig,badrows,dataSummary,extraArgs] = prepareData(X, varargin)

            % P is the minkowski exponent, parsed here to avoid conflict
            % with PredictorNames
            args = {{'p','exponent'},'predictornames'};
            defs = {              []               []};
            [pmink,prednames,~,extraArgs1] = internal.stats.parseArgs(args,defs,varargin{:});

            % Call prepareData() from base class
            [Xmat, dataSummary, extraArgs2] = anomaly.AnomalyDetector.prepareData(X, extraArgs1{:},...
                'predictornames',prednames);

            args = {'ordinalpredictors'};
            defs = {                [] };
            [ordpreds,~,extraArgs3] = internal.stats.parseArgs(args,defs,extraArgs2{:});

            % Extra arguments needed as output
            extraArgs = [extraArgs3,{'exponent',pmink}];

            % Error out if Ordinal Predictors are present since distance
            % cannot be calculated for ordinals.
            if any(ordpreds)
                error(message('stats:anomaly:lof:OrdinalsNotSupported'));
            end

            % pdist2 and knnsearch can handle sparse data
            internal.stats.checkSupportedNumeric('X',Xmat,false,true);

            D = size(Xmat,2);
            if ~ (isempty(dataSummary.CategoricalPredictors) ||...
                    (length(dataSummary.CategoricalPredictors) == D &&...
                    all(dataSummary.CategoricalPredictors == (1:D))))

                if dataSummary.WasInputTable
                    error(message('stats:ClassificationKNN:ClassificationKNN:BadCategoricalTable'));
                else
                    error(message('stats:ClassificationKNN:ClassificationKNN:BadCategoricalPre'));
                end
            end

            % Remove rows containing NaN values
            [Xmat,badrows] = LocalOutlierFactor.removeBadRows(Xmat);

            % Remove rows containing duplicates
            [Xuniq,W,unqtoorig] = LocalOutlierFactor.removeDuplicates(Xmat);
        end

        function modelParams = validateModelParameters(X,dataSummary,varargin)

            filename = 'lof';
            [modelParams,varargin] = anomaly.AnomalyDetector.validateModelParameters(filename,varargin{:});

            % 1.1 Set parameter defaults.
            dfltNumberNeighbors        = 20;
            dfltDistance               = [];
            dfltNSMethod               = [];
            dfltIncludeTies            = false;
            dfltExponent               = [];
            dfltCov                    = [];
            dfltBucketSize             = 50;
            dfltCacheSize              = 1e+3;

            % 1.2 Parse optional name/value pairs.
            paramNames = {'numneighbors', 'distance', {'nsmethod','searchmethod'},...
                'includeties', {'p','exponent'}, 'cov', 'bucketsize', 'cachesize'};
            paramDflts = {dfltNumberNeighbors, dfltDistance, dfltNSMethod,...
                dfltIncludeTies, dfltExponent, dfltCov, dfltBucketSize, dfltCacheSize};

            [numNN,distype,nsmethod,incties,p,cov,bsize,cachesize] = ...
                internal.stats.parseArgs(paramNames,paramDflts,varargin{:});

            % 1.3 Validate optional name/value pairs.

            % NumNeighbors
            if isempty(numNN)
                numNN = dfltNumberNeighbors;
            else
                validateattributes(numNN,{'single','double'},{'scalar','real','positive','integer'},...
                    filename,'NumNeighbors')
            end

            % IncludeTies
            if isempty(incties)
                incties = dfltIncludeTies;
            else
                incties = internal.stats.parseOnOff(incties,'IncludeTies');
            end

            % BucketSize
            if isempty(bsize)
                bsize = dfltBucketSize;
            else
                validateattributes(bsize,{'single','double'},{'scalar','real','positive','integer'},...
                    filename,'BucketSize')
            end

            % Distance
            iscatpresent = any(dataSummary.CategoricalPredictors);
            if isempty(distype)
                if ~iscatpresent
                    distype ='euclidean';
                else
                    distype ='hamming';
                end
            elseif ischar(distype)
                distype = validatestring(distype,{'euclidean','fasteuclidean',...
                    'mahalanobis','cityblock','minkowski',...
                    'chebychev','cosine','correlation','spearman',...
                    'hamming','jaccard'},filename,'Distance');
            else
                error(message('stats:anomaly:lof:BadDistanceType'));
            end

            % CacheSize
            if ismember(distype, {'fasteuclidean'})
                distype = statslib.internal.checkArgsForFastEuc(distype,issparse(X),cachesize,size(X,1),class(X));
            end

            % Exponent and Cov - complete validation of value is done by knnsearch
            if ~isempty(p)
                if ~strcmpi(distype,'minkowski')
                    error(message('stats:ExhaustiveSearcher:knnsearch:InvalidMinExp'));
                end
            end

            if ~isempty(cov)
                if ~strcmpi(distype,'mahalanobis')
                    error(message('stats:ExhaustiveSearcher:knnsearch:InvalidMahaCov'));
                end
            end

            % SearchMethod
            D = size(X,2);
            if isempty(nsmethod)
                distList = {'euclidean'; 'fasteuclidean'; 'cityblock'; 'chebychev'; 'minkowski'};
                isminkdist = any(strcmpi(distype,distList));
                % If the distance belongs to the Minkowski distance family and the
                % dimension is less than 10 a KDTreeSearcher is used,
                % otherwise a ExhaustiveSearcher is used
                if isminkdist && D <= 10  && ~issparse(X)
                    % KDTreeSearcher with regular Euclidean distance metric is faster
                    % than ExhaustiverSearcher with fast Euclidean distance in this
                    % case
                    if strcmpi(distype, 'fasteuclidean')
                        warning(message('stats:anomaly:lof:IgnoreFastEuc',distype,'euclidean'));
                        distype = 'euclidean';
                    end
                    nsmethod = 'kdtree';
                else
                    nsmethod = 'exhaustive';
                end
            elseif ischar(nsmethod)
                nsmethod = validatestring(nsmethod,{'kdtree','exhaustive'},...
                    filename,'SearchMethod');
            else
                error(message('stats:anomaly:lof:BadSearchMethodType'));
            end

            if strcmpi(nsmethod, 'kdtree')  && ~ismember(distype, {'euclidean', 'cityblock', 'chebychev', 'minkowski'}) 
                error(message('stats:anomaly:lof:InappropriateDistForKdtree', distype));
            end

            % create modelParams structure
            modelParams.BucketSize               = bsize;
            modelParams.Cov                      = cov;
            modelParams.Distance                 = distype;
            modelParams.Exponent                 = p;
            modelParams.IncludeTies              = incties;
            modelParams.SearchMethod             = nsmethod;
            modelParams.NumNeighbors             = numNN;
            modelParams.CacheSize                = cachesize;
        end
    end

    methods (Static=true,Hidden=true) % Helper methods
        function [xUnique,numOccurrences,uniqueToOriginal] = removeDuplicates(X)
            %removeDuplicates - Remove duplicate rows from X
            %   [XUnique,NumOccurances,UniqueToOriginal] = removeDuplicates(X)
            %   takes a N-by-P matrix X and removes duplicate rows from X.
            %   Xunique is a Nunique-by-1 matrix with unique rows of X.
            %
            %   NumOccurances is a Nunique-by-1 vector indicating the number of
            %   occurrences of each row of Xunique in X.
            %
            %   UniqueToOriginal is a N-by-1 logical vector marking the group
            %   of duplicate rows.

            N = size(X,1);
            % Remove duplicates and capture number of occurrences of unique
            % observations
            numOccurrences = ones(N,1,'like',X);
            [xUnique,~,uniqueToOriginal] = unique(X,'rows','stable');
            if size(xUnique,1) < size(X,1) % Duplicates exist
                numOccurrences = accumarray(uniqueToOriginal,numOccurrences);
            end
        end

        function [idxNN,distNN] = findKNeighbors(idxNN,distNN,Wdup,incties,k)
            N = size(Wdup,1);
            if ~incties
                removeidx = (1:N)';
                % Keep the queried observation as a neighbor if its Weight
                % is not 1 i.e if the observation has duplicates.
                keepmask   = ~((idxNN == removeidx) & (Wdup == 1));

                % If you keep the queried observation, remove the k+1th
                % neighbor.
                keepdup               = all(keepmask,2);
                keepmask(keepdup,end) = false;

                % Transpose index and distance for correct column-wise
                % reshape. This is significantly faster than using a for
                % loop when N > 10,000 (keeping K = 5000)
                idxNNtp  = idxNN';
                distNNtp = distNN';
                masktp   = keepmask';
                idxNN    = reshape(idxNNtp(masktp),k,N)';
                distNN   = reshape(distNNtp(masktp),k,N)';
            else
                % Using a for loop here is faster than using cellfun
                kval = k+1;
                for i = 1:N
                    % ith observation is the queried observation
                    % Find the k+1th neighbor and all its ties
                    numidx = length(idxNN{i,:});
                    tieidx = 1:numidx >= kval;

                    keepidx   = true(1,numidx);
                    % Don't keep the queried observation as a neighbor if its weight is 1
                    queryidx  = idxNN{i,:} == i;
                    if (Wdup(i) == 1)
                        keepidx(queryidx) = false;
                    end

                    % If queried observation is kept as a neighbor:
                    % Remove the k+1th neighbor and its ties if kth neighbor's distance
                    % is not equal to k+1th neighbor's distance. If kth neighbor's distance and k+1th
                    % neighbor's distance are equal, k+1th neighbor and its ties will be included as
                    % ties of the kth neighbor.
                    keptquery = all(keepidx);
                    if keptquery && ~isequal(distNN{i}(kval),distNN{i}(kval-1))
                        keepidx(tieidx) = false;
                    end

                    idxNN{i}  = idxNN{i}(keepidx);
                    distNN{i} = distNN{i}(keepidx);
                end
            end
        end

        function  [X,badrows] = removeBadRows(X)
            %removeBadRows - Remove unusable rows from X
            %   [X,badrows] = removeBadRows(X) takes a N-by-P matrix X and
            %   and removes some rows from X. Row k is removed from X if:
            %
            %   * X(k,:) has a NaN value.
            %
            %   badrows is a N-by-1 logical vector marking the location of bad
            %   rows.

            badrows = any(isnan(X),2);
            if any(badrows)
                warning(message('stats:anomaly:lof:IgnoreMissingRows'));
                X(badrows,:) = [];
            end
        end
    end

    methods % Define getters for dependent properties
        function bsize = get.BucketSize(obj)
            if isequal(obj.SearchMethod,'kdtree')
                bsize = obj.NSObject.BucketSize;
            else
                bsize = [];
            end
        end

        function distparam = get.DistanceParameter(obj)
            distparam = obj.NSObject.DistParameter;
        end

        function xout = get.X(obj)
            xout = obj.NSObject.X;
            % Decode categorical variables into original values
            if obj.WasInputTable
                xtbl = array2table(xout,'VariableNames',obj.PredictorNames);
                for j = 1:size(xout,2)
                    vrange = obj.VariableRange{j};
                    newx = anomaly.lof.decodeCategorical(xout(:,j),vrange);
                    xtbl.(obj.PredictorNames{j}) = newx;
                end
                xout = xtbl;
            elseif ~isempty(obj.VariableRange)
                for j = 1:size(xout,2)
                    vrange = obj.VariableRange{j};
                    newx = anomaly.lof.decodeCategorical(xout(:,j),vrange);
                    xout(:,j) = newx;
                end
            end
            if ~isempty(obj.UniqueToOriginal)
                xout = xout(obj.UniqueToOriginal,:);
            end
        end
    end
end

% LocalWords:  KDistance ISANOM LOFOBJ includeties Wdist wdist reachdist kdist includenan LRD lrd
% LocalWords:  lrdneigh Xmat Xuniq nsobject heteropdist ordinalpredictors numneighbors nsmethod
% LocalWords:  searchmethod bucketsize Maha XUnique Occurances Xunique Nunique badrows
