5
\$\begingroup\$

I was trying to rewrite the Python code in MATLAB.

The result is consistent.

enter image description here

But, the MATLAB code is so slow. Any help would be appreciated.

Ref python code link

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);

enter image description here

\$\endgroup\$
0

    2 Answers 2

    Reset to default
    6
    \$\begingroup\$

    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.

    \$\endgroup\$
      2
      \$\begingroup\$

      Avoid Symbolic Math for Numerical Evaluations:

      • Replace symbolic computations with MATLAB's numeric gamma function. 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
      \$\endgroup\$

        Start asking to get answers

        Find the answer to your question by asking.

        Ask question

        Explore related questions

        See similar questions with these tags.