bzr branch
http://darksoft.org/webbzr/ani/mrses
1
by Suren A. Chilingaryan
Initial import |
1 |
function res=mrses_hw_distance(A,B,k,Niter,Ncycle,distmod,block) |
2 |
if (nargin<7) |
|
3 |
block=256; |
|
4 |
end
|
|
5 |
if (nargin<6) |
|
6 |
distmod=1; |
|
7 |
end
|
|
8 |
if (nargin<5) |
|
9 |
Ncycle=1000; |
|
10 |
end
|
|
11 |
if (nargin<4) |
|
12 |
Niter=500; |
|
13 |
end
|
|
14 |
if (nargin<3) |
|
15 |
k=5; |
|
16 |
end
|
|
17 |
if (nargin<2) |
|
18 |
error('As minimum two matrixes needed for MRSES'); |
|
19 |
end
|
|
20 |
if (nargin>6) |
|
21 |
error('Too much parameters'); |
|
22 |
end
|
|
23 |
||
24 |
sa=size(A);sb=size(B); |
|
25 |
if (sa(2)==sb(2)) |
|
26 |
genes=sa(2); |
|
27 |
else
|
|
28 |
error('Features dimension mismatch'); |
|
29 |
end
|
|
30 |
||
31 |
nA=sa(1); nB=sb(1); |
|
32 |
||
33 |
%optki=zeros(Ncycle,k);
|
|
34 |
||
35 |
ctx = mrses_hw(); |
|
36 |
mrses_hw(ctx, 1, k, block, single(A), single(B), distmod); |
|
37 |
||
38 |
for icycle=0:block:(Ncycle-1) |
|
39 |
block_size = min(block, Ncycle - icycle); |
|
40 |
% SELECT k GENES {ki} FOR TEST AND EXCLUDE THEM FROM ALL GENES {ke}
|
|
41 |
||
42 |
ki=int16([]); ke=int16([]); |
|
43 |
for i=1:block_size |
|
44 |
tt=randperm(genes);% randomizing genes |
|
45 |
||
46 |
ki(:,i)=tt(1:k); % selecting first k |
|
47 |
ke(:,i)=tt(k+1:end); % the rest unuzed |
|
48 |
end
|
|
49 |
||
50 |
cur_dist = mrses_hw(ctx, 10, block_size, ki); |
|
51 |
%{
|
|
52 |
check_dist = [];
|
|
53 |
for i=1:block_size
|
|
54 |
check_dist(i) = bmc(A(:,ki(:,i)),B(:,ki(:,i)), distmod);
|
|
55 |
end
|
|
56 |
dist_diff = abs(cur_dist - check_dist);
|
|
57 |
allowed = abs(check_dist)/10000;
|
|
58 |
find ( dist_diff > allowed)
|
|
59 |
%}
|
|
60 |
||
61 |
for iter=1:Niter |
|
62 |
xki=ceil(rand(1,block_size)*k); % selecting random gen from selected |
|
63 |
xke=ceil(rand(1,block_size)*(genes-k)); % selected random gen from non-selected |
|
64 |
||
65 |
idx_i = sub2ind(size(ki), xki, 1:block_size); |
|
66 |
idx_e = sub2ind(size(ke), xke, 1:block_size); |
|
67 |
||
68 |
t=ki(idx_i); |
|
69 |
ki(idx_i)=ke(idx_e); |
|
70 |
ke(idx_e)=t; |
|
71 |
||
72 |
dist = mrses_hw(ctx, 10, block_size, ki); |
|
73 |
%{
|
|
74 |
check_dist = [];
|
|
75 |
for i=1:block_size
|
|
76 |
check_dist(i)=bmc(A(:,ki(:,i)),B(:,ki(:,i)),distmod); % compute distance between A and B with currently selected genes
|
|
77 |
end
|
|
78 |
find(dist ~= check_dist);
|
|
79 |
dist_diff = abs(dist - check_dist);
|
|
80 |
allowed = abs(dist)/10000;
|
|
81 |
find ( dist_diff > allowed)
|
|
82 |
%}
|
|
83 |
||
84 |
bad = find(dist < cur_dist); |
|
85 |
idx_i = idx_i(bad); |
|
86 |
idx_e = idx_e(bad); |
|
87 |
||
88 |
t=ki(idx_i); |
|
89 |
ki(idx_i)=ke(idx_e); |
|
90 |
ke(idx_e)=t; |
|
91 |
||
92 |
cur_dist = max(dist, cur_dist); |
|
93 |
end
|
|
94 |
optki(:,(icycle+1):(icycle+block_size))=ki; % save finally selected genes |
|
95 |
end
|
|
96 |
mrses_hw(ctx); |
|
97 |
||
98 |
optki=reshape(optki,1,[]); |
|
99 |
[n,g]=hist(optki,1:genes); |
|
100 |
H=[n./Ncycle;g]; |
|
101 |
res=flipud(sortrows(H')); |
|
102 |
||
103 |
||
104 |
function dist=bmc(x1,x2,distmod) |
|
105 |
c1=cov1(x1); |
|
106 |
c2=cov1(x2); |
|
107 |
c=(c1+c2)./2; |
|
108 |
||
109 |
if (distmod~=2) |
|
110 |
rcorr=2.*log(det(c)./sqrt(det(c1).*det(c2))); |
|
111 |
end
|
|
112 |
||
113 |
if (distmod~=3) |
|
114 |
m1=mean(x1); |
|
115 |
m2=mean(x2); |
|
116 |
rmahal=((m2-m1)/c)*(m2-m1)'; |
|
117 |
end
|
|
118 |
||
119 |
if (distmod==1) |
|
120 |
dist = rmahal./8+rcorr./4; |
|
121 |
elseif (distmode==2) |
|
122 |
dist=rmahal; |
|
123 |
else
|
|
124 |
dist=rcorr; |
|
125 |
end
|
|
126 |
||
127 |
||
128 |
function c = cov1(x, m) |
|
129 |
[rows, cols] = size(x); |
|
130 |
%rows = size(x(:,1))
|
|
131 |
||
132 |
if (nargin<2) |
|
133 |
nX = x - ones([rows,1]) * mean(x); |
|
134 |
else
|
|
135 |
nX = x - ones([rows,1]) * m'; |
|
136 |
end
|
|
137 |
||
138 |
c = nX' * nX / rows; |