function res=mrses_hw_distance(A,B,k,Niter,Ncycle,distmod,block) if (nargin<7) block=256; end 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 nA=sa(1); nB=sb(1); %optki=zeros(Ncycle,k); ctx = mrses_hw(); mrses_hw(ctx, 1, k, block, single(A), single(B), distmod); for icycle=0:block:(Ncycle-1) block_size = min(block, Ncycle - icycle); % SELECT k GENES {ki} FOR TEST AND EXCLUDE THEM FROM ALL GENES {ke} ki=int16([]); ke=int16([]); for i=1:block_size tt=randperm(genes);% randomizing genes ki(:,i)=tt(1:k); % selecting first k ke(:,i)=tt(k+1:end); % the rest unuzed end cur_dist = mrses_hw(ctx, 10, block_size, ki); %{ check_dist = []; for i=1:block_size check_dist(i) = bmc(A(:,ki(:,i)),B(:,ki(:,i)), distmod); end dist_diff = abs(cur_dist - check_dist); allowed = abs(check_dist)/10000; find ( dist_diff > allowed) %} for iter=1:Niter xki=ceil(rand(1,block_size)*k); % selecting random gen from selected xke=ceil(rand(1,block_size)*(genes-k)); % selected random gen from non-selected idx_i = sub2ind(size(ki), xki, 1:block_size); idx_e = sub2ind(size(ke), xke, 1:block_size); t=ki(idx_i); ki(idx_i)=ke(idx_e); ke(idx_e)=t; dist = mrses_hw(ctx, 10, block_size, ki); %{ check_dist = []; for i=1:block_size check_dist(i)=bmc(A(:,ki(:,i)),B(:,ki(:,i)),distmod); % compute distance between A and B with currently selected genes end find(dist ~= check_dist); dist_diff = abs(dist - check_dist); allowed = abs(dist)/10000; find ( dist_diff > allowed) %} bad = find(dist < cur_dist); idx_i = idx_i(bad); idx_e = idx_e(bad); t=ki(idx_i); ki(idx_i)=ke(idx_e); ke(idx_e)=t; cur_dist = max(dist, cur_dist); end optki(:,(icycle+1):(icycle+block_size))=ki; % save finally selected genes end mrses_hw(ctx); optki=reshape(optki,1,[]); [n,g]=hist(optki,1:genes); H=[n./Ncycle;g]; res=flipud(sortrows(H')); function dist=bmc(x1,x2,distmod) c1=cov1(x1); c2=cov1(x2); c=(c1+c2)./2; if (distmod~=2) rcorr=2.*log(det(c)./sqrt(det(c1).*det(c2))); end if (distmod~=3) m1=mean(x1); m2=mean(x2); rmahal=((m2-m1)/c)*(m2-m1)'; end if (distmod==1) dist = rmahal./8+rcorr./4; elseif (distmode==2) dist=rmahal; else dist=rcorr; end function c = cov1(x, m) [rows, cols] = size(x); %rows = size(x(:,1)) if (nargin<2) nX = x - ones([rows,1]) * mean(x); else nX = x - ones([rows,1]) * m'; end c = nX' * nX / rows;