How to port Matlab/Python's multivariate FoxH implementation in Mathematica?

Question

Ref code

My implementation

Are there some mistakes in my Mathematica code? Any help would be greatly appreciated.

Attempt Version 1 (❌)

Clear["Global`*"];

(*Function definitions*)
compMultiFoxH[params_, nsubdivisions_, boundaryTol_ : 0.0001] := 
 Module[{boundaries, dim, signs, inputs, code, quad, volume},
  boundaries = detBoundaries[params, boundaryTol];
  Print["boundaries=", boundaries];
  dim = Length[boundaries];
  signs = Tuples[{-1, 1}, dim];
  code = Tuples[Range[0, Floor[nsubdivisions/2] - 1], dim];
  quad = 0;
  Do[Module[{points, integrandVals}, 
    points = signs[[i]]  (code + 0.5)  (boundaries*2/nsubdivisions);
    integrandVals = compMultiFoxHIntegrand[points, params];
    quad += Total[integrandVals];], {i, 1, Length[signs]}];
  volume = Times @@ (2*boundaries/nsubdivisions);
  quad*volume]

detBoundaries[params_, tol_] := 
 Module[{boundaryRange, dims, boundaries, points, absIntegrand, 
   index}, boundaryRange = Range[0, 50, 0.05];
  dims = Length[params[[1]]];
  boundaries = ConstantArray[0, dims];
  Do[points = ConstantArray[0, {Length[boundaryRange], dims}];
   points[[All, dimL]] = boundaryRange;
   absIntegrand = Abs[compMultiFoxHIntegrand[points, params]];
   index = 
    Last[Flatten[
      Position[absIntegrand, _?(# > tol*absIntegrand[[1]] &), {1}, 
       Heads -> False]]];
   boundaries[[dimL]] = boundaryRange[[index]];, {dimL, 1, dims}];
  boundaries]

compMultiFoxHIntegrand[y_, params_] := 
 Module[{z, mn, pq, c, d, a, b, m, n, p, q, npoints, dims, s, lower, 
   upper, mindist, sigs, num, cnorm, newdist, s1, prodGamNum, 
   prodGamDenom, zs}, {z, mn, pq, c, d, a, b} = params;
  m = mn[[All, 1]];
  n = mn[[All, 2]];
  p = pq[[All, 1]];
  q = pq[[All, 2]];
  npoints = Length[y];
  dims = Length[First[y]];
  s = I  y;
  lower = ConstantArray[0, dims];
  upper = ConstantArray[0, dims];
  Do[If[b[[dimL]] =!= {}, 
    Module[{bj, Bj}, bj = b[[dimL, All, 1]][[;; m[[dimL + 1]]]];
     Bj = b[[dimL, All, 2]][[;; m[[dimL + 1]]]];
     lower[[dimL]] = -Min[bj/Bj];], lower[[dimL]] = -100];
   If[a[[dimL]] =!= {}, 
    Module[{aj, Aj}, aj = a[[dimL, All, 1]][[;; n[[dimL + 1]]]];
     Aj = a[[dimL, All, 2]][[;; n[[dimL + 1]]]];
     upper[[dimL]] = Min[(1 - aj)/Aj];], upper[[dimL]] = 1];, {dimL, 
    1, dims}];
  mindist = Norm[upper - lower];
  sigs = 0.5  (upper + lower);
  Do[num = 1 - c[[j, 1]] - Total[c[[j, 2 ;;]]  lower];
   cnorm = Norm[c[[j, 2 ;;]]];
   newdist = Abs[num]/(cnorm + $MachineEpsilon);
   If[newdist < mindist, mindist = newdist;
    sigs = lower + 0.5  num  c[[j, 2 ;;]]/(cnorm^2);], {j, 1, n[[1]]}];
  s = s + sigs;
  s1 = Join[ConstantArray[1, {npoints, 1}], s, 2];
  prodGamNum = 1 + 0  I;
  prodGamDenom = 1 + 0  I;
  Do[prodGamNum *= Gamma[(1 - s1 . c[[j]])];, {j, 1, n[[1]]}];
  Do[prodGamDenom *= Gamma[(1 - s1 . d[[j]])];, {j, 1, q[[1]]}];
  Do[prodGamDenom *= Gamma[(s1 . c[[j]])];, {j, n[[1]] + 1, p[[1]]}];
  Do[Do[prodGamNum *= 
      Gamma[(1 - a[[dimL, j, 1]] - 
         a[[dimL, j, 2]]  s[[All, dimL]])];, {j, 1, n[[dimL + 1]]}];
   Do[prodGamNum *= 
      Gamma[(b[[dimL, j, 1]] + b[[dimL, j, 2]]  s[[All, dimL]])];, {j,
      1, m[[dimL + 1]]}];
   Do[prodGamDenom *= 
      Gamma[(a[[dimL, j, 1]] + a[[dimL, j, 2]]  s[[All, dimL]])];, {j,
      n[[dimL + 1]] + 1, p[[dimL + 1]]}];
   Do[prodGamDenom *= 
      Gamma[(1 - b[[dimL, j, 1]] - 
         b[[dimL, j, 2]]  s[[All, dimL]])];, {j, m[[dimL + 1]] + 1, 
     q[[dimL + 1]]}];, {dimL, 1, dims}];
  zs = z^-s;
  (prodGamNum/prodGamDenom)  Product[zs, {2}]/(2  Pi)^dims // N]

(*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];
Print[result];

Attempt Version 2 (❌)

(*Define functions for gamma product \
computation*)ClearAll["Global`*"]
gammaProdNum[s1_, c_, n_] := 
 Times @@ (Gamma[1 - c[[#]] . s1] & /@ Range[n[[1]]])
gammaProdDenom[s1_, c_, p_, q_] := 
 Times @@ (Gamma[1 - c[[#]] . s1] & /@ Range[q[[1]]])*
  Times @@ (Gamma[c[[#]] . s1] & /@ Range[n[[1]] + 1, p[[1]]])

(*Compute boundaries*)
detBoundaries[params_, tol_] := 
 Module[{boundaryRange, dims, boundaries, points, absIntegrand, 
   index}, boundaryRange = Range[0, 50, 0.05];
  dims = Length[params[[1]]];
  boundaries = ConstantArray[0, dims];
  Do[points = ConstantArray[0, {Length[boundaryRange], dims}];
   points[[All, dimL]] = boundaryRange;
   absIntegrand = Abs[compMultiFoxHIntegrand[points, params]];
   index = 
    Max[FirstPosition[
        UnitStep[absIntegrand - tol*First[absIntegrand]], 1][[1]] - 1];
   boundaries[[dimL]] = boundaryRange[[index]];, {dimL, dims}];
  boundaries]

(*Compute complex integrand of the multivariate Fox-H function*)
compMultiFoxHIntegrand[y_, params_] := 
 Module[{z, mn, pq, c, d, a, b, m, n, p, q, s, lower, upper, dims, s1,
    prodGamNum, prodGamDenom, zs, result}, {z, mn, pq, c, d, a, b} = 
   params;
  {m, n} = Transpose[mn];
  {p, q} = Transpose[pq];
  {npoints, dims} = Dimensions[y];
  s = I*y;
  lower = ConstantArray[0, dims];
  upper = ConstantArray[0, dims];
  Do[If[Length[b[[dimL]]] > 0, 
    lower[[dimL]] = -Min[b[[dimL, All, 1]]/b[[dimL, All, 2]]], 
    lower[[dimL]] = -100];
   If[Length[a[[dimL]]] > 0, 
    upper[[dimL]] = Min[(1 - a[[dimL, All, 1]])/a[[dimL, All, 2]]], 
    upper[[dimL]] = 1];, {dimL, dims}];
  mindist = Norm[upper - lower];
  sigs = 0.5*(upper + lower);
  Do[num = 1 - c[[j, 1]] - Total[c[[j, 2 ;;]]*lower];
   cnorm = Norm[c[[j, 2 ;;]]];
   newdist = Abs[num]/(cnorm + $MachineEpsilon);
   If[newdist < mindist, mindist = newdist;
    sigs = lower + 0.5*num*c[[j, 2 ;;]]/(cnorm*cnorm)];, {j, 
    n[[1]]}];
  s += sigs;
  s1 = Transpose[Prepend[Transpose[s], ConstantArray[1, npoints]]];
  prodGamNum = gammaProdNum[s1, c, n];
  prodGamDenom = gammaProdDenom[s1, c, p, q];
  Do[Do[prodGamNum *= 
      Gamma[1 - a[[dimL, j, 1]] - 
        a[[dimL, j, 2]]*s[[All, dimL]]];, {j, n[[dimL + 1]]}];
   Do[prodGamNum *= 
      Gamma[b[[dimL, j, 1]] + b[[dimL, j, 2]]*s[[All, dimL]]];, {j, 
     m[[dimL + 1]]}];
   Do[prodGamDenom *= 
      Gamma[a[[dimL, j, 1]] + a[[dimL, j, 2]]*s[[All, dimL]]];, {j, 
     n[[dimL + 1]] + 1, p[[dimL + 1]]}];
   Do[prodGamDenom *= 
      Gamma[1 - b[[dimL, j, 1]] - 
        b[[dimL, j, 2]]*s[[All, dimL]]];, {j, m[[dimL + 1]] + 1, 
     q[[dimL + 1]]}];, {dimL, dims}];
  zs = z^-s;
  result = (prodGamNum/prodGamDenom)*Apply[Times, zs, {1}]/(2*Pi)^dims;
  result]

(*Compute multivariate Fox-H function*)
compMultiFoxH[params_, nsubdivisions_, boundaryTol_ : 0.0001] := 
 Module[{boundaries, dim, signs, code, quad, points, volume, result}, 
  boundaries = detBoundaries[params, boundaryTol];
  dim = Length[boundaries];
  signs = Tuples[{1, -1}, dim];
  code = Tuples[Range[0, Floor[nsubdivisions/2] - 1], dim];
  quad = 0;
  Do[points = 
    DiagonalMatrix[sign]*((code + 0.5)*boundaries*2/nsubdivisions);
   quad += Total[Re[compMultiFoxHIntegrand[points, params]]];, {sign, 
    signs}];
  volume = Apply[Times, 2*boundaries/nsubdivisions];
  result = quad*volume;
  result]

(*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}}}};
Print[compMultiFoxH[params1, 20]]

There are several typos in your code. Have you translated your code from Python and Matlab to Mathematica? — Alex Trounev, Commented Jul 21 at 7:27
What is the definition of Multivariate FoxH? A example paper ? — Mariusz Iwaniuk, Commented Jul 21 at 11:36
@MariuszIwaniuk The definition can be found at Eqn(1.1) of AN INTEGRAL INVOLVING GENERALIZED FUNCTION OF TWO VARIABLES and Eqn（A.1) of The H-Function Theory and Applications i.sstatic.net/D0lyUU4E.png — 138 Aspen, Commented Jul 21 at 12:37

138 Aspen · Accepted Answer · 2024-07-21 11:03:03Z

Now it works. The result is consistent with Python/MATLAB.

Any help for accelerating the Mathematica code would be greatly appreciated.

(*Define functions for gamma product computation*)ClearAll["Global`*"]

(*Compute boundaries*)
detBoundaries[params_, tol_] := 
 Module[{boundaryRange, dims, boundaries, points, absIntegrand, 
   index}, boundaryRange = Range[0, 50, 0.05];
  dims = Length[params[[1]]];
  boundaries = ConstantArray[0, dims];
  Do[points = ConstantArray[0, {Length[boundaryRange], dims}];
   points[[All, dimL]] = boundaryRange;
   absIntegrand = Abs[compMultiFoxHIntegrand[points, params]];
   index = 
    Max[Select[Range[Length[absIntegrand]], 
      absIntegrand[[#]] > tol*absIntegrand[[1]] &]];
   boundaries[[dimL]] = boundaryRange[[index]];, {dimL, dims}];
  boundaries]

(*Compute complex integrand of the multivariate Fox-H function*)
compMultiFoxHIntegrand[y_, params_] := 
 Module[{z, mn, pq, c, d, a, b, m, n, p, q, npoints, dims, s, lower, 
   upper, mindist, sigs, num, cnorm, newdist, s1, prodGamNum, 
   prodGamDenom, zs}, {z, mn, pq, c, d, a, b} = params;
  m = mn[[All, 1]];
  n = mn[[All, 2]];
  p = pq[[All, 1]];
  q = pq[[All, 2]];
  {npoints, dims} = Dimensions[y];
  s = I   y;
  lower = ConstantArray[0, dims];
  upper = ConstantArray[0, dims];
  Do[If[b[[dimL]] =!= {}, 
    Module[{bj, Bj}, bj = b[[dimL, All, 1]][[;; m[[dimL + 1]]]];
     Bj = b[[dimL, All, 2]][[;; m[[dimL + 1]]]];
     lower[[dimL]] = -Min[bj/Bj];], lower[[dimL]] = -100];
   If[a[[dimL]] =!= {}, 
    Module[{aj, Aj}, aj = a[[dimL, All, 1]][[;; n[[dimL + 1]]]];
     Aj = a[[dimL, All, 2]][[;; n[[dimL + 1]]]];
     upper[[dimL]] = Min[(1 - aj)/Aj];], upper[[dimL]] = 1];, {dimL, 
    1, dims}];
  mindist = Norm[upper - lower];
  sigs = 0.5   (upper + lower);
  Do[num = 1 - c[[j, 1]] - Total[c[[j, 2 ;;]]   lower];
   cnorm = Norm[c[[j, 2 ;;]]];
   newdist = Abs[num]/(cnorm + $MachineEpsilon);
   If[newdist < mindist, mindist = newdist;
    sigs = lower + 0.5   num   c[[j, 2 ;;]]/(cnorm^2);], {j, 1, 
    n[[1]]}];
  s = Table[s + sigEle, {sigEle, sigs}];
  s = Flatten[s, 1];
  s = ArrayReshape[s, {Length@s/dims, dims}];
  s1 = Table[{1.00}~Join~sEle, {sEle, s}];
  prodGamNum = 1 + 0   I;
  prodGamDenom = 1 + 0   I;
  Do[prodGamNum *= Gamma[(1 - s1 . c[[j]])];, {j, 1, n[[1]]}];
  Do[prodGamDenom *= Gamma[(1 - s1 . d[[j]])];, {j, 1, q[[1]]}];
  Do[prodGamDenom *= Gamma[(s1 . c[[j]])];, {j, n[[1]] + 1, p[[1]]}];
  Do[Do[prodGamNum *= 
      Gamma[(1 - a[[dimL, j, 1]] - 
         a[[dimL, j, 2]]   s[[All, dimL]])];, {j, 1, n[[dimL + 1]]}];
   Do[prodGamNum *= 
      Gamma[(b[[dimL, j, 1]] + 
         b[[dimL, j, 2]]   s[[All, dimL]])];, {j, 1, m[[dimL + 1]]}];
   Do[prodGamDenom *= 
      Gamma[(a[[dimL, j, 1]] + 
         a[[dimL, j, 2]]   s[[All, dimL]])];, {j, n[[dimL + 1]] + 1, 
     p[[dimL + 1]]}];
   Do[prodGamDenom *= 
      Gamma[(1 - b[[dimL, j, 1]] - 
         b[[dimL, j, 2]]   s[[All, dimL]])];, {j, m[[dimL + 1]] + 1, 
     q[[dimL + 1]]}];, {dimL, 1, dims}];
  zs = Table[z^-sEle, {sEle, s}];
  result = (prodGamNum/prodGamDenom)  *
     Table[Times @@ zsEle, {zsEle, zs}]/(2   Pi)^dims // N;
  result
  ]


(*Compute multivariate Fox-H function*)
compMultiFoxH[params_, nsubdivisions_, boundaryTol_ : 0.0001] := 
 Module[{boundaries, dim, signs, code, quad, points, volume, result}, 
  boundaries = detBoundaries[params, boundaryTol];
  dim = Length[boundaries];
  signs = Tuples[{1, -1}, dim];
  code = Tuples[Range[0, Floor[nsubdivisions/2] - 1], dim];
  quad = 0;
  Do[points = 
    Table[sign*((codeEle + 0.5)*boundaries*2/nsubdivisions), {codeEle,
       code}];
   quad += Total[compMultiFoxHIntegrand[points, params]];, {sign, 
    signs}];
  volume = Apply[Times, 2*boundaries/nsubdivisions];
  result = quad*volume;
  result]

(*Example usage*)
params1 = {
   {16.2982237081499, 16.2982237081499, 16.2982237081499, 16.2982237081499}, (*z*)
   {{0, 0}, {2, 1}, {2, 1}, {2, 1},{2, 1}}, (*mn*)
   {{0, 1}, {1, 2}, {1, 2}, {1, 2}, {1, 2}}, (*pq*)
   {}, (*c*)
   {{0, 1, 1, 1, 1}}, (*d*)
   {{{1, 2}}, {{1, 2}}, {{1, 2}}, {{1, 2}}},  (*a*)
   {{{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}}} (*b*)
   };
compMultiFoxH[params1, 20] // NumberForm[#, 20] & // Print

$\begingroup$ This is good version (+1). $\endgroup$
– Alex Trounev
Commented Jul 21 at 14:18 — Alex Trounev, Commented Jul 21 at 14:18

Stack Exchange Network

How to port Matlab/Python's multivariate FoxH implementation in Mathematica?

Ref code

My implementation

Attempt Version 1 (❌)

Attempt Version 2 (❌)

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
numerics
special-functions
numerical-value
or ask your own question.

Hot Network Questions

How to port Matlab/Python's multivariate FoxH implementation in Mathematica?

Ref code

My implementation

Attempt Version 1 (❌)

Attempt Version 2 (❌)

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged numericsspecial-functionsnumerical-value or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
numerics
special-functions
numerical-value
or ask your own question.