/ani/mrses

To get this branch, use:
bzr branch http://suren.me/webbzr/ani/mrses

« back to all changes in this revision

Viewing changes to mrses_software.m

  • Committer: Suren A. Chilingaryan
  • Date: 2010-04-28 04:30:08 UTC
  • Revision ID: csa@dside.dyndns.org-20100428043008-vd9z0nso9axezvlp
Initial import

Show diffs side-by-side

added added

removed removed

Lines of Context:
mrses_software.m
 
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;