load dryer2;
data_n = length(y2);
output = y2;
input = [[0; y2(1:data_n-1)] ...
        [0; 0; y2(1:data_n-2)] ...
        [0; 0; 0; y2(1:data_n-3)] ...
        [0; 0; 0; 0; y2(1:data_n-4)] ...
        [0; u2(1:data_n-1)] ...
        [0; 0; u2(1:data_n-2)] ...
        [0; 0; 0; u2(1:data_n-3)] ...
        [0; 0; 0; 0; u2(1:data_n-4)] ...
        [0; 0; 0; 0; 0; u2(1:data_n-5)] ...
        [0; 0; 0; 0; 0; 0; u2(1:data_n-6)]];
data = [input output];
data(1:6, :) = [];
input_name = str2mat('y(k-1)','y(k-2)','y(k-3)','y(k-4)','u(k-1)','u(k-2)','u(k-3)','u(k-4)','u(k-5)','u(k-6)');
trn_data_n = 300;
index = 1:100;
figure;
subplot(2,1,1);
plot(index, y2(index), '-', index, y2(index), 'o');
ylabel('y(k)');
subplot(2,1,2);
plot(index, u2(index), '-', index, u2(index), 'o');
ylabel('u(k)');





figure;
group1 = [1 2 3 4]; % y(k-1), y(k-2), y(k-3), y(k-4)
group2 = [1 2 3 4]; % y(k-1), y(k-2), y(k-3), y(k-4)
group3 = [5 6 7 8 9 10]; % u(k-1) through y(k-6)

anfis_n = 6*length(group3);
index = zeros(anfis_n, 3);
trn_error = zeros(anfis_n, 1);
chk_error = zeros(anfis_n, 1);
% ======= Training options 
mf_n = 2;
mf_type = 'gbellmf';
epoch_n = 1;
ss = 0.1;
ss_dec_rate = 0.5;
ss_inc_rate = 1.5;
% ====== Train ANFIS with different input variables
fprintf('\nTrain %d ANFIS models, each with 3 inputs selected from 10 candidates...\n\n',...
    anfis_n);
model = 1;
for i=1:length(group1),
    for j=i+1:length(group2),
        for k=1:length(group3),
            in1 = deblank(input_name(group1(i), :));
            in2 = deblank(input_name(group2(j), :));
            in3 = deblank(input_name(group3(k), :));
            index(model, :) = [group1(i) group2(j) group3(k)];
            trn_data = data(1:trn_data_n, [group1(i) group2(j) group3(k) size(data,2)]);
            chk_data = data(trn_data_n+1:trn_data_n+300, [group1(i) group2(j) group3(k) size(data,2)]);
            in_fismat = genfis1(trn_data, mf_n, mf_type);
            [trn_out_fismat t_err step_size chk_out_fismat c_err] = ...
                anfis(trn_data, in_fismat, ...
                [epoch_n nan ss ss_dec_rate ss_inc_rate], ...
                [0 0 0 0], chk_data, 1);
            trn_error(model) = min(t_err);
            chk_error(model) = min(c_err);
            fprintf('ANFIS model = %d: %s %s %s', model, in1, in2, in3);
            fprintf(' --> trn=%.4f,', trn_error(model));
            fprintf(' chk=%.4f', chk_error(model));
            fprintf('\n');
            model = model+1;
        end
    end
end

% ====== Reordering according to training error
[a b] = sort(trn_error);
b = flipud(b);  % List according to decreasing trn error
trn_error = trn_error(b);
chk_error = chk_error(b);
index = index(b, :);

% ====== Display training and checking errors
x = (1:anfis_n)';
subplot(2,1,1);
plot(x, trn_error, '-', x, chk_error, '-', ...
    x, trn_error, 'o', x, chk_error, '*');
tmp = x(:, ones(1, 3))';
X = tmp(:);
tmp = [zeros(anfis_n, 1) max(trn_error, chk_error) nan*ones(anfis_n, 1)]';
Y = tmp(:);
hold on; 
plot(X, Y, 'g'); 
hold off;
axis([1 anfis_n -inf inf]);
set(gca, 'xticklabel', []);

% ====== Add text of input variables
for k = 1:anfis_n,
    text(x(k), 0, ...
        [input_name(index(k,1), :) ' ' ...
            input_name(index(k,2), :) ' ' ...
            input_name(index(k,3), :)]);
end
h = findobj(gcf, 'type', 'text');
set(h, 'rot', 90, 'fontsize', 11, 'hori', 'right');
drawnow

% ====== Generate input_index for bjtrain.m
[a b] = min(trn_error);
input_index = index(b,:);
title('Training (Circles) and Checking (Asterisks) Errors');
ylabel('RMSE');


figure;
ss = 0.01;
ss_dec_rate = 0.5;
ss_inc_rate = 1.5;

trn_data = data(1:trn_data_n, [input_index, size(data,2)]);
chk_data = data(trn_data_n+1:600, [input_index, size(data,2)]);

% generate FIS matrix
in_fismat = genfis1(trn_data);

[trn_out_fismat trn_error step_size chk_out_fismat chk_error] = ...
    anfis(trn_data, in_fismat, [1 nan ss ss_dec_rate ss_inc_rate], ...
    nan, chk_data, 1);

y_hat = evalfis(data(1:600,input_index), chk_out_fismat);
%figure;
subplot(2,1,1);
index = 1:trn_data_n;
plot(index, data(index, size(data,2)), '-', ...
    index, y_hat(index), '.');
rmse = norm(y_hat(index)-data(index,size(data,2)))/sqrt(length(index));
title(['Training Data (Solid Line) and ANFIS Prediction (Dots) with RMSE = ' num2str(rmse)]);
xlabel('Time Index'); ylabel('');

subplot(2,1,2);
index = trn_data_n+1:600;
plot(index, data(index, size(data,2)), '-', index, y_hat(index), '.');
rmse = norm(y_hat(index)-data(index,size(data,2)))/sqrt(length(index));
title(['Checking Data (Solid Line) and ANFIS Prediction (Dots) with RMSE = ' num2str(rmse)]);
xlabel('Time Index'); ylabel('');