function disD = discretizeDistribution(disC, varargin)
% discretize a continuous distribution
%
% Syntax: disD = discretizeDistribution(disC)
%
% The standard method samples 100 points by a uniform distribution in the 
%   unit box with lower corner [0, .., 0] and upper corner [1, .., 1]. The 
%   input disC is generated by obtainDistributionContinuous.
%
% discretizeDistribution(..,  'argument', value) specifies additional
%   options of the procedure. The method itsef is selected by 'method' and
%   available values, including additional parameters are:
%     'uniform', 'box uniform' or 'uniform box' - samples uniformly inside a
%        rectangular box. The parameters 'regionLower' and 'regionUpper' 
%        specify the lower and upper coordinate of the box.
%   'normal' - uses importance sampling from a normal distribution. The
%        parameter 'mean' is the mean value and 'cov' represents the
%        covariance matrix. 'cov' might be a vector that is transformed to
%        a diagonal covariance matrix.
%
% The output structure contains all fields from the input distribution, 
%   while 'type' is set to 'disC'. The following fields are added:
%     parDiscrete - the method and the parameters used for generation
%     positionMatrix - the coordinates (column index) of each sample 
%       point (column vectors)
%     valueVector - the assigned value for the sample points from
%       disC.fhandle
%     probabilityVector - the probability density for selection of 
%       each sample point.
% 
% See also: obtainDistributionContinuous, obtainDistributionDiscrete

% 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/>.

%% Yet unsupported sampling:
%   'uniform sphere' or 'sphere uniform' - samples uniformly in a ball. The
%       parameters 'ballCenter' and 'ballRadius' define the ball.
%   'importance', 'box importance' or 'importance box' - uses importance
%       sampling in a rectangular box, where the importance is based on the
%       density function. Parameters as for 'uniform'. Note that a lot more
%       points are sampled to obtain the total nSample points.
%   'importance sphere' or 'sphere importance' - uses importance sampling
%       for a ball
%   'quasi box', 'box quasi' - uses the Sobolov series (a low discrepancy
%       series) to sample from a rectangular box.
%
% Importance sampling does not improve performance for Monte Carlo 
%   integrations as rejection sampling from an original uniform
%   distribution is used.

%% input handling
[o.nSample, o.method, o.regionLower, o.regionUpper, o.regionCenter, o.radius, ...
    o.mean, o.cov, optionStructIn] = mapOption(varargin, ...
    'nSample', 100, ...
    'method', 'box uniform', ...
    'regionLower', zeros([disC.nDim, 1]), ...
    'regionUpper', ones([disC.nDim, 1]), ...
    'regionCenter', zeros([disC.nDim, 1]), ...
    'regionRadius', 1, ...
    'mean', [], ...
    'cov', [], ...
    'options', struct());
optionStruct = mergeStruct(optionStructIn, o);

optionStruct.regionCenter = optionStruct.regionCenter(:);

%% initialization
disD = disC;
disD = rmfield(disD, 'fhandle');
disD.type = 'discrete distribution';
disD.nSample = optionStruct.nSample;

if strcmpi(optionStruct.method, 'uniform') || ...
   strcmpi(optionStruct.method, 'box uniform') || ...
   strcmpi(optionStruct.method, 'uniform box')

    %% case: uniform box
    % parameters
    regionLower = optionStruct.regionLower(:);
    regionUpper = optionStruct.regionUpper(:);
    disD.parDiscrete = struct('regionLower', regionLower, 'regionUpper', ...
        regionUpper);

    % sample point positions
    disD.positionMatrix = bsxfun(@plus, regionLower, ...
        diag(regionUpper-regionLower) * ...
        rand([disC.nDim, optionStruct.nSample]));
    
    % values at positions
    disD.valueVector = disC.fhandle(disD.positionMatrix);
    
    % probability at position
    disD.probabilityVector = disD.valueVector * 0 + ...
        1/prod(optionStruct.regionUpper - optionStruct.regionLower);
elseif strcmpi(optionStruct.method, 'normal')
    
    %% case: normal distribution
    % parameters
    [mean, cov] = ensureMeanAndCov(disC.nDim, optionStruct.mean, ...
        optionStruct.cov);
    
    disD.parDiscrete = struct('mean', mean, 'cov', cov);
    
    % assign position & probability
    [disD.positionMatrix disD.probabilityVector] = sampleNormal(disC.nDim, ...
        optionStruct.nSample, mean, cov);
    
    % values at positions
    disD.valueVector = disC.fhandle(disD.positionMatrix);    
    
else
    %% case: unknown >> error
    error('elk:distribution:wrongInput', ...
        ['The given method name (' method ') is not supported.']);
end