% DEMOHILOMOTDOECLASS demonstrates the functioning of the HILOMOTDOE toolbox for a nonlinear
% static regression example.

% HILOMOTDOE Nonlinear System Identification Toolbox for Design of Experiments
% Benjamin Hartmann, 25-November-2011
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2011 by Prof. Dr.-Ing. Oliver Nelles

clear all
clear java
clear
close all
clc

% Add path of HILOMOT class files
% restoredefaultpath
% hilomotDirectory = fileparts(which(mfilename));
% addpath(hilomotDirectory);

% Initialize hilomot object
design = hilomotDoE;

% User defined options
design.numberOfInitialPoints = 2;
design.smoothnessQuery = 1.0;

% Option for implementation of variable increments
% design.variableIncrements = [0 0.02 1; 0 0.03 1];

% Example selection 1, ..., 13
global example

example = 4;

% Standard deviation of white noise
noiseLevel = 0.01;


% [u1g u2g] = meshgrid(linspace(0,1,30), linspace(0,1,30));
% u1Plot = u1g(:);
% u2Plot = u2g(:);
% yPlot = demoFunctions([u1Plot u2Plot],example) + noiseLevel*rand(size(u1Plot,1),1);
%
% figure(24)
% surf(u1g,u2g,reshape(yPlot,30,30))
% set(gca,'fontsize',13)
% title('Process that has to be measured','fontsize',13)
% xlabel('u_1','fontsize',13)
% ylabel('u_2','fontsize',13)
% zlabel('y','fontsize',13)
% pause(1)
% close


% Generate candidate data set
% load demoCandidatePoints.mat
% design.candidatePoints = candidatePoints;
% u1 = design.candidatePoints(:,1);
% u2 = design.candidatePoints(:,2);
% NCand = size(u1,1);

% [u1gx u2gx] = meshgrid(linspace(0,1,50), linspace(0,1,50));
% u1 = u1gx(:);
% u2 = u2gx(:);
% design.candidatePoints = [u1 u2];
% u1 = u1Plot;
% u2 = u2Plot;
% NCand = size(u1,1);

u = linspace(0,1,100)';
design.candidatePoints = [u];
NCand = size(u,1);


% Simulate measurement with example function:
% y = demoFunctions([u1 u2], example) + noiseLevel*randn(size(u1,1),1);
rng(7,'twister')
y = demoFunctions([u], example);
y = y + noiseLevel*y.*randn(size(y,1),1);

% Initialization phase
[query] = design.DoEUpdate;

% Simulate measurement
yQuery = demoFunctions(query, example) + noiseLevel*randn(size(query,1),1);
yQuery = yQuery + noiseLevel*yQuery.*randn(size(yQuery,1),1);

% Generate data set: data = [u_ist y_ist u_soll]
design.inputTarget   = [design.inputTarget; query];
design.inputMeasured = [design.inputMeasured; query];
design.output        = [design.output; yQuery];

% Enable Active Learning procedure
design.flagSpaceFilling = false;


h = figure(123);
set(h,'Position',[522 446 700 650])


% Active learning phase
for k = 1:10
    
    % Perform active learning
    [query] = design.DoEUpdate;
    
    % Simulate measurement
    yQuery = demoFunctions(query, example) + noiseLevel*randn(size(query,1),1);
    yQuery = yQuery + noiseLevel*yQuery.*randn(size(yQuery,1),1);
    
    
    % Generate data set: data = [u_ist y_ist u_soll]
    design.inputTarget   = [design.inputTarget; query];
    design.inputMeasured = [design.inputMeasured; query];
    design.output        = [design.output; yQuery];
    
    % Example for failed measurement
    %     if k == 11 || k == 23
    %         fprintf('\n\n -> Simulation of failed measurement! \n\n')
    %         pause(0.1)
    %         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %         % Important note for user, if measurement failed:
    %         % Store query as target input. Store measured input and
    %         % output as NaN!
    %         design.inputTarget(end,:)   = query;
    %         design.inputMeasured(end,:) = [NaN NaN];
    %         design.output(end)          = NaN;
    %         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %     end
    
    %     figure(23)
    %     plot(design.inputMeasured(:,1),design.inputMeasured(:,2),'k.','markersize',11,'linewidth',2); hold on
    %     plot(query(:,1),query(:,2),'ro','markersize',11,'linewidth',2); hold off
    %     set(gca,'fontsize',13,'XTick',[0 0.5 1],'YTick',[0 0.5 1])
    %     title('Selected measurement points (black) and new query point (red)','fontsize',13)
    %     xlabel('u_1','fontsize',13)
    %     ylabel('u_2','fontsize',13)
    %     axis([0 1 0 1])
    
    drawnow
    plot(u, y,'color',[0.9 0.9 0.9],'linewidth',2); hold on
    plot(design.inputMeasured,design.output,'k.','markersize',11,'linewidth',2); hold on
    plot(query,yQuery,'ro','markersize',11,'linewidth',2); hold off
    set(gca,'fontsize',13,'XTick',[0 0.5 1],'YTick',[0 0.5 1])
    title('Selected measurement points (black) and new query point (red)','fontsize',13)
    xlabel('u_1','fontsize',13)
    ylabel('u_2','fontsize',13)
    %     axis([0 1 0 1])
    
end


% Plots

% close all
%
% figure(24)
% surf(u1g,u2g,reshape(yPlot,30,30))
% alpha(0.3); hold on
% plot3(design.inputMeasured(:,1),design.inputMeasured(:,2),design.output,'ro','markersize',11,'linewidth',2)
% set(gca,'fontsize',13)
% title('Process and selected measurement points','fontsize',13)
% xlabel('u_1','fontsize',13)
% ylabel('u_2','fontsize',13)
% zlabel('y','fontsize',13)
%
% figure(25)
% contour(u1g,u2g,reshape(yPlot,30,30))
% hold on
% plot3(design.inputMeasured(:,1),design.inputMeasured(:,2),design.output,'ro','markersize',11,'linewidth',2)
% set(gca,'fontsize',13,'XTick',[0 0.5 1],'YTick',[0 0.5 1])
% title('Process and selected measurement points','fontsize',13)
% xlabel('u_1','fontsize',13)
% ylabel('u_2','fontsize',13)





LMN = hilomot;
LMN.xRegressorDegree = 1;
LMN.scaleInput = [];
LMN.input = design.inputMeasured;
LMN.output = design.output;
LMN.LOOCV = true;
LMN.optLOOCV = true;
LMN.smoothness = 2;

LMN = LMN.train(false);

u = linspace(0.001,0.999,100)';
y = demoFunctions(u,example);
yM = LMN.calculateModelOutput(u);

Xneu = [ones(size(u)) u];

Phi = LMN.calculatePhi(u,LMN.leafModels);

% errorbar = 1*LMN.calcErrorbar(u);

minVal = min(y);
maxVal = max(y);


h = figure(123);
set(h,'Position',[522 446 550 650])

clf
% subplot(3,1,[1 2])
hold on
% f = [yM+errorbar; flipdim(yM-errorbar,1)];
% h = fill([u; flipdim(u,1)], f, [7 7 7]/8);
% set(h,'EdgeColor','None')

plot(u, y,'color',[0.9 0.9 0.9],'linewidth',2)
plot(u, yM,'k','linewidth',1)

% plot(u,yM+errorbar,'k--','linewidth',1)
% plot(u,yM-errorbar,'k--','linewidth',1)
plot(LMN.input,LMN.output,'kx','linewidth',2,'markersize',12)
% plot(query,demoFunctions(query, example),'ro','linewidth',2,'markersize',12)
hold off
xlabel('u','fontsize',13)
ylabel('y','fontsize',13)
set(gca,'fontsize',13)
axis([-0.01 1 -0.25 1.25*maxVal])
box on
% subplot(3,1,3)
% hold on
% LMN.plotPartition
% hold off
