add kernel dmd class

examples/todo/hamlyn_heart_video_kernel_dmd.m → examples/hamlyn_heart_video_kernel_dmd.m

@@ -26,7 +26,7 @@
 
 %% Try DMD on it
 clear variables
-load data/HamlynHeartImages.mat;
+load ../data/HamlynHeartImages.mat;
 
 disp('Data preparation...')
 
@@ -92,8 +92,8 @@ toc;
 Sfull = sqrtm(S2full);
 
 % truncate
-r = size(Ghat,1)-1; % number of dimensions to keep
-% r = floor(size(Ghat,1)/2);
+% 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);
 

examples/kernel_dmd.m

+classdef kernel_dmd < handle
+    % kernel_dmd
+    %   States are ROWS instead of columns, like they were in DMD.
+    %   Mostly follows notation from Williams et al. "Kernel-based methods
+    %   for Koopman spectral analysis." Journal of Computation Dynamics
+    %   (2015).
+    
+    properties (SetAccess = private)
+        % Kernel DMD variables
+        X; 
+        X2;
+        Ahat;
+        Ghat;
+        r;
+        Qfull; % left singular vectors of X
+        Sfull; % singular values of X
+        Q; % truncated left singular vectors of X
+        S; % truncated singular values of X
+        Khat; % low-dimensional "weight-space" dynamics
+        Vhat; % eigenvectors of weight-space dynamics
+        Lambda; % eigenvalues of weight-space dynamics
+        Phi; % Koopman eigenfunctions. projects weight-space to data-space (e.g. images)
+        Xi; % Koopman modes
+        
+        % kernel variables
+        kernel_type;
+        sigma2;
+        
+        % variables for preprocessing
+        ProcessSettings; % used for removing and re-inserting constant rows
+        scaling_method;
+        scaleFactor;
+        scaleOffset;
+        
+        % configuration
+        verbose; % controls verbosity
+    end
+    
+    methods
+        function obj = kernel_dmd(data,scaling_method,kernel_type, varargin)
+            %DMD Construct an instance of this class
+            %   The state at each timestep is a column vector)
+            
+            % parse some inputs
+            if(isempty(kernel_type) || kernel_type == "gaussian")
+                obj.kernel_type = kernel_type;
+            else
+                error('Only kernel_type = ''gaussian'' supported right now.')
+            end
+            if nargin == 2
+                obj.verbose = varargin{1};
+            else
+                obj.verbose = false; 
+            end
+            
+            % remove constant rows and scale data
+            [DataProcessed,obj.ProcessSettings] = removeconstantrows(data);
+            if(scaling_method ~= "none")
+                if(isempty(scaling_method))
+                    scaling_method = 'zscore';
+                end
+                obj.scaling_method = scaling_method;
+                [Xall,obj.scaleFactor, obj.scaleOffset] = ScaleData(DataProcessed',scaling_method);
+            end
+            
+            %store processed data in member variables
+            obj.X = Xall(1:end-1,:); % states are ROWS
+            obj.X2 = Xall(2:end,:);
+            obj.r = size(X,2); % no truncation by default
+        end
+        
+        function regression(obj,sigma2,varargin) % argument is truncation rank r
+            if nargin == 2
+                obj.r = varargin{1};
+            end
+            obj.sigma2 = sigma2;
+            
+            tic;
+            % compute left singular values Qfull and singular values Sfull
+            obj.Ahat = GaussianKernel(obj.X,obj.X2);
+            obj.Ghat = GaussianKernel(obj.X,obj.X);
+            [obj.Qfull, S2full] = eigs(obj.Ghat,size(obj.Ghat,1),'largestabs','Tolerance',1e-19); % i had problems with eigs before, needed to increase tolerance
+            obj.Sfull = sqrtm(S2full); % could potentially speed up by using the fact S2full is diagonal
+            
+            % truncate
+            obj.Q = obj.Qfull(:,1:obj.r);
+            obj.S = obj.Sfull(1:obj.r,1:obj.r);
+            
+            % compute the low-dimensional approximation of the Koopman
+            % operator
+            obj.Khat = (obj.S\obj.Q')*(obj.Ahat)*(obj.Q/obj.S);
+            [obj.Vhat, obj.Lambda] = eig(obj.Khat); 
+            obj.Phi = obj.Q*obj.S*obj.Vhat; % projection from weight-space to  data-space
+            obj.Xi = obj.Phi\obj.X; 
+            
+            % tell user how long it took
+            regressionTime = toc;
+            if(obj.verbose)
+                fprintf('DMD regression() took %f seconds.\n',regressionTime);
+            end
+        end
+        
+        function X2hat = reconstruct(obj)
+            % reconstructs an estimate of X2 using the low-rank
+            % "weight-space" dynamics x_r(k+1) = Atilde*x_r(k)
+            % 
+            % useful for getting a feel for how accurate the low-rank
+            % approximation is.
+            
+            if(isempty(obj.Qfull) || isempty(obj.Sfull) || isempty(obj.Vhat))
+                error('Need to run regression() before running reconstruct().');
+            end
+            
+            X2hatComplex = obj.Phi*obj.Lambda*obj.Xi;
+            assert(norm(imag(X2hatComplex)) <= 1e-6,'X2hat has non-negligible imaginary component.');
+            X2hat = real(X2hatComplex);
+            
+        end
+        
+        function xhat = predict(obj,initial_condition,num_timesteps)
+            % predicts num_timesteps timesteps from an initial condition
+            % like reconstruct, but more general.
+            % can be used to predict more timesteps than were in the data
+            
+            tic;
+            if(isempty(obj.U0) || isempty(obj.S0) || isempty(obj.V0))
+                error('Need to run regression() before running predict().');
+            end
+            
+            assert(all(size(initial_condition) == [1,size(obj.X,2)]),'Initial condition must be a ROW vector with the same dimension as a row in the data matrix X, which has length %d',size(obj.X,1)); 
+            xhat = nan(num_timesteps,length(initial_condition));
+            xhat(1,:) = initial_condition;
+            
+            for k = 2:num_timesteps
+                prev_x = xhat(k-1,:);
+                kernelmat = GaussianKernel(prev_x,obj.X);
+                kernelmat = kernelmat(1:num_modes)'; % truncate and turn into row vector
+                phi_at_t = kernelmat*obj.Q/obj.S*obj.Vhat;
+                xhat(k,:) = real(xnext);
+            end
+            
+            assert(norm(imag(xhat)) <= 1e-6,'Predicted xhat has non-negligible imaginary part!');
+            xhat = real(xhat);
+            
+            predictTime = toc;
+            if(obj.verbose)
+                fprintf('Predict() took %f seconds.\n',predictTime);
+            end
+        end
+        
+         function F = GaussianKernel(Xin,Yin)
+            x = Xin./sqrt(obj.sigma2);
+            y = Yin./sqrt(obj.sigma2);
+
+            % 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
+        
+    end
+    
+    methods (Access = private)
+        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
+    end
+end
+

examples/kernel_dmd_heart.m

+%% 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