Merge pull request #211 from mldiego/master

Fix precision error on caused by MaxPooling reachability

code/nnv/engine/nn/layers/MaxPooling2DLayer.m

@@ -305,14 +305,14 @@ function set_padding(obj, padding)
                 % input volume has only one channel                
                 h = n(1); % height of input
                 w = n(2); % width of input 
-                padded_I = zeros(t + h + b, l + w + r);
+                padded_I = cast(zeros(t + h + b, l + w + r), 'like', input);
                 padded_I(t+1:t+h,l+1:l+w) = input; % new constructed input volume
             elseif length(n) > 2
                 % input volume may has more than one channel
                 h = n(1); % height of input
                 w = n(2); % width of input 
                 d = n(3); % depth of input (number of channel)
-                padded_I = zeros(t + h + b, l + w + r, d); % preallocate 3D matrix
+                padded_I = cast(zeros(t + h + b, l + w + r, d), 'like', input); % preallocate 3D matrix
                 for i=1:d
                     padded_I(t+1:t+h, l+1:l+w, i) = input(:, :, i); % new constructed input volume
                 end              
@@ -336,7 +336,7 @@ function set_padding(obj, padding)
                 c = obj.get_zero_padding_input(ims.V(:,:,:,1));
                 k = size(c);
                 n = ims.numPred;
-                V1(:,:,:,n+1) = zeros(k);
+                V1(:,:,:,n+1) = cast(zeros(k), 'like', ims.V);
                 for i=1:n
                     V1(:,:,:,i+1) = obj.get_zero_padding_input(ims.V(:,:,:,i+1));
                 end
@@ -371,7 +371,7 @@ function set_padding(obj, padding)
 
             % compute feature map for each cell of map do it in parallel using cpu or gpu
             % TODO: explore power of using GPU for this problem
-            maxMap = zeros(1, h*w); % using single vector for parallel computation
+            maxMap = cast(zeros(1, h*w), 'like', input); % using single vector for parallel computation
 
             for l=1:h*w
                 a = mod(l, w);
@@ -466,7 +466,7 @@ function set_padding(obj, padding)
             end
 
             [h, w] = obj.get_size_maxMap(input.V(:,:,1,1));
-            new_V(:,:,input.numChannel, input.numPred+1) = zeros(h,w);
+            new_V(:,:,input.numChannel, input.numPred+1) = cast(zeros(h,w), 'like', input.V);
 
             % get max points for each channel
             channel_maxPoints(:,:,input.numChannel) = zeros(3, h*w);
@@ -540,7 +540,7 @@ function set_padding(obj, padding)
 
             % check max_id first           
             max_index = cell(h, w, pad_image.numChannel);
-            maxMap_basis_V(:,:,pad_image.numChannel, pad_image.numPred+1) = zeros(h,w); % pre-allocate basis matrix for the maxmap
+            maxMap_basis_V(:,:,pad_image.numChannel, pad_image.numPred+1) = cast(zeros(h,w), 'like', in_image.V); % pre-allocate basis matrix for the maxmap
             split_pos = [];
 
             % compute max_index and split_index when applying maxpooling operation
@@ -813,7 +813,7 @@ function set_padding(obj, padding)
                 fprintf('\n%d new variables are introduced\n', l);
             end                   
             % update new basis matrix
-            new_V(:,:,pad_image.numChannel,np+1) = zeros(h,w);
+            new_V(:,:,pad_image.numChannel,np+1) = cast(zeros(h,w), 'like', in_image.V);
             new_pred_index = 0;
             for k=1:pad_image.numChannel
                 for i=1:h
@@ -835,10 +835,10 @@ function set_padding(obj, padding)
 
             %update constraint matrix C and constraint vector d
             N = obj.PoolSize(1) * obj.PoolSize(2);            
-            new_C = zeros(new_pred_index*(N + 1), np);
-            new_d = zeros(new_pred_index*(N + 1), 1);
-            new_pred_lb = zeros(new_pred_index, 1); % update lower bound and upper bound of new predicate variables
-            new_pred_ub = zeros(new_pred_index, 1); 
+            new_C = cast(zeros(new_pred_index*(N + 1), np), 'like', in_image.V);
+            new_d = cast(zeros(new_pred_index*(N + 1), 1), 'like', in_image.V);
+            new_pred_lb = cast(zeros(new_pred_index, 1), 'like', in_image.V); % update lower bound and upper bound of new predicate variables
+            new_pred_ub = cast(zeros(new_pred_index, 1), 'like', in_image.V);
             new_pred_index = 0;
             for k=1:pad_image.numChannel
                 for i=1:h
@@ -849,14 +849,14 @@ function set_padding(obj, padding)
                             new_pred_index = new_pred_index + 1;                           
                             startpoint = startPoints{i,j};
                             points = pad_image.get_localPoints(startpoint, obj.PoolSize);
-                            C1 = zeros(1, np);
+                            C1 = cast(zeros(1, np), 'like', in_image.V);
                             C1(pad_image.numPred + new_pred_index) = 1;
                             [lb, ub] = pad_image.get_localBound(startpoint, obj.PoolSize, k, lp_solver);
                             new_pred_lb(new_pred_index) = lb;
                             new_pred_ub(new_pred_index) = ub;
                             d1 = ub;                           
-                            C2 = zeros(N, np);
-                            d2 = zeros(N, 1);
+                            C2 = cast(zeros(N, np),'like', in_image.V);
+                            d2 = cast(zeros(N, 1), 'like', in_image.V);
                             for g=1:N
                                 point = points(g,:);
                                 % add new predicate constraint: xi - y <= 0
@@ -876,7 +876,7 @@ function set_padding(obj, padding)
             end
 
             n = size(pad_image.C, 1);
-            C = [pad_image.C zeros(n, new_pred_index)];
+            C = [pad_image.C cast(zeros(n, new_pred_index), 'like', in_image.V)];
             new_C = [C; new_C];
             new_d = [pad_image.d; new_d];
             new_pred_lb = [pad_image.pred_lb; new_pred_lb];

code/nnv/examples/NN/medmnist/debugInconsistenciesNodulemnist.m

@@ -0,0 +1,56 @@
+%% Debug inconsistencies for nodulemnist verification results
+
+% get current dataset to verify
+% dataset = medmnist_path + datasets(i).name;
+dataset = "data/mat_files/nodulemnist3d.mat";
+
+
+disp("Begin verification of nodulemnist3d");
+
+% Load data
+load(dataset);
+
+% data to verify (test set)
+test_images = permute(test_images, [2 3 4 5 1]);
+test_labels = test_labels + 1;
+
+% load network
+load("models/model_nodulemnist3d.mat");
+matlabNet = net;
+net = matlab2nnv(net);
+
+% adversarial attacks
+names = ["dark"];
+max_pixels = [50;100;200];
+noise_vals = [1];
+
+% select volumes to verify
+N = 200;
+inputs = test_images(:,:,:,:,1:N);
+targets = test_labels(1:N);
+
+% Initialize results 
+results = zeros(2, N, length(names), length(max_pixels), length(noise_vals));
+outputSets = cell(3,1);
+
+% verify volumes with all attack combos
+for a=1:length(names)
+    for b=1:length(max_pixels)
+        for c=1:length(noise_vals)
+            % create attack from variables
+            adv_attack = struct;
+            adv_attack.Name = names(a);
+            adv_attack.max_pixels = max_pixels(b);
+            adv_attack.noise_de = noise_vals(c);
+            adv_attack.threshold = 150;
+
+            % Compute verification
+            [outputSets{b}, results(:,:,a,b,c)] = verify_medmnist3d_extraInfo(net, inputs, targets, adv_attack);
+        end
+    end
+end
+
+% save results
+save("results/verification_multipleAttacks_nodulemnist3d_debug.mat", "results", "outputSets");
+
+

code/nnv/examples/NN/medmnist/process_debugNodule.m

@@ -0,0 +1,93 @@
+%% Let's take a look into this
+
+load('results/verification_multipleAttacks_nodulemnist3d_debug.mat');
+
+% Only get verification results
+results = results(1,:,:,:);
+results = squeeze(results);
+results = results';
+outputSets = outputSets';
+
+% first discrepancy (5)
+R1 = outputSets(1,5); % 
+[lb1, ub1] = R1.getRanges;
+R2 = outputSets(2,5);
+[lb2, ub2] = R1.getRanges;
+R3 = outputSets(3,5);
+[lb3, ub3] = R1.getRanges;
+
+% So the ranges are the same, same results should be getting...
+% What is happening?
+
+% Nodulemnist is a binary classification problem
+% The bounds are clearly apart, so no unkown should be happening...
+% Are we representing the robustness halfspace incorrectly?
+
+% Verification property
+G = [1,-1]; g = 0;
+Hs = HalfSpace(G,g);
+
+res1 = verify_specification(R1, Hs);
+res2 = verify_specification(R2, Hs);
+res3 = verify_specification(R3, Hs);
+
+% bounds are the same, but the verification results are not...
+% why is there no intersections? Does the intersection method have a bug on
+% it?
+
+% Let's visualize the output sets (exact)
+% figure;
+% Star.plot(R1,'r');
+% figure;
+% Star.plot(R2,'r');
+% figure;
+% Star.plot(R3,'r');
+
+% Are the predicate bounds incorrect? They show 254 and 255... 
+% This seems just wrong... Can we simply remove the predicate bounds in 
+% all cases? I don't think so...
+% Then we need to fix either the how the predicate bounds get computed, or
+% do not use them for certain operations such as plotting or computing
+% intersections/empty sets.
+
+% Remove predicate bounds
+R1.predicate_lb = [];
+R1.predicate_ub = [];
+R2.predicate_lb = [];
+R2.predicate_ub = [];
+R3.predicate_lb = [];
+R3.predicate_ub = [];
+
+% Check verification results now
+res11 = verify_specification(R1, Hs);
+res22 = verify_specification(R2, Hs);
+res33 = verify_specification(R3, Hs);
+
+% We get the same results, so I guess the predicate bounds do not affect
+% this verification result (intersection of Star sets)
+
+% The difference between R3 and the other 2 is that there is an extra
+% constraint. Is this messing things up?
+% Let's manually visualize the stars...
+
+% get centers
+c1 = R1.V(:,1);
+c2 = R2.V(:,1);
+c3 = R3.V(:,1);
+
+% get basis vectors
+V1 = R1.V(:,2:end);
+V2 = R2.V(:,2:end);
+V3 = R3.V(:,2:end);
+
+% get constraints
+C1 = R1.C;
+d1 = R1.d;
+C2 = R2.C;
+d2 = R2.d;
+C3 = R3.C;
+d3 = R3.d;
+
+% can we get the bounds of the stars from here?
+R4 = Star(lb1,ub1);
+Star.plot(R4, 'r');

code/nnv/examples/NN/medmnist/process_results.m

@@ -1,73 +1,115 @@
 %% Process some of the results
 
 % for now, we just want to process the multiple attack 3D results
-% The results files are 5D ([#images, results, 
+% The results files are 5D 
+%   [# images, results, attacktype, # pixels changed, perturbation size]
 
-res_files = [...
-    "results/verification_multipleAttacks_adrenalmnist3d.mat";
-    "results/verification_multipleAttacks_fracturemnist3d.mat";
-    "results/verification_multipleAttacks_nodulemnist3d.mat";
-    "results/verification_multipleAttacks_organmnist3d.mat";
-    "results/verification_multipleAttacks_synapsemnist3d.mat";
-    "results/verification_multipleAttacks_vesselmnist3d.mat"];
+% res_files = [...
+%     "results/verification_multipleAttacks_adrenalmnist3d.mat";
+%     "results/verification_multipleAttacks_fracturemnist3d.mat";
+%     "results/verification_multipleAttacks_nodulemnist3d.mat";
+%     "results/verification_multipleAttacks_organmnist3d.mat";
+%     "results/verification_multipleAttacks_synapsemnist3d.mat";
+%     "results/verification_multipleAttacks_vesselmnist3d.mat"];
+
+res_files = "results/verification_multipleAttacks_nodulemnist3d.mat";
 
 nFiles = length(res_files);
+
+% Initialize var to summarize results
 res_summary = zeros(nFiles, 5);
+
 % Write results to a file
-diary multipleAttacks_3D_results.txt
+% diary multipleAttacks_3D_results.txt
+
+% Begin processing
 for i=1:nFiles
+
+    % Load data
     load(res_files(i));
+
+    % Results title
     disp("*******************************************************");
     disp("================= PROCESSING RESULTS: " + res_files(i) + " ...");
+
+    % Process results for all atacks in file
     process_multiple_attacks(results);
+
+    % Add end of table
     disp("*******************************************************");
+
 end
-diary off;
+
+% Close results file
+% diary off;
 
 
 %% Helper functions
 
 function process_multiple_attacks(results)
-    % results: 5D array
+
+    % ensure results: 5D array
     n = size(results);
     if length(n) ~= 5
         error("Wrong input.")
     end
-    % attack combinations
+
+    % remember attack combinations
     names = ["dark";"bright"];
     max_pixels = [50;100;200];
     noise_vals = [1;2;3];
-    for i = 1:n(3)
-        for j = 1:n(4)
-            for k = 1:n(5)
+
+    % Begin processing results one attack at a time
+    for i = 1:n(3)         % names (dark/bright)
+        for j = 1:n(4)     % max_pixels (50/100/200)
+            for k = 1:n(5) % noise_val (perturbation size)
+
                 % Print what results we are looking into
                 disp("... processing " + names(i) + " attack  with " ...
                     +string(max_pixels(j)) + " pixels perturbed with noise of "+ string(noise_vals(k)));
+
                 % process individual results
                 print_results(results(:,:,i,j,k));
+
+                % add empty line for readability
                 disp(" ");
-            end
-        end
-    end
+
+            end % k (perturbation size)
+        end     % j (# pixels)
+    end         % i (attack type)
+
 end
 
 % Print results to screen
 function res_summary = print_results(results)
+
     % print results for every attack
     N = size(results,2); 
     disp("----------- ROBUSTNESS RESULTS -------------");
     disp("Verification results of " + string(N) + " images.");
-    rob = sum(results(1,:) == 1); % robust
-    disp("Number of robust images          =  " + string(rob));
-    not_rob = sum(results(1,:) == 0); % not robust
-    disp("Number of not robust images      =  " + string(not_rob));
-    unk = sum(results(1,:) == 2); % unknown
-    disp("Number of unknown images         =  " + string(unk));
-    missclass = sum(results(1,:) == -1); % misclassified
-    disp("Number of missclassified images  =  " + string(missclass));
-    avgTime = sum(results(2,:))/N; % average computation time
-    disp("Average computation time of "         + string(avgTime));
+
+    % 1) Robust
+    rob = sum(results(1,:) == 1); 
+    disp("Number of robust images         =  " + string(rob));
+
+    % 0) Not Robust
+    not_rob = sum(results(1,:) == 0); 
+    disp("Number of not robust images     =  " + string(not_rob));
+
+    % 2) Unknown
+    unk = sum(results(1,:) == 2);
+    disp("Number of unknown images        =  " + string(unk));
+
+    % -1) Misclassified
+    missclass = sum(results(1,:) == -1); 
+    disp("Number of missclassified images =  " + string(missclass));
+
+    % Average computation time
+    avgTime = sum(results(2,:))/N; 
+    disp("Average computation time of        " + string(avgTime));
+
     % Get summary results 
     res_summary = [rob, not_rob, unk, missclass];
+
 end