-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnnmfEMG1122.m
107 lines (81 loc) · 2.01 KB
/
nnmfEMG1122.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
function [W_best, S_best,M_rec] = nnmfEMG1122(M,syn,option)
%% Normalize to Unit Variance
M(option.mus,:)=0;
if option.NormalizeToUnitVariance
std_val = std(M,0,2)+eps;
else
std_val = ones(size(M,1),1);
end
M = M ./ std_val;
%%
clear vaf rsq VAF RSQ r2_bootstat_lb r2_bootstat_ub vaf_bootstat_lb vaf_bootstat_ub
W_all = {};
S_all = {};
name={};
u2=1;
while u2 <= option.rep
opt = statset('MaxIter',1000,'Display','off');
[W0,S0] = nnmf(M,syn,'replicates',10,'options',opt,'algorithm','mult');
opt = statset('MaxIter',1000,'Display','off');
[W,S] = nnmf(M,syn,'w0',W0,'h0',S0,'options',opt,'replicates',20,'algorithm','als');
rec = (W*S);
if min(any(W,1)) == 0
continue;
end
% Saving the results ...
c = max(W);
W = W ./ c;
S = S .* repmat(c',1,size(S,2));
% if u2==1
% continue;
% end
W_all{u2} = W;
S_all{u2} = S;
name{u2}=num2str(u2);
u2=u2+1;
end
%% cluster W
[id,Wtype] = kmeanW(W_all, syn,0.6);
% calculate mean of each type
for i=1:syn
meanOfType{i} = mean(Wtype{i},2) ;
end
W_avg = cell2mat(meanOfType);
meanRepCorr=[];
% choose the best corr with W_avg
for i=1:option.rep
if min(id{i})==0
continue;
end
for k=1:syn
out{i}(k) = corr(W_all{i}(:,k), meanOfType{ id{i}(k) });
end
meanRepCorr(i) = mean(out{i});
end
if isempty(meanRepCorr)
W_best = W_all{1};
S_best = S_all{1};
else
[~,bestRepID] = max(meanRepCorr);
W_best = W_all{bestRepID};
S_best = S_all{bestRepID};
end
M_rec = (W_best .* std_val) * S_best;
%% plots
% for i=1:syn
% figure;
% subplot(2,1,1);hold on;
% bar(Wtype{1,i})
% subplot(2,1,2);hold on;
% bar([meanOfType{i}, W_best(:,i)])
%
%
%
% end
%%
% for i=1:5
% subplot(5,2,2*i-1)
% bar(W_all{33}(:,i))
% subplot(5,2,2*i)
% bar(W_best(:,i))
% end