function [membership,objectives,permmember,clusmeasure,...
    cluscoord,residsumsq] = orderpartitionfnd(prox,lincon)

% ORDERPARTITIONFND uses dynamic programming to
% construct a linearly constrained cluster analysis that
% consists of a collection of partitions with from 1 to
% n ordered classes.
%
% syntax: [membership,objectives,permmember,clusmeasure,...
%    cluscoord,residsumsq] = orderpartitionfnd(prox,lincon)
%
% PROX is the input proximity matrix (with a zero main diagonal
% and a dissimilarity interpretation); LINCON is the given
% constraining linear order (a permutation of the integers from
% 1 to n).
% MEMBERSHIP is the n x n matrix indicating cluster membership,
% where rows correspond to the number of ordered clusters, 
% and the columns are in the identity permutation input order 
% used for PROX. PERMMEMBER uses LINCON to reorder the columns
% of MEMBERSHIP.
% OBJECTIVES is the vector of merit values maximized in the
% construction of the ordered partitions; RESIDSUMSQ is the 
% vector of residual sum of squares obtained for the ordered
% partition construction.  CLUSMEASURE is the n x n matrix 
% (upper-triangular) containing the cluster measures for contiguous 
% object sets; the appropriate values in CLUSMEASURE are added
% to obtain the values optimized in OBJECTIVES; CLUSCOORD is also
% an n x n (upper-triangular) matrix but now containing the coordinates
% that would be would be used for all the (ordered) 
% objects within a class.

%initializing arrays

n = size(prox,1);
xi = zeros(n,n);
nxi = zeros(n,n);
clus = zeros(n,n);
proxpermut = prox(lincon,lincon);
cluscoord = zeros(n,n);


%initial single subset processing and
%homogeneity calculations

for nlb = 1:n
    for nub = nlb:n
        
                 
                subsetsize = (nub - nlb) + 1;              
               
               if((nlb == 1) & (nub < n))
                   
                   index = ...
                       ((sum(sum(proxpermut(nlb:nub,(nub+1):n))))^2)/subsetsize;
                   
                        coordindex = ...
                       -(sum(sum(proxpermut(nlb:nub,(nub+1):n))))/(subsetsize*n);
               end
               
               if((nlb == 1) & (nub == n))
                   
                   index = 0;
                   coordindex = 0;
                   
               end
               
               if ((nlb > 1) & (nub < n ))
                  
                   index = ...
                   ((sum(sum(proxpermut(nlb:nub,1:(nlb-1)))) - ...
                       sum(sum(proxpermut(nlb:nub,(nub+1):n))))^2)/subsetsize;
                   
                       coordindex = ...
                   (sum(sum(proxpermut(nlb:nub,1:(nlb-1)))) - ...
                       sum(sum(proxpermut(nlb:nub,(nub+1):n))))/(subsetsize*n);
                   
               end
               
               if((nlb > 1) & (nub == n))
                   
                   index = ...
                       ((sum(sum(proxpermut(nlb:nub,1:(nlb-1)))))^2)/subsetsize;
                   
                       coordindex = ...
                       (sum(sum(proxpermut(nlb:nub,1:(nlb-1)))))/(subsetsize*n);
               end
                  
     
        
        if (nlb == 1)
            xi(1,nub) = index;
            nxi(1,nub) = 1;
        end
        
        clus(nlb,nub) = index;
        cluscoord(nlb,nub) = coordindex;
    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



%retriving membership

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

permmember =  membership(:,lincon);

clusmeasure = clus;

residsumsq = ((1/2)*sum(sum(prox.^2))) - ((1/n)*objectives);



