% Sample Taguchi design array experiment (goes along with a class example):
% Barukyah Shaparenko and John M. Cimbala

clear all
close all

% parameter a = spring constant
% parameter b = weight of bullet
% parameter c = diameter of bullet
% X = response to be optimized, a function of a, b, and c

% First, input the data. We provide three different ways to do this.
% Comment out the method that you DON'T want to use.

% Method 1: Explicit data values entered manually
%{
a = [1; 1; 2; 2; 3; 3; 4; 4]; 
b = [1; 2; 3; 4; 1; 2; 3; 4]; 
c = [1; 4; 3; 2; 2; 3; 1; 4]; 
% Levels indicated as 1, 2, 3, and 4 for simplicity.

X = [5.5; 6.2; 6.0; 7.3; 7.6; 8.1; 9.7; 7.2];
%}

% ---OR---

% Method 2: Import from Excel
% NOTE: You must save the Excel file in the same folder as the Matlab m-file
% Matlab will search the working directory for a file with this name.
% Use xlsread or xlsxread

% xls_filename = 'Example_Taguchi_array.xls';
% sheet_number = 1;
% data_range_1 = 'C11:C18';
% data_range_2 = 'D11:D18';
% data_range_3 = 'E11:E18';
% data_range_4 = 'F11:F18';
% a = xlsread(xls_filename,sheet_number,data_range_1);
% b = xlsread(xls_filename,sheet_number,data_range_2);
% c = xlsread(xls_filename,sheet_number,data_range_3);
% X = xlsread(xls_filename,sheet_number,data_range_4); 


% ---OR---

% Method 3: Import from a text file
% NOTE: You must save the text file in the same folder as the Matlab m-file
% Matlab will search the working directory for a file with this name.
fid1=fopen('Example_Taguchi_array.txt');
num_lines = 8;
%fseek(fid1,1,'bof'); % Use if the data file contains a header.
read_data = textscan(fid1,'%f %f %f %f',num_lines);    
data = [read_data{1,:}];

a = data(:,1);
b = data(:,2);
c = data(:,3);
X = data(:,4);

% Find and sort unique values of the parameters:
%{
The A = unique(sort(a)) is to find the unique levels for the parameter a.
We use it to compute the level averages easily by averaging only the values
of X that have parameter a equal to each of the values in A.
Finding the unique values is necessary in order to have something to
compare with a. This is equivalent to automatically picking the values to
average as we do manually in Excel.
This command is especially useful when you do not have just level numbers,
but actual parameter values at these levels.
%}
A = unique(sort(a));
B = unique(sort(b));
C = unique(sort(c));

%Create arrays containing all zeroes to start with:
X_a_bar = zeros(length(A),1);
X_b_bar = zeros(length(B),1);
X_c_bar = zeros(length(C),1);

% Calculate the level averages:
for i = 1:length(A)
    X_a_bar(i) = mean(X(a==A(i)));
    X_b_bar(i) = mean(X(b==B(i)));
    X_c_bar(i) = mean(X(c==C(i)));
end

% Plot the level averages:
scrsz = get(0,'ScreenSize');
figure('Position',[scrsz(3)/10 scrsz(4)*3/8 scrsz(3)*4/5 scrsz(4)/2])

subplot(1,3,1)
hold on;
a_smooth = interp1(A,X_a_bar,min(A):(max(A)-min(A))/25:max(A),'cubic');
plot(min(A):(max(A)-min(A))/25:max(A),a_smooth,'r-','LineWidth',2)
plot(A,X_a_bar,'o','MarkerFaceColor','b','MarkerEdgeColor','b','MarkerSize',9)
xlim([1 4])
ylim([5 9])
xlabel('level for parameter \ita','FontWeight','Bold')
ylabel('level average for \itX_a','FontWeight','Bold')
pbaspect([1 1 1])
grid on;
hold off;

subplot(1,3,2)
hold on;
b_smooth = interp1(B,X_b_bar,min(B):(max(B)-min(B))/25:max(B),'cubic');
plot(min(B):(max(B)-min(B))/25:max(B),b_smooth,'r-','LineWidth',2)
plot(B,X_b_bar,'o','MarkerFaceColor','b','MarkerEdgeColor','b','MarkerSize',9)
xlim([1 4])
ylim([5 9])
xlabel('level for parameter \itb','FontWeight','Bold')
ylabel('level average for \itX_b','FontWeight','Bold')
pbaspect([1 1 1])
grid on;
hold off;

subplot(1,3,3)
hold on;
c_smooth = interp1(C,X_c_bar,min(C):(max(C)-min(C))/25:max(C),'cubic');
plot(min(C):(max(C)-min(C))/25:max(C),c_smooth,'r-','LineWidth',2)
plot(C,X_c_bar,'o','MarkerFaceColor','b','MarkerEdgeColor','b','MarkerSize',9)
xlim([1 4])
ylim([5 9])
xlabel('level for parameter \itc','FontWeight','Bold')
ylabel('level average for \itX_c','FontWeight','Bold')
pbaspect([1 1 1])
grid on;
hold off;

% Print the level averages:
fprintf('\n----------------------\n');
fprintf('\t\tLevel averages:\n');
fprintf('level\tname\tvalue\t\tname\tvalue\t\tname\tvalue\n');
for i = 1:length(A)
    fprintf('%g\t\tX_a%g =\t%g\t\tX_b%g =\t%g\t\tX_c%g =\t%g\t\t\n',i,i,X_a_bar(i),i,X_b_bar(i),i,X_c_bar(i))
end
