function [y name] = exampleFunctions(u,bsp)
% EXAMPLEFUNCTIONS calculates the output for a defined example for a 
% defined input space.
%
% y = polymotExamples(u,bsp)
%
% OUTPUT:
%   u:    (N x n)   Input Matrix u; N data samples - n dimensions.
%
%   bsp:  scalar    Indicates which example to choose.
%
% INPUT:
%   y:    (N x 1)   Output vektor for N data samples.
%
%
% PolyMoT Nonlinear System Identification Toolbox, Version 2.1
% Torsten Fischer, 08-Dec-2011
% Automatic Control - Mechatronics, University of Siegen, Germany
% Copyright (c) 2011 by Torsten Fischer

%% Get the right dimension and example
n = size(u,2); % number of input dimensions

if n==1
    switch bsp
        case 1 % hyperbel
            y = 1 ./ (0.1 + u);
            name = 'hyperbel';
        case 2 % sin superposition with straight line
            y = sin(2*pi*u) + 2.*u;
            name = 'sin superposition with straight line';
        case 3 % triangle (hard switch of grad)
            y = zeros(length(u),1);
            b = 0.55;
            for k = 1:length(u)
                if u(k) < b
                    y(k) = 1/b*u(k);
                else
                    y(k) = -1/(1-b) * u(k) + 1/(1-b);
                end
            end
            name = 'triangle (hard switch of grad)';
        case 4 % triangle (hard switch of grad)
            y = zeros(length(u),1);
            b = 0.25;
            if b >= 0.5
                error('choosen value not allowed!')
            end
            for k = 1:length(u)
                if u(k) < b
                    y(k) = 1/b*u(k);
                elseif u(k) < 1-b
                    y(k) = -1/(1-2*b) * u(k) + (1-b)/(1-2*b);
                else
                    y(k) = 1/b*u(k) + 1 - 1/b;
                end
            end
        otherwise
            error('invalid example')
    end
elseif n==2
    u1=u(:,1); u2=u(:,2);
    switch bsp
        case 1 % hyperbel (symmetric)
            y = 1./(0.1+0.5*(1-u1)+0.5*(1-u2));
            name = 'hyperbel (symmetric)';
        case 2 % hyperbel (not symmetric)
            y = 1./(0.1+0.3*(1-u1)+0.7*(1-u2));
            name = 'hyperbel (not symmetric)';
        case 14 % quad-root-combination
            y = 10*u1.^2.*sqrt(u2);
            name = 'quad-root-combination';
        case 4 % hyperbel (symmetric; very steep)
            y = 0.01./((1.1-u1).*(1.1-u2));
            name = 'hyperbel (symmetric; very steep)';
        case 5 % cup-shaped; very smooth
            y = 2-2*exp(-0.5*((u1-0.25).^2+(u2-0.25).^2));
            name = 'cup-shaped; very smooth';
        case 6 % similar to a moutain ridge
            y = 2*exp(8*((u1-0.5).^2+(u2-0.5).^2))./(exp(8*((u1-0.2).^2+(u2-0.7).^2))+exp(8*((u1-0.7).^2+(u2-0.2).^2)));
        case 7 % sinus combine with plain
            y = 2*sin(2*pi*u1).*(1-u2);
        case 8 % Gaussian (middle)
            y = exp(-50*((u1-0.5).^2+(u2-0.5).^2));
            name = 'Gaussian (middle)';
        case 9 % chrip signal combined with plain
            y = 0.25*cos(3*pi/2*u1.*(1+pi/2*u1.^2))+pi/2*u2;
        case 10 % sinus combined with quad
            y = 2*sin(2*pi*u1).*(u2-0.5).^2;
        case 11 % Branin function
            u1f = u1*15-5; u2f = u2*15; y = (u2f-(5.1/(4*pi^2))*u1f.^2 + 5*u1f/pi-6).^2 + 10*(1-1/(8*pi))*cos(u1f)+10;
        case 12 % Radcos function
            y = cos(9*sqrt(u1.^2+u2.^2) + 2) + 0.5*cos(11*u1+2) + 15*( (u1-0.4).^2 + (u2-0.4).^2 ).^2;
        case 13 % Gaussian with different standard deviation
            y = exp( -( 15*(u1-0.3).^2 + 38*(u2-0.5).^2 - 32*(u1-0.3).*(u2-0.5) ) );
        case 3 % parabola (symmetric)
            y = 1 - (u1-0.5).^2 - (u2-0.5).^2;
            name = 'parabola (symmetric)';
        otherwise
            error('invalid example')
    end
elseif n>=3
    switch bsp
        case 1 % hyperbel
            y = hyperbel(u,n);
            name = 'hyperbel';
        case 2 % parabola
            y = parabola(u,n);
            name = 'parabola';
        case 3 % perm function
            u = u*2*n-n;
            name = 'perm function';
            y = perm(u,n);
        case 4 % trid function
            u = u*2*n^2-n^2;
            name = 'trid function';
            y = trid(u,n);
        case 5 % Zakharov function
            u = u*15-5;
            y = zakh(u,n);
        case 6 % Radcos function
            y = Radcos(u,n);
        otherwise
            error('invalid example')
            
    end
else
    error('invalid dimension of data sample')
end

% Normalization to interval [0,1]
y = (y-min(y))./(max(y)-min(y));
end

%% function expressions

function y = hyperbel(u,n)
%
% Hyperbel function y = 0.1/(0.1+(1-u(1))/n+(1-u(2))/n+ ... +(1-u(n))/n)
% Matlab Code by B. Hartmann (Oct. 13, 2010)
%
den = 0.1;
for i = 1:n
    den = den+(1-u(:,i))/n;
end
y = 0.1./den;
end

function y = parabola(u,n)
%
% Parabola function y = 1-sum (u(i)-0.5)^2
% Matlab Code by B. Hartmann (Oct. 13, 2010)
%
y = 1;
for i = 1:n
    y = y-(u(:,i)-0.5).^2;
end
end

function y = perm(u,n)
%
% Perm function
% Matlab Code by A. Hedar (Nov. 23, 2005).
% The number of variables n should be adjusted below.
% The default value of n = 4.
%
b = 0.5;
s_out = 0;
for k = 1:n;
    s_in = 0;
    for j = 1:n
        s_in = s_in+(j^k+b)*((u(:,j)/j).^k-1);
    end
    s_out = s_out+s_in.^2;
end
y = s_out;
end

function y = trid(u,n)
%
% Trid function
% Matlab Code by A. Hedar (Nov. 23, 2005).
% The number of variables n should be adjusted below.
% The default value of n = 10.
%
s1 = 0;
s2 = 0;
for j = 1:n;
    s1 = s1+(u(:,j)-1).^2;
end
for j = 2:n;
    s2 = s2+u(:,j).*u(:,j-1);
end
y = s1-s2;
end

function y = Radcos(u,n)
F1 = ones(size(u,1),1);
F2 = zeros(size(u,1),1);
for d = 3:n
    F1 = F1 + u(:,d);
    b = 1+5*(d-3);
    a = 0.6/(1+b);
    F2 = F2 + (u(:,d) - a*(b*u(:,1) + 1 - u(:,2)) -0.2).^2;
end

y = cos(9*sqrt(u(:,1).^2+u(:,2).^2) + 2) + 0.5*cos(8*u(:,1).*( F1 ) +2) + 15*( (u(:,1)-0.4).^2 + (u(:,2)-0.4).^2 ).^2 + 15./(sqrt(n-2))*F2;

end