# using interpolation functions in ParallelMap

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}];

\[Phi]1 = RandomReal[{0., 2. \[Pi]}];
\[Theta]1 = RandomReal[{0., \[Pi]}];
y = RandomReal[{0., rc/l}];
d = l y;
dx = d Cos[\[Phi]1] Sin[\[Theta]1];
dy = d Sin[\[Phi]1] Sin[\[Theta]1];
dz = d Cos[\[Theta]1];

{y, \[Phi]1, \[Theta]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}];

\[Phi]1 = RandomReal[{0., 2. \[Pi]}];
\[Theta]1 = RandomReal[{0., \[Pi]}];
y = RandomReal[{0., rc/l}];
d = l y;
dx = d Cos[\[Phi]1] Sin[\[Theta]1];
dy = d Sin[\[Phi]1] Sin[\[Theta]1];
dz = d Cos[\[Theta]1];

{y, \[Phi]1, \[Theta]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}];

\[Phi]1 = RandomReal[{0., 2. \[Pi]}];
\[Theta]1 =  RandomReal[{0., \[Pi]}];
y = RandomReal[{0., rc/l}];
d = l y;
dx = d Cos[\[Phi]1] Sin[\[Theta]1];
dy = d Sin[\[Phi]1] Sin[\[Theta]1];
dz = d Cos[\[Theta]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?

-
You really need to work on your questions. To get good answers, it is the best, when you reduce your problem to a simple example. What can I do with your code if it does not run? I can't inspect anything, so it's useless unless I know the answer upfront. Therefore: Background information: yes, Long, not executable code: NO, small example reproducing the problem: YES! –  halirutan Oct 28 '12 at 23:52
@halirutan Thank you for your constructive criticism. I have edited the question with your points in mind. –  BeauGeste Oct 29 '12 at 14:32
Did you try the solution to execute all definitions inside ParallelEvaluate? It works here. –  halirutan Oct 30 '12 at 0:37

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}} *)

-
Hmm, would you call this a bug with DistributeDefinitions? –  Ajasja Oct 29 '12 at 8:46
@Ajasja This is long know, yes and I would call this a bug. I assume you have tried to simply distribute the InterpolatingFunction which does not work? –  halirutan Oct 29 '12 at 9:16
Yes, I've tried that and no, it does not work. I'd like to add the bugs tag (this is why I'm asking). Also do you have any idea why this workaround doesn't work? With[{aip = ip}, ParallelEvaluate[ip = aip]]; –  Ajasja Oct 29 '12 at 12:49
This is possibly related to this? –  Szabolcs Mar 21 '14 at 16:29