using interpolation functions in ParallelMap

Question

Edit

I am editing the question to give a simple workable example. It is at the bottom under the heading Example. The original question is now headed withe Background.

Background

I am trying to calculate large tables with values taken from a normal distribution and an interpolation function (as well as a spherical coordinate). The interpolation function is b[j] below where j is an index denoting a certain interpolation function. The interpolation functions come from large arrays.

I would like to calculate for different j in parallel since I have 30 different j to work with (j = 1..30). The problem I am seeing is that the parallel calculation is in fact much slower when the interpolation functions are present. The code is (which you will not be able to run since you don't have the bx,by,bz functions):

 l = 0.2;
 xEdge = 1990.;
 yEdge = 1990.;
 zBottom = 36.;
 zTop = 126.;


 dataMod[j_, rc_] := Module[{},

 X = RandomReal[{rc, xEdge - rc}];
 Y = RandomReal[{rc, yEdge - rc}];
 Z = RandomReal[{zBottom + rc, zTop - rc}];

 ϕ1 = RandomReal[{0., 2. π}];
 θ1 = RandomReal[{0., π}];
 y = RandomReal[{0., rc/l}];
 d = l y;
 dx = d Cos[ϕ1] Sin[θ1];
 dy = d Sin[ϕ1] Sin[θ1];
 dz = d Cos[θ1];

 {y, ϕ1, θ1,
 RandomReal[NormalDistribution[0., 1./2.]] + bx[j][X, Y, Z],
 RandomReal[NormalDistribution[0., 1./2.]] + bx[j][X + dx, Y + dy, Z + dz],
 RandomReal[NormalDistribution[0., 1./2.]] + by[j][X, Y, Z],
 RandomReal[NormalDistribution[0., 1./2.]] + by[j][X + dx, Y + dy, Z + dz],
 RandomReal[NormalDistribution[0., 1./2.]] + bz[j][X, Y, Z],
 RandomReal[NormalDistribution[0., 1./2.]] + bz[j][X + dx, Y + dy, Z + dz]}
 ];

 data[j_, rc_, nMax_] := Table[dataMod[j, rc], {nMax}];

I would like to calculate the table for several different values of j. Here I use only 2 though in actuality, I will use 30:

 list = Map[data[#, 1., 1*10^4] &, Range@2]; // AbsoluteTiming

This took 6.6 seconds. On two kernels,

 listParallel = 
 ParallelMap[data[#, 1., 1*10^4] &, Range@2]; // AbsoluteTiming

took 119 seconds on two kernels. Why is this so much slower? Are the (10000) elements of the table being calculated in parallel or are j = 1 and j = 2 being calculated in parallel (which is what I was aiming for). I thought ParallelMap would calculate j=1 and j=2 in parallel (i.e. simultaneously so I would expect the time to be $\sim 3.3$ seconds).

If I get rid of the interpolation functions,

 dataMod[j_, rc_] := Module[{},

 X = RandomReal[{rc, xEdge - rc}];
 Y = RandomReal[{rc, yEdge - rc}];
 Z = RandomReal[{zBottom + rc, zTop - rc}];

 ϕ1 = RandomReal[{0., 2. π}];
 θ1 = RandomReal[{0., π}];
 y = RandomReal[{0., rc/l}];
 d = l y;
 dx = d Cos[ϕ1] Sin[θ1];
 dy = d Sin[ϕ1] Sin[θ1];
 dz = d Cos[θ1];

 {y, ϕ1, θ1,
 RandomReal[NormalDistribution[0., 1./2.]],
 RandomReal[NormalDistribution[0., 1./2.]],
 RandomReal[NormalDistribution[0., 1./2.]],
 RandomReal[NormalDistribution[0., 1./2.]],
 RandomReal[NormalDistribution[0., 1./2.]],
 RandomReal[NormalDistribution[0., 1./2.]]}
 ];

then

 list = Map[data[#, 1., 1*10^4] &, Range@2]; // AbsoluteTiming

takes 0.99 seconds. The parallel version had a modest speed-up

 listParallel = 
 ParallelMap[data[#, 1., 1*10^4] &, Range@2]; // AbsoluteTiming

to 0.92 seconds which is at least not a lot slower as seen with the interpolation functions present. So why is it that the interpolation functions seem to present a problem for ParallelMap? Is there anything else obvious that I can do do speed this up (in reality nMax must be 1-100 million)?

Example

My question concerns using interpolation functions on multiple kernels with ParallelMap. Here is an example that I think illustrates my problem.

Define the interpolation functions:

 l = 0.2;
 xEdge = 4.;
 yEdge = 4.;
 zBottom = 0.;
 zTop = 4.;

 {cx[1], cx[2], cx[3], cx[4]} =  Table[Interpolation[
  Flatten[Table[{x, y, z, j*Sin[x y z]}, {x, 0, xEdge, 0.1}, {y, 0, 
   yEdge, 0.1}, {z, zBottom, zTop, 0.1}], 2]], {j, 1, 4}];

I now evaluate the interpolation functions at specific points governed by the following:

 ffMod[j_, rc_] := Module[{},

  X = RandomReal[{rc, xEdge - rc}];
  Y = RandomReal[{rc, yEdge - rc}];
  Z = RandomReal[{zBottom + rc, zTop - rc}];

  ϕ1 = RandomReal[{0., 2. π}];
  θ1 =  RandomReal[{0., π}];
  y = RandomReal[{0., rc/l}];
  d = l y;
  dx = d Cos[ϕ1] Sin[θ1];
  dy = d Sin[ϕ1] Sin[θ1];
  dz = d Cos[θ1];

  {cx[j][X, Y, Z], cx[j][X + dx, Y + dy, Z + dz]}
  ];

I then want to make a large table of values (20000 here, but actually want to do for 1 million or more):

 ffData[j_, rc_, nMax_] := Table[ffMod[j, rc], {nMax}];

for several different j.

 ffList = Map[ffData[ #, 1., 2*10^4] &, Range@4]; // AbsoluteTiming

takes 8.9 seconds. In parallel on 4 kernels:

 ffListP = ParallelMap[ffData[ #, 1., 2*10^4] &, Range@4]; // AbsoluteTiming

takes 17.6 seconds. How can I speed up the parallel version of this as I would expect that each kernel could operate on the 4 different j's?

Answer 1

I assume your InterpolatingFunction is not available on the parallel subkernels. Simple expample

ip = Interpolation[Table[{x, Sin[x]}, {x, 0, 2 Pi, Pi/100.}]];
dataMod[] := Module[{},
  Sum[Sum[ip[t], {t, 0, 2 Pi, 0.0001}], {10}]
  ]

Evaluating this in on the main kernel

AbsoluteTiming[dataMod[]]

(* {1.032286, -0.0000626701} *)

And now doing this on 4 subkernels which should take a bit longer since we have to transfer the data but this shows

AbsoluteTiming[ParallelEvaluate[dataMod[]]]

(* {4.024674, {-0.0000626701, -0.0000626701, -0.0000626701, -0.0000626701}} *)

That something is wrong here. Indeet, no ip on the subkernels

ParallelEvaluate[With[{ip = ip},
  Hold[ip]]]

(* {Hold[ip], Hold[ip], Hold[ip], Hold[ip]} *)

When you define it on the parallel kernels (after a Quit[] please)

ParallelEvaluate[
 ip = Interpolation[Table[{x, Sin[x]}, {x, 0, 2 Pi, Pi/100.}]];
 dataMod[] := Module[{},
   Sum[Sum[ip[t], {t, 0, 2 Pi, 0.0001}], {10}]
   ]
];
AbsoluteTiming[ParallelEvaluate[dataMod[]]]

(* {1.168965, {-0.0000626701, -0.0000626701, -0.0000626701, -0.0000626701}} *)