Результаты Affinity Propagation не совпадают

Я пытаюсь реализовать алгоритм кластеризации Affinity Propagation на C++. В рамках тестирования я хочу сравнить свои результаты с хорошо зарекомендовавшими себя реализациями алгоритма в Matlab (Ссылка) и в R (пакет apcluster). К сожалению, кластеризации не совпадают.

Чтобы быть более точным, (тестовый) набор данных:

0.9411760 0.9702140  
0.9607826 0.9744693  
0.9754896 0.9574479  
0.9852929 0.9489372  
0.9950962 0.9234050  
1.0000000 0.8936175  
1.0000000 0.8723408  
0.9852929 0.8595747  
1.0000000 0.8893622  
1.0000000 0.9191497

В R я набрал:

S<-negDistMat(data)  
A<-apcluster(S,maxits=1000,convits=100, lam=0.9,q=0.5) 

и получил:

> A@idx 
2 2 2 5 5 9 9 9 9 5  
2 2 2 5 5 9 9 9 9 5

В Matlab я просто набрал:

[idx,netsim,dpsim,expref]=apcluster(S,diag(S));

Из файла apcluster.m, реализующего apcluster (строка 77):

maxits=1000; convits=100; lam=0.9; plt=0; details=0; nonoise=0;

Это объясняет параметры для R, в Matlab они являются значениями по умолчанию. Поскольку мне удобнее использовать R в отношении Affinity Propagation, по причинам сравнения я придерживался значений по умолчанию Matlab, просто чтобы не испортить что-то непреднамеренно.

..но получил:

>> idx'  
ans =  
3     3     3     3     5     9     9     9     9     5

В обоих случаях матрицы сходства совпали. Что я мог упустить?

Обновление:
Я также реализовал код Matlab, предложенный Frey & Dueck в их оригинальной публикации. (Вы можете заметить, что я опустил шум), и хотя я могу воспроизвести индексы, предоставленные прежней реализацией Matlab, матрицы доступности и ответственности различаются по некоторым значениям. Ошибка меньше 0,01, но это существенно.

Их код:

function [idx,A,R]=frey(S);

N=size(S,1);
A=zeros(N,N);
R=zeros(N,N);
lam=0.9; % Set damping factor
for iter=1:122
    % Compute responsibilities
    Rold=R;
    AS=A+S;
    [Y,I]=max(AS,[],2);
    for i=1:N
            AS(i,I(i))=-realmax;
    end;
    [Y2,I2]=max(AS,[],2);
    R=S-repmat(Y,[1,N]);
    for i=1:N
            R(i,I(i))=S(i,I(i))-Y2(i);
    end;
    R=(1-lam)*R+lam*Rold; % Dampen responsibilities
    % Compute availabilities
    Aold=A;
    Rp=max(R,0);
    for k=1:N
            Rp(k,k)=R(k,k);
    end;
    A=repmat(sum(Rp,1),[N,1])-Rp;
    dA=diag(A);
    A=min(A,0);
    for k=1:N
            A(k,k)=dA(k);
    end;
    A=(1-lam)*A+lam*Aold; % Dampen availabilities
end;

E=R+A; % Pseudomarginals
I=find(diag(E)>0); K=length(I); % Indices of exemplars
[tmp c]=max(S(:,I),[],2); c(I)=1:K; idx=I(c); % Assignments