RSS

(root)/ani/mrses : /mrses_hw_debug.m (revision 2)

To get this branch, use:
bzr branch /webbzr/ani/mrses
Line Revision Contents
1 1
function res=mrses_hw_debug(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
  mean_a = mean(A,1);
39
  mean_b = mean(B,1);
40
  mdiff = mean_a - mean_b;
41
42
  dA = (A - ones([nA,1]) * mean_a) / sqrt(nA);
43
  dB = (B - ones([nB,1]) * mean_b) / sqrt(nB);
44
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}
48
    
49
    ki=int16([]); ke=int16([]); 
50
    for i=1:block_size
51
	tt=randperm(genes);% randomizing genes
52
	
53
	ki(:,i)=tt(1:k);	% selecting first k
54
	ke(:,i)=tt(k+1:end);	% the rest unuzed
55
    end
56
    
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)
63
64
%    for i=1:block_size
65
%	check_dist(i) = bmc(A(:,ki(:,i)),B(:,ki(:,i)), distmod); 
66
%    end
67
%    find(cur_dist ~= check_dist)
68
69
    for iter=1:Niter
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
72
73
	idx_i = sub2ind(size(ki), xki, 1:block_size);
74
	idx_e = sub2ind(size(ke), xke, 1:block_size);
75
	
76
	t=ki(idx_i);
77
	ki(idx_i)=ke(idx_e);
78
	ke(idx_e)=t;
79
	
80
	dist = multi_bmc(dA, dB, mdiff, ki, distmod);
81
%        for i=1:block_size
82
%	    check_dist(i)=bmc(A(:,ki(:,i)),B(:,ki(:,i)),distmod); % compute distance between A and B with currently selected genes
83
%	end
84
	%find(dist ~= check_dist)
85
%	dist_diff = abs(dist - check_dist);
86
%	allowed = abs(dist)/1000000;
87
%	find ( dist_diff > allowed)
88
89
	
90
	bad = find(dist < cur_dist);
91
	idx_i = idx_i(bad);
92
	idx_e = idx_e(bad);
93
	
94
	t=ki(idx_i);
95
	ki(idx_i)=ke(idx_e);
96
	ke(idx_e)=t;
97
	
98
	cur_dist = max(dist, cur_dist);
99
    end
100
    optki(:,(icycle+1):(icycle+block_size))=ki;	% save finally selected genes
101
  end
102
  mrses_hw(ctx);
103
  
104
  optki=reshape(optki,1,[]);
105
  [n,g]=hist(optki,1:genes);
106
  H=[n./Ncycle;g];
107
  res=flipud(sortrows(H'));
108
109
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);
114
    
115
    %ki(:,1) = [5,7,9,11,13];
116
    for i=1:block_size
117
	x = dA(:, ki(:,i));
118
	%x(1:5,:)'
119
	c1 = x'*x;
120
	x = dB(:, ki(:,i));
121
	c2 = x'*x;
122
	
123
	c=(c1+c2)./2;
124
	
125
	[L,p] = chol(c);
126
	%detc = prod(diag(L,0))^2;
127
	%tmp = diag(L,0);
128
	%detc = tmp' * tmp;
129
	if p > 0
130
	    detc = 0;
131
	else
132
	    detc = det(L)^2;
133
	end
134
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)));
138
139
	if detc == 0
140
	    mdiff = mean_diff(ki(:,i));
141
      	    rmahal(i)=(mdiff*pinv(c))*mdiff';
142
	else
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';
148
	end
149
    end
150
151
152
    if (distmod==1) 
153
	dist = rmahal./8 + rcorr./4;
154
    elseif (distmode==2) 
155
	dist = rmahal;
156
    else 
157
	dist = rcorr;
158
    end
159
160
function dist=bmc(x1,x2,distmod)
161
    c1=cov1(x1);
162
    c2=cov1(x2);
163
    c=(c1+c2)./2;
164
165
    if (distmod~=2)
166
	rcorr=2.*log(det(c)./sqrt(det(c1).*det(c2)));
167
    end
168
169
    if (distmod~=3)
170
	m1=mean(x1);
171
	m2=mean(x2);
172
	rmahal=((m2-m1)/c)*(m2-m1)';
173
    end
174
    
175
    if (distmod==1) 
176
	dist = rmahal./8+rcorr./4;
177
    elseif (distmode==2) 
178
	dist=rmahal;
179
    else 
180
	dist=rcorr;
181
    end
182
183
184
function c = cov1(x, m)
185
    [rows, cols] = size(x);
186
    %rows = size(x(:,1))
187
    
188
    if (nargin<2) 
189
	nX = x - ones([rows,1]) * mean(x);
190
    else
191
	nX = x - ones([rows,1]) * m';
192
    end
193
194
    c = nX' * nX / rows;
195
196
function c = cov2(x)
197
    c = x' * x;

Loggerhead 1.18.1 is a web-based interface for Bazaar branches