function [Corr]=pointmatchfun(points,labels,maxpoints)

npoints=length(points);
nlabels=length(labels);

maxpoints=1;

%computation of the affinity matrix for the point matching starts here
Diff1 = distance(points',points');
Diff2 = distance(labels',labels');

redmat=ones(npoints*nlabels); %it is used to reduce the affinity between in consistent labelings. 
plmat=ones(npoints*nlabels);  % matrix of compatibility of assigments of size (npoints x nlabels)^2, where an assigment is a pair (point, label)
mutmat=zeros(npoints*nlabels); %mutex matrix
for l1=1:nlabels
    for p1=1:npoints
        for l2=1:nlabels
            for p2=1:npoints
                plmat(p1+(l1-1)*npoints,p2+(l2-1)*npoints)= abs(Diff1(p1,p2)-Diff2(l1,l2));
                if (p1==p2 && l1~=l2) || (p1~=p2 && l1==l2) %we try here to reduce the chance of a given point getting two labels
                    redmat(p1+(l1-1)*npoints,p2+(l2-1)*npoints)=0;
                    mutmat(p1+(l1-1)*npoints,p2+(l2-1)*npoints)=1; %mutex matrix that enforces 1-1 constraints
                end
            end
        end
    end
end

spaMutmat = spatialMutex(labels, points, 0.3);
mutmat=max(spaMutmat,spaMutmat')+mutmat;

[SimMatrix]=affinityMatrix(plmat);
A=SimMatrix.*redmat; %we prefer that a site has at most one label

% %adding absorbing state
% meanA=mean(A(:));
% A=[A zeros(npoints*nlabels, 1)];
% A=[A' zeros(npoints*nlabels+1, 1)+meanA];
% A(end,end)=1;

D=diag(sum(A,1));
A=inv(D)*A; %normalization to a stochastic matrix
%A(end,:)=[]; A(:,end)=[]; %getting back to the original matrix size
sizeA=npoints*nlabels;
A=min(A, ones(sizeA)-eye(sizeA)); %set diagonal entries to zero
figure; imagesc(A); title('nromalized matrix A');

%%%%%%%%%%%%%%%%%%%%% Key Step!!!!!!!!!
%random walk compputation, input A is a stochastic matrix
R=inv(eye(sizeA) - A);

% %uncoment for iterative computation of R as alternative to the line above
% invmutmat= ones(npoints*nlabels) - mutmat; %invmutmat has zero in mutex pairs and ones everywhere else
% P=eye(npoints*nlabels); S=zeros(npoints*nlabels);
% for i=1:20
%    P=(P*A).*invmutmat;
%    S=(S+P).*invmutmat;
% end
% R=S;
% %iterative computation ends here

%figure; imagesc(R); title('matrix R');

matf=zeros(sizeA);
matX=zeros(sizeA);
invmutmat= ones(sizeA) - mutmat; %invmutmat has zero in mutex pairs and ones everywhere else
f=R(1,:); %the solution f is the likelihood of belonging to MWS under the assumpion that node 1 belongs to MWS
%%%computes connected subgraph that does not violate mutex with max sum of f values
for i=1:1%sizeA
    f=R(i,:);
    [fv, fi] = sort(f,'descend');
    X=zeros(1,sizeA);
    X(fi(1))=1;
    for i=1:sizeA-1 %computes connected subgraph that does not violate mutex with max sum of f values
        Z=find(X==1);
        nozero=1;
        for k=1:length(Z)
            if invmutmat(Z(k),fi(i+1))==0 || A(Z(k),fi(i+1))==0
                nozero=0;
            end
        end
        if nozero>0 %a new node is added if it does not violate any mutex constraints with current nodes
            X(fi(i+1))=1;
        end
    end
    sumf(i)=sum(f.*X);
    matf(i,:)=f;
    matX(i,:)=X;
end
[sv, si]=max(sumf);
f = matf(si,:);
X=matX(si,:);


% %vizualization
% CorrProb = reshape(f,npoints,nlabels);
% figure;
% imagesc(CorrProb); title('f values of random walk');
% %vizualization
% CorrProb=X; 
% CorrProb=reshape(CorrProb,npoints,nlabels);
% figure; imagesc(CorrProb); title('corr by random walk');


%%%%%%%%%%%%% Solution of our quadratic program with mutex constraints
w = A/max(A(:)) - mutmat;
 x = zeros(sizeA, 1);
    x(1) = 1;
    [fvalue,x_star] = MIQP_solver(w, x ,200);
   
 

%Corr=reshape(x_star,npoints,nlabels);
    %vizualization 
CorrProb=reshape(x_star,npoints,nlabels)
figure; imagesc(CorrProb); title('corr by MIPQ');
CorrProb
Corr = [];
for i=1:size(CorrProb,2)
    Corr = [Corr find(CorrProb(:,i))];
end;
Corr


function d = distance(a,b)
if (nargin ~= 2)
   error('Not enough input arguments');
end

if (size(a,1) ~= size(b,1))
   error('A and B should be of same dimensionality');
end

aa=sum(a.*a,1); bb=sum(b.*b,1); ab=a'*b; 
d = sqrt(abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab));




%made by Tianyang Ma

function  spaMutmat = spatialMutex(pts1, pts2, thresh)

npts1 = size(pts1, 1);
npts2 = size(pts2, 1);

spaMutmat = zeros(npts1*npts2, npts1*npts2, 'double');
pts_dist1 = sqrt(distance(pts1, pts1));
pts_dist2 = sqrt(distance(pts2, pts2));

for i=1:npts1
    for m=1:npts2
        for j=1:npts1
            for n=1:npts2
                if(~(i==j&&m==n))
                    lr_diff1 = pts1(i, 1) - pts1(j, 1);
                    lr_diff2 = pts2(m, 1) - pts2(n, 1);
                    ud_diff1 = pts1(i, 2) - pts1(j, 2);
                    ud_diff2 = pts2(m, 2) - pts2(n, 2);
                    
                    flag1 = 0;
                    flag2 = 0;
                    
                    if(abs(lr_diff1) > thresh && abs(lr_diff2)>thresh)
                        flag1 = lr_diff1*lr_diff2<0;
                    end
                    if(abs(ud_diff1) > thresh && abs(ud_diff2)>thresh)
                        flag2 = ud_diff1*ud_diff2<0;
                    end
                    if(flag1||flag2)
                        spaMutmat((i-1)*npts2 + m, (j-1)*npts2 + n) = 1;
                    end                                    
                end
            end
        end
    end
end


function [W]=affinityMatrix(Diff);
% computes an affinity matrix for a given distance matrix

Diff=(Diff+Diff')/2;
[T,INDEX]=sort(Diff,2);
[m,n]=size(Diff);
W=zeros(m,n);
K=15;

for i=1:m
    for j=1:n
        SIGMA=mean([T(i,2:K),T(j,2:K)]);
        W(i,j)=normpdf(Diff(i,j),0,SIGMA);
    end
end

% W=normpdf(Real_Diff,0,14.2);% the highest is 13.8, the accuracy is 89.78%
% [m,n]=size(W);
   % Norm=repmat(sum(W')',1,n);
   % P=W./Norm;% normalization
   return
%0.35*SIGMA for ASC (the newest inner distance) matrix for K=17
%0.25*SIGMA for inner distance matrix for K=14
%0.4*SIGMA for shape part matrix for K=10
%0.2*SIGMA for Face(All) for K=5

% MIQP solver for maximum subgraph with mutex constraints
% Tianyang Ma and Longin Jan Latecki
% maximize x^Twx 
% s.t. 0<=x_i<=1 



%% MIQP solver for maximum subgraph with mutex constraints
% Tianyang Ma and Longin Jan Latecki
% maximize x^Twx 
% s.t. 0<=x_i<=1 

function [fvalue,x_star] = MIQP_solver(w, x0, MaxIter)
% w:  A - \gamma M, A is affinity matrix, and M is mutex matrix
% x0: the initial solution (has to satisfy mutex constraints)
% MaxIter: maximum num of iterations
thres = -inf;
%function definitions
f = @(x) x'*w*x; %the target function
itercount = 0;
% initialization
y_t = x0;
x_star = x0;
gofy = inf;
% iteration to update
while  (gofy > thres) && (itercount <MaxIter)
    itercount = itercount+1;
    x_bar = w'*y_t>=0; % it is discrete
    if(norm(x_bar- y_t)==0)
        x_star = x_bar;
        break;
    end
    if f(x_bar)>f(y_t)
        y_t1 = x_bar;
        x_star = x_bar;
    else
        eta = -((w'*y_t)'*(x_bar-y_t))/(((x_bar-y_t)'*w*(x_bar-y_t))+eps);
        y_t1 = y_t + eta*(x_bar - y_t);
    end
    gofy = norm(y_t1 - y_t);
    y_t    = y_t1;
end
% figure, plot(y_t);
% fprintf('\nitercount = %d,  gofy = %f, sum(x) = %d\n', itercount, gofy, sum(x_star));
fvalue = f(x_star);



