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.

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 to rewrite it to make it fastetr?

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=%dn",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);
disp(result);