% lolimotDemo6 Demo 6. Estimate a model given an initial structure

% LoLiMoT - Nonlinear System Identification Toolbox
% Torsten Fischer, 20-October-2011
% Institute of Mechanics & Automatic Control, University of Siegen, Germany
% Copyright (c) 2012 by Prof. Dr.-Ing. Oliver Nelles


% Training
lolimotObject = lolimot;                   % Generate an empty net and data set structure

% lolimotObject.scaleInput = []; % do not scale the inputs

[u1g u2g] = meshgrid(linspace(0,1,25), linspace(0,1,25));
u1 = u1g(:);
u2 = u2g(:);
y = 1./(0.1+u1) + (2*u2).^2;
lolimotObject.input = [u1 u2];
lolimotObject.output = y;


% Define initial structure
lolimotObject.leafModels = true(1,6);
lowerLeftCorner = [0 0; 0.0625 0; 0.125 0; 0.25 0; 0.5 0.5; 0.5 0];

%lowerLeftCorner = bsxfun(@minus,lowerLeftCorner,lolimotObject.scaleInput.parameter.mean);
%lowerLeftCorner = bsxfun(@rdivide,lowerLeftCorner,lolimotObject.scaleInput.parameter.std);

upperRightCorner = [0.0625 1; 0.125 1; 0.25 1; 0.5 1; 1 1; 1 0.5];
%upperRightCorner = bsxfun(@minus,upperRightCorner,lolimotObject.scaleInput.parameter.mean);
%upperRightCorner = bsxfun(@rdivide,upperRightCorner,lolimotObject.scaleInput.parameter.std);

for k = 1:size(lowerLeftCorner,1)
    lolimotObject.localModels = [lolimotObject.localModels gaussianOrthoLocalModel(lowerLeftCorner(k,:),upperRightCorner(k,:),lolimotObject.smoothness)];
end


lolimotObject.maxNumberOfLM = 15;                           % Termination criterion for maximal number of LLMs
lolimotObject.minError = 0.05;                              % Termination criterion for minimal error
lolimotObject.kStepPrediction = 0;                          % Static model
lolimotObject.history.displayMode = true;                   % display information

% Train net
lolimotObject = lolimotObject.train;

% Find best performance/complexity trade-off
[~, idxBest] = min(lolimotObject.history.penaltyLossFunction);
lolimotObject.leafModels = lolimotObject.history.leafModelIter{idxBest};


% Generalization
[u1g u2g] = meshgrid(linspace(0,1,30), linspace(0,1,30));
u1G = u1g(:);
u2G = u2g(:);
yG = 1./(0.1+u1G) + (2*u2G).^2;
uG = [u1G u2G];

% Simulate net
yGModel = calculateModelOutput(lolimotObject, uG, yG);
JG = calcGlobalLossFunction(lolimotObject ,yG, yGModel);

figure
lolimotObject.plotModel

figure
lolimotObject.plotPartition

figure
lolimotObject.plotMSFValue