%% Comparison of two second degree designs
% (c) jean-marie.furbringer@epfl.ch  sept 2021
% Table of Contents
% The reference model
%     Setup of the reference model
%     Visualization of the reference model
% Strategy 1: composite design
%     Set up of the model matrix
%     Measurement simulation
%     Computing of the coefficients
%     Canonical analysis
%     Visualisation of the second degree fonction
% Strategy 2: Doehlert design
% Strategy 3 : Box-Behnken design
%
%
%% The reference model
%% Setup of the reference model
%  A reference model is chosen to simulate the experiments. It is a quadratic model at 3 dimensions
% y=a_o+a_1x_1+a_2x_2+a_3x_3+a_{12}x_1x_2+a_{13}x_1x_3+a_{23}x_2x_3+a_{11}x_1^2+a_{22}x_2^2+a_{33}x_3^2
% Coefficients of the reference model are chosen at random. This allows a honest comparison between the designs.

%kappa=[.4 -.5 .3 -.2 .7 -.6 -.5 .8 -.9 -.75];
 kappa=rand(1,10);

% An experimental variance is used for the simulation of the experiment:

sigma=0.05;
disp(['sigma = ',num2str(sigma)])

%% Visualization of the reference model
% To visualize the reference model that we attemp to fit, let's draw a slice view and an isovalue representation. To do that it is necessary to
% 1) declare the symbolic variables x_1, x_2, x_3,

syms x_1 x_2 x_3
assume(x_1,'real')
assume(x_2,'real')
assume(x_3,'real')

% 2) build  the model coordinate as symbolic vector ,
% f(x)=\left( \matrix{1 \cr x_1 \cr x_2 \cr x_3 \cr x_1x_2   \cr x_1x_3  \cr x_2x_3   \cr x_1^2   \cr x_2^2  \cr x_3^2  } \right)

f=[1;x_1;x_2;x_3;x_1*x_2;x_1*x_3;x_2*x_3;x_1^2;x_2^2;x_3^2];

% 3) build  the model as the product of the coefficients and the model coordinate,y=\alpha \cdot f(x), ce qui donne:

y=kappa*f

% 4) define the domain,

% domain
min_x=-1;
max_x=1;

% 5) define the number of points per interval

% number of points per dimension
N=21;
disp(['N = ',num2str(N)])

% 6) compute the values of the model on a grid defined with 3 3D-matrices X_1, X_2 and X_3 (one per dimension),

% factor values
[X1,X2,X3]=meshgrid(linspace(min_x,max_x,N));

% position of the planes for the slices

p1=-1;
p2=0.9;
p3=0;

Xs=[p1 p2 p3];

% evaluation of the model
Y_ref=double(subs(y,{x_1,x_2,x_3},{X1,X2,X3}));

% 7) visualize the model with the function viz() defined at the end of this file

viz(X1,X2,X3,Y_ref,Xs)

% This ends the setup of the experimental "reality". Obviously, in the real situation, those informations are not known by the experimenter.
% Now let's start with the experimenter mimicking.
%
%% Strategy 1: composite design
% Let's built a factorial design for 3 factors:

E1=(ff2n(3)-.5)*2;
disp(array2table(E1,"VariableNames",{'X_1','X_2','X_3'}))

% Let's built a star design for 3 factors:

star= [1 0 0
      -1 0 0
       0 1 0
       0 -1 0
       0 0 1
       0 0 -1];

% Ratio of the distance of the star points in relation with the factorial points

alfa=1.215; 
E2=alfa*star;
disp(array2table(E2,"VariableNames",{'X_1','X_2','X_3'}))

% Set up the central composite design with one central experiment

E= [E1;E2;0 0 0];
disp(array2table(E,"VariableNames",{'X_1','X_2','X_3'}))

%
%% Set up of the model matrix
% Model matrix

M=x2fx(E,'quadratic');

% Dispertion matrix D=(X'X)^{-1}

D=inv(M'*M);

% Representation of the model matrix

figure
imagesc(M),colorbar
set(gca, 'Xtick',.5:1:10.5,...
'Ytick',.5:1:15.5,...
'FontSize',16,...
'XTickLabel',{' ' 'I' '1' '2' '3' '12' '13' '23' '11' '22' '33'},...
'YTickLabel',{ ' ' '1' '2' '3' '4' '5' '6' '7' '8' '9' '10' '11' '12' '13' '14' '15'})
colormap(winter(10))
axis ij
axis square
title('Matrix of the model','FontSize',16)
xlabel('Coefficients','FontSize',16)
ylabel('Experiments','FontSize',16)

%  Representation of the dispertion matrix

figure
imagesc(D),colorbar
set(gca, 'Xtick',.5:1:10.5,...
'Ytick',.5:1:10.5,...
'FontSize',16,...
'XTickLabel',{' ' 'I' '1' '2' '3' '12' '13' '23' '11' '22' '33'},...
'YTickLabel',{ ' ' 'I' '1' '2' '3' '12' '13' '23' '11' '22' '33'})
%c with olormap(pink(10)) 
axis ij
axis square
title('Matrix of dispersion','FontSize',16) 
xlabel('Coefficients','FontSize',16)
ylabel('Coefficients','FontSize',16)

% The structure of the dispersion matrix is typical from a quadratic model orthogonalyty for the first main effects and the interactions and the correlation with the constant for the quadratic terms. The computing of the alias matrix between the linear model with interactions and the quadratic terms shows that if the quadratic terms are not evaluated their values will biased the constant with 73 percent of the quadratic curvature.

M1=M(:,1:7);
M2=M(:,8:10);
M_Alias=(M1'*M1)\(M1'*M2);
disp(M_Alias)

%
%% Measurement simulation
% The measurement is simulated by adding a random error to the real model
% Y=X \kappa + \epsilon_i

Nexp=size(M,1); % number of experiments
 Y = M*kappa'.*(ones(Nexp,1)+randn(Nexp,1)*sigma);

% Graphical representation of the experimental results

figure
stem(Y,"filled")
set(gca,'FontSize',14)
title('Experimental results','FontSize',16)
xlabel('Runs','FontSize',16)
grid on

%
%% Computing of the coefficients
% The coefficients are computed by a least square fit algorithm
% $\hat{\beta}=(X'X)^{-1} X' Y$

beta=D*M'*Y;

% and then they can be compared with the original ones

figure
bar([beta/beta(1), kappa'/kappa(1)])
set(gca,'FontSize',16,...
    'XTickLabel',{'a_o' 'a_1' 'a_2' 'a_3' ...
    'a_{12}' 'a_{13}' 'a_{23}' 'a_{11}' 'a_{22}' 'a_{33}'})
title('Model','FontSize',16)
xlabel('coefficients','FontSize',16)
ylabel('aij/ao','FontSize',16)
legend("Estimated", "True")
grid on
mdl_1=fitlm(E,Y,"quadratic")
anova(mdl_1,"summary",2)

%
%% Canonical analysis
% The first step consists in expressing the model in terms of constant a_o, linear vector \vec{a} and quadratic matrix A:

% Routine belonging to LISA library
ao=beta(1);
a=beta(2:4);
A=[beta(8) beta(5)/2 beta(6)/2
beta(5)/2 beta(9) beta(7)/2
beta(6)/2 beta(7)/2 beta(10)];

% The quadratic matric is:

disp(A)

% Then is is possible to use the user-defined function analyse_canonique() provided in Moodle which computes the fix point and the eigen vectors \vec{V}_i and eigen values \lambda_i :

[Xs,Ys,V,lamdda]=analyse_canonique(ao,a,A);

% The fix point X_s is

disp(array2table(Xs',"VariableNames",{'x_1','x_2','x_3'}))

% that we want to check to know if it is placed inside the domain.
% The eigen vectors are:

disp(array2table(V,"VariableNames",{'V_1','V_2','V_3'}))

% The eigen values are

disp(array2table(diag(lamdda)',"VariableNames",{'lambda_1','lambda_2','lambda_3'}))

% Now it is possible to draw the position of the new reference frame (in red the experimental points are also indicated)


% Setup of the new referential frame
base=[Xs,Xs+V(:,1),Xs,Xs+V(:,2),Xs,Xs+V(:,3)];

figure
plot3(Xs(1),Xs(2),Xs(3),'*',base(1,:),base(2,:),base(3,:),'-',...
"LineWidth",2,"Color","blue")
hold on
plot3(E(:,1),E(:,2),E(:,3),'ok',...
"MarkerFaceColor","red")
hold off
box on
xlabel('x_1')
ylabel('x_2')
zlabel('x_3')
grid on

% coordinates of the labels of the eigen vectors
dl=1.2;
coord=[Xs-0.2*V(:,1),Xs+dl*V(:,1),Xs+dl*V(:,2),Xs+dl*V(:,3)]';
text(coord(:,1),coord(:,2),coord(:,3),["X_s","V_1","V_2", "V_3"])

%%
%% Visualisation of the second degree fonction
% Compute the value of the fonction in a parallelepiped around the origin

y_hat=beta'*f

% 4) Computing values of the model on a grid

% number of points per dimension
N=21;

% domain
min_x=-1;
max_x=1;

% factor values
[X1,X2,X3]=meshgrid(linspace(min_x,max_x,N));
Xs=zeros(1,3);

% evaluation of the model
Y_hat=double(subs(y_hat,{x_1,x_2,x_3},{X1,X2,X3}));

% 5) Visualization

viz(X1,X2,X3,Y_hat,Xs)

% 6) analyse the residuals

plotResiduals(mdl_1)
plotAdded(mdl_1)
plotAdjustedResponse(mdl_1,"x1")
plotAdjustedResponse(mdl_1,"x2")
plotAdjustedResponse(mdl_1,"x3")

%
%
%% Strategy 2: Doehlert design
% Let's see how the same situation results when using a Doehlert design

E_doehlert=doehlert(3);
disp(array2table(E_doehlert,"VariableNames",{'X_1','X_2','X_3'}))
M_doehlert=x2fx(E_doehlert,"quadratic");
Nexp=size(M_doehlert,1); % number of experiments
Y_doehlert = M_doehlert*kappa'.*(ones(Nexp,1)+randn(Nexp,1)*sigma);

% Graphical representation of the experimental results

figure
stem(Y_doehlert)
set(gca,'FontSize',14)
title('Experimental results (Doehlert)','FontSize',16)
xlabel('Runs','FontSize',16)
grid on
mdl_2=fitlm(E_doehlert,Y_doehlert,"quadratic")
anova(mdl_2,"summary",2)
figure
gamma=mdl_2.Coefficients.Estimate;
bar([gamma/gamma(1), beta/beta(1), kappa'/kappa(1)])
set(gca,'FontSize',16,...
    'XTickLabel',{'a_o' 'a_1' 'a_2' 'a_3' ...
    'a_{12}' 'a_{13}' 'a_{23}' 'a_{11}' 'a_{22}' 'a_{33}'})
title('Model','FontSize',16)
xlabel('coefficients','FontSize',16)
ylabel('aij/ao','FontSize',16)
legend("Doehlert","composite", "True")
grid on
plotResiduals(mdl_2)
plotAdded(mdl_2)
plotAdjustedResponse(mdl_2,"x1")
plotAdjustedResponse(mdl_2,"x2")
plotAdjustedResponse(mdl_2,"x3")

%
%% Strategy 3 : Box-Behnken design
% Let's see how the same situation results when using a Doehlert design

E_BB=bbdesign(3,"center",3);
disp(array2table(E_BB,"VariableNames",{'X_1','X_2','X_3'}))
M_BB=x2fx(E_BB,"quadratic");
Nexp=size(M_BB,1); % number of experiments
Y_BB = M_BB*kappa'.*(ones(Nexp,1)+randn(Nexp,1)*sigma);

% Graphical representation of the experimental results

figure
stem(Y_BB)
set(gca,'FontSize',14)
title('Experimental results (Box-Behnken)','FontSize',16)
xlabel('Runs','FontSize',16)
grid on


mdl_3=fitlm(E_BB,Y_BB,"quadratic")
anova(mdl_3,"summary",2)


figure
iota=mdl_3.Coefficients.Estimate;
bar([iota/iota(1),gamma/gamma(1), beta/beta(1), kappa'/kappa(1)])
set(gca,'FontSize',16,...
    'XTickLabel',{'a_o' 'a_1' 'a_2' 'a_3' ...
    'a_{12}' 'a_{13}' 'a_{23}' 'a_{11}' 'a_{22}' 'a_{33}'})
title('Model','FontSize',16)
xlabel('coefficients','FontSize',16)
ylabel('aij/ao','FontSize',16)
legend("BB","Doehlert","composite", "True")
grid on


plotResiduals(mdl_3)
plotAdded(mdl_3)
plotAdjustedResponse(mdl_3,"x1")
plotAdjustedResponse(mdl_3,"x2")
plotAdjustedResponse(mdl_3,"x3")



%

function viz(X1,X2,X3,Y,Xs)
% viz(ValuesMatrix,X_s)
% Y: matrix with the values of the function
% Xi: factors values

% The 3 planes go by x_s.

slice(X1,X2,X3,Y,Xs(1),Xs(2),Xs(3))
title('Values of the function in the experimental space')
xlabel('x1')
ylabel('x2')
zlabel('x3')
colorbar
view([-24.70 40.40])

% Isosurface drawing

% save the colormap for defining automatiquely the color of the surface
map=colormap;

% Compute the limit of Y
Ymin=min(min(min(Y)));
Ymax=max(max(max(Y)));

% number of surfaces
Ns=7;

% define the levels of the isosurfaces
Levels= linspace(Ymin,Ymax,Ns); 

figure


for k=2:6
    % compute the iso surface and built the patches
    p=patch(isosurface(X1,X2,X3,Y,Levels(k)));
    isonormals(X1,X2,X3,Y,p)

    % define the color of the surface from the actual color map
    p.FaceColor=map(round(120+120/7*k),1:3); 
    p.EdgeColor='none';
    hold on
end

colorbar

% close the drawing loop
hold off

% define a cubique zone for the 3D la representation
daspect([1 1 1]) 
%view(3);
view([-48 37])
axis tight
camlight
lighting gouraud % lighting adéquat pour surface courbée
view([-48.30 8.40])
title('iso surface de la fonction y')
xlabel('x_1')
ylabel('x_2')
zlabel('x_3')


end