% C1 and C2 are K x M dimensional, with K the number of clusters and M the
% number of spatio-temporal features
function [RealCorrs,NullCors] = MOT_ANA_CompareClusters(C1,C2,dist)

    
    if size(C2,1) >= size(C1,1)
        
        switch dist
            case 'SpatialCorr'
        
                
                
                Dis = 1-corr(C1',C2');

                ExtraDis1 = 1-corr(C1',C1');
                ExtraDis2 = 1-corr(C2',C2');
                
                % First, we use Hungarian algorithm to match
                [assignment] = munkres(Dis);
            case 'MSE'
                
                for i = 1:size(C1,1)
                    for j = 1:size(C2,1)
                        Dis(i,j) = sqrt(sum((C1(i,:)-C2(j,:)).^2));
                    end
                end
                
                for i = 1:size(C1,1)
                    for j = 1:size(C1,1)
                        ExtraDis1(i,j) = sqrt(sum((C1(i,:)-C1(j,:)).^2));
                    end
                end

                for i = 1:size(C2,1)
                    for j = 1:size(C2,1)
                        ExtraDis2(i,j) = sqrt(sum((C2(i,:)-C2(j,:)).^2));
                    end
                end

                [assignment] = munkres(Dis);
                
            otherwise
                
        end
        
        % We reorder
        C2 = C2(assignment,:);

        
    else
        
        switch dist
            case 'SpatialCorr'
                Dis = 1-corr(C2',C1');

                ExtraDis1 = 1-corr(C1',C1');
                ExtraDis2 = 1-corr(C2',C2');
                
                % First, we use Hungarian algorithm to match
                [assignment] = munkres(Dis);
        
            case 'MSE'
                
                
                for i = 1:size(C2,1)
                    for j = 1:size(C1,1)
                        Dis(i,j) = sqrt(sum((C2(i,:)-C1(j,:)).^2));
                    end
                end
                
                
                for i = 1:size(C1,1)
                    for j = 1:size(C1,1)
                        ExtraDis1(i,j) = sqrt(sum((C1(i,:)-C1(j,:)).^2));
                    end
                end

                for i = 1:size(C2,1)
                    for j = 1:size(C2,1)
                        ExtraDis2(i,j) = sqrt(sum((C2(i,:)-C2(j,:)).^2));
                    end
                end
                
                [assignment] = munkres(Dis);
                
            otherwise
                
        end
                
        % We reorder
        C1 = C1(assignment,:);
    end
    
    switch dist
        case 'SpatialCorr'
    
            % We compute the diagonal and off diagonal elements
            NewDisMat = 1-corr(C1',C2');

            EXTRA1 = jUpperTriMatToVec(ExtraDis1);
            EXTRA2 = jUpperTriMatToVec(ExtraDis2);
            
            
        case 'MSE'
            
            for i = 1:size(C1,1)
                for j = 1:size(C2,1)
                    NewDisMat(i,j) = sqrt(sum((C1(i,:)-C2(j,:)).^2));
                end
            end
            
            
            
            EXTRA1 = jUpperTriMatToVec(ExtraDis1);
            EXTRA2 = jUpperTriMatToVec(ExtraDis2);
    end
            
    RealCorrs = diag(NewDisMat);
    
    NullCors = jUpperTriMatToVec(NewDisMat);
    NullCors = [NullCors;jUpperTriMatToVec(NewDisMat')];
    NullCors = [NullCors;EXTRA1;EXTRA2];
    
    figure;
    set(gcf,'color','w');
    
    COL{1} = [0.8,0.8,0.8];
    COL{2} = [0 0 0];
    MakeViolin(NullCors',gca,{''},'',COL,1,1);
    
    hold on;
    plot(ones(1,size(C1,1)),RealCorrs,'rd');

    MAXVAL = max([NullCors;RealCorrs])+0.05*max([NullCors;RealCorrs]);
    MINVAL = min([NullCors;RealCorrs])-0.05*min([NullCors;RealCorrs]);
    
    ylim([MINVAL,MAXVAL]);
    
end