RSS

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

To get this branch, use:
bzr branch /webbzr/ani/mrses
Line Revision Contents
1 1
function res=mrses_software(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_cl(ctx, 1, block, A, B, k);
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=[]; ke=[]; 
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 = multi_bmc(dA, dB, mdiff, ki, distmod);
58
%    for i=1:block_size
59
%	check_dist(i) = bmc(A(:,ki(:,i)),B(:,ki(:,i)), distmod); 
60
%    end
61
%    find(cur_dist ~= check_dist)
62
63
    for iter=1:Niter
64
	xki=ceil(rand(1,block_size)*k);			% selecting random gen from selected
65
	xke=ceil(rand(1,block_size)*(genes-k));		% selected random gen from non-selected
66
67
	idx_i = sub2ind(size(ki), xki, 1:block_size);
68
	idx_e = sub2ind(size(ke), xke, 1:block_size);
69
	
70
	t=ki(idx_i);
71
	ki(idx_i)=ke(idx_e);
72
	ke(idx_e)=t;
73
	
74
	dist = multi_bmc(dA, dB, mdiff, ki, distmod);
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)/1000000;
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
% DISTANCE CALCULATOR
105
% ifs are taking to much time, do a separate functions for each distance
106
function dist=multi_bmc(dA, dB, mean_diff, ki, distmod)
107
    block_size = size(ki,2);
108
    
109
    %ki(:,1) = [5,7,9,11,13];
110
    for i=1:block_size
111
	x = dA(:, ki(:,i));
112
	c1 = x'*x;
113
	x = dB(:, ki(:,i));
114
	c2 = x'*x;
115
	
116
	c=(c1+c2)./2;
117
	
118
	[L,p] = chol(c);
119
	%detc = prod(diag(L,0))^2;
120
	%tmp = diag(L,0);
121
	%detc = tmp' * tmp;
122
	if p > 0
123
	    detc = 0;
124
	else
125
	    detc = det(L)^2;
126
	end
127
128
%	rcorr(i)=log((detc.*detc)./(det(c1).*det(c2)));
129
	rcorr(i)=2.*log(detc./sqrt(det(c1).*det(c2)));
130
%	rcorr(i)=2.*log(detc./sqrt(det(c1*c2)));
131
132
	if detc == 0
133
	    mdiff = mean_diff(ki(:,i));
134
      	    rmahal(i)=(mdiff*pinv(c))*mdiff';
135
	else
136
%      	    rmahal(i)=(mdiff/c)*mdiff';
137
%           rmahal(i)=((mdiff/L)/L')*mdiff';
138
%	    rmahal(i) = prod(mdiff/L)^2;
139
	    tmp = mean_diff(ki(:,i))/L;
140
	    rmahal(i) = tmp * tmp';
141
	end
142
    end
143
144
    if (distmod==1) 
145
	dist = rmahal./8 + rcorr./4;
146
    elseif (distmode==2) 
147
	dist = rmahal;
148
    else 
149
	dist = rcorr;
150
    end
151
152
function dist=bmc(x1,x2,distmod)
153
    c1=cov1(x1);
154
    c2=cov1(x2);
155
    c=(c1+c2)./2;
156
157
    if (distmod~=2)
158
	rcorr=2.*log(det(c)./sqrt(det(c1).*det(c2)));
159
    end
160
161
    if (distmod~=3)
162
	m1=mean(x1);
163
	m2=mean(x2);
164
	rmahal=((m2-m1)/c)*(m2-m1)';
165
    end
166
    
167
    if (distmod==1) 
168
	dist = rmahal./8+rcorr./4;
169
    elseif (distmode==2) 
170
	dist=rmahal;
171
    else 
172
	dist=rcorr;
173
    end
174
175
176
function c = cov1(x, m)
177
    [rows, cols] = size(x);
178
    %rows = size(x(:,1))
179
    
180
    if (nargin<2) 
181
	nX = x - ones([rows,1]) * mean(x);
182
    else
183
	nX = x - ones([rows,1]) * m';
184
    end
185
186
    c = nX' * nX / rows;
187
188
function c = cov2(x)
189
    c = x' * x;

Loggerhead 1.18.1 is a web-based interface for Bazaar branches