function varargout = mcint(h, n, opt)
%MCINT Monte Carlo integrator for multiple functions
%
% Satisfy your wildest integration fantasies!
%
%
%
% written by Lee Ferchoff, 2006
% Do not modify this file in any way unless for your own personal
% use. Do not re-release any part of the code contained in this file.
% Do not attach your own name to any part of the code contained in this
% file. Please reference me and e-mail to inform me if you use any part
% of this code, or original ideas from the algorithm. Thanks! Hope you find
% it helpful.
% E-mail: umferch1 at cc DOT umanitoba DOT ca
%
%
%
%
% A real-valued integral is calculated in arbitrary n-dimensional
% coordinates by averaging the value of a
% function over a large number of randomly selected points within the
% hypervolume to be integrated. Multiple functions can be integrated
% simultaneously over the same domain. The hypervolume can be of
% arbitrary shape, as long as it can be expressed as a series of
% logical conditions on the coordinates.
%
% Examples of the usage of mcint can be found in learnmcint.m
%
% Throughout, x(i,j) contains the matrix of "darts", or random
% points used to evaluate the function to be integrated.
% The rows i, refer to the ith coordinate of a
% point, where each column j, is a separate point. In standard
% notation, the coordinate vectors of all points 'n' are given by:
%
% 1-coordinate vector is x(1,:)
% 2-coordinate vector is x(2,:)
% 3-coordinate vector is x(3,:)
% ...
% mth-coordinate is x(m,:)
%
% In Cartesian coordinates, you might define the 1-dimension to be x,
% the 2-dimension to be y, or if you are using spherical coordinates,
% you could define the 1-dimension to be r, the 2-dimension to be
% theta, and the 3-dimension to be phi.
%
% This is the format that must be used when passing functions
% of coordinates to mcint. Operations on these coordinates must
% be written as vector operations (i.e. .* and .^, not * and ^).
%
% To integrate in a non-Cartesian coordinate system, in your
% function definition, multiply by the Jacobian from
% the Cartesian to 'system' transformation, for whatever coordinate
% system you are using. Jacobians for common coordinate systems
% can be obtained from jacobian.m. For example, to integrate
% f(r,theta,phi) in spherical coordinates, on paper you would write
%
% triple integral ( f(r,theta,phi) * r^2 * sin(theta) dr d(theta) d(phi) )
%
% So the integrand should be defined as:
%
% "f(x)" .* x(1,:).^2 + sin(x(2,:))
% or making use of jacobian.m;
% "f(x)" .* jacobian(x,'spherical')
%
% You can write your own Jacobian in, if you choose not to call jacobian.m.
% This flexibility allows you to mix coordinate systems. For example, you
% can do a 6-dimensional integral where three of the dimensions are in
% spherical coordinates and the other three are in elliptical. Furthermore,
% you can choose whatever order you want for the coordinate
% defintions, as long as you write the Jacobian accordingly. For example,
% x(1,:) could be theta, x(2,:) r, x(3,:) the angular variable in the
% elliptic coordinates, and so on.
%
%
%INPUT
%
% h -- data structure defining a hypervolume integration problem
%
% h.b(i,j) -- bounding parallepiped of the hypervolume to be
% integrated.
% size(b) = [dim 2] where dim is the number of dimensions
% of the hypervolume. (i,1) is the minimum bounding coordinate
% in the ith dimension and (i,2) is the maximum bounding coordinate
% in the ith dimension.
%
% Example: To integrate x from -1 to 1, and y from 0 to Inf,
% b(:,1) = [-1; 0]
% b(:,2) = [1; Inf]
%
% h.cond{1,j} (OPTIONAL) -- cell array which contains the functions
% expressing additional logical conditional statements that define
% a hypervolume within its bounding box. These conditional
% statements, indexed by j, are functions.
% h.cond may be empty, if one wishes to integrate
% over a rectangular parallepiped and doesn't
% require any further logical conditions on the
% hypervolume.
%
% Example: cond = {@outCircle @aboveLine}
% where the functions are defined as:
%
% function tf = outcircle(x, par)
% tf = x(1,:).^2 + x(2,:).^2 >= 1;
%
% function tf = aboveLine(x, par)
% tf = x(2,:) > x(1,:);
%
% In this example, our hypervolume of integration
% is restricted to the region outside of circle of
% radius 1 at the origin, and to the upper left
% of the line y=x.
%
% All h.cond function definitions should take a structure par as
% their second input, even if no parameters from 'par' are used in
% the funciton.
%
% h.funcs{1,j} -- cell array which contains the functions to be
% integrated. When the jth function is evaluated, it returns
% the evaluation of the jth integrand. The usual usage is
% to enter a function that returns an row vector for the integrand
% f(1,k), where the kth column is the integrand evaluated at the
% kth dart point. mcint tests the conditions in the order you enter
% them in the cell array, therefore it is most efficient to list
% them from most likely to fail to least likely to fail.
%
% Example: function f = simpleExample(x, par)
% f = x(1,:).^2 .* 2*exp(-x(2,:));
%
% In this case, the integral returned will be a scalar, I(1,1)
% containing the value of the integral.
%
% If h.funcs contains several functions, I(1,m) will be the
% value of the integral for the mth function in h.funs.
%
% Example: function f = simpleExample(x, par)
% f = x(1,:).^2 .* 2*exp(-x(2,:));
% function f = Id(x)
% f = ones( 1, size(x,2));
%
% In this example, I(1,1) will be the integral of @simpleExample
% and I(1,2) will be the integral of @Id. Note that @Id is the
% identity function; to obtain the hypervolume itself, simply
% integrate this identity function.
%
% To allow greater flexibility, a function handle may contain
% more than one function. If h.funcs contains one function handle
% which returns a function with m rows, then I(m,1) will be
% the integral for the integrand defined by the mth row returned
% by the function.
%
% Example: function f = energy(x, par)
% p = calculatePsi(x);
% H = getHamiltonian(x);
% f(1,:) = conj(p).*H.*p;
% f(2,:) = conj(p).*p;
%
% This example shows you how you might calculate the expectation
% value of energy for a quantum system. Here we've calculated
% p (the wavefunction psi) and H (the Hamiltonian) on the random
% dart matrix (x). Row one of the function will be the
% evaluation of the integrand for the expectation value of the
% Hamiltonian, and row two will be the probability density.
% mcint will give I(1,1) which is the expectation value of the
% Hamiltonian, and I(2,1) which is the total probability. Then
% the user can simply find the expectation value of energy as
% I(1,1)/I(2,1).
%
% This method of putting multiple integrands in one function
% handle is preferable when they are calculated with the same
% ingredients. It would be more time consuming to make a separate
% function handle for both expectation values calculated above
% since the wavefunction psi would have to be calculated in
% BOTH of the function handles with a call to @calculatePsi:
% p = calculatePsi(x);
%
% This inefficiency can thus be avoided. In Quantum Energy Solver
% (Lee Ferchoff), I use this method to calculate the expectation
% value of the Hamiltonian, and the identity operator for ALL of
% the individuals in the population at at once under one function
% handle.
%
% NOTE: As of now, if multiple function handles are used, each one
% must contain the same number of rows. Support will be added
% for varying numbers of rows in a future version.
%
% In general, if h.funcs contains j function handles, then mcint
% returns I, where I(i,j) is the integral of the function in the ith
% row of the jth function handle.
%
% All h.funcs function definitions should take a structure 'par' as
% their second input, even if no parameters from 'par' are used in
% the funciton.
%
%
% n -- the number of randomly chosen sampling points
%
%
% opt (OPTIONAL) -- integration options
%
% opt.maxArray -- adjustable max array size parameter. mcint has
% built-in memory management. This is the largest number
% of elements that will be stored in memory for
% use by mcint at any given time.
%
% NOTE: The default max array size of 1e4 has been tested
% extensively, this gives the shortest run time for a fixed amount
% of points. This has been tested for integrals of dimensions 1 and
% 3 and n of 1e5 and 5e6, and it was found that run time increased
% if the max array size was either increased or decreased from 1e4
% in all cases.
% If it is decreased from 1e4, the run time sharply increases. This
% is because we lose most of the vectorization of the code,
% since the function evaluations are broken up into vectorized
% groups of the maximum array size. If the max array size is
% increased from 1e4, the run time gradually increases by about
% 25% until memory management is not being used and the whole
% integration is done in one shot. I suppose this is because
% doing stuff with really large matricies is unwieldy and in
% this case breaking up a bit of the vectorization at the benefit
% of using smaller matricies seems to be a good thing.
%
% 1e-4 is optimum for the linux boxes in Room 523.
% Anything above 1e-3 seems to be handled equally well by
% hactar.
% opt.noerror -- if true, statistical errors are not calculated and
% only the integral value is returned. Note that this mode cannot
% be used when a tolerance is requested with opt.tol.
% opt.tol -- TOLERANCE STOPPING CONDITION
% absolute tolerance requested (integration will stop
% when reached)
% opt.tol(i,j) is the absolute tolerance
% for the ith row of the jth function in h.funcs. If a tolerance
% is set to zero, then it is taken to mean that no tolerance is
% requested for that particular integral.
% opt.time -- TIME STOPPING CONDITION
% maximum time allowed in seconds. Integration
% will end normally after this time.
% opt.warnOff -- if true, turn off printed warnings, such as not meeting a
% requested tolerance
% opt.mcintPar -- a struture containing parameters that is passed to the
% functions in h.funcs, which can be used in function evaluation
%
%OUTPUT
%
% varargout(1) is
% I(i,j,k) -- integral (with optional error)
% I(i,j,1) -- value of the integral for the integrand defined
% in the ith row of the jth function handle in h.funcs.
% I(i,j,2) -- value of the statistical error of the integral
% for the integrand defined in the ith row of the jth function
% handle in h.funcs. The error is calculated from random fluctations
% of the integrand about its mean value.
%
% Errors are calculated by default, so k=1:2 by default.
% If user passes opt.noerror=1 to mcint, then no errors are
% calculated, and so k=1 only.
%
% varargout(2) is
% info (OPTIONAL) -- technical data about the integration. The data
% can be used for error checking.
%
% info.maxArray -- maximum array size used
% info.memoryManaged -- if true, the requested number of points
% n was large enough, that the integration had to be
% performed with chunks of size maxArray
% info.noerror -- if true, no statistical error calculated
% info.extraParameters -- if true, extra parameters were passed
% to the functions to be integrated in opt.mcintPar
%
% info.dim -- number of dimensions
% info.n -- number of random points requested
% info.points -- number of random points used
% info.pointdiff -- info.n-info.points
%
% info.tolReq(i,j) -- tolerance requested by user. The abs tolerance
% for the ith row of the jth function in h.funcs. If a tolerance
% is set to zero, then it is taken to mean that no tolerance is
% requested for that particular integral. This is what the user
% entered in opt.tol.
% info.someTolReq -- if true, at least one integral tolerance has
% been requested
% info.tolMet(i,j) -- true in the indicies for integrals which
% have had their requested tolerance met. (i,j) labels the ith
% row of the jth function in h.funcs
% info.tolMet=-1 in the indicies representing integrals for which
% a specific tolerance was not requested. If no tolerance was
% required for any integral, info.tolMet is a scalar -1.
% info.allTolMet -- if true, the requested tolerance for every
% tolerance was met
% info.timeReq -- max time requested by the user. info.timeReq=-1
% if a specific time limit was not requested
% info.timeMet -- if true, the integration was within the requested
% time. info.timeMet=-1 if a specific time limit was not
% requested
% info.totalTime -- total time required for integration
%
%
%FEATURES TO BE ADDED IN THE NEXT VERSION:
%
% - relative tolerance stopping condition
% - differening amounts of rows in the functions in h.funcs
% - protection against overflow from a chance hit near a singularity
%
%Lee Ferchoff
t0 = clock; % record initial time
if ~isfield(h,'cond')
h.cond = {};
end
% ===========================
% SET OPTIONS
% ===========================
% --------------------------------------
% mcint is memory managed, and will not
% use an array that has more elements
% than the number maxArray
p.info.maxArray = 1e4; % default max array size
% will be changed if opt.maxArray is a field
% --------------------------------------
% -1 means no request
p.info.tolReq = -1;
p.info.someTolReq = -1;
p.info.tolMet = -1;
p.info.allTolMet = -1;
p.info.timeReq = -1;
p.info.timeMet = -1;
p.info.noerror = 0; % if true, statistical errors are not calculated
p.info.extraParameters = 0; % if true, the user has provided opt.mcintPar
warnOff = 0; % if true, no warnings will be printed to screen
p.par = [];
if nargin==3 % optional parameter structure opt has been provided
if isfield(opt,'maxArray')
p.info.maxArray = opt.maxArray;
end
if isfield(opt,'noerror')
p.info.noerror = opt.noerror;
end
if isfield(opt,'mcintPar') % user has provided extra parameters
p.info.extraParameters = 1;
p.par = opt.mcintPar;
end
if isfield(opt,'tol') % user requested tolerances
p.info.tolReq = abs(opt.tol);
p.info.someTolReq = 1;
p.info.tolMet = 0;
p.info.allTolMet = 0;
end
if isfield(opt,'time') % user requested max integration time
p.info.timeReq = opt.time;
p.info.timeMet = 0;
end
if isfield(opt,'warnOff')
warnOff = opt.warnOff;
end
end
% error: can't not calculate error and still request an error tolerance
if p.info.noerror==1 && p.info.someTolReq==1
if warnOff==0
disp('WARNING: IF a tolerance is requested, error must be calculated.');
disp('Integration will continue WITH error calculation.');
end
p.info.noerror = 0;
end
% ===========================
% SETUP INTEGRATION PROBLEM
% ===========================
p.dim = size(h.b, 1);
p.nf = length(h.funcs); % # functions
points = n; % running counter that keeps track of how many of the n darts have been thrown
% Transform any improper integral bounds:
% true in rows of b having at least one improper integral limit
p.infBound = sum(isinf(h.b),2) ~= 0;
% perform inverse tangent coordinate transformation on affected rows.
% This will allow us to calculate the improper integral in the old
% coordinate system as a proper integral in the new coordinate system.
h.b(p.infBound,:) = atan(h.b(p.infBound,:));
% p.infDims isthe numbers of the dimensions with improper limits
p.infDims = find(p.infBound==1);
p.improperProblem = ~isempty(p.infDims); % if true, at least one infinite dim
V = prod(h.b(:,2) - h.b(:,1)); % bounding volume (after transformation)
% ===========================
% INTEGRATE
% ===========================
if p.dim*n <= p.info.maxArray
p.info.memoryManaged = 0;
p.cols = n;
p.xMin=repmat(h.b(:,1),1,p.cols);
p.xMax=repmat(h.b(:,2),1,p.cols);
s = updateSums(h,p);
else % memory management needed
p.info.memoryManaged = 1;
p.cols = floor(p.info.maxArray/p.dim);
p.xMin=repmat(h.b(:,1),1,p.cols);
p.xMax=repmat(h.b(:,2),1,p.cols);
% # of runs which use the full maxArray size
fullruns = ceil(p.dim*n/p.info.maxArray)-1;
for j=1:fullruns
newS = updateSums(h,p);
if ~exist('s') % first run through loop; initialize s, s2
s = zeros(size(newS));
end
s = s + newS;
points = j*p.cols;
% Check for stopping conditions
if (p.info.timeReq~=-1 && etime(clock,t0)>p.info.timeReq)
% executes if max time has elapsed
p.info.timeMet = 1;
break;
elseif p.info.someTolReq==1
% checks if all tolerances have been met
E = getError(s,points,V);
tolMatrix = zeros(size(E));
tolMatrix(:) = Inf;
tolMatrix = p.info.tolReq(find(p.info.tolReq~=0));
p.info.tolMet = E < tolMatrix;
p.info.allTolMet = prod(+p.info.tolMet(:));
if p.info.allTolMet==1
break;
end
end
end
% if integration wasn't stopped by some condition, finish the leftover
% points
if ~(p.info.timeMet==1 || p.info.allTolMet==1)
% The last update covers the leftover points on the exact
% number of points the user requested
p.cols = n - points;
p.xMin=repmat(h.b(:,1),1,p.cols);
p.xMax=repmat(h.b(:,2),1,p.cols);
newS = updateSums(h,p);
s = s + newS;
points = points + p.cols;
end
end
% ===========================
% REPORT RESULT
% ===========================
% return integral and possibly error
I(:,:,1) = V/points * s(:,:,1); % integral
if p.info.noerror==0
I(:,:,2) = getError(s,points,V); % error
end
varargout(1) = {I};
% some tolerance was requested and wasn't met
if p.info.allTolMet==0 && warnOff==0
disp('***');
disp('Max number of requested darts thrown.');
disp('Requested tolerance was not achieved.');
disp('***');
end
% build info structure to pass to user
if nargout==2
p.info.dim = p.dim;
p.info.n = n;
p.info.points = points;
p.info.pointdiff = n - points;
p.info.totalTime = etime(clock,t0);
varargout(2) = {p.info};
end
function E = getError(s,points,V)
%GETERROR Find the current statistical error of the integral
%
%INPUT
% s -- s(:,:,1) function sums
% -- s(:,:,2) function squared sums
% points -- number of points in sum
% V -- hypervolume
%
%OUTPUT
%
% E(i,j) -- the current error of the ith row and jth column using the
% function sums in s
E = V/points * sqrt(s(:,:,2) - s(:,:,1).*s(:,:,1)/points);
function s = updateSums(h, p)
%UPDATESUMS Calculate function sums for mcint
%
%INPUT
% h -- data structure defining the hypervolume integration problem
% p -- data structure containing the current parameters of the
% integration
%
%OUTPUT
% s(i,j,k) -- k=1, sum of the functions over all darts
% k=2, sum of the functions squared over all darts
% for the ith row of the jth function handle in h.funcs
%
% if p.info.noerror=1, then the k=2 matrix is not returned
% p.xMin, p.xMax passed into function rather than being calculated
% in it so that p.xMin, p.xMax only need to be calculated once
% for many memory blocks of size info.maxArray
% ===========================
% GENERATE RANDOM DARTS
% ===========================
% x(i,j) is the ith Cartesian coordinate of the jth random point
x = p.xMin + (p.xMax-p.xMin) .* rand(p.dim,p.cols); % random dart matrix
% ===========================
% COORDINATE TRANSFORMATIONS
% ===========================
% perform coordinate transformations if there is at least one improper
% integral. Transform only the infinite dimensions, since the
% transformation costs time and is not needed for the finite dimensions
if p.improperProblem
x(p.infDims,:) = tan( x(p.infDims,:) );
end
% remove all points from x which aren't in hypervolume
% by using the extra restrictions given on the hypervolume
for i=1:length(h.cond)
x = x( :, feval(h.cond{i},x));
end
% we calculate the Jacobian of the arctan transformation AFTER removing
% the points that aren't in the hypervolume so we don't waste time
% calculating the Jacobian for points that aren't in the integration
if p.improperProblem
jacobian = prod( x(p.infDims,:).^2 + 1, 1 );
end
% ===========================
% INTEGRAND EVALUATION
% ===========================
for i = p.nf: -1: 1
f = feval(h.funcs{i}, x, p.par);
if p.improperProblem % multiply by the Jacobian of the arctan transform
jacobian = repmat(jacobian(1,:), [size(f,1) 1]);
f = f .* jacobian;
end
s(:,i,1) = sum(f,2);
if p.info.noerror==0
s(:,i,2) = sum(f.*f,2);
end
end
