function stuff = kmtrain(Xin,nmean,init,yin)

if nargin<2  nmean = 8;  end; % make eight means

[n,d] = size(Xin);

% -- initialize the means somehow ...

% ... using a scaled random set
if init == 1
    xscale = max(Xin,[],1);

    means=[];
    for i=1:nmean
        means(i,:) = rand(1,d) .* xscale;
    end;
end;

% ... picking points at random
if init == 2
    for i=1:nmean
        means(i,:) = Xin(floor(n*rand),:);
    end;
end;


% ... picking one point from each class
if init == 3
    dx=1;
    for i=min(yin):max(yin)
        idxl = find(yin==i);
        if(~isempty(idxl))
            means(dx,:) = ...
                 Xin(idxl(ceil(length(idxl)*rand)),:);
        end;
        dx = dx+1;
    end;
end;


% ... taking the raw means of each class
if init == 4
    dx=1;
    for i=min(yin):max(yin)
        idxl = find(yin==i);
        if(~isempty(idxl))
            means(dx,:) = mean(Xin(idxl,:),1);
        end;
        dx = dx+1;
    end;
end;


% -- do the training
lmeans = zeros(nmean,d);
cont = 1;
soft_assn = [];
lsoft_assn = ones(1,n);

while(cont)

    %Estep
    for i=1:n
        for j=1:nmean
            dist(j) = distance(Xin(i,:),means(j,:));
        end;
        [y,idx] = min(dist);
        soft_assn(i) = idx(1);
    end;
    
    if (sum(abs(lsoft_assn-soft_assn)) == 0)
        cont = 0;
    else
        fprintf('%d points changed\n',...
            sum(abs(lsoft_assn-soft_assn)));
        lsoft_assn = soft_assn;
        
        %Mstep
        for j=1:nmean
            jdx = find(soft_assn==j);
            tassn(j) = length(jdx);
            if (~isempty(jdx))
                means(j,:) = mean(Xin(jdx,:),1);
            else
                fprintf('none were assigned to %d\n',j); 
                means(j,:) = xscale .* rand(1,d);
            end;
        end;
        tassn
    end;        
end;

stuff{1} = means;
stuff{2} = soft_assn;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% -- find the distance between a review and a mean
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dist = distance(rev,mean)
dist = sqrt(sum((rev-mean).^2));