/ani/mrses

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;