function X = testfunctionInput(numberOfInputs, numberOfSamples, dataDistribution)

% Generiert Matrix mit Eingaengen input = [u1 u2 u3 ... u_numberOfInputs] jeweils der Laenge numberOfSamples.
% Auswahl Anzahl der Eingangsgroessen: numberOfInputs = 1,2,...,10
% Auswahl Anzahl der Datenpunkte: numberOfSamples
%
% Benjamin Hartmann, 08. Oktober 2008

% if nargin<3
%    uMin = 0;
%    uMax = 1;
% end

% u = linspace(uMin,uMax,numberOfSamples)';
% switch numberOfInputs
%    case 1
%       input = [u];
%    case 2
%       [u1g u2g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:);
%       input = [u1 u2];
%    case 3
%       [u1g u2g u3g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:); u3=u3g(:);
%       input = [u1 u2 u3];
%    case 4
%       [u1g u2g u3g u4g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:); u3=u3g(:); u4=u4g(:);
%       input = [u1 u2 u3 u4];
%    case 5
%       [u1g u2g u3g u4g u5g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:); u3=u3g(:); u4=u4g(:);  u5=u5g(:);
%       input = [u1 u2 u3 u4 u5];
%    case 6
%       [u1g u2g u3g u4g u5g u6g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:); u3=u3g(:); u4=u4g(:);  u5=u5g(:); u6=u6g(:);
%       input = [u1 u2 u3 u4 u5 u6];
%    case 7
%       [u1g u2g u3g u4g u5g u6g u7g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:); u3=u3g(:); u4=u4g(:);  u5=u5g(:); u6=u6g(:); u7=u7g;
%       input = [u1 u2 u3 u4 u5 u6 u7];
%    case 8
%       [u1g u2g u3g u4g u5g u6g u7g u8g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:); u3=u3g(:); u4=u4g(:);  u5=u5g(:); u6=u6g(:); u7=u7g(:); u8=u8g(:);
%       input = [u1 u2 u3 u4 u5 u6 u7 u8];
%    case 9
%       [u1g u2g u3g u4g u5g u6g u7g u8g u9g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:); u3=u3g(:); u4=u4g(:);  u5=u5g(:); u6=u6g(:); u7=u7g(:); u8=u8g(:); u9=u9g(:);
%       input = [u1 u2 u3 u4 u5 u6 u7 u8 u9];
%    case 10
%       [u1g u2g u3g u4g u5g u6g u7g u8g u9g u10g] = ndgrid(u);
%       u1=u1g(:); u2=u2g(:); u3=u3g(:); u4=u4g(:);  u5=u5g(:); u6=u6g(:); u7=u7g(:); u8=u8g(:); u9=u9g(:); u10=u10g(:);
%       input = [u1 u2 u3 u4 u5 u6 u7 u8 u9 u10];
% end

inputDimension = numberOfInputs;
if nargin < 3
    dataDistribution = 'linear';
end
switch dataDistribution
    
    case 'linear'
        
        % 1 Dimension
        X = linspace(0,1,ceil(numberOfSamples^(1/inputDimension)))'; % 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 'logaritmic'
        
        % 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
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

