function [fitted,vaf,weights,end_condition] = partitionfit(prox,member)

% PARTITIONFIT provides a least-squares approximation to a proximity 
% matrix based on a given collection of partitions.
%
% syntax: [fitted,vaf,weights,end_condition] = partitionfit(prox,member)
%
% PROX is the n x n input proximity matrix (with a zero main diagonal
% and a dissimilarity interpretation); MEMBER is the m x n matrix
% indicating cluster membership, where each row corresponds to a specific
% partition (there are m partitions in general); the columns of MEMBER
% are in the same input order used for PROX.
% FITTED is an n x n matrix fitted to PROX (through least-squares)
% constructed from the nonnegative weights given in the m x 1 WEIGHTS
% vector corresponding to each of the partitions.  VAF is the variance-
% accounted-for in the proximity matrix PROX by the fitted matrix FITTED.
% END_CONDITION should be zero for a normal termination of the optimization
% process.

n = size(prox,1);
nsize = n*(n-1)/2;
npart = size(member,1);

partvec = zeros(npart,nsize);
datavec = zeros(1,nsize);
fitted = zeros(n,n);
proxave = zeros(n,n);

index = 0;

for i = 1:(n-1)
    for j = (i+1):n
        
        index = index + 1;        
        datavec(index) = prox(i,j);
        
        for k = 1:npart
            
            if(member(k,i) ~= member(k,j))
                 partvec(k,index) = 1;
            end
            
        end
    end
end

identity = eye(nsize);
equality_flag = zeros(npart,1);
constraint_constant = zeros(npart,1);

[solution,kuhn_tucker,iterations,end_condition] = ...
    dykstra(datavec,identity,partvec,constraint_constant,equality_flag);

weights = (.5)*kuhn_tucker;

fittedvec = datavec - solution';

index = 0;

for i = 1:(n-1)
    for j = (i+1):n
        
        index = index + 1;
        fitted(i,j) = fittedvec(index);
        fitted(j,i) = fitted(i,j);
        
    end
end

aveprox = sum(sum(prox))/(n*(n-1));

for i = 1:n
    for j = 1:n
        
        if(i ~= j)
            proxave(i,j) = aveprox;
        else
            proxave(i,j) = 0;
        end
        
    end
end

diff = sum(sum((prox - fitted).^2));
denom = sum(sum((prox - proxave).^2));

vaf = 1 - (diff/denom);