function [centroids, costAssignment, iter] = kModesFunction(k, maxIter, choices, nPrespec, mapNonfit, numNonfitSubj, exoSubj)


[nPoints, nProfiles] = size(choices);
%nObs = nPoints - nPrespec;
points = 1:nPoints;



%% Randomly choose initial clusters
%  among subjects who do not follow a prespecified rule exactly

if numel(exoSubj) == 0
    randSubj = randsample(numNonfitSubj, k);

    selectedSubj = zeros(numNonfitSubj,1);
    selectedSubj(randSubj) = 1;

    % centroidsVec is a 0/1 vector that denotes for each observation
    % whether or not it is a centroid
    
    % centroidsIndex is an integer vector that denotes the positions of the
    % centroid observations
    
    centroidsVec = mapNonfit*selectedSubj;
    centroidsVec(1:nPrespec) = 1;

    centroidsIndex = find(centroidsVec);
    centroidsIndex = centroidsIndex(nPrespec+1: numel(centroidsIndex));
else
    centroidsIndex = nPrespec + exoSubj;
    centroidsVec = zeros(nPoints, 1);
    centroidsVec(1:nPrespec) = 1;
    centroidsVec(centroidsIndex) = 1;
end




%% Assign each observation to the closest centroid

assignment = points;

for subj = nPrespec+1:nPoints
    subjChoices = repmat(choices(subj,:),[nPoints 1]);
    dist = subjChoices ~= choices; 
    dist = sum(dist,2);
    dist(~centroidsVec) = 1000;     % Exclude all non-centroid or prespecified observations from considerations by making the distance very high
    [~, assignment(subj)] = min(dist);   
end



centroids = zeros(k, nProfiles);

% Calculate mode for each cluster

for i = 1 : k
    clusterObservations = points(assignment == centroidsIndex(i));
    centroids(i,:) = mode(choices(clusterObservations, :), 1);
end

centroids = [choices(1:nPrespec,:); centroids];


%% Iterate to improve centroids

% Assign each observation to the closest centroid

assignment = points;
cost = zeros(size(points));

iter = 0; 
change = 1;

while change ~= 0 && iter < maxIter

    for subj = nPrespec+1:nPoints
        subjChoices = repmat(choices(subj,:),[nPrespec + k, 1]);
        dist = subjChoices ~= centroids; 
        dist = sum(dist,2);
        [cost(subj), assignment(subj)] = min(dist);   
    end

    % Calculate mode for each cluster and update centroids

    centroidsIndex = [nPrespec+1: nPrespec+k];
    centroidsNew = centroids;
    for i = 1 : k
        clusterObservations = points(assignment == centroidsIndex(i));
        centroidsNew(nPrespec + i,:) = mode(choices(clusterObservations, :), 1);
    end

    change = +max(centroidsNew ~= centroids, [], 'all');

    iter = iter + 1;

    % disp(iter) 

    costAssignment = sum(cost(nPrespec+1:nPoints));

    centroids = centroidsNew;

end