Parallelization with xAct
Edu Serna
to run on a single core. I am experimenting with parallelizing it
using the nooby-friendly parallel commands in Mathematica (such as
ParallelMap) which are supposed to distribute variables among the
kernels and such stuff. However I am getting the following warning
message as if somehow the parallel computation didn't know about some
of the xAct definitions:
RuleDelayed::rhs: "Pattern xAct`xTensor`Private`a_ appears on the
right-hand side of rule
HoldPattern[xAct`xTensor`Private`makepattern[xAct`xTensor`Private`a_Symbol]]:>xAct`xTensor`Private`a_.
"
Nonetheless I don't seem to get any problems.
Thanks in advance
JMM
This is just a warning message. There is no actual error.
In general you do not want to have definitions with patterns on the
right hand side, and hence Mathematica sends this warning message to
report the possibility of an error. But sometimes you do need to have
patterns on the right hand side, for example if you are constructing
patterns programmatically, something frequent in xAct. In fact there
are a few Off[ RuleDelayed::rhs ] ... On[ RuleDelayed::rhs ] in
xTensor.m to avoid having this warning message at loading time.
I don't know why parallelism is making these warnings reappear. Can
you turn them off by writing your own Off[ RuleDelayed::rhs ] ?
Automatic parallelization is another future project for xAct. Tell us
about your experiences, and whether you saw a big difference in
parallelizing anything in particular. I guess the most important thing
would be parallelizing ToCanonical on sums. There are two definitions
to touch: MapIfPlus[f_,expr_Plus] (in xCore.m) and
ToCanonicalObject[expr_Plus,options___] in xTensor.m.
Cheers,
Jose.
Edu Serna
Thanks so much for your reply, I will describe my current setup in
further detail. First thing I load
Needs["xAct`xTensor`"]
ParallelNeeds["xAct`xTensor`"]
Off[RuleDelayed::rhs]
Off[RuleDelayed::rhs] has taken care of the warning message as when I
apply mathematica core functions to xact objects, I am not sure if
this is what you meant, or I have to change the definition in
xTensor.m.
ParallelNeeds["xAct`xTensor`"] is needed in order to have the
replacement rules I have defined with MakeRule work in parallel.
With this setup however xact's functions are not recognised by the
other kernels, in particular I am trying ContractMetric, in order to
get this to parallelize I have tried (with success)
DistributeDefinitions[ContractMetric]. This produce the same warning
message again from each Kernel:
RuleDelayed::rhs: Pattern xAct`xTensor`Private`a_ appears on the right-
hand side of rule
At the moment I am trying to set up this to work on a dual Core laptop
before trying the actual run in a 24 core remote server.
Thanks everyone for their time.
Thomas Bäckdahl
I just wonder which part of xTensor that takes so long time for you that
you have to parallelize.
Is it MakeRule, or is it many things?
I myself have problems with MakeRule took long time when I had
expressions with many indices with a lot of contractions and symmetries.
This was because very many rules had to be generated. I wrote a new
MakeCompareRule that generates only one rule that cares about the index
configuration only at "apply time". With this I managed to cut my
computational time for some things 100 times or so. It is still not very
well written, and my current implementation can only handle tensors or
product of tensors in the left hand side. I will have to write something
more to handle derivatives, but this can be done.
I am happy to post the code if you think this could solve some of your
problems.
Still it would be very nice if everything can be parallelized anyway.
Regards
Thomas
Edu Serna
I doubt my MakeRules are there the problem as I mostly replace the
contraction of two vectors with scalar quantities (I imagine i could
write my own cheaper rules but I doubt its worth the effort) as I
calculate scattering amplitudes. The thing is just that my expressions
are VERY long (hundreds of thousands of terms) so I need to do it in
parallel.
I would like to see that code just to see how you did it but I doubt
that is my problem
Thanks so much anyways.
JMM
computation as much as you can. A few recommendations:
As Thomas suggests, the most important thing is to find out which part
of the computation takes more time, to focus on improving that part.
The main idea is to use ToCanonical only when it comes to reorder
indices of an individual product of tensors, and do the rest by hand.
Let us start with an arbitrary tensor expression texpr:
1) First use Expand[ texpr ] if there are nonexpanded products of sums
of tensors.
2) Use SameDummies on the result of 1). This will try to simplify the
expression by minimizing the number of different dummies, without
actually calling group-theoretical algorithms.
3) Map ToCanonical to each term, or use a Do loop to monitor timings
of individual canonicalizations. If the result is much smaller than
the input, it might be useful to add intermediate results to an
accumulator variable, instead of doing a big sum at the end.
Each call to ToCanonical takes a couple milliseconds in the MathLink
transfers, and then a few more milliseconds in the actual
computations, unless it happens to be a case of very high symmetry in
the contractions, which could potentially take very long (in the order
of seconds or worse). If you identify in advance the structure of the
terms with very high symmetry, perhaps you can canonicalize them
separately, with an adapted MyToCanonical.
Jose.
Thomas Bäckdahl
Here I post the code for my MakeCompareRule, even though this thread
might not be the right place for this.
It is unfortunately not well written and is a bit complicated.
It is based on EqualExpressionsQ to determine if two tensor expressions
are equal if one would change the indices.
SpecialEqualExpressionsQ does basically the same thing as
EqualExpressionsQ, but it also gives information about how expr1 and
expr2 are related. (sign changes, index changes etc.) I think this code
can be improved a bit, but I have not looked at it in a while.
SpecialEqualExpressionsQ[expr1_, expr2_] :=
Module[{inds1, frees1, frees2, dummies1, dummies2, sortedfrees2,
dummies1slots, frees1slots, slotrules1, expr1SGS, stabilizerGS,
transversal1, indstoslotnumbers, freesym1, perms1, arrexpr2,
newexpr2,
count, result = False, perm = {}, xxx, equation,
solution = 0},
(* Indices of expr1 *)
inds1 = FindIndices[expr1];
dummies1 = List @@ xAct`xTensor`Private`TakeEPairs@inds1;
frees1 = List @@ xAct`xTensor`Private`TakeFrees@inds1;
(* Indices of expr2 *)
{frees2, dummies2} =
List @@ xAct`xTensor`Private`FindFreeAndDummyIndices[expr2];
If[Length[frees1] == Length[frees2],
sortedfrees2 = List @@ IndexSort[frees2];
dummies2 = List @@ dummies2;
inds1 = List @@ inds1;
{slotrules1, expr1SGS} = List @@ Take[SymmetryOf@expr1, -2];
indstoslotnumbers = (Reverse /@ slotrules1) /.
xAct`xTensor`Private`slot[x_] -> x;
dummies1slots = dummies1 /. indstoslotnumbers;
frees1slots = frees1 /. indstoslotnumbers;
xAct`xTensor`Private`checkChangeFreeIndices[frees1, sortedfrees2];
(* Computing permutations of the free indices of expr1 *)
(* One coud use perms1=Permutations[List@@frees1],
but this can be reduced by considering the symmetries of the \
free indices of expr1 *)
(* The stabilizer of the dummies for the symmetry of expr1 gives \
a subset of the symmetry on dummies1.
The numbers in the group correspond to the position of the index \
in frees1. *)
stabilizerGS =
Stabilizer[dummies1slots, Last@expr1SGS] /.
Thread@Rule[frees1slots, Range@Length@frees1];
(* If the stabilizer is trivial, we will have all permutations. *)
If[Length@stabilizerGS == 0,
perms1 = Permutations[List @@ frees1],
(* The code for the transversal is probably not optimal. *)
transversal1 = RightTransversal[stabilizerGS, Length@frees1];
(* Translate this to permutations of the same kind as \
Permutations[List@@frees1] gives. *)
perms1 =
TranslatePerm[transversal1, {Perm, Length@frees1}] /.
Perm[{x___}] -> {x} /.
Thread@Rule[Range@Length@frees1, frees1];
];
(* Remove collisions between frees1 and dummies2 *)
arrexpr2 =
xAct`xTensor`Private`arrangedummies[expr2, dummies2, frees1];
DefConstantSymbol[xxx, DefInfo -> False];
Do[
newexpr2 =
xAct`xTensor`Private`changeFreeIndices[arrexpr2, sortedfrees2,
perms1[[count]]];
equation = ToCanonical[expr1 - xxx newexpr2];
solution = First[xxx /. Solve[equation == 0, xxx]];
If[ConstantQ[solution], result = True;
perm = Inner[Rule, perms1[[count]], sortedfrees2, List];
Break[]],
{count, Length[perms1]}]];
{result, solution, sortedfrees2, perm}
];
SetNumberOfArguments[SpecialEqualExpressionsQ, 2];
Protect[SpecialEqualExpressionsQ];
Make a rule that looks for a tensors and replaces it with a properly
reformulated RHS, if tensora[indsa] is equal to tensora[indsaeq].
The expressions are compared with SpecialEqualExpressionsQ
SetAttributes[MakeCompareRule, HoldFirst];
For a single tensor:
MakeCompareRule[{tensora_?xTensorQ[indsaeq___], RHS_}, ___] :=
With[{dummieseq = xAct`xTensor`Private`TakeEPairs[{indsaeq}]},
With[{numdummiesaeq = Length@Intersection[dummieseq, {indsaeq}]},
RuleDelayed[tensora[indsa___],
With[{dummies = xAct`xTensor`Private`TakeEPairs[{indsa}],
frees = xAct`xTensor`Private`TakeFrees[{indsa}]},
Module[{equalexpr, sign, sortedfrees,
perm}, {equalexpr, sign, sortedfrees, perm} =
SpecialEqualExpressionsQ[tensora[indsaeq], tensora[indsa]];
With[{arrexpr2 = ReplaceDummies[RHS]}, (sign*arrexpr2 /. perm)] /;
equalexpr] /;
Length@Intersection[dummies, {indsa}] == numdummiesaeq]]]]
A similar thing for a product of tensors:
MakeCompareRule[{LHS : Times[_?xTensorQ[___], _?xTensorQ[___] ...],
RHS_}, ___] :=
With[{indslist =
ToExpression[StringJoin["tinds", ToString[#]] & /@ Range[Length@LHS]]},
With[{rulelhs = MapIndexed[#1[[0]][First@indslist[[#2]]] &, LHS],
rulelhspattern =
MapIndexed[#1[[0]][
With[{patterninds = First@indslist[[#2]]},
Pattern[patterninds, BlankNullSequence[]]]] &, LHS]},
RuleDelayed[rulelhspattern,
With[{rulelhsinds = FindIndices[rulelhs]},
With[{dummieslhs = xAct`xTensor`Private`TakeEPairs[rulelhsinds],
freeslhs = xAct`xTensor`Private`TakeFrees[rulelhsinds]},
Module[{equalexpr, sign, sortedfrees,
perm}, {equalexpr, sign, sortedfrees, perm} =
SpecialEqualExpressionsQ[LHS, rulelhs];
Module[{arrexpr2 = ReplaceDummies[RHS]}, (sign*arrexpr2 /.
perm)] /;
equalexpr]]]]]]
A code that handles derivatives is missing. I should probably write that
some time.
ApplyCompareRule is used to make rules from equations.
ApplyCompareRule[expr_] := MakeCompareRule[Evaluate[List @@ expr]]
Please tell me if you find the code useful.
Regards
Thomas B�ckdahl
On 2012-02-24 20:45, Edu Serna wrote:
> Hi Thomas,
>
> I doubt my MakeRules are there the problem as I mostly replace the
> contraction of two vectors with scalar quantities (I imagine i could
> write my own cheaper rules but I doubt its worth the effort) as I
> calculate scattering amplitudes. The thing is just that my expressions
> are VERY long (hundreds of thousands of terms) so I need to do it in
> parallel.
>
> I would like to see that code just to see how you did it but I doubt
> that is my problem
>
> Thanks so much anyways.
>
Guillaume Faye
Here is a small addendum to José María's post.
>
> 1) First use Expand[ texpr ] if there are nonexpanded products of sums
> of tensors.
> 2) Use SameDummies on the result of 1). This will try to simplify the
> expression by minimizing the number of different dummies, without
> actually calling group-theoretical algorithms.
> 3) Map ToCanonical to each term, or use a Do loop to monitor timings
> of individual canonicalizations. If the result is much smaller than
> the input, it might be useful to add intermediate results to an
> accumulator variable, instead of doing a big sum at the end.
It is very important, in item 3, when using a Do loop (which seems to be less
memory consuming in general), to store the result of ToCanonical acting on
each term into a list, and not a Plus-type expression. Indeed, the command
Plus tries to perform operations that can take time (reordering, flattening,
etc.) if a sum is very lengthy. It is thus much more efficient not to add
canonicalized terms to partials sum at each step, but to achieve the summation
only at the end.
I also advise to prepare the expressions with as many Scalar as possible.
For very long sums of monomials made of a lot of factors, this can make the
computation less memory consuming, specially when the tensors that are involved
in the expressions have only one or two indices (scalar products of vectors).
Using the strategy proposed by José María and the two tricks described above,
I can manipulate hundreds of thousand tensorial terms in a few minutes. I need
to make parallel computation only when I have to apply to each term of my sum
a very time consuming operation (typically an integration).
Yours,
Guillaume.
Edu Serna
Thanks everyone for all of this, I now have faith that I will be able
to run this reasonably fast.
I am gonna spell out exactly what I am doing so I can clarify some of
the points raised, and hopefully for other people's future reference.
First of all I need to expand products of reasonably long expressions.
This I have to parallelize or it takes for ever so I take the longest
factor in my expression and make a list of it where each term is of
some reasonable length:
split =
Append[ ParallelTable[
Take[LONGFACTOR, {Floor[L/M] i + 1, Floor[L/M] (i + 1)}], {i, 0,
M - 1}], Take[LONGFACTOR, {Floor[L/M] M + 1, L}]];
with L=Length[LONGFACTOR] and M so that each term is of some
reasonable length.
I then calculate my amplitude as
DistributeDefinitions[ContractMetric];
EXPANDEDEXPR=ParallelMap[(ContractMetric[Expand[#]]) &, (LONGFACTOR
RESTEXPR)]
I gather that according to Gillaume I should now do
Apply[List, EXPANDEDEXPR]
so that each term of the sum becomes each term of the list
What I want to do now is map each term that will have contracted terms
of the form (in parallel I assume, I hope the overheads of treating
each small term of the list in parallel are outweighed):
v1[a] v2[b]..v1[-a] v2[-b]...vi[m] vj[n]
using /.MakeRule[{v1[a] v1[-a], Scalar1[] }, MetricOn -> All,
ContractMetrics -> True];
so that my term end up something like
Scalar1[].. ScalarN[] vi[m] vj[n]
I was wondering whether I should do
> 2) Use SameDummies on the result of 1). This will try to simplify the
> expression by minimizing the number of different dummies, without
> actually calling group-theoretical algorithms.
also is there a cheaper equivalent to /.MakeRule[{v1[a] v1[-a],
Scalar1[] }, MetricOn -> All, ContractMetrics -> True]; considering
that there no symmetries or just using something like /.v1[_] v1[-_] -
> Scalar1[]
I think I should avoid using ToCanonical altogether?
Guillaume Faye
Dear Edu,
>
> using /.MakeRule[{v1[a] v1[-a], Scalar1[] }, MetricOn -> All,
> ContractMetrics -> True];
>
> so that my term end up something like
>
> Scalar1[].. ScalarN[] vi[m] vj[n]
>
> I was wondering whether I should do
This is a possible solution. You may also want to use the general command
PutScalar to introduce Scalar objects. It is worthy when your monomials
involve factors Scalar1[]^q1 ... ScalarN[]^qn with q1, ... or qn much greater
than 1. You may alternatively prefer having your own faster version of
PutScalar specifically adapted to your problem (I chose the latter approach).
>
>> 2) Use SameDummies on the result of 1). This will try to simplify the
>> expression by minimizing the number of different dummies, without
>> actually calling group-theoretical algorithms.
As for me, I combine both strategies (Scalar approach + avoiding to use
ToCanonical except at the very final stage, if needed). This makes my
computations *at least* twice faster (for lengthy expressions). [The time
reduction factor can be significantly greater than 2.]
Yours,
Guillaume.
Edu Serna
Thanks so much for all this, I will put in practice and post my
experience here.
Edu
Filippo
Many thanks,
Thomas Bäckdahl
The problem here is that the different kernels have their own lists of defined tensors, abstract indices, covds etc. These are created when you load xTensor in each kernel, and they are not updated from the data in the main kernel. If you don't give an extra argument to SortCovDs, it will look in the list $CovDs, and do the sorting for all covariant derivatives it finds there. In your parallel kernels it will not find CD so it will not sort it.
One way to fix this is to remove the local variables in each parallel kernel. This means that when some code is evaluated in a parallel kernel, the main kernel variables are used. I think the main names and setting should be in the $ variables in the xTensor context. Hence, the following code should help:
xTensorDollarNames = Names["xAct`xTensor`$*"]
ParallelEvaluate[Clear @@ #, Kernels[]] &@xTensorDollarNames
If it turns out that this becomes slow because these variable names have to be fetched from the main kernel, one can copy the current values of these variables into each kernel.
ParallelEvaluate[Set[Evaluate@Symbol@#[[1]], #[[2]]];,
Kernels[]] & /@ ({#, Symbol[#]} & /@ xTensorDollarNames);
Observe that if you define a new tensor or makes something else that changes the $ variables, you have to clear the local variables again, and copy the data one more. If you don't copy the data, but just clear the local variables, everything should update automatically.
There where some annoying messages that can be turned off by:
Off[RuleDelayed::rhs]
PrintAsCharacter[] = "";
I don't know where these messages came from though.
Regards
Thomas
--
You received this message because you are subscribed to the Google Groups "xAct Tensor Computer Algebra" group.
To unsubscribe from this group and stop receiving emails from it, send an email to xact+uns...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Filippo
Filippo
The method is in fact normally slower than a single kernel calculation. If I understand well the names are taken from the main Kernel every time a call is sent to a subKernel.
However I did not understand how to use the second method, i.e. to copy names definitions to each kernel.
If I do
xTensorDollarNames = Names["xAct`xTensor`$*"]
ParallelEvaluate[Clear @@ #, Kernels[]] &@xTensorDollarNames
And also by evaluating the second ParallelEvaluate leads to a lot of errors.
Cheers,
Filippo
Thomas Bäckdahl
My code was not very good, but it actually worked for me. An alternative code to set the $ variables in each kernel would be
ParallelEvaluate[SetFromName[str_String, val_] := ToExpression[str, InputForm, Function[name, name = val, HoldAll]], Kernels[]]
ParallelEvaluate[SetFromName @@ #;, Kernels[]] & /@ ({#, Symbol[#]} & /@ xTensorDollarNames);
I think a better alternative might be to run all Define commands in all kernels as in the attached example. If you have many definitions for objects that are needed in your parallel code, it might be a good idea to make an .m file with these definitions, and load that .m file in each kernel. With this philosophy you don't distribute any data. Let us say that you want to compute a function F operating on an expression E involving A. Then you define both F and A in each kernel separately and give them all the properties they should have before you make the parallel function call. This way only the expression E needs to be passed from the main kernel to one of the parallel kernels. No properties needs to be passed. Because of this we can use the setting $DistributedContexts = None.
Regards
Thomas
