function out = integrateCuba(fhandle, xmin, xmax, varargin)
% integrateCuba - multidimensional integration using the Cuba library
%
% Syntax: res = integrateCuba(fhandle, xmin, xmax, ...)
%
% This function wraps the Cuba library, licenced under GNU lesser general
%   public licence, version 3.0. It is programmed and maintained by Thomas
%   Hahn, Max Planck Institute, Munich and the code is available online:
%   http://www.feynarts.de/cuba/.
%
% integrateCuba integrates a given function handle of the form @(x) in the
%   rectangular boundaries given by the corners xmin and xmax. Further
%   options can be added by name, value pairs. In particular, the algorithm
%   is selected by: integrateCuba(fhandle, xmin, xmax, 'algorithm',
%   'vegas').
% Relative (epsrel) and absolute tolerances (epsabs) can be provide, while
%   only one will be fulfilled. In the standard parameters, epsabs is set
%   to 1e-20 so that epsrel governs the error of the result in most of the
%   cases.
%
% All options are described in the documentation of the Cuba package.
%   Therefore, the names of the parameters are unaltered. Note, that the
%   following general parameters are computed automatically from the input:
%     * ndim, ncomp, userdata
%   some more parameters are general to all algorithms:
%     * epsrel (1e-2), epsabs (1e-20), flags (16), seed (0), mineval (0),
%       maxeval (5e4)
% All output variables are put in the structure res. These are:
%     * value, error, prob, fail, nregions (not in vegas)
%   Additional output is added with the field relError = value ./error.
% A warning is thrown when [fail] is not equal to zero, or [prob] is not 
%   below 0.95 for all components.
%     
% Algorithm: VEGAS. Vegas supports the following parameters:
%     * nstart (1e2), nincrease (5e2), nbatch (1e3)
%   while the following parameters are not supported by this MATLAB
%   wrapper:
%     * gridno (0), statefile (NULL)
%
% Algorithm: DIVONNE. Divonne supports the following parameters:
%     * key1 (7), key2 (-1), key3 (0), ngiven (automatic), ldxgiven
%       (automatic), xgiven ([])
%   the following parameters are not supported:
%     * maxpass (5), border (0.0), maxchisq (10), mindeviation (0.25),
%       nextra (0), peakfinder_t (NULL)
%
% Selegting CLEANUP for the algorithm deletes all compiled files, while 
%   selecting COMPILE compiles all MEX-files, including a testrun.

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 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/>.

%% ADDING A NEW CUBA LIBRARY
% To enable verbosity from the Cuba algorithms, you need to find the
%   line(s): "#define Print(s) .." and replace it/them by:
%     #define Print(s) mexPrintf(s);
% Currently, the version 3.1 does not run while it is suspected that it
%   might be according to the added parallel computing support that does
%   not work well with MATLAB.

%% Input handling

% ensure fhandle, xmin, xmax - this is done to give the right input
%   parameters for the testAlgorithm cases. And also to provide provide
%   generic cases for the compile option.
if ~exist('fhandle', 'var') || isempty(fhandle)
    fhandle = @(x)(sum(x, 1));
end
if ~exist('xmin', 'var') || isempty(xmin)
    xmin = [0 0 0]';
end
if ~exist('xmax', 'var') || isempty(xmax)
    xmax = xmin*0 + 1;
end

% parameters and options
optStruct = mapOptionStruct(varargin, ...
    'algorithm', 'vegas', ...
    'epsrel', 1e-2, 'epsabs', 1e-20, ...
    'seed', 0, 'flags', 16, ...
    'mineval', 0, 'maxeval', 5e4, ...
    'nstart', 1e3', 'nincrease', 5e2, 'nbatch', 1e3, ...
    'key', 7, 'key1', 7, 'key2', 7, 'key3', 0, 'xgiven', [], ...
    'flatness', 25, 'nnew', 1000, ...
    'cubaDirectory', 'Cuba-2.1', ...
    'gccCall', 'gcc-4.4 -fPIC');

% automatically generated parameters
optStruct.ndim = length(xmin);
optStruct.ncomp = length(fhandle(xmin));
if ~isempty(optStruct.xgiven) && ...
        size(optStruct.xgiven, 1) ~= optStruct.ndim
    error('elk:externCuba:wrongInput', ['The provided peakMatrix ' ...
        'must have the size [' num2str(optStruct.xdim) ']x[ngiven].']);
end
optStruct.ngiven = size(optStruct.xgiven, 2);
% standard output
out = struct('value', [], 'error', [], 'prob', [], 'fail', 0, ...
    'relError', [], 'time', []);

%% Input handling
switch lower(optStruct.algorithm)
    case 'cleanup'
        warnStruct = warning('off', 'MATLAB:DELETE:FileNotFound');
        try
            myDirectory = fileparts(which(mfilename));
            cubaDirectory = [myDirectory filesep optStruct.cubaDirectory];
            eval(['delete ' cubaDirectory filesep 'libcuba.a']);
            eval(['delete ' cubaDirectory filesep 'makefile']);
            eval(['delete ' cubaDirectory filesep 'config.h']);
            eval(['delete ' cubaDirectory filesep 'config.log']);
            eval(['delete ' cubaDirectory filesep 'config.status']);
            eval(['delete ' myDirectory filesep 'libcuba.a']);
            eval(['delete ' myDirectory filesep 'cuba.h']);
            eval(['delete ' myDirectory filesep '*.' mexext]);
            warning(warnStruct);
        catch thisError
            % ensure warning is on again
            warning(warnStruct);
            rethrow(thisError);
        end
    case 'compile'
        out = integrateCuba([], [], [], 'algorithm', 'vegas')
        out = integrateCuba([], [], [], 'algorithm', 'suave')
        out = integrateCuba([], [], [], 'algorithm', 'divonne')
        out = integrateCuba([], [], [], 'algorithm', 'cuhre')
        
    case 'testvegas'
        generateMexFile('testCubaVegas', optStruct);
        testCubaVegas;
        
    case 'vegas'
        generateMexFile('integrateCubaVegas', optStruct);
%         setenv('CUBACORES', '2');
        tic
        out = integrateCubaVegas(fhandle, xmin, xmax, optStruct);
        out.time = toc;
        
    case 'suave'
        generateMexFile('integrateCubaSuave', optStruct);
%         setenv('CUBACORES', '2');
        tic
        out = integrateCubaSuave(fhandle, xmin, xmax, optStruct);
        out.time = toc;
        
    case 'divonne'
        generateMexFile('integrateCubaDivonne', optStruct);
%         setenv('CUBACORES', '2');
        tic
        out = integrateCubaDivonne(fhandle, xmin, xmax, optStruct);
        out.time = toc;
        
    case 'cuhre'
        generateMexFile('integrateCubaCuhre', optStruct);
%         setenv('CUBACORES', '2');
        tic
        out = integrateCubaCuhre(fhandle, xmin, xmax, optStruct);
        out.time = toc;
    otherwise
        error(['The specified algorithm ''' algo '''is not supported.']);
end

%% general postprocessing
out.relError = out.error ./ out.value;
if out.fail ~= 0
    warning('elk:externCuba:integrationFailed', ...
        'Integration failed, likely because tolerance cannot be reached');
end
if any(out.prob > 0.95)
    warning('elk:externCuba:integrationFailed', ['You might not want to ' ...
        'trust the given error bounds. See the Cuba documentation, ' ...
        'concerning the prob output']);
end
        
end
        
function generateMexFile(file, optStruct)
        % identify my location
        originalDirectory = cd;
        myDirectory = fileparts(which(mfilename));
        
        % nothing to do, if mexfile already exists
        if exist([myDirectory filesep file '.' mexext], 'file')
            return
        end

        % build cuba library
        if isunix
            try
                %% build cuba library object
                disp('--- BUILD CUBA LIRARY ---');
                cd([myDirectory filesep optStruct.cubaDirectory])
                % makefile should be adopted to use gcc4.1 (gcc4.2 throws
                % errors for divonne, gcc4.3 is not yet supported by MATLAB
                % and might have issues)
                
                % configure
                eval(['!./configure CC=''' optStruct.gccCall '''']);
                !make lib
                !cp cuba.h ../
                !cp libcuba.a ../
                
                %% build mexfile
                disp('--- BUILD MEXFILE ---');
                disp(['Name: ' file '.' mexext]);
                cd(myDirectory)
                eval(['mex CC="' optStruct.gccCall '" ' file '.c libcuba.a']);
                
                % restore original path
                cd(originalDirectory);
                
                disp('--- BUILD DONE ---');
            catch thisError
                % ensure we are back in the original directory
                cd(originalDirectory);
                rethrow(thisError);
            end
        %! DEV: another operating system might be added, via:
%         elseif ispc
%             % .. windows specific code
%         elseif ismac
%             % .. mac specific code
        %/
        
        else
            error(['I cannot automatically build the cuba library from' ...
                'source for a non-unix operating system. Please ' ...
                'adopt the implementation to your system.']);
        end
end