function [membership,objectives,decrease,permmember] = partition_lincon_kmeans(prox,lincon,criterion)

% PARTITION_LINCON_KMEANS carries out a linearly constrained K-means
% cluster analysis based on dynamic programming.


% INPUT
% ------------------------------------------------------------------
% PROX       = the (n x n) proximity matrix of inter-object squared
%              Euclidean distances
% LINCON     = the linear constraint of the objects
% CRITERION  = if criterion == 1 the within-cluster sum-of-squared-
%               error is minimized. if criterion == 2 the sum of the
%               diameters of the clusters are minimized.

% OUTPUT
% ------------------------------------------------------------------
% MEMBERSHIP = the (n x n) matrix where the rows indicate the num-
%              ber of clusters and the columns indicate each
%               object, where the objects are in the original order
%               as indicated by the ordering in the data matrix.
% OBJECTIVES = the (n x 1) vector of objective function values,
%              where the nth values corresponds to the objective
%              function when partitioning the data into N clusters
%              based on LINCON
% DECREASE   = the (n x 1) vector indicating the percent decrease
%              in the objective function when moving from the ith
%              to the ith + 1 position of OBJECTIVES
% PERMMEMBER = the (n x n) matrix that is the result of permuting
%              the columns of MEMBERSHIP based on LINCON.
% ------------------------------------------------------------------


%initializing arrays
tic;
n = size(prox,1);
xi = (1.0e+20)*ones(n,n);
nxi = zeros(n,n);
clus = zeros(n,n);
proxpermut = prox(lincon,lincon);


%initial single subset processing and
%homogeneity calculations

for nlb = 1:n
    for nub = nlb:n
        
        index = 0.0;
        
        for i = nlb:nub
            for j = nlb:nub
                
                subsetsize = (nub - nlb) + 1;
                
%                 if  (index < proxpermut(i,j))
%                     index = proxpermut(i,j);
%                 end
                if criterion == 1
                    index = sum(sum(proxpermut(nlb:nub,nlb:nub)))/(2*subsetsize);
                end
                if criterion == 2
                    index = max(max(proxpermut(nlb:nub,nlb:nub)));
                end           
            end
        end
        
        if (nlb == 1)
            xi(1,nub) = index;
            nxi(1,nub) = 1;
        end
        
        clus(nlb,nub) = index;
    end
end

%processing the recursion

for k = 2:n
    
    for i = (k-1):(n-1)
        for jone = (i+1):n
            
            temp = xi(k-1,i) + clus(i+1,jone);
            
            if(temp < xi(k,jone))
                
                xi(k,jone) = temp;
                nxi(k,jone) = i + 1;
                
            end
        end
    end
end


%this does the membership for a given cluster size (not all subsets of that size)
membership=[];
objectives=[];
for j = 1:n
for k = j:j
    membershipvec=zeros(1,n);
    cluster = 0;

    
    nend = n;
    nk = k;
    objective = xi(nk,nend);
    
          
        
        add=lincon(nxi(nk,nend):nend);
        cluster = cluster+1;
        a = length(add);
        for i=1:a
            membershipvec(add(i)) = cluster;
        end
      
        
      
     while(nxi(nk,nend) ~= 1)
        
        nend = nxi(nk,nend) - 1;
        nk = nk - 1;
                        
        
        add=lincon(nxi(nk,nend):nend);
        cluster = cluster+1;
        a = length(add);
        for i=1:a
            membershipvec(add(i)) = cluster;
        end

        
    end
    cluster = cluster+1;
    
end
membership = [membership; membershipvec];

objectives = [objectives;objective];
end
plot(objectives)
permmember =  membership(:,lincon);

decrease=zeros(n,1);
for i=2:n
    if i<n
    decrease(i) = ((objectives(i-1) - objectives(i)))/(objectives(i-1));
    end
    if i==n
        decrease(i) = 1;
    end
end

toc