The problem
Consider some pre-generated list with coordinates x1,x2,x3,x4 and values of some function func at these points. For fixed x1 belonging to the x1 range from list but generally not equal to any of the grid points, I need to interpolate it somehow and integrate it with some pre-factor function prefactor[x1,x2,x3,x4] over x2, x3, x4.
The problem is that a brute-force integral using Interpolation is very slow, caused by interpolation built in Mathematica. I attempted to improve the integration speed by ~5-10 times (depending on the values of the parameter and x1) by keeping an accuracy of ~10% - by "mapping" the data onto a grid if assuming the given value of x1 and summing it over x2,x3,x4 (the code is given below). However, the implementation is very ugly - I am sure it is not optimized and contains rough approximations, while the speedup may be higher.
Could you please help me speed up the integration in a proper way?
Test data
This is the test list:
func[x1_, x2_, x3_, x4_] = If[x3 > x1, Exp[-(Sqrt[x1^2 + x3^2]/5)]*Sin[x2]*(Cos[20*x2])^x4*((100 - x4)/50)^2, 10^-90];
gridx1 = Table[x, {x, 0.05, 5.1, (5.1 - 0.05)/30}];
gridx2 = Table[10^x, {x, -5, Log10[0.05], (Log10[0.05] + 5)/20}];
gridx3 = Table[10^x, {x, Log10[0.051], Log10[350], (Log10[350] - Log10[0.051])/50}];
gridx4 = Table[x, {x, 38., 88., 2.}];
GridIn1 = {gridx1, gridx2, gridx3, gridx4};
GridIn = Tuples[GridIn1];
list = Join[GridIn, Partition[func @@@ GridIn, 1], 2];
Below, there are the interpolation funcInt, the pre-factor prefactor and the integral intv:
funcInt[x1_,x2_,x3_,x4_] = 10^(Interpolation[Log10[list],InterpolationOrder->1][Log10[x1],Log10[x2],Log10[x3],Log10[x4]]);
prefactor[x1_, x2_, x3_, x4_, ct_] = Exp[-x4/(Cos[x2]*ct*x3/x1)]/(Cos[x2]*ct*x3/x1);
intv[x1_, ct_] := NIntegrate[funcInt[x1, x2, x3, x4]*prefactor[x1, x2, x3, x4, ct], {x2, Min[gridx2], Max[gridx2]}, {x3, Min[gridx3], Max[gridx3]}, {x4, Min[gridx4], Max[gridx4]}, Method -> "AdaptiveMonteCarlo"]
Here, ct is a positive parameter with meaningful values between 0.01 and 10^8. Like everything that I do in Mathematica, it is very slow:
intv[4, 30] // AbsoluteTiming
{0.708241, 0.0000166047}
I interpolated over the logarithmized data to accurately integrate into the region where the integrand drops exponentially (the case of small ct).
The potential speedup may be (compared to the integration of the initial function func) 10 times:
intv0[x1_, ct_] :=
NIntegrate[
func[x1, x2, x3, x4]*prefactor[x1, x2, x3, x4, ct], {x2,
Min[gridx2], Max[gridx2]}, {x3, Min[gridx3], Max[gridx3]}, {x4,
Min[gridx4], Max[gridx4]}, Method -> "AdaptiveMonteCarlo"]
intv0[4, 30] // AbsoluteTiming
{0.0670157, 0.0000184466}
It may also be imlroved further, depending of the performance of the built-in integration methods in Mathematica.
My attempt
I made the following attempt. For the given value of x1, I "map" Log10[list] onto a grid {{x1},gridx21,gridx31, gridx41}, where gridx21,gridx31,gridx41 are denser than the initial grid gridx1,.... Then, having the grid I compute the size of the intervals Delta x, multiply the mapped list with the intervals and the prefactor, and then sum over it.
It gives me a $\sim 5$x speedup per list of ct, with ~10% accuracy compared to intv.
However, I am pretty sure that the implementation is very dumb. First, it is not defined correctly mathematically (there are some rough approximations) - instead of using, e.g., Darboux sum $\sum_{i}\text{sup}(f_{x\in \{x_{i}\}})\Delta x_{i}$ I just sum $\sum_{i}f_{x_{i}}\Delta x_{i}$, where $f(x_{i})$ is just a value of the function at a point $x_{i}$, with no adequate relation between $x_{i}$ and $\Delta x_{i}$. Second, in particular, it becomes very inconvenient if there are several lists (with different grids) that I need to interpolate and integrate. I would appreciate it if you comment on this implementation.
So, first, let us define the "mapping" function - one may use any of the beautiful solutions presented in this post:
nd = 4;
cf4 = Module[{xgvars = Unique["xg"] & /@ slist @@ Range@nd,
igvars = Unique["ig"] & /@ slist @@ Range@nd,
tgvars = Unique["tg"] & /@ slist @@ Range@nd,
ivars = Unique["i"] & /@ slist @@ Range@nd,
svars = Unique["s"] & /@ slist @@ Range@nd,
tvars = Unique["t"] & /@ slist @@ Range@nd,
jvars = Unique["j"] & /@ slist @@ Range@nd},
Inactivate[
Compile[{seq@{xgvars, _Real, 1}, {y, _Real, nd}},
Module[{seq@igvars, seq@tgvars, seq@ivars, seq@svars,
seq@tvars},
seq[igvars =
Floor[xgvars] - UnitStep[xgvars - indexed@Dimensions@y]];
seq[tgvars = xgvars - igvars];
Table[seq[tvars = Compile`GetElement[tgvars, jvars]];
seq[ivars = Compile`GetElement[igvars, jvars]];
seq[svars = 1. - tvars];
eval@Total[
Times @@ #[[All, 1]] Compile`GetElement[y,
Sequence @@ #[[All, 2]]] & /@
Tuples@Transpose[{{svars, ivars}, {tvars,
ivars + 1}}, {2, 3, 1}]],
seq@{jvars, Length@igvars}]], CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"], Except[seq | eval | indexed]] /.
seq[expr_] :>
RuleCondition@(Sequence @@
Table[Inactivate[expr,
Except[slist | indexed]] /. {l_slist :> l[[i]],
indexed[l_] :> Inactive[Compile`GetElement][l, i]}, {i,
nd}]) /. eval@expr_ :>
RuleCondition[Activate[expr /. slist -> List]]] // Activate;
Next, let us define the out grid and the intervals:
(*Slightly different maximal value of x2max since cf4 gets stuck for x2max=0.05*)
gridx21 = Table[x, {x, 10^-5., 0.05, (0.05 - 10^-5.)/40}];
(*Out grid for x3. Will be computed once x1 is fixed*)
gridx31temp[x1_] := Table[x, {x, x1, 350., (350. - x1)/200}];
gridx41 = Table[x, {x, 38., 88., 0.5}];
(*Integvals*)
DxVals[list_] :=
Join[Table[
list[[i]] - list[[i - 1]], {i, 2, Length[list], 1}], {list[[-1]] -
list[[-2]]}];
(*Delta x2 and Delta x4*)
Dx2vals = DxVals[gridx21];
Dx4vals = DxVals[gridx41];
(*Reshaped array of logarithmized values of the function from list. Needs for mapping*)
vals = ArrayReshape[
Log10[list[[All, 5]]], {Length[gridx1], Length[gridx2],
Length[gridx3], Length[gridx4]}];
(*Compiled code which computes the product of Delta x*)
DxValsTotComp =
Hold@Compile[{{Dx2vals, _Real, 1}, {Dx3vals, _Real,
1}, {Dx4vals, _Real, 1}},
Table[{Dx2vals[[i]]*Dx3vals[[j]]*Dx4vals[[k]]}, {i, 1,
Length[Dx2vals]}, {j, 1, Length[Dx3vals]}, {k, 1,
Length[Dx4vals]}], CompilationTarget -> "C",
RuntimeOptions -> "Speed"] // ReleaseHold;
Next, let us define the block returning the mapped list:
(*Block producing the mapped list*)
BlockDiscreteIntegral[x1_] := Block[{},
gridx31 = gridx31temp[x1];
Dx3vals = DxVals[gridx31];
GridOut1 = Log10[{{x1 // N}, gridx21, gridx31, gridx41}];
GridOut = 10^Tuples[{GridOut1[[2]], GridOut1[[3]], GridOut1[[4]]}];
Dxvals =
Flatten[DxValsTotComp[Dx2vals, Dx3vals, Dx4vals], {1, 2, 3}];
(*Mapping routine*)
xig = MapThread[
Interpolation[Transpose@{#, Range@Length@#},
InterpolationOrder -> 1][#2] &, {Log10[GridIn1], GridOut1}];
Join[GridOut, 10^Partition[Flatten@cf4[Sequence @@ xig, vals], 1],
Dxvals, 2]
]
Finally, this is the block returning the integral for the list ctlist:
(*Compiled summed integral*)
IntegralCompiled =
Hold@Compile[{{tab, _Real, 2}, {ct, _Real}, {x1, _Real}},
Total[((*prefactor[x1,#[[1]],#[[2]],#[[3]],ct]**)
Exp[-#[[3]]/(Cos[#[[1]]]*ct*#[[2]]/x1)]/(
Cos[#[[1]]]*ct*#[[2]]/x1) #[[4]]*#[[5]]) & /@ tab],
CompilationTarget -> "C", RuntimeOptions -> "Speed"] // ReleaseHold
BlockIntegralList[x1_, ctlist_] := Block[{},
dat = BlockDiscreteIntegral[x1];
Table[{ctlist[[i]], IntegralCompiled[dat, ctlist[[i]], x1]}, {i, 1,
Length[ctlist], 1}]
]
Let us compare the integral and the discrete integral:
mtest = 4;
ctlist = {0.5, 1, 2, 3, 4, 5, 6, 10, 50, 100, 1000, 10^4} // N;
t1 = BlockIntegralList[mtest, ctlist]; // AbsoluteTiming // First
t2 = Table[{ctlist[[i]], intv[mtest, ctlist[[i]]]}, {i, 1,
Length[ctlist], 1}]; // AbsoluteTiming // First
t2[[All, 2]]^-1 (t1[[All, 2]] - t2[[All, 2]])
1.79184
10.8851
{0.0916109, 0.0606622, 0.0491823, 0.0590111, 0.0686524, 0.113182,
0.0597432, 0.0117492, -0.0875845, -0.0977589, -0.032049, -0.0811554}
Another approach would be to generate values of the interpolated function for random points inside the integration domain and sum them. In principle, I may use cf4 for this by calling it for any particular point, but the problem is that it is not optimized for this purpose.
