classdef SemiSupervisedGraphModel < semisupervised.Model
    
    %   SemisupervisedGraphModel properties:
    %   FittedLabels          - Labels fitted to the unlabeled data
    %   LabelScores           - Scores for the fitted labels
    %   Method                - Labeling technique used by the algorithm
    %   ClassNames            - Names of classes in Y
    %   CategoricalPredictors - Indices of categorical predictors
    %   PredictorNames        - Names of predictors used for this model
    %   ResponseName          - Name of the response variable
    
    
    %   SemisupervisedGraphModel methods:
    %   predict      - Predict labels and scores on unseen data
    
    %   Copyright 2020 The MathWorks, Inc.
    
    properties (SetAccess = protected, GetAccess = public)
        Method
    end
    
    % properties needed for prediction
    properties (SetAccess = protected, GetAccess= protected)
        PrivX
        LabeledDataScore
        Distance
        KernelScale
    end
    
    methods (Hidden)
        function mdl = SemiSupervisedGraphModel(labels, scores, method,...
                trainX, labeleddatascrore ,dataSummary, classSummary, distance, kernelscale)
            % do not protect access to the constructor, in case future
            % models want to invoke it but, try to guide the user
            if nargin~=9
                error(message('stats:SemiSupervised:DoNotUseConstructorGraph'));
            end
            
            % store data & class summaries on the object, for (1)
            % consistency in appearance with classreg objects, and (2) use
            % in prediction
            mdl = mdl@semisupervised.Model(dataSummary, classSummary);
            
            % most important user-visible properties
            mdl.FittedLabels = labels;
            mdl.LabelScores = scores;
            mdl.Method =  method;
            
            % store training data and scores for prediction
            mdl.PrivX = trainX;
            mdl.LabeledDataScore = labeleddatascrore;
            
            % store similarity parameters, needed for prediction
            mdl.Distance = distance;
            mdl.KernelScale = kernelscale;
        end
    end
    
    methods (Hidden, Static)
        function mdl = fit(X, Y, XU, varargin)
            Y = classreg.learning.internal.ClassLabel(Y); % cast all Ys to a single type
            [exponent, userProvided, varargin] = internal.stats.parseArgs({'P'},{2.0}, varargin{:}); % parse the minkowski exponent before predictor names are parsed 
            
            [X,Y,XU, dataSummary,userclassnames,nonzeroClassNames, varargin] = semisupervised.Model.prepareData(X, Y, XU, varargin{:});
            numLabeled = size(X,1);
            numClasses = numel(userclassnames); % size of scores will still be determined by user classes, similar to what we do for classreg predict()
            
            isemptyXU = isempty(XU);
                      
            trainX = [X;XU];
            
            [~,classidx] = ismember(Y, userclassnames);
            dummyY= dummyvar(classidx); % dummify as needed by the algorithm
            
            if isempty(dataSummary.CategoricalPredictors)
                distdflt = 'euclidean';
            else % for heterogeneous and categorical data, select goodall3 as the default
                distdflt = 'goodall3';
            end
            
            % parse training algorithm specific parameters
            algoparams = {'ordinalpredictors',... % ordinal predictors need to be parsed because they are populated from table conversion
                'method', ...
                'alpha',...
                'iterationlimit',...
                'tolerance',...
                'distance',...
                'kernelscale',...
                } ;
            
            algodflts = {[],...
                'labelpropagation',...
                0.01, ...
                1e3,...
                1e-3,...
                distdflt,...
                1.0,...  % kernelscale should use a default of 1 consistent with other functions in the toolbox
                }; 
            
            % parse the params
            [~, fittermethod,alpha, iterationlimit, tolerance, distance, kernelscale,...
                ~,varargin] = internal.stats.parseArgs(algoparams, algodflts, varargin{:});
            
            % validate params, for partial matches, prefer iterative methods
            validMethods = {'labelpropagation', 'labelspreading', ...
                'labelpropagationexact', 'labelspreadingexact'};
            fittermethod = validatestring(fittermethod, validMethods);
            
            fName = 'fitsemigraph';
            validateattributes(alpha,{'double', 'single'},...
                {'scalar','>',0,'<',1}, fName, 'Alpha');  % alpha  must be in the set (0,1)
            validateattributes(iterationlimit,'numeric', ...
                {'scalar', 'positive', 'integer'}, fName,'IterationLimit');
            validateattributes(tolerance, {'double', 'single'}, ...
                {'scalar', 'positive', 'real'}, fName, 'Tolerance');
            
            % allowed distances in graph construction
            if ischar(distance)
                validDists = {'euclidean'; 'seuclidean'; 'cityblock'; 'chebychev'; ...
                    'mahalanobis'; 'minkowski'; 'cosine'; 'correlation'; ...
                    'spearman'; 'hamming'; 'jaccard'; 'goodall3'; 'ofd'};
                distance = validatestring(distance,validDists,fName,'Distance');
            end
            
            % mixed data support, add a default for a mix of categorical and continuous predictors
            if strcmpi(distance, 'goodall3') || strcmpi(distance, 'ofd')  % for mixed types
                dist = @(y,X) internal.stats.heteropdist2(X,y, distance,dataSummary.CategoricalPredictors);
            else
                dist = distance; % all other distance strings, supported in knnsearch
            end
            
            % If KernelScale is set to 'auto', compute value based on heuristic
            if ischar(kernelscale)
                validatestring(kernelscale,{'auto'},fName,'KernelScale');
                if ~(strcmp(dist,'euclidean') || strcmp(dist,'seuclidean'))
                    error(message('stats:similarity:AutoInvalidDist'));
                end
                
                if strcmp(dist,'seuclidean')
                    % If the distance is seculidean standardize X before calculating
                    % optimalKernelScale
                    if  isempty(dataSummary.Mu) % if not standardized already standardize X just before calculating optimalKernelScale
                        Xstd = zscore(trainX); % treat all variables as continuous here
                    else
                        Xstd = trainX;
                    end
                    % Optimal Kernel Scale for seuclidean distance
                    kernelscale = classreg.learning.svmutils.optimalKernelScale(Xstd,[],1);
                else
                    % Optimal Kernel Scale for euclidean distance
                    kernelscale = classreg.learning.svmutils.optimalKernelScale(trainX,[],1);
                end
            end % KernelScale will be validated in the call to similarity
            
            % if P was provided by the user, append it to the list of arguments to similarity
            if userProvided.P
                varargin = [varargin, {'P', exponent}]; % if the user provides the minkowski exponent, add it to the list of arguments
            end
            
            if numClasses==1 % when the user specifies only one class, do not do any training, return outputs of appropriate sizes
                numUnlabeled = size(XU,1);
                fittedLabels = repmat(labels(userclassnames(1)),numUnlabeled,1);
                fittedLabelScores = ones(numUnlabeled,1);
                labeledDataScores = ones(numLabeled, 1);
                % make a dummy call to similarity to validate P, even for a
                % single class case, update this dummy call to account for
                % a new assert added in heteropdist2 in 21a
                x = X(1,:); % X is always non-empty
                internal.stats.similarity(x,'Distance',dist,'KernelScale', kernelscale, varargin{:});
            else % numClasses > 1 => can train
                % for label propagation do not train for empty unlabeled
                % data, memorize the labeled data and other properties
                % needed for prediction
                % for label spreading, the algorithm still needs to learn
                % soft labels based on alpha, for such cases do train
                if isemptyXU && ~strncmpi(fittermethod, 'labelspread', 11)
                    fittedLabels = repmat(labels(userclassnames(1,:)),0,1);
                    fittedLabelScores = nan(0,numClasses);
                    labeledDataScores = dummyY; % label propagation does hard label assignments for labeled data
                    % make a dummy call to similarity to validate P, even for
                    % an empty XU case/label prop
                    x = X(1,:); % X is always non-empty
                    internal.stats.similarity(x,'Distance',dist,'KernelScale', kernelscale, varargin{:});
                else
                    % compute adjacency matrix - rest of the similarity parameters
                    % are validated in similarity calculation
                    adjmat = internal.stats.similarity(trainX,'Distance',dist,'KernelScale', kernelscale, varargin{:});
                    
                    numTotal = size(adjmat,1);
                    
                    inputs = {adjmat, dummyY, numLabeled, numTotal, numClasses}; % inputs common to all methods except exact label propagation
                    fitterfunc = str2func(fittermethod);
                    
                    hasSoftLabels = false; % label spreading assigns soft labels to labeled data as well instead of a hard clamping of the user-provided labels
                    
                    switch fittermethod % add additional inputs
                        case 'labelpropagation'
                            inputs = [inputs {iterationlimit,tolerance}];
                            fittedLabelScores = fitterfunc(inputs{:});
                        case 'labelspreading'
                            inputs = [inputs {alpha, iterationlimit,tolerance}];
                            F = fitterfunc(inputs{:});
                            hasSoftLabels = true;
                            labeledDataScores = F(1:numLabeled,:);
                            fittedLabelScores = F(numLabeled+1:end,:);
                        case 'labelspreadingexact'
                            inputs = [inputs {alpha}];
                            F = fitterfunc(inputs{:});
                            hasSoftLabels = true;
                            labeledDataScores = F(1:numLabeled,:);
                            fittedLabelScores = F(numLabeled+1:end,:);
                        otherwise % 'labelpropagationexact'
                            inputs = inputs(1:4); % exact label propagation needs only the first 4 inputs - adjmat, dummyY, numLabeled, numTotal
                            fittedLabelScores = fitterfunc(inputs{:});
                    end
                    
                    normalizingSum = sum(fittedLabelScores,2);
                    if normalizingSum~=0
                        fittedLabelScores = fittedLabelScores./normalizingSum;
                    end
                    [~,maxscoreidx] = max(fittedLabelScores,[],2);
                    classLabels  = userclassnames(maxscoreidx);
                    fittedLabels = labels(classLabels); % return the same type passed by the user
                    
                    if ~hasSoftLabels % for hard labels
                        labeledDataScores = dummyY;
                    else
                        normalizingSum = sum(labeledDataScores,2);
                        if normalizingSum~=0
                            labeledDataScores = labeledDataScores./normalizingSum; % normalize
                        end
                    end
                end
            end
            classSummary.ClassNames = userclassnames;
            classSummary.NonzeroProbClasses = nonzeroClassNames;
            mdl = SemiSupervisedGraphModel(fittedLabels, fittedLabelScores, fittermethod,...
                trainX, labeledDataScores, ...
                dataSummary, ...
                classSummary, distance, kernelscale);
        end
    end % end hidden static methods
    
    
    methods
        function [predlabels, predscores] = predict(this,X)
            %   PREDICT Predict response of the SemiSupervisedGraphModel.
            %   [LABEL,SCORE]=PREDICT(MDL,X) returns predicted class
            %   labels and scores from the model, MDL and predictors X. X
            %   must be a table if MDL was originally trained on a table,
            %   or a numeric matrix if MDL was originally trained on a
            %   matrix. If X is a table, it must contain all the predictors
            %   used for training this model. If X is a matrix, it must
            %   have P columns, where P is the number of predictors used
            %   for training this model. Classification labels LABEL have
            %   the same type as Y used for training. Scores SCORE are an
            %   N-by-K numeric matrix for N observations and K classes. The
            %   predicted label is assigned to the class with the largest
            %   score.
            
            X = classreg.learning.internal.table2PredictMatrix(X,[],[],...
                this.VariableRange,...
                this.CategoricalPredictors,this.PredictorNames);
            % check size post-conversion
            D = numel(this.PredictorNames);
            Dpassed = size(X,2);
            str = getString(message('stats:classreg:learning:classif:ClassificationModel:predictEmptyX:columns'));
            if Dpassed~=D
                error(message('stats:classreg:learning:classif:ClassificationModel:predictEmptyX:XSizeMismatch', D, str));
            end
            
            if isempty(X) % delegate to base class
                [predlabels, predscores] = predictEmptyX(this);
                return;
            end
            
            standardizeBeforePredict = ~isempty(this.Mu);
            if  standardizeBeforePredict
                X = semisupervised.Model.standardizePredictData(X,this.Mu,this.Sigma);
            end
            
            dist = this.Distance;
            
            if strcmpi(dist, 'goodall3') ...
                    || strcmpi(dist, 'ofd')
                % use heteropdist2 for a mix of categorical and continuous
                % predictors
                % edge cases covered for categorical predictors: inf and
                % nan
                %
                % inf: finite ismember indices are used as categorical
                % levels because heteropdist2 does not know how to handle
                % inf as a category because it expects positive integers as
                % categorical levels.
                
                % nan:  if a category is not seen during training, it will
                % be encoded as nan by table2FitMatrix- heteropdist2 uses
                % Laplace correction to calculate frequency of such unseen
                % categories as 1/N+1.
                %
                distmat = internal.stats.heteropdist2(this.PrivX, X,...
                    dist, this.CategoricalPredictors); % training data has to be passed first due to the convention adopted by heteropdist2
                distmat = distmat';
            else
                % for distances that pdist2 supports, delegate distance
                % calculation to pdist2
                distmat = pdist2(X, this.PrivX, dist);
            end
            
            kernelscale = this.KernelScale;
            % in the future, even the weighting function can be custom
            weightingfunction = @(d) exp(-(d./kernelscale).^2); % use the same weighting function used during training
            W = weightingfunction(distmat);
            F = [this.LabeledDataScore; this.LabelScores];
            predscores = (W*F)./sum(W,2); % weighted sum of training scores
            [~,maxscoreidx] = max(predscores,[],2);
            cnames = this.ClassSummary.ClassNames;
            predlabels = labels(cnames(maxscoreidx)); % call the labels method
        end
    end % end methods
    
    methods (Access = protected) % add properties to display
        function propStruct = propsForDisplay(this)
            propStruct = propsForDisplay@semisupervised.Model(this); % get common properties
            propStruct.('Method') = this.Method;  % add method to the list of displayed properties
        end
    end
end


function FU = labelpropagation( W, dummyY, numLabeled, numTotal, numClasses, iterationlimit, tolerance)
FU = zeros(numTotal-numLabeled, numClasses);
FUOld = zeros(numTotal-numLabeled, numClasses);
FL = dummyY(1:numLabeled,:);
iter = 0;
T = W./sum(W,2);
Tul = T(numLabeled+1:end, 1:numLabeled);
Tuu = T(numLabeled+1:end, numLabeled+1:end);

while iter < iterationlimit
    FU =  Tuu*FU + Tul*FL;
    criterion = sum(abs(FUOld - FU),2);
    if criterion < tolerance
        break;
    end
    FUOld = FU;
    iter = iter+1;
end
end

function YU = labelpropagationexact(W, dummyY, numLabeled, numTotal)
T = W./sum(W,2); % (D^-1)W
Tul = T(numLabeled+1:end, 1:numLabeled);
Tuu = T(numLabeled+1:end, numLabeled+1:end);
n = numTotal - numLabeled; % number of unlabeled points
Tuu(1:n+1:end) = Tuu(1:n+1:end) - 1;
YU = -Tuu\(Tul*dummyY);
end

function F = labelspreading(W, dummyY, numLabeled, numTotal, numClasses, alpha, iterationlimit, tolerance)
F = zeros(numTotal, numClasses);
Fold = F;
Ytotal = zeros(numTotal,numClasses);
Ytotal(1:numLabeled,:) = dummyY;
rowsum = sum(W,2);
sqrtD = spdiags(1./sqrt(rowsum),0, numTotal, numTotal);
W = (sqrtD)*W*(sqrtD);
W(1:numTotal+1:end) = 0; % set diagonal to zero
iter = 0;
while iter < iterationlimit
    F = alpha.*(W*F) + (1 - alpha).*Ytotal;
    criterion = sum(abs(Fold - F),2);
    if criterion < tolerance
        break;
    end
    Fold = F;
    iter = iter+1;
end
end

function F = labelspreadingexact(W, dummyY, numLabeled, numTotal, numClasses, alpha)
Ytotal = zeros(numTotal,numClasses);
Ytotal(1:numLabeled,:) = dummyY;
rowsum = sum(W,2);
sqrtD = spdiags(1./sqrt(rowsum),0, numTotal, numTotal);
W = (sqrtD)*W*(sqrtD);
W(1:numTotal+1:end) = 0; % set diagonal to zero
AlphaWminusI = alpha*W;
AlphaWminusI(1:numTotal+1:end) = AlphaWminusI(1:numTotal+1:end) -1;  
F = -AlphaWminusI\Ytotal; % I-alpha*W\Ytotal
end

