function label=init_methods(data, K, choose)
% 输入：无标签数据，聚类数，选择方法
% 输出：聚类标签
if choose==1
    %随机初始化，随机选K行作为聚类中心，并用欧氏距离计算其他点到其聚类，将数据集分为K类，输出每个样例的类标签
    [X_num, ~]=size(data);
    rand_array=randperm(X_num);    %产生1~X_num之间整数的随机排列   
    para_miu=data(rand_array(1:K), :);  %随机排列取前K个数，在X矩阵中取这K行作为初始聚类中心
    %欧氏距离，计算（X-para_miu）^2=X^2+para_miu^2-2*X*para_miu\'，矩阵大小为X_num*K
    distant=repmat(sum(data.*data,2),1,K)+repmat(sum(para_miu.*para_miu,2)\',X_num,1)-2*data*para_miu\';
    %返回distant每行最小值所在的下标
    [~,label]=min(distant,[],2);
elseif choose==2
    %用kmeans进行初始化聚类，将数据集聚为K类，输出每个样例的类标签
    label=kmeans(data, K);
elseif choose==3
    %用FCM算法进行初始化
    options=[NaN, NaN, NaN, 0];
    [~, responsivity]=fcm(data, K, options);   %用FCM算法求出隶属度矩阵
    [~, label]=max(responsivity\', [], 2);
elseif choose==4
    label = litekmeans(data, K,\'Replicates\',20);
end

function [label, center, bCon, sumD, D] = litekmeans(X, k, varargin)
% [IDX, C] = litekmeans(data, K,\'Replicates\',20);
%LITEKMEANS K-means clustering, accelerated by matlab matrix operations.
%
%   label = LITEKMEANS(X, K) partitions the points in the N-by-P data matrix
%   X into K clusters.  This partition minimizes the sum, over all
%   clusters, of the within-cluster sums of point-to-cluster-centroid
%   distances.  Rows of X correspond to points, columns correspond to
%   variables.  KMEANS returns an N-by-1 vector label containing the
%   cluster indices of each point.
%
%   [label, center] = LITEKMEANS(X, K) returns the K cluster centroid
%   locations in the K-by-P matrix center.
%
%   [label, center, bCon] = LITEKMEANS(X, K) returns the bool value bCon to
%   indicate whether the iteration is converged.  
%
%   [label, center, bCon, SUMD] = LITEKMEANS(X, K) returns the
%   within-cluster sums of point-to-centroid distances in the 1-by-K vector
%   sumD.    
%
%   [label, center, bCon, SUMD, D] = LITEKMEANS(X, K) returns
%   distances from each point to every centroid in the N-by-K matrix D. 
%
%   [ ... ] = LITEKMEANS(..., \'PARAM1\',val1, \'PARAM2\',val2, ...) specifies
%   optional parameter name/value pairs to control the iterative algorithm
%   used by KMEANS.  Parameters are:
%
%   \'Distance\' - Distance measure, in P-dimensional space, that KMEANS
%      should minimize with respect to.  Choices are:
%            {\'sqEuclidean\'} - Squared Euclidean distance (the default)
%             \'cosine\'       - One minus the cosine of the included angle
%                              between points (treated as vectors). Each
%                              row of X SHOULD be normalized to unit. If
%                              the intial center matrix is provided, it
%                              SHOULD also be normalized.
%
%   \'Start\' - Method used to choose initial cluster centroid positions,
%      sometimes known as "seeds".  Choices are:
%         {\'sample\'}  - Select K observations from X at random (the default)
%          \'cluster\' - Perform preliminary clustering phase on random 10%
%                      subsample of X.  This preliminary phase is itself
%                      initialized using \'sample\'. An additional parameter
%                      clusterMaxIter can be used to control the maximum
%                      number of iterations in each preliminary clustering
%                      problem.
%           matrix   - A K-by-P matrix of starting locations; or a K-by-1
%                      indicate vector indicating which K points in X
%                      should be used as the initial center.  In this case,
%                      you can pass in [] for K, and KMEANS infers K from
%                      the first dimension of the matrix.
%
%   \'MaxIter\'    - Maximum number of iterations allowed.  Default is 100.
%
%   \'Replicates\' - Number of times to repeat the clustering, each with a
%                  new set of initial centroids. Default is 1. If the
%                  initial centroids are provided, the replicate will be
%                  automatically set to be 1.
%
% \'clusterMaxIter\' - Only useful when \'Start\' is \'cluster\'. Maximum number
%                    of iterations of the preliminary clustering phase.
%                    Default is 10.  
%
%
%    Examples:
%
%       fea = rand(500,10);
%       [label, center] = litekmeans(fea, 5, \'MaxIter\', 50);
%
%       fea = rand(500,10);
%       [label, center] = litekmeans(fea, 5, \'MaxIter\', 50, \'Replicates\', 10);
%
%       fea = rand(500,10);
%       [label, center, bCon, sumD, D] = litekmeans(fea, 5, \'MaxIter\', 50);
%       TSD = sum(sumD);
%
%       fea = rand(500,10);
%       initcenter = rand(5,10);
%       [label, center] = litekmeans(fea, 5, \'MaxIter\', 50, \'Start\', initcenter);
%
%       fea = rand(500,10);
%       idx=randperm(500);
%       [label, center] = litekmeans(fea, 5, \'MaxIter\', 50, \'Start\', idx(1:5));
%
%
%   See also KMEANS
%
%    [Cite] Deng Cai, "Litekmeans: the fastest matlab implementation of
%           kmeans," Available at:
%           http://www.zjucadcg.cn/dengcai/Data/Clustering.html, 2011. 
%
%   version 2.0 --December/2011
%   version 1.0 --November/2011
%
%   Written by Deng Cai (dengcai AT gmail.com)

if nargin < 2
    error(\'litekmeans:TooFewInputs\',\'At least two input arguments required.\');
end

[n, p] = size(X);

pnames = {   \'distance\' \'start\'   \'maxiter\'  \'replicates\' \'onlinephase\' \'clustermaxiter\'};
dflts =  {\'sqeuclidean\' \'sample\'       []        []        \'off\'              []        };
[eid,errmsg,distance,start,maxit,reps,online,clustermaxit] = getargs(pnames, dflts, varargin{:});
if ~isempty(eid)
    error(sprintf(\'litekmeans:%s\',eid),errmsg);
end

if ischar(distance)
    distNames = {\'sqeuclidean\',\'cosine\'};
    j = strcmpi(distance, distNames);
    j = find(j);
    if length(j) > 1
        error(\'litekmeans:AmbiguousDistance\', ...
            \'Ambiguous \'\'Distance\'\' parameter value:  %s.\', distance);
    elseif isempty(j)
        error(\'litekmeans:UnknownDistance\', ...
            \'Unknown \'\'Distance\'\' parameter value:  %s.\', distance);
    end
    distance = distNames{j};
else
    error(\'litekmeans:InvalidDistance\', ...
        \'The \'\'Distance\'\' parameter value must be a string.\');
end

center = [];
if ischar(start)
    startNames = {\'sample\',\'cluster\'};
    j = find(strncmpi(start,startNames,length(start)));
    if length(j) > 1
        error(message(\'litekmeans:AmbiguousStart\', start));
    elseif isempty(j)
        error(message(\'litekmeans:UnknownStart\', start));
    elseif isempty(k)
        error(\'litekmeans:MissingK\', ...
            \'You must specify the number of clusters, K.\');
    end
    if j == 2
        if floor(.1*n) < 5*k
            j = 1;
        end
    end
    start = startNames{j};
elseif isnumeric(start)
    if size(start,2) == p
        center = start;
    elseif (size(start,2) == 1 || size(start,1) == 1)
        center = X(start,:);
    else
        error(\'litekmeans:MisshapedStart\', ...
            \'The \'\'Start\'\' matrix must have the same number of columns as X.\');
    end
    if isempty(k)
        k = size(center,1);
    elseif (k ~= size(center,1))
        error(\'litekmeans:MisshapedStart\', ...
            \'The \'\'Start\'\' matrix must have K rows.\');
    end
    start = \'numeric\';
else
    error(\'litekmeans:InvalidStart\', ...
        \'The \'\'Start\'\' parameter value must be a string or a numeric matrix or array.\');
end

% The maximum iteration number is default 100
if isempty(maxit)
    maxit = 100;
end

% The maximum iteration number for preliminary clustering phase on random
% 10% subsamples is default 10 
if isempty(clustermaxit)
    clustermaxit = 10;
end


% Assume one replicate
if isempty(reps) || ~isempty(center)
    reps = 1;
end

if ~(isscalar(k) && isnumeric(k) && isreal(k) && k > 0 && (round(k)==k))
    error(\'litekmeans:InvalidK\', ...
        \'X must be a positive integer value.\');
elseif n < k
    error(\'litekmeans:TooManyClusters\', ...
        \'X must have more rows than the number of clusters.\');
end

bestlabel = [];
sumD = zeros(1,k);
bCon = false;

for t=1:reps
    switch start
        case \'sample\'
            center = X(randsample(n,k),:);
        case \'cluster\'
            Xsubset = X(randsample(n,floor(.1*n)),:);
            [dump, center] = litekmeans(Xsubset, k, varargin{:}, \'start\',\'sample\', \'replicates\',1 ,\'MaxIter\',clustermaxit);
        case \'numeric\'
    end
    
    last = 0;label=1;
    it=0;
    
    switch distance
        case \'sqeuclidean\'
            while any(label ~= last) && it<maxit
                last = label;
                
                bb = full(sum(center.*center,2)\');
                ab = full(X*center\');
                D = bb(ones(1,n),:) - 2*ab;
                
                [val,label] = min(D,[],2); % assign samples to the nearest centers
                ll = unique(label);
                if length(ll) < k
                    %disp([num2str(k-length(ll)),\' clusters dropped at iter \',num2str(it)]);
                    missCluster = 1:k;
                    missCluster(ll) = [];
                    missNum = length(missCluster);
                    
                    aa = sum(X.*X,2);
                    val = aa + val;
                    [dump,idx] = sort(val,1,\'descend\');
                    label(idx(1:missNum)) = missCluster;
                end
                E = sparse(1:n,label,1,n,k,n);  % transform label into indicator matrix
                center = full((E*spdiags(1./sum(E,1)\',0,k,k))\'*X);    % compute center of each cluster
                it=it+1;
            end
            if it<maxit
                bCon = true;
            end
            if isempty(bestlabel)
                bestlabel = label;
                bestcenter = center;
                if reps>1
                    if it>=maxit
                        aa = full(sum(X.*X,2));
                        bb = full(sum(center.*center,2));
                        ab = full(X*center\');
                        D = bsxfun(@plus,aa,bb\') - 2*ab;
                        D(D<0) = 0;
                    else
                        aa = full(sum(X.*X,2));
                        D = aa(:,ones(1,k)) + D;
                        D(D<0) = 0;
                    end
                    D = sqrt(D);
                    for j = 1:k
                        sumD(j) = sum(D(label==j,j));
                    end
                    bestsumD = sumD;
                    bestD = D;
                end
            else
                if it>=maxit
                    aa = full(sum(X.*X,2));
                    bb = full(sum(center.*center,2));
                    ab = full(X*center\');
                    D = bsxfun(@plus,aa,bb\') - 2*ab;
                    D(D<0) = 0;
                else
                    aa = full(sum(X.*X,2));
                    D = aa(:,ones(1,k)) + D;
                    D(D<0) = 0;
                end
                D = sqrt(D);
                for j = 1:k
                    sumD(j) = sum(D(label==j,j));
                end
                if sum(sumD) < sum(bestsumD)
                    bestlabel = label;
                    bestcenter = center;
                    bestsumD = sumD;
                    bestD = D;
                end
            end
        case \'cosine\'
            while any(label ~= last) && it<maxit
                last = label;
                W=full(X*center\');
                [val,label] = max(W,[],2); % assign samples to the nearest centers
                ll = unique(label);
                if length(ll) < k
                    missCluster = 1:k;
                    missCluster(ll) = [];
                    missNum = length(missCluster);
                    [dump,idx] = sort(val);
                    label(idx(1:missNum)) = missCluster;
                end
                E = sparse(1:n,label,1,n,k,n);  % transform label into indicator matrix
                center = full((E*spdiags(1./sum(E,1)\',0,k,k))\'*X);    % compute center of each cluster
                centernorm = sqrt(sum(center.^2, 2));
                center = center ./ centernorm(:,ones(1,p));
                it=it+1;
            end
            if it<maxit
                bCon = true;
            end
            if isempty(bestlabel)
                bestlabel = label;
                bestcenter = center;
                if reps>1
                    if any(label ~= last)
                        W=full(X*center\');
                    end
                    D = 1-W;
                    for j = 1:k
                        sumD(j) = sum(D(label==j,j));
                    end
                    bestsumD = sumD;
                    bestD = D;
                end
            else
                if any(label ~= last)
                    W=full(X*center\');
                end
                D = 1-W;
                for j = 1:k
                    sumD(j) = sum(D(label==j,j));
                end
                if sum(sumD) < sum(bestsumD)
                    bestlabel = label;
                    bestcenter = center;
                    bestsumD = sumD;
                    bestD = D;
                end
            end
    end
end

label = bestlabel;
center = bestcenter;
if reps>1
    sumD = bestsumD;
    D = bestD;
elseif nargout > 3
    switch distance
        case \'sqeuclidean\'
            if it>=maxit
                aa = full(sum(X.*X,2));
                bb = full(sum(center.*center,2));
                ab = full(X*center\');
                D = bsxfun(@plus,aa,bb\') - 2*ab;
                D(D<0) = 0;
            else
                aa = full(sum(X.*X,2));
                D = aa(:,ones(1,k)) + D;
                D(D<0) = 0;
            end
            D = sqrt(D);
        case \'cosine\'
            if it>=maxit
                W=full(X*center\');
            end
            D = 1-W;
    end
    for j = 1:k
        sumD(j) = sum(D(label==j,j));
    end
end


function [eid,emsg,varargout]=getargs(pnames,dflts,varargin)
%GETARGS Process parameter name/value pairs 
%   [EID,EMSG,A,B,...]=GETARGS(PNAMES,DFLTS,\'NAME1\',VAL1,\'NAME2\',VAL2,...)
%   accepts a cell array PNAMES of valid parameter names, a cell array
%   DFLTS of default values for the parameters named in PNAMES, and
%   additional parameter name/value pairs.  Returns parameter values A,B,...
%   in the same order as the names in PNAMES.  Outputs corresponding to
%   entries in PNAMES that are not specified in the name/value pairs are
%   set to the corresponding value from DFLTS.  If nargout is equal to
%   length(PNAMES)+1, then unrecognized name/value pairs are an error.  If
%   nargout is equal to length(PNAMES)+2, then all unrecognized name/value
%   pairs are returned in a single cell array following any other outputs.
%
%   EID and EMSG are empty if the arguments are valid.  If an error occurs,
%   EMSG is the text of an error message and EID is the final component
%   of an error message id.  GETARGS does not actually throw any errors,
%   but rather returns EID and EMSG so that the caller may throw the error.
%   Outputs will be partially processed after an error occurs.
%
%   This utility can be used for processing name/value pair arguments.
%
%   Example:
%       pnames = {\'color\' \'linestyle\', \'linewidth\'}
%       dflts  = {    \'r\'         \'_\'          \'1\'}
%       varargin = {{\'linew\' 2 \'nonesuch\' [1 2 3] \'linestyle\' \':\'}
%       [eid,emsg,c,ls,lw] = statgetargs(pnames,dflts,varargin{:})    % error
%       [eid,emsg,c,ls,lw,ur] = statgetargs(pnames,dflts,varargin{:}) % ok

% We always create (nparams+2) outputs:
%    one each for emsg and eid
%    nparams varargs for values corresponding to names in pnames
% If they ask for one more (nargout == nparams+3), it\'s for unrecognized
% names/values

%   Original Copyright 1993-2008 The MathWorks, Inc. 
%   Modified by Deng Cai (dengcai@gmail.com) 2011.11.27

% Initialize some variables
emsg = \'\';
eid = \'\';
nparams = length(pnames);
varargout = dflts;
unrecog = {};
nargs = length(varargin);

% Must have name/value pairs
if mod(nargs,2)~=0
    eid = \'WrongNumberArgs\';
    emsg = \'Wrong number of arguments.\';
else
    % Process name/value pairs
    for j=1:2:nargs
        pname = varargin{j};
        if ~ischar(pname)
            eid = \'BadParamName\';
            emsg = \'Parameter name must be text.\';
            break;
        end
        i = strcmpi(pname,pnames);
        i = find(i);
        if isempty(i)
            % if they\'ve asked to get back unrecognized names/values, add this
            % one to the list
            if nargout > nparams+2
                unrecog((end+1):(end+2)) = {varargin{j} varargin{j+1}};
                % otherwise, it\'s an error
            else
                eid = \'BadParamName\';
                emsg = sprintf(\'Invalid parameter name:  %s.\',pname);
                break;
            end
        elseif length(i)>1
            eid = \'BadParamName\';
            emsg = sprintf(\'Ambiguous parameter name:  %s.\',pname);
            break;
        else
            varargout{i} = varargin{j+1};
        end
    end
end

varargout{nparams+1} = unrecog;

function data = normlization(data, choose)
% 数据归一化 
if choose==0
    % 不归一化
    data = data;
elseif choose==1
    % Z-score归一化
    data = bsxfun(@minus, data, mean(data));
    data = bsxfun(@rdivide, data, std(data));
elseif choose==2
    % 最大-最小归一化处理
    [data_num,~]=size(data);
    data=(data-ones(data_num,1)*min(data))./(ones(data_num,1)*(max(data)-min(data)));
end