RSS

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

To get this branch, use:
bzr branch /webbzr/ani/mrses
Line Revision Contents
1 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;

Loggerhead 1.18.1 is a web-based interface for Bazaar branches