Code cleanup

README.md

@@ -2,3 +2,11 @@ This repository contains data and code to generate results from the paper
 
 _Residual Selection for Consistency Based Diagnosis Using Machine Learning Models_
 by Erik Frisk and Mattias Krysander, Department of Electrical Engineering, Linköping University, Sweden
+
+The repository includes code in both Matlab and Python (in the form of a Jupyter notebook).
+
+Note that the generated plots are not identical to the results in the paper where a specific Matlab implementation of the machine learning algorithms were used. However, the methodology is the same and the results are similar.
+
+If you have any questions, you are welcome to contact any of the authors
+Erik Frisk (erik.frisk@lliu.se) or
+Mattias Krysander (mattias.krysander@liu.se)

code/main.m

-% TODO
-%
-% 3. Flytta viss kod till förklarande funktioner
-% 5. Uppdatera readme
-% 6. Kommentarer om version av matlab och figurer. Klura ut vilken version
-%    som använts i artikeln.
-% 7. Vilken version av Matlab krävs?
-
 % This script produces the results shown in:
 %
 % Erik Frisk and Mattias Krysander (2018). Residual Selection for
@@ -70,7 +62,7 @@ figure(16);
 IsolabilityAnalysisFSM(data.fsm(:, 2:8), data.modes(2:8), 'permute', 0);
 title('Structural isolability matrix using all available tests');
 ylabel('Injected fault', 'FontWeight', 'bold');
-xlabel('Diagnsoed fault', 'FontWeight', 'bold');
+xlabel('Diagnosed fault', 'FontWeight', 'bold');
 
 %% Create thresholded residuals
 thdata = data;
@@ -153,40 +145,16 @@ xlabel('Predictors','FontWeight', 'bold');
 clear imp
 
 %% Performance evaluation of selected tests using consistency based diagnosis.
-
-no_samples = size(data.res, 1); % n 
-no_modes = numel(data.modes);  % m 
-no_of_samples_per_mode = no_samples/no_modes; % k
-
-structuralIsolationMatrix = ...
-  IsolabilityAnalysisFSM(data.fsm(:, 2:8), data.modes(2:8), 'permute', 0);
-
-stim2 =  sum(((2.^(0:6)).*structuralIsolationMatrix), 2); 
-probabilityMaximumIsolation = zeros(42, no_modes-1);
-
-C = cell(length(residualRanking), 1);
 isolationerror = zeros(length(residualRanking), 1);
 mean_md = zeros(length(residualRanking), 1);
 fa = zeros(length(residualRanking), 1);
+probabilityMaximumIsolation = zeros(42, numel(data.modes)-1);
+
 for i = 1:length(residualRanking)
-    selected_residuals = residualRanking(1:i); % select the i best residuals
-    dx = SingleFaultIsolation(data.res(:, selected_residuals),...
-        data.fsm(selected_residuals, (2:8))); % apply CBD.
-    % Compute for each mode the probability of maximum isolation
-    dx2 = sum(((2.^(0:6)).*dx(:, (2:end))), 2);
-    for j = 1:no_modes-1 % not NF
-        probabilityMaximumIsolation(i,j) = sum(dx2((1:no_of_samples_per_mode)+...
-            j*no_of_samples_per_mode)==stim2(j))/no_of_samples_per_mode; 
-    end
-    % Compute confusion matrix of CBD using the i best residuals
-    C{i} = CBDPlots(data, selected_residuals);
-    % Compute performance indicators
-    isolationerror(i) = sum(sum(abs(structuralIsolationMatrix-C{i}(2:8,2:8)./((1-C{i}(2:end,1))*ones(1,7)))))/7^2;
-    mean_md(i) = mean(C{i}(2:end, 1));
-    fa(i) = 1-C{i}(1, 1);
+  selected_residuals = residualRanking(1:i); % select the i best residuals
+  [mean_md(i), fa(i), isolationerror(i), probabilityMaximumIsolation(i, :)] = ...
+    PerformanceEval(data, selected_residuals);
 end
-clear stim2 dx dx2 no_samples no_modes no_of_samples_per_mode selected_residuals C
-clear structuralIsolationMatrix i j 
 
 %% Plot Fig 11: False alarm probability, mean missed detection probability,
 % and aggregated isolation error as a function of selected tests.

code/src/CBDPlots.m

-function C = CBDPlots(data, resIdx)
+function C = CBDPlots(data, resIdx, varargin)
 % CBDPlots Plot result from consistency based diagnosis based on residuals
 %
 %  C = CBDPlots(data, resIdx)
 %
 %  Inputs:
-%    data    - residual data structure
-%    resIdx  - Index to residuals to evaluate 
-%  
+%    data         - residual data structure
+%    resIdx       - Index to residuals to evaluate 
+%    compute_only - only compute the confusion matrix, do not plot
+%                   (default: false)
 %  Outputs:
 %    C       - Confusion matrix
 %
 
+  if nargin < 3
+    compute_only = false;
+  else
+    compute_only = varargin{1};
+  end
+
   FSM42 = data.fsm(:, 2:end);
 
   %% Evaluate selected set of residuals
@@ -28,29 +35,31 @@ function C = CBDPlots(data, resIdx)
       C(fj,fi) = sum(dx(fjIdx,fi))/sum(fjIdx);
     end
   end
-
+  
   %% Plot confusion matrix in a 3D-bar plot
-  b=bar3(C);
-  set(gca, 'XTickLabel', data.modes, 'TickLabelInterpreter', 'none')
-  set(gca, 'YTickLabel', data.modes, 'TickLabelInterpreter', 'none')
-  xlabel('Diagnosed fault', 'FontWeight', 'bold');
-  ylabel('Injected fault', 'FontWeight', 'bold');
-  zlabel('P(fi diag|fj)');
-  title('Fault Isolation Performance Matrix')
+  if ~compute_only
+    b=bar3(C);
+    set(gca, 'XTickLabel', data.modes, 'TickLabelInterpreter', 'none')
+    set(gca, 'YTickLabel', data.modes, 'TickLabelInterpreter', 'none')
+    xlabel('Diagnosed fault', 'FontWeight', 'bold');
+    ylabel('Injected fault', 'FontWeight', 'bold');
+    zlabel('P(fi diag|fj)');
+    title('Fault Isolation Performance Matrix')
 
-  for finjIdx=1:nf
-    for fdiagIdx=1:nf
-      text(fdiagIdx,finjIdx,C(finjIdx,fdiagIdx)+0.05,...
-        sprintf('%.1f',C(finjIdx,fdiagIdx)*100),'HorizontalAlignment', 'center')
+    for finjIdx=1:nf
+      for fdiagIdx=1:nf
+        text(fdiagIdx,finjIdx,C(finjIdx,fdiagIdx)+0.05,...
+          sprintf('%.1f',C(finjIdx,fdiagIdx)*100),'HorizontalAlignment', 'center')
+      end
     end
+    for k = 1:length(b)
+        zdata = b(k).ZData;
+        b(k).CData = zdata;
+        b(k).FaceColor = 'interp';
+    end
+    colormap('Summer')
+    view(0,90)
   end
-  for k = 1:length(b)
-      zdata = b(k).ZData;
-      b(k).CData = zdata;
-      b(k).FaceColor = 'interp';
-  end
-  colormap('Summer')
-  view(0,90)
 end
 
 function [dx] = EvaluateResiduals(data, FSM, resIdx)

code/src/PerformanceEval.m

+function [mean_md, fa, fi, pma] = PerformanceEval(data, residuals)
+  no_samples = size(data.res, 1); % n 
+  no_modes = numel(data.modes);  % m 
+  no_of_samples_per_mode = no_samples/no_modes; % k
+
+  structuralIsolationMatrix = ...
+    IsolabilityAnalysisFSM(data.fsm(:, 2:8), data.modes(2:8), 'permute', 0);
+
+  stim2 =  sum(((2.^(0:6)).*structuralIsolationMatrix), 2); 
+  pma = zeros(1, no_modes-1);
+
+
+  dx = SingleFaultIsolation(data.res(:, residuals),...
+      data.fsm(residuals, (2:8))); % apply CBD.
+    
+  % Compute for each mode the probability of maximum isolation
+  dx2 = sum(((2.^(0:6)).*dx(:, (2:end))), 2);
+  for j = 1:no_modes-1 % not NF
+      pma(j) = sum(dx2((1:no_of_samples_per_mode)+...
+          j*no_of_samples_per_mode)==stim2(j))/no_of_samples_per_mode; 
+  end
+  % Compute confusion matrix of CBD using the i best residuals
+  C = CBDPlots(data, residuals, true);
+  
+  % Compute performance indicators
+  fi = sum(sum(abs(structuralIsolationMatrix-C(2:8,2:8)./((1-C(2:end,1))*ones(1,7)))))/7^2;
+  mean_md = mean(C(2:end, 1));
+  fa = 1-C(1);
+end
\ No newline at end of file