/ani/mrses

function res=mrses_hw_debug(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);

 

  mean_a = mean(A,1);

  mean_b = mean(B,1);

  mdiff = mean_a - mean_b;

 

  dA = (A - ones([nA,1]) * mean_a) / sqrt(nA);

  dB = (B - ones([nB,1]) * mean_b) / sqrt(nB);

 

  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_hw = mrses_hw(ctx, 10, block_size, ki);

    cur_dist = multi_bmc(dA, dB, mdiff, ki, distmod);

%    find(cur_dist ~= cur_dist_hw)

    dist_diff = abs(cur_dist_hw - cur_dist);

    allowed = abs(cur_dist)/100000;

    find ( dist_diff > allowed)

 

%    for i=1:block_size

%       check_dist(i) = bmc(A(:,ki(:,i)),B(:,ki(:,i)), distmod); 

%    end

%    find(cur_dist ~= check_dist)

 

    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 = multi_bmc(dA, dB, mdiff, ki, distmod);

%        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)/1000000;

%       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'));

 

 

% DISTANCE CALCULATOR

% ifs are taking to much time, do a separate functions for each distance

function dist=multi_bmc(dA, dB, mean_diff, ki, distmod)

    block_size = size(ki,2);

    

    %ki(:,1) = [5,7,9,11,13];

    for i=1:block_size

        x = dA(:, ki(:,i));

        %x(1:5,:)'

        c1 = x'*x;

        x = dB(:, ki(:,i));

        c2 = x'*x;

        

        c=(c1+c2)./2;

        

        [L,p] = chol(c);

        %detc = prod(diag(L,0))^2;

        %tmp = diag(L,0);

        %detc = tmp' * tmp;

        if p > 0

            detc = 0;

        else

            detc = det(L)^2;

        end

 

%       rcorr(i)=log((detc.*detc)./(det(c1).*det(c2)));

        rcorr(i)=2.*log(detc./sqrt(det(c1).*det(c2)));

%       rcorr(i)=2.*log(detc./sqrt(det(c1*c2)));

 

        if detc == 0

            mdiff = mean_diff(ki(:,i));

            rmahal(i)=(mdiff*pinv(c))*mdiff';

        else

%           rmahal(i)=(mdiff/c)*mdiff';

%           rmahal(i)=((mdiff/L)/L')*mdiff';

%           rmahal(i) = prod(mdiff/L)^2;

            tmp = mean_diff(ki(:,i))/L;

            rmahal(i) = tmp * tmp';

        end

    end

 

 

    if (distmod==1) 

        dist = rmahal./8 + rcorr./4;

    elseif (distmode==2) 

        dist = rmahal;

    else 

        dist = rcorr;

    end

 

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;

 

function c = cov2(x)

    c = x' * x;