Commit 3758dadc authored by Felix Kong's avatar Felix Kong

add data and scripts (to be replaced)

scripts will be replaced with classes
parent 2453440b
% scales each column (the time-series for each state variable) to be
% between 0 and 1, and returns the scaling transform for each time series.
function [processedData, scale, offset] = ScaleData(X,method)
% scale each column
scale = nan(1,size(X,2));
offset = nan(1,size(X,2));
processedData = nan(size(X));
for k = 1:size(X,2)
if method == "range" % scale everything to be between 0 and 1
xmin = min(X(:,k));
xmax = max(X(:,k));
scale(k) = xmax - xmin;
if(scale(k) < 1e-3) % all numbers are the same
scale(k) = xmax;
offset(k) = 0;
processedData(:,k) = ones(size(X,1),1);
else
% normal case
offset(k) = xmin;
end
elseif method == "zscore" % scale to fit standard normal (centered at zero, std deviation of 1)
scale(k) = std(X(:,k));
offset(k) = mean(X(:,k));
else
error('second argument "method" must be either ''range'' or ''zscore''.')
end
processedData(:,k) = (X(:,k) - offset(k))/scale(k);
end
end
\ No newline at end of file
% convert hamlyn heart video to image sequence
% workingDir = tempname;
% mkdir(workingDir)
% mkdir(workingDir,'images')
%
% shuttleVideo = VideoReader('~/Downloads/hamlyn_heart.avi');
%
% ii = 1;
% while hasFrame(shuttleVideo)
% img = readFrame(shuttleVideo);
% filename = [sprintf('%03d',ii) '.jpg'];
% fullname = fullfile(workingDir,'images',filename);
% imwrite(img,fullname) % Write out to a JPEG file (img1.jpg, img2.jpg, etc.)
% ii = ii+1;
% end
%
% imageNames = dir(fullfile(workingDir,'images','*.jpg'));
% imageNames = {imageNames.name}';
%
%
% for ii = 1:length(imageNames)
% img = imread(fullfile(workingDir,'images',imageNames{ii}));
% grayImages(:,:,ii) = rgb2gray(img);
% end
%% Try DMD on it
clear variables
load data/HamlynHeartImages.mat;
disp('Data preparation...')
num_images = 100;
% for k = 1:size(grayImages,3)
for k = 1:num_images % first 100 images
image_k = grayImages(:,:,k);
Xall(:,k) = double(image_k(:));
end
X = Xall(:,1:end-1);
X2 = Xall(:,2:end);
disp('DMD Regression...')
[U0,S0,V0] = svd(X, 'econ');
figure(2);
clf
semilogy(diag(S0),'o');
%% Compute DMD (Phi are eigenvectors)
% r = size(X,2); % truncate at 5 modes
r = 50;
U = U0(:,1:r);
S = S0(1:r,1:r);
V = V0(:,1:r);
Atilde = U'*X2*(V/S);
[W,Lambda] = eig(Atilde);
Phi = X2*V*inv(S)*W;
eigs = diag(Lambda);
b0 = Phi\X(:,1);
%% Calculate DMD output
disp('Reconstructing...')
steps = size(X,2);
xkvec = zeros(steps,size(X,1));
% calculate the output states
for j = 1:steps
xkvec(j,:) = real(Phi*Lambda^(j-1)*b0)';
if(norm(imag(xkvec(j,:))) > 1e-3)
error('xkvec(j,:) isnt real (has non-negligble imaginary component)')
end
reconstruction(:,:,j) = uint8((reshape(xkvec(j,:),size(grayImages,1),size(grayImages,2)))); % im2uint8 doesnt work
end
%% replay
reconsructionHdl = implay(reconstruction);
set(reconsructionHdl.Parent,'Name',sprintf('Reconstruction with %d modes.',r));
originalHdl = implay(grayImages(:,:,1:num_images));
set(originalHdl.Parent,'Name','Original data');
%%
figure(3);
clf
mode_number = 1;
% dont use uint8 and imshow, U0 is unitary so columns have very small
% values
singleMode = surf(flipud(reshape(U0(:,mode_number),size(grayImages,1),size(grayImages,2))),'LineStyle','none');
colormap(flipud(gray))
modeEig = Lambda(mode_number,mode_number);
title(sprintf('Mode %d, eig %.3f + %.3fi', mode_number, real(modeEig), imag(modeEig)));
% view from above
axis([0 360 -0.001 288 -1 0])
view(0,90)
%% convert hamlyn heart video to image sequence
% workingDir = tempname;
% mkdir(workingDir)
% mkdir(workingDir,'images')
%
% shuttleVideo = VideoReader('~/Downloads/hamlyn_heart.avi');
%
% ii = 1;
% while hasFrame(shuttleVideo)
% img = readFrame(shuttleVideo);
% filename = [sprintf('%03d',ii) '.jpg'];
% fullname = fullfile(workingDir,'images',filename);
% imwrite(img,fullname) % Write out to a JPEG file (img1.jpg, img2.jpg, etc.)
% ii = ii+1;
% end
%
% imageNames = dir(fullfile(workingDir,'images','*.jpg'));
% imageNames = {imageNames.name}';
%
%
% for ii = 1:length(imageNames)
% img = imread(fullfile(workingDir,'images',imageNames{ii}));
% grayImages(:,:,ii) = rgb2gray(img);
% end
%% Try DMD on it
clear variables
load data/HamlynHeartImages.mat;
disp('Data preparation...')
image_start = 1;% frame number to start from
num_images = 100;
allImages = [];
% for k = 1:size(grayImages,3)
for k = image_start:image_start+num_images
image_k = grayImages(:,:,k);
% allImages(:,k) = double(image_k(:));
allImages = [allImages, double(image_k(:))];
end
XallPrescaled = [allImages(:,2:end);
allImages(:,1:end-1)];
% XallPrescaled = allImages;
% remove constant rows which make XallPrescaled lose rank
[XallPrescaledProcessed,ProcessSettings] = removeconstantrows(XallPrescaled);
XallPrescaledProcessed = XallPrescaledProcessed';
disp('DMD Regression...')
%% kernel dmd
disp('Scaling data...')
[Xall,scaleFactor,scaleOffset] = ScaleData(XallPrescaledProcessed,'zscore'); % zscore may cause problems with body_length, which is constant
% Xall = XallPrescaledProcessed;
% scaleFactor = 1;
% scaleOffset = 0;
X = Xall(1:end-1,:);
X2 = Xall(2:end,:);
%% Kernel DMD
disp('Doing Kernel DMD...')
% kernel function
% inputs: x1, x2 are vectors of the same size.
% kernel = @(x1,x2) (1 + dot(x1,x2))^10; % polynomial kernel, doesnt work that well. often fails during reconstruction
sigma2 = 1e6; % 36 works well for noise = 1e-6. but so does 169. Noise dependent?
kernel = @(x1,x2) exp(-norm(x1-x2)^2/sigma2); % Gausssian kernel
% for j = 1:size(X,1)
% disp(j)
% for k = 1:size(X2,1)
% Ahat(j,k) = kernel(X(k,:),X2(j,:)); % Y'Y2 is Ahat in Williams et al.
% Ghat(j,k) = kernel(X(k,:),X(j,:)); % Y*Y is the Gramian Ghat in Williams et al.
% end
% end
Ahat = GaussianKernel(X,X2,sqrt(sigma2));
Ghat = GaussianKernel(X,X,sqrt(sigma2));
%%
% TODO: truncate to get low-order model.
% use Ghat to compute Sigma (called 'S') and V
% [Qfull,S2full] = eig(Ghat); % since it's true that Ghat*Q = Q*Sigma^2
tic;
[Qfull, S2full] = eigs(Ghat,size(Ghat,1),'largestabs','Tolerance',1e-19); % use increased tolerance with GaussianKernel
toc;
% S2full = flipud(fliplr(S2full)); % make largest singular values first instead of largest singular values last
% Qfull = fliplr(Qfull);
Sfull = sqrtm(S2full);
% truncate
r = size(Ghat,1)-1; % number of dimensions to keep
% r = floor(size(Ghat,1)/2);
Q = Qfull(:,1:r);
S = Sfull(1:r,1:r);
figure(2);
clf;
semilogy(diag(Sfull),'o');
hold on
semilogy(diag(S),'.');
% Sinv = diag(ones(size(S,1),1)./diag(S)); % better numerical stability?
% Khat = (Sinv*Q')*Ahat*(Q*Sinv);
Khat = (S\Q')*(Ahat)*(Q/S); % Khat, an approximation of the Koopman operator
[Vhat, Lambda] = eig(Khat);
if(rank(Vhat) < size(Vhat,1))
error('Vhat is not full rank (rank = %d). I didnt code up the case where Vhat isnt.', rank(Vhat));
end
Phi = Q*S*Vhat; % eigenfunctions of Khat
Xi = Phi\X;
X2hat = Phi*Lambda*Xi;
if(norm(imag(X2hat)) > 1e-3)
error('X2hat has non-negligible imaginary component.')
end
X2hat = real(X2hat); % this is quite an accurate reconstruction of X2.
disp('Finished computing kernel-DMD regression. Attempting simulation from IC.');
%% simulate forward
num_modes = size(Q,1); % i THINK this is the correct one
tic;
xhat(1,:) = X(1,:);
% xhat(1,1) = xhat(1,1) + 0.1;
% xhat(1,:) = xhat(1,:) + 0.001*randn(1,size(xhat,2));
for t = 2:2*size(X,1)
% calculate phi given x
prev_x = xhat(t-1,:);
% kernelmat = nan(1,num_modes);
% for M = 1:num_modes
% kernelmat(M) = kernel(prev_x,X(M,:));
% end
kernelmat = GaussianKernel(prev_x,X, sqrt(sigma2));
kernelmat = kernelmat(1:num_modes)';
phi_at_t = kernelmat*Q/S*Vhat;
% estimate
xnext = phi_at_t*Lambda*Xi;
if(norm(imag(xnext)) > 1e-3)
error('xnext isnt real (has non-negligble imaginary component)')
end
xhat(t,:) = real(xnext);
end
toc;
%% Check
figure(1)
clf
subplot(2,1,1)
plot(Xall(:,1));
hold on
plot([nan(1);X2hat(:,1)],'--');
plot(xhat(:,1),'--');
legend('true','1-step recon', 'n-step recon');
title('Reconstruction element 1')
subplot(2,1,2)
plot(Xall(:,6));
hold on
plot([nan(1);X2hat(:,6)],'--');
plot(xhat(:,6),'--');
legend('true','1-step recon', 'n-step recon');
title('Reconstruction element 6');
disp('Simulated from IC. Unscaling and casting data.');
%%
% offsetx = repmat(scaleOffset,size(xhat,1),1);
% reconstruction = xhat*diag(scaleFactor) + offsetx;
for k = 1:size(xhat,1)
xhatRescaled(k,:) = xhat(k,:).*scaleFactor + scaleOffset;
end
xhatRescaled = xhatRescaled';
xhatRescaled = removeconstantrows('reverse',xhatRescaled,ProcessSettings); % add removed rows back in.
% reconstruction = reshape(xhatRescaled,size(grayImages,1),size(grayImages,2),num_images);
for k = 1:size(xhatRescaled,2)
reconstruction(:,:,k) = reshape(xhatRescaled(1:size(allImages,1),k),size(grayImages,1),size(grayImages,2));
end
disp('Done.')
%% replay
reconsructionHdl = implay(uint8(reconstruction));
set(reconsructionHdl.Parent,'Name',sprintf('Reconstruction: %d modes, sigma2 = %.2e.',r,sigma2));
originalHdl = implay(grayImages(:,:,image_start:image_start+num_images));
set(originalHdl.Parent,'Name','Original data');
%%
% figure(3);
% clf
% mode_number = 1;
% % dont use uint8 and imshow, U0 is unitary so columns have very small
% % values
% singleMode = surf(flipud(reshape(Qfull(:,mode_number),size(grayImages,1),size(grayImages,2))),'LineStyle','none');
% % singleMode = surf(flipud(reshape(U0(:,mode_number),size(grayImages,1),size(grayImages,2))),'LineStyle','none');
% colormap(flipud(gray))
% modeEig = Lambda(mode_number,mode_number);
% title(sprintf('Mode %d, eig %.3f + %.3fi', mode_number, real(modeEig), imag(modeEig)));
% % view from above
% axis([0 360 -0.001 288 -1 0])
% view(0,90)
function F = GaussianKernel(Xin,Yin,sigma)
X = Xin./sigma;
Y = Yin./sigma;
% each term is |xi - yj|^2. Expanding:
% |x-y|^2 = x'x - 2*x'y + y'y. these three lines are the three terms.
X2 = sum(X.^2,2); %x'x as a column vector
Y2 = sum(Y.^2,2); %y'y as a column vector
XY = bsxfun(@minus, X2, (2*X)*Y.'); % add 2x'y terms
total = bsxfun(@plus, Y2.', XY); % add y'y terms
F = exp(-total)';
end
clear all, close all, clc
load data/VortexData.mat % loads variables VORTALL, nx, ny
X = VORTALL(:,1:end-1);
X2 = VORTALL(:,2:end);
[U0,S0,V0] = svd(X,'econ');
%% Compute DMD (Phi are eigenvectors)
r = 21; % truncate at 21 modes
U = U0(:,1:r);
S = S0(1:r,1:r);
V = V0(:,1:r);
Atilde = U'*X2*V*inv(S);
[W,Lambda] = eig(Atilde);
Phi = X2*V*inv(S)*W;
eigs = diag(Lambda);
b0 = Phi\X(:,1);
%% prepare objects for animating the raw data
figure;
subplot(3,1,1);
xIC = reshape(VORTALL(:,1),nx,ny); % vortex field first timestep
im = imagesc(xIC); % plot vorticity field for first timestep
load CCcool.mat
colormap(CC); % use custom colormap
% clean up axes
set(gca,'XTick',[1 50 100 150 200 250 300 350 400 449],'XTickLabel',{'-1','0','1','2','3','4','5','6','7','8'})
set(gca,'YTick',[1 50 100 150 199],'YTickLabel',{'2','1','0','-1','-2'});
set(gcf,'Position',[100 100 300 130])
axis equal
hold on
% plot cylinder
theta = (1:100)/100'*2*pi;
x = 49+25*sin(theta);
y = 99+25*cos(theta);
fill(x,y,[.3 .3 .3]) % place cylinder
plot(x,y,'k','LineWidth',1.2) % cylinder boundary
title('Data')
colorbar;
% add contour lines (positive = solid, negative = dotted);
[~,cneg] = contour(xIC,[-5.5:.5:-.5 -.25 -.125],':k','LineWidth',1.2);
[~,cpos] = contour(xIC,[.125 .25 .5:.5:5.5],'-k','LineWidth',1.2);
%% prepare objects to animate the DMD reconstruction
subplot(3,1,2);
imDmd = imagesc(xIC);
load CCcool.mat
colormap(CC); % use custom colormap
% clean up axes
set(gca,'XTick',[1 50 100 150 200 250 300 350 400 449],'XTickLabel',{'-1','0','1','2','3','4','5','6','7','8'})
set(gca,'YTick',[1 50 100 150 199],'YTickLabel',{'2','1','0','-1','-2'});
set(gcf,'Position',[100 74 604 713])
axis equal
hold on
% plot cylinder
theta = (1:100)/100'*2*pi;
x = 49+25*sin(theta);
y = 99+25*cos(theta);
fill(x,y,[.3 .3 .3]) % place cylinder
plot(x,y,'k','LineWidth',1.2) % cylinder boundary
% add contour lines (positive = solid, negative = dotted);
[~,cnegDmd] = contour(xIC,[-5.5:.5:-.5 -.25 -.125],':k','LineWidth',1.2);
[~,cposDmd] = contour(xIC,[.125 .25 .5:.5:5.5],'-k','LineWidth',1.2);
titlename = sprintf('DMD reconstruction (%d dimensional, sorta)', r);
title(titlename)
colorbar;
%% Reconstruction error
subplot(3,1,3);
imErr = surf(zeros(size(xIC)));
set(gca,'XTick',[1 50 100 150 200 250 300 350 400 449],'XTickLabel',{'-1','0','1','2','3','4','5','6','7','8'})
set(gca,'YTick',[1 50 100 150 199],'YTickLabel',{'2','1','0','-1','-2'});
set(gcf,'Position',[100 74 604 713])
axis equal
hold on
imErr.EdgeAlpha = 0;
zlim([-0.01 0.01]);
% plot cylinder
theta = (1:100)/100'*2*pi;
x = 49+25*sin(theta);
y = 99+25*cos(theta);
fill(x,y,[.3 .3 .3]) % place cylinder
plot(x,y,'k','LineWidth',1.2) % cylinder boundary
title('Reconstruction error')
colorbar;
%% Run the animation
for tstep = 1:size(VORTALL,2) % every column is a timestep
% update ground truth
xTrue = reshape(VORTALL(:,tstep),nx,ny);
im.CData = xTrue;
cneg.ZData = xTrue;
cpos.ZData = xTrue;
% update dmd reconstruction
xkvec = Phi*Lambda^(tstep-1)*b0;
assert(all(abs(imag(xkvec)) <= 1e-6));
xk = reshape(real(xkvec),nx,ny);
imDmd.CData = xk;
cnegDmd.ZData = xk;
cposDmd.ZData = xk;
% update error plot
errork = xTrue-xk;
fprintf('Max abs error: %f \n', max(abs(errork),[],'all'));
imErr.ZData = errork;
view([0 90]);
drawnow;
% pause(0.1);
end
clear all, close all, clc
load data/VortexData.mat % loads variables VORTALL, nx, ny
X = VORTALL(:,1:end-1)';
X2 = VORTALL(:,2:end)';
Xall = VORTALL;
% [Xall,scaleFactor,scaleOffset] = ScaleData(XallPrescaled,'zscore'); % zscore may cause problems with body_length, which is constant
% % Xall2 = normalize(XallPrescaled');
% % Xall = normalize([dmd_data(i).X(:,1), dmd_data(i).X_dash],2,'range');
% X = Xall(1:end-1,:);
% X2 = Xall(2:end,:);
%% Kernel DMD
disp('Doing Kernel DMD...')
% kernel function
% inputs: x1, x2 are vectors of the same size.
% kernel = @(x1,x2) (1 + dot(x1,x2))^10; % polynomial kernel, doesnt work that well. often fails during reconstruction
sigma2 = 100000; % 36 works well for noise = 1e-6. but so does 169. Noise dependent?
kernel = @(x1,x2) exp(-norm(x1-x2)^2/sigma2); % Gausssian kernel
% for j = 1:size(X,1)
% for k = 1:size(X2,1)
% Ahat(j,k) = kernel(X(k,:),X2(j,:)); % Y'Y2 is Ahat in Williams et al.
% Ghat(j,k) = kernel(X(k,:),X(j,:)); % Y*Y is the Gramian Ghat in Williams et al.
% end
% end
Ahat = GaussianKernel(X,X2,sqrt(sigma2));
Ghat = GaussianKernel(X,X,sqrt(sigma2));
%%
% TODO: truncate to get low-order model.
% use Ghat to compute Sigma (called 'S') and V
% [Qfull,S2full] = eig(Ghat); % since it's true that Ghat*Q = Q*Sigma^2
tic;
[Qfull, S2full] = eigs(Ghat,size(Ghat,1),'largestabs','Tolerance',1e-19); % use increased tolerance with GaussianKernel
toc;
% S2full = flipud(fliplr(S2full)); % make largest singular values first instead of largest singular values last
% Qfull = fliplr(Qfull);
Sfull = sqrtm(S2full);
% truncate
r = size(Ghat,1)-1; % number of dimensions to keep
Q = Qfull(:,1:r);
S = Sfull(1:r,1:r);
figure(2);
clf;
semilogy(diag(Sfull),'o');
hold on
semilogy(diag(S),'.');
% Sinv = diag(ones(size(S,1),1)./diag(S)); % better numerical stability?
% Khat = (Sinv*Q')*Ahat*(Q*Sinv);
Khat = (S\Q')*(Ahat)*(Q/S); % Khat, an approximation of the Koopman operator
[Vhat, Lambda] = eig(Khat);
if(rank(Vhat) < size(Vhat,1))
error('Vhat is not full rank (rank = %d). I didnt code up the case where Vhat isnt.', rank(Vhat));
end
Phi = Q*S*Vhat; % eigenfunctions of Khat
Xi = Phi\X;
X2hat = Phi*Lambda*Xi;
if(norm(imag(X2hat)) > 1e-3)
error('X2hat has non-negligible imaginary component.')
end
X2hat = real(X2hat); % this is quite an accurate reconstruction of X2.
disp('Finished computing kernel-DMD regression. Attempting simulation from IC.');
%% simulate forward
num_modes = size(Q,1); % i THINK this is the correct one
tic;
xhat(1,:) = X(1,:);
% xhat(1,1) = xhat(1,1) + 0.1;
% xhat(1,:) = xhat(1,:) + 0.001*randn(1,size(xhat,2));
for t = 2:2*size(X,1)
% calculate phi given x
prev_x = xhat(t-1,:);
% kernelmat = nan(1,num_modes);
% for M = 1:num_modes
% kernelmat(M) = kernel(prev_x,X(M,:));
% end
kernelmat = GaussianKernel(prev_x,X, sqrt(sigma2));
kernelmat = kernelmat(1:num_modes)';
phi_at_t = kernelmat*Q/S*Vhat;
% estimate
xnext = phi_at_t*Lambda*Xi;
if(norm(imag(xnext)) > 1e-3)
error('xnext isnt real (has non-negligble imaginary component)')
end
xhat(t,:) = real(xnext);
end
toc;
%% Check
figure(1)
clf
subplot(2,1,1)
plot(Xall(:,1));
hold on
plot([nan(1);X2hat(:,1)],'--');
plot(xhat(:,1),'--');
legend('true','1-step recon', 'n-step recon');
title('Reconstruction element 1')
subplot(2,1,2)
plot(Xall(:,6));
hold on
plot([nan(1);X2hat(:,6)],'--');
plot(xhat(:,6),'--');
legend('true','1-step recon', 'n-step recon');
title('Reconstruction element 6');
disp('Simulated from IC.');