1
function res=mrses_hw_debug(A,B,k,Niter,Ncycle,distmod,block)
18
error('As minimum two matrixes needed for MRSES');
21
error('Too much parameters');
24
sa=size(A);sb=size(B);
28
error('Features dimension mismatch');
33
%optki=zeros(Ncycle,k);
36
mrses_hw(ctx, 1, k, block, single(A), single(B), distmod);
40
mdiff = mean_a - mean_b;
42
dA = (A - ones([nA,1]) * mean_a) / sqrt(nA);
43
dB = (B - ones([nB,1]) * mean_b) / sqrt(nB);
45
for icycle=0:block:(Ncycle-1)
46
block_size = min(block, Ncycle - icycle);
47
% SELECT k GENES {ki} FOR TEST AND EXCLUDE THEM FROM ALL GENES {ke}
49
ki=int16([]); ke=int16([]);
51
tt=randperm(genes);% randomizing genes
53
ki(:,i)=tt(1:k); % selecting first k
54
ke(:,i)=tt(k+1:end); % the rest unuzed
57
cur_dist_hw = mrses_hw(ctx, 10, block_size, ki);
58
cur_dist = multi_bmc(dA, dB, mdiff, ki, distmod);
59
% find(cur_dist ~= cur_dist_hw)
60
dist_diff = abs(cur_dist_hw - cur_dist);
61
allowed = abs(cur_dist)/100000;
62
find ( dist_diff > allowed)
65
% check_dist(i) = bmc(A(:,ki(:,i)),B(:,ki(:,i)), distmod);
67
% find(cur_dist ~= check_dist)
70
xki=ceil(rand(1,block_size)*k); % selecting random gen from selected
71
xke=ceil(rand(1,block_size)*(genes-k)); % selected random gen from non-selected
73
idx_i = sub2ind(size(ki), xki, 1:block_size);
74
idx_e = sub2ind(size(ke), xke, 1:block_size);
80
dist = multi_bmc(dA, dB, mdiff, ki, distmod);
82
% check_dist(i)=bmc(A(:,ki(:,i)),B(:,ki(:,i)),distmod); % compute distance between A and B with currently selected genes
84
%find(dist ~= check_dist)
85
% dist_diff = abs(dist - check_dist);
86
% allowed = abs(dist)/1000000;
87
% find ( dist_diff > allowed)
90
bad = find(dist < cur_dist);
98
cur_dist = max(dist, cur_dist);
100
optki(:,(icycle+1):(icycle+block_size))=ki; % save finally selected genes
104
optki=reshape(optki,1,[]);
105
[n,g]=hist(optki,1:genes);
107
res=flipud(sortrows(H'));
110
% DISTANCE CALCULATOR
111
% ifs are taking to much time, do a separate functions for each distance
112
function dist=multi_bmc(dA, dB, mean_diff, ki, distmod)
113
block_size = size(ki,2);
115
%ki(:,1) = [5,7,9,11,13];
126
%detc = prod(diag(L,0))^2;
135
% rcorr(i)=log((detc.*detc)./(det(c1).*det(c2)));
136
rcorr(i)=2.*log(detc./sqrt(det(c1).*det(c2)));
137
% rcorr(i)=2.*log(detc./sqrt(det(c1*c2)));
140
mdiff = mean_diff(ki(:,i));
141
rmahal(i)=(mdiff*pinv(c))*mdiff';
143
% rmahal(i)=(mdiff/c)*mdiff';
144
% rmahal(i)=((mdiff/L)/L')*mdiff';
145
% rmahal(i) = prod(mdiff/L)^2;
146
tmp = mean_diff(ki(:,i))/L;
147
rmahal(i) = tmp * tmp';
153
dist = rmahal./8 + rcorr./4;
160
function dist=bmc(x1,x2,distmod)
166
rcorr=2.*log(det(c)./sqrt(det(c1).*det(c2)));
172
rmahal=((m2-m1)/c)*(m2-m1)';
176
dist = rmahal./8+rcorr./4;
184
function c = cov1(x, m)
185
[rows, cols] = size(x);
189
nX = x - ones([rows,1]) * mean(x);
191
nX = x - ones([rows,1]) * m';