function [clusterlbs,idxbestcandpivots] = cluster_candpivots(candpivots,scores,adjmat_cp2cxpts,varargin)
% clusterlbs = cluster label for each candpivots
% idxbestcandpivots = index of the highest score candpivot for each cluster

% distance is weighted by scores to make it easier to cluster bad matches
% together, giving more clusters for good matches

if isempty(candpivots)
    clusterlbs = [];
    idxbestcandpivots = [];
    return;
end

pp = inputParser;
pp.addParamValue('distfunc',@srec_distfunc_norm2);
pp.addParamValue('usescoredist',false); %modify distance by score
pp.addParamValue('usecenterdist',false); %modify distance by candpivot distance
pp.addParamValue('useiter',false); %use iterative clustering?
pp.addParamValue('nclust',GLOBALCONST.THRES_NCLUSTER_PIVOT_PER_LABEL);
pp.parse(varargin{:});

distfunc = pp.Results.distfunc;
usescoredist = pp.Results.usescoredist;
nmaxcluster = pp.Results.nclust;
use_simple_clustering = ~pp.Results.useiter;

npivots = size(candpivots,2);
clusterlbs = 1:npivots;
datapts = double(full(adjmat_cp2cxpts));

ncluster_final = nmaxcluster;
cluster_ratio = 0.5;
ncluster_current = size(datapts,1);
idxactive = 1:npivots;
distmat = distfunc(adjmat_cp2cxpts);

%shorten the distances between bad matches so that we have more clusters in
%good matches
if usescoredist
    s1 = repmat(scores(:),1,length(scores));
    s2 = repmat(scores(:)',length(scores),1);
    smax = max(s1,s2);
    distmat = distmat.*smax;
end

if pp.Results.usecenterdist
    centerdistmat = squareform(pdist(candpivots'));
    distmat = distmat.*centerdistmat;
end

if use_simple_clustering
    z = linkage(squareform(distmat));
    clusterlbs = cluster(z,'maxclust',nmaxcluster);
else
    while ncluster_current > ncluster_final
        n = round(ncluster_current * cluster_ratio);
        n = max(n,ncluster_final);
        thisdatapts = datapts(idxactive,:);
        thiscores = scores(idxactive);
        prev_cluster_lbs = clusterlbs(idxactive);
        prev_max_cluster_lbs = max(prev_cluster_lbs);
        
        thisdistmat = distmat(idxactive,idxactive);
        z = linkage(squareform(thisdistmat));
        t = cluster(z,'maxclust',n);
        %     t = clusterdata(thisdatapts,'maxclust',n,'distance','cosine');
        uids = unique(t);
        idxkeep = zeros(1,length(uids));
        maskdone = false(size(clusterlbs));
        for i=1:length(uids)
            mask = t==uids(i);
            tmpscores = thiscores;
            tmpscores(~mask) = -inf;
            [~,idxmax] = max(tmpscores);
            idxkeep(i) = idxactive(idxmax);
            
            clbs_to_merge = prev_cluster_lbs(mask);
            tf = ismember(clusterlbs,clbs_to_merge);
            clusterlbs(tf) = prev_max_cluster_lbs + i;
            maskdone(tf) = 1;
        end
        
        ncluster_current=n;
        idxactive = idxkeep;
    end
end


% clusterlbs = clusterdata(datapts,'maxclust',ncluster_final,'distance','cosine');
[~,~,clusterlbs] = unique(clusterlbs);

ncluster = max(clusterlbs);
idxbestcandpivots = zeros(1,ncluster);
for i=1:ncluster
    mask = clusterlbs == i;
    thiscores = scores;
    thiscores(~mask) = -inf;
    [~,idxmax] = max(thiscores);
    idxbestcandpivots(i) = idxmax;
end

end

