forked from andrewssobral/nway
-
Notifications
You must be signed in to change notification settings - Fork 0
/
tuckdemo.m
122 lines (100 loc) · 3.52 KB
/
tuckdemo.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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
%tuckdemo.m
%
% Made for the Self-Study application
%
% REQUIRES THE DATA SETS !!!!!
%
% Copyright (C) 1995-2006 Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, [email protected]
%
close all
clear all
load('data\dataset1');
load('data\dataset1res');
fprintf('\n 1 Pak ### Inspect raw data - press a key..\n');pause
figure(1);set(gcf,'Position',[-1 31 804 534]);
for i=1:28,
m=reshape(X(i,:),R(2),R(3));
mesh(EmAx,ExAx,m);
title(['Raw data. Time ' int2str(i)]);
xlabel('Excitation [nm]')
ylabel('Emission [nm]')
axis([EmAx(1) EmAx(311) ExAx(1) ExAx(20) 0 1000]);
grid on
drawnow
end;
fprintf('\n 2 Pak ### Inspect calibration data - press a key..\n');pause
%Remove wrong/disturbing observations
for i=1:28,
m=reshape(X(i,:),R(2),R(3));
mn=SetNaNs1(m,ExAx(1),ExAx(20),EmAx(1),EmAx(311),18,20,NaN,2,0);
mnr=cmatrep(mn,'gt',999.99,NaN);
mesh(EmAx,ExAx,mnr);
title(['Corrected data. Time ' int2str(i)]);
xlabel('Excitation [nm]')
ylabel('Emission [nm]')
axis([EmAx(1) EmAx(311) ExAx(1) ExAx(20) 0 1000]);
grid on
drawnow
end;
fprintf('\n 3 Pak ### Calculate all possible 1-4 models and list - press a key\n');pause
for i=1:37,
fprintf(' %2i [%i %i %i] %f \n',i,list(i,:));
end;
fprintf('\n 4 Pak ### Calculate all possible 1-4 models and list sorted and plot - press a key\n');pause
[i j]=sort(list(:,4));
lists=list(j,:);
for i=1:37,
fprintf(' %2i [%i %i %i] %f \n',i,lists(i,:));
end;
figure(1);set(gcf,'Position',[-1 31 804 534]);
plot(lists(:,4));
grid on
ylabel('Explained variation');
xlabel('Model ranking');
fprintf('\n 5 Pak ### Calculate a [1 3 3] model, list core and facplot - press a key\n');pause
format short
format compact
W=[1 3 3];
load('data\dataset3res') %[Factors,G,XExpl,Xf]=tucker(X,R,W);
int2str(G)
figure(1);set(gcf,'Position',[-1 31 804 534]);
% Convert factors to new format
clear ff,id1 = 0;
for i = 1:length(DimX)
id2 = sum(DimX(1:i).*W(1:i));
ff{i} = reshape(Factors(id1+1:id2),DimX(i),W(i));
id1 = id2;
end
Factors = ff;
plotfac(Factors,[],[1 28;250 440;250 560]);drawnow;legend('1','2','3')
set(gcf,'Position',[-1 31 804 534]);
fprintf('\n 6 Pak ### Continue with a [3 3 3] model, list core and facplot - press a key\n');pause
W=[3 3 3];
disp('[3 3 3] Before rotation')
int2str(G)
figure(1);set(gcf,'Position',[-1 31 804 534]);
plotfac(Factors,[],[1 28;250 440;250 560]);legend('1','2','3')
[A B C]=fac2let(Factors);
fprintf('\n 7 Pak ### Continue with a [3 3 3] model with VOS rotation and facplot - press a key\n');pause
disp('[3 3 3] After rotation to optimal variance of squares')
G = reshape(G,W);
[Gv,Ov1,Ov2,Ov3]=maxvar3(G,1e-5);
Gv = reshape(Gv,size(Gv,1),size(Gv,2)*size(Gv,3));
int2str(Gv)
Av=A*Ov1;Bv=B*Ov2;Cv=C*Ov3;SSres=sum(sum( (X-Av*Gv*kron(Cv',Bv')).^2 ));
figure(2);set(gcf,'Position',[-1 31 804 534]);
plotfac({Av;Bv;Cv},[],[1 28;250 440;250 560]);legend('1','2','3')
fprintf('\n 8 Pak ### The [3 3 3] model with DIA rotation and facplot - press a key\n');pause
disp('[3 3 3] After rotation to optimal diagonality')
[Gd,Od1,Od2,Od3]=maxdia3(G,1e-5);
Gd = reshape(Gd,size(Gd,1),size(Gd,2)*size(Gd,3));
int2str(Gd)
Ad=A*Od1;Bd=B*Od2;Cd=C*Od3;SSres=sum(sum( (X-Ad*Gd*kron(Cd',Bd')).^2 ));
figure(3);set(gcf,'Position',[-1 31 804 534]);
plotfac({Ad; Bd;Cd},[],[1 28;250 440;250 560]);legend('1','2','3')
fprintf('\n 9 Pak ### Compare cores - press a key\n');pause
fprintf('\nNot rotated - Fig 1\n');disp(int2str(G));
fprintf('\nVOS rotated - Fig 2\n');disp(int2str(Gv));
fprintf('\nDIA rotated - Fig 3\n');disp(int2str(Gd));
format