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.0001 function out = integrateCuba(fhandle, xmin, xmax, varargin) 0002 % integrateCuba - multidimensional integration using the Cuba library 0003 % 0004 % Syntax: res = integrateCuba(fhandle, xmin, xmax, ...) 0005 % 0006 % This function wraps the Cuba library, licenced under GNU lesser general 0007 % public licence, version 3.0. It is programmed and maintained by Thomas 0008 % Hahn, Max Planck Institute, Munich and the code is available online: 0009 % http://www.feynarts.de/cuba/. 0010 % 0011 % integrateCuba integrates a given function handle of the form @(x) in the 0012 % rectangular boundaries given by the corners xmin and xmax. Further 0013 % options can be added by name, value pairs. In particular, the algorithm 0014 % is selected by: integrateCuba(fhandle, xmin, xmax, 'algorithm', 0015 % 'vegas'). 0016 % Relative (epsrel) and absolute tolerances (epsabs) can be provide, while 0017 % only one will be fulfilled. In the standard parameters, epsabs is set 0018 % to 1e-20 so that epsrel governs the error of the result in most of the 0019 % cases. 0020 % 0021 % All options are described in the documentation of the Cuba package. 0022 % Therefore, the names of the parameters are unaltered. Note, that the 0023 % following general parameters are computed automatically from the input: 0024 % * ndim, ncomp, userdata 0025 % some more parameters are general to all algorithms: 0026 % * epsrel (1e-2), epsabs (1e-20), flags (16), seed (0), mineval (0), 0027 % maxeval (5e4) 0028 % All output variables are put in the structure res. These are: 0029 % * value, error, prob, fail, nregions (not in vegas) 0030 % Additional output is added with the field relError = value ./error. 0031 % A warning is thrown when [fail] is not equal to zero, or [prob] is not 0032 % below 0.95 for all components. 0033 % 0034 % Algorithm: VEGAS. Vegas supports the following parameters: 0035 % * nstart (1e2), nincrease (5e2), nbatch (1e3) 0036 % while the following parameters are not supported by this MATLAB 0037 % wrapper: 0038 % * gridno (0), statefile (NULL) 0039 % 0040 % Algorithm: DIVONNE. Divonne supports the following parameters: 0041 % * key1 (7), key2 (-1), key3 (0), ngiven (automatic), ldxgiven 0042 % (automatic), xgiven ([]) 0043 % the following parameters are not supported: 0044 % * maxpass (5), border (0.0), maxchisq (10), mindeviation (0.25), 0045 % nextra (0), peakfinder_t (NULL) 0046 % 0047 % Selegting CLEANUP for the algorithm deletes all compiled files, while 0048 % selecting COMPILE compiles all MEX-files, including a testrun. 0049 0050 % The elk-library: convex geometry applied to crystallization modeling. 0051 % Copyright (C) 2012 Alexander Reinhold 0052 % 0053 % This program is free software: you can redistribute it and/or modify it 0054 % under the terms of the GNU General Public License as published by the 0055 % Free Software Foundation, either version 3 of the License, or (at your 0056 % option) any later version. 0057 % 0058 % This program is distributed in the hope that it will be useful, but 0059 % WITHOUT ANY WARRANTY; without even the implied warranty of 0060 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0061 % General Public License for more details. 0062 % 0063 % You should have received a copy of the GNU General Public License along 0064 % with this program. If not, see <http://www.gnu.org/licenses/>. 0065 0066 %% ADDING A NEW CUBA LIBRARY 0067 % To enable verbosity from the Cuba algorithms, you need to find the 0068 % line(s): "#define Print(s) .." and replace it/them by: 0069 % #define Print(s) mexPrintf(s); 0070 % Currently, the version 3.1 does not run while it is suspected that it 0071 % might be according to the added parallel computing support that does 0072 % not work well with MATLAB. 0073 0074 %% Input handling 0075 0076 % ensure fhandle, xmin, xmax - this is done to give the right input 0077 % parameters for the testAlgorithm cases. And also to provide provide 0078 % generic cases for the compile option. 0079 if ~exist('fhandle', 'var') || isempty(fhandle) 0080 fhandle = @(x)(sum(x, 1)); 0081 end 0082 if ~exist('xmin', 'var') || isempty(xmin) 0083 xmin = [0 0 0]'; 0084 end 0085 if ~exist('xmax', 'var') || isempty(xmax) 0086 xmax = xmin*0 + 1; 0087 end 0088 0089 % parameters and options 0090 optStruct = mapOptionStruct(varargin, ... 0091 'algorithm', 'vegas', ... 0092 'epsrel', 1e-2, 'epsabs', 1e-20, ... 0093 'seed', 0, 'flags', 16, ... 0094 'mineval', 0, 'maxeval', 5e4, ... 0095 'nstart', 1e3', 'nincrease', 5e2, 'nbatch', 1e3, ... 0096 'key', 7, 'key1', 7, 'key2', 7, 'key3', 0, 'xgiven', [], ... 0097 'flatness', 25, 'nnew', 1000, ... 0098 'cubaDirectory', 'Cuba-2.1', ... 0099 'gccCall', 'gcc-4.4 -fPIC'); 0100 0101 % automatically generated parameters 0102 optStruct.ndim = length(xmin); 0103 optStruct.ncomp = length(fhandle(xmin)); 0104 if ~isempty(optStruct.xgiven) && ... 0105 size(optStruct.xgiven, 1) ~= optStruct.ndim 0106 error('elk:externCuba:wrongInput', ['The provided peakMatrix ' ... 0107 'must have the size [' num2str(optStruct.xdim) ']x[ngiven].']); 0108 end 0109 optStruct.ngiven = size(optStruct.xgiven, 2); 0110 % standard output 0111 out = struct('value', [], 'error', [], 'prob', [], 'fail', 0, ... 0112 'relError', [], 'time', []); 0113 0114 %% Input handling 0115 switch lower(optStruct.algorithm) 0116 case 'cleanup' 0117 warnStruct = warning('off', 'MATLAB:DELETE:FileNotFound'); 0118 try 0119 myDirectory = fileparts(which(mfilename)); 0120 cubaDirectory = [myDirectory filesep optStruct.cubaDirectory]; 0121 eval(['delete ' cubaDirectory filesep 'libcuba.a']); 0122 eval(['delete ' cubaDirectory filesep 'makefile']); 0123 eval(['delete ' cubaDirectory filesep 'config.h']); 0124 eval(['delete ' cubaDirectory filesep 'config.log']); 0125 eval(['delete ' cubaDirectory filesep 'config.status']); 0126 eval(['delete ' myDirectory filesep 'libcuba.a']); 0127 eval(['delete ' myDirectory filesep 'cuba.h']); 0128 eval(['delete ' myDirectory filesep '*.' mexext]); 0129 warning(warnStruct); 0130 catch thisError 0131 % ensure warning is on again 0132 warning(warnStruct); 0133 rethrow(thisError); 0134 end 0135 case 'compile' 0136 out = integrateCuba([], [], [], 'algorithm', 'vegas') 0137 out = integrateCuba([], [], [], 'algorithm', 'suave') 0138 out = integrateCuba([], [], [], 'algorithm', 'divonne') 0139 out = integrateCuba([], [], [], 'algorithm', 'cuhre') 0140 0141 case 'testvegas' 0142 generateMexFile('testCubaVegas', optStruct); 0143 testCubaVegas; 0144 0145 case 'vegas' 0146 generateMexFile('integrateCubaVegas', optStruct); 0147 % setenv('CUBACORES', '2'); 0148 tic 0149 out = integrateCubaVegas(fhandle, xmin, xmax, optStruct); 0150 out.time = toc; 0151 0152 case 'suave' 0153 generateMexFile('integrateCubaSuave', optStruct); 0154 % setenv('CUBACORES', '2'); 0155 tic 0156 out = integrateCubaSuave(fhandle, xmin, xmax, optStruct); 0157 out.time = toc; 0158 0159 case 'divonne' 0160 generateMexFile('integrateCubaDivonne', optStruct); 0161 % setenv('CUBACORES', '2'); 0162 tic 0163 out = integrateCubaDivonne(fhandle, xmin, xmax, optStruct); 0164 out.time = toc; 0165 0166 case 'cuhre' 0167 generateMexFile('integrateCubaCuhre', optStruct); 0168 % setenv('CUBACORES', '2'); 0169 tic 0170 out = integrateCubaCuhre(fhandle, xmin, xmax, optStruct); 0171 out.time = toc; 0172 otherwise 0173 error(['The specified algorithm ''' algo '''is not supported.']); 0174 end 0175 0176 %% general postprocessing 0177 out.relError = out.error ./ out.value; 0178 if out.fail ~= 0 0179 warning('elk:externCuba:integrationFailed', ... 0180 'Integration failed, likely because tolerance cannot be reached'); 0181 end 0182 if any(out.prob > 0.95) 0183 warning('elk:externCuba:integrationFailed', ['You might not want to ' ... 0184 'trust the given error bounds. See the Cuba documentation, ' ... 0185 'concerning the prob output']); 0186 end 0187 0188 end 0189 0190 function generateMexFile(file, optStruct) 0191 % identify my location 0192 originalDirectory = cd; 0193 myDirectory = fileparts(which(mfilename)); 0194 0195 % nothing to do, if mexfile already exists 0196 if exist([myDirectory filesep file '.' mexext], 'file') 0197 return 0198 end 0199 0200 % build cuba library 0201 if isunix 0202 try 0203 %% build cuba library object 0204 disp('--- BUILD CUBA LIRARY ---'); 0205 cd([myDirectory filesep optStruct.cubaDirectory]) 0206 % makefile should be adopted to use gcc4.1 (gcc4.2 throws 0207 % errors for divonne, gcc4.3 is not yet supported by MATLAB 0208 % and might have issues) 0209 0210 % configure 0211 eval(['!./configure CC=''' optStruct.gccCall '''']); 0212 !make lib 0213 !cp cuba.h ../ 0214 !cp libcuba.a ../ 0215 0216 %% build mexfile 0217 disp('--- BUILD MEXFILE ---'); 0218 disp(['Name: ' file '.' mexext]); 0219 cd(myDirectory) 0220 eval(['mex CC="' optStruct.gccCall '" ' file '.c libcuba.a']); 0221 0222 % restore original path 0223 cd(originalDirectory); 0224 0225 disp('--- BUILD DONE ---'); 0226 catch thisError 0227 % ensure we are back in the original directory 0228 cd(originalDirectory); 0229 rethrow(thisError); 0230 end 0231 %! DEV: another operating system might be added, via: 0232 % elseif ispc 0233 % % .. windows specific code 0234 % elseif ismac 0235 % % .. mac specific code 0236 %/ 0237 0238 else 0239 error(['I cannot automatically build the cuba library from' ... 0240 'source for a non-unix operating system. Please ' ... 0241 'adopt the implementation to your system.']); 0242 end 0243 end