function [U] = GridDef(lbmub,varargin)

% This functions defines a grid for sampling using three different
% approaches:
%  - Regular grids
%  - Latin Hypercube sampling with the function lhsdesign.m
%  - Orthogonal sampling with the function orthogonal_sample.m
%  - Sobol sampling with the function sobolset.m
%
% The input of this function is given by
%       lbmub:	3 x ntot matrix, with first row determining the minimum
%               value, the second row the reference value and the third row
%               the maximum value. This matrix includes all independent
%               variables, also the unvaried.
%
%       indin:  Indices ofthe lbmub matrix, which should be varied
%
%       opt:	Structure containing the options of the system, includes
%               .type:  String containing CG, LHS or OD
%               .Np:    CG:     Points per variable
%                       LHS:    Total number of points
%                       OD:     Number of points per subdivision
%                       RND:    Random numbers generated
%                       SOB:    Sobol sequence sampling
%               .Ns     OD:     Number of subdivision per variable

% Initialize the global variable SOB giving the sobol sequence defined in
% the first section
global SOB

% Safety function to use all available input points and a LHS with 1000 points in
% the case no inputs or grids are defined.
if nargin < 2
    indin = 1:size(lbmub,2);
    opt.type = 'LHS';
    opt.Np = 1e3;
    opt.val = 'y';
elseif nargin < 3
    indin = varargin{1};
    opt.type = 'LHS';
    opt.Np = 1e3;
    opt.val = 'y';
elseif nargin < 4
    indin = varargin{1};
    opt = varargin{2};
end

% Definition of the number of varied independent variables and the total
% number of variables
nvar = length(indin);
ntot = size(lbmub,2);
nsamp = opt.Np;

switch opt.type
    case 'CG'
        
        % Predefinition of the cell vector for the varied values
        vec = cell(1,ntot);

        for k = indin
            vec{k} = linspace(0,1,opt.Np);
        end

        grid_comb = cell(1,ntot);
        [grid_comb{indin}] = ndgrid(vec{indin});

        for k = indin
            grid_comb{k} = grid_comb{k}(:);
        end

        U.grid = cell2mat(grid_comb(indin));

        % Definition of the number of samples
        nsamp = opt.Np^nvar;

        U.var  = U.grid.*(ones(nsamp,1)*diff(lbmub([1 3],indin)))+ones(nsamp,1)*lbmub(1,indin);
        
        % Definition of the vectors for the feed composition without any change
        U.out = ones(nsamp,1)*lbmub(2,:);

        % Definition of the adjusted matrix
        U.out(:,indin) = U.var;
        
        % Rearragnement of the variables which need to be changed
        if isfield(opt,'dep')
            for k = 1:size(opt.dep,1)
                U.out(:,opt.dep(k,2)) = U.out(:,opt.dep(k,1))*opt.dep(k,3).*...
                                       (1+(U.grid(:,opt.dep(k,2))-.5)*2*opt.dep(k,4));%+opt.dep(k,4);
                
            end
        end
                
                     
    case 'RND'
        % Performing the random sampling
        U.init = rand(nsamp,nvar);
        
        if opt.val == 'n'
            % Predefinition of the cell vector for the varied values
            vec = cell(1,ntot);

            for k = indin
                vec{k} = linspace(0,1,2);
            end

            grid_comb = cell(1,ntot);
            [grid_comb{indin}] = ndgrid(vec{indin});

            for k = indin
                grid_comb{k} = grid_comb{k}(:);
            end

            U.CG = cell2mat(grid_comb(indin));
            
            % Definition of the number of samples
            nsamp = opt.Np+2^nvar;
        else
            U.CG = [];
            
            % Definition of the number of samples
            nsamp = opt.Np;
        end

        % Combination to a grid
        U.grid = [U.CG;...
                  U.init];
        
        U.var  = U.grid.*(ones(nsamp,1)*diff(lbmub([1 3],indin)))+ones(nsamp,1)*lbmub(1,indin);
        
        % Definition of the vectors for the feed composition without any change
        U.out = ones(nsamp,1)*lbmub(2,:);

        % Definition of the adjusted matrix
        U.out(:,indin) = U.var;
        
        % Rearragnement of the variables which need to be changed
        if isfield(opt,'dep')
            for k = 1:size(opt.dep,1)
                U.out(:,opt.dep(k,2)) = U.out(:,opt.dep(k,1))*opt.dep(k,3).*...
                                       (1+(U.grid(:,opt.dep(k,2))-.5)*2*opt.dep(k,4));%+opt.dep(k,4);
                
            end
        end
                
                
    case 'LHS'
        % Performing the LHS
        U.init = lhsdesign(nsamp,nvar,'criterion','maximin');
        
        if opt.val == 'n'
            % Predefinition of the cell vector for the varied values
            vec = cell(1,ntot);

            for k = indin
                vec{k} = linspace(0,1,2);
            end

            grid_comb = cell(1,ntot);
            [grid_comb{indin}] = ndgrid(vec{indin});

            for k = indin
                grid_comb{k} = grid_comb{k}(:);
            end

            U.CG = cell2mat(grid_comb(indin));
            
            % Definition of the number of samples
            nsamp = opt.Np+2^nvar;
        else
            U.CG = [];
            
            % Definition of the number of samples
            nsamp = opt.Np;
        end

        % Combination to a grid
        U.grid = [U.CG;...
                  U.init];
        
        U.var  = U.grid.*(ones(nsamp,1)*diff(lbmub([1 3],indin)))+ones(nsamp,1)*lbmub(1,indin);
        
        % Definition of the vectors for the feed composition without any change
        U.out = ones(nsamp,1)*lbmub(2,:);

        % Definition of the adjusted matrix
        U.out(:,indin) = U.var;
        
        % Rearragnement of the variables which need to be changed
        if isfield(opt,'dep') 
            for k = 1:size(opt.dep,1)
                if any(indin == opt.dep(k,2))
                    U.out(:,opt.dep(k,2)) = U.out(:,opt.dep(k,1))*opt.dep(k,3).*...
                                           (1+(U.grid(:,opt.dep(k,2))-.5)*2*opt.dep(k,4));%+opt.dep(k,4);
                end
            end
        end
                   
                
    case 'SOB'
        % Safety loops for the sobol sequence
        if ~isfield(opt.SOB,'leap')
            opt.SOB.leap = 0;
        end
        if ~isfield(opt.SOB,'skip')
            opt.SOB.skip = 1;
        end
        % Generating the sobol sequence
        SOB    = sobolset(nvar,'Skip',opt.SOB.skip,'Leap',opt.SOB.leap);
        U.init = net(SOB,nsamp);
        
        if opt.val == 'n'
            % Predefinition of the cell vector for the varied values
            vec = cell(1,ntot);

            for k = indin
                vec{k} = linspace(0,1,2);
            end

            grid_comb = cell(1,ntot);
            [grid_comb{indin}] = ndgrid(vec{indin});

            for k = indin
                grid_comb{k} = grid_comb{k}(:);
            end

            U.CG = cell2mat(grid_comb(indin));
            
            % Definition of the number of samples
            nsamp = opt.Np+2^nvar;
        else
            U.CG = [];
            
            % Definition of the number of samples
            nsamp = opt.Np;
        end

        % Combination to a grid
        U.grid = [U.CG;...
                  U.init];
        
        U.var  = U.grid.*(ones(nsamp,1)*diff(lbmub([1 3],indin)))+ones(nsamp,1)*lbmub(1,indin);
        
        % Definition of the vectors for the feed composition without any change
        U.out = ones(nsamp,1)*lbmub(2,:);

        % Definition of the adjusted matrix
        U.out(:,indin) = U.var;
        
        % Rearragnement of the variables which need to be changed
        if isfield(opt,'dep') 
            for k = 1:size(opt.dep,1)
                if any(indin == opt.dep(k,2))
                    U.out(:,opt.dep(k,2)) = U.out(:,opt.dep(k,1))*opt.dep(k,3).*...
                                           (1+(U.grid(:,opt.dep(k,2))-.5)*2*opt.dep(k,4));%+opt.dep(k,4);
                end
            end
        end     
end
