MorrisScreeningTestFunction.m

%% Tutorial: Performing Morris-screening on Ishigami function
% By Adem R.N. Aouichaoui @DTU
% revised by GSI / 8.8.2024

% Initializing
clear;close all; clc         
format short

% Model: Ishigami function [https://www.sfu.ca/~ssurjano/ishigami.html]
f = @(x) sin(x(:,1)) + 7.*sin(x(:,2)).^2 + 0.1.*x(:,3).^4.*sin(x(:,1));  

% Uniform Input Space
pars = {'x1','x2','x3'}; % input parameter names
lbs  = -pi.*ones(1,3);   % lower bounds of input parameters
ubs  =  pi.*ones(1,3);   % upper bounds of input parameters

%% Step 1: perform Morris sampling
disp('Step 1: perform Morris sampling')

% Morris sampling parameters
k = length(pars) ; % no of parameters or factors
p = 4 ; % number of levels {4,6,8}
dt = p/(2*(p-1)) ; % perturbation factor .
r = 50; % number of repetion for calculating the EEi, e.g. 4 - 15

%Morris sampling will produce discrete uniform probabilities for each
% factor.
X = morris(p,dt,k,r);
Xmean = mean(X) ;

% from uniform distribution [0 1] to real values
Xval = ones(r*(k+1),1)*lbs + (ones(r*(k+1),1)*(ubs-lbs)).*X ;
%a=icdf('uniform',X(:,1),lbs(1),ubs(1));

% Plot Morris sampling results for the visualization
figure(1)
[~,ax,~,~] = plotmatrix(X) ;
set(ax,'FontSize',14,'FontWeight','bold')
for i=1:k
ylabel(ax(i,1),pars(i))
xlabel(ax(k,i),pars(i))
end
title(['Morris Sampling with r=', num2str(r),', p=',num2str(p),' levels,',' and \Delta =',num2str(dt,'%5.2f'),': unit range'])
saveas(1,'MorrisSamples in P','tiff')

figure(2)
[h,ax,bax,P] = plotmatrix(Xval) ;
set(ax,'FontSize',14,'FontWeight','bold')
for i=1:k
ylabel(ax(i,1),pars(i))
xlabel(ax(k,i),pars(i))
end
title(['Morris Sampling with r=', num2str(r),', p=',num2str(p),' levels,',' and \Delta =',num2str(dt,'%5.2f'),': real values'])
saveas(2,'MorrisSamples','tiff')


%% step 2: perform simulations with Morris Samples
disp('step 2: perform simulations with Morris Samples')

% Evlauate / simulations of function f
y = f(Xval);

%% Step 3: visualize the simulation results
disp('step 3: visualize the simulation results')

fs=12;
figure % Create a new figure window
plot(y,'b','LineWidth',1.5)
xlabel('iteration','FontSize',fs,'FontWeight','bold') 
ylabel('f - Ishigami ','FontSize',fs,'FontWeight','bold') 
set(gca,'LineWidth',2,'FontSize',12,'FontWeight','bold') 
title('Results from simulating the Ishigami function')

saveas(3,'visualizethesimulations','tiff')

%% Step Calculate EEi
disp('step 4: compute Elementary Effects, EEi')
sig_y = std(y);
sig_x = std(Xval);


[n, k] = size(Xval) ;
mr = n / (k+1); % repetition
[mm,ll] = size(y) ;

for i=1:r
    for j=1:k
        m1 = (k+1)*(i-1);
        r1 = m1 + j ;
        ix = find (X(r1,:) - X(r1+1,:)) ; % find the non-zero value
        if length(ix) > 1
            warning('there is more than one factor changed')
            return
        end
        m2 = k*(i-1)+j; % quality check: this is to track which params changed
        idx(m2) = ix ;
        dtheta = Xval(r1+1,ix) - Xval(r1,ix) ;
        EE(i,ix) = (y(r1+1,1) - y(r1,1))/dtheta * (sig_x(ix)/sig_y(1))  ; % SEEi = EEi* sigx/sigy where EEi = (Ytheta - Ytheta+dtheta) / dtheta *
    end
end

%Calculate mean and std from Fi = EEi
mu  = mean(EE);
sig = std(EE);
mua =mean(abs(EE));

%% Step Plots results 
disp('step 5: Plot Elementary Effects, EEi')
figure
subplot(3,1,1)
hist(EE(:,1))
xlabel(['EEi of ',pars{1}])
ylabel('Frequency')
title('Ishigami')
subplot(3,1,2)
hist(EE(:,2))
xlabel(['EEi of ',pars{2}])
ylabel('Frequency')
subplot(3,1,3)
hist(EE(:,3))
xlabel(['EEi of ',pars{3}])
ylabel('Frequency')

saveas(4,'HistogramEEi','tiff')

%% Plotting
 figure

for j=1:k
    sem(j) = 2*sig(j)/sqrt(r);
    plot(mu(j),sig(j),'ko',2*sem(j),sig(j),'k',-2*sem(j),sig(j),'k')
    hold on
    text(mu(j),sig(j),[pars{j},'=',num2str(j)])
    ylim([0 max(sig)])
end
    title('Ishigami function')
    xlabel('mean, \mu')
    ylabel('stdev, \sigma')

  
 figure 
 bar(pars,mua)   
ylabel('absolute mean')

saveas(6,'ImportanceRanking','tiff')

%% plot of muA
T = table;
T.mua =mua;
T