%Implementation of RRWM, the codes are written based on
%ˇ°Reweighted Random Walks for Graph Matching, Minsu Cho et al, 2010ˇ±
addpath('RRWM\Methods\RRWM');
addpath('RRWM\utils');

[pair_4,pair_2,names]=pair_2_4(write_folder);

paircellI=0;
for i=1:size(pair_2,1)
    candis{i}.name1=pair_2{i,1};
    candis{i}.name2=pair_2{i,2};
    
    candis{i}.pt=load([write_folder...
        pair_2{i,1} '-' pair_2{i,2} ...
        '-ptcandi.txt']);
    candis{i}.line=load([write_folder...
        pair_2{i,1} '-' pair_2{i,2}...
        '-linecandi.txt']);
    maxl1=max(candis{i}.line(:,1));
    maxl2=max(candis{i}.line(:,2));
    
    candis{i}.all=[candis{i}.line;
        candis{i}.pt(:,1)+maxl1  candis{i}.pt(:,2)+maxl2];
    [group1,group2] = make_group12(candis{i}.all);
    
    [idx1,ID1] = make_groups(group1);
    [idx2,ID2] = make_groups(group2);
    
    if ID1(end) < ID2(end)
        [idx1, ID1, idx2, ID2, dumVal, dumSize] = make_groups_slack(idx1, ID1, idx2, ID2);
        dumDim = 1;
    elseif ID1(end) > ID2(end)
        [idx2,ID2,idx1,ID1,dumVal, dumSize] = make_groups_slack(idx2, ID2, idx1, ID1);
        dumDim = 2;
    else
        dumDim = 0; dumVal = 0; dumSize = 0;
    end
    
    idx1 = idx1-1; idx2 = idx2-1;
    
    candis{i}.dumDim=dumDim;
    candis{i}.dumVal=dumVal;
    candis{i}.dumSize=dumSize;
    candis{i}.ID1=ID1;
    candis{i}.ID2=ID2;
    candis{i}.idx1=idx1;
    candis{i}.idx2=idx2;
    
    candis{i}.stid=paircellI+1;
    candis{i}.edid=paircellI+size( candis{i}.all,1);
    
    
    paircellI=paircellI+size( candis{i}.all,1);
    
    pairwise{i}.pt=load([write_folder...
        pair_2{i,1} '-' pair_2{i,2}...
        '-pairwise_pt.txt']);
    pairwise{i}.line=load([write_folder...
        pair_2{i,1} '-' pair_2{i,2}...
        '-pairwise_line.txt']);
end
str=[write_folder 'pairwise.mat'];
save(str,'pairwise');


%Resconstruct the index
lineIst=size(candis{1}.line,1);
index_add=[0 lineIst];
lineIst=lineIst+size(candis{1}.pt,1);
for i=2:size(pair_2,1)
    index_add=[index_add;...
        lineIst lineIst+size(candis{i}.line,1)];
    lineIst=size(candis{i}.line,1)+size(candis{i}.pt,1)+lineIst;
end

%Triplet candidates
for i=1:size(pair_4,1)
    triple{i}.pt=load([write_folder...
        pair_4{i,1} '-' pair_4{i,2} '-' pair_4{i,4}...
        '-pm.txt']);
    triple{i}.line=load([write_folder...
        pair_4{i,1} '-' pair_4{i,2} '-' pair_4{i,4}...
        '-lm.txt']);
    triple{i}.name1=pair_4{i,1};
    triple{i}.name2=pair_4{i,2};
    triple{i}.name3=pair_4{i,3};
    triple{i}.name4=pair_4{i,4};
    
    for j=1:size(pair_2,1)
        
        if strcmp(pair_4{i,1},pair_2{j,1})&&...
           strcmp(pair_4{i,2},pair_2{j,2})
            triple{i}.line(:,1)=index_add(j,1)+triple{i}.line(:,1);
            triple{i}.pt(:,1)=index_add(j,2)+triple{i}.pt(:,1);
        end
        
        if strcmp(pair_4{i,3},pair_2{j,1})&&...
           strcmp(pair_4{i,4},pair_2{j,2})
            triple{i}.line(:,2)=index_add(j,1)+triple{i}.line(:,2);
            triple{i}.pt(:,2)=index_add(j,2)+triple{i}.pt(:,2);
        end
    end
    
end

%The sum of each column
colsum=zeros(lineIst,1);
for i=1:size(pair_2,1)
    %Update the index
    pairwise{i}.line(:,1:2)=pairwise{i}.line(:,1:2)+index_add(i,1);
    pairwise{i}.pt(:,1:2)=pairwise{i}.pt(:,1:2)+index_add(i,1);
    pairwise{i}.pt(:,3)=pairwise{i}.pt(:,3)+index_add(i,2);
    %
    for j=1:size(pairwise{i}.line,1)
        
        colid=pairwise{i}.line(j,2)+1;
        colsum(colid)=...
            colsum(colid)+...
            pairwise{i}.line(j,3);
        
        colid=pairwise{i}.line(j,1)+1;
        colsum(colid)=...
            colsum(colid)+...
            pairwise{i}.line(j,3);
    end
    
    for j=1:size(pairwise{i}.pt,1)
        
        colid=pairwise{i}.pt(j,1)+1;
        value=(pairwise{i}.pt(j,4)+pairwise{i}.pt(j,5))/2;
        colsum(colid)=...
            colsum(colid)+value;
        
        colid=pairwise{i}.pt(j,2)+1;
        colsum(colid)=...
            colsum(colid)+value;
        
        colid=pairwise{i}.pt(j,3)+1;
        colsum(colid)=...
            colsum(colid)+value;
        
       
        pairwise{i}.pt(j,4)=value;
    end
    
end

%With triplet candidates
for i=1:size(pair_4,1)
    for j=1:size(triple{i}.line,1)
        colid=triple{i}.line(j,2)+1;
        colsum(colid)=...
            colsum(colid)+...
            triple{i}.line(j,3);
        
        colid=triple{i}.line(j,1)+1;
        colsum(colid)=...
            colsum(colid)+...
            triple{i}.line(j,3);
    end
    
    for j=1:size(triple{i}.pt,1)
        
        value=triple{i}.pt(j,3);
        
        colid=triple{i}.pt(j,1)+1;
        colsum(colid)=...
            colsum(colid)+value;
        
        colid=triple{i}.pt(j,2)+1;
        colsum(colid)=...
            colsum(colid)+value;
      
    end
end
%Normalization with the maxmum
maxD = max(colsum);
for i=1:size(pair_2,1)
    pairwise{i}.line(:,3)=pairwise{i}.line(:,3)/maxD;
    pairwise{i}.pt(:,4)=pairwise{i}.pt(:,4)/maxD;
end

for i=1:size(pair_4,1)
    triple{i}.line(:,3)=triple{i}.line(:,3)/maxD;
    triple{i}.pt(:,3)=triple{i}.pt(:,3)/maxD;
end

%Parameters of random walks
c= 0.2;
amp_max= 30;
iterMax= 300;
thresConvergence=1e-15;
tolC= 1e-3 ;

%Initialization the random walk
nMatch = length(colsum);
prev_score = ones(nMatch,1)/nMatch; % buffer for the previous score
prev_score2 = prev_score;         % buffer for the two iteration ahead
prev_assign = prev_score ; % buffer for Sinkhorn result of prev_score
bCont = 1;  iter_i = 0;

subweights=zeros(size(pair_2,1),1);
%% Walks
while bCont && iter_i < iterMax
    iter_i = iter_i + 1;
    %% Reweight
    rwscore=( c*prev_score + (1-c)*prev_assign );
    cur_score =zeros(nMatch,1);
    

    for i=1:size(pair_2,1)
        for j=1:size(pairwise{i}.line,1)
            
            rowid=pairwise{i}.line(j,1)+1;
            colid=pairwise{i}.line(j,2)+1;
            
            cur_score(rowid)=cur_score(rowid)+...
                pairwise{i}.line(j,3)*rwscore(colid);
            
            cur_score(colid)=cur_score(colid)+...
                pairwise{i}.line(j,3)*rwscore(rowid);
            
        end
    end
    %points weight
    for i=1:size(pair_2,1)
        for j=1:size(pairwise{i}.pt,1)
         
            rowid=pairwise{i}.pt(j,1)+1;
            colid=pairwise{i}.pt(j,3)+1;
            
            cur_score(rowid)=cur_score(rowid)+...
                pairwise{i}.pt(j,4)*rwscore(colid);
            
            cur_score(colid)=cur_score(colid)+...
                pairwise{i}.pt(j,4)*rwscore(rowid);
            
            %l2 pt
            rowid=pairwise{i}.pt(j,2)+1;
            colid=pairwise{i}.pt(j,3)+1;
            
            cur_score(rowid)=cur_score(rowid)+...
                pairwise{i}.pt(j,4)*rwscore(colid);
            
            cur_score(colid)=cur_score(colid)+...
                pairwise{i}.pt(j,4)*rwscore(rowid);
            
        end
    end
    
    %Triplet line candidates
    for i=1:size(pair_4,1)
        for j=1:size(triple{i}.line,1)
            
            rowid=triple{i}.line(j,1)+1;
            colid=triple{i}.line(j,2)+1;
            
            cur_score(rowid)=cur_score(rowid)+...
                triple{i}.line(j,3)*rwscore(colid);
            
            cur_score(colid)=cur_score(colid)+...
                triple{i}.line(j,3)*rwscore(rowid);
            
        end
    end
    
    %Triplet point candidates
    for i=1:size(pair_4,1)
        for j=1:size(triple{i}.pt,1)
         
           
            rowid=triple{i}.pt(j,1)+1;
            colid=triple{i}.pt(j,2)+1;
            
            cur_score(rowid)=cur_score(rowid)+...
                triple{i}.pt(j,3)*rwscore(colid);
            
            cur_score(colid)=cur_score(colid)+...
                triple{i}.pt(j,3)*rwscore(rowid);
        end
    end
    
    sumCurScore = sum(cur_score); % normalization of sum 1
    if sumCurScore>0, cur_score = cur_score./sumCurScore; end
    
    cur_assign = cur_score;
    % attenuate small values and amplify large values
    amp_value = amp_max/ max(cur_assign);  % compute amplification factor
    cur_assign = exp( amp_value*cur_assign );  % amplification 
    
   
    for i=1:size(pair_2,1)
        subweights(i)=sum(cur_assign(candis{i}.stid:candis{i}.edid));
    end
    subweights=subweights./sum(subweights);
    
    %Sinkhorn method of iterative bistocastic normalizations seperately for
    %each stereo
    for i=1:size(pair_2,1)
        X_slack = [cur_assign(candis{i}.stid:candis{i}.edid);...
                   candis{i}.dumVal*ones(candis{i}.dumSize,1)];
        X_slack = mexBistocNormalize_match_slack(X_slack, int32(candis{i}.idx1),...
        int32(candis{i}.ID1), int32(candis{i}.idx2), int32(candis{i}.ID2),...
        tolC, candis{i}.dumDim, candis{i}.dumVal, int32(1000));
    
        cur_assign(candis{i}.stid:candis{i}.edid)=subweights(i)*...
            X_slack(1:(candis{i}.edid-candis{i}.stid)+1);
    end
    
     sumCurAssign = sum(cur_assign); % normalization of sum 1
    if sumCurAssign>0, cur_assign = cur_assign./sumCurAssign; end
    
    diff1 = sum((cur_score-prev_score).^2)
    diff2 = sum((cur_score-prev_score2).^2); % to prevent oscillations
    
    if diff1 < thresConvergence
        bCont = 0;
    end
   
    prev_score2 = prev_score;
    prev_score = cur_score;
    prev_assign = cur_assign;
    
end

X = cur_score;
 save([write_folder...
        'X.txt'],'X','-ascii');
















