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;
