Slight cleanup of code and added SingleFaultIsolation.m to src directory

README.md

 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 Univeristy, Sweden
\ No newline at end of file
+by Erik Frisk and Mattias Krysander, Department of Electrical Engineering, Linköping University, Sweden

code/main.m

 % This script produces the results shown in:
 %
-% Erik Frisk and Mattias Krysander (2018). Residaul Selection for
+% Erik Frisk and Mattias Krysander (2018). Residual Selection for
 % Consistency Based Diagnosis Using Machine Learning Models. In Safeprocess
 % 2018, Warsaw, Polen.
 %
@@ -52,47 +52,48 @@ plotResiduals(absdata)
 clear absdata
 
 %% Plot Fig. 5: Probability of residual alarm given a given mode
-no_residuals = size(data.res,2); 
-all_residuals =[1:no_residuals]; % select all residuals
-plotDecisionMatrix(all_residuals,data);
+no_residuals = size(data.res, 2);
+all_residuals = 1:no_residuals; % select all residuals
+plotDecisionMatrix(all_residuals, data);
 clear no_residuals ans
 
 %% Plot Fig 6: Structural isolability of all available residuals
 load model.mat;
-figure;
-model.IsolabilityAnalysisFSM(data.fsm(:,[2:8]),'permute',0);
+figure(16);
+model.IsolabilityAnalysisFSM(data.fsm(:, 2:8), 'permute', 0);
 title('Structural isolability matrix using all available tests');
-ylabel('Injected fault','FontWeight', 'bold');
-xlabel('Diagnsoed fault','FontWeight', 'bold');
+ylabel('Injected fault', 'FontWeight', 'bold');
+xlabel('Diagnsoed fault', 'FontWeight', 'bold');
 clear ans
+
 %% Create thresholded residuals
 thdata = data;
-thdata.res = abs(data.res)>=1;
+thdata.res = abs(data.res) >= 1;
 
 %% Plot Fig 7: Confusion matrix obtained by consistency-based diagnosis (CBD) using all thresholded residuals
-figure;
-CBDPlots(thdata,model,all_residuals);
+figure(17);
+CBDPlots(thdata, model, all_residuals);
 clear all_residuals ans
 
 %% Train a random forest of 100 classification trees using the all data.
 trees = 100;
 rng(1); % For reproducibility
-t = templateTree('PredictorSelection','interaction-curvature');
+t = templateTree('PredictorSelection', 'interaction-curvature');
 
 fprintf('Building random forest classifier ... ');
   tic;
-Mdl = fitcensemble(single(thdata.res),thdata.mode,'Method','bag',...
-    'NumLearningCycles',trees,'Learners',t);
+Mdl = fitcensemble(single(thdata.res), thdata.mode, 'Method', 'bag',...
+    'NumLearningCycles', trees, 'Learners', t);
 fprintf(' Finished in %.2f sec.\n', toc);
 clear t trees
 
 %% Plot Fig 8: Random forest classification performance using all tests on out-of-bag data.
 Y = oobPredict(Mdl);
-Ival = confusionmat(categorical(thdata.mode),categorical(Y));
+Ival = confusionmat(categorical(thdata.mode), categorical(Y));
 nf = 8;
-C = Ival./(sum(Ival')'*ones(1,nf));
+C = Ival./(sum(Ival, 2)*ones(1, nf));
 
-figure;
+figure(18);
 b = bar3(C);
 title('Confusion matrix - Validation data set')
 set(gca, 'XTickLabel', data.modes, 'TickLabelInterpreter', 'none')
@@ -104,8 +105,9 @@ 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')
+    text(fdiagIdx, finjIdx, C(finjIdx, fdiagIdx) + 0.05,...
+      sprintf('%.1f', C(finjIdx,fdiagIdx)*100), ...
+              'HorizontalAlignment', 'center')
   end
 end
 for k = 1:length(b)
@@ -114,15 +116,15 @@ for k = 1:length(b)
     b(k).FaceColor = 'interp';
 end
 colormap('Summer')
-view(0,90)
+view(0, 90)
 clear Y Ival nf C b fdiagIdx finjIdx k zdata thdata
 
 %% Plot Fig 9: Random forest out-of-bag classification error as a function of the number of trees.
-figure;
-oobErr = oobLoss(Mdl,'mode','cumulative');
-plot(oobErr,'LineWidth',2)
-xlabel( 'Number of grown trees' , 'FontWeight', 'bold');
-ylabel( 'Out-of-bag classification error' , 'FontWeight', 'bold');
+figure(19);
+oobErr = oobLoss(Mdl, 'mode', 'cumulative');
+plot(oobErr, 'LineWidth', 2)
+xlabel('Number of grown trees', 'FontWeight', 'bold');
+ylabel('Out-of-bag classification error', 'FontWeight', 'bold');
 clear oobErr
 
 %% Estimate residual importance by permuting out-of-bag observations.  
@@ -132,13 +134,14 @@ importance = oobPermutedPredictorImportance(Mdl);
 fprintf(' Finished in %.2f sec.\n', toc);
 
 %% Plot Fig 10: Residual importance
-[imp,residualRanking] = sort(importance,'descend');
+[imp,residualRanking] = sort(importance, 'descend');
 
-figure;
-plot(imp,'LineWidth',2)
-xticks([1:length(residualRanking)])
+figure(20);
+plot(imp,'LineWidth', 2)
+xticks(1:length(residualRanking))
 axis([0 43 -0.4 1.6])
-set(gca, 'XTickLabel', residualRanking, 'TickLabelInterpreter', 'none', 'FontSize',8)
+set(gca, 'XTickLabel', residualRanking, 'TickLabelInterpreter', 'none', ...
+    'FontSize', 8)
 title('Out-of-Bag Permuted Predictor Importance Estimates');
 ylabel('Estimates','FontWeight', 'bold');
 xlabel('Predictors','FontWeight', 'bold');
@@ -146,80 +149,87 @@ clear imp
 
 %% Performance evaluation of selected tests using consistency based diagnosis.
 
-no_samples = size(data.res,1); % n 
+no_samples = size(data.res, 1); % n 
 no_modes = numel(data.modes);  % m 
 no_of_samples_per_mode = no_samples/no_modes; % k
 
-structuralIsolationMatrix = model.IsolabilityAnalysisFSM(data.fsm(:,[2:8]),'permute',0);
+structuralIsolationMatrix = ...
+  model.IsolabilityAnalysisFSM(data.fsm(:, 2:8),'permute', 0);
+
+stim2 =  sum(((2.^(0:6)).*structuralIsolationMatrix), 2); 
+probabilityMaximumIsolation = zeros(42, no_modes-1);
 
-stim2 =  sum(((2.^[0:6]).*structuralIsolationMatrix)'); 
-probabilityMaximumIsolation = zeros(42,no_modes-1);
-C ={};
+C = cell(length(residualRanking), 1);
+isolationerror = zeros(length(residualRanking), 1);
+mean_md = zeros(length(residualRanking), 1);
+fa = zeros(length(residualRanking), 1);
 for i = 1:length(residualRanking)
-    selected_residauls = residualRanking(1:i); % select the i best residuals
-    dx = SingleFaultIsolation(data.res(:,selected_residauls),...
-        data.fsm(selected_residauls,[2:8])); % apply CBD.
+    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]))');
+    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]+...
+        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,model,selected_residauls);
+    C{i} = CBDPlots(data, model, 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);
+    mean_md(i) = mean(C{i}(2:end, 1));
+    fa(i) = 1-C{i}(1, 1);
 end
-clear stim2 dx dx2 no_samples no_modes no_of_samples_per_mode selected_residauls C
+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.
-figure;
-plot([mean_md;fa;isolationerror]','LineWidth',2)
+figure(21);
+plot([mean_md, fa, isolationerror], 'LineWidth', 2)
 hold
-selected_no_of_residauls = [12 14 26 27];
-plot(selected_no_of_residauls,isolationerror(selected_no_of_residauls),'x',...
-    'LineWidth',2,'MarkerSize',8,'Color','k')
+selected_no_of_residuals = [12 14 26 27];
+plot(selected_no_of_residuals,isolationerror(selected_no_of_residuals), ...
+     'x', 'LineWidth', 2, 'MarkerSize', 8, 'Color', 'k')
 
-legend('Missed detection probability','False alarm probability','Fault isolation errors')
-xlabel('Number of selected residauls','FontWeight', 'bold')
-ylabel('Probability','FontWeight', 'bold')
+legend('Missed detection probability', 'False alarm probability', ...
+       'Fault isolation errors')
+xlabel('Number of selected residuals','FontWeight', 'bold')
+ylabel('Probability', 'FontWeight', 'bold')
 grid on
-clear mean_md isolationerror selected_no_of_residauls
+clear mean_md isolationerror selected_no_of_residuals
 
 %% Plot Fig 12: Consistency based isolation performance per mode as a function of number of selected residuals.
-figure;
-plot([1-fa; probabilityMaximumIsolation']','LineWidth',2)
+figure(22);
+plot([1-fa, probabilityMaximumIsolation], 'LineWidth', 2)
 l = legend(data.modes{1:end});
 set(l, 'Interpreter', 'none')
-xlabel('Number of selected residuals','FontWeight', 'bold')
-ylabel('Probability of maximum isolation performance','FontWeight', 'bold')
+xlabel('Number of selected residuals', 'FontWeight', 'bold')
+ylabel('Probability of maximum isolation performance', 'FontWeight', 'bold')
 grid on
 clear probabilityMaximumIsolation fa l
 
 %% Plot Fig 13: Consistency-based diagnosis results obtained when running different numbers k of most important residuals.
 % Plots the cases when using 10, 12, 26, and 27 residuals:
-selected_no_of_residauls = [10 12 26 27];
-figure;
-C ={};
-for i = 1:length(selected_no_of_residauls) 
-    subplot(2,2,i);
-    selected_residauls = residualRanking(1:selected_no_of_residauls(i));
-    C{i} = CBDPlots(data,model,selected_residauls);
-    title(['No of tests: ' num2str(selected_no_of_residauls(i))]); 
+selected_no_of_residuals = [10 12 26 27];
+figure(23);
+C = cell(length(selected_no_of_residuals), 1);
+for i = 1:length(selected_no_of_residuals) 
+    subplot(2, 2, i);
+    selected_residuals = residualRanking(1:selected_no_of_residuals(i));
+    C{i} = CBDPlots(data, model, selected_residuals);
+    title(['No of tests: ' num2str(selected_no_of_residuals(i))]); 
 end
 clear C i 
 
 %% Plot structural isolability using the k best residuals. Compare with previous plot.
-figure;
-for i = 1:length(selected_no_of_residauls)
-    subplot(2,2,i)
-    selected_residauls = residualRanking(1:selected_no_of_residauls(i));
-    model.IsolabilityAnalysisFSM(data.fsm(selected_residauls,[2:8]),'permute',0); % full str isolability
-    title(['No of tests: ' num2str(selected_no_of_residauls(i))]) 
+figure(24);
+for i = 1:length(selected_no_of_residuals)
+    subplot(2, 2, i)
+    selected_residuals = residualRanking(1:selected_no_of_residuals(i));
+    model.IsolabilityAnalysisFSM(data.fsm(selected_residuals, (2:8)), ...
+                                 'permute', 0); % full str isolability
+    title(['No of tests: ' num2str(selected_no_of_residuals(i))]) 
 end
-%clear selected_no_of_residauls selected_residauls residualRanking i ans
+%clear selected_no_of_residuals selected_residuals residualRanking i ans
 

code/src/SingleFaultIsolation.m

+function isol=SingleFaultIsolation( res, FSM )
+  N = size(res,1);
+  nf = size(FSM,2);
+  
+  isol = zeros(N,nf + 1);
+  for k=1:N
+    alarms = abs(res(k,:))>=1;
+    if any(alarms)
+      isol(k,:) = [false all(FSM(alarms,:),1)];
+    else
+      isol(k,:) = [true boolean(zeros(1,nf))];
+    end
+  end
+end
\ No newline at end of file