I was trying to rewrite the Python code in MATLAB.
The result is consistent.
But, the MATLAB code is so slow. Any help would be appreciated.
The MATLAB code written by me is as follows: so slow, whereas definitely correct.
I guess the critical time-costing step is compMultiFoxHIntegrand function. How can I rewrite it to make it faster?
function result = compMultiFoxH(params, nsubdivisions, boundaryTol)
if nargin < 3
boundaryTol = 0.0001;
end
boundaries = detBoundaries(params, boundaryTol);
dim = numel(boundaries);
signs = dec2bin(0:2^dim-1) - '0';
signs(signs == 0) = -1;
signs=signs.*(-1);
inputs = repmat({(0:floor(nsubdivisions/2)-1)}, 1, dim);
code = cartesian(inputs{:});
quad = 0;
% res = [];
for i = 1:size(signs, 1)
fprintf("i=%d\n",i);
points = signs(i, :) .* (bsxfun(@plus, code, 0.5)) .* (boundaries * 2 / nsubdivisions);
integrandVals = (compMultiFoxHIntegrand(points, params));
% res = [res; integrandVals];
quad = quad + sum(integrandVals);
end
volume = prod(2 * boundaries / nsubdivisions);
result = quad * volume;
end
function boundaries = detBoundaries(params, tol)
boundary_range = 0:0.05:50;
dims = numel(params{1});
boundaries = zeros(1, dims);
for dim_l = 1:dims
points = zeros(length(boundary_range), dims);
points(:, dim_l) = boundary_range;
abs_integrand = abs(compMultiFoxHIntegrand(points, params));
index = find(abs_integrand > tol * abs_integrand(1), 1, 'last');
boundaries(dim_l) = boundary_range(index);
end
end
function result = compMultiFoxHIntegrand(y, params)
z = params{1};
mn = params{2};
pq = params{3};
c = params{4};
d = params{5};
a = params{6};
b = params{7};
m = cellfun(@(x) x(1), mn);
n = cellfun(@(x) x(2), mn);
p = cellfun(@(x) x(1), pq);
q = cellfun(@(x) x(2), pq);
npoints = size(y, 1);
dims = size(y, 2);
s = 1j * y;
lower = zeros(1, dims);
upper = zeros(1, dims);
for dim_l = 1:dims
if ~isempty(b{dim_l})
bj = cellfun(@(x) x(1), b{dim_l}(1:m(dim_l+1)));
Bj = cellfun(@(x) x(2), b{dim_l}(1:m(dim_l+1)));
lower(dim_l) = -min(bj ./ Bj);
else
lower(dim_l) = -100;
end
if ~isempty(a{dim_l})
aj = cellfun(@(x) x(1), a{dim_l}(1:n(dim_l+1)));
Aj = cellfun(@(x) x(2), a{dim_l}(1:n(dim_l+1)));
upper(dim_l) = min((1 - aj) ./ Aj);
else
upper(dim_l) = 1;
end
end
mindist = norm(upper - lower);
sigs = 0.5 * (upper + lower);
for j = 1:n(1)
num = 1 - c{j}(1) - sum(c{j}(2:end) .* lower);
cnorm = norm(c{j}(2:end));
newdist = abs(num) / (cnorm + eps);
if newdist < mindist
mindist = newdist;
sigs = lower + 0.5 * num * c{j}(2:end) / (cnorm * cnorm);
end
end
s = bsxfun(@plus, s, sigs);
s1 = [ones(npoints, 1), s];
prod_gam_num = 1 + 0j;
prod_gam_denom = 1 + 0j;
for j = 1:n(1)
prod_gam_num = prod_gam_num .* gamma(sym(1 - s1 * c{j}'));
end
for j = 1:q(1)
prod_gam_denom = prod_gam_denom .* gamma(sym(1 - s1 * d{j}'));
end
for j = n(1)+1:p(1)
prod_gam_denom = prod_gam_denom .* gamma(sym(s1 * c{j}'));
end
for dim_l = 1:dims
for j = 1:n(dim_l+1)
prod_gam_num = prod_gam_num .* gamma(sym(1 - a{dim_l}{j}(1) - a{dim_l}{j}(2) * s(:, dim_l)));
end
for j = 1:m(dim_l+1)
prod_gam_num = prod_gam_num .* gamma(sym(b{dim_l}{j}(1) + b{dim_l}{j}(2) * s(:, dim_l)));
end
for j = n(dim_l+1)+1:p(dim_l+1)
prod_gam_denom = prod_gam_denom .* gamma(sym(a{dim_l}{j}(1) + a{dim_l}{j}(2) * s(:, dim_l)));
end
for j = m(dim_l+1)+1:q(dim_l+1)
prod_gam_denom = prod_gam_denom .* gamma(sym(1 - b{dim_l}{j}(1) - b{dim_l}{j}(2) * s(:, dim_l)));
end
end
zs = z .^ -s;
result = double((prod_gam_num ./ prod_gam_denom) .* prod(zs, 2) ./ (2 * pi)^dims);
end
function C = cartesian(varargin)
args = varargin;
n = nargin;
[F{1:n}] = ndgrid(args{:});
for i=n:-1:1
G(:,i) = F{i}(:);
end
C = unique(G , 'rows');
end
clear;clc;
% Example usage
params1 = {...
[16.2982237081499, 16.2982237081499, 16.2982237081499, 16.2982237081499], ...
{[0, 0], [2, 1], [2, 1], [2, 1], [2, 1]}, ...
{[0, 1], [1, 2], [1, 2], [1, 2], [1, 2]}, ...
{}, ...
{[0, 1, 1, 1, 1]}, ...
{{[1, 2]}, {[1, 2]}, {[1, 2]}, {[1, 2]}}, ...
{{[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}} ...
};
result = compMultiFoxH(params1, 20);
format longG
disp(result);
2 Answers 2
There are some obvious inefficiencies in the code that I can point out.
The profiler shows that most time is spent doing symbolic math. You use symbolic math to compute the product of a series of gamma functions. There is no reason to compute this symbolically. MATLAB has a numeric gamma function. Simply leaving out the call to sym in your code should reduce the computational cost of your code by a few orders of magnitude.
Symbolic math is useful in a few circumstances. Mostly when doing symbolic manipulation, for example when finding the zeros of a function, which yields a new equation that can be evaluated many times over. The complicated part is done once, then the numerical evaluation that needs to be done millions of times becomes cheaper. In your case, all you do with the symbolic math is evaluate a numerical result, there is no symbolic manipulation. Numerical evaluation should be done numerically, it is a lot cheaper that way.
The other obvious inefficiency is all the unnecessary calls to cellfun. If you define the params data structure to be more convenient, then you won't need to copy and move around all these values every time you call the inner function compMultiFoxHIntegrand.
For example,
mn = params{2};
% ...
m = cellfun(@(x) x(1), mn);
n = cellfun(@(x) x(2), mn);
can be simplified by defining params{2} to be a 2D array instead of a cell array of 1D arrays. Your current implementation is an obvious 1-to-1 translation of the Python list of lists, but this is not very efficient at all. It is also not efficient in Python, the original code would also be more efficient using a 2D array.
mn = params{2};
m = mn(:, 1);
n = mn(:, 2);
% ...
params1 = {...
% ...
[0, 0; 2, 1; 2, 1; 2, 1; 2, 1], ...
% ...
};
Also, since you don't use mn in any other way but to separate it into its columns (at least not in the code you show), you could define params to have the two arrays separately. Note that m = mn(:, 1) means the data is being copied over.
If you use a struct instead of a cell array, this would be more convenient too:
params1 = [];
params1.m = [0, 2, 2, 2, 2];
params1.n = [0, 1, 1, 1, 1];
When specific indices into a cell array have specific meanings, it becomes very easy to mix things up. Making a change to the definition of params (e.g. by adding a value) could mean having to update all the code that uses it, and making mistakes there is very easy. A struct in MATLAB is very much like a cell array, but the indexing is done with a name instead of a number. This makes it much less likely that indices get mixed up, makes the code easier to read (e.g. you don't need to do mn = params{2}; just to give the value a readable name in the remainder of the code, you can just use params.mn everywhere, this indexing is very cheap), and makes it trivial to make changes to the definition of params.
Avoid Symbolic Math for Numerical Evaluations:
- Replace symbolic computations with MATLAB's numeric
gammafunction. Symbolic math should only be used for symbolic manipulations, not for numerical evaluations.
It does work. Now it's much faster.
clear;clc;
% Example usage
params1 = {...
[16.2982237081499, 16.2982237081499, 16.2982237081499, 16.2982237081499], ...
{[0, 0], [2, 1], [2, 1], [2, 1], [2, 1]}, ...
{[0, 1], [1, 2], [1, 2], [1, 2], [1, 2]}, ...
{}, ...
{[0, 1, 1, 1, 1]}, ...
{{[1, 2]}, {[1, 2]}, {[1, 2]}, {[1, 2]}}, ...
{{[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}, {[1, 0.6666666666666666], [3.5, 0.5]}} ...
};
result = compMultiFoxH(params1, 20);
format longG
disp(result);
function result = compMultiFoxH(params, nsubdivisions, boundaryTol)
if nargin < 3
boundaryTol = 0.0001;
end
boundaries = detBoundaries(params, boundaryTol);
dim = numel(boundaries);
signs = dec2bin(0:2^dim-1) - '0';
signs(signs == 0) = -1;
signs=signs.*(-1);
inputs = repmat({(0:floor(nsubdivisions/2)-1)}, 1, dim);
code = cartesian(inputs{:});
quad = 0;
% res = [];
for i = 1:size(signs, 1)
% fprintf("i=%d\n",i);
points = signs(i, :) .* (bsxfun(@plus, code, 0.5)) .* (boundaries * 2 / nsubdivisions);
integrandVals = (compMultiFoxHIntegrand(points, params));
% res = [res; integrandVals];
quad = quad + sum(integrandVals);
end
volume = prod(2 * boundaries / nsubdivisions);
result = quad * volume;
end
function boundaries = detBoundaries(params, tol)
boundary_range = 0:0.05:50;
dims = numel(params{1});
boundaries = zeros(1, dims);
for dim_l = 1:dims
points = zeros(length(boundary_range), dims);
points(:, dim_l) = boundary_range;
abs_integrand = abs(compMultiFoxHIntegrand(points, params));
index = find(abs_integrand > tol * abs_integrand(1), 1, 'last');
boundaries(dim_l) = boundary_range(index);
end
end
function result = compMultiFoxHIntegrand(y, params)
z = params{1};
mn = params{2};
pq = params{3};
c = params{4};
d = params{5};
a = params{6};
b = params{7};
m = cellfun(@(x) x(1), mn);
n = cellfun(@(x) x(2), mn);
p = cellfun(@(x) x(1), pq);
q = cellfun(@(x) x(2), pq);
npoints = size(y, 1);
dims = size(y, 2);
s = 1j * y;
lower = zeros(1, dims);
upper = zeros(1, dims);
for dim_l = 1:dims
if ~isempty(b{dim_l})
bj = cellfun(@(x) x(1), b{dim_l}(1:m(dim_l+1)));
Bj = cellfun(@(x) x(2), b{dim_l}(1:m(dim_l+1)));
lower(dim_l) = -min(bj ./ Bj);
else
lower(dim_l) = -100;
end
if ~isempty(a{dim_l})
aj = cellfun(@(x) x(1), a{dim_l}(1:n(dim_l+1)));
Aj = cellfun(@(x) x(2), a{dim_l}(1:n(dim_l+1)));
upper(dim_l) = min((1 - aj) ./ Aj);
else
upper(dim_l) = 1;
end
end
mindist = norm(upper - lower);
sigs = 0.5 * (upper + lower);
for j = 1:n(1)
num = 1 - c{j}(1) - sum(c{j}(2:end) .* lower);
cnorm = norm(c{j}(2:end));
newdist = abs(num) / (cnorm + eps);
if newdist < mindist
mindist = newdist;
sigs = lower + 0.5 * num * c{j}(2:end) / (cnorm * cnorm);
end
end
s = bsxfun(@plus, s, sigs);
s1 = [ones(npoints, 1), s];
prod_gam_num = 1 + 0j;
prod_gam_denom = 1 + 0j;
for j = 1:n(1)
prod_gam_num = prod_gam_num .* gammac( (1 - s1 * c{j}'));
end
for j = 1:q(1)
prod_gam_denom = prod_gam_denom .* gammac( (1 - s1 * d{j}'));
end
for j = n(1)+1:p(1)
prod_gam_denom = prod_gam_denom .* gammac( (s1 * c{j}'));
end
for dim_l = 1:dims
for j = 1:n(dim_l+1)
prod_gam_num = prod_gam_num .* gammac( (1 - a{dim_l}{j}(1) - a{dim_l}{j}(2) * s(:, dim_l)));
end
for j = 1:m(dim_l+1)
prod_gam_num = prod_gam_num .* gammac( (b{dim_l}{j}(1) + b{dim_l}{j}(2) * s(:, dim_l)));
end
for j = n(dim_l+1)+1:p(dim_l+1)
prod_gam_denom = prod_gam_denom .* gammac( (a{dim_l}{j}(1) + a{dim_l}{j}(2) * s(:, dim_l)));
end
for j = m(dim_l+1)+1:q(dim_l+1)
prod_gam_denom = prod_gam_denom .* gammac( (1 - b{dim_l}{j}(1) - b{dim_l}{j}(2) * s(:, dim_l)));
end
end
zs = z .^ -s;
result = double((prod_gam_num ./ prod_gam_denom) .* prod(zs, 2) ./ (2 * pi)^dims);
end
function C = cartesian(varargin)
args = varargin;
n = nargin;
[F{1:n}] = ndgrid(args{:});
for i=n:-1:1
G(:,i) = F{i}(:);
end
C = unique(G , 'rows');
end
function [f] = gammac(z)
% GAMMA Gamma function valid in the entire complex plane.
% Accuracy is 15 significant digits along the real axis
% and 13 significant digits elsewhere.
% This routine uses a superb Lanczos series
% approximation for the complex Gamma function.
%
% z may be complex and of any size.
% Also n! = prod(1:n) = gamma(n+1)
%
%usage: [f] = gamma(z)
%
%tested on versions 6.0 and 5.3.1 under Sun Solaris 5.5.1
%
%References: C. Lanczos, SIAM JNA 1, 1964. pp. 86-96
% Y. Luke, "The Special ... approximations", 1969 pp. 29-31
% Y. Luke, "Algorithms ... functions", 1977
% J. Spouge, SIAM JNA 31, 1994. pp. 931-944
% W. Press, "Numerical Recipes"
% S. Chang, "Computation of special functions", 1996
% W. J. Cody "An Overview of Software Development for Special
% Functions", 1975
%
%see also: GAMMA GAMMALN GAMMAINC PSI
%see also: mhelp GAMMA
%
%Paul Godfrey
%[email protected]
%http://winnie.fit.edu/~gabdo/gamma.txt
%Sept 11, 2001
siz = size(z);
z=z(:);
zz=z;
f = 0.*z; % reserve space in advance
p=find(real(z)<0);
if ~isempty(p)
z(p)=-z(p);
end
% 15 sig. digits for 0<=real(z)<=171
% coeffs should sum to about g*g/2+23/24
g=607/128; % best results when 4<=g<=5
c = [ 0.99999999999999709182;
57.156235665862923517;
-59.597960355475491248;
14.136097974741747174;
-0.49191381609762019978;
.33994649984811888699e-4;
.46523628927048575665e-4;
-.98374475304879564677e-4;
.15808870322491248884e-3;
-.21026444172410488319e-3;
.21743961811521264320e-3;
-.16431810653676389022e-3;
.84418223983852743293e-4;
-.26190838401581408670e-4;
.36899182659531622704e-5];
%Num Recipes used g=5 with 7 terms
%for a less effective approximation
z=z-1;
zh =z+0.5;
zgh=zh+g;
%trick for avoiding FP overflow above z=141
zp=zgh.^(zh*0.5);
ss=0.0;
for pp=size(c,1)-1:-1:1
ss=ss+c(pp+1)./(z+pp);
end
%sqrt(2Pi)
sq2pi= 2.5066282746310005024157652848110;
f=(sq2pi*(c(1)+ss)).*((zp.*exp(-zgh)).*zp);
f(z==0 | z==1) = 1.0;
%adjust for negative real parts
if ~isempty(p)
f(p)=-pi./(zz(p).*f(p).*sin(pi*zz(p)));
end
%adjust for negative poles
p=find(round(zz)==zz & imag(zz)==0 & real(zz)<=0);
if ~isempty(p)
f(p)=Inf;
end
f=reshape(f,siz);
return
endfunction [outputArg1,outputArg2] = untitled2(inputArg1,inputArg2)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
outputArg1 = inputArg1;
outputArg2 = inputArg2;
end