% gdama13_01
% mixing of endmembers, visualized as vectors in 3D space
clear all;
fprintf('gdama13_01\n');
% box in 3D (x,y,z) space
xmin = 0;
xmax = 1;
ymin = 0;
ymax = 1;
zmin = 0;
zmax = 1;
% set up 3D plot
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-1, 1, -1, 1, -1, 1]);
% factors
% factor 1
f1 = [1, 0.2, 0.8]';
norm = 1.2*sqrt(f1'*f1);
f1=f1/norm;
% factor 2
f2 = [0.2, 1, 0.8]';
norm = 0.8*sqrt(f2'*f2);
f2=f2/norm;
gda_arrow3(f1,'r-',3);
gda_arrow3(f2,'r-',3);
% samples
s1=0.8*f1+(1-0.8)*f2;
s2=0.6*f1+(1-0.6)*f2;
s3=0.4*f1+(1-0.4)*f2;
s4=0.2*f1+(1-0.2)*f2;
gda_arrow3(s1,'k-',2);
gda_arrow3(s2,'k-',2);
gda_arrow3(s3,'k-',2);
gda_arrow3(s4,'k-',2);
% improvise outline of 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('E_1');
ylabel('E_2');
zlabel('E_3');
% view angle
view([-7,56]);
fprintf("Caption: Samples (black arrows) are linear combinations of end-members (red arrows)\n");
% gdama13_02
% factors visualized as vectors in 3D space
clear all;
fprintf('gdama13_02\n');
% define 3D region of (x,y,z) space
xmin = 0;
xmax = 1;
Lx = xmax-xmin;
ymin = 0;
ymax = 1;
Ly = ymax-ymin;
zmin = 0;
zmax = 1;
Lz = zmax-zmin;
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-1, 1, -1, 1, -1, 1]);
% factors
% factor 1
f1 = [1, 0.2, 0.8]';
norm = 1.2*sqrt(f1'*f1);
f1=f1/norm;
% factor 2
f2 = [0.2, 1, 0.8]';
norm = 0.8*sqrt(f2'*f2);
f2=f2/norm;
gda_arrow3(f1,'r-',3);
gda_arrow3(f2,'r-',3);
% samples
s1=0.8*f1+(1-0.8)*f2;
s2=0.6*f1+(1-0.6)*f2;
s3=0.4*f1+(1-0.4)*f2;
s4=0.2*f1+(1-0.2)*f2;
gda_arrow3(s1,'k-',2);
gda_arrow3(s2,'k-',2);
gda_arrow3(s3,'k-',2);
gda_arrow3(s4,'k-',2);
if( 1 )
% use SVD. The minus signs are just to get the vectors
% in quadrants that look good in the plot
S=zeros(4,3);
S(1,:)=s1';
S(2,:)=s2';
S(3,:)=s3';
S(4,:)=s4';
[U,LAMBDA,V] = svd(S);
v1 = -V(:,1);
v2 = -V(:,2);
v3 = -V(:,3);
else
% just use the mean vector, and two vectors perpendicular
% to it, one of which is in the F1-F2 plane
v1 = f1+f2;
norm = sqrt(v1'*v1);
v1=v1/norm;
v3 = cross(f1,f2);
norm = sqrt(v3'*v3);
v3=v3/norm;
v2 = cross(v1,v3);
norm = sqrt(v2'*v2);
v2=v2/norm;
end
gda_arrow3(v1,'b-',3);
gda_arrow3(v2,'b-',3);
gda_arrow3(v3,'b:',3);
% improvise 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('E_1');
ylabel('E_2');
zlabel('E_3');
% set view
view([-7,56]);
fprintf("Caption: Samples (black arrows) are linear combinations of factors (red arrows)\n");
fprintf('Two factors (solid blue) lie in the plane of the samples, one (blue dashed) is normal to it.\n')
% gdama13_03
% visualizatiov of two different sets of factors
% was vectors in 3D space
clear all;
fprintf('gdama13_03\n');
% 3D box in (x,y,z) space
xmin = 0;
xmax = 1;
Lx = xmax-xmin;
ymin = 0;
ymax = 1;
Ly = ymax-ymin;
zmin = 0;
zmax = 1;
Lz = zmax-zmin;
% plot
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-1, 1, -1, 1, -1, 1]);
% factors for Case 1
% factor 1
f1 = [1, 0.2, 0.8]';
norm = 1.2*sqrt(f1'*f1);
f1=f1/norm;
% factor 2
f2 = [0.2, 1, 0.8]';
norm = 0.8*sqrt(f2'*f2);
f2=f2/norm;
gda_arrow3(f1,'r-',3);
gda_arrow3(f2,'r-',3);
% samples
s1=0.8*f1+(1-0.8)*f2;
s2=0.6*f1+(1-0.6)*f2;
s3=0.4*f1+(1-0.4)*f2;
s4=0.2*f1+(1-0.2)*f2;
gda_arrow3(s1,'k-',2);
gda_arrow3(s2,'k-',2);
gda_arrow3(s3,'k-',2);
gda_arrow3(s4,'k-',2);
% improvise 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('E_1');
ylabel('E_2');
zlabel('E_3');
view([-7,56]);
fprintf("Caption: Samples (black arrows) and factors (red arrows)\n");
% plot
figure(2);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-1, 1, -1, 1, -1, 1]);
% factors for Case 2
f1a = [1, 0.2, 0.8]';
norm = 1.2*sqrt(f1a'*f1a);
f1a=f1a/norm;
f2a = [0.2, 1, 0.8]';
norm = 0.8*sqrt(f2a'*f2a);
f2a=f2a/norm;
f1 = f1a+0.2*f2a;
norm = 1.2*sqrt(f1'*f1);
f1=f1/norm;
f2 = f2a + 0.2*f1a;
norm = 0.8*sqrt(f2'*f2);
f2=f2/norm;
gda_arrow3(f1,'r-',3);
gda_arrow3(f2,'r-',3);
gda_arrow3(s1,'k-',2);
gda_arrow3(s2,'k-',2);
gda_arrow3(s3,'k-',2);
gda_arrow3(s4,'k-',2);
% improvise 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('E_1');
ylabel('E_2');
zlabel('E_3');
view([-7,56]);
fprintf("Caption: Samples (black arrows) and an alternative set of factors (red arrows)\n");
% gdama13_04
% factor analysis on Atlantic Rocks dataset
% using singular value decomposition
clear all;
fprintf('gdama13_04\n');
% load data
D = load('../data/rocks.txt');
sio2 = D(:,1); % SiO2
tio2 = D(:,2); % TiO2
als03 = D(:,3); % Al203
feot = D(:,4); % FeO-total
mgo = D(:,5); % MgO
cao = D(:,6); % CaO
na20 = D(:,7); % Na2O
k20 = D(:,8); % K2O
Ns = size(D);
N = Ns(1);
M = Ns(2);
% compute factors and factor loadings using singular value decompostion
S=D;
[U, LAMBDA, V] = svd(S,0);
lambda = diag(LAMBDA);
Ns = length(lambda);
F = V';
C = U*LAMBDA;
% plot singular values
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot( [1:Ns], lambda, 'k-', 'LineWidth', 2 );
plot( [1:Ns], lambda, 'ko', 'LineWidth', 2 );
% title('singular values, s(i)');
xlabel('index, i');
ylabel('lambda(i)');
fprintf('Caption: Singular-values.\n');
% display first five factors
for j = (1:5)
f1=F(j,:);
fprintf('factor %d\n', j);
fprintf('SiO2 %f\n', f1(1));
fprintf('TiO2 %f\n', f1(2));
fprintf('Al203 %f\n', f1(3));
fprintf('FeO-total %f\n', f1(4));
fprintf('MgO %f\n', f1(5));
fprintf('CaO %f\n', f1(6));
fprintf('Na2O %f\n', f1(7));
fprintf(' \n');
end
% plot loadings of factors (f2,f3,f4) in 3D (c2,c3,c4) space
cmin=(-50);
cmax=50;
figure(2);
clf;
set(gca, 'LineWidth', 3);
set(gca, 'FontSize', 14);
hold on;
axis( [cmin, cmax, cmin, cmax, cmin, cmax] );
plot3( C(:,2), C(:,3), C(:,4), 'k.', 'LineWidth', 2 );
plot3( [cmin,cmin], [cmin,cmin], [cmin,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmin,cmin], [cmin,cmax], [cmin,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmin,cmax], [cmin,cmin], [cmin,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmax], [cmax,cmax], [cmax,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmax], [cmax,cmin], [cmax,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmin], [cmax,cmax], [cmax,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmin], [cmin,cmin], [cmax,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmax], [cmin,cmin], [cmax,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmin,cmin], [cmax,cmin], [cmax,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmin,cmin], [cmax,cmax], [cmax,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmax], [cmax,cmin], [cmin,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmin], [cmax,cmax], [cmin,cmin], 'k-', 'LineWidth', 2 );
xlabel('C(2)');
ylabel('C(3)');
zlabel('C(4)');
view(3);
fprintf('Csption: Scatter plot of the loadings of factors 2, 3, and 4.\n');
% gdama13_05
% factor analysis on Atlantic Rocks dataset
% using singular value decomposition
% and the varimax procedure
clear all;
fprintf('gdama13_05');
% load data
D = load('../data/rocks.txt');
sio2 = D(:,1); % SiO2
tio2 = D(:,2); % TiO2
als03 = D(:,3); % Al203
feot = D(:,4); % FeO-total
mgo = D(:,5); % MgO
cao = D(:,6); % CaO
na20 = D(:,7); % Na2O
k20 = D(:,8); % K2O
Ns = size(D);
N = Ns(1);
M = Ns(2);
% compute factors and factor loadings using singular value decompostion
[U, S, V] = svd(D,0);
% keep only first five singular values
P=5;
F = V(:,1:P)';
C = U(:,1:P)*S(1:P,1:P);
% initial rotated factor matrix, FP, abd rotation matrix, MR
MR=eye(P,P);
FP=F;
% spike these factors using the varimax procedure
k = [2, 3, 4, 5]';
Nk = length(k);
% varimax is an iterative procedure, but converges very rapidly,
% so do only 3 iterations. Within each iteration, one applies
% the procedure to every pair of factors.
for iter = (1:3) % iterations
for ii = (1:Nk) % pairs of factors
for jj = (ii+1:Nk)
% spike factors i and j
i=k(ii);
j=k(jj);
% copy factors from matrix to vectors
fA = FP(i,:)';
fB = FP(j,:)';
% standard varimax procedure to determine rotation angle q
u = fA.^2 - fB.^2;
v = 2* fA .* fB;
A = 2*M*u'*v;
B = sum(u)*sum(v);
top = A - B;
CC = M*(u'*u-v'*v);
DD = (sum(u)^2)-(sum(v)^2);
bot = CC - DD;
q = 0.25 * atan2(top,bot);
cq = cos(q);
sq = sin(q);
fAp = cq*fA + sq*fB;
fBp = -sq*fA + cq*fB;
% put back into factor matrix, FP
FP(i,:) = fAp';
FP(j,:) = fBp';
% accumulate rotation
mA = MR(i,:)';
mB = MR(j,:)';
mAp = cq*mA + sq*mB;
mBp = -sq*mA + cq*mB;
MR(i,:) = mAp';
MR(j,:) = mBp';
end
end
end
% new factor loadings
CP=C*MR';
% display first five factors
for j = (1:5)
f1=FP(j,:);
fprintf('factor %d\n', j);
fprintf('SiO2 %f\n', f1(1));
fprintf('TiO2 %f\n', f1(2));
fprintf('Al203 %f\n', f1(3));
fprintf('FeO-total %f\n', f1(4));
fprintf('MgO %f\n', f1(5));
fprintf('CaO %f\n', f1(6));
fprintf('Na2O %f\n', f1(7));
fprintf('K2O %f\n', f1(8));
fprintf(' \n');
end
gda_draw(' ', abs(F(2,:))', 'caption f2', abs(F(3,:))', 'caption f3', abs(F(4,:))', 'caption f4', abs(F(5,:))', 'caption f5', ...
' ', abs(FP(2,:))', 'caption f2p', abs(FP(3,:))', 'caption f3p', abs(FP(4,:))', 'caption f4p', abs(FP(5,:))', 'caption f5p');
fprintf('Caption: (Left) Factors 2-5 before varimax rotation. (Right) After rotation.\n');
% This is a test that the rotated matrices are correct
% (that is, that they reproduce the original sample matrix)
e = C*F-CP*FP;
E = sqrt(sum(sum(e.^2))) / sqrt(sum(sum((C*F).^2)));
fprintf('Fractional error in reconstruction: %.3f percent\n', 100*E );
% gdama13_06
% factor analysis with simple shapes
% that allows these shapes to be classified
clear all;
fprintf('gdama13_06\n');
% figure 1 is plot of samples
figure(1);
clf;
% define samples. These are simple shapes
% that might represent cross-sections through
% a mountain range
N=15;
M=11;
S=zeros(N,M);
S(1,:) = [0, 1, 2, 3, 4, 5, 4, 3, 2, 1, 0];
S(2,:) = [0, 2, 3, 4, 5, 6, 5, 4, 3, 1, 0];
S(3,:) = [0, 1, 2, 4, 5, 6, 5, 4, 2, 1, 0];
S(4,:) = [0, 1, 3, 5, 5, 4.5, 4, 3.5, 3, 2.5, 0];
S(5,:) = fliplr(S(4,:));
S(6,:) = [0, 6, 6, 6, 5, 5, 5, 3.5, 2, 1, 0];
S(7,:) = fliplr(S(6,:));
S(8,:) = [0, 0.25, 0.5, 1, 3, 6, 3, 1, 0.5, 0.25, 0];
S(9,:) = [0, 1, 2.5, 3, 4, 6, 4, 3, 2.5, 1, 0];
S(10,:) = [0, 1, 2, 4, 6, 4, 2, 1, 1, 0.5, 0];
S(11,:) = fliplr(S(10,:));
S(12,:) = [0, 0, 0.5, 0, 1, 1, 2, 4, 6, 3, 0];
S(13,:) = fliplr(S(12,:));
S(14,:) = [0, 0, 0.5, 1, 3, 6, 5.5, 5.5, 4.5, 3.5, 0];
S(15,:) = fliplr(S(14,:));
% plot samples
for i=(1:N)
subplot(5,3,i);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, M-1, 0, 7] );
plot( [0:M-1], S(i,:), 'k-', 'LineWidth', 2 );
plot( [0:M-1], S(i,:), 'k.', 'LineWidth', 2 );
xlabel('i');
ylabel(sprintf('s_{%d}',i));
end
fprintf('Caption: Mountain profiles 1 through 15.\n');
% factor analysis
[U, LAMBDA, V] = svd(S);
C = U*LAMBDA;
F = V';
fprintf('First 3 singular values: %f %f %f\n', LAMBDA(1,1), LAMBDA(2,2), LAMBDA(3,3) );
% figure 2 is plot of the first three factors
figure();
clf;
% plot first 3 factors. SIgn flip is just for aesthetic plot
for i=(1:3)
subplot(1,3,i);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 10, -1, 1] );
plot( [0:M-1], -F(i,:), 'k-', 'LineWidth', 2 );
plot( [0:M-1], -F(i,:), 'k.', 'LineWidth', 2 );
plot( [0:M-1], zeros(1,M), 'k:', 'LineWidth', 2 );
xlabel('i');
ylabel(sprintf('F_{%d}',i));
end
fprintf('Caption: First three EOFs (factors).\n');
% figure 3 is a plot of the profiles arranged by loadings 2 and 3
% figure 2 is plot of the first three factors
figure(3);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-8, 8, -8, 8] );
plot( [-8, 8], [0, 0], 'k:', 'LineWidth', 2);
plot( [0, 0], [-8, 8], 'k:', 'LineWidth', 2);
for i=(1:N)
% loadings
c2 = C(i,2);
c3 = C(i,3);
% demean the sample and rescale
scale1 = 0.4;
x = [0:M-1];
s = S(i,:);
s = s - mean(s);
s = s*scale1;
[smax, ismax] = max(s);
xmax = x(ismax);
xbar = mean(x);
x = x-xbar;
xmax = xmax-xbar;
x = x*scale1;
xmax = xmax*scale1;
% offset by scaled version of loadings
x = x+c2;
s = s+c3;
fill( x, s, [0, 0.5, 1.0], 'FaceAlpha', 0.5 );
plot( x, s, 'k-', 'Linewidth', 2);
plot( x, s, 'k.', 'Linewidth', 2);
plot( [x(1), x(M)], [s(1), s(M)], 'k-', 'LineWidth', 2 );
text( c2+xmax, c3, sprintf('%d',i), 'HorizontalAlignment', 'center', 'FontSize', 20 );
end
xlabel('Factor loading, C_{i2}','FontName','Cambria Math');
ylabel('Factor loading, C_{i3}','FontName','Cambria Math');
fprintf('Caption: zMountsin profiles (blue) arranged according to loadings of\n');
fprintf('EOFs (factors) 2 and 3');
% gdama13_07
% Synthetic EOF example using 2D images
% that are a mis of spatial patterns with
% time varying amplitudes
clear all;
fprintf('gdama13_07\n');
% time
Nt=25;
dt=1.0;
t = dt*[0:Nt-1]';
% space Nx by Nx grid
Nx=20;
Nx2=Nx*Nx;
dx=1.0;
L = dx*(Nx-1);
x = dx*[0:Nx-1]';
% cook up data, Nt images, each Nx by Nx
% Each image contains three spatial modes
% (patterns) of variation
% define the modes
m0 = zeros(Nx,Nx);
m1 = zeros(Nx,Nx);
m2 = zeros(Nx,Nx);
mn = 0.1 * random('norm',0,1,Nx,Nx);
for p = (1:Nx)
for q = (1:Nx)
m0(p,q) = sin(pi*x(p)/L).*sin(pi*x(q)/L);
m1(p,q) = sin(pi*x(p)/L).*sin(2*pi*x(q)/L);
m2(p,q) = sin(2*pi*x(p)/L).*sin(3*pi*x(q)/L);
end
end
% build data with modes with these amplitudes at different times
% These loadings were hand-tuned to look aesthetic
c0 = [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1 ]';
c1 = [1, 2, 3, 4, 5, 5, 4, 3, 2, 1, 1, 2, 3, 4, 5, 5, 4, 3, 2, 1, 1, 2, 3, 4, 5 ]';
c2 = [0, 0, 0, 1, 2, 3, 2, 1, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 1, 2, 3, 2, 1 ]';
% data is sum of modes with prescribed coefficients, plus random noise
figure();
set(gcf,'pos',[10 10 700 700]);
clf;
colormap('jet');
hold on;
D = zeros( Nx, Nx, Nt );
for p = (1:Nt)
D(:,:,p) = c0(p)*m0 + c1(p)*m1 + c2(p)*m2 + 0.7*random('norm',0,1,Nx,Nx);
subplot(5,5,p);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [0 20 0 20] );
axis off;
hold on;
imagesc(D(:,:,p));
text(0,-1.5,sprintf('%d',p));
hold on;
end
fprintf('Caption: Set of 25 sequential images (time shown below eacch image).')
% rearrange Nx by Nx image Nx*Nx as row vector
S=zeros(Nt,Nx2);
for r = (1:Nt)
for p = (1:Nx)
for q = (1:Nx)
s = Nx*(p-1)+q;
S(r,s) = D(p,q,r);
end
end
end
% test of code to rebuild Nx by Nx image froom Nx*Nx row vector
if( 0 )
D2 = zeros( Nx, Nx, Nt );
for r = (1:Nt)
for s = (1:Nx2)
p = floor( (s-1)/Nx+0.01 ) + 1;
q = s - Nx*(p-1);
D2(p,q,r) = S(r,s);
end
end
figure();
set(gcf,'pos',[10 10 700 700]);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
for p = (1:Nt)
subplot(5,5,p);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [0 20 0 20] );
axis off;
hold on;
imagesc(D2(:,:,p));
text(0,-1.5,sprintf('%d',p));
hold on;
end
end % end if test
% basic SVD analysis
[U,LAMBDA,V] = svd(S);
C = U*LAMBDA;
S2 = C*(V');
evalue = spdiags(LAMBDA,0);
% plot singular-values
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot([1:Nt], evalue, 'k-', 'LineWidth',3);
plot([1:Nt], evalue, 'ko', 'LineWidth',3);
fprintf('Caption: Singular-values.\n');
% test of code to rebuild D from S2=U*LAMBDA*V'=C*V'
if( 0 )
D3 = zeros( Nx, Nx, Nt );
for r = (1:Nt)
for s = (1:Nx2)
p = floor( (s-1)/Nx+0.01 ) + 1;
q = s - Nx*(p-1);
D3(p,q,r) = S2(r,s);
end
end
figure();
clf;
set(gcf,'pos',[10 10 700 700]);
for p = (1:Nt)
subplot(5,5,p);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis( [0 20 0 20] );
axis off;
imagesc(D3(:,:,p));
text(0,-1.5,sprintf('%d',p));
hold on;
end
end % end if test
% select largest eigenvalues, reconstruct data from them
Nf=3;
Ua=U(:,1:Nf);
LAMBDAa=LAMBDA(1:Nf,1:Nf);
Ca = Ua*LAMBDAa;
Va = V(:,1:Nf);
Sa = Ca * (Va');
% retrieve images of modes from V matrix
M = zeros( Nx, Nx, Nt );
for r = (1:Nt)
for s = (1:Nx2)
p = floor( (s-1)/Nx+0.01 ) + 1;
q = s - Nx*(p-1);
M(p,q,r) = V(s,r);
end
end
% plot reconstructed data
D4 = zeros( Nx, Nx, Nt );
for r = (1:Nt)
for s = (1:Nx2)
p = floor( (s-1)/Nx+0.01 ) + 1;
q = s - Nx*(p-1);
D4(p,q,r) = Sa(r,s);
end
end
figure(5);
set(gcf,'pos',[10 10 700 700]);
clf;
for p = (1:Nt)
subplot(5,5,p);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis( [0 20 0 20] );
axis off;
imagesc(D4(:,:,p));
text(0,-1.5,sprintf('%d',p));
end
fprintf('Caption: Reconstruction of set of 25 sequential images using three EOFs\n');
% plot first Nf modes
figure();
set(gcf,'pos',[10 10 700 700]);
clf;
for p = (1:Nf)
subplot(5,5,p);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis( [0 20 0 20] );
axis off;
imagesc(M(:,:,p));
text(0,-1.5,sprintf('%d',p));
end
fprintf('Caption: First 3 EOFs (factors).\n');
% plot timeseries of coefficients of modes
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
for p = (1:Nf)
subplot(Nf,1,p);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
plot( [1:Nt], Ca(:,p), 'k-', 'LineWidth', 3 );
plot( [1:Nt], Ca(:,p), 'ko', 'LineWidth', 2 );
xlabel('time');
ylabel(sprintf('C%d',p));
end
fprintf('Caption: Loadings C of first 3 EOFs (factors) as time series.\n');