function disD = obtainDistributionDiscrete(name, varargin)
% obtain a discrete probability distribution
%
% Syntax: disD = obtainDistributionDiscrete(name, par)
%   OR    disD = obtainDistributionDiscrete(name, 'option', value, ..)
%
% This function returns a discrete realization of a probability
%   distribution. Common option/value pairs are: 
%      'nSample' - number of sample points {100}
%      'verbose' - verbosity level (0 is off, 1 is on) {0}
%
% The input "name" is one out of the list below. The input "par" is a 
%   structure containing the parameters of the distribution. The supported
%   distributions are:
%     uniform polytope: uniform distribution in a given polytope; 
%       parameters: poly (the polygon in V- or H-representation)
%     uniform ray: uniform distribution in a cone around a given ray,
%       parameters: ray (the center ray); radius (radius of the base at the 
%       end of the ray)
%     normal: multivariate Gaussian normal distribution with parameters: 
%       mean (vector with mean values); cov (matrix with covariances)
%     metropolis: the points are sampled by the Metropolis algorithm and
%       the first parameter set considers the Metropolis algorithm itself:
%       metroHanlde (distribution density @(x)), metroInit (an initial
%       point x for which metroHandle(x) should be large), metroCov (a
%       guess for the covariance matrix of sampled points), metroSpacing
%       (measure of the steps waited before a new sample is taken taken
%       relative to the time the autocorellation reaches a negative value
%       the first time), nSampleBatch (amount of points sampled in parallel
%       for speedup); the second set of parameters is used to estimate
%       the integral of the given distribution to calculate the probability
%       density: metroIntegral (the integral value that is directly used),
%       metroIntegralBatch (points to evaluate in each step of the Monte
%       Carlo Integration), metroRelTol (the coefficient of variation to
%       be obtained by the MC integration)
%
% The output disD is a structure with the fields:
%     name - the name of the distribution
%     nDim - the dimension of the distribution
%     par - the input parameters of the distribution
%     pivotPositionMatrix - the coordinates the sampled points (each is a 
%       column vector)
%     pivotProbabilityVector - the assigned value for the sample point
%
%
% See also: obtainDistributiondisContinuous, discretizeDistribution

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>

%% input handling
par = mapOptionStruct(varargin, 'nSample', 100, 'verbose', 0, ...
    'ray', 1, 'radius', 1, ...
    'mean', 1, 'cov', 1, ...
    'poly', obtainVrep('cube'), ...
    'metroInit', 1, 'metroCov', 1, 'metroHandle', @(x)(1./(x.^2+1)), ...
    'metroSpacing', 2, 'metroNonAdaptive', 10, 'metroNoHistory', 10, ...
    'metroHistory', 1000, 'metroLag', 100, 'metroShowAutoCor', 0, ...
    'metroRelTol', 0.01, 'metroIntegral', [], 'metroIntegralBatch', 1e5, ...
    'metroSampleBatch', 100 ...
    );

if strcmpi(name, 'ray') || strcmpi(name, 'ray uniform') || ...
        strcmpi(name, 'uniform ray')
    
    % basic information
    disD.name = 'uniform ray';
    disD.nDim = length(par.ray);
    
    
    [disD.positionMatrix disD.probabilityVector] = sampleRay(...
        par.nSample, par.ray(:),  par.radius);
    
    % parameter mapping
    disD.parDiscrete.ray = par.ray(:);
    disD.parDiscrete.radius = par.radius;

elseif strcmpi(name, 'poly') || strcmpi(name, 'uniform poly') || ...
        strcmpi(name, 'polytope') || strcmpi(name, 'uniform polytope') || ...
        strcmpi(name, 'poly uniform') || strcmpi(name, 'polytope uniform')
    
    % basic information
    disD.name = 'uniform polygon';
    infoPoly = identifyPolytope(par.poly);
    disD.nDim = infoPoly.dim;
    
    [disD.positionMatrix disD.probabilityVector] = samplePolygon(...
        par.nSample, par.poly);
    
    % parameter mapping
    disD.parDiscrete.poly = par.poly;
    
elseif strcmpi(name, 'normal') || strcmpi(name, 'normal distribution')
    % basic information
    disD.name = 'normal';
    disD.nDim = length(par.mean);
    [disD.parDiscrete.mean, disD.parDiscrete.cov] = ensureMeanAndCov(...
        disD.nDim, par.mean, par.cov);
       
    % assign position & probability
    [disD.positionMatrix disD.probabilityVector] = sampleNormal(disD.nDim, ...
        par.nSample, disD.parDiscrete.mean, disD.parDiscrete.cov);
   
elseif strcmpi(name, 'metro') || strcmpi(name, 'metropolis')
    % basic information
    disD.name = 'metropolis';
    
    [disD.positionMatrix disD.probabilityVector disD.resMetro] = ...
        sampleMetropolis(par, par.nSample);
    
    % parameter mapping
    disD.parDiscrete = par;
    
else
    
    error('elk:distribution:wrongInput', ...
        ['The given distribution name (' name ') is not supported.']);
end
