function res=mrses_orig(A,B,k,Niter,Ncycle,distmod)

if (nargin<6)
    distmod=1;
end
if (nargin<5)
    Ncycle=1000;
end
if (nargin<4)
    Niter=500;
end
if (nargin<3)
    k=5;
end
if (nargin<2)
    error('As minimum two matrixes needed for MRSES');
end
if (nargin>6)
    error('Too much parameters');
end


sa=size(A);sb=size(B);
if (sa(2)==sb(2))
    genes=sa(2);
else
    error('Features dimension mismatch');
end


optki=zeros(Ncycle,k);
for icycle=1:Ncycle
    %   SELECT k GENES {ki} FOR TEST AND EXCLUDE THEM FROM ALL GENES {ke}
    tt=randperm(genes);

    ki=tt(1:k);
    ke=tt(k+1:end);

for iter=1:Niter
dist1=bmc(A(:,ki),B(:,ki));

xki=ceil(rand(1)*k);
xke=ceil(rand(1)*(genes-k));

t=ki(xki);
ki(xki)=ke(xke);
ke(xke)=t;

dist2=bmc(A(:,ki),B(:,ki));

if(dist2(distmod)<dist1(distmod))           % COMPARES BHATA DISTANCES 1-bhata, 2-mahal, 3-corr
    t=ki(xki);
    ki(xki)=ke(xke);
    ke(xke)=t;
end

end
optki(icycle,:)=ki;
end
optki=reshape(optki,1,[]);
[n,g]=hist(optki,1:genes);
H=[n./Ncycle;g];
res=flipud(sortrows(H'));


% DISTANCE CALCULATOR
function [rbhata,rmahal,rcorr]=bmc(x1,x2)

c1=cov(x1,1);
c2=cov(x2,1);
c=(c1+c2)./2;

m1=mean(x1);
m2=mean(x2);

rmahal=((m2-m1)/c)*(m2-m1)';
rcorr=2.*log(det(c)./sqrt(det(c1).*det(c2)));
rbhata=rmahal./8+rcorr./4;
