function [X, Y, modelStruct, noise, variance] = exampleFunctions(name, inputDimension, dataDistribution, numberOfSamples, noiseLevel, modelStruct, dataGap)
% exampleFunctions is a simple script to create training and test data samples
%
% [X, Y, modelStruct, noise,variance] = exampleFunctions(name, inputDimension,
% dataDistribution, numberOfSamples, noiseLevel, modelStruct, dataGap)
%
% INPUTS:
%
%   name
%   (string) name of the examples function or data set. Options are:
%       'polynom'
%       'doppler'
%       'stenman'
%       'mars1'
%       'sinus1'
%       'sinus2'
%       'step'
%       'quadratic'
%       'linear'
%       'exponential'
%       'gaussian'
%
%   inputDimension
%   (scalar) specify the number of dimensions the input shall posses
%
%   dataDistribution
%   (string) specify how the data samples should be distributed. Options
%   are:
%       'grid' for data on a grid 0...1
%       'random' for random data 0...1
%       'spacefill' for pseudo monte-carlo distribution 0...1
%       'logaritmic' for a logarithmic grid
%
%   numberOfSamples
%   (scalar) the number of samples
%
%   noiseLevel
%   specify a noise level in relation to the variance of the data. 0 for no
%   noise, 0.05 is for low noise, 1 would totaly obscure the underlying
%   function
%
%   modelStruct
%   you may hand an old example over to conveniently create test data for
%   the training data. Example:
%
%       [trainingInput, trainingOutput, modelStruct] = exampleFunctions(...);
%
%       [testInput, testOutput] = exampleFunctions
%           (...,<differentDataDistribution>,...,modelStruct);
%
%   dataGap
%   if you want a gap or "hole" in your data, you may give the start and
%   end point here. For example [0.4 0.6] would induce a cube of emptiness
%   in a 3D input space
%
%
% OUTPUTS:
%   X
%   (n x d) the input data for training or testing
%
%   Y
%   (n x 1) the output data for training or testing
%
%   modelStruct
%   a struct with a specifications of the example
%
%   noise
%   (n x 1) if you specified a noiselevel for the input, the absulute noise is
%   given here
%
%   variance
%   (scalar) the true variance of the noise
%
%
% EXAMPLES
%
%   example of a exponential function ('exponential') in 2D (2) on a grid
%   ('linear') with 200 samples (200) without noise (0):
%
%   [trainingInput, trainingOutput] = exampleFunctions('exponential', 2, 'linear', 200, 0)
%
%
%   if you want to same example as above, but you also want 33 random
%   distributed test data samples with low noise (0.05), you could to it
%   like this:
%
%   % create the training data:
%   [trainingInput, trainingOutput, exampleStruct] = exampleFunctions('exponential', 2, 'linear', 200, 0)
%   % create the test data
%   [testInput, testOutput] = exampleFunctions([], [], 'random', 33, 0.05, exampleStruct)
%
%--------------------------------------------------------------------------

% Tobias Ebert, 10-January-2011
% Institut für Mess- und Regelungstechnik, Universität Siegen, Germany
% Copyright (c) 2011 by Tobias Ebert

% Update: Added new example functions. 31.01.2012, Benjamin Hartmann


%% create model

% test if an old modelStruct is given
if exist('modelStruct','var') && ~isempty(modelStruct)
    % use specifications from old example
    
    if ~isstruct(modelStruct)
        error('createTestfata:notAnOldModel','The modelStruct you have given is not useable (not a struct)')
    end
    
    % an old model exists, use it
    if ~isempty(name) || ~isempty(inputDimension)
        error('createTestdata:oldModel','an old model was given, leave <name> and <inputDimension> empty!')
    end
    
    % load old model
    inputDimension = modelStruct.inputDimension;
    datafun = modelStruct.datafun;
    name = modelStruct.name;
    
    if ~exist('dataDistribution','var') || isempty(dataDistribution)
        dataDistribution = modelStruct.dataDistribution;
    end
    
    if ~exist('numberOfSamples','var') || isempty(numberOfSamples)
        numberOfSamples = modelStruct.numberOfSamples;
    end
    
    if ~exist('noiseLevel','var') || isempty(noiseLevel)
        noiseLevel = modelStruct.noiseLevel;
    end
    
else
    % build a totaly new example
    
    if ~exist('dataDistribution','var') || isempty(dataDistribution)
        dataDistribution = 'linear';
    end
    if ~exist('numberOfSamples','var') || isempty(numberOfSamples)
        numberOfSamples = 100;
    end
    if ~exist('noiseLevel','var') || isempty(noiseLevel)
        noiseLevel = 0;
    end
    
    % create new testmodel
    if inputDimension<=2 && ~ischar(name)
        datafun = @lowDimExamples;% lowDimExamples(x,name);
    else
        switch name
            case 'polynom'
                % random number of points for polynom
                if inputDimension==1
                    nop = ceil(rand(1)*8) + 2;
                    x = rand(nop,1);
                    y = rand(nop,1);
                    P = polyfit(x,y,nop-1);
                    datafun = @(x) polyval(P,x);
                else % dimension > 1
                    % number of points for polynom
                    nop = ceil(rand(1)*8^inputDimension) + 2*inputDimension;
                    x = rand(nop,inputDimension);
                    y = rand(nop,1);
                    P = polyfitn(x,y,ceil(sqrt(nop))-1);
                    datafun = @(x) polyvaln(P,x);
                end
            case 'doppler'
                %datafun = @(x) 20.*sqrt(x.*(1-x)) .* sin(2.*pi.*(1.05 ./ (x+0.05))); % Doppler
                datafun = @(x) 20.*sqrt((prod(x,2).^(1./inputDimension)).*(1-(prod(x,2).^(1./inputDimension)))) .* sin(2.*pi.*(1.05 ./ ((prod(x,2).^(1./inputDimension))+0.05))); % Doppler
            case 'stenman'
                % datafun = @(x) 3*exp(-x.^2/0.3^2) + 2*exp(-(x-1).^2/0.7^2); % Stenman
                datafun = @(x) 3*exp(-(prod(x,2).^(1./inputDimension)).^2/0.3^2) + 2*exp(-((prod(x,2).^(1./inputDimension))-1).^2/0.7^2); % Stenman
            case 'mars1' % 2D
                if inputDimension ~= 2
                    error('data function <mars1> is for 2D only')
                end
                datafun = @(u) 2*exp(8*((u(:,1)-0.5).^2 + (u(:,2)-0.5).^2)) ./ ...
                    (exp(8*((u(:,1)-0.2).^2 + (u(:,2)-0.7).^2)) + exp(8*((u(:,1)-0.7).^2 + (u(:,2)-0.2).^2)));
            case 'sinus1'
                datafun = @(x) (sin(2*pi*(prod(x,2).^(1./inputDimension)))+2*(prod(x,2).^(1./inputDimension)))/2; % extent do 2d!
            case 'sinus2'
                datafun = @(x) sin(2*pi*10*(prod(x,2).^(1./inputDimension)));
            case 'step'
                datafun = @(x) round((prod(x,2).^(1./inputDimension)));
            case 'quadroot'
                %datafun = @(x) (prod(x-0.5,2).^(1./dimension)).^2;
                datafun = @(x) real((prod(x-0.5,2).^(1./inputDimension)).^2); % 14.07.2011
            case 'linroot'
                datafun = @(x) real((prod(x,2).^(1./inputDimension))+1);
            case 'hyperbola'
                % datafun = @(x) 1./(0.1+size(x,2)-(1-sum(x,2)));
                datafun = @(x) 0.1./(0.1+1-(1/size(x,2)*sum(x,2)));
            case 'gaussian'
                datafun = @gaussian;
            case 'gaussian2'
                datafun = @gaussian2;
            case 'parabola'
                datafun = @parabola;
            case 'perm'
                datafun = @perm;
            case 'trid'
                datafun = @trid;
            case 'hypGauss'
                datafun = @hypGauss;
            case 'hypSin'
                datafun = @hypSin;
            case 'zakh'
                datafun = @zakh
            otherwise
                error('name not known')
        end
    end
end


%% create input data

switch dataDistribution
    
    case 'grid'
        
        % 1 Dimension
        X = linspace(0,1,ceil(numberOfSamples^(1/inputDimension)))'; % model input
        
        %X = [linspace(0.1,0.6,ceil(numberOfSamples^(1/inputDimension)))'; linspace(0.85,0.9,2)']; % model input
        
        % n Dimensions
        if inputDimension>1 % the following part is from ndgrid
            siz(1:inputDimension) = length(X);
            newX = cell(1,inputDimension);
            for i=1:inputDimension,
                s = siz; s(i) = []; % Remove i-th dimension
                preX = reshape(X(:,ones(1,prod(s))),[length(X) s]); % Expand x
                newX{i} = permute(preX,[2:i 1 i+1:inputDimension]); % Permute to i'th dimension
            end
            X = zeros(size(newX{1}(:),1),inputDimension);
            % do a nice reshaping
            for k = 1:inputDimension
                X(:,k) = newX{k}(:);
            end
        end
        
        if size(X,1) ~= numberOfSamples
            warning('%1.0f instead of %2.0f number of samples were created',size(X,1),numberOfSamples)
        end
        
    case 'random'
        
        X = rand( numberOfSamples ,inputDimension);
        %X = sortrows(X);
        
    case 'spacefill' % good for high dimensions
        
        X = randomSpaceFill(zeros(1,inputDimension),ones(1,inputDimension),numberOfSamples,numberOfSamples*5);
        
    case 'logarithmic'
        
        % one input dimension
        X = logspace(-2,0,ceil(numberOfSamples^(1/inputDimension)))'; % model input
        
        % for n input dimensions
        if inputDimension>1 % the following part is from ndgrid
            siz(1:inputDimension) = length(X);
            newX = cell(1,inputDimension);
            for i=1:inputDimension,
                s = siz; s(i) = []; % Remove i-th dimension
                preX = reshape(X(:,ones(1,prod(s))),[length(X) s]); % Expand x
                newX{i} = permute(preX,[2:i 1 i+1:inputDimension]); % Permute to i'th dimension
            end
            X = zeros(size(newX{1}(:),1),inputDimension);
            % do a nice reshaping
            for k = 1:inputDimension
                X(:,k) = newX{k}(:);
            end
        end
        
    otherwise
        
        error('<distribution> is unknown')
        
end



%% delete points in gap

if exist('dataGap','var') && ~isempty(dataGap)
    X(logical(floor(sum((X>=dataGap(1) & X<=dataGap(2)),2)/inputDimension)),:) = [];
end


%% calculate data output
if ~ischar(name)
    Y = datafun(X,name);
else
    Y = datafun(X); % model output
end

%% add noise if noiselevel is given

if ~isempty(noiseLevel)
    
    if iscell(noiseLevel)
        
        switch noiseLevel{1}
            case 'rising' % left zero right maximum
                %noise = X(:,1) * noiselevel{2}  * (max(Y(:))-min(Y(:))) .* randn(size(X,1),1);
                variance = (X(:,1) * noiseLevel{2}  * (max(Y(:))-min(Y(:)))).^2;
            case 'quadratic' % left zero right maximum
                %noise = X(:,1) * noiselevel{2}  * (max(Y(:))-min(Y(:))) .* randn(size(X,1),1);
                variance = (X(:,1).^2 * noiseLevel{2}  * (max(Y(:))-min(Y(:)))).^2;
            case 'block'
                %noise = noiselevel{2}  * (max(Y(:))-min(Y(:)))  * randn(size(X,1),1);
                %noise(1:floor(length(noise)/2)) = 0;
                variance = (noiseLevel{2}  * (max(Y(:))-min(Y(:))) *ones(size(Y,1),1)).^2;
                variance(1:floor(length(variance)/2)) = 0;
            otherwise
                warning('not known')
        end
        
    elseif isnumeric(noiseLevel)
        
        % calculate the true variance of the noise
        variance = (noiseLevel * (max(Y(:))-min(Y(:))) * ones(size(Y,1),1) ).^2;
    end
    noise = sqrt(variance) .*randn(size(Y,1),1);
    Y = Y+noise; % add noise to model output
else
    noise = [];
    variance=0;
end

%% create modelStruct
modelStruct.datafun = datafun;
modelStruct.name = name;
modelStruct.inputDimension = inputDimension;
modelStruct.numberOfSamples = numberOfSamples;
modelStruct.dataDistribution = dataDistribution;
modelStruct.noiseLevel = noiseLevel;
if exist('P','var')
    modelStruct.P = P;
end


end


%% subfunction to create pseudo monte-carlo data distribution
function u = randomSpaceFill(Min,Max,N,Nrand)


% Generate homogeneous distributed data samples.
% Procedure description:
% For all random samples (candidate points) the distances to all existing points are calculated and the minimal distance dMin is determined.
% Next, the point is added to the data set that has maximum dMin.

% Candidate points
dim = size(Min,2); % dimension of the Problerm
% uRand = rand(Nrand,dim) * diag(Max-Min);
z_qmc = sobolset(dim);
uRand = net(z_qmc,Nrand) * diag(Max-Min);
for k = 1:dim
    uRand(:,k) = uRand(:,k)+Min(k); % add offset
end

% Initiatlize vectors
index = [];
newPoints = [];
existingPoints = uRand(1,:);
dMin = inf(1,size(uRand,1));
for k = 1:N
    % Evaluate the distances of candidate points to the lastly added point u(k,:)
    distanceOfk = sqrt(sum(bsxfun(@minus,uRand,existingPoints(k,:)).^2,2))';
    % Take over the smallest distances in the distance vector dMin
    dMin = min([distanceOfk;dMin],[],1);
    
    if k >= size(existingPoints,1)
        % Select the point with the largest dMin-value
        [~, idxMax] = max(dMin);
        % Update existing points with selected candidate point
        existingPoints = [existingPoints; uRand(idxMax,:)];
        % Update index vector
        index = [index; idxMax];
        newPoints = [newPoints; uRand(idxMax,:)];
        
    end
end
u = newPoints;

end





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

%% perm
function y = perm(u)
%
% 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.
%
n = size(u,2);
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

%% trid
function y = trid(u)
%
% 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.
%
n = size(u,2);
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

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

%% gauss
function y = gaussian2(u,c,sigma)
[N n] = size(u);
if nargin == 1
    c     = ones(1,n)/2; % Centers for Gaussian
    sigma = 0.25;        % Standard deviation (equally for all inputs)
end
x = zeros(N,1);
for i = 1:n
    x = x + ((u(:,i)-c(i))/sigma).^2;
end
y = exp(-0.5*x);
end

%% hypGauss
function y = hypGauss(u)
% hyperbel + gauss
n       = size(u,2);
c       = zeros(1,n); % Centers for Gaussian
sigma   = 0.25;        % Standard deviation (equally for all inputs)
hheight = 1;          % Maximum height of the Gaussian function
gheight = 0.75;        % Maximum height of the hyperbola function
y = hheight*hyperbel(u) + (gheight-hyperbel(c))*gaussian2(u,c,sigma);
end

%% hyperbel
function y = hyperbel(u)
%
% 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)
%
n = size(u,2);
den = 0.1;
for i = 1:n
    den = den+(1-u(:,i))/n;
end
y = 0.1./den;
end

%% hypSin
function y = hypSin(u)
p = size(u,2);
y = (1+sin(3*pi*sum(u,2)/p +1*pi/2))/1.5.*(1-sum(u,2)/p);
% y = sin(1./(0.15+sum(u,2)/p));
end

%% 1D and 2D example functions
function y = lowDimExamples(u,bsp)
n = size(u,2);
if n==1
    switch bsp
        case 1
            y = sqrt(2)*sin(u) + sqrt(2)*cos(u) + 1/(2*sqrt(2))*sin(2*u) + 1/(sqrt(2))*cos(2*u) - sqrt(2)*sin(3*u) + sqrt(2)*cos(3*u);
        case 2
            y = zeros(size(u,1),1);
            y(u<0.3)=8-30*u(u<0.3).^3;
            y(u>=0.3) = 5+5*u(u>=0.3).^2;
        case 3
            v1 = 1;
            v0 = 0.65;
            smoothness = 3;
            sigmoide = 1./(1+exp((-v0+v1*u)/v1/smoothness*20));
            y2 = 0.1./(0.1+u);
            y1 = 1+sin(4*pi*u-pi/8);
            y = y1.*sigmoide + y2.*(1-sigmoide);
        case 4
            v1 = 1;
            v0 = 0.65;
            smoothness = 3;
            sigmoide = 1./(1+exp((-v0+v1*u)/v1/smoothness*20));
            y2 = 0.1./(0.1+u);
            y1 = 1+sin(3*pi*u-pi/8);
            y = y1.*sigmoide + y2.*(1-sigmoide)+1;
    end
else
    u1=u(:,1); u2=u(:,2);
    switch bsp
        case 1
            y = 0.1./(0.1+0.5*(1-u1)+0.5*(1-u2));
        case 2
            y = 0.1./(0.1+0.3*(1-u1)+0.7*(1-u2));
        case 3
            y = u1.^2.*sqrt(u2);
        case 4
            y = 0.01./((1.1-u1).*(1.1-u2));
        case 5
            y = 2-2*exp(-0.5*((u1-0.25).^2+(u2-0.25).^2));
        case 6 %Mars
            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
            y = 2*sin(2*pi*u1).*(u2-0.5).^2;
        case 8
            y = exp(-25*((u1-0.5).^2+(u2-0.5).^2));
        case 9
            y = 0.25*cos(3*pi/2*u1.*(1+pi/2*u1.^2))+pi/2*u2;
        case 10
            y = 2*sin(2*pi*u1).*(u2-0.5).^2;
        case 11
            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; % Branin
        case 12
            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
            y = exp( -( 15*(u1-0.3).^2 + 38*(u2-0.5).^2 - 32*(u1-0.3).*(u2-0.5) ) );
            
        case 14
            idx = u1>1-u2;
            u11 = u1(idx);  u21 = u2(idx);
            u12 = u1(~idx); u22 = u2(~idx);
            y = zeros(size(u,1),1);
            y(idx)  = 5 + 3*u11 + 7*u21.^3;
            %y(~idx) = 1 + u12.^2 + 0.5*u22;
        case 15
            y  = 5 + 3*u1 + 7*u2.^3;
        case 16
            y = sin(4*pi*(u1 + u2.^2)/2);
        case 17
            y = 0.5./(0.5+0.5*(1-u1)+0.5*(1-u2.^2));
    end
end
end