wip: kernel_dmd running, unverified.

examples/kernel_dmd.m

@@ -31,6 +31,8 @@ classdef kernel_dmd < handle
         scaling_method;
         scaleFactor;
         scaleOffset;
+        imageWidth;
+        imageHeight;
         
         % configuration
         verbose; % controls verbosity
@@ -47,14 +49,28 @@ classdef kernel_dmd < handle
             else
                 error('Only kernel_type = ''gaussian'' supported right now.')
             end
-            if nargin == 2
+            if nargin == 4
                 obj.verbose = varargin{1};
             else
                 obj.verbose = false; 
             end
             
+            % vectorize 3D matrix of 2D image data
+            obj.imageHeight = size(data,1);
+            obj.imageWidth = size(data,2);
+            
+            vectorizedImages = [];
+            for k = 1:size(data,3)
+                image_k = data(:,:,k);
+                vectorizedImages = [vectorizedImages, double(image_k(:))];
+            end
+            
+            % construct "second order state"
+            data2 = [vectorizedImages(:,2:end);
+                     vectorizedImages(:,1:end-1)];
+            
             % remove constant rows and scale data
-            [DataProcessed,obj.ProcessSettings] = removeconstantrows(data);
+            [DataProcessed,obj.ProcessSettings] = removeconstantrows(data2);
             if(scaling_method ~= "none")
                 if(isempty(scaling_method))
                     scaling_method = 'zscore';
@@ -64,21 +80,21 @@ classdef kernel_dmd < handle
             end
             
             %store processed data in member variables
-            obj.X = Xall(1:end-1,:); % states are ROWS
+            obj.X = Xall(1:end-1,:); % states are ROWS now
             obj.X2 = Xall(2:end,:);
-            obj.r = size(X,2); % no truncation by default
+            obj.r = size(obj.X,1); % no truncation by default
         end
         
         function regression(obj,sigma2,varargin) % argument is truncation rank r
-            if nargin == 2
+            if nargin == 3
                 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.Ahat = obj.GaussianKernel(obj.X,obj.X2);
+            obj.Ghat = obj.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
             
@@ -100,7 +116,7 @@ classdef kernel_dmd < handle
             end
         end
         
-        function X2hat = reconstruct(obj)
+        function [X2hat,imageSeq] = reconstruct(obj)
             % reconstructs an estimate of X2 using the low-rank
             % "weight-space" dynamics x_r(k+1) = Atilde*x_r(k)
             % 
@@ -113,8 +129,10 @@ classdef kernel_dmd < handle
             
             X2hatComplex = obj.Phi*obj.Lambda*obj.Xi;
             assert(norm(imag(X2hatComplex)) <= 1e-6,'X2hat has non-negligible imaginary component.');
-            X2hat = real(X2hatComplex);
+            X2hatScaled = real(X2hatComplex);
             
+            X2hat = obj.ReapplyScaling(X2hatScaled);
+            imageSeq = obj.VecSeq2ImSeq(X2hat);
         end
         
         function xhat = predict(obj,initial_condition,num_timesteps)
@@ -136,7 +154,7 @@ classdef kernel_dmd < handle
                 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);
+                xhat(k,:) = real(phi_at_t*obj.Lambda*obj.Xi);
             end
             
             assert(norm(imag(xhat)) <= 1e-6,'Predicted xhat has non-negligible imaginary part!');
@@ -148,54 +166,36 @@ classdef kernel_dmd < handle
             end
         end
         
-         function F = GaussianKernel(Xin,Yin)
+         function F = GaussianKernel(obj,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
+            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
         
-    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);
+         function rescaledData = ReapplyScaling(obj,data)
+            for k = 1:size(data,1) % for each state, which is a ROW
+                xhatRescaled(k,:) = data(k,:).*obj.scaleFactor + obj.scaleOffset;
             end
-
-        end
+            xhatRescaled = xhatRescaled';
+            rescaledData = removeconstantrows('reverse',xhatRescaled,obj.ProcessSettings); % add removed rows back in.
+         end
+         
+         % Vector sequence to image sequence
+         function imageSequence = VecSeq2ImSeq(obj,vectorSequence)
+             for k = 1:size(vectorSequence,2) % now each state/image is a column, to comply with removeconstantrows()
+                 imageSequence(:,:,k) = reshape(vectorSequence(1:obj.imageHeight*obj.imageWidth,k),obj.imageHeight,obj.imageWidth);
+             end
+         end
+        
     end
 end
 

examples/kernel_dmd_heart.m

@@ -26,7 +26,7 @@
 
 %% Try DMD on it
 clear variables
-load ../data/HamlynHeartImages.mat;
+load data/HamlynHeartImages.mat;
 
 disp('Data preparation...')
 
@@ -49,6 +49,12 @@ XallPrescaled = [allImages(:,2:end);
 [XallPrescaledProcessed,ProcessSettings] = removeconstantrows(XallPrescaled);
 XallPrescaledProcessed = XallPrescaledProcessed';
 
+
+%% do Kernel DMD via object implementation
+D = kernel_dmd(grayImages(:,:,image_start:image_start+num_images),'','gaussian','true');
+D.regression(1e6,49);
+rec = D.reconstruct();
+
 disp('DMD Regression...')
 %% kernel dmd
 disp('Scaling data...')
@@ -64,7 +70,7 @@ 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?
+sigma2 = 1e6; 
 kernel = @(x1,x2) exp(-norm(x1-x2)^2/sigma2); % Gausssian kernel
 
 % for j = 1:size(X,1)
@@ -75,7 +81,7 @@ kernel = @(x1,x2) exp(-norm(x1-x2)^2/sigma2); % Gausssian kernel
 %     end 
 % end
 
-
+tic;
 Ahat = GaussianKernel(X,X2,sqrt(sigma2));
 Ghat = GaussianKernel(X,X,sqrt(sigma2));
 
@@ -84,9 +90,9 @@ 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);
@@ -114,7 +120,7 @@ end
 
 Phi = Q*S*Vhat; % eigenfunctions of Khat 
 Xi = Phi\X; 
-
+toc;
 X2hat = Phi*Lambda*Xi;
 if(norm(imag(X2hat)) > 1e-3)
     error('X2hat has non-negligible imaginary component.')
@@ -177,8 +183,8 @@ for k = 1:size(xhat,1)
 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)
+
+for k = 1:size(xhatRescaled,2) % now each state/image is a column, to comply with removeconstantrows()
     reconstruction(:,:,k) = reshape(xhatRescaled(1:size(allImages,1),k),size(grayImages,1),size(grayImages,2));
 end