ParallelTable behaves differently when it is run inside a package

Question

I've been butting my head against some very weird behaviour by ParallelTable, and its interaction with packages, over the past few weeks, and making relatively little headway towards even reproducing the weird behaviour in a stable way. I just managed to crystallize some of this behaviour in a clean format and I would like some help understanding what's going on.

Consider, as an example of a relatively heavy calculation that one might want to parallelize, the following sum:

AbsoluteTiming[
Table[
  Sum[
   BesselJ[0, 10^-9 k]/(n + 1.6^k), {k, 0, 10000}
   ]
  , {n, 0, 12}]
]

(* {5.67253, {etc.}} *)

If I parallelize this, even for something this small, it gets faster:

AbsoluteTiming[
ParallelTable[
  Sum[
   BesselJ[0, 10^-9 k]/(n + 1.5^k), {k, 0, 10000}
   ]
  , {n, 0, 12}]
]

(* {1.89187, {etc.}} *)

(Minor change in the denominator to avoid what looks like caching.) OK, so far so good. Now, suppose I wish to make this calculation into part of a package, which might look like this:

BeginPackage["package`"];
function::usage = "function[x] is a function to calculate stuff";
RunInParallel::usage = "RunInParallel is an option for function which determines whether it runs in parallel or not.";
Begin["Private`"];

Options[function] = {RunInParallel -> False};

function[x_, OptionsPattern[]] := Block[{TableCommand, SumCommand},
  Which[
   OptionValue[RunInParallel] === False,
    TableCommand = Table; SumCommand = Sum;,
   OptionValue[RunInParallel] === True,
   TableCommand = ParallelTable; SumCommand = Sum;,
   True, TableCommand = OptionValue[RunInParallel][[1]]; 
   SumCommand = OptionValue[RunInParallel][[2]];
   ];
  TableCommand[
   SumCommand[
    BesselJ[0, 10^-9 k]/(n + x^k), {k, 0, 50000}
    ]
   , {n, 0, 12}]
  ]

End[];
EndPackage[];

In particular, I have given it the option RunInParallel to decide whether to use a normal Table or a parallelized one. If I run it like this, however, I get much worse timings:

AbsoluteTiming[function[1.1, RunInParallel -> True]]
AbsoluteTiming[function[1.2, RunInParallel -> False]]

(* {31.465, {etc.}} *)
(* {34.5198, {etc.}} *)

Note here that (i) both versions are much slower than their non-packaged cousins, and (ii) all the speedup from the parallelization is gone.

To try and probe this a bit further, I tried to add some functionality to let me extract the calculation and then run it separately. That is, running

function[1.3, RunInParallel -> {Inactive[ParallelTable], Inactive[Sum]}]]

returns the calculation that it would have run, but with the Table and Sum wrapped in Inactive statements:

Inactive[ParallelTable][
 Inactive[Sum][
  BesselJ[0, Private`k/1000000000]/(1.3^Private`k + Private`n)
 , {Private`k, 0, 50000}]
, {Private`n, 0, 12}]

I can then simply pop them open with a corresponding Activate statement. However, when I do this,

AbsoluteTiming[Activate[function[1.9, RunInParallel -> {Inactive[ParallelTable], Inactive[Sum]}]]]
AbsoluteTiming[Activate[function[1.8, RunInParallel -> {Inactive[Table], Inactive[Sum]}]]]

(* {11.7112, {etc.}} *)
(* {35.7969, {etc.}} *)

the timings come out as something else entirely, yet again. I'm a bit baffled about why the calculation is slower through the package than outside it, but mostly it's the parallelization that bothers me: why isn't the in-package parallelization able to work as well as the Activate[Inactive] route? Why was the parallelization lost in the first place? Did I fall into a bug or something?

Any help in understanding this will be welcome.

(All of this run, by the way, on a 4-core 4-thread Intel Core i5-2500 with 4GB RAM, MM v10.4 over Ubuntu 15.10.)

Answer 1

In short, the reason is that SumCommand's (temporary) definition won't automatically get distributed to the parallel kernels because now SumCommand lives in the package`Private` context, not Global`. This means that SumCommand won't get evaluated to Sum on the subkernels. It gets returned as-is to the main kernel, where now SumCommand does have a definition, gets evaluated to Sum, which in turn gets evaluated to the desired result. But all the evaluation happens on the main kernel.

Aside: Note that Begin["Private`"] should be Begin["`Private`"] so that private symbols will go into package`Private` and not into Private`.

Why is the function slow when in a package?

Most parallel functions, such as ParallelTable, automatically distribute the definitions of symbols. (The notable exception is ParallelEvaluate.)

When you put the function into a package, the context of all the auxiliary symbols (such as SumCommand and TableCommand) change. This prevents the distribution mechanism from kicking in. Package symbols, which reside in contexts other than Global`, do not get distributed by default. This is to prevent the distribution of package symbol definitions to subkernels, which would break packages which must also do initialization (such as loading LibraryFunctions) in addition to issuing definitions. Instead packages should be properly loaded on each subkernel, which can be automated with ParallelNeeds.

Unfortunately I do not fully understand the distribution rules for contexts though ... you can read more at DistributedContexts and links within that page.

This theory that the problem is that SumCommand doesn't get distributed can be verified by adding DistributeDefinitions[SumCommand] right after TableCommand = ParallelTable; SumCommand = Sum;. This will make it run fast again (but it is not a good workaround, see below).

How to fix it?

I have never written packages which use parallelization, so I have no experience with this, and no experience with what the major pitfalls are. But I would not use DistributeDefinitions the way I suggested above (which was only for testing the theory that the problem is with distribution).

One problem with the way Block is used here is that Block won't have any effect across kernels. It only works on the main kernel. Thus if we simply insert DistributeDefinitions[SumCommand] inside of the Block body, the definition will get distributed to the subkernels, but it won't get cleared on the subkernels when the Block finishes. Instead it will persist even after function finishes. You can verify this with ParallelEvaluate[package`Private`SumCommand].

Instead I suggest never sending the symbol SumCommand to the subkernel in the first place. Just send Sum instead. One way to achieve this is with a With-definition (instead of a Block-definition), which does a direct replacement of SumCommand within the body of With.

Here's the final code:

BeginPackage["package`"];
function::usage = "function[x] is a function to calculate stuff";
RunInParallel::usage = "RunInParallel is an option for function which determines whether it runs in parallel or not.";
Begin["`Private`"];

Options[function] = {RunInParallel -> False};

function[x_, OptionsPattern[]] := Block[{TableCommand, SumCommand},
  Which[
   OptionValue[RunInParallel] === False,
    TableCommand = Table; SumCommand = Sum;,
   OptionValue[RunInParallel] === True,
   TableCommand = ParallelTable; SumCommand = Sum;,
   True, TableCommand = OptionValue[RunInParallel][[1]]; 
   SumCommand = OptionValue[RunInParallel][[2]];
   ];
  With[{SumCommand=SumCommand},
   TableCommand[
     SumCommand[
       BesselJ[0, 10^-9 k]/(n + x^k), {k, 0, 50000}
      ]
     , {n, 0, 12}]
   ]
  ]

End[];
EndPackage[];

(To avoid the red colouring you might consider using a different name for the With variable.)

This version is robust and runs fast.

Answer 2

As pointed out by @szabolcs part of the problem is associated with the distribution over the kernels. But I think that is not only problem.

We can generalize the foo function to run in parallel and achieve the same performance as the unpackage code.

The code is much slower in the packaged case, because it is different and lot more complex. If you create a function with the same instructions as the unpackaged code, you should get the equivalent timing and the same results.

For instance, consider only the serial case

BeginPackage["package`"];
foo::usage = "foo[x, test] is a function to calculate stuff";
Begin["Private`"];
foo[ x_Real, test_:False ] := If[
  test == True
, ParallelTable[
    Sum[
      BesselJ[0, 10^-9 k]/(n + x^k), {k, 0, 10000}
    ]
  , {n, 0, 12}
  ]
, Table[
    Sum[
      BesselJ[0, 10^-9 k]/(n + x^k), {k, 0, 10000}
    ]
  , {n, 0, 12}
  ]
]

End[];
EndPackage[];

which runs roughly the same time as the unpackaged code (using ClearSystemCache[] to ensure correct timing)

ClearSystemCache[]
res1 = AbsoluteTiming[
  Table[
    Sum[
      BesselJ[0, 10^-9 k]/(n + 1.6^k)
    , {k, 0, 10000}
    ]
  , {n, 0, 12}]
];
ClearSystemCache[]
res3a = AbsoluteTiming[foo[1.6,True]];
ClearSystemCache[]
res3b = AbsoluteTiming[foo[1.6,False]];
ClearSystemCache[]
res3c = AbsoluteTiming[foo[1.6]];

with comparable timing and exactly the same results

{ res1[[ 1 ]], res3a[[ 1 ]], res3b[[ 1 ]], res3c[[ 1 ]] }
res1[[ 2 ]] == res3a[[ 2 ]] == res3b[[ 2 ]] == res3c[[ 2 ]]

{4.5524, 1.92163, 4.38698, 4.38055}

True

Now, if we define an unpackaged function equivalent to your packaged one

Options[bar] = {RunInParallel -> False};
bar[x_, OptionsPattern[]] := Block[{TableCommand, SumCommand},
  Which[
   OptionValue[RunInParallel] === False
 , TableCommand = Table; SumCommand = Sum;
 , OptionValue[RunInParallel] === True
 , TableCommand = ParallelTable; SumCommand = Sum;
 , True
 , TableCommand = OptionValue[RunInParallel][[1]]; 
   SumCommand = OptionValue[RunInParallel][[2]];
 ];
 TableCommand[
   SumCommand[
     BesselJ[0, 10^-9 k]/(n + x^k), {k, 0, 50000}
   ]
 , {n, 0, 12}]
 ]

evaluating it and comparing with you packaged function (res2)

ClearSystemCache[]
res4 = AbsoluteTiming[ bar[1.6, RunInParallel -> False] ]

{ res2[[1]], res4[[1]] }
res2[[2]]==res4[[2]]
{42.4399, 42.0656}
True

All results are equivalent

res1[[ 2 ]] == res2[[ 2 ]] == res3[[ 2 ]] == res4[[ 2 ]]
True