Hello Wolfram Community,
I am interested in the parallel evaluation integrals where the integrand is very expensive to compute. Such a feature can be found in some implementations of “cubature” where the whole vector of points needed for the current recursive subdivision are made available so that they can be handled in whatever batch seems appropriate.
I have not found such a capability in Wolfram Language/Mathematica, but I have attached a notebook where I demonstrate it is possible “trick” NIntegrate into yielding its points by supplying a fake/dummy function that is rapid to evaluate, such as f[x]:=0. Those are Reap/Sowed, and then sent to ParallelMap. (Note that previous recursion must be successively provided to this trick so that it knows how to subdivide the current level!)
After memoizing those values, they can be supplied to a second invocation of NIntegrate, but this time using the memoized variable as the integrand.
(Here, I use “LocalAdaptive” because I believe that it optimally-reuses lower levels of recursion.)
Disclaimer: this is a rough cut of an idea, and not well tested. But I am finding it useful in higher dimensions, and the technique is applicable to other functions which can be expensive an adaptive, like ContourPlot3D as well.
I am interested in the parallel evaluation integrals where the integrand is very expensive to compute. Such a feature can be found in some implementations of “cubature” where the whole vector of points needed for the current recursive subdivision are made available so that they can be handled in whatever batch seems appropriate.
I have not found such a capability in Wolfram Language/Mathematica, but I have attached a notebook where I demonstrate it is possible “trick” NIntegrate into yielding its points by supplying a fake/dummy function that is rapid to evaluate, such as f[x]:=0. Those are Reap/Sowed, and then sent to ParallelMap. (Note that previous recursion must be successively provided to this trick so that it knows how to subdivide the current level!)
After memoizing those values, they can be supplied to a second invocation of NIntegrate, but this time using the memoized variable as the integrand.
(Here, I use “LocalAdaptive” because I believe that it optimally-reuses lower levels of recursion.)
Disclaimer: this is a rough cut of an idea, and not well tested. But I am finding it useful in higher dimensions, and the technique is applicable to other functions which can be expensive an adaptive, like ContourPlot3D as well.
In[]:=
LaunchKernels[]
Out[]=
KernelObject,KernelObject,KernelObject,KernelObject
In[]:=
f[x_?NumericQ]:=Cos[Log[x]/x];
1
x
In[]:=
Plot[f[x],{x,0.1,0.2}]
Out[]=
NIntegrate Reference Values
NIntegrate Reference Values
Establish reference values for integration using standard way:
MaxRecursion 0
MaxRecursion 0
In[]:=
NIntegrate[f[x],{x,0.1,0.2},Method->{"LocalAdaptive","SymbolicProcessing"->0},MinRecursion->0,MaxRecursion0]
Out[]=
-0.0955463
MaxRecursion 1
MaxRecursion 1
In[]:=
NIntegrate[f[x],{x,0.1,0.2},Method->{"LocalAdaptive","SymbolicProcessing"->0},MinRecursion->1,MaxRecursion1]
Out[]=
-0.100505
MaxRecursion 2
MaxRecursion 2
In[]:=
NIntegrate[f[x],{x,0.1,0.2},Method->{"LocalAdaptive","SymbolicProcessing"->0},MinRecursion->2,MaxRecursion2]
Out[]=
-0.100476
MaxRecursion 3
MaxRecursion 3
In[]:=
NIntegrate[f[x],{x,0.1,0.2},Method->{"LocalAdaptive","SymbolicProcessing"->0},MinRecursion->3,MaxRecursion3]
-6
10
-7
10
Out[]=
-0.100476
ParallelNIntegrate Prototype
ParallelNIntegrate Prototype
In[]:=
ParallelNIntegrate[f_,{xmin_,xmax_},maxr_]:=Module[{memoize,mtemp,res,fake,mem,xnew,known,minr,integrand},(*keepalistofknownxvalues*)known={};(*ensurenumericevaluationinsideNIntegrate*)mem[x_?NumericQ]:=memoize[x];fake[x_?NumericQ]:=0;(*incrementMaxRecursionfrom0tomaxr*)Do[(*FirstInvocation:"trick"NIntegratesothatiteitherlooksupvaluesfrompreviouslowerrecursivelevels,orquicklyyieldsafakevalue*)(*Print["Recursion r = "<>ToString[r]];Print["---------"];*)xnew=Quiet[Reap[NIntegrate[If[MemberQ[known,x],mem[x],fake[x]],{x,xmin,xmax},Method->{"LocalAdaptive","SymbolicProcessing"->0},MinRecursion->r,MaxRecursionr,EvaluationMonitor:>((*Print[{"1st: ",x,If[MemberQ[known,x],mem[x],fake[x]],MemberQ[known,x]}];*)Sow[x])]][[2,1]]];(*returnthepreviousSecondInvocationresultifnonewdata*)If[Length[xnew]==0,Return[res]];(*forvaluesthatwehavenotseen,evaluateinparallel*)xnew=Complement[xnew,known];mtemp=ParallelMap[f[#]&,xnew,Method->"FinestGrained"];(*memoizethosevaluesthatwerecomputedinparallel*)Do[memoize[xnew[[i]]]=mtemp[[i]],{i,Length[mtemp]}];(*appendthosex-valuestotheknownlist*)known=Join[known,xnew];(*SecondInvocation:NIntegratewithmemoizedvalues*)res=Quiet[NIntegrate[mem[x],{x,xmin,xmax},Method->{"LocalAdaptive","SymbolicProcessing"->0},MinRecursion->r,MaxRecursionr(*,EvaluationMonitor:>Print[{"2nd: ",x,mem[x]}]*)]];(*incrementuptomaxr*),{r,0,maxr}];res];
Sleep Test
Sleep Test
Insert a Pause statement to simulate a long running integration
Best regards,
Eric
Eric