classdef DriftDetectionMethod < incremental.drift.ConceptDriftDetector
    %DRIFTDETECTIONMETHOD Incremental drift detector that utilizes
    % Drift Detection Method (DDM).
    %
    %   DriftDetectionMethod properties:
    %       WarningThreshold         - Number of standard deviations to determine if drift might exist.
    %       DriftThreshold           - Number of standard deviations to determine if drift exists.
    %       Mean                     - Weighted average of all input data used for training.
    %       OptimalMean              - Optimal weighted average observed up to the most recent data point.
    %       OptimalStandardDeviation - Optimal weighted standard deviation observed up to the most recent data point.
    %       StandardDeviation        - Weighted standard deviation of all input data used for training.
    %       Alternative              - Type of hypothesis test.
    %       InputType                - Type of the input for which drift is tracked.
    %       NumTrainingObservations  - Nonnegative integer specifying the number of observations 
    %                                  used for training drift detector.
    %       IsWarm                   - Flag indicating if drift detector has gone through the warmup period.
    %       WarmupPeriod             - Nonnegative integer specifying the number of observations used for
    %                                  drift detector warmup.
    %       PreviousDriftStatus      - Status of drift detector prior to training with the most 
    %                                  recent data.
    %       DriftStatus              - Current status of drift detector.
    %       DriftDetected            - Flag indicating whether DriftStatus is
    %                                  'Drift'.
    %       WarningDetected          - Flag indicating whether DriftStatus is
    %                                  'Warning'.
    %
    %   DriftDetectionMethod methods:
    %       detectdrift              - Update drift detector states and drift status with new data.
    %       reset                    - Reset incremental concept drift detector.
    %
    %
    %   Example 1: Generate random bit stream and monitor for any changes.
    %
    %     % Set random seed for reproducibility
    %     rng(0)
    %     numObservations = 3000;
    %     switchPeriod = 1000;
    %     detector = incrementalConceptDriftDetector('ddm');
    %     % Value of 1 indicates a failure.
    %     % For the first 1000 observations, failure rate is 0.1.
    %     % After 1000, failure rate increases to 0.6.
    %     for i = 1:numObservations
    %         if i <= switchPeriod
    %             failurerate = 0.1;
    %         else
    %             failurerate = 0.6;
    %         end
    %         X = rand() < failurerate;
    %         detector = detectdrift(detector,X);
    %         if detector.DriftDetected
    %             sprintf("Drift detected at Observation #%d. Detector reset.",i)
    %             detector = reset(detector);
    %         elseif detector.WarningDetected
    %             sprintf("Warning detected at Observation #%d",i)
    %         end
    %     end
    %
    % References:
    %
    %   [1] Gama, J., Medas, P., Castillo, G., & Rodrigues, P. (2004, September).
    %       Learning with drift detection. In Brazilian symposium on artificial
    %       intelligence (pp. 286-295). Springer, Berlin, Heidelberg.
    %
    %   See also incrementalConceptDriftDetector.

    %   Copyright 2021-2022 The MathWorks Inc.
    properties  (GetAccess = public, SetAccess = protected, Dependent)
        %WARNINGTHRESHOLD Number of standard deviations to determine if drift might exist.
        %   WarningThreshold determines how many standard deviations
        %   away the overall test statistic is allowed to be
        %   from the optimal test statistic before DriftStatus is set to
        %   'Warning'.
        %
        %   See also DriftThreshold, DriftStatus, PreviousDriftStatus.
        WarningThreshold

        %DRIFTTHRESHOLD Number of standard deviations to determine if drift exists.
        %   DriftThreshold determines how many standard deviations
        %   away the overall test statistic is allowed to be
        %   from the optimal test statistic before DriftStatus is set to
        %   'Drift'.
        %
        %   See also WarningThreshold, DriftStatus, PreviousDriftStatus.
        DriftThreshold
    end

    properties  (GetAccess = public, SetAccess = protected)
        %MEAN Weighted average of all input data used for training.
        %
        %   See also OptimalMean, StandardDeviation.
        Mean

        %OPTIMALMEAN Optimal weighted average observed up to the most recent data point.
        %       * When Alternative is 'greater', OptimalMean and OptimalStandardDeviation
        %         are updated if Mean + StandardDeviation is less than or
        %         equal to OptimalMean + OptimalStandardDeviation.
        %       * When Alternative is 'less', OptimalMean and OptimalStandardDeviation
        %         are updated if Mean - StandardDeviation is greater than
        %         or equal to OptimalMean - OptimalStandardDeviation.
        %
        %   See also Mean, OptimalStandardDeviation.
        OptimalMean = inf

        %OPTIMALSTANDARDDEVIATION Optimal weighted standard deviation
        %   observed up to the most recent data point.
        %       * When Alternative is 'greater', OptimalMean and OptimalStandardDeviation
        %         are updated if Mean + StandardDeviation is less than or
        %         equal to OptimalMean + OptimalStandardDeviation.
        %       * When Alternative is 'less', OptimalMean and OptimalStandardDeviation
        %         are updated if Mean - StandardDeviation is greater than
        %         or equal to OptimalMean - OptimalStandardDeviation.
        %
        %   See also OptimalMean, StandardDeviation.
        OptimalStandardDeviation = 0
    end

    properties (SetAccess = protected, Dependent)
        %StandardDeviation Weighted standard deviation of all input data used
        %   for training.
        %
        %   See also Mean, OptimalStandardDeviation.
        StandardDeviation
    end

    properties (Access = protected)
        DefaultWarningThreshold = 2
        DefaultDriftThreshold = 3
        Length = 0
        Var
        WeightsSum = 0
        DriftThresholdPriv
        WarningThresholdPriv
    end

    methods (Hidden)
        function this = DriftDetectionMethod(varargin)
            %DRIFTDETECTIONMETHOD Construct an incremental concept drift detector that utilizes
            % Drift Detection Method (DDM).
            %
            %   DETECTOR=incremental.drift.DriftDetectionMethod() returns an incremental concept
            %   drift detector using DDM [1] configured to detect drift for a binary input that
            %   consists 0's and 1's. By default, drift detector employs a one-sided test and tracks
            %   any changes in test statistics that imply a drift in the positive
            %   direction; i.e., given 2 test statistics F1(x) and F2(x), null
            %   hypothesis states F1(x)<=F2(x). This corresponds to
            %   'greater' for 'Alternative' Name-Value argument.
            %
            %   DETECTOR=incremental.drift.DriftDetectionMethod('PARAM1',val1,'PARAM2',val2,...) specifies optional
            %   parameter name/value pairs:
            %
            %       'Alternative'             - String indicating the type
            %                                   of hypothesis test. 'Alternative'
            %                                   refers to the alternative
            %                                   hypothesis test.
            %                                   Can be one of:
            %                                       * 'greater' -- "F1(x) > F2(x)" (one-sided test)
            %                                       * 'less' -- "F1(x) < F2(x)" (one-sided test).            
            %       'DriftThreshold'           - Number of standard deviations to determine if drift exists.
            %                                   'DriftThreshold' indicates one of the two levels used for determining drift status.
            %                                   'DriftThreshold' must be strictly greater than 'WarningThreshold'.
            %                                   'DriftThreshold' determines how many standard deviations
            %                                   away the overall test statistic is allowed to be
            %                                   from the optimal test statistic before DriftStatus is set to
            %                                   'Drift'. Default is 3 ([1]), which corresponds to 99.7 %
            %                                   confidence level.
            %       'InputType'               - Type of the input to the drift
            %                                   detector. Can be 'binary' or
            %                                   'continuous'. Default: 'binary'.
            %       'WarmupPeriod'            - Nonnegative integer specifying the number of
            %                                   observations that are used for the drift detector
            %                                   warmup. Until drift detector is warm,
            %                                   drift detector trains with incoming input and updates internal
            %                                   statistics, but does not check to see if 'DriftStatus'
            %                                   has changed. After drift detector is warm, drift detector
            %                                   starts checking for changes in 'DriftStatus'.
            %                                   Default: 30.
            %       'WarningThreshold'        - Number of standard deviations to determine if drift might exist.
            %                                   'WarningThreshold' indicates one of the two levels used for determining drift status.
            %                                   'WarningThreshold' must be strictly less than 'DriftThreshold'.
            %                                   'WarningThreshold' determines how many standard deviations
            %                                   away the overall test statistic is allowed to be
            %                                   from the optimal test statistic before DriftStatus is set to
            %                                   'Warning'. Default is 2 ([1]), which corresponds to 95 %
            %                                   confidence level.
            %
            % References:
            %
            %   [1] Gama, J., Medas, P., Castillo, G., & Rodrigues, P. (2004, September).
            %       Learning with drift detection. In Brazilian symposium on artificial
            %       intelligence (pp. 286-295). Springer, Berlin, Heidelberg.
            %
            %   See also incrementalConceptDriftDetector.

             % Pass values for supported and default values that need to be
             % modified
             newArguments = {'SupportedInputTypes',{'binary','continuous'}};
             modifiedArgsList = [varargin newArguments];
             this = this@incremental.drift.ConceptDriftDetector(modifiedArgsList{:});
        end
    end

    methods
        function std = get.StandardDeviation(this)
            std = sqrt(this.Var);
        end

        function this = set.WarningThreshold(this,warningthreshold)
            validateattributes(warningthreshold,{'double','single'},{'scalar','real','finite','positive'},mfilename,'WarningThreshold');
            this.WarningThresholdPriv = warningthreshold;
        end

        function this = set.DriftThreshold(this,driftthreshold)
            validateattributes(driftthreshold,{'double','single'},{'scalar','real','finite','positive'},mfilename,'DriftThreshold');
            if (driftthreshold <= this.WarningThresholdPriv)
                error(message('stats:incremental:drift:ConceptDriftDetector:InvalidThreshold','DriftThreshold','WarningThreshold'));
            end 
            this.DriftThresholdPriv = driftthreshold;
        end

        function warningthreshold = get.WarningThreshold(this)
              warningthreshold = this.WarningThresholdPriv;
        end

        function driftthreshold = get.DriftThreshold(this)
              driftthreshold = this.DriftThresholdPriv;
        end
    end

    methods (Access = protected)


        function this = resetDetector(this,inputClassType)
            this = initializeOptimalMean(this,inputClassType);
            this.OptimalStandardDeviation = cast(0,inputClassType);
            this.Mean = cast(0,inputClassType);
            this.Var = cast(0,inputClassType);
            this.WeightsSum = cast(0,inputClassType);
            this.Length = cast(0,inputClassType);
        end

        function this = initializeOptimalMean(this,inputClassType)
            if isequal(this.Alternative,'less')
                this.OptimalMean = cast(-inf,inputClassType);
            else
                this.OptimalMean = cast(inf,inputClassType);
            end
        end

        function [this,drift,warning] = updateDetectorStatus(this)
            drift = false;
            warning = false;
            if isequal(this.Alternative,'greater')
                term = this.Mean + this.StandardDeviation;
                updateStatisticsCondition = term <= this.OptimalMean + this.OptimalStandardDeviation;
                driftCondition = this.OptimalMean + this.DriftThreshold*this.OptimalStandardDeviation < term;
                warningCondition = this.OptimalMean + this.WarningThreshold*this.OptimalStandardDeviation < term;
            else %isequal(this.Alternative,'less')
                term = this.Mean - this.StandardDeviation;
                updateStatisticsCondition = this.OptimalMean - this.OptimalStandardDeviation <= term;
                driftCondition = term < this.OptimalMean - this.DriftThreshold*this.OptimalStandardDeviation;
                warningCondition = term < this.OptimalMean - this.WarningThreshold*this.OptimalStandardDeviation;
            end
            if updateStatisticsCondition
                this.OptimalMean = this.Mean;
                this.OptimalStandardDeviation = this.StandardDeviation;
            end
            % drift if term is nan
            if  driftCondition || isnan(term) 
                drift = true;
                warning = false;
            elseif warningCondition
                drift = false;
                warning = true;
            end
        end

        function this = updateDetectorWithOneObservation(this,X,W)

            if isempty(X)
                return;
            end

            [this.Mean,this.Var,this.WeightsSum,this.Length] = this.updateMoments(X,W,this.Mean,this.Var,this.WeightsSum,this.Length,this.Binary);

        end

        function this = initializeDetector(this,varargin)

            [~,extraArgs] = internal.stats.parseArgs({},{[]}, varargin{:});
            if ~isempty(extraArgs)
                error(message('stats:classreg:learning:FitTemplate:fillIfNeeded:UnknownBaseFitObjectParameter', extraArgs{1}));
            end
            this = initializeOptimalMean(this,this.InputClassType);
        end

    end
    methods(Access = protected,Static)
        function [newMean,newVar,newWeight,newSize] = updateMoments(X,W,oldMean,oldVar,oldWeight,oldSize,binary)
            newSize = oldSize;
            chunkMean = classreg.learning.internal.wnanmean(X,W);
            chunkWeight = sum(W);
            chunkVar = classreg.learning.internal.wnanvar(X,W,0);
            nonzeroChunkWeight = 0 < chunkWeight;
            [newMean,newVar,newWeight] = incremental.utils.combineMoments(oldMean,chunkMean,oldVar,chunkVar,oldWeight,chunkWeight);
            if binary
                newVar = oldVar;    
                if 1 < newSize  %newSize is incremented after variance calculation. Denominator is N-1; i.e., we use Bessel's correction for unbiased variance estimation
                    newVar(nonzeroChunkWeight,:) =  (newMean(nonzeroChunkWeight,:)).*((1-newMean(nonzeroChunkWeight,:)));
                    newVar(nonzeroChunkWeight,:) = newVar(nonzeroChunkWeight,:)./newSize;
                end
            end
            newSize = oldSize + numel(nonzeroChunkWeight);
        end
    end
end




