Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this community
I need to analyze very large datasets, and they obviously don't fit in a computer memory. I need to produce some statistics on these datasets, and this brings me here:
How can I make some statistical analysis on distributed datasets?
I imagine Mathematica may split (hopefully automagically) data over many kernels, being local or remote.
Anyone had any experience?
This is not exactly an answer to the question as posed, but rather a suggestion for how one can potentially cause the separate subkernels to be able to communicate with one another. I have previously posted this on MathGroup, but there seems not to be any good reason not to present it here as well.
The approach demonstrated is based on the idea that MathLink connections can be established between subkernels on a peer-to-peer basis. Then, one has a message passing mechanism that could in principle be employed similarly to MPI. However, the code shown here does not include most essential elements of MPI, e.g. collective communication primitives, which would need to be implemented first.
Firstly, we start the subkernels between which MathLink connections will later be established. Since the environment will no longer be stateless after the subkernels are able to communicate with one another, we must set the appropriate options so that there will be no attempt at recovery in case any subkernel fails.
SetSystemOptions["ParallelOptions" -> {
"RecoveryMode" -> "Abandon", "RelaunchFailedKernels" -> False
}];
Quiet[CloseKernels[]; LaunchKernels[4], LaunchKernels::nodef];
Internally, the Parallel` package is based on anonymous kernel objects signifying locations to which expressions may be submitted in a queue for evaluation. The $KernelIDs correspond to these abstract objects and may be duplicated or reassigned arbitrarily; they are consequently unsuitable to our present purpose, necessitating the use of $RankID instead.
ranks = Range@Length[
kernels = Kernels[]
];
With[{rankAssignments = Unevaluated[$RankID = #] & /@ ranks},
Parallel`Developer`ParallelDispatch[rankAssignments, kernels]
];
Next, we define the Rank object, which will store (as downvalues) the properties of each of the ranks and implement (through upvalues) the context-sensitive substitution of these properties into various expressions for which they are required. Note that KernelObjects and Parallel` evaluation functions are not ordinarily propagated to the subkernels since they refer to and operate on the subkernels' $ParentLink. Therefore, we define the following relations only in the master kernel.
ClearAll[Rank];
SetAttributes[Rank, NHoldFirst];
Rank[ranks_List?(VectorQ[#, IntegerQ] &), "KernelObject"] :=
Rank[#, "KernelObject"] & /@ ranks;
MapThread[Set, {Rank[ranks, "KernelObject"], kernels}];
Rank[_, "KernelObject"] = $Failed;
Rank /: ParallelEvaluate[cmd_, Rank[ranks : _Integer | _List?(VectorQ[#, IntegerQ] &)]] :=
ParallelEvaluate[cmd, Rank[ranks, "KernelObject"]];
Rank /: Parallel`Developer`ParallelDispatch[
cmds_, Rank[ranks : _Integer | _List?(VectorQ[#, IntegerQ] &)]
] /; Length[cmds] == Length[ranks] :=
Parallel`Developer`ParallelDispatch[cmds, Rank[ranks, "KernelObject"]];
Having created a mapping between subkernels and ranks, we are now in a position to open MathLink connections between them. For each pair of ranks we require that one endpoint creates a link, to which its partner subsequently connects. A convenient way to divide this work (and that does not depend on the number of ranks) is to partition the matrix of connections into upper and lower triangular portions. Here we form the lower triangle and hence determine by which ranks the links are to be created.
With[{
allRanks = $AllRanks = ranks,
rankCount = $RankCount = Length[ranks],
machineRules = ParallelEvaluate[$RankID -> $MachineID]
},
ParallelEvaluate[
$AllRanks = allRanks;
$RankCount = rankCount;
locallyLinkedRanks = Take[$AllRanks, $RankID - 1];
localLinks = If[#1 === #2,(* Same $MachineID on both link endpoints? *)
(* Optimization for intra-machine links. *)
LinkCreate[LinkProtocol -> "SharedMemory"],
(* Inter-machine link; must use TCP/IP. *)
LinkCreate[LinkProtocol -> "TCPIP"]
] & @@@ Thread[{$RankID, locallyLinkedRanks} /. machineRules];
(* We also create a loopback link. Although this is technically unnecessary,
it may prove useful for testing or help to simplify the construction of
various collective communication algorithms. Loopback links do not have
any LinkProtocol. *)
ownLink = LinkCreate[LinkMode -> Loopback];
]
];
With the necessary links created, all that remains is for the remaining rank in each pair to connect to the link hosted by its partner. After the links have been opened from both ends, the partners each call LinkActivate on the connection in order to acknowledge each other.
With[{
linkRules = ParallelEvaluate[
(* Results must be returned in rank order or otherwise links will
be established between incorrectly permuted pairs of ranks. *)
Thread[locallyLinkedRanks -> localLinks] /.
linkObject : LinkObject[linkName_String, __] :> {
linkName,
LinkProtocol -> "LinkProtocol" /. MathLink`LinkDeviceInformation[linkObject]
},
Rank[$AllRanks]
]
},
ParallelEvaluate[
remoteLinks = LinkConnect @@@ Flatten[ReplaceList[$RankID, linkRules], 1];
Scan[
(* LinkActivate is blocking at both endpoints and so the calls
should be made in a particular order so as to avoid deadlock.
The rotation is not required for correctness but rather overlaps
pairwise calls that have no serial dependence. *)
LinkActivate,
RotateRight[$RankLinks = Flatten[{localLinks, ownLink, remoteLinks}], $RankID - 1]
];
]
];
The links having been established, it is now useful to be able to call MathLink functions with a rank ID rather than a LinkObject, so here we overload Rank in order to accomplish this. This time, the definitions are are rank-specific, and meaningful only to the subkernels.
ParallelEvaluate[
ClearAll[Rank];
SetAttributes[Rank, NHoldFirst];
Rank[ranks_List?(VectorQ[#, IntegerQ] &), "LinkObject"] :=
Rank[#, "LinkObject"] & /@ ranks;
MapThread[Set, {Rank[$AllRanks, "LinkObject"], $RankLinks}];
Rank[_, "LinkObject"] = $Failed;
Rank /: LinkWrite[Rank[rank_Integer], expr_] :=
LinkWrite[Rank[rank, "LinkObject"], Unevaluated[expr]];
Rank /: LinkWrite[Rank[ranks_List?(VectorQ[#, IntegerQ] &)], expr_] :=
LinkWrite[#, Unevaluated[expr]] & /@ Rank[ranks, "LinkObject"];
Rank /: LinkFlush@Rank[rank_Integer] :=
LinkFlush@Rank[rank, "LinkObject"];
Rank /: LinkFlush@Rank[ranks_List?(VectorQ[#, IntegerQ] &)] :=
LinkFlush /@ Rank[ranks, "LinkObject"];
Rank /: LinkReadyQ@Rank[rank_Integer] :=
LinkReadyQ@Rank[rank, "LinkObject"];
Rank /: LinkReadyQ@Rank[ranks_List?(VectorQ[#, IntegerQ] &)] :=
LinkReadyQ /@ Rank[ranks, "LinkObject"];
Rank /: LinkRead[Rank[rank_Integer], h_: Identity] :=
LinkRead[Rank[rank, "LinkObject"], h];
Rank /: LinkRead[Rank[ranks_List?(VectorQ[#, IntegerQ] &)], h_: Identity] :=
LinkRead[#, h] & /@ Rank[ranks, "LinkObject"];
];
To demonstrate the general principle, let's perform a simple all-to-all collective communication:
ParallelEvaluate[
LinkWrite[Rank[#],
"Message from rank " ~~ ToString[$RankID] ~~ " to rank " ~~ ToString[#]
] & /@ $AllRanks;
LinkReadyQ@Rank[$AllRanks],
Rank[$AllRanks]
]
(* -> {{True, False, False, False}, {True, True, False, False},
{True, True, True, False}, {True, True, True, True}} *)
ParallelEvaluate[
LinkFlush@Rank[$AllRanks];
LinkReadyQ@Rank[$AllRanks],
Rank[$AllRanks]
]
(* -> {{True, True, True, True}, {True, True, True, True},
{True, True, True, True}, {True, True, True, True}} *)
ParallelEvaluate[
LinkRead@Rank[$AllRanks],
Rank[$AllRanks]
]
(* -> {{"Message from rank 1 to rank 1", "Message from rank 2 to rank 1",
"Message from rank 3 to rank 1", "Message from rank 4 to rank 1"},
{"Message from rank 1 to rank 2", "Message from rank 2 to rank 2",
"Message from rank 3 to rank 2", "Message from rank 4 to rank 2"},
{"Message from rank 1 to rank 3", "Message from rank 2 to rank 3",
"Message from rank 3 to rank 3", "Message from rank 4 to rank 3"},
{"Message from rank 1 to rank 4", "Message from rank 2 to rank 4",
"Message from rank 3 to rank 4", "Message from rank 4 to rank 4"}} *)
So, that's very nice, but what can we do productively? As an example, let's try a slightly more complicated communication pattern: the dissemination barrier, which can be used for synchronization between ranks. Each rank waits at the barrier for as long as necessary for all ranks to catch up; they then exit the barrier simultaneously.
ParallelEvaluate[
BarrierMessage;
(* This implementation of the dissemination barrier algorithm will
deadlock if any of the ranks are duplicated because a blocking read
is used for synchronization. The ordering of the participating ranks
must be consistent between calls from different ranks for the same reason. *)
Barrier@Rank[ranks_List?(VectorQ[#, IntegerQ] && MemberQ[#, $RankID] &)] :=
MapThread[(
LinkWrite[#1, BarrierMessage[ranks]]; LinkFlush[#1];
Block[{msg},
(* Here we simply discard irrelevant messages.
A receive queue would be a considerable improvement over this. *)
While[msg =!= BarrierMessage[ranks], msg = LinkRead[#2]]
]
) &,
With[{
rankCount = Length[ranks],
thisRankPosition = Position[ranks, $RankID][[1, 1]]
},
Map[
Rank,
ranks[[#]] & /@ Mod[
0~Range~Log[2, rankCount] //
{thisRankPosition + 2^#, thisRankPosition - 2^#} &,
rankCount, 1
],
{2}
]
]
];
Barrier@Rank[_] = Null,
Rank[$AllRanks]
];
And, the result:
ParallelEvaluate[
First@AbsoluteTiming[
Through@{Print, Pause}@RandomReal[10];
Barrier@Rank[$AllRanks]
],
Rank[$AllRanks]
] // AbsoluteTiming
(* -> (prints) 6.43307
(prints) 6.14898
(prints) 7.92761
(prints) 1.00838
{7.9375000, {7.9375000, 7.9375000, 7.9375000, 7.9375000}} *)
Finally, to clean up, we close all of the links and clear the symbols defined above.
ParallelEvaluate[
ClearAll[BarrierMessage, Barrier];
ClearAll[Rank];
LinkClose /@ $RankLinks;
$RankLinks =.;
$RankCount =.;
$AllRanks =.;
$RankID =.;
];
ClearAll[Rank];
$RankCount =.;
$AllRanks =.;
More an extended comment rather than a complete answer (you've posed a rather general question).
I like this question.
Mathematica certainly includes the components and operations to do this. Unless some feature exists that hasn't come to mind, you'll need to combine those components and operations as a bespoke solution.
That said...
Kernels and RAM - You need more than just kernels. You need kernels on different machines (physical or virtual) each with dedicated RAM.
Assign data responsibility - Consider linking a machine ID ($MachineID) to a specific chunk of data readable from some central file or database. You could put all of these into a function that you distribute (DistributeDefinitions) to the remote kernels. To make this more easily scalable, you could create a table of machine ID's paired with ranges of data (positions in the file or specific records in a database). Place the table in a .csv file or some equivalent in a central location and have your distributed function read it into each kernel. This would allow you to add or delete kernels and reassign their responsibility without needing to rewrite your code.
None of this unusual or particularly insightful, just thinking through one way to accomplish what you want to do.
Read the data - Look at:
In the documentation to read in discrete chunks of the data to each remote kernel.
Collect the results - Run your analytics then collect the results in the main kernel.
Seems very workable to me and maybe not as daunting a task as @Oleksandr R. suggests.
Sorry no Mathematica automagic comes to mind ;-)
I just don't know enough about your specific problem to suggest more at this point. Any additional detail you can provide would help.
Still, a very interesting question.
