SymManipulator

Symmetrized tensor expressions

Thomas Bäckdahl     thomas.backdahl@chalmers.se
Mathematical Sciences,
Chalmers University of Technology and University of Gothenburg
Gothenburg, Sweden

(c) 2011-2021, under GPL

http://www.xAct.es/
http://groups.google.com/group/xAct

SymManipulator is a package for handling symmetrized tensor expressions and making irreducible decompositions of spinor expressions.

SymManipulator is distributed under the GNU General Public License, and runs on top of xTensor, a free package for fast manipulation of abstract tensor expressions. Both can be downloaded from http://www.xAct.es/

DateList[]

xAct`SymManipulator`$Version={"0.9.5",{2021,09,14}}

xAct`SymManipulator`$xTensorVersionExpected = {"1.1.5", {2021, 8, 29}};

1. Initialization

1.0. TODO list

*The canonicalization of general groups don’t not move symmerizations exactly as we would like. A new algorithm is desperately needed to speed it up.
*Allow non-spinor VBundles in the fundamental spinor operators.

1.1. GPL

(* SymManipulator: Symmetrized tensor expressions *)

(* Copyright (C) 2011-2021 Thomas Bäckdahl *)

(* This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License,or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place-Suite 330, Boston, MA 02111-1307, USA. 
*)

1.2. Info package

(* :Title: SymManipulator *)

(* :Author: Thomas Bäckdahl *)

(* :Summary: Symmetrized tensor expressions and irreducible decomposition 
    of spinor expressions *)

(* :Brief Discussion:
   - Handles symmetrized tensor expressions. It also handles irreducible decompositions of spinor expressions, and fundamental spinor operators. *)
  
(* :Context: xAct`SymManipulator` *)

(* :Package Version: 0.9.5 *)

(* :Copyright: Thomas Bäckdahl (2011-2021) *)

(* :History: See SymManipulator.History *)

(* :Keywords: *)

(* :Source: SymManipulator.nb *)

(* :Warning: *)

(* :Mathematica Version: 9.0 and later *)

(* :Limitations:
	- Only abstract indices can be used. *)
	
(* :Acknowledgements:
	This work was partly funded by scholarships from the Albert Einstein 
	Institut and the Wenner-Gren foundarions. Furthermore the autor would like to 
	thank J. M. Martín-García, A. García-Parrado, Teake Nutma, Leonard Soicher 
	and Peter Cameron for helpful comments. *)

1.3. BeginPackage

If[Unevaluated[xAct`xCore`Private`$LastPackage]===xAct`xCore`Private`$LastPackage,xAct`xCore`Private`$LastPackage="xAct`SymManipulator`"];

xAct`xCore`Private`$LastPackage

Explicit (not hidden) import of xTensor, xPerm and xCore:

BeginPackage["xAct`SymManipulator`",{"xAct`xTensor`","xAct`xPerm`","xAct`xCore`"}]

Print[xAct`xCore`Private`bars]
Print["Package xAct`SymManipulator`  version ",$Version[[1]],", ",$Version[[2]]];
Print["CopyRight (C) 2011-2021, Thomas Bäckdahl, under the General Public License."];

Check if we have the correct version of xTensor.

If[Not@OrderedQ@Map[Last, {xAct`SymManipulator`$xTensorVersionExpected, xAct`xTensor`$Version}], 
	Message[General::versions, "xTensor", xAct`xTensor`$Version, xAct`SymManipulator`$xTensorVersionExpected];
	Abort[]
];

We specify the context xAct`SymManipulator` to avoid overriding the Disclaimer of the other packages. However we need to turn off the message General:shdw temporarily:

Off[General::shdw]
xAct`SymManipulator`Disclaimer[]:=Print["These are points 11 and 12 of the General Public License:\n\nBECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM `AS IS´ WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.\n\nIN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES."]
On[General::shdw]

If SymManipulator is not being called from other package then write this GPL short disclaimer:

If[xAct`xCore`Private`$LastPackage==="xAct`SymManipulator`",
Unset[xAct`xCore`Private`$LastPackage];
Print[xAct`xCore`Private`bars];
Print["These packages come with ABSOLUTELY NO WARRANTY; for details type Disclaimer[]. This is free software, and you are welcome to redistribute it under certain conditions. See the General Public License for details."];
Print[xAct`xCore`Private`bars]]

Note that symbols in the Global` context cannot be accessed right now.

$ContextPath

$Context

Established connection to external executable?

$xpermQ

1.4. Non-standard setup

Screen all dollar indices:

$PrePrint=ScreenDollarIndices;

Timings will be shown if they are above 0.1 Second (Only for this notebook; this is not included in the package):

<<xAct/ShowTime1.m

$ShowTimeThreshold=0.1

We also read the package (this is not automatic in xCore anymore since its version 0.6.2):

<<xAct/ExpressionManipulation.m

1.5. Usage messages

SymH::usage="SymH[headlist, sym, label][indices] represents a symmetrized tensor expression if the product of the tensors in headlist. The imposed symmetry is sym, and label displayed under the symbol. Observe that traces are taken after the symmetrization. These operations do not generally commute unless they work on disjoint sets of indices.";
ImposeSym::usage="ImposeSym[expr, inds, sym] creates a symmetrized expression equivalent to ImposeSymmetry[expr, inds, sym]. The result is displayed as SymH objects. If expr, already has the symmetry sym, expr is returned as is is. The output is 0 if the result due to single terms symmetries can be reduced to zero. Warning if expr contains traces, the resulting expression will be an expression where the traces are taken after the symmetization. The trace and symmetriation operations do not generally commute unless they work on disjoint sets of indices. The user needs to check that symmetrizations are not taken on contracted indices. Warning, ImposeSym requires expr to be expanded. For general expressions use ImposeSym[Expand@expr, inds, sym]";
ExpandSym::usage ="ExpandSym[expr] expands the SymH objects in expr. If the option SmartExpand is used, a smaller number of terms are produced by using the interplay between the internal symmetry and the imposed symmetry.";
RemoveSym::usage ="RemoveSym[expr] removes the SymH objects in the sence that the first term of ExpandSym[expr] is computed. Observe that traces of the symmetrized expressions are taken after the symmetrization. The trace and symmetriation operations do not generally commute unless they work on disjoint sets of indices. The user needs to check that symmetrizations are not taken on contracted indices.";
RemoveSuperfluousSym::usage ="RemoveSuperfluousSym[expr] if an expression inside a SymH object already has the imposed symmetry, the extra symmetrization will be removed.";
RemoveSuperfluousInnerSym::usage ="RemoveSuperfluousInnerSym[expr] if an expression has nested SymH symbols, and an inner symmetrization is a part of an outer symmetrization, the inner symmetrization will be removed. For instance if an expression is symmetrized over slot 1 and 2, and then symmetrized over slots 1, 2, and 3, then the symmetrization over slot 1 and 2 can be removed.";
RemoveTrivialSym::usage = "RemoveTrivialSym[expr] replaces SymH objects in expr with 0 if they vanish due to single terms symmetries.";
MoveTensorsOutsideSym::usage ="MoveTensorsOutSideSym[expr] moves all tensors that are not part of the symmetry outside the SymH object.";
MoveTensorsInsideSym::usage = "MoveTensorsInsideSym[expr] moves all tensors that inside the SymH object.";
CanonicalizeGroupInSym::usage ="CanonicalizeGroupInSym[expr] finds a canonical representaion of all SymH objects. Observe that the indices are not canonicalized, so the correct usage for complete canonicalization is ToCanonical@CanonicalizeGroupInSym[expr].";
SubgroupQ::usage ="SubgroupQ[sym1, sym2] returns true iff sym1 is a subgroup of sym2.";
CompatibleSymmetric::usage ="CompatibleSymmetric[inds] splits inds into sets accoring to vbundle. The result is then the disjoint union of the symmetric groups on these sets.";
ExpandSymOneIndex::usage ="ExpandSymOneIndex[expr, ExpansionIndex] expands expr so that ExpansionIndex is not part of the symmetry.";
InertHeadHead::usage="InertHeadHead[inerthead, tensorh] a representaiton of a tensor head such that InertHeadHead[H1,T][inds]=H1[T[inds]].";
CovarD::usage ="CovarD[CD, tensorh ,slotlist] a representaiton of a tensor head such that CovarD[CD,T,{-a,-b,...}][inds]=CD[-a, -b,...][T[inds]]].";
ZeroTensor::usage ="ZeroTensor[slotlist_] is a representation of the zero tensor with slots defined by slotlist. That is ZeroTensor[{TangentM, ...}][-a, ...]=0.";
TensorToZeroRule::usage="TensorToZeroRule[tensorhead] makes a replacement rule tensorhead->ZeroTensor[slotlist] with the aproproate slotlist.";
MoveSymIndicesDown::usage ="MoveSymIndicesDown[expr] moves the indices in the SymH objects so that the symmetrized indices are in a down position. This position is required for a good TexPrint output. Warning, this function requires that all indices can be moved with the metric. This check is not done automatically.";
SortTensorsInSym::usage = "SortTensorsInSym[expr] sorts the tensors in the SymH objects alphabetically.";
IrrDecompose::usage = "IrrDecompose[expr] computes the irreducible decomposition of over all spinor indices.";
CompleteIrrDecompose::usage = "CompleteIrrDecompose[expr] computes the irreducible decomposition of over all spinor indices. The symmetries of expr is then imposed on the resulting expression. This can also be done with the option UseSym in IrrDecompose.";
UseSym::usage="UseSym is an option for IrrDecompose. If this option is set to True, then the symmetries of the expression is omposed on the decomposed expression. This often simplifies the expression.";
IrrDecImposeNewMethod::usage ="An option for IrrDecompose.";
ResultType::usage="An option for IrrDecompose. The default is\"Expression\". Other possibilities are \"Equation\", \"Rule\" or \"CompareRule\".";
InternalCommutingSymmetry::usage = "InternalCommutingSymmetry is a tag that describes the part of the internal symmetry of a SymH object that commutes with the imposed symmetry.";
SmartExpand::usage ="An option for ExpandSym.";
ImposeLargerSym::usage ="An option for ImposeSym.";
ImposeSuperfluousSym::usage ="An option for ImposeSym.";
deltaH::usage ="deltaH[vbundle] is the delta on the vector bundle vbundle.";
ContractMetricsInsideSym::usage ="ContractMetricsInsideSym[expr] checks if any SymH object has a metric inside that can be contracted, and contracting this.";
SmartSymmetrize::usage ="Does the same thing as Symmetrize, but reducing the number of terms due to the internal symmetry.";
SmartAntisymmetrize::usage ="Does the same thing as Antisymmetrize, but reducing the number of terms due to the internal symmetry.";
ExpandInternalSym::usage="Applies ExpandSym to the internal expression of a SymH object.";
ToCanonicalSym::usage="Runs CanonicalizeGroupInSym followed by ToCanonical.";
RuleInSym::usage ="Applies a rule inside a SymH object. Observe that this can decanonicalize all SymH objects.";
FunctionInSym::usage="Applies a function inside a SymH object. Observe that this can decanonicalize all SymH objects.";
$MixedSymSpecialCode::usage = "Default True, General code False. Old behaviour SetStabilizer. Observe that SetStabilizer gives only a subgroup of the correct group, but it is faster.";
TransversalInSymmetricGroup::usage ="TransversalInSymmetricGroup[H_StrongGenSet,G_Symmetric] computes the transversal of H in G. If H is not a subgroup, an error is thrown. (Not tested well.)";
$MixedSymVerbose::usage ="Display timing info from MixedSymmetry calculations.";
SymmetryOfExpression::usage ="SymmetryOfExpression[expr] computes the symmetry group of the free indices of expr and returns a strong generating set and a list of the free indices.";
MakeCompareRule::usage = "MakeCompareRule is used in the same way as MakeRule, but only one rule is generated. When applied to an expression, it is compared to the LHS. If a match is found it is replaced by the corresponding RHS. At the moment the LHS can only be a single tensor, a product of tensors, or a covariant derivative of a single tensor.";
CompatibleSymQ::usage ="CompatibleSymQ[vbundles, sym] returns true if the total symmetry is a subgroup of sym. Each vbundle gives a symmetric group. The total symmetry if the union of these groups. Signs of vbundles are dropped before the rest of the computation.";
RemoveSuperfluousPartialSym::usage = "Tries to factor the imposed group, and remove superfluous factors. It is also an option for ToCanonicalSym.";
IrrDecSpinBundles::usage ="An option for IrrDecompose and CompleteIrrDecompose.";
ExpandBox::usage ="ExpandBox[expr,CDe] converts the BoxCDe in expr to curvature.";
DefFundSpinOperators::usage ="DefFundSpinOperators[CDe] defines the operators DivCDe, CurlCDe, CurlDgCDe and TwistCDe what acts on symmetric spinors. Also commutator rules like CommuteDivCurlCDe are defined.";
HideCovDSymbolReference::usage="HideCovDSymbolReference is an option for DefFundSpinOperators. If set to True, the default PrintAs command for the derived symbols will not contain a bracket with the covd symbol.";
ShowValenceInfo::usage="ShowValenceInfo is an option for DefFundSpinOperators. If set to True, the default PrintAs command for the operators will contain numbers indicating the valence of the spinor it acts on.";
UndefFundSpinOperators::usage ="UndefFundSpinOperators[CDe] undefines all objects defined by DefFundSpinOperators[CDe].";
FundSpinOpQ::usage="FundSpinOpQ[op] returns True if the symbol op is a fundamental spinor operator.";
CovDOfFundSpinOp::usage="CovDOfFundSpinOp[op] returns the covariant derivative associated with the fundamental spinor operator op.";
ImposeFundSpinOp::usage="ImposeFundSpinOp[expr, op, {inds}] applies the fundamental spinor operator op to expr and gives the result the indices inds.";
ExpandFundSpinOp::usage="ExpandFundSpinOp[expr, covd] expands the fundamental spinor operators related to covd into symmetrized derivatives.";
FunctionInFundSpinOp::usage="FunctionInFundSpinOp[expr, func, l, outerfunc] applies the function func, inside the fundamental spinor operators at level l. The optional function outerfunc is applied to the results at all lower levels. The default level is 1 and the default outerfunc is Identity.";
IrrDecomposeCovD::usage="IrrDecomposeCovD[covd[...]@T[...]] produces the irreducible decomposition of covd[...]@T[...] using the corresponding fundamental operators.
IrrDecomposeCovD[T, covd] does the same, but you do not have to give any indices.";
ToFundSpinOp::usage="ToFundSpinOp[expr, covd] tries to represent all covariant spinor derivatives in expr in terms of the fundamental spinor operators. When several nested derivatives appear, one might need to apply the function several times.";
CommuteFundSpinOp::usage="CommuteFundSpinOp[op1,op2] gives the commutator rule for the fundamental spinor operators op1, op2. Observe that this function only gives the proper commutator relations. There are other commutator like relations with two derivative terms in the right hand side. CommuteFundSpinOp[op1,op2,op3,op4] gives a commutator like rule where op1[op2[expr]] is turned into an expression with a term op3[op4[expr]] if possible.";
GetIndexRange::usage="GetIndexRange[k, vbundle] gives the first k indices of vbundle. New indices are constructed if needed.";
GiveIndicesToTensor::usage ="GiveIndicesToTensor[TT] gives the tensor TT down indices in the correct vbundles. A list of vbundles can be given as a second argument. Indices in these vbundles will then be up indices.";
ChangeFreeIndicesToDefault::usage ="ChangeFreeIndicesToDefault[expr] changes the free indices of expr. The indices are taken from the beginning of the list of indices for the corresponding VBundle. A list of VBundles can be given as a second argument. Free indices belonging to these VBundles are given up indices. The default is that all new free indices are down. The dummy indices are also changed.";
MoveCovDInsideSym::usage=="MoveCovDInsideSym is a rule that moves covariant derivatives inside SymH objects.";
FindTensorCoefficients::usage=="FindTensorCoefficients[expr,{T1[...], T2[...]T3[...]}] finds the coefficients of the tensor expressions T1[...] and T2[...]T3[...] in expr. The result is a list of the form {{rest, 1},{c1, T1[...]},{c2, T2[...]T3[...]}}. Plus@@Times@@@FindTensorCoefficients[expr, {...}] is equivalent to expr, but with tensors collected." ;
SpatialSpinCovDQ::usage ="SpatialSpinCovDQ Reserved for spatial spinors.";
PlaceIndicesInTensor::usage ="PlaceIndicesInTensor[T,{inds}] places the indices inds in their correct positions in the tensor T, where the correct positions are given by SlotsOfTensor.";
ValenceNumbersOfSpinor::usage="ValenceNumbersOfSpinor[T, Spin], gives a list {k,l} where k is the number of Spin slots and l is the number of Spin† slots of the spinor T.";
IndepIrrDecompose::usage="IndepIrrDecompose[expr, options] computes the irreducible deomposition of expr, tries to find relations between the irreducible parts, and reduces the expression with respect to these relations. If all such relations wever found the remainign irreducible parts should be independent. However, the algorithm is not analysed for completeness, so this can not be guanranteed of several irreducible pieces have the same valence.";
FindSpinorsInIrrDec::usage="FindSpinorsInIrrDec[expr, {ϵ, ϵ†}] gives the irreducible parts of expr assuming that expr is the result of IrrDecompose or IndepIrrDecompose.";
DefIrreducibleSpinors::usage="DefIrrDecSpinors[irrdecspinors, basename, baseprintas ,spin]";
EqToCompareRule::usage="EqToCompareRule[eq] takes an equation eq and tranforms it into a replacement rule with MakeCompareRule.";
ExtractUnprimedAndPrimed::usage="ExtractUnprimedAndPrimed[Spin] returns a function that extract the unprimed and primed indices from a list.";
NumOfUnprimedVBundles::usage="NumOfUnprimedVBundles[Spin] returns a function that counts the number of unprimed spin bundles in a list.";
NumOfPrimedVBundles::usage="NumOfPrimedVBundles[Spin] returns a function that counts the number of primed spin bundles in a list.";
NumOfUnprimedSlots::usage="NumOfUnprimedSlots[Spin] returns a function that counts the number of unprimed slots of a spinor.";
NumOfPrimedSlots::usage="NumOfPrimedSlots[Spin] returns a function that counts the number of primed slots of a spinor.";
NumOfUnprimedSlotsDyad::usage="";
NumOfPrimedSlotsDyad::usage="";
SymmetricSpinorOfValenceQ::usage="";
SymmetricSpinorOfArbitraryValenceQ::usage="";
ExpandSymS::usage="ExpandSymS[expr] is a shorthand for ExpandSym[expr, SmartExpand->True]";
EnoughIndsQSpin::usage="EnoughIndsQSpin[T,k,l,Spin] returns True if the spinor T has at least k unprimed and l primed indices with respect to the spin bundle Spin.";
SymmetricOfValence::usage="SymmetricOfValence[k,l,spin,ovb] gives the slot VBundles of a symmetric spinor of valence k,l, possibly with extra VBundles ovb.";
ZeroTensorOfTensor::usage="";
EqToRule::usage="EqToRule[lhs == rhs] transforms the equation into an rule using MakeRule. This is used for indexed expressions assuming that lhs only contains one term.";
FlipEquation::usage="FlipEquation[lhs == rhs] flips the lhs and rhs. If rhs contains a constant factor, the equation is divided by this constant.";
Isolate::usage="Isolate[eq, p] tries to rewrite the equation eq so that the pattern p becomes the left hand side.";

1.6. Begin private

Begin["`Private`"]

The reserved words in version 0.9.1 :

Names["xAct`SymManipulator`*"]

{"CanonicalizeGroupInSym","ChangeFreeIndicesToDefault","CompatibleSymmetric","CompatibleSymQ","CompleteIrrDecompose","ContractMetricsInsideSym","CovarD","CovDOfFundSpinOp","DefFundSpinOperators","deltaH","Disclaimer","ExpandBox","ExpandFundSpinOp","ExpandInternalSym","ExpandSym","ExpandSymOneIndex","FindTensorCoefficients","FunctionInFundSpinOp","FunctionInSym","FundSpinOpQ","GetIndexRange","GiveIndicesToTensor","HideCovDSymbolReference","ImposeFundSpinOp","ImposeLargerSym","ImposeSuperfluousSym","ImposeSym","InertHeadHead","InternalCommutingSymmetry","IrrDecImposeNewMethod","IrrDecompose","IrrDecomposeCovD","IrrDecSpinBundles","MakeCompareRule","MoveCovDInsideSym","MoveSymIndicesDown","MoveTensorsInsideSym","MoveTensorsOutsideSym","RemoveSuperfluousInnerSym","RemoveSuperfluousPartialSym","RemoveSuperfluousSym","RemoveSym","RemoveTrivialSym","ResultType","RuleInSym","SmartAntisymmetrize","SmartExpand","SmartSymmetrize","SortTensorsInSym","SubgroupQ","SymH","SymmetryOfExpression","TensorToZeroRule","ToCanonicalSym","ToFundSpinOp","TransversalInSymmetricGroup","UndefFundSpinOperators","UseSym","ZeroTensor","$MixedSymSpecialCode","$MixedSymVerbose","$Version","$xTensorVersionExpected"}

%//Length

$ContextPath

$Context

Now this is not being used. Deactivated cell :

BeginExamples[]:=(Set[$Context,"Global`"];Global`symbols=Names["Global`*"];Null)
Global`EndExamples[]:=(If[#=!={},Print["There are new symbols in Global`:",#]];Set[$Context,"xAct`SymManipulator`Private`"];Null)&@Complement[Names["Global`*"],Global`symbols]

2. Group theory

2.0. MathLink

MathLinkBaseChangeStabilizerChain[sgs_,newbase_]:=MathLinkBaseChangeStabilizerChain[sgs,newbase,Max[PermLength@sgs,Max@newbase]]

MathLinkBaseChangeStabilizerChain[sgs_,newbase_,len_]:=
Module[{tmp},
tmp=Drop[xAct`xPerm`Private`TMPHead[
len+2,
xAct`xPerm`Private`tosgslist[sgs,len,False],
newbase
],{2}];
TranslatePerm[xAct`xPerm`Private`ToSign[#,len],NotationOfPerm[sgs]]&/@Apply[xAct`xPerm`Private`MLBaseChangeStabilizerChain,tmp]];

MathLinkBaseChange[sgs_,newbase_,len_,options___]:=
Module[{tmp},
tmp=Drop[xAct`xPerm`Private`TMPHead[
len+2,
xAct`xPerm`Private`tosgslist[sgs,len,False],
newbase
],{2}];
TranslatePerm[xAct`xPerm`Private`ToSign[Apply[xAct`xPerm`Private`MLBaseChange,tmp],len],NotationOfPerm[sgs]]
];

Options[BaseChange2]={MathLink->MemberQ[LinkPatterns[$xpermLink]/.HoldForm->Head,xAct`xPerm`Private`MLBaseChange]};

BaseChange2[sgs_StrongGenSet,newbase_]:=BaseChange2[sgs,newbase,Max[PermLength@sgs,Max@newbase]]

BaseChange2[StrongGenSet[base_, gs_],base_,len_]:=StrongGenSet[base, gs];

BaseChange2[sgs_,newbase_,len_,options:OptionsPattern[]]:=If[OptionValue[MathLink],MathLinkBaseChange,BaseChange][sgs,newbase,len,options];

Observe that the current MathLink code keeps the order of the rest of the base as the Mathematica impementation.

MathLinkStabilizerSGS[pts_,sgs_,len_]:=
Module[{tmp},
tmp=Drop[xAct`xPerm`Private`TMPHead[
len+2,
xAct`xPerm`Private`tosgslist[sgs,len,False],
pts
],{2}];
TranslatePerm[xAct`xPerm`Private`ToSign[Apply[xAct`xPerm`Private`MLStabilizerSGS,tmp],len],NotationOfPerm[sgs]]
];

Options[StabilizerSGS]={MathLink->MemberQ[LinkPatterns[$xpermLink]/.HoldForm->Head,xAct`xPerm`Private`StabilizerSGS]};

StabilizerSGS[points_List,SGS_,options:OptionsPattern[]]:=StabilizerSGS[points,SGS,PermLength[SGS],options];
StabilizerSGS[ps_List,Symmetric[list_],options:OptionsPattern[]]:=Symmetric[Complement[list,ps]];
StabilizerSGS[ps_List,Antisymmetric[list_],options:OptionsPattern[]]:=Antisymmetric[Complement[list,ps]];
StabilizerSGS[points_List,SGS:StrongGenSet[base_,GS_],len_Integer,options:OptionsPattern[]]:=If[OptionValue[MathLink],MathLinkStabilizerSGS[points,SGS,len],Stabilizer[points,SGS,len]];

2.1. Basic group and set theory

We need the internal function PartitionRagged. For old versions of Mathematica, we need to supply an implementation.

oldPartitionRagged[l_,p_]:=Last/@Rest@FoldList[{Drop[#1[[1]],#2],Take[#1[[1]],#2]}&,{l,{}},p];

If[System`$VersionNumber<8.,partitionRagged=oldPartitionRagged,partitionRagged=Internal`PartitionRagged];

Check if a set is a subset of the other. Observe that a new function in Mathematica 10 has the same name and meaning, but a different order of the arguments.

If[System`$VersionNumber<10.,subsetQ[inds1_,inds2_]:=And@@(MemberQ[List@@inds2,#]&/@(List@@inds1)),subsetQ[inds1_,inds2_]:=SubsetQ[inds2,inds1]];

Check if a group is a subgroup of the other.

Symmetric cases.

SubgroupQ[sym1:(_Symmetric|_Antisymmetric),sym2_StrongGenSet]:=SubgroupQ[GenSet@sym1,sym2]

SubgroupQ[sym1_,sym2:(_Symmetric|_Antisymmetric)]:=SubgroupQ[sym1,Append[sym2,Cycles]]

For general cases we need to check that the generators of the first group is in the other group.

SubgroupQ[sym1_StrongGenSet,sym2_StrongGenSet]:=SubgroupQ[Last@sym1,sym2]

SubgroupQ[gs1_GenSet,sym2_StrongGenSet]:=And@@(PermMemberQ[#,sym2]&/@List@@gs1)

Computation of the group support, i. e. the points that are moved by some element in the group.

GroupSupport[sym:(_Symmetric|_Antisymmetric)]:=Sort@First@sym

GroupSupport[gs1_GenSet]:=xAct`xTensor`Private`support[gs1]

GroupSupport[SGS1_StrongGenSet]:=xAct`xTensor`Private`support[Last@SGS1]

GroupSupport[sym_]:=xAct`xTensor`Private`support[Last@xAct`xTensor`Private`SGSofsym@sym]

Find a nicer representation of the group.

TidySGS[sym_]:=DeleteSomeRedundantGenerators@BaseChange2[xAct`xTensor`Private`SGSofsym@sym,GroupSupport@sym]

TidyGroup[sym_]:=Module[{unstablepoints=GroupSupport@sym},
Which[
Length@unstablepoints==0,StrongGenSet[{},GenSet[]],
SubgroupQ[Symmetric@unstablepoints,sym],Symmetric@unstablepoints,
SubgroupQ[Antisymmetric@unstablepoints,sym],Antisymmetric@unstablepoints,
True,TidySGS@sym]]

Try to factor the group.

GroupFactorization1[SGS_StrongGenSet]:=Module[{orbits=Orbits@SGS,factors,prodgroup},
factors=Stabilizer[Flatten@Delete[orbits,#],SGS]&/@Range@Length@orbits;
prodgroup=JoinSGS@@factors;
(* If the group factors completely into ints orbits, return the factorization, otherwise return the group as it is. If we have several orbits the second case can be improved... *)
If[And[SubgroupQ[SGS,prodgroup],SubgroupQ[prodgroup,SGS]],factors,{SGS}]
]

2.2. More group theory

We implement a new code for the stabilizer chain which is much quicker.

StabilizerChain2[StrongGenSet[base_List,GS_GenSet]]:=Thread@StrongGenSet[Drop[base,#]&/@Range[0,Length@base],FoldList[Stabilizer[{#2},#1]&,GS,base]]

Here we extend the base change algorithm to work directly on the stabilizer chain instead. This can actually make it quicker.

BasicOrbitStabChain[j_Integer,k_Integer,chain_List]:=Orbit[chain[[1,1,j]],chain[[k+1]]]/;k<j;

This is algorithm 12.3 of Butler with some slight modifications so it now works on the stabilizer chain.

InterchangeStabChain[chain_List,j_Integer,len_Integer]:=
With[{base=First@First@chain},Module[{Deltaj,Deltajp1,LDeltaBarjp1,T,Gamma,Delta,gamma,p,g1,g2,newT, newbase,newGSj},
Deltaj=BasicOrbitStabChain[j,j-1,chain];
Deltajp1=BasicOrbitStabChain[j+1,j,chain];
LDeltaBarjp1=Length[Deltajp1]Length[Deltaj]/Length[BasicOrbitStabChain[j+1,j-1,chain]];
T=chain[[j+2,2]];
Gamma=Complement[Deltaj,base[[{j,j+1}]]];
Delta={base[[j]]};
While[Length[Delta]≠LDeltaBarjp1,
gamma=First[Gamma];
g1=TraceSchreier[gamma,SchreierOrbit[base[[j]],chain[[j]],len]];
p=OnPoints[base[[j+1]],InversePerm[g1]];
If[MemberQ[Deltajp1,p],
g2=TraceSchreier[p,SchreierOrbit[base[[j+1]],chain[[j+1]],len]];
(* Do we need to check if PermProduct[g2,g1] is already in T before we add it? *)
AppendTo[T,PermProduct[g2,g1]];
Delta=Orbit[base[[j]],T];
Gamma=Complement[Gamma,Delta],
Gamma=Complement[Gamma,Orbit[gamma,T]]
]
];
newT=DeleteCases[Drop[T,Length@chain[[j+2,2]]],Alternatives@@chain[[j,2]]];
newGSj=Stabilizer[{base[[j+1]]},Union[chain[[j,2]],newT]];
newbase=Join[Take[base,j-1],{base[[j+1]],base[[j]]},Drop[base,j+1]];
Join[Thread@StrongGenSet[Drop[newbase,#]&/@Range[0,j],Append[Join[#,newT]&/@chain[[1;;j,2]],newGSj]],Drop[chain,j+1]]
]]

ConjugateStabChain[chain_List,g_?PermQ]:=Module[{newbase=OnPoints[First@First@chain,g],oldGS=Last@First@chain,newGS,l=Length@chain},
newGS=PermProduct[InversePerm[g],#,g]&/@oldGS;
Thread@StrongGenSet[Drop[newbase,#]&/@Range[0,l-1],(#[[2]]&/@chain)/.Thread@Rule[List@@oldGS,List@@newGS]]]

AppendBasePointStabChain[chain_List,p_]:=Append[MapAt[Append[#,p]&,chain,{#,1}&/@Range@Length@chain],StrongGenSet[{},GenSet[]]]

Options[BaseChangeStabChain]={MathLink->MemberQ[LinkPatterns[$xpermLink]/.HoldForm->Head,xAct`xPerm`Private`MLBaseChangeStabilizerChain]};

BaseChangeStabChain[SGS_StrongGenSet,base_,options:OptionsPattern[]]:=BaseChangeStabChain[SGS,base,PermLength[SGS],options];
BaseChangeStabChain[SGS_StrongGenSet,newbase_List,len_,options:OptionsPattern[]]:=If[OptionValue[MathLink],MathLinkBaseChangeStabilizerChain[SGS,newbase,len],BaseChangeStabChain2[StabilizerChain2[SGS],newbase,len]];
BaseChangeStabChain[chain_List,base_,options:OptionsPattern[]]:=BaseChangeStabChain[chain,base,PermLength@First[chain],options];
BaseChangeStabChain[{StrongGenSet[{},GenSet[]]},base_,len_,options:OptionsPattern[]]:=StrongGenSet[#,GenSet[]]&/@(Drop[base,#]&/@Range[0,Length[base]]);
BaseChangeStabChain[chain_List,newbase_List,len_,options:OptionsPattern[]]:=If[OptionValue[MathLink],MathLinkBaseChangeStabilizerChain[First@chain,newbase,len],BaseChangeStabChain2[chain,newbase,len]];

BaseChangeStabChain2[chain_List,newbase_List,len_]:=Module[{i=0,j,g=ID,more,B=First@First@chain,newchain=chain,pos,sorbits={}},
more=If[Length[newbase]>0,
MemberQ[BasicOrbitStabChain[1,0,chain],newbase[[1]]],
False];
If[more,sorbits=SchreierOrbits[Last@First@chain,len]];
While[more,
i=i+1;
g=PermProduct[TraceSchreier[OnPoints[newbase[[i]],InversePerm[g]],sorbits],g];
more=(i+1)≤Min[Length[B],Length[newbase]];
If[more,more=MemberQ[OnPoints[newbase[[i+1]],InversePerm[g]],BasicOrbitStabChain[i+1,i,chain]]]
];
If[Not@PermEqual[g,ID],
newchain=ConjugateStabChain[chain,g];
B=First@First@newchain];
For[j=i+1,j≤Length[newbase],j++,
If[MemberQ[B,newbase[[j]]],pos=Position[B,newbase[[j]]][[1,1]],AppendTo[B,newbase[[j]]];newchain=AppendBasePointStabChain[newchain,newbase[[j]]];pos=Length[B]];
While[pos≠j &&pos>1,newchain=InterchangeStabChain[newchain,pos-1,len];B=First@First@newchain;pos=pos-1];
];
newchain
];

Instead of using the slow function DeleteRedundantGenerators, we use a quicker version, which only removes the redundant generators that are easy to find.

DeleteSomeRedundantGenerators[StrongGenSet[base_,GS_GenSet]]:=Module[{Si,Sip1,i,T=GenSet[],orbit,gens,check,toadd},
Sip1=GenSet[];
check[i_,gen_,t_]:=FreeQ[Orbit[base[[i]],t],OnPoints[base[[i]],gen]];
toadd[i_,gs_,T_]:=Module[{TT=T},If[check[i,#,TT],AppendTo[TT,#]]&/@gs;Complement[TT,T]];
For[i=Length[base],i≥1,i--,
Si=Stabilizer[Take[base,i-1],GS];
gens=Complement[Si,Sip1];
T=Join[T,toadd[i,gens,T]];
Sip1=Si;
];
StrongGenSet[base,T]
];

2.3. New code for SetStabilizer

This is not used at the moment, but it might be needed for larger examples later.

SetStabilizerProperty[points_List,perm_]:=Complement[OnPoints[points,perm],points]==={};

Options[Search2]={xPermVerbose->False};

SearchSetStabilizer[chain_List,s_Integer, l_Integer,stabset_,SGSK_,options___]:=Module[{base=chain[[1,1]],newSGSK=SGSK,k=Length[chain]-1,Deltas,gammas,Korbit,verb},
verb=xPermVerbose/.CheckOptions[options]/.Options[Search];
If[s==l,
newSGSK=chain[[l]],
newSGSK=SearchSetStabilizer[chain,s+1,l,stabset,SGSK,options];
Deltas=Orbit[base[[s]],chain[[s,2]]];
If[verb,Print["Branching over points ",Deltas]];
(* Avoid rechecking permutations. Not in Butler's algorithm *)
(*If[s=!=k,Deltas=Drop[Deltas,1]];*)
Do[
gammas=Deltas[[j]];
Korbit=Orbit[gammas,newSGSK[[2]]];
If[gammas===MinB[Korbit,base],
newSGSK=GenerateSetStabilizer[chain,s,l, s+1,Append[Take[base,s-1],gammas], stabset,newSGSK,verb];
If[Head[newSGSK]===Times,newSGSK=-newSGSK]]
,{j,1,Length[Deltas]}];
];
newSGSK
];

GenerateSetStabilizer[chain_List,s_,l_,i_,list_,stabset_,SGSK_,verb_]:=Module[{base=chain[[1,1]],k=Length[chain]-1,g,newSGSK=SGSK,Deltag,gammai,otherSGSK},
(*If[verb,Print["GenerateSetStabilizer[",chain,",",s,",",l,",",i,",",list,",",stabset,",",SGSK,",",verb,"]"]];*)
g=FromBaseImage[list,ReplacePart[chain[[1]],Take[base,Length[list]],1]];
If[i==l,
If[verb,Print["Generate at level ",s," with i=",i,". We have list ",list," and permutation ",g]];
If[Not@PermEqual[g,ID]&&SetStabilizerProperty[stabset,g],
Module[{},newSGSK=-StrongGenSet[base,Append[Last[SGSK],g]]; (* Is this a good base? *)
If[verb,Print["  Added permutation ",g,", newSGSK: ",-newSGSK, ", SGSK: ",SGSK]]]],
Deltag=Intersection[OnPoints[Orbit[base[[i]],chain[[i,2]]],g],stabset];
If[verb,Print["Generating over points ",Deltag]];
Do[
gammai=Deltag[[j]];
newSGSK=GenerateSetStabilizer[chain,s,l,i+1,Append[list,gammai],stabset,newSGSK,verb];
If[Head[newSGSK]===Times,Break[]],
{j,1,Length[Deltag]}];
];
newSGSK
];

SetStabilizer2[points_,sgs:StrongGenSet[base_,GS_],options___]:=
Module[{newbasesgschain=StabilizerChain2@BaseChange[sgs,points]},
SearchSetStabilizer[newbasesgschain,1, Length[points]+1, points,newbasesgschain[[Length[points]+1]],options]];

2.4. Smallest image of a set

This is an algorithm for computing the smallest image of a set Delta when a group G is acting on it. 
The algorithm can be found in the paper: Linton, S., “Finding the smallest image of a set”, ISSAC ‘04, 229-234 (2004) http://dx.doi.org/10.1145/1005285.1005319

FindSmallestImage[Delta_List,G_StrongGenSet,n_]:=Module[{MinOrbit,FindFirstAndAppend,MapPoints,GpChain=StabilizerChain2@G,Cands={{{},Delta}},Pass,Mp={},k=Length@Delta,m,i,MRs,FirstSchreierOrbit,xs,gs,DRjs,DRgs},
MinOrbit[expr_]:=Min@Orbit[#,First@GpChain]&/@expr;
FindFirstAndAppend[elements_,cand_]:=Append[cand[[1]],First@First@Position[cand[[2]],#,1,1]]&/@elements;
MapPoints[DR_,gs_]:=OnPoints[DR,#]&/@gs;
For[i=1,i≤k,i+=1,
m=n+1;
MRs=Min[MinOrbit[#]]&/@(Part[#[[2]],Complement[Range[1,k],#[[1]]]]&/@Cands);
While[Length@MRs>0,
If[First@MRs<m,
  {Pass={First@Cands},m=First@MRs},
If[First@MRs==m,Pass=Append[Pass,First@Cands]]];
MRs=Rest@MRs;
Cands=Rest@Cands;];
GpChain=BaseChangeStabChain[GpChain,{m}];
FirstSchreierOrbit=SchreierOrbit[m,First@GpChain,n];
xs=Intersection[First@FirstSchreierOrbit,#[[2]]]&/@Pass;
gs=InversePerm@TraceSchreier[#,FirstSchreierOrbit]&/@#&/@xs;
DRjs=MapThread[FindFirstAndAppend,{xs,Pass}];
DRgs=MapThread[MapPoints[#1[[2]],#2]&,{Pass,gs}];
Cands=MapThread[List,{Flatten[DRjs,1],Flatten[DRgs,1]}];
GpChain=Rest@GpChain;
Mp=Append[Mp,m];
];
Mp]

The same thing but keep track of the permutation that maps to the smallest image with slightly streamlined programming. Can we use a better base change algorithm here?

FindSmallestImageAndMapping[Delta_List, G_StrongGenSet, n_] := Module[{MinOrbit, FindFirstAndAppend, MapPoints, PermProducts, Gp = G, Cands = {{{}, Delta, ID}}, Pass, Mp = {}, k = Length@Delta, m, i, MRs, FirstSchreierOrbit, gs, CandFromPass},
MinOrbit[expr_] := Min@Orbit[#, Gp] & /@ expr;CandFromPass[{JR_, DR_, gperm_}] :=With[{xs = Intersection[First@FirstSchreierOrbit, DR]},gs = InversePerm@TraceSchreier[#, FirstSchreierOrbit] & /@ xs;Sequence @@ MapThread[List, {Append[JR, First@First@Position[DR, #, 1, 1]] & /@ xs, OnPoints[DR, #] & /@ gs, PermProduct[gperm, #] & /@ gs}]
];
For[i = 1, i <= k, i += 1,
m = n + 1;MRs = Min[MinOrbit[#]] & /@ (Part[#[[2]], Complement[Range[1, k], #[[1]]]] & /@ Cands);While[Length@MRs > 0,
If[First@MRs < m,
Pass = {First@Cands}; m = First@MRs,If[First@MRs == m, Pass = Append[Pass, First@Cands]]];
MRs = Rest@MRs;
Cands = Rest@Cands;];
Gp = BaseChange2[Gp, {m}];
FirstSchreierOrbit = SchreierOrbit[m, Gp, n];
Cands = CandFromPass /@ Pass;Gp = StrongGenSet[Rest@First@Gp, Stabilizer[{m}, Gp[[2]]]];
Mp = Append[Mp, m];
];
{Mp, Cands[[1, 3]]}]

FindSmallestImageAndMapping2[Delta_List, G_StrongGenSet, n_] := Module[{MinOrbit, FindFirstAndAppend, MapPoints, PermProducts, GpChain = StabilizerChain2@G, Cands = {{{}, Delta, ID}}, Pass, Mp = {}, k = Length@Delta, m, i, MRs, FirstSchreierOrbit, gs, CandFromPass},
      MinOrbit[expr_] := Min@Orbit[#, First@GpChain] & /@ expr;
      CandFromPass[{JR_, DR_, gperm_}] := With[{xs = Intersection[First@FirstSchreierOrbit, DR]},
            gs = InversePerm@TraceSchreier[#, FirstSchreierOrbit] & /@ xs;
Sequence @@ MapThread[List, {Append[JR, First@First@Position[DR, #, 1, 1]] & /@ xs,OnPoints[DR, #] & /@ gs, PermProduct[gperm, #] & /@ gs}]
];
For[i = 1, i <= k, i += 1,
m = n + 1;MRs = Min[MinOrbit[#]] & /@ (Part[#[[2]], Complement[Range[1, k], #[[1]]]] & /@ Cands);While[Length@MRs > 0,
If[First@MRs < m,
Pass = {First@Cands}; 
m = First@MRs,
If[First@MRs == m, Pass = Append[Pass, First@Cands]]];
MRs = Rest@MRs;
Cands = Rest@Cands;];
GpChain = BaseChangeStabChain[GpChain, {m}];
        FirstSchreierOrbit = SchreierOrbit[m, First@GpChain, n];
        Cands = CandFromPass /@ Pass;
        GpChain = Rest@GpChain;
        Mp = Append[Mp, m];
];
{Mp, Cands[[1, 3]]}]

2.5. Comparing groups

This section contains algorithms for comparing groups. In theory we compare two groups by computing the lexicographically ordered lists of elements, and then compare these lists lexicographically. This is not practical so we use algorithms that make an equivalent comparison.
The algorithms can be found in the paper: Hulpke, A., Linton, S., “Total Ordering on Subgroups and Cosets”, ISSAC ‘03, 156-160 (2003) http://dx.doi.org/10.1145/860854.860892

SmallestGenSet[G_StrongGenSet]:=SmallestGenSet[SmallestBaseStabChain@G]

SmallestGenSet[SmallBaseStabChainG_List]:=Module[{n=PermLength[First@SmallBaseStabChainG],L={},i, Gi,betai,orbbetai,orb,t,x},
i=Length[SmallBaseStabChainG]-1;
While[i>0,
Gi=SmallBaseStabChainG[[i]];
betai=First@First@Gi;
orbbetai=SchreierOrbit[betai,Gi[[2]]];
orb=Complement[First@orbbetai,Orbit[betai,GenSet@@L]];
While[Length[orb]>0,
t=TraceSchreier[Min[orb],orbbetai];
x=SmallestCoset[t,SmallBaseStabChainG[[i+1;;Length[SmallBaseStabChainG]]]];
L=Append[L,x];
orb=Complement[First@orbbetai,Orbit[betai,GenSet@@L]];
];
i=i-1;];
L]

SmallestCoset[p_,G_StongGenSet]:=SmallestCoset[p,SmallestBaseStabChain@G]

SmallestCoset[p_,StabChainG_List]:=Module[{i=1,rep=p, Gi,betai,orb, img, mu, omega,t,x},
While[i< Length[StabChainG],
Gi=StabChainG[[i]];
betai=First@First@Gi;
(* Print["i, betai, rep: ",i,", ",betai,", ",rep];*)
orb=SchreierOrbit[betai,Gi[[2]]];
img=OnPoints[First@orb,rep];
mu=Min[img];
omega=OnPoints[mu,InversePerm[rep]];
t=TraceSchreier[omega,orb];
(* Print["img, mu, omega, omegatest, t: ",img, ", ",mu,", ",omega,", ",OnPoints[betai,t],", ",t ];*)
rep=PermProduct[t,rep];
i=i+1;];
rep]

LargestElementInGroup[G_StongGenSet]:=LargestElementInGroup[SmallestBaseStabChain@G]

LargestElementInGroup[StabChainG_List]:=Module[{i=1,rep=ID, Gi,betai,orb, img, mu, omega,t,x},
While[i< Length[StabChainG],
Gi=StabChainG[[i]];
betai=First@First@Gi;
orb=SchreierOrbit[betai,Gi[[2]]];
img=OnPoints[First@orb,rep];
mu=Max[img];
omega=OnPoints[mu,InversePerm[rep]];
t=TraceSchreier[omega,orb];
rep=PermProduct[t,rep];
i=i+1;];
rep]

RemoveSuperfluousBasePointInStabChain[StabChain_List]:=Module[{firstelem=First@StabChain,tempchain},
If[Length[StabChain]==1,
StabChain,
If[firstelem[[2]]=== StabChain[[2,2]],
RemoveSuperfluousBasePointInStabChain[Rest@StabChain],
tempchain=RemoveSuperfluousBasePointInStabChain[Rest@StabChain];
Prepend[tempchain,StrongGenSet[Prepend[tempchain[[1,1]],firstelem[[1,1]]],firstelem[[2]]]]]]]

SmallestBaseStabChain[G_StrongGenSet]:=RemoveSuperfluousBasePointInStabChain@BaseChangeStabChain[G,Range[PermLength[G]]]

GroupLessOrEqualQ[U_StrongGenSet,V_StrongGenSet]:=Module[{
SmallChainU=SmallestBaseStabChain@U,
SmallChainV=SmallestBaseStabChain@V,
SmallGenSetU,SmallGenSetV,lenU,lenV,MaxElem,n},
n=Max[PermLength[U],PermLength[V]];
SmallGenSetU=SmallestGenSet@SmallChainU;
SmallGenSetV=SmallestGenSet@SmallChainV;
lenU=Length[SmallGenSetU];
lenV=Length[SmallGenSetV];
Which[lenU<lenV,
MaxElem=LargestElementInGroup@SmallChainU;
SmallGenSetU=Join[SmallGenSetU,Table[MaxElem,{lenV-lenU}]],
lenU>lenV,
MaxElem=LargestElementInGroup@SmallChainV;
SmallGenSetV=Join[SmallGenSetV,Table[MaxElem,{lenU-lenV}]],
True,Null];
SmallGenSetU=xAct`xPerm`Private`FromSign[TranslatePerm[#,{Images,n}],n]&/@SmallGenSetU;
SmallGenSetV=xAct`xPerm`Private`FromSign[TranslatePerm[#,{Images,n}],n]&/@SmallGenSetV;
(*Print["SmallGenSetU: ",SmallGenSetU];
Print["SmallGenSetV: ",SmallGenSetV];*)
OrderedQ[{SmallGenSetU,SmallGenSetV}]]

SmallestGroup[GroupList_List]:=Module[{
SmallChains=SmallestBaseStabChain/@GroupList,
SmallGenSets,SetLenths,MaxLength,MaxElem,i=1,n, SmallestIndex,SmallestCanonicalGroup},
n=Max@@(PermLength/@GroupList);
SmallGenSets=SmallestGenSet/@SmallChains;
SetLenths=Length/@SmallGenSets;
MaxLength=Max@@SetLenths;
While[i≤Length@GroupList,
If[SetLenths[[i]]<MaxLength,
MaxElem=LargestElementInGroup[SmallChains[[i]]];
SmallGenSets[[i]]=Join[SmallGenSets[[i]],Table[MaxElem,{MaxLength-SetLenths[[i]]}]]];
i=i+1;];
SmallGenSets=Map[xAct`xPerm`Private`FromSign[TranslatePerm[#,{Images,n}],n]&,SmallGenSets,{2}];
SmallestIndex=First@Ordering[SmallGenSets,1];
SmallestCanonicalGroup=GenSet@@Take[SmallGenSets[[SmallestIndex]],SetLenths[[SmallestIndex]]];
(* TranslatePerm[#,Cycles] *)
SmallestCanonicalGroup=Map[xAct`xPerm`Private`ToSign[#,n]&,SmallestCanonicalGroup,{1}];
{SmallestIndex,StrongGenSet[ SmallChains[[SmallestIndex,1,1]],SmallestCanonicalGroup]}]

2.6. Group canonicalization (old version)

This is an old code that is not used at the moment.

We canonicalize the representation of a group by checking if it has been used before. If it has been used, use the representation in the list $UsedGroups. If not, add it to $UsedGroups. 
All groups seen by ToCanonicalSGS can be found in $UsedGroups, along with the label for this group.
Hence, if you want to give a group sym, (not in $UsedGroups) a specific label, call ToCanonicalSGS[sym, label].

$UsedGroups={};

ToCanonicalSGS[sym_StrongGenSet]:=ToCanonicalSGS[sym,Null]

ToCanonicalSGS[sym_StrongGenSet,label_]:=Module[{groupcount=1,found=False},
While[(groupcount<=Length[$UsedGroups])&&!found,
found=SubgroupQ[$UsedGroups[[groupcount,1]],sym]&&SubgroupQ[sym,$UsedGroups[[groupcount,1]]];
If[!found,groupcount=groupcount+1]];
If[!found,
If[label===Null,
$UsedGroups=Append[$UsedGroups,{sym,DefaultLabelSym[sym]}],
$UsedGroups=Append[$UsedGroups,{sym,label}]]];
$UsedGroups[[groupcount]]]

CanonicalizeGroupsInSym[expr_]:=expr/.SymH[headlist_,sym_,label_]:>Module[{canonicalsym,newlabel},
{canonicalsym,newlabel}=ToCanonicalSGS[sym,label];
SymH[headlist,canonicalsym,newlabel]]

2.7. Product of groups

Here follows an algorithm that checks if perm is such that perm is in H.G where G is given by GS1 and H is given GS2. Observe that we need the same base for the two strong generating sets.
The idea is to write perm=PermProduct[u2,newperm,u1], where u1 is in G, u2 in H and newperm don't move the first point in the base. 
This means First@base = OnPoints[First@base,PermProduct[InversePerm[u2], perm, InversePerm[u1]]]. But this is equivalent to 
First@base = OnPoints[First@base,PermProduct[u1, InversePerm[perm], u2]] = OnPoints[OnPoints[First@base, u1], PermProduct[InversePerm[perm]], u2]].
We can rewrite this as OnPoints[First@base,PermProduct[InversePerm[u2]], perm]] = OnPoints[First@base, u1].
Let orbit1 and orbit2 denote the orbits of the first base point in G and H respectively.
Define point1=OnPoints[First@base, u1], then point1 must be in orbit1, but also in OnPoints[orbit2,perm].

In reality we don't know u1 and u2, but we can compute them from a TraceSchreier operation using point.
Hence, we consider all possible choices of point1 in the intersection between orbit1 and OnPoints[orbit2,perm].

From a specific choice of point1, we compute u1, and u2, and newperm.
We can then iterate with the first point removed from the base.

The different choices of point1 gives a tree. We check the first possible choice first, if that turns out to not give a solution, we remove that choice, and try the next one.

Sifting[perm_,StrongGenSet[base_List,GS1_GenSet],StrongGenSet[base_List,GS2_GenSet],word1_List:{},word2_List:{}]:=
If[Length[base]===0,
{PermEqual[perm,ID],word1,word2},
Module[{
Sorbit1=SchreierOrbit[First[base],GS1,Max[PermDeg[perm],PermDeg[GS1],PermDeg[GS2]]],
Sorbit2=SchreierOrbit[First[base],GS2,Max[PermDeg[perm],PermDeg[GS1],PermDeg[GS2]]],
mappedorbit2,possiblepoints,
point1,u1,u2,newperm, foundID=False,newword1={},newword2={}},
mappedorbit2=OnPoints[First@Sorbit2,perm];
possiblepoints=Cases[First@Sorbit1,Alternatives@@mappedorbit2];
While[(Length[possiblepoints]>0) && Not[foundID],
point1=First@possiblepoints;
possiblepoints=Rest@possiblepoints;
u1=TraceSchreier[point1,Sorbit1];
u2=InversePerm@TraceSchreier[OnPoints[point1,InversePerm[perm]],Sorbit2];
newperm=PermProduct[InversePerm[u2],perm,InversePerm[u1]];
{foundID,newword1,newword2}=Sifting[newperm,StrongGenSet[Rest[base],Stabilizer[{First[base]},GS1]],StrongGenSet[Rest[base],Stabilizer[{First[base]},GS2]],
Append[word1,u2],Prepend[word2,u1]];
];
{foundID,newword1,newword2}
]];

A somewhat streamlined version of the above algorithm: 
Here follows an algorithm that checks if perm is such that perm is in H.G where G is given by GS1 and H is given GS2. 
Observe that we need the same base for the two strong generating sets.
The idea is to write perm=PermProduct[InversePerm[u2], newperm, u1]
or equivalently newperm = PermProduct[u2, perm, InversePerm[u1]], where u1 is in G, u2 in H and newperm don't move the first point in the base. 
This means First@base = OnPoints[First@base,PermProduct[u2, perm, InversePerm[u1]]]. But this is equivalent to 
First@base = OnPoints[First@base,PermProduct[u1, InversePerm[perm], InvesePerm[u2]]] 
= OnPoints[OnPoints[First@base, u1], PermProduct[InversePerm[perm]], InvesePerm[u2]]]
= OnPoints[OnPoints[First@base, u1], InversePerm[PermProduct[u2,perm]]].
We can rewrite this as OnPoints[First@base,PermProduct[u2, perm]] = OnPoints[First@base, u1].
Let orbit1 and orbit2 denote the orbits of the first base point in G and H respectively.
Define point1=OnPoints[First@base, u1], then point1 must be in orbit1, but also in OnPoints[orbit2,perm].

In reality we don't know u1 and u2, but we can compute them from a TraceSchreier operation using point1.
Hence, we consider all possible choices of point1 in the intersection between orbit1 and OnPoints[orbit2,perm].

From all possible choices of point1, we compute u1, and u2, and newperm. 
All posibilities of newperm are computed, and sorted and duplicates removed. Call this list possiblenewperm.
We can then iterate with the first point removed from the base, and the first element in possiblenewperm as the permutation.
If the result is false, we try the next element in possiblenewperm.

Sifting2[perm_,base_List,GS1_GenSet,GS2_GenSet,n_]:=
If[Length[base]===0,
PermEqual[perm,ID],
Module[{
(* Compute the Schreierorbits of the first base point *)
Sorbit1=SchreierOrbit[First[base],GS1,n],
Sorbit2=SchreierOrbit[First[base],GS2,n],
mappedorbit2,possiblepoints1,
foundID=False, possibleu1,possibleu2,possiblenewperm},
(* Which points does perm map Sorbit2 to? *)
mappedorbit2=OnPoints[First@Sorbit2,perm];
(* point1 must be in orbit1 and in mappedorbit2. *)
possiblepoints1=Cases[First@Sorbit1,Alternatives@@mappedorbit2];
(* Compute the corresponding u1 and u2 for the points in possiblepoints1 *)
possibleu1=TraceSchreier[#,Sorbit1]&/@possiblepoints1;
possibleu2=TraceSchreier[OnPoints[#,InversePerm[perm]],Sorbit2]&/@possiblepoints1;
(* Compute all possible newperm, remove duplicates and sort the list, so we find ID as quickly as possible.*)
possiblenewperm=PermSort@Union@Inner[PermProduct[#2,perm,InversePerm[#1]]&,possibleu1,possibleu2,List];
While[(Length[possiblenewperm]>0) && Not[foundID],
foundID=Sifting2[First@possiblenewperm,Rest[base],Stabilizer[{First[base]},GS1],Stabilizer[{First[base]},GS2],n];
possiblenewperm=Rest@possiblenewperm;
];
foundID
]];

The same thing, but we take the the shreierorbits of the first points in the base for all groups in the stabilizer chain.

Sifting2b[perm_,base_List,StabChainSchOrbitsG1_List,StabChainSchOrbitsG2_List,n_]:=
If[Length[base]===0,
PermEqual[perm,ID],
Module[{
(* Compute the Schreierorbits of the first base point *)
Sorbit1=First@StabChainSchOrbitsG1,
Sorbit2=First@StabChainSchOrbitsG2,
mappedorbit2,possiblepoints1,
foundID=False, possibleu1,possibleu2,possiblenewperm},
(* Which points does perm map Sorbit2 to? *)
mappedorbit2=OnPoints[First@Sorbit2,perm];
(* point1 must be in orbit1 and in mappedorbit2. *)
possiblepoints1=Cases[First@Sorbit1,Alternatives@@mappedorbit2];
(* Compute the corresponding u1 and u2 for the points in possiblepoints1 *)
possibleu1=TraceSchreier[#,Sorbit1]&/@possiblepoints1;
possibleu2=TraceSchreier[OnPoints[#,InversePerm[perm]],Sorbit2]&/@possiblepoints1;
(* Compute all possible newperm, remove duplicates and sort the list, so we find ID as quickly as possible.*)
possiblenewperm=PermSort@Union@Inner[PermProduct[#2,perm,InversePerm[#1]]&,possibleu1,possibleu2,List];
While[(Length[possiblenewperm]>0) && Not[foundID],
foundID=Sifting2b[First@possiblenewperm,Rest[base],Rest@StabChainSchOrbitsG1,Rest@StabChainSchOrbitsG2,n];
possiblenewperm=Rest@possiblenewperm;
];
foundID
]];

The sifting would be much quicker if the points moved by the permutation are the first points in the base.
Sifting3 changes base so this is the case and then calles Sifting2. The base change can probably be done in a better and safer way. Can we guarantee that the bases of the groups are the same?

Sifting3[perm_,SGS1_StrongGenSet,SGS2_StrongGenSet]:=Sifting3[perm,SGS1,SGS2,Max[PermDeg@perm,PermLength@SGS1,PermLength@SGS2]]

Sifting3[perm_,SGS1_StrongGenSet,SGS2_StrongGenSet, len_]:=Module[{base=Join[First@SGS1,Complement[First@SGS2,First@SGS1]],movedpoints,GS1,GS2},
movedpoints:=Complement[Range[1,PermDeg@perm],StablePoints@perm];
base=Join[movedpoints,Complement[base,movedpoints]];
GS1=Last@BaseChange2[SGS1,base];
GS2=Last@BaseChange2[SGS2,base];
Sifting2[perm,base,GS1,GS2,len]
]

Check if - ID is in the product of the groups: (Improved 2014-09-14)

MinusIDInProductQ[H0_Symmetric,G0_Symmetric]:=False

MinusIDInProductQ[H0_Antisymmetric,G0_Antisymmetric]:=False

MinusIDInProductQ[H0_Symmetric,G0_Antisymmetric]:=Length@Intersection[First@H0,First@G0]≥2

MinusIDInProductQ[H0_Antisymmetric,G0_Symmetric]:=Length@Intersection[First@H0,First@G0]≥2

MinusIDInProductQ[H0_StrongGenSet,G0_Symmetric]:=Not@FreeQ[sign/@List@@Stabilizer[Complement[GroupSupport@H0,First@G0],H0[[2]]],-1]

MinusIDInProductQ[H0_Symmetric,G0_StrongGenSet]:=Not@FreeQ[sign/@List@@Stabilizer[Complement[GroupSupport@G0,First@H0],G0[[2]]],-1]

Old code for checking if -ID is in the product.

MinusIDInProductQold[H0_StrongGenSet,G0_StrongGenSet]:=
Module[{G1,H1},
(* Arrange so that the groups have the same bases and cycle notation *)
H1=BaseChange2[H0,First@G0];
G1=BaseChange2[G0,First@H1];
Sifting2[-ID,G1[[1]],G1[[2]],H1[[2]],Max[PermLength@G1,PermLength@H1]]
];

Slightly improved code. Now the Schreier vectors are only computed once. Changed 2013-10-09.

MinusIDInProductQ[H0_StrongGenSet,G0_StrongGenSet]:=
Module[{G1,H1,n},
(* Arrange so that the groups have the same base *)
H1=BaseChange2[H0,First@G0];
G1=ReplacePart[G0,1-> First@H1];
n=Max[PermLength@G1,PermLength@H1];
Sifting2b[-ID,G1[[1]],SchreierOrbit[First@First@#,#,n]&/@(Most@StabilizerChain2[G1]),SchreierOrbit[First@First@#,#,n]&/@(Most@StabilizerChain2[H1]),n]
];

MinusIDInProductQ[H0_,G0_]:=MinusIDInProductQ[xAct`xTensor`Private`SGSofsym@H0,xAct`xTensor`Private`SGSofsym@G0]

Old version of almost brute force search for the symmetry of the product of two groups.

GpinHGQBruteForceOld[p_,G_,H_,n_,Gelements_]:=
If[Sifting2[p,First@G,Last@G,Last@H,n],
Module[
{TestList=Gelements,
result=True},
While[(Length[TestList]>0)&&result,
result=And@@(Sifting2[PermProduct[First@TestList,p],First@G,Last@G,Last@H,n]);
TestList=Rest@TestList;];
result
],
False]

SymmetryOfProductOfGroupsBruteForceOld[G0_StrongGenSet,H0_StrongGenSet]:=
Module[{G1,H1,permdeg,K1,F1, GElements, unstablepoints,time},
(* Arrange so that the groups have the same bases and cycle notation *)
 (* Print["SymmetryOfProductOfGroupsBruteForce[",G0,",",H0,"]"];*)
If[Length@Last@H0===0,{G0,H0},
H1=BaseChange2[H0,First@G0];
G1=BaseChange2[G0,First@H1];
permdeg=Max[Max@@First[H1],PermDeg[Last@G1],PermDeg[Last@H1]];
(*Print["H1: ",H1];
Print["G1: ",G1];
Print["permdeg: ",permdeg];*)
(* GElements=Rest[List@@Dimino[G1]];*)
(* It should be enough to test one from each coset H.p where p is in G *)
{time,GElements}=AbsoluteTiming@Rest[Union[First[RightCosetRepresentative[#,permdeg,H1]]&/@List@@Dimino@G1]];
(* Print["GElements timing: ",time, " length: ", Length@GElements];*)
(*Print["GElements: ",GElements];*)
K1=Search[H1,GpinHGQBruteForce[#,G1,H1,permdeg,GElements]&,1,StrongGenSet[{},GenSet[]]];
(* Print["K1: ",K1];*)
F1=SchreierSims[Range@permdeg,Join[Last@K1,Last@G1]];
F1=DeleteSomeRedundantGenerators[F1];
unstablepoints=Complement[Range@permdeg,StablePoints[Last@F1,permdeg]];
F1=ReplacePart[F1,1->unstablepoints];
(*Print["F1: ",F1];*)
{F1,K1}]]

Newer version where the Schreier vectors are computed only once. Changed 2013-10-09. This is still not efficient, so a new algorithm is needed.

GpinHGQBruteForce[p_,base_,StabChainSchOrbG_,StabChainSchOrbH_,n_,Gelements_]:=
If[xAct`SymManipulator`Private`Sifting2b[p,base,StabChainSchOrbG,StabChainSchOrbH,n],
Module[
{TestList=Gelements,
result=True},
While[(Length[TestList]>0)&&result,
result=xAct`SymManipulator`Private`Sifting2b[PermProduct[First@TestList,p],base,StabChainSchOrbG,StabChainSchOrbH,n];
TestList=Rest@TestList;];
result
],
False]

SymmetryOfProductOfGroupsBruteForce[G0_StrongGenSet,H0_StrongGenSet]:=
Module[{G1,H1,permdeg,K1,F1, GElements, unstablepoints,time,StabChainSchOrbG, StabChainSchOrbH },
(* Arrange so that the groups have the same bases and cycle notation *)
 (* Print["SymmetryOfProductOfGroupsBruteForce[",G0,",",H0,"]"];*)
If[Length@Last@H0===0,{G0,H0},
H1=BaseChange2[H0,First@G0];
G1=ReplacePart[G0,1-> First@H1];
permdeg=Max[Max@@First[H1],PermDeg[Last@G1],PermDeg[Last@H1]];
StabChainSchOrbG=SchreierOrbit[First@First@#,#,permdeg]&/@(Most@StabilizerChain2[G1]);
StabChainSchOrbH=SchreierOrbit[First@First@#,#,permdeg]&/@(Most@StabilizerChain2[H1]);
(*Print["H1: ",H1];
Print["G1: ",G1];
Print["permdeg: ",permdeg];*)
(* GElements=Rest[List@@Dimino[G1]];*)
(* It should be enough to test one from each coset H.p where p is in G *)
{time,GElements}=AbsoluteTiming@Rest[Union[First[RightCosetRepresentative[#,permdeg,H1]]&/@List@@Dimino@G1]];
(* Print["GElements timing: ",time, " length: ", Length@GElements];*)
(*Print["GElements: ",GElements];*)
K1=Search[H1,GpinHGQBruteForce[#,First@G1,StabChainSchOrbG,StabChainSchOrbH,permdeg,GElements]&,1,StrongGenSet[{},GenSet[]]];
(* Print["K1: ",K1];*)
F1=SchreierSims[Range@permdeg,Join[Last@K1,Last@G1]];
F1=xAct`SymManipulator`Private`DeleteSomeRedundantGenerators[F1];
unstablepoints=Complement[Range@permdeg,StablePoints[Last@F1,permdeg]];
F1=ReplacePart[F1,1->unstablepoints];
(*Print["F1: ",F1];*)
{F1,K1}]]

SymmetryOfProductOfGroupsBruteForce[G0_,H0_]:=SymmetryOfProductOfGroupsBruteForce[xAct`xTensor`Private`SGSofsym@G0,xAct`xTensor`Private`SGSofsym@H0]

2.8. Transversals

SchreierTransversal[x_,base_]:=TraceSchreier[#,x]&/@(SortB[First@x,base])

SchreierTransversals[SGS_]:=
Module[{Gchain=StabilizerChain2@SGS,SchreierOrbits1},
SchreierOrbits1=SchreierOrbit[Gchain[[#,1,1]],Gchain[[#]],
PermLength@Gchain[[#]]]&/@Range[Length[Gchain]-1];
SchreierTransversal[#,First@SGS]&/@SchreierOrbits1
]

The following algorithm is taken from “Handbook of computational group theory” by Derek F. Holt page 128.
This implementation has to be verified (does it produce the correct thing even if the base is not sorted?):

TransversalComputation[H0_,G0_]:=
Module[{base=First@G0,baselength=Length@First@G0,StabChainG0=StabilizerChain2@G0,
StabChainH0=StabilizerChain2@H0,
STransG0=SchreierTransversals@G0,
TransversalLengths,RestBaseOrbitsH,TList,TBar,i,j,Checku,ElementsNeeded},
Checku[u_,n_]:=
Module[{tuList=PermProduct[#,u]&/@TList,FirstBaseImageUnderu=OnPoints[StabChainH0[[n,1,1]],u],
SelectionList},
(*Print["Checku[",u,",",n,"]"];*)
SelectionList=And@@@Outer[Not[xAct`xPerm`Private`LessB[{OnPoints[#2,#1],FirstBaseImageUnderu},base]]&,tuList,RestBaseOrbitsH[[n]]];
(*Print["SelectionList:",SelectionList];*)
Pick[tuList,SelectionList]
];
TransversalLengths=(OrderOfGroup[StabChainG0[[#]]]/OrderOfGroup[StabChainH0[[#]]])&/@Range[baselength+1];
RestBaseOrbitsH=Rest@Orbit[StabChainH0[[#,1,1]],StabChainH0[[#]]]&/@Range[baselength];
TList={ID};
i=baselength;
While[i≥1,
TBar={};
j=2;
ElementsNeeded=TransversalLengths[[i]]-Length[TList];
While[And[Length[TBar]<ElementsNeeded,
j≤Length[STransG0[[i]]]],
TBar=Join[TBar,Checku[STransG0[[i,j]],i]];
j=j+1;
];
TList=Join[TList,TBar];
i=i-1;];
TList
]

TransversalInSymmetricGroup[H0_StrongGenSet,G0_Symmetric]:=
Module[{suppH0=GroupSupport@H0,H0partofbase, newbase,G1},If[Not@subsetQ[suppH0,First@G0],Throw[Message[General::error,"TransversalInSymmetricGroup: Not a subgroup."]]];
(* Arrange so both groups have the same base. *)
(* Elliminate points in the base that is not in G0.*)
H0partofbase=DeleteCases[First@H0,Not@MemberQ[#,First@G0]&];
(* construct the new base. *)
newbase=Join[H0partofbase,Complement[First@G0,H0partofbase]];
(* Construct symmetric group with correct base. *)
G1=StrongGenSet[newbase,GenSet@@(Cycles/@Partition[DeleteCases[newbase,Not@MemberQ[#,First@G0]&],2,1])];
(* Compute the transversal. *)
TransversalComputation[ReplacePart[H0,1->newbase],G1]
]

TransversalInCompatibleSymmGroup[H0_StrongGenSet,G0_Symmetric]:=TransversalInSymmetricGroup[H0,G0]

TransversalInCompatibleSymmGroup[H0_StrongGenSet,G0_StrongGenSet]:=
Module[{suppH0=GroupSupport@H0,suppG0=GroupSupport@G0,H0partofbase, newbase,G1},If[Not@subsetQ[suppH0,suppG0],Throw[Message[General::error,"TransversalInCompatibleSymmGroup: Not a subgroup."]]];
(* Arrange so both groups have the same base. *)
(* Elliminate points in the base that is not in G0.*)
H0partofbase=DeleteCases[First@H0,Not@MemberQ[#,suppG0]&];
(* construct the new base. *)
newbase=Join[H0partofbase,Complement[suppG0,H0partofbase]];
(* Construct symmetric group with correct base. *)
G1=BaseChange2[G0,newbase];
(* Compute the transversal. *)
TransversalComputation[ReplacePart[H0,1->newbase],G1]
]

3. Definition of Sym and development of some tools

3.1. Definition of ZeroTensor and deltaH

xTensorQ[ZeroTensor[slotlist_]]^:=True;
PrintAs[ZeroTensor[slotlist_]]^:="0";
xAct`TexAct`Tex[ZeroTensor[slotlist_]]^:="0";
SlotsOfTensor[ZeroTensor[slotlist_]]^:=slotlist;
DependenciesOfTensor[ZeroTensor[slotlist_]]^:={};
DefInfo[ZeroTensor[slotlist_]]^:={"tensor",StringJoin["ZeroTensor with slots ",ToString/@slotlist," ."]};
HostsOf[ZeroTensor[slotlist_]]^:={};
TensorID[ZeroTensor[slotlist_]]^:={};
Dagger[ZeroTensor[slotlist_List]]^:=ZeroTensor[Dagger/@slotlist]

ZeroTensor[slotlist_][___]:=0

TensorToZeroRule[T_?xTensorQ]:=T->ZeroTensor[SlotsOfTensor@T]

xTensorQ[deltaH[vbundle_]]^:=True;
PrintAs[deltaH[vbundle_]]^:="δ";
SlotsOfTensor[deltaH[vbundle_]]^:={-vbundle,vbundle};
SymmetryGroupOfTensor[deltaH[vbundle_]]^=StrongGenSet[{},GenSet[]];
DependenciesOfTensor[deltaH[vbundle_]]^:={};
DefInfo[deltaH[vbundle_]]^:={"tensor",StringJoin[" delta on ",ToString/@vbundle," ."]};
HostsOf[deltaH[vbundle_]]^:={};
TensorID[deltaH[vbundle_]]^:={};

deltaHTodelta=deltaH[vbundle_]:>delta;

MetricOrdeltaHOfVBundle[vbundle_]:=
Module[{metrics=MetricsOfVBundle[vbundle]},
If[Length@metrics>0,First@metrics,deltaH[vbundle]]];

3.2. Definition of the SymH object

nosign[-x_]:=x;
nosign[x_]:=x;
sign[-x_]:=-1;
sign[x_]:=1;

slot=xAct`xTensor`Private`slot;

The SymH symbol is a place holder for ImposeSymmetry.

First general definition of the SymH symbol.

xTensorQ[SymH[headlist_,sym_,label_]]^:=True;
PrintAs[SymH[headlist_,sym_,label_]]^:=RowBox[{UnderscriptBox["Sym",StyleBox[label,Tiny]],"[",Sequence@@(PrintAs/@headlist),"]"}];
SlotsOfTensor[SymH[headlist_,sym_,label_]]^:=Join@@(SlotsOfTensor/@headlist);
DependenciesOfTensor[SymH[headlist_,sym_,label_]]^:=Union@@(DependenciesOfTensor/@headlist);
DefInfo[SymH[headlist_,sym_,label_]]^:={"tensor",StringJoin["Symmetrized ",Sequence@@(ToString/@headlist)," ."]};
(* JMM: These two should also be generalized, but are not relevant for us *)
HostsOf[SymH[headlist_,sym_,label_]]^:={};
TensorID[SymH[headlist_,sym_,label_]]^:={};

If the ZeroTensor appears in the headlist, the entire object is a ZeroTensor.

SymH[headlist:{___,ZeroTensor[_],___},sym_,label_]:=ZeroTensor[Join@@(SlotsOfTensor/@headlist)]

3.3. Expanding the symmetry

Extracting the internal expression. Observe that CovarD is not expanded.

InternalExprSym[SymH[headlist_,sym_,label_]]:=
Module[{numindices=Length/@SlotsOfTensor/@headlist,slotlist,partitionedslots,internalexpr},
slotlist=IndexList@@slot/@Range[Plus@@numindices];
(* Extract indices belonging to the different tensors. Can this be done in a simpler way? *)
(* partitionedslots=IndexList@@@Last/@Rest@FoldList[{Drop[#1[[1]],#2],Take[#1[[1]],#2]}&,{List@@slotlist,{}},numindices];*)
partitionedslots=IndexList@@@partitionRagged[List@@slotlist,numindices];
internalexpr=Inner[(#1@@#2)&,headlist,partitionedslots,Times]]

Extracting the internal symmetry.
Unfortunately you can not take SymmetryOf@InternalExprSym[x], because Identify requires specification of VBundle. Therefore, we introduce temporary indices.

InternalSymmetrySym[x:(SymH[headlist_,sym_,label_])]:=
Module[{internalexpr=InternalExprSym[x],tempinds,internalsym,slotrules,inverseslotrules, n,tmpinternal,internalslotrules, slotpermutations},
tempinds=DummyIn/@SlotsOfTensor[x];
slotrules=Thread@Rule[slot/@Range[Length@tempinds],tempinds];
inverseslotrules=Thread@Rule[tempinds,slot/@Range[Length@tempinds]];
{n,tmpinternal,internalslotrules,internalsym}=List@@SymmetryOf[internalexpr/.slotrules];
slotpermutations=(internalslotrules/.inverseslotrules);
internalsym/.(slotpermutations/.slot[nn_]:>nn)
]

Expanding the imposed symmetry.

ExpandSym1[x:(SymH[headlist_,sym_,label_][inds___])]:=Module[{indlist=List[inds],slotlist=slot/@Range@Length@List[inds],internalexpr, slotexpansion},
internalexpr=InternalExprSym[Head@x];
slotexpansion=ImposeSymmetry[internalexpr,IndexList@@slotlist,sym];
slotexpansion/.Thread@Rule[slotlist,indlist]]

The following expands the imposed symmetry, but the internal symmetry is used to reduce the number of terms. 
Instead of expanding the imposed symmetry, the complete mixed symmetry is used instead. The terms is then taken as a transversal of InternalCommutingSymmetry as a subgroup of this mixed symmetry.

ExpandSym2[x:(SymH[headlist_,sym_,label_][inds___])]:=Module[{indlist=List[inds],slotlist=slot/@Range@Length@List[inds],internalexpr, slotexpansion, bigsym,smallsym,transversal},
internalexpr=xAct`SymManipulator`Private`InternalExprSym[Head@x];
bigsym=SymmetryGroupOfTensor[Head@x];
smallsym=ReplacePart[BaseChange2[xAct`xTensor`Private`SGSofsym@InternalCommutingSymmetry[Head@x],First@bigsym],1->First@bigsym];
(* If[Not@SubgroupQ[smallsym,bigsym],Print["Not a subgroup..."]]; *)
transversal=InversePerm/@TransversalComputation[smallsym,bigsym];
slotexpansion=Plus@@(PermuteIndices[internalexpr,IndexList@@slotlist,#]&/@transversal)/Length[transversal];
slotexpansion/.Thread@Rule[slotlist,indlist]]

Removing the imposed symmetry. This is the same thing as InternalExprSym, but with the correct indices, not some dummy indices. Observe that this can be dangerous if there are contractions on any of the slots in the symmetry. The contractions should be taken after the symmetry is imposed.

RemoveSym1[x:(SymH[headlist_,sym_,label_][inds___])]:=Module[{indlist=List[inds],slotlist=IndexList@@slot/@Range@Length@List[inds]},
InternalExprSym[Head@x]/.Thread@Rule[List@@slotlist,indlist]]

Options[ExpandSym]={SmartExpand->False};

User versions of the functions above. These work on any expression, it does not have to be just a SymH object. It also expands CovarD.

ExpandSym[expr_,options:OptionsPattern[]]:=ExpandCovarDAndInertHeadHead[
If[OptionValue[SmartExpand],
expr/.x:(SymH[___][inds___]):>ExpandSym2[x],
expr/.x:(SymH[___][inds___]):>ExpandSym1[x]]]

ExpandSymS[LL_]:=ExpandSym[LL,SmartExpand->True]

RemoveSym[expr_]:=ExpandCovarDAndInertHeadHead[expr/.x:(SymH[___][inds___]):>RemoveSym1[x]]

3.4. Labels

Here are some functions for generating a label that can sometimes describe the symmetry. If the group is not recognized, the unstable slots are listed as the description. One can recognize a disjoint union of symmetric and antisymmetric groups by studying the group orbit by orbit. This is probably not true for other groups.

DefaultOrbitLabel[points_,sym_]:=
Module[{orbitlabel=StringJoin[ToString/@points]},
Which[SubgroupQ[Symmetric[points],sym],(* Symmetric *)
StringJoin["(",orbitlabel,")"],
SubgroupQ[Antisymmetric[points],sym],(* Antisymmetric *)
StringJoin["[",orbitlabel,"]"],
True,(* Everything else *)
orbitlabel]]

DefaultLabelSym[sym_StrongGenSet]:=Module[{n=PermLength@Last[sym], longorbits,unstableslots, orbitlabels},
longorbits=Sort/@Select[Orbits[sym,n],Length[#]>1&];
orbitlabels=DefaultOrbitLabel[#,sym]&/@longorbits;
StringJoin[ToString/@orbitlabels]]

DefaultLabelSym[sym_Symmetric]:=StringJoin["(",StringJoin[ToString/@First@sym],")"]

DefaultLabelSym[sym_Antisymmetric]:=StringJoin["[",StringJoin[ToString/@First@sym],"]"]

3.5. Computation of the total symmetry of the symmetrized object

This computation is only works if the imposed symmetry is Symmetric[ _ ] or Antisymmetric[ _ ].
In all other cases it returns the imposed symmetry only.
Observe that we compute the symmetry only once.

$MixedSymSpecialCode=True;
$MixedSymVerbose=False;

SymmetryGroupOfTensor[x : SymH[headlist_, sym_, label_]] ^:=SymmetryGroupOfTensor[x]^=Module[{n=Plus@@(Length/@SlotsOfTensor/@headlist),unstableslots,newlabel,symmetrygroup, internalcommutingsym},
unstableslots=GroupSupport@sym;
(* In the special cases when the symmetry is 
Symmetric[unstableslots], or Antysymmetric[unstableslots], we compute the mixed symmetry. *)
{symmetrygroup, internalcommutingsym}=
Which[And[$MixedSymSpecialCode===True,Head@sym===Symmetric],(* Symmetric *)
SymmetryOfProductOfGroupsSymmetricCase2b[unstableslots,InternalSymmetrySym[SymH[headlist, sym, label]],False],
And[$MixedSymSpecialCode===SetStabilizer,Head@sym===Symmetric],(* Symmetric *)
SymmetricMixedSymmetry[unstableslots,InternalSymmetrySym[SymH[headlist, sym, label]],sym,n],
And[$MixedSymSpecialCode===SetStabilizer,Head@sym===Antisymmetric],(* Antisymmetric *)
AntisymmetricMixedSymmetry[unstableslots,InternalSymmetrySym[SymH[headlist, sym, label]],sym,n],
True,(* Everything else *)
GeneralMixedSymmetry[unstableslots,InternalSymmetrySym[SymH[headlist, sym, label]],sym,n]];
InternalCommutingSymmetry[x]^= internalcommutingsym;
symmetrygroup]

ClearRememberedSymHSymmetry[tensor_]:=UpValues[SymH]=DeleteCases[UpValues[SymH],_?((#[[1,0]]==HoldPattern)&&Or[#[[1,1,0]]==SymmetryGroupOfTensor,#[[1,1,0]]==InternalCommutingSymmetry]&&(#[[1,1, 1,0]]==SymH)&&MemberQ[#[[1,1,1,1]],tensor,Infinity]&)];

3.6. Changing UndefTensor to clear SymmetryGroupOfTensor results

xTension["SymManipulator`",UndefTensor,"Beginning"]=ClearRememberedSymHSymmetry;

3.7. InertHeadHead and CovarD

Unprotect[InertHeadHead];

xTensorQ[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ]]^:=True;
PrintAs[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ]]^:=RowBox[{PrintAs[H1],PrintAs[T]}];
SlotsOfTensor[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ]]^:=SlotsOfTensor[T];
DependenciesOfTensor[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ]]^:=DependenciesOfTensor[T];
SymmetryGroupOfTensor[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ]]^:=SymmetryGroupOfTensor[T];
DefInfo[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ]]^:={"tensor",StringJoin[ToString[H1],ToString[T]]};
HostsOf[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ]]^:={};
TensorID[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ]]^:={};

Unprotect[CovarD];

xTensorQ[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]^:=True;
PrintAs[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]^:=RowBox[{Last@SymbolOfCovD[D1],PrintAs[T]}];
SlotsOfTensor[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]^:=Join[list,SlotsOfTensor[T]];
(* Is this correct? *)
DependenciesOfTensor[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]^:=DependenciesOfTensor[T];
SymmetryGroupOfTensor[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]^:=
JoinSGS[SymmetryGroupOfCovD[D1],SymmetryGroupOfTensor[T]/.Thread@Rule[Range@Length@SlotsOfTensor[T],Range[1+Length[list],Length[list]+Length@SlotsOfTensor[T]]]];
DefInfo[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]^:={"tensor",StringJoin[ToString[D1],ToString[T]]};
(* These two should also be generalized... *)
HostsOf[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]^:={};
TensorID[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]^:={};

CovarD[D1_?CovDQ,ZeroTensor[slotlist_],list_]:=ZeroTensor[Join[list,slotlist]]

InertHeadHead[H1_?InertHeadQ,ZeroTensor[slotlist_]]:=ZeroTensor[slotlist]

ExpandInertHeadHead[InertHeadHead[H1_?InertHeadQ,T_?xTensorQ][inds___]]:=H1[T[inds]]

ExpandCovarD[CovarD[D1_?CovDQ,T_?xTensorQ,list_][inds___]]:=D1[Sequence@@Take[List[inds],Length@list]][T[Sequence@@Drop[List[inds],Length@list]]]

ExpandCovarDAndInertHeadHead[expr_]:=expr//.{x:(InertHeadHead[___][inds___]):>ExpandInertHeadHead[x],(x:CovarD[___][inds___]):>ExpandCovarD[x]}

ToInertHeadHead[H1_?InertHeadQ[T_?xTensorQ[inds___]]]:=InertHeadHead[H1,T][inds]

ToInertHeadHead[H1_?InertHeadQ[H2_?InertHeadQ[expr_]]]:=ToInertHeadHead[H1@ToInertHeadHead[H2[expr]]]

ToInertHeadHead[H1_?InertHeadQ[D1_?CovDQ[covdinds___][expr_]]]:=ToInertHeadHead[H1@ToCovarD[D1[covdinds][expr]]]

Always use DownIndex as standard position. It is important that we do not have any freedom in the choice. Otherwise we need to deal with this ambiguity during canonicalization.

ToCovarD[D1_?CovDQ[covdinds___][T_?xTensorQ[tensorinds___]]]:=
CovarD[D1,T,(-VBundleOfIndex[#])&/@List[covdinds]][Sequence@@Join[List[covdinds],List[tensorinds]]]

ToCovarD[D1_?CovDQ[covdinds___][D2_?CovDQ[covdinds2___][expr_]]]:=ToCovarD[D1[covdinds]@ToCovarD[D2[covdinds2][expr]]]

ToCovarD[D1_?CovDQ[covdinds___][H1_?InertHeadQ[expr_]]]:=ToCovarD[D1[covdinds]@ToInertHeadHead[H1[expr]]]

3.8. Dagger

Unprotect[Dagger];

This overloads the dagger symbol for ordinary tensors, so we have to remove that first, and add it later.

Dagger[xAct`xTenosr`Private`tensor_?xTensorQ[xAct`xTenosr`Private`inds___]] =.

To avoid problems with contracted indices, we replace the indices with temporary ones.

Dagger[x:(SymH[headlist_,sym_,label_]?xTensorQ[inds___])]:=Module[{tempinds=DummyAs/@IndexList[inds]},
ImposeSym[Dagger@RemoveSym1[Apply[Head@x,tempinds]],DaggerIndex/@tempinds,sym]/.Thread@Rule[DaggerIndex/@(List@@tempinds),DaggerIndex/@List[inds]]]

Dagger[tensor_?xTensorQ[inds___]] := Apply[Dagger[tensor], DaggerIndex@IndexList[inds]];

This might not be correct for all cases

Dagger[CovarD[D1_?CovDQ,T_?xTensorQ,list_]]:=CovarD[D1,Dagger[T],Dagger/@list]

Protect[Dagger];

3.9. Imposing the symmetry on expressions

First some tools

SeparateVBundles[inds_IndexList]:=Module[{VBundleOfIndexQ,VBlist,sepatatedinds},
(* Temporary function for checking if an index belongs to VBundle1 *)
VBundleOfIndexQ[VBundle1_]:=VBundle1==VBundleOfIndex[#]&;
(* List the VBundles in inds *)
VBlist=List@@Union@Sort[VBundleOfIndex/@inds];
(* Separate the indices according to VBundles *)
sepatatedinds=Select[List@@inds,VBundleOfIndexQ[#]]&/@VBlist;
(* Find the corresponding position in the list inds and retun that.*)
First@First@Position[inds,#,1]&/@#&/@sepatatedinds
]

The following code computes the strong generating set for a totally symmetric tensor/spinor. The symmetries are taken over the different VBundles separately.
This is important if the spinor has both primed and unprimed indices.

CompatibleSymmetric[inds_IndexList]:=Module[{VBundles,SGSlist},
(* Compute the strong generating sets for the separate symmetries *)
VBundles=Select[SeparateVBundles[inds],(Length[#]≥2)&];
If[Length@VBundles==1,SGSlist=Symmetric/@VBundles,
SGSlist=Symmetric[#,Cycles]&/@VBundles;
(* Join the symmetries *)
While[Length[SGSlist]≥2,
SGSlist[[2]]=JoinSGS[SGSlist[[1]],SGSlist[[2]]];
SGSlist=Delete[SGSlist,1];
];];
(* Return the remaining SGS or the trivial SGS if the list is empty. *)
If[Length[SGSlist]==1,First@SGSlist,StrongGenSet[{}, GenSet[]]]
]

The same thing but just give it the VBundles. Perhaps the code above can be replaced by this.

CompatibleSym[vbundles_List]:=With[{sgs=JoinSGS@@(Symmetric[#,Cycles]&/@GatherBy[Range@Length@vbundles,(UpIndex/@vbundles)[[#]]&])},
If[Last@sgs===GenSet[],
StrongGenSet[{},GenSet[]],
sgs]];

A second argument can be given to specify which VBundles we want to symmetrize over.

CompatibleSym[vbundles_List,selectionVBs_List]:=With[{sgs=JoinSGS@@(Symmetric[#,Cycles]&/@(Select[GatherBy[Range@Length@vbundles,(UpIndex/@vbundles)[[#]]&],MemberQ[selectionVBs,UpIndex@vbundles[[First[#]]]]&]))},
If[Last@sgs===GenSet[],
StrongGenSet[{},GenSet[]],
sgs]];

Check if sym is CompatibleSym.

CompatibleSymQ[{},sym_]:=OrderOfGroup@sym==1

CompatibleSymQ[vbundles_List,sym_]:=Module[{n=Length@vbundles,upvbundles=UpIndex/@vbundles,disjsubgroups},disjsubgroups=Symmetric/@GatherBy[Range@n,upvbundles[[#]]&];
SubgroupQ[Join@@(GenSet/@disjsubgroups),sym]]

A third argument can be given to specify which VBundles we want to symmetrize over.

CompatibleSymQ[vbundles_List,sym_,selectionVBs_List]:=Module[{n=Length@vbundles,disjsubgroups},disjsubgroups=(Symmetric[#]&/@(Select[GatherBy[Range@n,(UpIndex/@vbundles)[[#]]&],MemberQ[selectionVBs,UpIndex@vbundles[[First[#]]]]&]));
If[Length@disjsubgroups==0,True,
SubgroupQ[Join@@(GenSet/@disjsubgroups),sym]]]

Refer back to the previous code for tensors or product of tensors

TimesToList[expr_Times]:=List@@expr

TimesToList[xAct`xTensor`Private`VerbatimProduct[Times][expr__]]:=(List[expr]/.xAct`xTensor`Private`Scalar1[Scalarexpr_]:>Scalarexpr)

TimesToList[T_?xTensorQ[tensorinds__]]:=List[T[tensorinds]]

Main code

Options[ImposeSym]={ImposeLargerSym->False,ImposeSuperfluousSym->False};

Shortcuts

ImposeSym[0]:=0
ImposeSym[expr_,options___?OptionQ]:=ImposeSym[expr,FindFreeIndices[expr],options]

ImposeSym[expr_,inds_,options___?OptionQ]:=ImposeSym[expr,inds,CompatibleSymmetric[inds],options]

ImposeSym[expr_,imposedinds_,sym_GenSet,options___?OptionQ]:=Module[{n=PermLength[sym],unstableslots,SGS},
unstableslots=Complement[Range@n,StablePoints[sym,n]];
SGS=SchreierSims[unstableslots,sym];
ImposeSym[expr,imposedinds,SGS,options]]

ImposeSym[expr_,IndexList[___],StrongGenSet[{___},GenSet[]],other___]:=expr

Step 0: Expand input.

ImposeSym[expr_,inds_,sym_,options___?OptionQ]:=ImposeSym1[xAct`xTensor`Private`MathInputExpand@expr,inds,sym,options]

Step 1: Reduce to a product of tensors

ImposeSym1[expr_Plus,other__]:=ImposeSym1[#,other]&/@expr
ImposeSym1[0,other__]:=0

ImposeSym1[x_?ConstantQ y_,other__]:=x ImposeSym1[y,other]
ImposeSym1[Scalar[x_] y_,other__]:=Scalar[x] ImposeSym1[y,other]
ImposeSym1[(F_?ScalarFunctionQ)[Scalar[x_]]y_,other__]:=F[Scalar[x]] ImposeSym1[y,other]
ImposeSym1[(F_?ScalarFunctionQ)[x_?xTensorQ[]]y_,other__]:=F[x[]] ImposeSym1[y,other]
ImposeSym1[x_?ScalarQ y_,other__]:=x ImposeSym1[y,other]
ImposeSym1[Times[x___,y_?xTensorQ[],z___],other__]:=y[]ImposeSym1[Times[x,z],other]
ImposeSym1[Times[x___,Power[y_?xTensorQ[],n_],z___],other__]:=Power[y[],n]ImposeSym1[Times[x,z],other]
ImposeSym1[Times[x___,Power[Scalar[y_],n_],z___],other__]:=Power[Scalar[y],n]ImposeSym1[Times[x,z],other]
ImposeSym1[Times[x___,Power[y_?ScalarQ ,n_],z___],other__]:=Power[y,n]ImposeSym1[Times[x,z],other]

ImposeSym1[z:(x_?InertHeadQ[___]),other__]:=ImposeSym1[ToInertHeadHead[z],other]

ImposeSym1[z:(x_?CovDQ[__][y_]),other__]:=ImposeSym1[ToCovarD[z],other]

ImposeSym1[LieD[dir_][expr_],other___]:=LieD[dir]@ImposeSym1[expr,other]

ImposeSym1[Times[z___,x:(_?InertHeadQ[___]),y___],other__]:=ImposeSym1[Times[z,ToInertHeadHead[x],y],other]

ImposeSym1[Times[z___,x:(_?CovDQ[___][___]),y___],other__]:=ImposeSym1[Times[z,ToCovarD[x],y],other]

ImposeSym1[expr:Times[_?xTensorQ[__],_?xTensorQ[__]...],other__]:=ImposeSym2[expr,other]

ImposeSym1[T_?xTensorQ[tensorinds__],other__]:=ImposeSym2[T[tensorinds],other]

ImposeSym1[other__]:=Throw@Message[ImposeSym::error,"Unable to impose symmetry: ImposeSym[",StringJoin@@Riffle[ToString/@{other},","],"]"]

Step 2: Impose symmetry on product of tensors.

ImposeSym2[expr_,imposedinds_,sym:(_StrongGenSet|_Symmetric|_Antisymmetric),options___?OptionQ]:=
Module[{LargerSymOption,SuperfluousSymOption,n,slotexpr,slotrules,internalsym,slotexprlist, slotlist1, slotlist2,HeadList, tensorinds,imposedorderinds,extendedtensorinds,perm1,perm2,newsym,newlabel, deltapos, deltareplacements},
{LargerSymOption,SuperfluousSymOption}={ImposeLargerSym,ImposeSuperfluousSym}/.CheckOptions[options]/.Options[ImposeSym];
If[And[Length@GroupSupport@sym==0,Not@SuperfluousSymOption],
expr//.{x:(CovarD[__][___]):>ExpandCovarD[x],x:(InertHeadHead[__][___]):>ExpandInertHeadHead[x]},
{n,slotexpr,slotrules,internalsym}=List@@SymmetryOf[expr];
(* Sort the tensors *)
slotexprlist=Sort@TimesToList[slotexpr];
slotlist1=Join@@List@@@slotexprlist; 
(* Remove slots that does not contain abstract indices *)
slotlist2=Select[slotlist1,Head[#]===xAct`xTensor`Private`slot&];
HeadList=Head/@slotexprlist;
(* We need to replace deltas with metrics or deltaH that carries its VBundle *)
deltapos=First/@Position[HeadList,delta,1];
If[Length@deltapos>0,
deltareplacements=MetricOrdeltaHOfVBundle/@VBundleOfIndex/@((First/@Part[slotexprlist,deltapos])/.slotrules);
HeadList=ReplacePart[HeadList,Thread@Rule[deltapos, deltareplacements]];
];
tensorinds=slotlist1/.slotrules;
(* If we have indices in the sym that are not in tensorinds *)
extendedtensorinds=Join[tensorinds,Complement[List@@imposedinds,tensorinds]];
imposedorderinds=Join[List@@imposedinds,Cases[tensorinds,Except[Alternatives@@imposedinds]]];
perm1=TranslatePerm[PermutationFromTo[imposedorderinds,extendedtensorinds],Rules];
newsym=sym/.List@@perm1;
(* The extra indices in the sym that are not in tensorinds, does not act on anything so stabilize these points.*) 
(* Changed to StabilizerSGS instead of Stabilizer *)
newsym=StabilizerSGS[Range[1+Length@tensorinds,Length@extendedtensorinds],newsym];
newsym=TidyGroup@newsym;
perm2=TranslatePerm[PermutationFromTo[slot/@Range@n,slotlist2],Rules];internalsym=internalsym/.List@@perm2;
internalsym=TidyGroup@internalsym;
(* Simplify if possible *)
If[MinusIDInProductQ[newsym,internalsym],
0,
If[And[Not@SuperfluousSymOption,SubgroupQ[newsym,internalsym]],
ExpandCovarDAndInertHeadHead[expr],
If[LargerSymOption,
Module[{LargerGroup, tempsymhead=SymH[HeadList,newsym,DefaultLabelSym[newsym]]},
LargerGroup=SymmetryGroupOfTensor@tempsymhead;
newlabel=DefaultLabelSym[LargerGroup];
SymH[HeadList,LargerGroup,newlabel][Sequence@@tensorinds]],
SymH[HeadList,newsym,DefaultLabelSym[newsym]][Sequence@@tensorinds]]]]]]

3.10. SmartSymmetrize and SmartAntisymmetrize

These are new tools to compute ToCanonical@Symmetrize[expr, inds] quicker. Observe that ToCanonical might still be needed, but the number of intermediate terms should be greatly reduced. In principle a general SmartImposeSymmetry should work, but my current code for computing the mixed symmetries for a general imposed symmetry is way too slow to be useful for this purpose.

SmartSymmetrize[expr_]:=SmartSymmetrize[expr,FindFreeIndices[expr]];

SmartSymmetrize[expr_,inds_]:=ExpandSym[ImposeSym[expr,IndexList@@inds,Symmetric@Range@Length@inds,ImposeSuperfluousSym->True],SmartExpand->True];

SmartAntisymmetrize[expr_]:=SmartAntisymmetrize[expr,FindFreeIndices[expr]];

SmartAntisymmetrize[expr_,inds_]:=ExpandSym[ImposeSym[expr,IndexList@@inds,Antisymmetric@Range@Length@inds,ImposeSuperfluousSym->True],SmartExpand->True];

3.11. More tools

Sorting the Heads in the SymH object

SortTensorsInSym1[SymH[headlist_,sym_,label_][inds___]]:=
Module[{sortedheads=Sort[headlist],
slotlist=IndexList@@slot/@Range@Length@List[inds],
numindices=Length/@SlotsOfTensor/@headlist,partitionedslots,symbolpermutation, newslots,slotpermutation,newinds,newsym,newlabel},
partitionedslots=IndexList@@@partitionRagged[List@@slotlist,numindices];(*partitionedslots=IndexList@@@Last/@Rest@FoldList[{Drop[#1[[1]],#2],Take[#1[[1]],#2]}&,{List@@slotlist,{}},numindices];*)
symbolpermutation=PermutationFromTo[headlist,sortedheads];
newslots=Flatten[List@@@PermuteList[partitionedslots,symbolpermutation]];
slotpermutation=PermutationFromTo[List@@slotlist,newslots];
newinds=PermuteList[List[inds],slotpermutation];
newsym=sym/.List@@TranslatePerm[slotpermutation,Rules];
(* Also sort the base of the sym *)
newsym=TidyGroup@newsym;
(* One should probably change the label too. *)
newlabel=DefaultLabelSym[newsym];
SymH[sortedheads,newsym,newlabel][Sequence@@newinds]]

SortTensorsInSym[expr_]:=
(expr/.x:SymH[___][___]:>SortTensorsInSym1[x])

More tools

SuperfluousSymQ[x:(SymH[headlist_,sym_,label_])]:=
SubgroupQ[sym,InternalSymmetrySym[x]]

RemoveSuperfluousSym[expr_]:=(expr/.x:(SymH[headlist_,sym_,label_][inds___]):>ExpandCovarDAndInertHeadHead[InternalExprSym[Head@x]/.Thread@Rule[slot/@Range[Length@List[inds]],List[inds]]]/;SuperfluousSymQ[Head@x])

RemoveSuperfluousInnerSym[expr_]:=ExpandCovarDAndInertHeadHead[expr/.SymH[{SymH[headlist_,innersym_,innerlabel_]},outersym_,outerlabel_]:> Module[{},
If[SubgroupQ[innersym,outersym],
RemoveTrivialSym@SymH[headlist,outersym,outerlabel],
SymH[{SymH[headlist,innersym,innerlabel]},outersym,outerlabel]]]]

RemoveTrivialSym[expr_]:=expr/.x:(SymH[headlist_,sym_,label_]):>ZeroTensor[SlotsOfTensor@x]/;MinusIDInProductQ[sym,InternalSymmetrySym@x]

The following code partially expands a Sym object. The symmetrizations in the expansion does not include "ExpansionIndex". This code needs to be verified theoretically.

ExpandSymOneIndex1[SymH[headlist_,sym_,label_][inds___],ExpansionIndex_]:=
If[Length@Position[List[inds],ExpansionIndex,1,1]==0,
SymH[headlist,sym,label][inds],
Module[{ExpansionSlot=First@First@Position[List[inds],ExpansionIndex,1,1], newsym},
newsym=BaseChange2[xAct`xTensor`Private`SGSofsym@sym,{ExpansionSlot}];
Module[{IndexSOrbit=SchreierOrbit[ExpansionSlot,newsym,PermLength@newsym],
PermsForOrbit,SlotChanges,GroupList,SignList,LabelList,SymHHeads,PermutedIndices,SymHList,FF},
PermsForOrbit=TraceSchreier[#,IndexSOrbit]&/@First[IndexSOrbit];
SignList=sign/@PermsForOrbit;
SlotChanges=List@@@TranslatePerm[nosign/@PermsForOrbit,Rules];
GroupList=TidyGroup/@(StabilizerChain2[newsym][[2]]/.SlotChanges);
SymHHeads=SymH[headlist,#,DefaultLabelSym@#]&/@GroupList;
PermutedIndices=FF@@PermuteList[List[inds],#]&/@PermsForOrbit;
SymHList=MapThread[(#1@@#2)&,{SymHHeads,PermutedIndices},1];
RemoveSuperfluousSym[Inner[Times,SignList,SymHList,Plus]/Length[First[IndexSOrbit]]]
]]]

ExpandSymOneIndexInternal[y:SymH[headlist_,sym_,label_][inds___],ExpansionIndex_]:=
Module[{SymHpositions, return=y},
SymHpositions=Position[headlist,SymH];
If[Length[SymHpositions]>0,
Module[{internalexpr=InternalExprSym@Head@y,expandedinternalexpr, dummylist,dummiestoinds,indstodummies},
dummylist=DummyAs/@List[inds];
dummiestoinds=Thread@Rule[dummylist,List[inds]];
indstodummies=Thread@Rule[List[inds],dummylist];
internalexpr=internalexpr/.Thread@Rule[slot/@Range@Length@List[inds],dummylist];
expandedinternalexpr=internalexpr/.x:SymH[___][___]:>ExpandSymOneIndex1[x,ExpansionIndex/.indstodummies]/.x:SymH[___][___]:>ExpandSymOneIndexInternal[x,ExpansionIndex/.indstodummies];
return=RemoveSuperfluousInnerSym[ImposeSym[Expand@expandedinternalexpr,IndexList@@dummylist,sym]/.dummiestoinds];
]];
return]

ExpandSymOneIndex[expr_,ExpansionIndex_]:=Module[{result},
result=(expr/.x:SymH[___][___]:>ExpandSymOneIndex1[x,ExpansionIndex]/.y:SymH[___][___]:>ExpandSymOneIndexInternal[y,ExpansionIndex]);
result]

ExpandSymOneSlot[SymH[headlist_,sym_,label_][inds___],ExpansionSlot_]:=
Module[{newsym=BaseChange2[xAct`xTensor`Private`SGSofsym@sym,{ExpansionSlot}]},
Module[{IndexSOrbit=SchreierOrbit[ExpansionSlot,newsym,PermLength@newsym],
PermsForOrbit,stabsym=TidyGroup[StabilizerChain2[newsym][[2]]], newinds, newlabel},
PermsForOrbit=InversePerm/@(TraceSchreier[#,IndexSOrbit]&/@First[IndexSOrbit]);
newinds=PermuteList[List[inds],#]&/@PermsForOrbit;
newlabel=DefaultLabelSym[stabsym];
1/Length[First[IndexSOrbit]]*Plus@@(SymH[headlist,stabsym,newlabel]@@@newinds)]]

There is still a problem if headlist contains a Sym object which symmetrizes over ExpansionIndex.....?

MoveCovDInsideSym=covd_?CovDQ[covdinds___]@SymH[headlist_, sym_, label_][inds___]:>Module[{indlist=List[inds],slotlist=xAct`SymManipulator`Private`slot/@Range@Length@List[inds],internalexpr,tempinds=DummyAs/@List[inds],newinternalexpr},
internalexpr=xAct`SymManipulator`Private`InternalExprSym[SymH[headlist,sym,label]];
newinternalexpr=Expand@xAct`SymManipulator`Private`ExpandCovarDAndInertHeadHead[internalexpr/.Thread@Rule[slotlist,tempinds]];
ImposeSym[covd[covdinds]@newinternalexpr,IndexList@@tempinds,sym]/.Thread@Rule[tempinds,indlist]];

MoveTensorsOutsideSym[expr_, onlymetrics_:False]:=
(expr/.x:(SymH[headlist_,sym_,label_][inds___]):>
Module[{stableslots,newlabel,tensornotinsymQ,
numindices=Length/@SlotsOfTensor/@headlist,slotlist,partitionedslots, tensornotinsymQpos,
symbolpermutation,newslots,slotpermutation, newinds,innerinds, innersym,innerheads, outerheads, newsym, newpartitionedslots, newnumindices, partitionednewinds,tensorinsymQpos},
slotlist=Range[Plus@@numindices];
(* Extract indices belonging to the different tensors. Can this be done in a simpler way? *)
(*partitionedslots=Last/@Rest@FoldList[{Drop[#1[[1]],#2],Take[#1[[1]],#2]}&,{slotlist,{}},numindices];*)
partitionedslots=partitionRagged[slotlist,numindices];
stableslots=Complement[Range@Length@List[inds],GroupSupport@sym];
tensornotinsymQ = subsetQ[#,stableslots]&/@partitionedslots;
(* If the flag onlymetrics is set we only move metrics outside. *)
If[onlymetrics,
tensornotinsymQ=Thread@And[tensornotinsymQ,MetricQ/@headlist];];
tensornotinsymQpos=Select[Range@Length@headlist,tensornotinsymQ[[#]]&];
outerheads=headlist[[tensornotinsymQpos]];
tensorinsymQpos=Complement[Range@Length@headlist,tensornotinsymQpos];
innerheads=headlist[[tensorinsymQpos]];
symbolpermutation=PermutationFromTo[Range@Length@headlist,Join[tensorinsymQpos, tensornotinsymQpos]];
newpartitionedslots=PermuteList[partitionedslots,symbolpermutation];
newslots=Flatten[newpartitionedslots];
slotpermutation=PermutationFromTo[slotlist,newslots];
newinds=PermuteList[List[inds],slotpermutation];
newnumindices=Length/@Take[newpartitionedslots, -Length@outerheads];
newnumindices=Prepend[newnumindices,Length[newinds]-Plus@@newnumindices];
partitionednewinds=IndexList@@@Last/@Rest@FoldList[{Drop[#1[[1]],#2],Take[#1[[1]],#2]}&,{List@@newinds,{}},newnumindices];
newsym=sym/.List@@TranslatePerm[slotpermutation,Rules];
(* Also sort the base of the sym *)
newsym=TidyGroup@newsym;
(* One should probably change the label too. *)
newlabel=DefaultLabelSym[newsym];
outerheads=Prepend[outerheads,SymH[innerheads,newsym,newlabel]];
Inner[(#1@@#2)&,outerheads,partitionednewinds,Times]]
)//.y:(CovarD[___][___]):>ExpandCovarD[y]

MoveTensorsInsideSym1[SymH[headlist_,sym_,label_][sinds___],hllist_List,indlist_List]:=
SymH[Join[headlist,hllist],sym,label][sinds,Sequence@@indlist]

MoveTensorsInsideSym[expr_]:=Module[{TempFF, intermediateexpr},
intermediateexpr=expr/.y:SymH[headlist_,sym_,label_][sinds___]:> TempFF[y,{},{}];
intermediateexpr=intermediateexpr//.Times[TT_?xTensorQ[tinds___],TempFF[y_,hllist_List,indlist_List],z___]:> Times[z,TempFF[y,Append[hllist,TT],Join[indlist,{tinds}]]];
intermediateexpr=intermediateexpr/.TempFF->MoveTensorsInsideSym1;
intermediateexpr/.x:SymH[___][___]:>SortTensorsInSym[x]]

RemoveDummiesFromSym[expr_Plus]:=RemoveDummiesFromSym/@expr

RemoveDummiesFromSym[expr_]:=Module[{dummies=FindDummyIndices[expr]},
dummies=Join[Times[-1,#]&/@dummies,dummies];
Fold[ExpandSymOneIndex,expr,dummies]]

MoveSymIndicesDown1[x:(SymH[headlist_,sym_,label_][sinds___]),expr2_]:= Module[{n=Length@List[sinds],unstableslots,unstableindices,needtomoveindices, metricfactors,indexchanges},
unstableslots=GroupSupport[sym];
unstableindices=List[sinds][[unstableslots]];
needtomoveindices=Select[unstableindices,DownIndex[#]=!=#&];
metricfactors=(First[MetricsOfVBundle@VBundleOfIndex[#]][#,DummyAs[#]])&/@needtomoveindices;
indexchanges=Rule[#[[1]],ChangeIndex[#[[2]]]]&/@metricfactors;
Times[ReplaceIndex[x,indexchanges],ContractMetric@Times[Times@@metricfactors,expr2]]
];

Needs to implement something that works with derivatives too ...

MoveSymIndicesDown[expr_]:=Module[{TempFF, intermediateexpr},
(* Derivatives *)
intermediateexpr=expr//.x:(D1_?CovDQ[___,-(A_Symbol),___][SymH[headlist_,sym_,label_][___,(A_Symbol),___]]):>xAct`xTensor`Private`SymmetryOfIndex[A]*ReplaceIndex[x,{A->-A,-A->A}];
(* Tensors *)
intermediateexpr=intermediateexpr/.x:SymH[headlist_,sym_,label_][sinds___]:> TempFF[x,1];
intermediateexpr=intermediateexpr//.Times[TempFF[x_,y_],z___]:> TempFF[x,Times[y,z]];
intermediateexpr=intermediateexpr/.TempFF->MoveSymIndicesDown1]

The following code contracts metrics inside the SymH object. This is new in release 0.8.2 and should initially be considered as experimental.

ContractMetricsInsideSym[expr_]:=
expr//.x:(SymH[(headlist:{___,TT_?MetricQ,___}),sym_,label_][inds___]):>
Module[{stableslots,newlabel,removetensor=0,
numindices=Length/@SlotsOfTensor/@headlist,slotlist,partitionedslots,
indlist=List[inds],metricpositions,metricslots,slottoind,contractionpos1,contractionpos2,firstcontractedmetric1,firstcontractedmetric2,contractionslots,contractionindex, stablecontraction,slotchange1,contractionsign=1,
symbolpermutation,newslots,slotpermutation, newinds,innerheads, newsym, newpartitionedslots},
slotlist=Range[Plus@@numindices];
metricpositions=First/@Position[MetricQ/@headlist,True];
(* If we have no metric just return *)
If[Length@metricpositions==0,x,
slottoind=Thread@Rule[slotlist,List[inds]];
(* Extract indices belonging to the different tensors. Can this be done in a simpler way? *)
partitionedslots=partitionRagged[slotlist,numindices];
stableslots=Complement[Range@Length@List[inds],GroupSupport@sym];
metricslots=Part[partitionedslots,metricpositions];
(* At which slot do we find an index contracted with the index in metricslots[[n,1]]? *)
contractionpos1=(First/@Position[indlist,-((#1)[[1]]/.slottoind),1])&/@metricslots;
(* At which slot do we find an index contracted with the index in metricslots[[n,2]]? *)
contractionpos2=(First/@Position[indlist,-((#1)[[2]]/.slottoind),1])&/@metricslots;
firstcontractedmetric1=First/@Position[(Length[#]>0)&/@contractionpos1,True,1];
firstcontractedmetric2=First/@Position[(Length[#]>0)&/@contractionpos2,True,1];
newpartitionedslots=partitionedslots;
If[Length@firstcontractedmetric1>0,
contractionslots={metricslots[[First@firstcontractedmetric1,1]],First@contractionpos1[[First@firstcontractedmetric1]]};
(* Print["Found contraction between the first index (at slot ",contractionslots[[1]],") of metric tensor ", headlist[[metricpositions[[First@firstcontractedmetric1]]]]," with the index at slot ", contractionslots[[2]]]; *)
stablecontraction=subsetQ[contractionslots,stableslots];
If[stablecontraction,
(* Print["Stable"]; *)
removetensor=metricpositions[[First@firstcontractedmetric1]];
slotchange1=partitionedslots[[removetensor,2]]->contractionslots[[2]];
contractionindex=contractionslots[[1]]/.slottoind;
contractionsign=If[UpIndexQ[contractionindex],xAct`xTensor`Private`SymmetryOfIndex[contractionindex],1];
If[contractionslots[[2]]==partitionedslots[[removetensor,2]],
(* Print["Self contraction"];*)
contractionsign=contractionsign*DimOfVBundle@VBundleOfIndex[contractionindex];]
,
(* Print["Not stable: ", x];*)
removetensor=0;
];
];
If[And[removetensor==0,Length@firstcontractedmetric2>0], 
contractionslots={metricslots[[First@firstcontractedmetric2,2]],First@contractionpos2[[First@firstcontractedmetric2]]};
(* Print["Found contraction between the second index (at slot ",contractionslots[[1]],") of metric tensor ", headlist[[metricpositions[[First@firstcontractedmetric2]]]]," with the index at slot ", contractionslots[[2]]]; *)
stablecontraction=subsetQ[contractionslots,stableslots];
If[stablecontraction,
(* Print["Stable"];*)
removetensor=metricpositions[[First@firstcontractedmetric2]];
slotchange1=partitionedslots[[removetensor,1]]->contractionslots[[2]];
contractionindex=contractionslots[[1]]/.slottoind;
contractionsign=If[UpIndexQ[contractionindex],1,xAct`xTensor`Private`SymmetryOfIndex[contractionindex]];
If[contractionslots[[2]]==partitionedslots[[removetensor,1]],Print["Self contraction"];
contractionsign=contractionsign*DimOfVBundle@VBundleOfIndex[contractionindex];]
,
(* Print["Not stable: ", x]; *)
removetensor=0;
];
];
If[removetensor==0,x,
innerheads=Delete[headlist,removetensor];
symbolpermutation=PermutationFromTo[Range@Length@headlist,Append[Delete[Range@Length@headlist,removetensor],removetensor]];
newpartitionedslots=PermuteList[partitionedslots,symbolpermutation];
newslots=Flatten[newpartitionedslots];
slotpermutation=PermutationFromTo[slotlist/.{slotchange1,Reverse@slotchange1},newslots];
newinds=Drop[PermuteList[List[inds],slotpermutation],-2];
newsym=sym/.List@@TranslatePerm[slotpermutation,Rules];
(* Also sort the base of the sym *)
newsym=TidyGroup@newsym;
(* One should probably change the label too. *)
newlabel=DefaultLabelSym[newsym];
RemoveSuperfluousSym[contractionsign*SymH[innerheads,newsym,newlabel][Sequence@@newinds]]]]]

ExpandInternalSym[expr_,options___]:=
expr/.x:(SymH[___][inds___]):>ExpandInternalSym1[x,options]

ExpandInternalSym1[x:(SymH[headlist_,sym_,label_][inds___]),options___]:=Module[{indlist=List[inds],slotlist=IndexList@@slot/@Range@Length@List[inds],slotrules},
slotrules=Thread@Rule[List@@slotlist,indlist];
ImposeSym[Expand@ExpandSym[InternalExprSym[Head@x]/.slotrules,options],indlist,sym,options]]

RuleInSym1[SymH[headlist_,sym_,label_][inds___],rule_]:=Module[{indlist=List[inds],slotlist=slot/@Range@Length@List[inds],internalexpr,tempinds=DummyAs/@List[inds],newinternalexpr},
internalexpr=InternalExprSym[SymH[headlist,sym,label]];
newinternalexpr=Expand[Expand@ExpandCovarDAndInertHeadHead[internalexpr/.Thread@Rule[slotlist,tempinds]]/.rule];
ImposeSym[newinternalexpr,IndexList@@tempinds,sym]/.Thread@Rule[tempinds,indlist]]

RuleInSym[expr_,rule_]:=expr/.x:(SymH[___][inds___]):>RuleInSym1[x,rule]

FunctionInSym1[SymH[headlist_,sym_,label_][inds___],func_]:=Module[{indlist=List[inds],slotlist=slot/@Range@Length@List[inds],internalexpr,tempinds=DummyAs/@List[inds],newinternalexpr},
internalexpr=InternalExprSym[SymH[headlist,sym,label]];
newinternalexpr=Expand@func[Expand[ExpandCovarDAndInertHeadHead[internalexpr/.Thread@Rule[slotlist,tempinds]]]];
ImposeSym[newinternalexpr,IndexList@@tempinds,sym]/.Thread@Rule[tempinds,indlist]]

FunctionInSym[expr_,func_]:=expr/.x:(SymH[___][inds___]):>FunctionInSym1[x,func]

Try to factor the imposed group, and remove factors that are not needed.

RemoveSuperfluousPartialSym1[x:(SymH[headlist_,sym_,label_])]:=
With[{intsym=InternalSymmetrySym[x],symfactors=GroupFactorization1@xAct`xTensor`Private`SGSofsym@sym},
With[{newsym=TidyGroup[JoinSGS@@(Select[symfactors,Not@SubgroupQ[#,intsym]&])]},
SymH[headlist,newsym,DefaultLabelSym@newsym]
]]

RemoveSuperfluousPartialSym[expr_]:=
(expr/.x:SymH[___]:>RemoveSuperfluousPartialSym1[x])

4. Computation of mixed symmetries

4.1. Symmetric case

Old code which is incomplete

SymmetricMixedSymmetry[unstableslots_,internalsym_,sym_,n_]:=Module[{internalstabilizer, newsym},
internalstabilizer=SetStabilizer[unstableslots,internalsym];
newsym=TidySGS@SchreierSims[{},Union[internalstabilizer[[2]],Append[sym,Cycles][[2]]]];
{newsym,internalstabilizer}]

Newer code. By setting antisym to True, one can also handle the antisymmetric case.

SymmetryOfProductOfGroupsSymmetricCase[Delta0_List,H0_StrongGenSet,antisym_:False]:=
Module[{Delta=Sort@Delta0,Deltacompl,G1,H1, chainH1,permdeg, mm,alpham,Hmm1,Kmm1,K1,F1, GElements, unstablepoints,time,testfunc},
(* Arrange so that the groups have the same bases and cycle notation *)
(* Print["SymmetryOfProductOfGroupsBruteForce[",G0,",",H0,"]"]; *)
If[subsetQ[GroupSupport@H0,Delta],{If[antisym,
Antisymmetric[Delta,Cycles],
Symmetric[Delta,Cycles]],H0},
permdeg=Max[Max@@First[H0],Max@Delta,PermDeg[Last@H0]];
mm=permdeg-Length@Delta;
Deltacompl=Complement[Range@permdeg,Delta];
H1=DeleteSomeRedundantGenerators@BaseChange2[H0,Join[Deltacompl,Delta]];
If[antisym,
G1=ReplacePart[Antisymmetric[Delta,Cycles],1->First@H1],
G1=ReplacePart[Symmetric[Delta,Cycles],1->First@H1]];
(*Print["H1: ",H1];
Print["G1: ",G1];*)
chainH1=StabilizerChain2[H1];
Hmm1=chainH1[[mm]];
alpham=Hmm1[[1,1]];
(*Print["permdeg: ",permdeg];
Print["mm: ",mm];
Print["Hm: ",StabilizerChain[H1][[1+mm]]];
Print["Hm-1: ",Hmm1];
Print["alpham: ",alpham];*)
If[Intersection[Orbit[alpham,Hmm1[[2]]],Delta]===Delta,Kmm1=Hmm1,Kmm1=ReplacePart[chainH1[[1+mm]],1->chainH1[[mm,1]]]];
 (*Print["Kmm1: ",Kmm1];*)
If[mm==1,K1=Kmm1,
testfunc=Length@Complement[Deltacompl,StablePoints[First@SpecialPermWord[#,H1,Length@Delta],permdeg]]==0&;
(* GElements=Rest[List@@Dimino[G1]];*)
(* It should be enough to test one from each coset Hm.p where p is in G1 *)
{time,GElements}=AbsoluteTiming[InversePerm/@Rest[TransversalComputation[chainH1[[1+mm]],ReplacePart[G1,1->chainH1[[1+mm,1]]]]]]; If[$MixedSymVerbose,Print["GElements timing: ",time, " length: ", Length@GElements]];
(*Print["GElements: ",GElements];
Print["testfunc: ",testfunc];
Print["Search ",{H1,GpinHGQBruteForceSymmetric[#,testfunc,GElements]&,1,Kmm1}];*)
K1=SearchSym1[chainH1,GpinHGQBruteForceSymmetric[#,testfunc,GElements]&,1,mm+1,Kmm1,Delta];];
(*Print["K1: ",K1];*)
F1=StrongGenSet[H1[[1]],Join[Complement[K1[[2]],StabilizerChain2[H1][[1+mm,2]]],G1[[2]]]];
F1=DeleteRedundantGenerators[F1];
unstablepoints=GroupSupport@F1;
(*Print["F1: ",F1];*)
{DeleteSomeRedundantGenerators@ReplacePart[BaseChange2[F1,unstablepoints],1->unstablepoints],DeleteSomeRedundantGenerators@ReplacePart[BaseChange2[K1,unstablepoints],1->unstablepoints]}]]

An improved version. Updated 2013-10-09 and 2015-04-26. Bug detected 2016-11-29.

SymmetryOfProductOfGroupsSymmetricCase2[Delta0_List,H0_StrongGenSet,antisym_:False]:=
Module[{Delta=Sort@Delta0,Deltacompl,G1,H1, chainH1,permdeg, mm,alpham,Hmm1,Kmm1,K1,F1, GElements, unstablepoints,time},
(* Arrange so that the groups have the same bases and cycle notation *)
(* Print["SymmetryOfProductOfGroupsBruteForce[",G0,",",H0,"]"]; *)
If[subsetQ[GroupSupport@H0,Delta],{If[antisym,
Antisymmetric[Delta,Cycles],
Symmetric[Delta,Cycles]],H0},
permdeg=Max[Max@@First[H0],Max@Delta,PermDeg[Last@H0]];
mm=permdeg-Length@Delta;
Deltacompl=Complement[Range@permdeg,Delta];
H1=DeleteSomeRedundantGenerators@BaseChange2[H0,Join[Deltacompl,Delta]];
If[antisym,
G1=ReplacePart[Antisymmetric[Delta,Cycles],1->First@H1],
G1=ReplacePart[Symmetric[Delta,Cycles],1->First@H1]];
(*Print["H1: ",H1];
Print["G1: ",G1];*)
chainH1=StabilizerChain2[H1];
Hmm1=chainH1[[mm]];
alpham=Hmm1[[1,1]];
(*Print["permdeg: ",permdeg];
Print["mm: ",mm];
Print["Hm: ",chainH1[[1+mm]]];
Print["Hm-1: ",Hmm1];
Print["alpham: ",alpham];*)
If[Intersection[Orbit[alpham,Hmm1[[2]]],Delta]===Delta,Kmm1=Hmm1,Kmm1=ReplacePart[chainH1[[1+mm]],1->chainH1[[mm,1]]]];
(* Print["Kmm1: ",Kmm1];*)
If[mm==1,K1=Kmm1,
(* GElements=Rest[List@@Dimino[G1]];*)
(* It should be enough to test one from each coset Hm.p where p is in G1 *)
{time,GElements}=AbsoluteTiming[InversePerm/@Rest[TransversalComputation[chainH1[[1+mm]],ReplacePart[G1,1->chainH1[[1+mm,1]]]]]]; If[$MixedSymVerbose,Print["GElements timing: ",time, " length: ", Length@GElements]];
(* Print["GElements: ",GElements];
Print["Search ",{chainH1,1,mm+1,Kmm1,Delta,GElements}]; *)
K1=SearchSym2[chainH1,SchreierOrbits[#,permdeg]&/@chainH1,1,mm+1,Kmm1,Delta,GElements];];
 (* Print["K1: ",K1];*)
F1=StrongGenSet[H1[[1]],Join[Complement[K1[[2]],chainH1[[1+mm,2]]],G1[[2]]]];
F1=DeleteSomeRedundantGenerators[F1];
unstablepoints=GroupSupport@F1;
(* Print["F1: ",F1];*)
{DeleteSomeRedundantGenerators@ReplacePart[BaseChange2[F1,unstablepoints],1->unstablepoints],DeleteSomeRedundantGenerators@ReplacePart[BaseChange2[K1,unstablepoints],1->unstablepoints]}]]

Attempted bugfix 2016-11-29. Only works with antisym=False.

SymmetryOfProductOfGroupsSymmetricCase2b[Delta0_List,H0_StrongGenSet,antisym_:False]:=
Module[{Delta=Sort@Delta0,Deltacompl,H0b,G1,H1, chainH1,origpermdeg,permdeg, mm,alpham,Hmm1,Kmm1,K1,F1, GElements, unstablepoints, unstablepointsnosign,time},
(* Arrange so that the groups have the same bases and cycle notation *)
If[antisym,Print["Warning: SymmetryOfProductOfGroupsSymmetricCase2b called with antisym=True"]];
If[subsetQ[GroupSupport@H0,Delta],{If[antisym,
Antisymmetric[Delta,Cycles],
Symmetric[Delta,Cycles]],H0},
origpermdeg=Max[Max@@First[H0],Max@Delta,PermDeg[Last@H0]];
H0b=xAct`xPerm`Private`FromSign[H0,origpermdeg];
permdeg=Max[Max@@First[H0b],Max@Delta,PermDeg[Last@H0b]];
mm=permdeg-Length@Delta;
Deltacompl=Complement[Range@permdeg,Delta];
H1=DeleteSomeRedundantGenerators@BaseChange2[H0b,Join[Deltacompl,Delta]];
If[antisym,
G1=DeleteSomeRedundantGenerators@BaseChange2[xAct`xPerm`Private`FromSign[Antisymmetric[Delta,Cycles],origpermdeg],First@H1],
G1=ReplacePart[Symmetric[Delta,Cycles],1->First@H1]];
(*Print["H1: ",H1];
Print["G1: ",G1];*)
chainH1=StabilizerChain2[H1];
Hmm1=chainH1[[mm]];
alpham=Hmm1[[1,1]];
(*Print["permdeg: ",permdeg];
Print["mm: ",mm];
Print["Hm: ",chainH1[[1+mm]]];
Print["Hm-1: ",Hmm1];
Print["alpham: ",alpham];
Print["Orbit alpham: ",Orbit[alpham,Hmm1[[2]]]];*)
(* The following line is probably wrong if Hmm1 has signed permutations...*)
If[Intersection[Orbit[alpham,Hmm1[[2]]],Delta]===Delta,Kmm1=Hmm1,Kmm1=ReplacePart[chainH1[[1+mm]],1->chainH1[[mm,1]]]];
 (*Print["Kmm1: ",Kmm1];*)
If[mm==1,K1=Kmm1,
(* GElements=Rest[List@@Dimino[G1]];*)
(* It should be enough to test one from each coset Hm.p where p is in G1 *)
{time,GElements}=AbsoluteTiming[InversePerm/@Rest[TransversalComputation[chainH1[[1+mm]],ReplacePart[G1,1->chainH1[[1+mm,1]]]]]]; If[$MixedSymVerbose,Print["GElements timing: ",time, " length: ", Length@GElements]];
(* Print["GElements: ",GElements];
Print["Search ",{chainH1,1,mm+1,Kmm1,Delta,GElements}]; *)
K1=SearchSym2[chainH1,SchreierOrbits[#,permdeg]&/@chainH1,1,mm+1,Kmm1,Delta,GElements];];
  (*Print["K1: ",K1];*)
F1=StrongGenSet[H1[[1]],Join[Complement[K1[[2]],chainH1[[1+mm,2]]],G1[[2]]]];
F1=DeleteSomeRedundantGenerators[F1];
unstablepoints=GroupSupport@F1;
unstablepointsnosign=DeleteCases[unstablepoints,origpermdeg+1|origpermdeg+2];
(*Print["F1: ",F1];
 Print["unstablepoints: ",unstablepoints];*)
{DeleteSomeRedundantGenerators@ReplacePart[xAct`xPerm`Private`ToSign[BaseChange2[F1,unstablepoints],origpermdeg],1->unstablepointsnosign],DeleteSomeRedundantGenerators@ReplacePart[xAct`xPerm`Private`ToSign[BaseChange2[K1,unstablepoints],origpermdeg],1->unstablepointsnosign]}]]

SearchSym1[chain_List,property_,s_Integer, l_Integer,SGSK_,Delta_,options___]:=Module[{base=chain[[1,1]],newSGSK=SGSK,k=Length[chain]-1,Deltas,gammas,Korbit,verb},
verb=xPermVerbose/.CheckOptions[options]/.Options[Search];
If[verb, Print["SearchSym1[",chain,",",property,",",s,",", l,",",SGSK,",",Delta,"]"]];
If[s==l-1,
newSGSK=SGSK,
newSGSK=SearchSym1[chain,property,s+1,l,SGSK,Delta,options];
Deltas=Orbit[base[[s]],chain[[s,2]]];
If[Intersection[Deltas,Delta]=!=Delta,Deltas=Complement[Deltas,Delta]];
If[verb,Print["Branching over points ",Deltas, " s=",s]];
(* Avoid rechecking permutations. Not in Butler's algorithm *)
(*If[s=!=k,Deltas=Drop[Deltas,1]];*)
If[s=!=l-2,Deltas=Drop[Deltas,1]];
Do[
gammas=Deltas[[j]];
Korbit=Orbit[gammas,newSGSK[[2]]];
If[verb,Print["gammas ", gammas, " Korbit ", Korbit]];
If[gammas===MinB[Korbit,base],
newSGSK=GenerateSym1[chain,property,s,l,s+1,Append[Take[base,s-1],gammas],newSGSK,Delta,verb];
If[Head[newSGSK]===Times,newSGSK=-newSGSK]]
,{j,1,Length[Deltas]}];
];
If[verb,Print["Return SGSK: ", newSGSK]];
newSGSK
];

GenerateSym1[chain_List,property_,s_,l_,i_,list_,SGSK_,Delta_,verb_]:=Module[{base=chain[[1,1]],k=Length[chain]-1,g,newSGSK=SGSK,Deltag,gammai,otherSGSK,schvecs=SchreierOrbits[#,PermLength[chain[[1]]]]&/@chain},
(* If[verb, Print["GenerateSym1[",chain,",",property,",",s,",",l,",",i,",",list,",",SGSK,",",Delta,",",verb,"]"]];*)
g=First@FromBaseImage2[list,schvecs];
If[verb,Print["Pregenerate at level ",s," with i=",i,". We have list ",list," and permutation ",g]];
If[i==l,
If[verb,Print["Generate at level ",s," with i=",i,". We have list ",list," and permutation ",g]];
If[Not@PermEqual[g,ID]&&Or[Length@Intersection[Delta,list]==0,property[g]]&&Not@PermMemberQ[g,SGSK],
newSGSK=-StrongGenSet[chain[[s,1]],Append[Last[SGSK],g]];
If[verb,Print["  Added permutation ",g," New SGSK: ", newSGSK]]],
Deltag=OnPoints[Orbit[base[[i]],chain[[i,2]]],g];
(* Here we should check if a point in Deltag is in Delta, 
then the entire Delta has to be in the orbit in chain[[s]]. 
Otherwise g can not be in the group. 
Can we refine this to the orbit in chain[[i]]??? *)
Deltag=Select[Deltag,Or[FreeQ[Delta,#],subsetQ[Delta,Orbit[#,chain[[s,2]]]]]&];
If[verb,Print["Generating over points ",Deltag, ", s=", s, ", i=",i,", list ",list]];
Do[
gammai=Deltag[[j]];
newSGSK=GenerateSym1[chain,property,s,l,i+1,Append[list,gammai],newSGSK,Delta,verb];
If[Head[newSGSK]===Times,Break[]],
{j,1,Length[Deltag]}];
];
newSGSK
];

SearchSym2[chain_List,schvecs_,s_Integer, l_Integer,SGSK_,Delta_,GElements_,options___]:=Module[{base=chain[[1,1]],newSGSK=SGSK,k=Length[chain]-1,Deltas,gammas,Korbit,verb},
verb=xPermVerbose/.CheckOptions[options]/.Options[Search];
If[verb, Print["SearchSym2[",chain,",",schvecs,",",s,",", l,",",SGSK,",",Delta,",",GElements,"]"]];
If[s==l-1,
newSGSK=SGSK,
newSGSK=SearchSym2[chain,schvecs,s+1,l,SGSK,Delta,GElements,options];
Deltas=Orbit[base[[s]],chain[[s,2]]];
If[Intersection[Deltas,Delta]=!=Delta,Deltas=Complement[Deltas,Delta]];
If[verb,Print["Branching over points ",Deltas, " s=",s]];
(* Avoid rechecking permutations. Not in Butler's algorithm *)
(*If[s=!=k,Deltas=Drop[Deltas,1]];*)
If[s=!=l-2,Deltas=Drop[Deltas,1]];
Do[
gammas=Deltas[[j]];
Korbit=Orbit[gammas,newSGSK[[2]]];
If[verb,Print["gammas ", gammas, " Korbit ", Korbit]];
If[gammas===MinB[Korbit,base],
newSGSK=GenerateSym2[chain,schvecs,s,l,s+1,Append[Take[base,s-1],gammas],newSGSK,Delta,GElements,verb];
If[Head[newSGSK]===Times,newSGSK=-newSGSK]]
,{j,1,Length[Deltas]}];
];
If[verb,Print["Return SGSK: ", newSGSK]];
newSGSK
];

GenerateSym2[chain_List, schvecs_,s_,l_,i_,list_,SGSK_,Delta_,GElements_,verb_]:=Module[{base=chain[[1,1]],k=Length[chain]-1,g,newSGSK=SGSK,Deltag,gammai,otherSGSK,Deltacompl},
(* If[verb, Print["GenerateSym1[",chain,",",property,",",s,",",l,",",i,",",list,",",SGSK,",",Delta,",",verb,"]"]];*)
Deltacompl=Drop[base,-Length@Delta];
g=First@FromBaseImage2[list,schvecs];
If[verb,Print["Pregenerate at level ",s," with i=",i,". We have list ",list," and permutation ",g]];
If[i==l,
If[verb,Print["Generate at level ",s," with i=",i,". We have list ",list," and permutation ",g]];
If[Not@PermEqual[g,ID]&&Or[Length@Intersection[Delta,list]==0,GpinHGQBruteForceSymmetric2[list,Deltacompl, schvecs,GElements]]&&Not@PermMemberQ[g,SGSK],
newSGSK=-StrongGenSet[chain[[s,1]],Append[Last[SGSK],g]];
If[verb,Print["  Added permutation ",g," New SGSK: ", newSGSK]]],
Deltag=OnPoints[Orbit[base[[i]],chain[[i,2]]],g];
(* Here we should check if a point in Deltag is in Delta, 
then the entire Delta has to be in the orbit in chain[[s]]. 
Otherwise g can not be in the group. 
Can we refine this to the orbit in chain[[i]]??? *)
Deltag=Select[Deltag,Or[FreeQ[Delta,#],subsetQ[Delta,Orbit[#,chain[[s,2]]]]]&];
If[verb,Print["Generating over points ",Deltag, ", s=", s, ", i=",i,", list ",list]];
Do[
gammai=Deltag[[j]];
newSGSK=GenerateSym2[chain,schvecs,s,l,i+1,Append[list,gammai],newSGSK,Delta,GElements,verb];
If[Head[newSGSK]===Times,Break[]],
{j,1,Length[Deltag]}];
];
newSGSK
];

SpecialPermWord[perm_,StrongGenSet[base_List,GS_GenSet],Deltalength_,word_List:{}]:=If[Length[base]≤Deltalength,
Prepend[word,perm],
Module[{
Sorbit=SchreierOrbit[First[base],GS,Max[PermDeg[perm],PermDeg[GS],Max[base]]],
point=OnPoints[First[base],perm],
u},
If[MemberQ[First[Sorbit],point],
u=TraceSchreier[point,Sorbit];
SpecialPermWord[PermProduct[perm,InversePerm[u]],StrongGenSet[Rest[base],Stabilizer[{First[base]},GS]],Deltalength,Prepend[word,u]],
Prepend[word,perm]]
]
];

SpecialPermMemberQ[p_,H_,Delta_,n_]:=Length@Complement[Complement[Range@n,Delta],StablePoints[First@SpecialPermWord[p,H,Length@Delta],n]]==0

GpinHGQBruteForceSymmetric[p_,testfunc_,Gelements_]:=
Module[
{TestList=Gelements,
result=True},
While[(Length[TestList]>0)&&result,
result=testfunc@PermProduct[p,First@TestList];
TestList=Rest@TestList;];
result
]

GpinHGQBruteForceSymmetric2[list_,deltacompl_, schvecs_,Gelements_]:=
Module[
{TestList=Gelements,
result=True},
While[(Length[TestList]>0)&&result,
result=deltacompl==Last@FromBaseImage2[OnPoints[list,First@TestList],schvecs];
(* Print["Tested ", First@TestList, ": ", result];*)
TestList=Rest@TestList;];
result
]

FromBaseImage2[images_List,schvecs_]:=Module[{i,u,g=ID,imgs=images},
For[i=1,i≤Length[imgs],i++,
If[MemberQ[schvecs[[i,1]],imgs[[i]]],
u=TraceSchreier[imgs[[i]],schvecs[[i]]];
g=PermProduct[u,g];
imgs=OnPoints[imgs,InversePerm[u]]
]];
{g,imgs}
]

4.2. Antisymmetric case

Old code which is incomplete

AntisymmetricMixedSymmetry[unstableslots_,internalsym_,sym_,n_]:=Module[{internalstabilizer, newsym},
internalstabilizer=SetStabilizer[unstableslots,internalsym];
newsym=TidySGS@SchreierSims[{},Union[internalstabilizer[[2]],Append[sym,Cycles][[2]]]];
{newsym,internalstabilizer}]

4.1. More general cases

GeneralMixedSymmetry[unstableslots_,internalsym_,sym_,n_]:=Module[{time,result,longorbits,joinedsymmetric},
If[SubgroupQ[internalsym,sym],{xAct`xTensor`Private`SGSofsym@sym,internalsym},
{time,result}=AbsoluteTiming@SymmetryOfProductOfGroupsBruteForce[sym,internalsym];
If[$MixedSymVerbose,
Print["GeneralMixedSymmetry: ", time, "s."]
];
result]]

5. Canonicalization of the group

5.1. Canonicalization for nested SymH objects

FindInternalCanonicalSymGroupAndPerm[SymH[headlist:{SymH[intheadlist_, intsym_, intlabel_]},sym_,label_]]:=Module[{newinternalSymH, intcanonperm,newsym,newSymH, canonicalizingperm, tmpexpr},
{newinternalSymH, intcanonperm}=FindInternalCanonicalSymGroupAndPerm[SymH[intheadlist, intsym, intlabel]];
newsym=sym/.List@@(TranslatePerm[nosign@intcanonperm,Rules]);
newsym=TidyGroup@newsym;
tmpexpr=SymH[{newinternalSymH},newsym,""];
{newSymH, canonicalizingperm}=FindCanonicalSymGroupAndPerm[tmpexpr];
{newSymH, PermProduct[intcanonperm,canonicalizingperm]}]

FindInternalCanonicalSymGroupAndPerm[SymH[headlist_,sym_,label_]]:=FindCanonicalSymGroupAndPerm[SymH[headlist,sym,label]]

FindCanonicalSymGroupAndPerm[SymH[headlist_,sym_,label_]]:=Module[{n=Plus@@(Length/@SlotsOfTensor/@headlist),unstableslots, newsymh, perm},
unstableslots=GroupSupport@sym;
{newsymh,perm}=Which[Head@sym===Symmetric,(* Symmetric *)
CanonicalizeSymmetricGroupInSym[SymH[headlist,sym,label],n,unstableslots],
Head@sym===Antisymmetric,(* Antisymmetric *)
CanonicalizeAntisymmetricGroupInSym[SymH[headlist,sym,label],n,unstableslots],
True,(* Everything else *)
CanonicalizeGeneralGroupInSym[SymH[headlist,sym,label],n]]]

CanonicalizeGroupInSym[expr_]:=expr/.SymH[headlist_,sym_,label_][inds___]:>Module[
{newSymH, canonicalizingperm},
{newSymH,canonicalizingperm}=FindInternalCanonicalSymGroupAndPerm[SymH[headlist,sym,label]];
sign[canonicalizingperm]*newSymH@@PermuteList[{inds},canonicalizingperm]]

Options[ToCanonicalSym]={RemoveSuperfluousPartialSym->False}

ToCanonicalSym[expr_,options:OptionsPattern[]]:=ToCanonical@With[{remsuppartsym=OptionValue[RemoveSuperfluousPartialSym]},If[remsuppartsym,CanonicalizeGroupInSym@RemoveSuperfluousPartialSym@RemoveSuperfluousSym@expr,CanonicalizeGroupInSym@expr]]

5.2. Canonicalization for Symmetric and Antisymmetric groups

CanonicalizeSymmetricGroupInSym[SymH[headlist_,sym_,label_],n_,unstableslots_]:=
Module[
{CanonicalSym=sym,newlabel=label,InternalSym, canonicalslots,canonicalizingperm},
InternalSym=InternalSymmetrySym[SymH[headlist,sym,label]];
 {canonicalslots,canonicalizingperm}=FindSmallestImageAndMapping[unstableslots,InternalSym,n];
CanonicalSym=Symmetric@Sort@canonicalslots;
newlabel=StringJoin["(",ToString/@canonicalslots,")"];
{SymH[headlist,CanonicalSym,newlabel],canonicalizingperm}]

CanonicalizeAntisymmetricGroupInSym[SymH[headlist_,sym_,label_],n_,unstableslots_]:=
Module[
{CanonicalSym=sym,newlabel=label,InternalSym, canonicalslots,canonicalizingperm},
InternalSym=InternalSymmetrySym[SymH[headlist,sym,label]];
 {canonicalslots,canonicalizingperm}=FindSmallestImageAndMapping[unstableslots,InternalSym,n];
CanonicalSym=Antisymmetric@Sort@canonicalslots;
newlabel=StringJoin["[",ToString/@canonicalslots,"]"];
{SymH[headlist,CanonicalSym,newlabel],canonicalizingperm}]

5.3. Canonicalization for general groups

This is experimental and very slow.

The original version pushes the symmetrizations to the right, but we need them pushed to the left. Therefore we reverse the points. This is still not exactly what we want...

Older code

CanonicalizeGeneralGroupInSym[SymH[headlist_,sym_,label_],n_]:=
Module[{time,result},
{time,result}=AbsoluteTiming@Module[{newlabel=label, newsym, InternalSym, InternalSymElements, PermIndex,newunstableslots,reverserule=Thread@Rule[Range@n,Reverse@Range@n],intsymel,reversesym},
(* Find a short representation of the reversed version of sym. This makes the following computations faster. *)
reversesym=StrongGenSet[#[[1,1]],GenSet@@xAct`SymManipulator`Private`SmallestGenSet@#]&@xAct`SymManipulator`Private`SmallestBaseStabChain[sym/.reverserule];
InternalSym=InternalSymmetrySym[SymH[headlist,sym,label]]/.reverserule;
InternalSymElements=List@@Dimino@InternalSym;
{PermIndex, newsym}=SmallestGroup[reversesym/.List@@@TranslatePerm[nosign/@InternalSymElements,Rules]];
newsym=TranslatePerm[StrongGenSet[#[[1,1]],GenSet@@xAct`SymManipulator`Private`SmallestGenSet@#]&@xAct`SymManipulator`Private`SmallestBaseStabChain[TranslatePerm[newsym,Cycles]/.reverserule],Cycles];
(*newsym=TranslatePerm[TranslatePerm[newsym,Rules]/.reverserule,Cycles];*)
newsym=MapAt[Sort,newsym,{2}];
intsymel=(InternalSymElements[[PermIndex]])/.reverserule;
newlabel=DefaultLabelSym[newsym];
newunstableslots=Complement[Range@n,StablePoints[Last@newsym,n]];
newunstableslots=DeleteDuplicates@Join[First@newsym,newunstableslots];
newsym=ReplacePart[newsym,1->newunstableslots];
{SymH[headlist,newsym,newlabel],intsymel}];
(*Print["CanonicalizeGeneralGroupInSym: ", time, "s."];*)
result]

Newer code which does not work properly yet.

FindCanonicalGroup[G1_,H1_,n_,perm_:Cycles[]]:=Module[{G2=StrongGenSet[#[[1,1]],GenSet@@xAct`SymManipulator`Private`SmallestGenSet@#]&@xAct`SymManipulator`Private`SmallestBaseStabChain[G1],longorbits=Select[Sort/@Orbits@H1,Length@#>1&],firstorbitpoint, firstschreierorbit, firstperms,PermIndex,G3},
If[Length@longorbits==0,{G2,perm},
firstorbitpoint=First@Sort[First/@longorbits];
firstschreierorbit=SchreierOrbit[firstorbitpoint,Last@H1,n];
firstperms=TraceSchreier[#,firstschreierorbit]&/@First@firstschreierorbit;
{PermIndex, G3}=xAct`SymManipulator`Private`SmallestGroup[G2/.List@@@TranslatePerm[xAct`SymManipulator`Private`nosign/@firstperms,Rules]];
FindCanonicalGroup[TranslatePerm[G3,Cycles],Stabilizer[{firstorbitpoint},H1],n,PermProduct[perm,firstperms[[PermIndex]]]]]]

CanonicalizeGeneralGroupInSym2[SymH[headlist_,sym_,label_],n_]:=
Module[{time,result},
{time,result}=AbsoluteTiming@Module[{newlabel=label, newsym, InternalSym, InternalSymElements, PermIndex,newunstableslots,reverserule=Thread@Rule[Range@n,Reverse@Range@n],intsymel,reversesym, firstorbitpoint},
(* Find a short representation of the reversed version of sym. This makes the following computations faster. *)
reversesym=StrongGenSet[#[[1,1]],GenSet@@xAct`SymManipulator`Private`SmallestGenSet@#]&@xAct`SymManipulator`Private`SmallestBaseStabChain[sym/.reverserule];
InternalSym=xAct`SymManipulator`Private`InternalSymmetrySym[SymH[headlist,sym,label]]/.reverserule;
{newsym,intsymel}=FindCanonicalGroup[reversesym,InternalSym,n];
newsym=TranslatePerm[TranslatePerm[newsym,Rules]/.reverserule,Cycles];
newsym=MapAt[Sort,newsym,{2}];
intsymel=intsymel/.reverserule;
newlabel=xAct`SymManipulator`Private`DefaultLabelSym[newsym];
newunstableslots=Complement[Range@n,StablePoints[Last@newsym,n]];
newunstableslots=DeleteDuplicates@Join[First@newsym,newunstableslots];
newsym=ReplacePart[newsym,1->newunstableslots];
{SymH[headlist,newsym,newlabel],intsymel}];
(* Print["CanonicalizeGeneralGroupInSym2: ", time, "s."];*)
result]

6. Irreducible decomposition

6.1. IrrDecompose

Check that the spinors package is loaded. Expand the expression before irreducible decomposition.

Options[IrrDecompose]^={Verbose->False,IrrDecSpinBundles->All, UseSym->True,IrrDecImposeNewMethod->False,TimeVerbose->False, ResultType->"Expression"};

IrrDecompose[expr_,options:OptionsPattern[]]:=Module[{verb,AllIrrDecSpinBundles,spinbundlesoption,SpinVBs,expandedexpr=Expand[expr], irrdecexpr1,MethodOption,TimeVerboseOption, usesymoption,time,resulttype},
If[Not@MemberQ[$ContextPath,"xAct`Spinors`"],
Throw@Message[IrrDecompose::error,"The Spinors package is needed for this. Please load the package Spinors first."]];
{MethodOption,TimeVerboseOption,spinbundlesoption,usesymoption,resulttype,verb}=OptionValue[{IrrDecImposeNewMethod,TimeVerbose,IrrDecSpinBundles, UseSym, ResultType,Verbose}];
AllIrrDecSpinBundles=xAct`Spinors`VBundleOfSolderingForm/@xAct`Spinors`$SolderingForms;AllIrrDecSpinBundles=Union[AllIrrDecSpinBundles,Dagger/@AllIrrDecSpinBundles];
SpinVBs=AllIrrDecSpinBundles;
If[Not[spinbundlesoption===All],
SpinVBs=Intersection[AllIrrDecSpinBundles,spinbundlesoption]];
{time, irrdecexpr1} = AbsoluteTiming@IrrDecompose2[expandedexpr,{SpinVBs,verb}];
If[TimeVerboseOption,Print["IrrDecompose step 1 timing:", time]];
If[And[usesymoption,Head@expandedexpr=!=Plus],
Module[{freeinds, freesym},
{freeinds, freesym}=SymmetryOfFreeIndices[expandedexpr,SpinVBs];
If[MethodOption,
{time, irrdecexpr1} = AbsoluteTiming@ToCanonical@ImposeSymmetryOnIrrdecTerms[irrdecexpr1,freeinds, freesym,SpinVBs],
{time, irrdecexpr1} = AbsoluteTiming@ToCanonical@ImposeSymmetry[irrdecexpr1,freeinds,freesym]];
If[TimeVerboseOption,Print["IrrDecompose step 2 timing:", time]];
]];
Switch[resulttype,
"Expression",irrdecexpr1,
"Equation",Equal[expr,irrdecexpr1],
"Rule",MakeRule[Evaluate[{expr,irrdecexpr1},MetricOn->All,ContractMetrics->True]],
"CompareRule",MakeCompareRule[Evaluate[{expr,irrdecexpr1}]],
_,Null]]

For backwards compatibility we also define

CompleteIrrDecompose[expr_,opt___]:=IrrDecompose[expr,UseSym->True,opt]

Thread over plus.

IrrDecompose2[expr_Plus,opt_]:=IrrDecompose2[#,opt]&/@expr

IrrDecompose2[s_?ScalarQ,opt_]:=s

IrrDecompose2[expr_,{SpinVBs_,verb_}]:=
Module[{frees=FindFreeIndices@expr,n,slotexpr,slotrules,sym,freeslots,totalsymoffrees,return},
(* Finds the free indices that are in SpinVBs *)
frees=Select[frees,MemberQ[SpinVBs,VBundleOfIndex@#]&];
{n,slotexpr,slotrules,sym}=List@@SymmetryOf[expr];
(* Which slots contains the free indices *)
freeslots=First/@List@@Replace[List@@frees, Reverse/@slotrules, {1}];
totalsymoffrees=CompatibleSymmetric[frees]/.Thread@Rule[Range@Length@frees,freeslots];
(* Check if expr is already irreducible, that is if it has maximal symmetry on its free indices *)
If[verb,Print["IrrDecompose2[",expr,",{",SpinVBs,",",verb,"}]"]];
return=If[SubgroupQ[totalsymoffrees,sym],
expr,
(* If it is not irreducible, rewrite ImposeSym[expr,frees] with DecomposeStep2[expr,frees]
  After that metrics are moved outside the sym if they are not symmetrized upon. *)
ToCanonical[
(expr+ToCanonicalSym@RemoveSuperfluousInnerSym@MoveTensorsOutsideSym[ImposeSym[expr,frees],True]-DecomposeStep2[expr,frees,{SpinVBs,verb}])]];
return]

DecomposeStep2[expr, inds] represents the same expression as ImposeSym[expr, inds] but written as expr + epsilon times traces.

DecomposeStep2[expr_, inds_,{SpinVBs_,verb_}] := Module[{frees = FindFreeIndices[expr], splitindices},
If[verb,Print["DecomposeStep2[",expr,",",inds,",{",SpinVBs,",",verb,"}]"]];
  If[Length[inds] < 2,
(* If we have less than two indices we can not rewrite the expression *)
expr,
splitindices = Select[inds, VBundleOfIndex[#] == VBundleOfIndex@First@inds &];
If[Length[splitindices] < 2,
(* If we only have one index in this VBundle, then we can not do anything with this index,
   so we continue the procedure with the rest of the indices. *)
DecomposeStep2[expr, Rest@inds,{SpinVBs,verb}],
(* If we more indices in this VBundle, then we do something more complicated. *)
Module[{tmpexpr, tmpexpr2, tmpexpr3, tmpexpr4, tmpexpr4b, tmpexpr5, tmpexpr6, finaltmpexpr, dummyindex, permutationrules},
(* First, make a recursive call to guarantee that tmpexpr represents ImposeSym[expr,Rest@inds]
   written as expr + epsilon times traces *)
tmpexpr = DecomposeStep2[expr, Rest@inds,{SpinVBs,verb}];
(* Let tmpexpr2 be the same thing but represented as a SymH object *)
tmpexpr2 = ImposeSym[expr, Rest@inds, CompatibleSymmetric[Rest@inds]];
(* We want to take a trace between splitindices[[1]] and splitindices[[2]].
 First we remove these two indices from the symmetry, 
   and move possible metrics outside the SymH object. *)
If[verb,Print["DecomposeStep2[",expr,",",inds,"...]: tmpexpr=",tmpexpr];
Print["DecomposeStep2[",expr,",",inds,"...]: tmpexpr2=",tmpexpr2]];
tmpexpr3 =MoveTensorsOutsideSym[Expand@ExpandSymOneIndex[ExpandSymOneIndex[tmpexpr2, splitindices[[1]]], splitindices[[2]]],True];
(* Generate a new dummy index for the trace. *)
dummyindex =DownIndex@DummyIn@VBundleOfIndex[splitindices[[2]]];
(* Make the contraction. *)
tmpexpr4 = ReplaceIndex[Evaluate@tmpexpr3, {splitindices[[1]] -> -dummyindex, splitindices[[2]] -> dummyindex}];
If[verb,Print["DecomposeStep2[",expr,",",inds,"...]: tmpexpr4=",tmpexpr4]];
(* If we now got a metric with one index contracted and the other in the symmetry,
   we can contract this metric. After that, the expression is canonicalized. *)
tmpexpr4b = ToCanonicalSym@ContractMetricsInsideSym@tmpexpr4;
(* Find the Irreducible decomposition of this symmetrized and contracted expression. *)
tmpexpr5 = IrrDecompose2[tmpexpr4b,{SpinVBs,verb}];
(* Multiply this with the epsilon of the first two splitindices *)
tmpexpr6 = Expand[(First@MetricsOfVBundle@VBundleOfIndex[splitindices[[1]]])[splitindices[[1]], splitindices[[2]]]*tmpexpr5];
(* The idea is now that expr can be written as tmpexpr + c*ImposeSym[tmpexpr6, Rest@splitindices].
The problem is that the second term will contain epsilons with symmetrizations. 
We want to avoid that. We notive that tmpexpr6 is symmetric in Rest@Rest@splitinices.
The symmetrization can therefore be represented by the permutations that interchange 
 splitindices[[2]] with Rest@splitindices *)
(* Construct the rules that interchange splitindices[[2]] with Rest@splitindices. *)
permutationrules = {Rule[splitindices[[2]], splitindices[[#]]], Rule[splitindices[[#]], splitindices[[2]]]}& /@ Range[2, Length[splitindices]];
(* The final expression therefore has this form. *)
finaltmpexpr = Expand[tmpexpr + 1/Length[splitindices]*Plus @@ (tmpexpr6 /. permutationrules)];
finaltmpexpr]]]]

SymmetryOfFreeIndices[expr_,SpinVBs_]:=Module[{freeinds=FindFreeIndices@expr, 
freeinvbslots, restslots, slotexpr,slotrules,sym,n,freefirstperm, freesym,freeinvb,restfrees},
{n,slotexpr,slotrules,sym}=List@@SymmetryOf[expr];
freeinvb=Select[freeinds,MemberQ[SpinVBs,UpIndex@VBundleOfIndex[#]]&];
freeinvbslots=First/@List@@Replace[List@@freeinvb, Reverse/@slotrules, {1}];
restslots=Complement[Range@n,freeinvbslots];
freefirstperm=List@@TranslatePerm[PermutationFromTo[Range@n,Join[freeinvbslots,restslots]],Rules];
(* The stabilizer only gives a subset of the symmetry, but might be good enough. *)
freesym=StabilizerSGS[restslots,sym];
freesym[[1]]=Complement[freesym[[1]],restslots];
freesym=freesym/.freefirstperm;
{IndexList@@freeinvb, freesym}]

When CompleteIrrDecImposeNewMethod is set to True, we don’t use the complete group G = freesym. The ideal method would be to use ImposeSym and ExpandSym with the option SmartExpand->True. This however, requires the computation of the mixed symmetry which takes a long time to compute. Instead we extract the indices of the irreducible part. If we stabilize all other points of G we get a subgroup P of H (the internal symmetry). All permutations in P will give expressions that canonicalizes to the same expression. We can therefore use a transversal of P in G instead of the complete group G.

FreeIndsOfIrrDecPart[expr_,spinmetrics_]:=With[{vbPMQs=xAct`xTensor`Private`VBundleIndexPMQ/@VBundleOfMetric/@spinmetrics},
FindFreeIndices@Evaluate[expr/.((#[A_,B_]->1)&/@spinmetrics)/.((delta[A_?#,B_?#]->1)&/@vbPMQs)]]

ImposeSymmetryOnIrrdecTerms[expr_,freeinds_, freesym_,SpinVBs_]:=Module[{spinmetrics},
spinmetrics=First/@MetricsOfVBundle/@SpinVBs;ImposeSymmetryOnIrrdecTerm[Expand@expr,spinmetrics,freeinds, freesym]]

ImposeSymmetryOnIrrdecTerm[expr_Plus,spinmetrics_,freeinds_, freesym_]:=ImposeSymmetryOnIrrdecTerm[#,spinmetrics,freeinds, freesym]&/@expr

ImposeSymmetryOnIrrdecTerm[expr_,spinmetrics_,freeinds_, freesym_]:=Module[{irrdecpartfrees,stabpoints,newfreesym,permlist},
irrdecpartfrees=FreeIndsOfIrrDecPart[expr,spinmetrics];
stabpoints=Select[Range@Length@freeinds,FreeQ[irrdecpartfrees,freeinds[[#]],1]&];
newfreesym=BaseChange2[freesym,stabpoints];
permlist=InversePerm/@xAct`SymManipulator`Private`TransversalComputation[StrongGenSet[newfreesym[[1]],Stabilizer[stabpoints,newfreesym[[2]]]],newfreesym];
Expand@ImposeSymmetry[expr,freeinds,Group@@permlist,Identity]
]

6.2. IndepIrrDecompose

DgSpinVBOfSpinVB[spin_]:=If[FreeQ[ServantsOf[spin],Dagger@spin],spin,Dagger@spin]

RemoveConstantFactors[c_?ConstantQ*expr_]:=expr;
RemoveConstantFactors[expr_]:=expr;

FindSpinorsInIrrDec[expr_,spinmetrics_]:=With[{vbPMQs=xAct`xTensor`Private`VBundleIndexPMQ/@VBundleOfMetric/@spinmetrics,termlist=xAct`xTensor`Private`ListOfTerms[expr],spindgs=Union[DgSpinVBOfSpinVB/@VBundleOfMetric/@spinmetrics]},
Union[ChangeFreeIndicesToDefault[RemoveConstantFactors[(#/.((#[A_,B_]->1)&/@spinmetrics)/.((delta[A_?#,B_?#]->1)&/@vbPMQs))],spindgs]&/@termlist]]

Options[IndepIrrDecompose]^={Verbose->False,IrrDecSpinBundles->All, UseSym->True,IrrDecImposeNewMethod->False,TimeVerbose->False, ResultType->"Expression"};

IndepIrrDecompose[expr_,options:OptionsPattern[]]:=Module[{irrdec,irrdecparts,AllIrrDecSpinBundles,SpinVBs,spinbundlesoption,SpinMetrics,valences,duplicatevalences,possilbydep,irrdecrule,relations,reductioneqs,indepirrdec,resulttype,verb},
{spinbundlesoption, resulttype,verb}=OptionValue[{IrrDecSpinBundles,ResultType,Verbose}];
AllIrrDecSpinBundles=xAct`Spinors`VBundleOfSolderingForm/@xAct`Spinors`$SolderingForms;
AllIrrDecSpinBundles=Union[AllIrrDecSpinBundles,Dagger/@AllIrrDecSpinBundles];
SpinVBs=AllIrrDecSpinBundles;
If[Not[spinbundlesoption===All],
SpinVBs=Intersection[AllIrrDecSpinBundles,spinbundlesoption]];
SpinMetrics=First/@MetricsOfVBundle/@SpinVBs;
If[verb,Print["Computing irreducible decomposition."]];
irrdec=IrrDecompose[expr,ResultType->"Expression",FilterRules[{options},(Options[IrrDecompose]/.ResultType->""/.Verbose->False)]];
If[verb,Print["Finding irrducible parts."]];
irrdecparts=FindSpinorsInIrrDec[irrdec,SpinMetrics];
If[verb,Print["Found ",Length@irrdecparts ," irreducible parts."]];
valences=Tally/@Map[VBundleOfIndex,List@@@(FindFreeIndices/@irrdecparts),{2}];
duplicatevalences=Select[Tally@valences,Last[#]>1&];
possilbydep=Function[x,Pick[irrdecparts,#===x&/@valences]]/@First/@duplicatevalences;
If[verb,Print["Possibly dependent ", possilbydep]];
irrdecrule=MakeCompareRule[Evaluate[{expr,irrdec}]];
relations=Map[#==ScreenDollarIndices@ToCanonical@ImposeSym@ToCanonical@ContractMetric[SeparateMetric[]@RemoveSym@#/.irrdecrule]&,possilbydep,{2}];
reductioneqs=Flatten[Solve/@relations]/.Rule->Equal;
If[verb,Print["Found relations ", reductioneqs]];
indepirrdec=Expand[irrdec/.Flatten[EqToCompareRule/@reductioneqs]];
Switch[resulttype,
"Expression",indepirrdec,
"Equation",Equal[expr,indepirrdec],
"Rule",MakeRule[Evaluate[{expr,indepirrdec},MetricOn->All,ContractMetrics->True]],
"CompareRule",MakeCompareRule[Evaluate[{expr,indepirrdec}]],
_,Null]
]

6.3. Naming irreducible parts

NameTensors1[struct_List,basename_String,baseprintas_String,spin_]:=If[Length@struct==1,{{Symbol[StringJoin[basename,ToString[struct[[1,1]]],ToString[struct[[1,2]]]]]@@Join[DownIndex/@GetIndexRange[struct[[1,1]],spin],GetIndexRange[struct[[1,2]],Dagger@spin]],struct[[1,3]],StringJoin["\!\(\*UnderscriptBox[\(",baseprintas,"\), \(",ToString[struct[[1,1]]],",", ToString[struct[[1,2]]],"\)]\)"],StringJoin["\\underset{",ToString[struct[[1,1]]],",", ToString[struct[[1,2]]],"}{",xAct`TexAct`Tex[baseprintas],"}{}"]}},
{Symbol[StringJoin[basename,ToString[struct[[#,1]]],ToString[struct[[#,2]]],FromCharacterCode[96+#]]]@@Join[DownIndex/@GetIndexRange[struct[[#,1]],spin],GetIndexRange[struct[[#,2]],Dagger@spin]],struct[[#,3]],StringJoin["\!\(\*UnderscriptBox[\(",baseprintas,FromCharacterCode[96+#],"\), \(",ToString[struct[[#,1]]],",", ToString[struct[[#,2]]],"\)]\)"],StringJoin["\\underset{",ToString[struct[[#,1]]],",", ToString[struct[[#,2]]],"}{",xAct`TexAct`Tex[baseprintas],FromCharacterCode[96+#],"}{}"]}&/@Range@Length@struct
]

GiveNamesToIrreducibleParts[irrdectensors_List,basename_String,baseprintas_String,spin_]:=With[{NumOfUnprimedVBundles=Length@Cases[UpIndex/@#,spin]&,
NumOfPrimedVBundles=Length@Cases[UpIndex/@#,Dagger[spin]]&},
Module[{sortedtensors=
Sort[{NumOfUnprimedVBundles[VBundleOfIndex/@FindFreeIndices@#],NumOfPrimedVBundles[VBundleOfIndex/@FindFreeIndices@#],#}&/@irrdectensors],numberofkind},
Flatten[NameTensors1[#,basename,baseprintas,spin]&/@GatherBy[sortedtensors,Most],1]
]]

DefIrreducibleSpinors[irrdectensors_List,basename_String,baseprintas_String,spin_]:=Module[{namestruct=GiveNamesToIrreducibleParts[irrdectensors,basename,baseprintas,spin],defcommands},defcommands={First[#],BaseOfVBundle[spin],CompatibleSymmetric[IndexList@@First[#]],PrintAs->#[[3]]}&/@namestruct;
xAct`Spinors`DefSpinor@@@defcommands;
(xAct`TexAct`Tex[Head[First[#]]]^=#[[4]])&/@namestruct;
Equal@@@(Take[#,2]&/@namestruct)]

7. More tools

7.1. SymmetryOfExpression

This function computed the symmetry of the free indices of any expression by a conducting a tree search. At the nodes the symmetry is tested with HasSym, which checks if a permutation is in the symmetry by checking if the permuted expression minus the original expression canonicalizes to 0. Signed permutations are handled as a permutation of two extra points as in xPerm.

SymmetryOfExpression[expr_]:=Module[{frees=FindFreeIndices@expr,HasSym,completesymmetricSGS},
With[{n=Length@frees},
HasSym[p_]:=ToCanonical[PermuteIndices[expr,frees,xAct`xPerm`Private`ToSign[p,n]]-expr]===0;
completesymmetricSGS=JoinSGS[Symmetric[Range@n,Cycles],Symmetric[{n+1,n+2},Cycles]];
completesymmetricSGS=ReplacePart[completesymmetricSGS,1->Append[Range[n-1],n+1]];
{xAct`xPerm`Private`ToSign[Search[completesymmetricSGS,HasSym,1,StrongGenSet[{},GenSet[]]],n],frees}]]

7.2. MakeCompareRule

MakeCompare rule is an alternative to MakeRule. Instead of generating all symmetry equivalent left hand sides it makes a general rule that decides at application time if the expression can be rewritten as the left hand side of the rule.
The comparison is made with a modified version of EqualExpressionsQ.

SpecialEqualExpressionsQ[0,expr2_]:={expr2===0,1,{},{}};

SpecialEqualExpressionsQ[expr1_,expr2_]:=
    Module[{inds1,frees1,frees2,dummies1,dummies2,sortedfrees1,sortedfrees2, sortrule1,dummies1slots, frees1slots,slotrules1, expr1SGS, stabilizerSGS, compatiblesymfrees1, transversal1, indstoslotnumbers,freesym1,perms1,arrexpr2,newexpr2,
        count,result=False,perm={},xxx,equation,solution=0},
      (* Indices of expr1 *)
(* Print["SpecialEqualExpressionsQ[",expr1,",",expr2,"]"];*)
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];
sortedfrees2=List@@(Last/@Sort[{VBundleOfIndex[#],#}&/@IndexSort[frees2]]);
sortedfrees1=List@@(Last/@Sort[{VBundleOfIndex[#],#}&/@IndexSort[frees1]]);
If[And[Length[frees1]==Length[frees2],
(VBundleOfIndex/@sortedfrees1)===(VBundleOfIndex/@sortedfrees2)],
dummies2=List@@dummies2;
inds1=List@@inds1;
{slotrules1,expr1SGS}=List@@Take[SymmetryOf@expr1,-2];
indstoslotnumbers=(Reverse/@slotrules1)/.xAct`xTensor`Private`slot[x_]->x;
sortrule1=Thread@Rule[Range@Length@frees1,sortedfrees1/.indstoslotnumbers];
dummies1slots=dummies1/.indstoslotnumbers;
frees1slots=frees1/.indstoslotnumbers;
xAct`xTensor`Private`checkChangeFreeIndices[sortedfrees1,sortedfrees2];
(* Computing permutations of the free indices of expr1 *)
(* One could 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. *)
stabilizerSGS=StabilizerSGS[dummies1slots,expr1SGS]/.Reverse/@sortrule1;
compatiblesymfrees1=CompatibleSymmetric[IndexList@@sortedfrees1];
(* If the stabilizer is trivial, we will have all permutations. *)
If[stabilizerSGS[[2]]===GenSet[],
transversal1 =List@@Dimino@compatiblesymfrees1,
(* The code for the transversal is probably not optimal. *)
transversal1 =TransversalInCompatibleSymmGroup[stabilizerSGS,compatiblesymfrees1];
];
(* 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,sortedfrees1];
(* 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[HoldPattern@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]},ReplaceIndex[Evaluate[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[HoldPattern@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]},ReplaceIndex[Evaluate[sign*arrexpr2],perm]]/;equalexpr]]]]]];

For a derivative of a single tensor:

MakeCompareRule[{covd_?CovDQ[indscdeq___][tensora_?xTensorQ[indsaeq___]],RHS_},___]:=With[{dummieseq=xAct`xTensor`Private`TakeEPairs[{indscdeq,indsaeq}]},With[{numdummiesaeq=Length@Intersection[dummieseq,{indscdeq,indsaeq}]},RuleDelayed[HoldPattern@covd[indscd___][tensora[indsa___]],With[{dummies=xAct`xTensor`Private`TakeEPairs[{indscd,indsa}],frees=xAct`xTensor`Private`TakeFrees[{indscd,indsa}]},Module[{equalexpr,sign,sortedfrees,perm},{equalexpr,sign,sortedfrees,perm}=SpecialEqualExpressionsQ[covd[indscdeq][tensora[indsaeq]],covd[indscd][tensora[indsa]]];With[{arrexpr2=ReplaceDummies[RHS]},ReplaceIndex[Evaluate[sign*arrexpr2],perm]]/;equalexpr]/;Length@Intersection[dummies,{indscd,indsa}]==numdummiesaeq]]]];

7.3. Giving indices to tensors

Get a range of k indices in vbundle.

GetIndexRange[k_,vbundle_]:=Module[{n=Length@Flatten[IndicesOfVBundle[vbundle]]},NewIndexIn[vbundle]&/@Range[k-n];Take[Flatten[IndicesOfVBundle[vbundle]],k]]

Give indices to a tensor. The second argument is a list of vbundles which should have up indices. The default is down indices for all vbundles.

GiveIndicesToTensor[TT_?xTensorQ,upvbundles_:{}]:=Module[{n=Length@SlotsOfTensor@TT,vbundles=UpIndex/@SlotsOfTensor@TT,gatheredslots,inds,slotrules},gatheredslots=GatherBy[Range@n,vbundles[[#]]&];
inds=Flatten[GetIndexRange@@@({Length@#,vbundles[[First@#]]}&/@gatheredslots)];
inds=If[MemberQ[upvbundles,VBundleOfIndex[#]],#,DownIndex[#]]&/@inds;
slotrules=Thread@Rule[Flatten@gatheredslots,inds];
TT@@(Range@n/.slotrules)]

ChangeFreeIndicesToDefault[expr_,upvbundles_:{}]:=
With[{freeinds=List@@FindFreeIndices@Evaluate[expr]},Module[{n=Length@freeinds,vbundles=UpIndex/@VBundleOfIndex/@freeinds,gatheredslots,inds,slotrules,indsrule},gatheredslots=GatherBy[Range@n,vbundles[[#]]&];
inds=Flatten[GetIndexRange@@@({Length@#,vbundles[[First@#]]}&/@gatheredslots)];
inds=If[MemberQ[upvbundles,VBundleOfIndex[#]],#,DownIndex[#]]&/@inds;
slotrules=Thread@Rule[Flatten@gatheredslots,inds];
indsrule=Thread@Rule[freeinds,Range@n/.slotrules];
ScreenDollarIndices[ReplaceIndex[Evaluate[ReplaceDummies@expr],indsrule]]]];

PlaceIndicesInTensor[TT_?xTensorQ,{inds___}]:=Module[{gatheredinds=Prepend[#,(VBundleOfIndex@First@#)]&/@GatherBy[{inds},VBundleOfIndex],pos,slotrules},
pos=Flatten[Position[DownIndex@SlotsOfTensor@TT,#,1]]&/@First/@gatheredinds;
If[(Length/@pos)=!=(Length/@(Rest/@gatheredinds)),Throw[Message[General::error,StringJoin["PlaceIndicesInTensor: Indices does not fit in tensor: ",ToString@TT,", ", ToString[{inds}]]]]];
slotrules=Flatten[MapThread[Thread[Rule[#1,#2]]&,{pos,Rest/@gatheredinds}]];
TT@@(Range@Length@SlotsOfTensor[TT]/.slotrules)
]

PlaceIndicesInSpinor[TT_?xTensorQ,{inds___},{spin_,spindg_}]:=With[{slotsTT=SlotsOfTensor@TT,unprimedinds=Select[{inds},VBundleOfIndex[#]==spin&],primedinds=Select[{inds},VBundleOfIndex[#]==spindg&]},
With[{unprimedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spin)&],primedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spindg)&]},TT[Sequence@@(Range@Length@slotsTT/.Join[Thread[primedpos->primedinds],Thread[unprimedpos->unprimedinds],SlotRulesOtherVBs[slotsTT,{inds},spin,spindg]])]]]

PlaceIndicesInSpinor[TT_?xTensorQ,{inds___},spin_?VBundleQ]:=PlaceIndicesInSpinor[TT,{inds},{spin,Dagger@spin}]

7.4. Valence numbers

ValenceNumbersOfSpinor[T_?xTensorQ,spin_]:={NumOfUnprimedVBundles[spin][SlotsOfTensor[T]],NumOfPrimedVBundles[spin][SlotsOfTensor[T]]}

ExtractUnprimedAndPrimed[spin_]:={Select[#,xAct`xTensor`Private`VBundleIndexPMQ@spin],Select[#,xAct`xTensor`Private`VBundleIndexPMQ@Dagger@spin]}&

SymmetricOfValence[k_,l_,spin_]:=SymmetricOfValence[k,l,spin,{}];
SymmetricOfValence[k_?NumericQ,l_?NumericQ,spin_,ovb_]:=Join[ovb,ConstantArray[-spin,k],ConstantArray[-Dagger@spin,l]]/;And[k≥0,l≥0]

ExtractOtherVBundles[spin_][vbs_List]:=DeleteCases[vbs,Alternatives[spin,-spin,Dagger@spin,-Dagger@spin]]

NumOfUnprimedVBundles[spin_][SymmetricOfValence[k_,l_,spin_,ovb_]]:=k;
NumOfPrimedVBundles[spin_][SymmetricOfValence[k_,l_,spin_,ovb_]]:=l;
ExtractOtherVBundles[spin_][SymmetricOfValence[k_,l_,spin_,ovb_]]:=ovb;

NumOfUnprimedVBundles[spin_][vbs_List]:=Length@Cases[UpIndex/@vbs,spin];
NumOfPrimedVBundles[spin_][vbs_List]:=Length@Cases[UpIndex/@vbs,Dagger@spin];
NumOfUnprimedSlots[spin_]:=NumOfUnprimedVBundles[spin][SlotsOfTensor[#]]&;
NumOfPrimedSlots[spin_]:=NumOfPrimedVBundles[spin][SlotsOfTensor[#]]&;
NumOfUnprimedSlotsDyad[dyad_]:=NumOfUnprimedSlots[xAct`xCoba`VBundleOfBasis@dyad];
NumOfPrimedSlotsDyad[dyad_]:=NumOfPrimedSlots[xAct`xCoba`VBundleOfBasis@dyad];

SymmetricSpinorOfValenceQ[TT_,k_,l_,spin_]:=And[xTensorQ[TT],
With[{slots=SlotsOfTensor[TT]},If[Head[slots]===SymmetricOfValence,True,CompatibleSymQ[slots,SymmetryGroupOfTensor@TT,{spin,Dagger@spin}]]],
ValenceNumbersOfSpinor[TT,spin]==={k,l}];

SymmetricSpinorOfArbitraryValenceQ[spin_][TT_?xTensorQ]:=SymmetricSpinorOfArbitraryValenceQ[spin][TT]=With[{k=NumOfUnprimedSlots[spin][TT],l=NumOfPrimedSlots[spin][TT]},SymmetricSpinorOfValenceQ[TT,k,l,spin]]

EnoughIndsQSpin=Function[{TTT,k,l,spin},And[Simplify[NumOfUnprimedSlots[spin][TTT]≥k],Simplify[NumOfPrimedSlots[spin][TTT]≥l]]];

ZeroTensorOfTensor[TT_]:=ZeroTensor@Block[{xAct`SymManipulator`Private`$SlotsOfTensorWarning=False},SlotsOfTensor@TT]

7.5. Finding tensor coefficients

It is often important to be able to find the coefficient of a tensor. The following idea works only if the complete expression does not have any free indices. However, if the expression has free indices, these can usually be contracted with an arbitrary tensor. It is also important that the tensor we are interested in only appears once per term. The tensor we are interested in can be a single tensor or a product, with or without contractions. Observe that different contractions will be treated as different tensors unless they canonicalize to the same expression.
The idea here is to change xSortPrecedence for the tensor we are interested in. After canonicalization, this tensor will then have the same index configuration in all terms where it appears. Unfortunately it might be difficult to know which index configuration will be due to possible dummy indices, therefore we construct some rules that looks for the tensor or product of tensors in question and records the actual index configuration it found in $IndexConfigs. The tensor coefficients can then be found by the Mathematica function Coefficient.

$IndexConfigs={};

AddToIndexConfigRule[tensora_?xTensorQ[indsaeq___]]:=With[{dummieseq=xAct`xTensor`Private`TakeEPairs[{indsaeq}]},With[{numdummiesaeq=Length@Intersection[dummieseq,{indsaeq}]},RuleDelayed[HoldPattern@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]];
If[equalexpr,AppendTo[$IndexConfigs,tensora[indsa]]];
tensora[indsa]/;equalexpr]/;Length@Intersection[dummies,{indsa}]==numdummiesaeq]]]];

AddToIndexConfigRule[LHS:Times[_?xTensorQ[___],_?xTensorQ[___]...]]:=
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[HoldPattern@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];
If[equalexpr,AppendTo[$IndexConfigs,rulelhs]];
rulelhs/;equalexpr]]]]]];

FullSymmetryOfFreeIndices[expr_]:=Module[{freeinds=FindFreeIndices@expr,abstractfree,actractfreeslots, restslots, slotexpr,slotrules,reverseslotrules,sym,n,freefirstperm, freesym,restfrees, tmpsym,inds,HasSym,completesymmetricSGS},
{n,slotexpr,slotrules,sym}=List@@SymmetryOf[expr];
If[And[OrderOfGroup[sym]==1, Length@freeinds==n],{StrongGenSet[{},GenSet[]],freeinds},
inds=(xAct`xTensor`Private`slot/@Range[n])/.slotrules;
reverseslotrules=Reverse/@slotrules;
abstractfree=Select[freeinds,AIndexQ];
actractfreeslots=First/@List@@Replace[List@@abstractfree,reverseslotrules , {1}];
restslots=Complement[Range@n,actractfreeslots];
(* Find a permutation that moves the free indices first. *)
freefirstperm=List@@TranslatePerm[PermutationFromTo[Range@n,Join[actractfreeslots,restslots]],Rules];
(* Change the base of sym so it begins with the dummy indices. *)
tmpsym=StabilizerSGS[restslots,sym]/.freefirstperm;
n=Length@freeinds;
If[CompatibleSymQ[List@@(UpIndex/@VBundleOfIndex/@freeinds),tmpsym],
freesym=tmpsym,
HasSym[p_]:=ToCanonical[PermuteIndices[expr,freeinds,xAct`xPerm`Private`ToSign[p,n]]-expr]===0;
completesymmetricSGS=JoinSGS[CompatibleSymmetric[freeinds],Symmetric[{n+1,n+2},Cycles]];
completesymmetricSGS=ReplacePart[completesymmetricSGS,1->Append[Range[n-1],n+1]];
freesym=xAct`xPerm`Private`ToSign[Search[completesymmetricSGS,HasSym,1,tmpsym],n]];
{freesym,freeinds}]]

FindTensorHeads[TT_?xTensorQ[inds___]]:={TT}

FindTensorHeads[expr_Times]:=FindTensorHeads/@(List@@expr)

FindTensorCoefficients[expr_?ScalarQ,tensors_List]:=Module[{AddToIndexConfigR1=AddToIndexConfigRule/@tensors,TensorHeads,canonicalizedexpr,tensorsasinexpr,ImposeCorrectSym, oldxsortprecendence, coefflist},
TensorHeads=Flatten[FindTensorHeads/@tensors];
ImposeCorrectSym[expr1_,expr2_]:=With[{fullsym=FullSymmetryOfFreeIndices[expr2]},ImposeSym[expr1,ChangeIndex/@fullsym[[2]],fullsym[[1]]]];
oldxsortprecendence=xSortPrecedence/@TensorHeads;
(xSortPrecedence[#]^=2)&/@TensorHeads;
$IndexConfigs={};
canonicalizedexpr=ScreenDollarIndices[ToCanonical[expr]]/.AddToIndexConfigR1;
tensorsasinexpr=Union@$IndexConfigs;
(xSortPrecedence[TensorHeads[[#]]]^=oldxsortprecendence[[#]])&/@Range@Length@TensorHeads;
 coefflist={ImposeCorrectSym[Coefficient[canonicalizedexpr,#],#],#}&/@tensorsasinexpr;
Join[{{ToCanonical@ExpandSym[expr-Plus@@Times@@@coefflist],1}},coefflist]
]

FindTensorCoefficients[___]:=Throw@Message[FindTensorCoefficients::error, "Only scalar expressions are allowed."]

7.6. General equation handling tools

EqToRule[expr_Equal]:=MakeRule[Evaluate[List@@expr],MetricOn->All,ContractMetrics->True]
EqToRule[True]={};

EqToCompareRule[expr_Equal]:=MakeCompareRule[Evaluate[List@@expr]];
EqToCompareRule[True]={};

FlipEquation[LHS_==Times[a_?NumberQ,RHS1_,RHSS2___]]:=Times[RHS1,RHSS2]==(LHS)/a
FlipEquation[LHS_==RHS_]:=RHS==LHS

Isolate[LHS_==RHS_,pattern_]:=Module[{eq=LHS-RHS,pos},
pos=Position[eq,pattern];
If[Length@pos==0,Throw@Message[Isolate::error, "Could not find "<>ToString@pattern]];
pos=First@First@pos;
FlipEquation[eq[[pos]]-eq==eq[[pos]]]]

8. Fundamental spinor operators

8.1. Common functions for the fundamental spinor operators

DivName[covd_?xAct`Spinors`SpinCovDQ]:=Symbol[StringJoin["Div",SymbolName[covd]]]; CurlName[covd_?xAct`Spinors`SpinCovDQ]:=Symbol[StringJoin["Curl",SymbolName[covd]]]; CurlDgName[covd_?xAct`Spinors`SpinCovDQ]:=Symbol[StringJoin["CurlDg",SymbolName[covd]]]; TwistName[covd_?xAct`Spinors`SpinCovDQ]:=Symbol[StringJoin["Twist",SymbolName[covd]]];
BoxName[covd_?xAct`Spinors`SpinCovDQ]:=Symbol[StringJoin["Box",SymbolName[covd]]];

ExtractUnprimedAndPrimedIndsGivenCovD[inds_List,covd_?xAct`Spinors`SpinCovDQ]:=With[{spin=VBundlesOfCovD[covd][[2]],spindg=VBundlesOfCovD[covd][[3]]},{Select[inds,xAct`xTensor`Private`VBundleIndexPMQ@spin],Select[inds,xAct`xTensor`Private`VBundleIndexPMQ@spindg],Select[inds,Not@Or[(xAct`xTensor`Private`VBundleIndexPMQ@spin)[#],(xAct`xTensor`Private`VBundleIndexPMQ@spindg)[#]]&]}];

SpatialSpinCovDQ[_]:=False;

SpacetimeSpinCovDQ=And[xAct`Spinors`SpinCovDQ[#],Not@SpatialSpinCovDQ[#]]&;

8.2. Definition of fundamental spinor operators and tools

FundSpinOpQ[_]:=False;
CommuteFundSpinOp[__]:={};

$FundSpinOpCovDs={};
$FundSpinVBundles={};

SymFromOtherVBs[TT_,spin_,spindg_]:=Module[{slotsTT=SlotsOfTensor@TT,spinpositions,oVBpositions},
spinpositions=Flatten@Select[GatherBy[Range@Length@slotsTT,(UpIndex/@slotsTT)[[#]]&],MemberQ[{spin,spindg},UpIndex@slotsTT[[First[#]]]]&];
oVBpositions=Complement[Range@Length@slotsTT,spinpositions];
StabilizerSGS[spinpositions,TranslatePerm[SymmetryGroupOfTensor@TT,Cycles]]/.Thread@Rule[oVBpositions,Range@Length@oVBpositions]
]

SlotRulesOtherVBs[slotsTT_,inds_,spin_,spindg_]:=Module[{gatheredinds=Prepend[#,(xAct`xTensor`Private`SignedVBundleOfIndex@First@#)]&/@GatherBy[inds,xAct`xTensor`Private`SignedVBundleOfIndex],otherinds,otherpos},
otherinds=Select[gatheredinds,And[First@#=!=-spin,First@#=!=-spindg,First@#=!=spin,First@#=!=spindg]&];
otherpos=Flatten[Position[slotsTT,#,1]]&/@First/@otherinds;
Flatten[MapThread[Thread[Rule[#1,#2]]&,{otherpos,Rest/@otherinds}]]
]

Options[DefFundSpinOperators]^={HideCovDSymbolReference->True,ShowValenceInfo->True};

SetAttributes[DefFundSpinOperators,ReadProtected]

Definition of the divergence, curl, curldagger and twist operators, and most tools that are needed to handle them.

$SlotsOfTensorWarning=True;

OpValenceNumberPrintFunc[TT_,opsymbol_,opsymbsuper_, covdsymb_,spin_]:=Join[{If[OptionValue[DefFundSpinOperators,ShowValenceInfo],SubsuperscriptBox[opsymbol,StringJoin[ToString[DisplayForm[NumOfUnprimedSlots[spin][TT]],StandardForm],",",ToString[DisplayForm[NumOfPrimedSlots[spin][TT]],StandardForm]],opsymbsuper],
If[StringLength[opsymbsuper]==0,opsymbol,SuperscriptBox[opsymbol,opsymbsuper]]]},If[OptionValue[DefFundSpinOperators,HideCovDSymbolReference],{},{"[",covdsymb,"]"}]]

TexValenceNumbers[TT_,spin_]:=If[OptionValue[DefFundSpinOperators,ShowValenceInfo],StringJoin["_{",xAct`TexAct`Tex[NumOfUnprimedSlots[spin][TT]/.xAct`TexAct`ToOrderedPlus],",",xAct`TexAct`Tex[NumOfPrimedSlots[spin][TT]/.xAct`TexAct`ToOrderedPlus],"} "]," "];

DefFundSpinOperators[covd_?SpacetimeSpinCovDQ,options:OptionsPattern[]]:=
If[Not@And[MemberQ[$ContextPath,"xAct`Spinors`"],MemberQ[$ContextPath,"xAct`TexAct`"]],Throw@Message[DefFundSpinOperators::error,"Please load the packages Spinors and TexAct first."],
If[Not@xAct`Spinors`SpinCovDQ@covd,Throw@Message[DefFundSpinOperators::error,"Not a SpinCovD."],
With[{
(* Names of the operators *)
divcd=DivName@covd, curlcd=CurlName@covd, curldgcd=CurlDgName@covd, twistcd=TwistName@covd,boxcd=BoxName@covd,
(* Names of lists *)
$irrdecrulescd=Symbol[StringJoin["$IrrDecRules",SymbolName[covd]]],
$expandfundspinoprulescd=Symbol[StringJoin["$ExpandFundSpinOpRules",SymbolName[covd]]],
(* Extract the spin bundles *)
spin=VBundlesOfCovD[covd][[2]],spindg=VBundlesOfCovD[covd][[3]]},
If[Not[Dagger@spin==spindg],Throw@Message[DefFundSpinOperators::error,"VBundles are not conjugated."]];
With[{
(* The spin metrics *)
eps=First@MetricsOfVBundle@spin,
epsdg=First@MetricsOfVBundle@spindg,
(* Functions for extraction of unprimed and primed indices and their numbers. *)
ExtractOtherInds=Select[#,Not@Or[(xAct`xTensor`Private`VBundleIndexPMQ@spin)[#],(xAct`xTensor`Private`VBundleIndexPMQ@spindg)[#]]&]&,
EnoughIndsQ=Function[{TTT,k,l},And[Length@Cases[UpIndex/@SlotsOfTensor@TTT,spin]≥k,Length@Cases[UpIndex/@SlotsOfTensor@TTT,spindg]≥l]]},
Module[{expanddivcd,expandcurlcd,expandcurldgcd,expandtwistcd},
(* Add to the list $FundSpinOpCovDs and $FundSpinVBundles *)
AppendTo[$FundSpinOpCovDs,covd];
$FundSpinVBundles=Union[{spin,spindg},$FundSpinVBundles];

(* FundSpinOpQ *)

FundSpinOpQ[divcd]^=True;
FundSpinOpQ[curlcd]^=True;
FundSpinOpQ[curldgcd]^=True;
FundSpinOpQ[twistcd]^=True;

(* Linking *)
CovDOfFundSpinOp[divcd]^=covd;
CovDOfFundSpinOp[curlcd]^=covd;
CovDOfFundSpinOp[curldgcd]^=covd;
CovDOfFundSpinOp[twistcd]^=covd;

(* Tensor definitions. *)
xTensorQ[divcd[TT_?xTensorQ]]^= True;
xTensorQ[curlcd[TT_?xTensorQ]]^= True;
xTensorQ[curldgcd[TT_?xTensorQ]]^= True;
xTensorQ[twistcd[TT_?xTensorQ]]^= True;

(* Links to other symbols *)
VisitorsOf[covd]^=Union[VisitorsOf[covd],{divcd,curlcd,curldgcd,twistcd}];

(* How to undefine *)
Undef[divcd]^:=UndefFundSpinOperators@covd;
Undef[curlcd]^:=UndefFundSpinOperators@covd;
Undef[curldgcd]^:=UndefFundSpinOperators@covd;
Undef[twistcd]^:=UndefFundSpinOperators@covd;

(* Usage messages *)
$expandfundspinoprulescd::usage=StringJoin@@(ToString/@{$expandfundspinoprulescd,"contains rules for expansion of the fundamental operators ", divcd, ", ", curlcd,", ", curldgcd," and ", twistcd," of expr in terms of symmetrized derivatives."});
$irrdecrulescd::usage=StringJoin@@(ToString/@{$irrdecrulescd," is a list of rules that are used by the function ", ToFundSpinOp," to convert derivatives to fundamental operators. Both the functions ", AddToIrrdecRules," and ",ToFundSpinOp ," adds rules to this list."});

(* Calculation of the slots of the tensors. *)
SlotsOfTensor[divcd[TT_?xTensorQ]]^:=Module[{k,l,slotsTT=SlotsOfTensor@TT,ovb},
k=NumOfUnprimedVBundles[spin]@slotsTT;
l=NumOfPrimedVBundles[spin]@slotsTT;
ovb=ExtractOtherVBundles[spin]@slotsTT;
If[And[$SlotsOfTensorWarning,Or[k<1, l<1]],Throw@Message[divcd::error,StringJoin["Wrong VBundles: ", ToString[slotsTT]]],
SymmetricOfValence[k-1,l-1,spin,ovb]]];
SlotsOfTensor[curlcd[TT_?xTensorQ]]^:=Module[{k,l,slotsTT=SlotsOfTensor@TT,ovb},
k=NumOfUnprimedVBundles[spin]@slotsTT;
l=NumOfPrimedVBundles[spin]@slotsTT;
ovb=ExtractOtherVBundles[spin]@slotsTT;
If[And[$SlotsOfTensorWarning,Or[k<0, l<1]],Throw@Message[curlcd::error,StringJoin["Wrong VBundles: ", ToString[slotsTT]]],
SymmetricOfValence[k+1,l-1,spin,ovb]]];
SlotsOfTensor[curldgcd[TT_?xTensorQ]]^:=Module[{k,l,slotsTT=SlotsOfTensor@TT,ovb},
k=NumOfUnprimedVBundles[spin]@slotsTT;
l=NumOfPrimedVBundles[spin]@slotsTT;
ovb=ExtractOtherVBundles[spin]@slotsTT;
If[And[$SlotsOfTensorWarning,Or[k<1, l<0]],Throw@Message[curldgcd::error,StringJoin["Wrong VBundles: ", ToString[slotsTT]]],
SymmetricOfValence[k-1,l+1,spin,ovb]]];
SlotsOfTensor[twistcd[TT_?xTensorQ]]^:=Module[{k,l,slotsTT=SlotsOfTensor@TT,ovb},
k=NumOfUnprimedVBundles[spin]@slotsTT;
l=NumOfPrimedVBundles[spin]@slotsTT;
ovb=ExtractOtherVBundles[spin]@slotsTT;
If[And[$SlotsOfTensorWarning,Or[k<0, l<0]],Throw@Message[twistcd::error,StringJoin["Wrong VBundles: ", ToString[slotsTT]]],
SymmetricOfValence[k+1,l+1,spin,ovb]]];

(* Quicker calculation of the number of slots *)
TagSetDelayed[divcd,NumOfUnprimedSlots[spin][divcd[TT_?xTensorQ]],NumOfUnprimedSlots[spin][TT]-1];
TagSetDelayed[divcd,NumOfPrimedSlots[spin][divcd[TT_?xTensorQ]],NumOfPrimedSlots[spin][TT]-1];
TagSetDelayed[curlcd,NumOfUnprimedSlots[spin][curlcd[TT_?xTensorQ]],NumOfUnprimedSlots[spin][TT]+1];
TagSetDelayed[curlcd,NumOfPrimedSlots[spin][curlcd[TT_?xTensorQ]],NumOfPrimedSlots[spin][TT]-1];
TagSetDelayed[curldgcd,NumOfUnprimedSlots[spin][curldgcd[TT_?xTensorQ]],NumOfUnprimedSlots[spin][TT]-1];
TagSetDelayed[curldgcd,NumOfPrimedSlots[spin][curldgcd[TT_?xTensorQ]],NumOfPrimedSlots[spin][TT]+1];
TagSetDelayed[twistcd,NumOfUnprimedSlots[spin][twistcd[TT_?xTensorQ]],NumOfUnprimedSlots[spin][TT]+1];
TagSetDelayed[twistcd,NumOfPrimedSlots[spin][twistcd[TT_?xTensorQ]],NumOfPrimedSlots[spin][TT]+1];

(* Calculation of the symmetry groups of the tensors *)

SymmetryGroupOfTensor[divcd[TT_?xTensorQ]]^:=If[Length@ExtractOtherVBundles[spin]@SlotsOfTensor@TT≤1, CompatibleSym@SlotsOfTensor[divcd[TT]],
JoinSGS[SymFromOtherVBs[TT,spin,spindg],CompatibleSym[SlotsOfTensor[divcd[TT]],{spin,spindg}]]];
SymmetryGroupOfTensor[curlcd[TT_?xTensorQ]]^:=If[Length@ExtractOtherVBundles[spin]@SlotsOfTensor@TT≤1, CompatibleSym@SlotsOfTensor[curlcd[TT]],
JoinSGS[SymFromOtherVBs[TT,spin,spindg],CompatibleSym[SlotsOfTensor[curlcd[TT]],{spin,spindg}]]];
SymmetryGroupOfTensor[curldgcd[TT_?xTensorQ]]^:=If[Length@ExtractOtherVBundles[spin]@SlotsOfTensor@TT≤1, CompatibleSym@SlotsOfTensor[curldgcd[TT]],
JoinSGS[SymFromOtherVBs[TT,spin,spindg],CompatibleSym[SlotsOfTensor[curldgcd[TT]],{spin,spindg}]]];
SymmetryGroupOfTensor[twistcd[TT_?xTensorQ]]^:=If[Length@ExtractOtherVBundles[spin]@SlotsOfTensor@TT≤1, CompatibleSym@SlotsOfTensor[twistcd[TT]],
JoinSGS[SymFromOtherVBs[TT,spin,spindg],CompatibleSym[SlotsOfTensor[twistcd[TT]],{spin,spindg}]]];

(* Quicker check if it is symmetric *)
TagSetDelayed[divcd,SymmetricSpinorOfArbitraryValenceQ[spin][divcd[TT_?xTensorQ]],SymmetricSpinorOfArbitraryValenceQ[spin][TT]];
TagSetDelayed[curlcd,SymmetricSpinorOfArbitraryValenceQ[spin][curlcd[TT_?xTensorQ]],SymmetricSpinorOfArbitraryValenceQ[spin][TT]];
TagSetDelayed[curldgcd,SymmetricSpinorOfArbitraryValenceQ[spin][curldgcd[TT_?xTensorQ]],SymmetricSpinorOfArbitraryValenceQ[spin][TT]];
TagSetDelayed[twistcd,SymmetricSpinorOfArbitraryValenceQ[spin][twistcd[TT_?xTensorQ]],SymmetricSpinorOfArbitraryValenceQ[spin][TT]];

(* Interaction with the ZeroTensor *)
divcd[ZeroTensor[vbundles_]]:=With[{k=NumOfUnprimedVBundles[spin]@vbundles,l=NumOfPrimedVBundles[spin]@vbundles,ovb=ExtractOtherVBundles[spin]@vbundles},
ZeroTensor[SymmetricOfValence[k-1,l-1,spin,ovb]]];
curlcd[ZeroTensor[vbundles_]]:=With[{k=NumOfUnprimedVBundles[spin]@vbundles,l=NumOfPrimedVBundles[spin]@vbundles,ovb=ExtractOtherVBundles[spin]@vbundles},
ZeroTensor[SymmetricOfValence[k+1,l-1,spin,ovb]]];
curldgcd[ZeroTensor[vbundles_]]:=With[{k=NumOfUnprimedVBundles[spin]@vbundles,l=NumOfPrimedVBundles[spin]@vbundles,ovb=ExtractOtherVBundles[spin]@vbundles},
ZeroTensor[SymmetricOfValence[k-1,l+1,spin,ovb]]];
twistcd[ZeroTensor[vbundles_]]:=With[{k=NumOfUnprimedVBundles[spin]@vbundles,l=NumOfPrimedVBundles[spin]@vbundles,ovb=ExtractOtherVBundles[spin]@vbundles},
ZeroTensor[SymmetricOfValence[k+1,l+1,spin,ovb]]];

(* Printing *)
PrintAs[divcd[TT_?xTensorQ]]^:=RowBox@Join[OpValenceNumberPrintFunc[TT,"𝒟","",Last@SymbolOfCovD[covd],spin],{PrintAs[TT]}];
PrintAs[curlcd[TT_?xTensorQ]]^:=RowBox@Join[OpValenceNumberPrintFunc[TT,"𝒞","",Last@SymbolOfCovD[covd],spin],{PrintAs[TT]}];
PrintAs[curldgcd[TT_?xTensorQ]]^:=RowBox@Join[OpValenceNumberPrintFunc[TT,"𝒞","†",Last@SymbolOfCovD[covd],spin],{PrintAs[TT]}];
PrintAs[twistcd[TT_?xTensorQ]]^:=RowBox@Join[OpValenceNumberPrintFunc[TT,"𝒯","",Last@SymbolOfCovD[covd],spin],{PrintAs[TT]}];

(* Tex output. *)
xAct`TexAct`$TexInitLatexPackages=DeleteDuplicates@Join[xAct`TexAct`$TexInitLatexPackages,{"{mathrsfs}","{eucal}"}];
xAct`TexAct`$TexInitLatexExtraCode=DeleteDuplicates@Join[xAct`TexAct`$TexInitLatexExtraCode,{"\\def\\sDiv{\\mathscr{D}}","\\def\\sCurl{\\mathscr{C}}","\\def\\sCurlDagger{\\mathscr{C}^\\dagger}","\\def\\sTwist{\\mathscr{T}}"}];
TexFundOp[divcd[TT_?xTensorQ]]^:=StringJoin["\\sDiv",TexValenceNumbers[TT,spin],TexFundOp[TT]];
TexFundOp[curlcd[TT_?xTensorQ]]^:=StringJoin["\\sCurl",TexValenceNumbers[TT,spin],TexFundOp[TT]];
TexFundOp[curldgcd[TT_?xTensorQ]]^:=StringJoin["\\sCurlDagger",TexValenceNumbers[TT,spin],TexFundOp[TT]];
TexFundOp[twistcd[TT_?xTensorQ]]^:=StringJoin["\\sTwist",TexValenceNumbers[TT,spin],TexFundOp[TT]];
TexFundOp[expr_]:=xAct`TexAct`Tex@expr;
xAct`TexAct`Tex[x:((divcd|curlcd|curldgcd|twistcd)[_])]:=StringJoin[xAct`TexAct`Private`TexOpen["("],TexFundOp@x,xAct`TexAct`Private`TexClose[")"]];

(* Complex conjugates. *)
HermitianQ[divcd[TT_?xTensorQ]]^:=HermitianQ[TT];
HermitianQ[twistcd[TT_?xTensorQ]]^:=HermitianQ[TT];
Dagger[divcd[TT_?xTensorQ]]^:=If[HermitianQ[TT],divcd[TT],divcd[Dagger[TT]]];
Dagger[curlcd[TT_?xTensorQ]]^:=If[HermitianQ[TT],curldgcd[TT],curldgcd[Dagger[TT]]];
Dagger[curldgcd[TT_?xTensorQ]]^:=If[HermitianQ[TT],curlcd[TT],curlcd[Dagger[TT]]];
Dagger[twistcd[TT_?xTensorQ]]^:=If[HermitianQ[TT],twistcd[TT],twistcd[Dagger[TT]]];
Dagger[divcd[TT_?xTensorQ][inds___]]^:=Dagger[divcd[TT]][Sequence@@ExtractOtherInds[DaggerIndex/@{inds}],Sequence@@(Join@@ExtractUnprimedAndPrimed[spin][DaggerIndex/@{inds}])];
Dagger[curlcd[TT_?xTensorQ][inds___]]^:=Dagger[curlcd[TT]][Sequence@@ExtractOtherInds[DaggerIndex/@{inds}],Sequence@@(Join@@ExtractUnprimedAndPrimed[spin][DaggerIndex/@{inds}])];
Dagger[curldgcd[TT_?xTensorQ][inds___]]^:=Dagger[curldgcd[TT]][Sequence@@ExtractOtherInds[DaggerIndex/@{inds}],Sequence@@(Join@@ExtractUnprimedAndPrimed[spin][DaggerIndex/@{inds}])];
Dagger[twistcd[TT_?xTensorQ][inds___]]^:=Dagger[twistcd[TT]][Sequence@@ExtractOtherInds[DaggerIndex/@{inds}],Sequence@@(Join@@ExtractUnprimedAndPrimed[spin][DaggerIndex/@{inds}])];

(* Expansion of the operators. *)
expanddivcd[divcd[TT_?xTensorQ][inds___]]:=Module[{slotsTT=SlotsOfTensor@TT,primedinds,unprimedinds,A=DummyIn@spin,Adg=DummyIn@spindg,primedpos,unprimedpos,TTinds},
{unprimedinds,primedinds}=ExtractUnprimedAndPrimed[spin][{inds}];
unprimedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spin)&];
primedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spindg)&];
TTinds=Range@Length@slotsTT/.Thread[primedpos->Append[primedinds,-Adg]]/.Thread[unprimedpos->Append[unprimedinds,-A]]/.SlotRulesOtherVBs[slotsTT,{inds},spin,spindg];
covd[A,Adg][TT@@TTinds]];
expandcurlcd[curlcd[TT_?xTensorQ][inds___]]:=Module[{slotsTT=SlotsOfTensor@TT,primedinds,unprimedinds,Adg=DummyIn@spindg,primedpos,unprimedpos,TTinds},
{unprimedinds,primedinds}=ExtractUnprimedAndPrimed[spin][{inds}];
unprimedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spin)&];
primedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spindg)&];
TTinds=Range@Length@slotsTT/.Thread[primedpos->Append[primedinds,-Adg]]/.Thread[unprimedpos->Rest@unprimedinds]/.SlotRulesOtherVBs[slotsTT,{inds},spin,spindg];
ToCanonicalSym@ImposeSym[covd[First@unprimedinds,Adg][TT@@TTinds],IndexList@@unprimedinds]];
expandcurldgcd[curldgcd[TT_?xTensorQ][inds___]]:=Module[{slotsTT=SlotsOfTensor@TT,primedinds,unprimedinds,A=DummyIn@spin,primedpos,unprimedpos,TTinds},
{unprimedinds,primedinds}=ExtractUnprimedAndPrimed[spin][{inds}];
unprimedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spin)&];
primedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spindg)&];
TTinds=Range@Length@slotsTT/.Thread[primedpos->Rest@primedinds]/.Thread[unprimedpos->Append[unprimedinds,-A]]/.SlotRulesOtherVBs[slotsTT,{inds},spin,spindg];
ToCanonicalSym@ImposeSym[covd[A,First@primedinds][TT@@TTinds],IndexList@@primedinds]];
expandtwistcd[twistcd[TT_?xTensorQ][inds___]]:=Module[{slotsTT=SlotsOfTensor@TT,primedinds,unprimedinds,primedpos,unprimedpos,TTinds},
{unprimedinds,primedinds}=ExtractUnprimedAndPrimed[spin][{inds}];
unprimedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spin)&];
primedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spindg)&];
TTinds=Range@Length@slotsTT/.Thread[primedpos->Rest@primedinds]/.Thread[unprimedpos->Rest@unprimedinds]/.SlotRulesOtherVBs[slotsTT,{inds},spin,spindg];
ToCanonicalSym@ImposeSym[covd[First@unprimedinds,First@primedinds][TT@@TTinds],IndexList@@Join[unprimedinds,primedinds]]];
$expandfundspinoprulescd={x:(divcd[TT_?xTensorQ][inds___]):>expanddivcd[x],x:(curlcd[TT_?xTensorQ][inds___]):>expandcurlcd[x],x:(curldgcd[TT_?xTensorQ][inds___]):>expandcurldgcd[x],x:(twistcd[TT_?xTensorQ][inds___]):>expandtwistcd[x]};

(* Imposing Twist *)
ImposeFundSpinOp[expr1_,twistcd, iinds_]:=
Module[{primedinds1,unprimedinds1,irrdecexpr,k,l, (* If the expression is an expanded symmetry, it is important that we do the contractions with the same index on all terms. Therefore, we sort the free indices so the contraction is done with the same index on all terms.*)
frees=Sort[List@@FindFreeIndices@Evaluate[expr1]], primedfrees, unprimedfrees,indexrules={}, intexpr},
{unprimedinds1,primedinds1}=ExtractUnprimedAndPrimed[spin][iinds];
{unprimedfrees,primedfrees}=ExtractUnprimedAndPrimed[spin][frees];
indexrules=Join[Thread@Rule[unprimedfrees,Rest@unprimedinds1],Thread@Rule[primedfrees,Rest@primedinds1]];
intexpr=covd[First@unprimedinds1,First@primedinds1][ReplaceIndex[expr1,indexrules]];
intexpr=ToFundSpinOp[intexpr,covd];
ImposeSym[intexpr,IndexList@@Join[unprimedinds1,primedinds1]]];

(* Imposing Curl *)
ImposeFundSpinOp[expr1_,curlcd, iinds_]:=
Module[{primedinds1,unprimedinds1,irrdecexpr,k,l, (* If the expression is an expanded symmetry, it is important that we do the contractions with the same index on all terms. Therefore, we sort the free indices so the contraction is done with the same index on all terms.*)
frees=Sort[List@@FindFreeIndices@Evaluate[expr1]], primedfrees, unprimedfrees,indexrules={}, intexpr,Adg=DummyIn@spindg},
{unprimedinds1,primedinds1}=ExtractUnprimedAndPrimed[spin][iinds];
{unprimedfrees,primedfrees}=ExtractUnprimedAndPrimed[spin][frees];
indexrules=Join[Thread@Rule[unprimedfrees,Rest@unprimedinds1],Thread@Rule[primedfrees,Append[primedinds1,-Adg]]];
intexpr=covd[First@unprimedinds1,Adg][ReplaceIndex[expr1,indexrules]];
intexpr=ToFundSpinOp[intexpr,covd];
ImposeSym[intexpr,IndexList@@Join[unprimedinds1,primedinds1]]];

(* Imposing CurlDg *)
ImposeFundSpinOp[expr1_,curldgcd, iinds_]:=Module[{primedinds1,unprimedinds1,irrdecexpr,k,l, (* If the expression is an expanded symmetry, it is important that we do the contractions with the same index on all terms. Therefore, we sort the free indices so the contraction is done with the same index on all terms.*)
frees=Sort[List@@FindFreeIndices@Evaluate[expr1]], primedfrees, unprimedfrees,indexrules={}, intexpr,A=DummyIn@spin},
{unprimedinds1,primedinds1}=ExtractUnprimedAndPrimed[spin][iinds];
{unprimedfrees,primedfrees}=ExtractUnprimedAndPrimed[spin][frees];
indexrules=Join[Thread@Rule[unprimedfrees,Append[unprimedinds1,-A]],Thread@Rule[primedfrees,Rest@primedinds1]];
intexpr=covd[A,First@primedinds1][ReplaceIndex[expr1,indexrules]];
intexpr=ToFundSpinOp[intexpr,covd];
ImposeSym[intexpr,IndexList@@Join[unprimedinds1,primedinds1]]];

(* Imposing Div *)
ImposeFundSpinOp[expr1_,divcd, iinds_]:=
Module[{primedinds1,unprimedinds1,irrdecexpr,k,l, (* If the expression is an expanded symmetry, it is important that we do the contractions with the same index on all terms. Therefore, we sort the free indices so the contraction is done with the same index on all terms.*)
frees=Sort[List@@FindFreeIndices@Evaluate[expr1]], primedfrees, unprimedfrees,indexrules={}, intexpr,A=DummyIn@spin,Adg=DummyIn@spindg},
{unprimedinds1,primedinds1}=ExtractUnprimedAndPrimed[spin][iinds];
{unprimedfrees,primedfrees}=ExtractUnprimedAndPrimed[spin][frees];
indexrules=Join[Thread@Rule[unprimedfrees,Append[unprimedinds1,-A]],Thread@Rule[primedfrees,Append[primedinds1,-Adg]]];
intexpr=covd[A,Adg][ReplaceIndex[expr1,indexrules]];
intexpr=ToFundSpinOp[intexpr,covd];
ImposeSym[intexpr,IndexList@@Join[unprimedinds1,primedinds1]]];

$irrdecrulescd={};

(* Commutators *)
(* DivCurl *)
CommuteFundSpinOp[divcd,curlcd]^=
divcd[curlcd[TT_?xTensorQ]][inds___]:>Module[{Bdg=DummyIn@spindg,Cdg=DummyIn@spindg,k=NumOfUnprimedSlots[spin]@TT},k/(k+1)curlcd[divcd[TT]][inds]-boxcd[-Bdg,-Cdg]@TT[inds,Bdg,Cdg]]/;EnoughIndsQ[TT,0,2];
CommuteFundSpinOp[divcd,curlcd,curlcd,divcd]^=CommuteFundSpinOp[divcd,curlcd];
(* CurlDiv *)
CommuteFundSpinOp[curlcd,divcd]^=curlcd[divcd[TT_?xTensorQ]][inds___]:>Module[{Bdg=DummyIn@spindg,Cdg=DummyIn@spindg,k=NumOfUnprimedSlots[spin]@TT},(k+1)/k*divcd[curlcd[TT]][inds]+(k+1)/k*boxcd[-Bdg,-Cdg]@TT[inds,Bdg,Cdg]]/;EnoughIndsQ[TT,1,2];
CommuteFundSpinOp[curlcd,divcd,divcd,curlcd]^=CommuteFundSpinOp[curlcd,divcd];
(* CurlTwist *)
CommuteFundSpinOp[curlcd,twistcd]^=curlcd[twistcd[TT_?xTensorQ]][inds___]:>With[{l=NumOfPrimedSlots[spin]@TT,unprimedinds=First@ExtractUnprimedAndPrimed[spin][{inds}]},
l/(l+1)twistcd[curlcd[TT]][inds]-ToCanonicalSym@ImposeSym[(boxcd@@Take[unprimedinds,2])[TT@@DeleteCases[{inds},Alternatives@@Take[unprimedinds,2]]],IndexList@@unprimedinds]];
CommuteFundSpinOp[curlcd,twistcd,twistcd,curlcd]^=CommuteFundSpinOp[curlcd,twistcd];
(* TwistCurl *)
CommuteFundSpinOp[twistcd,curlcd]^=twistcd[curlcd[TT_?xTensorQ]][inds___]:>With[{l=NumOfPrimedSlots[spin]@TT,unprimedinds=First@ExtractUnprimedAndPrimed[spin][{inds}]},(l+1)/l*curlcd[twistcd[TT]][inds]+(l+1)/l*ToCanonicalSym@ImposeSym[(boxcd@@Take[unprimedinds,2])[TT@@DeleteCases[{inds},Alternatives@@Take[unprimedinds,2]]],IndexList@@unprimedinds]]/;EnoughIndsQ[TT,0,1];
CommuteFundSpinOp[twistcd,curlcd,curlcd,twistcd]^=CommuteFundSpinOp[twistcd,curlcd];
(* DivCurlDg *)
CommuteFundSpinOp[divcd,curldgcd]^=divcd[curldgcd[TT_?xTensorQ]][inds___]:>Module[{B=DummyIn@spin,C=DummyIn@spin,l=NumOfPrimedSlots[spin]@TT,k=NumOfUnprimedSlots[spin]@TT},
l/(l+1)curldgcd[divcd[TT]][inds]-boxcd[-B,-C]@TT[Sequence@@Drop[{inds},-k-l+2],B,C,Sequence@@Take[{inds},-k-l+2]]]/;EnoughIndsQ[TT,2,0];
CommuteFundSpinOp[divcd,curldgcd,curldgcd,divcd]^=CommuteFundSpinOp[divcd,curldgcd];
(* CurlDgDiv *)
CommuteFundSpinOp[curldgcd,divcd]^=curldgcd[divcd[TT_?xTensorQ]][inds___]:>Module[{B=DummyIn@spin,C=DummyIn@spin,l=NumOfPrimedSlots[spin]@TT,k=NumOfUnprimedSlots[spin]@TT},(l+1)/l*divcd[curldgcd[TT]][inds]+(l+1)/l*boxcd[-B,-C]@TT[Sequence@@Drop[{inds},-k-l+2],B,C,Sequence@@Take[{inds},-k-l+2]]]/;EnoughIndsQ[TT,2,1];
CommuteFundSpinOp[curldgcd,divcd,divcd,curldgcd]^=CommuteFundSpinOp[curldgcd,divcd];
(* CurlDgTwist *)
CommuteFundSpinOp[curldgcd,twistcd]^=curldgcd[twistcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT},k/(k+1)twistcd[curldgcd[TT]][inds]-ToCanonicalSym@ImposeSym[(boxcd@@Take[{inds},-2])[TT@@Drop[{inds},-2]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]];
CommuteFundSpinOp[curldgcd,twistcd,twistcd,curldgcd]^=CommuteFundSpinOp[curldgcd,twistcd];
(* TwistCurlDg *)
CommuteFundSpinOp[twistcd,curldgcd]^=twistcd[curldgcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT},(k+1)/k*curldgcd[twistcd[TT]][inds]+(k+1)/k*ToCanonicalSym@ImposeSym[(boxcd@@Take[{inds},-2])[TT@@Drop[{inds},-2]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]]/;EnoughIndsQ[TT,1,0];
CommuteFundSpinOp[twistcd,curldgcd,curldgcd,twistcd]^=CommuteFundSpinOp[twistcd,curldgcd];
(* CurlCurlDgToCurlDgCurl *)
CommuteFundSpinOp[curlcd,curldgcd,curldgcd,curlcd]^=
curlcd[curldgcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT,B=DummyIn@spin,Bdg=DummyIn@spindg},
curldgcd[curlcd[TT]][inds]+
(1/(k+1)-1/(l+1))*twistcd[divcd[TT]][inds]-ToCanonicalSym@ImposeSym[(boxcd[First[Take[{inds},-k-l]],B])[TT[Sequence@@Drop[{inds},-k-l],Sequence@@Prepend[Rest[Take[{inds},-k-l]],-B]]],IndexList@@First@ExtractUnprimedAndPrimed[spin][{inds}]]+ToCanonicalSym@ImposeSym[(boxcd[Bdg,Last[{inds}]])[TT@@Append[Most[{inds}],-Bdg]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]
]/;EnoughIndsQ[TT,1,1];
(* CurlDgCurlToCurlCurlDg *)
CommuteFundSpinOp[curldgcd,curlcd,curlcd,curldgcd]^=curldgcd[curlcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT,B=DummyIn@spin,Bdg=DummyIn@spindg},
curlcd[curldgcd[TT]][inds]-
(1/(k+1)-1/(l+1))*twistcd[divcd[TT]][inds]+ToCanonicalSym@ImposeSym[(boxcd[First[Take[{inds},-k-l]],B])[TT[Sequence@@Drop[{inds},-k-l],Sequence@@Prepend[Rest[Take[{inds},-k-l]],-B]]],IndexList@@First@ExtractUnprimedAndPrimed[spin][{inds}]]-ToCanonicalSym@ImposeSym[(boxcd[Bdg,Last[{inds}]])[TT@@Append[Most[{inds}],-Bdg]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]
]/;EnoughIndsQ[TT,1,1];
(* DivTwistToCurlCurlDg *)
CommuteFundSpinOp[divcd,twistcd,curlcd,curldgcd]^=divcd[twistcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT,B=DummyIn@spin,Bdg=DummyIn@spindg},
-(1/(k+1)+1/(l+1))*curlcd[curldgcd[TT]][inds]+
l*(l+2)/(l+1)^2*twistcd[divcd[TT]][inds]-(l+2)/(l+1)*ToCanonicalSym@ImposeSym[(boxcd[First[Take[{inds},-k-l]],B])[TT[Sequence@@Drop[{inds},-k-l],Sequence@@Prepend[Rest[Take[{inds},-k-l]],-B]]],IndexList@@First@ExtractUnprimedAndPrimed[spin][{inds}]]-l/(l+1)*ToCanonicalSym@ImposeSym[(boxcd[Bdg,Last[{inds}]])[TT@@Append[Most[{inds}],-Bdg]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]
]/;EnoughIndsQ[TT,1,0];
(* DivTwistToCurlDgCurl *)
CommuteFundSpinOp[divcd,twistcd,curldgcd,curlcd]^=divcd[twistcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT,B=DummyIn@spin,Bdg=DummyIn@spindg},
-(1/(k+1)+1/(l+1))*curldgcd[curlcd[TT]][inds]+
k*(k+2)/(k+1)^2*twistcd[divcd[TT]][inds]-k/(k+1)*ToCanonicalSym@ImposeSym[(boxcd[First[Take[{inds},-k-l]],B])[TT[Sequence@@Drop[{inds},-k-l],Sequence@@Prepend[Rest[Take[{inds},-k-l]],-B]]],IndexList@@First@ExtractUnprimedAndPrimed[spin][{inds}]]-(k+2)/(k+1)*ToCanonicalSym@ImposeSym[(boxcd[Bdg,Last[{inds}]])[TT@@Append[Most[{inds}],-Bdg]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]
]/;EnoughIndsQ[TT,0,1];
(* CurlCurlDgToDivTwist *)
CommuteFundSpinOp[curlcd,curldgcd,divcd,twistcd]^=curlcd[curldgcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT,B=DummyIn@spin,Bdg=DummyIn@spindg},
-divcd[twistcd[TT]][inds]/(1/(k+1)+1/(l+1))+
((1+k)*l*(2+l))/((1+l)*(2+k+l))*twistcd[divcd[TT]][inds]-((1+k)*(2+l))/(2+k+l)*ToCanonicalSym@ImposeSym[(boxcd[First[Take[{inds},-k-l]],B])[TT[Sequence@@Drop[{inds},-k-l],Sequence@@Prepend[Rest[Take[{inds},-k-l]],-B]]],IndexList@@First@ExtractUnprimedAndPrimed[spin][{inds}]]-((1+k)*l)/(2+k+l)*ToCanonicalSym@ImposeSym[(boxcd[Bdg,Last[{inds}]])[TT@@Append[Most[{inds}],-Bdg]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]
]/;EnoughIndsQ[TT,1,0];
CommuteFundSpinOp[curlcd,curldgcd,twistcd,divcd]^=CommuteFundSpinOp[curlcd,curldgcd,divcd,twistcd];
(* CurlDgCurlToDivTwist *)
CommuteFundSpinOp[curldgcd,curlcd,divcd,twistcd]^=curldgcd[curlcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT,B=DummyIn@spin,Bdg=DummyIn@spindg},
-divcd[twistcd[TT]][inds]/(1/(k+1)+1/(l+1))+
(k*(2+k)*(1+l))/((1+k)*(2+k+l))*twistcd[divcd[TT]][inds]-(k*(1+l))/(2+k+l)*ToCanonicalSym@ImposeSym[(boxcd[First[Take[{inds},-k-l]],B])[TT[Sequence@@Drop[{inds},-k-l],Sequence@@Prepend[Rest[Take[{inds},-k-l]],-B]]],IndexList@@First@ExtractUnprimedAndPrimed[spin][{inds}]]-((2+k)*(1+l))/(2+k+l)*ToCanonicalSym@ImposeSym[(boxcd[Bdg,Last[{inds}]])[TT@@Append[Most[{inds}],-Bdg]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]
]/;EnoughIndsQ[TT,0,1];
CommuteFundSpinOp[curldgcd,curlcd,twistcd,divcd]^=CommuteFundSpinOp[curldgcd,curlcd,divcd,twistcd];
(* TwistDivToCurlCurlDg *)
CommuteFundSpinOp[twistcd,divcd,curlcd,curldgcd]^=twistcd[divcd[TT_?xTensorQ]][inds___]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT,B=DummyIn@spin,Bdg=DummyIn@spindg},
+(1+l)^2/(l*(2+l))*divcd[twistcd[TT]][inds]
+((1+l)*(2+k+l))/((1+k)*l*(2+l))*curlcd[curldgcd[TT]][inds]+(1+l^(-1))*ToCanonicalSym@ImposeSym[(boxcd[First[Take[{inds},-k-l]],B])[TT[Sequence@@Drop[{inds},-k-l],Sequence@@Prepend[Rest[Take[{inds},-k-l]],-B]]],IndexList@@First@ExtractUnprimedAndPrimed[spin][{inds}]]+(1+l)/(2+l)*ToCanonicalSym@ImposeSym[(boxcd[Bdg,Last[{inds}]])[TT@@Append[Most[{inds}],-Bdg]],IndexList@@Last@ExtractUnprimedAndPrimed[spin][{inds}]]
]/;EnoughIndsQ[TT,1,1];
(* Evaluate extensions *)
MakexTensions[DefFundSpinOperators,"End",covd,options];
]]]]];

UndefFundSpinOperators[covd_?CovDQ]:=
With[{
(* Names of the operators *)
divcd=DivName@covd, curlcd=CurlName@covd, curldgcd=CurlDgName@covd, twistcd=TwistName@covd,
(* Names of lists *)
$irrdecrulescd=StringJoin["$IrrDecRules",SymbolName[covd]],
$expandfundspinoprulescd=StringJoin["$ExpandFundSpinOpRules",SymbolName[covd]]},

(* Removing special code *)
ImposeFundSpinOp[expr1_,twistcd, iinds_]=.;
ImposeFundSpinOp[expr1_,curlcd, iinds_]=.;
ImposeFundSpinOp[expr1_,curldgcd, iinds_]=.;
ImposeFundSpinOp[expr1_,divcd, iinds_]=.;

(* Unlinking *)
VisitorsOf[covd]^=DeleteCases[VisitorsOf[covd],Alternatives[divcd,curlcd,curldgcd,twistcd]];

(* Undefining objects *)
Remove@$irrdecrulescd;
Remove@$expandfundspinoprulescd;

Remove@divcd;
Remove@curlcd;
Remove@curldgcd;
Remove@twistcd;

$FundSpinOpCovDs=DeleteCases[$FundSpinOpCovDs,covd];
];

8.3. Basic Tools

(* Converting the box to curvature *)
ExpandBox[expr_,covd_?SpacetimeSpinCovDQ]:=ToCanonical@ContractMetric[xAct`Spinors`Decomposition[xAct`Spinors`BoxToCurvature[expr,BoxName@covd],xAct`Spinors`Chi]];
ExpandBox[expr_,_]:=expr;
ExpandBox[expr_]:=If[Head@xAct`Spinors`Private`$SpinCovDs==List,Fold[ExpandBox,expr,xAct`Spinors`Private`$SpinCovDs],expr];

(* Expanding the fundamental spin operators *)
ExpandFundSpinOp[expr_,covd_?xAct`Spinors`SpinCovDQ]:=expr/.Symbol[StringJoin["$ExpandFundSpinOpRules",SymbolName[covd]]];
ExpandFundSpinOp[expr_]:=Fold[ExpandFundSpinOp,expr,$FundSpinOpCovDs];

(* Imposing the fundamental operators. These functions are general for all operators. Special code is added at the definition of the operators.*)
ImposeFundSpinOp[expr_Plus,op_?FundSpinOpQ, iinds_]:=ImposeFundSpinOp[#,op, iinds]&/@expr;
ImposeFundSpinOp[0,op_?FundSpinOpQ, iinds_]:=0;
ImposeFundSpinOp[c_?ConstantQ*expr_,op_?FundSpinOpQ, iinds_]:=c*ImposeFundSpinOp[expr,op, iinds];
ImposeFundSpinOp[TT_?xTensorQ[tinds___],op_?FundSpinOpQ, iinds_]:=(op[TT]@@Join[Select[{tinds},Not@MemberQ[$FundSpinVBundles,UpIndex@VBundleOfIndex@#]&],iinds])/;CompatibleSymQ[SlotsOfTensor[TT],SymmetryGroupOfTensor@TT,$FundSpinVBundles];

(* Applying a function inside a fundamental operator *)
FunctionInFundSpinOp[expr1_,func1_]:=FunctionInFundSpinOp[expr1,func1,1];
FunctionInFundSpinOp[expr1_,func1_,1,extrafunc_:Identity]:=expr1/.(op1_Symbol?FundSpinOpQ)[TT2_?xTensorQ][tinds1___]:>
Module[{unprimedinds,primedinds,otherinds},
{unprimedinds,primedinds,otherinds}=ExtractUnprimedAndPrimedIndsGivenCovD[{tinds1},CovDOfFundSpinOp@op1];
ImposeFundSpinOp[func1[TT2[Sequence@@otherinds,Sequence@@Drop[List@@GiveIndicesToTensor@TT2,Length@otherinds]]],op1,Join[unprimedinds,primedinds]]];
FunctionInFundSpinOp[expr1_,func1_,level_Integer,extrafunc_:Identity]:=FunctionInFundSpinOp[expr1,extrafunc@FunctionInFundSpinOp[#,func1,level-1]&,1];

8.4. Irreducible decompositions of covariant derivatives

Irreducible decompositions of covariant derivatives of symmetric spinors.

Options[IrrDecomposeCovD]^={ResultType->"Expression",ExpandSym->False};

(* Needs fixing for other VBundles -- Check that the order is right *)
IrrDecomposeCovD[covd_?SpacetimeSpinCovDQ[A_,Adg_]@TT_?xTensorQ[inds___],options:OptionsPattern[]]:=If[CompatibleSymQ[SlotsOfTensor@TT,SymmetryGroupOfTensor@TT]=!=True,Throw@Message[IrrDecomposeCovD::error,"Only decompositions of derivatives of symmetric tensors are implemented."],
With[{eps=First@MetricsOfVBundle@VBundlesOfCovD[covd][[2]],
epsdg=First@MetricsOfVBundle@VBundlesOfCovD[covd][[3]]},
Module[{primedinds,unprimedinds,irrdecexpr,k,l, resulttype,expandsym,otherinds},
{resulttype,expandsym}=OptionValue[{ResultType,ExpandSym}];
{unprimedinds,primedinds,otherinds}=ExtractUnprimedAndPrimedIndsGivenCovD[{inds},covd];
k=Length@unprimedinds;
l=Length@primedinds;
irrdecexpr=TwistName[covd][TT][Sequence@@otherinds,A,Sequence@@unprimedinds,Adg,Sequence@@primedinds];
If[l>0,irrdecexpr=irrdecexpr-l/(l+1)*ImposeSym[epsdg[Adg,First@primedinds]*CurlName[covd][TT][Sequence@@otherinds,A,Sequence@@unprimedinds,Sequence@@Rest@primedinds],IndexList@@primedinds]];
If[k>0,irrdecexpr=irrdecexpr-k/(k+1)*ImposeSym[eps[A,First@unprimedinds]*CurlDgName[covd][TT][Sequence@@otherinds,Sequence@@Rest@unprimedinds,Adg,Sequence@@primedinds],IndexList@@unprimedinds]];
If[And[k>0,l>0],irrdecexpr=irrdecexpr+l*k/(l+1)/(k+1)*ImposeSym[epsdg[Adg,First@primedinds]*eps[A,First@unprimedinds]*DivName[covd][TT][Sequence@@otherinds,Sequence@@Rest@unprimedinds,Sequence@@Rest@primedinds],IndexList@@Join[unprimedinds,primedinds]]];
If[expandsym,irrdecexpr=ToCanonical@ExpandSym[irrdecexpr,SmartExpand->True]];
Switch[resulttype,
"Expression",irrdecexpr,
"Equation",Equal[covd[A,Adg]@TT[inds],irrdecexpr],
"Rule",MakeRule[Evaluate[{covd[A,Adg]@TT[inds],irrdecexpr},MetricOn->All,ContractMetrics->True]],
"CompareRule",MakeCompareRule[Evaluate[{covd[A,Adg]@TT[inds],irrdecexpr}]],
_,Null]
]]];

We can also just give the name of the tensor and the covd.

IrrDecomposeCovD[TT_?xTensorQ,covd_?SpacetimeSpinCovDQ,options:OptionsPattern[]]:=Module[{numindsinvb,generatedinds,unprimedinds,primedinds,cdeexpr,slotsTT=SlotsOfTensor@TT,spin=VBundlesOfCovD[covd][[2]],spindg=VBundlesOfCovD[covd][[3]],unprimedpos,primedpos,otherinds,otherpos,otherrules,TTinds},numindsinvb=Tally[UpIndex/@SlotsOfTensor@TwistName[covd][TT]];
generatedinds=Rule[#1,GetIndexRange[#2,#1]]&@@@(numindsinvb);
unprimedinds=DownIndex/@(spin/.generatedinds);
primedinds=DownIndex/@(spindg/.generatedinds);
unprimedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spin)&];
primedpos=Select[Range@Length@slotsTT,(UpIndex[slotsTT[[#]]]==spindg)&];
otherinds=List@@@Select[generatedinds,And[First@#=!=spin,First@#=!=spindg]&];
otherpos=Flatten[Position[UpIndex/@slotsTT,#,1]]&/@First/@otherinds;
otherrules=Flatten[MapThread[Thread[Rule[#1,#2]]&,{otherpos,Last/@otherinds}]];
otherrules=Rule[#1,#2*sign@slotsTT[[#1]]]&@@@otherrules;
TTinds=Range@Length@slotsTT/.Thread[primedpos->Rest@primedinds]/.Thread[unprimedpos->Rest@unprimedinds]/.otherrules;
cdeexpr=covd[First@unprimedinds,First@primedinds][TT@@TTinds];
IrrDecomposeCovD[cdeexpr,options]];

IrrDecomposeCovD[TT_?xTensorQ,options:OptionsPattern[]]:=IrrDecomposeCovD[TT,First@Select[$CovDs,SpinCovDQ],options]

Collecting rules for irreducible decompositions.

Options[AddToIrrDecRules]^={ExpandSym->True};

AddToIrrDecRules[TT_?xTensorQ, covd_?xAct`Spinors`SpinCovDQ,options:OptionsPattern[]]:=
With[{$irrdecrulescd=StringJoin["$IrrDecRules",SymbolName[covd]]},
ToExpression[StringJoin["$IrrDecRules",SymbolName[covd]],InputForm,Function[name,name=DeleteDuplicates@Join[Evaluate@Symbol[$irrdecrulescd],IrrDecomposeCovD[TT,covd,ExpandSym->OptionValue[ExpandSym],ResultType->"Rule"]],HoldAll]]];

Options[ToFundSpinOp]^={ExpandSym->True};

ToFundSpinOp[expr_,covd_?xAct`Spinors`SpinCovDQ,options:OptionsPattern[]]:=Module[{newtensors={},intermexpr},intermexpr=ToCanonical@ContractMetric[expr//.Symbol[StringJoin["$IrrDecRules",SymbolName[covd]]]]/.x:(covd[A_,Adg_][TTT_?xTensorQ[inds___]]):>Module[{},AppendTo[newtensors,TTT];x];
newtensors=DeleteDuplicates@newtensors;
newtensors=Select[newtensors,CompatibleSymQ[SlotsOfTensor@#,SymmetryGroupOfTensor@#]&];
If[Length@newtensors>0,Print["Adding rules for ", newtensors]];
AddToIrrDecRules[#,covd,ExpandSym->OptionValue[ExpandSym]]&/@newtensors;
ToCanonical@ContractMetric[expr//.Symbol[StringJoin["$IrrDecRules",SymbolName[covd]]]]
];
ToFundSpinOp[expr_]:=Fold[ToFundSpinOp,expr,$FundSpinOpCovDs];

9. End private and package

Protect[InertHeadHead];

End[];
EndPackage[];

10. Examples

The user context should be empty :

?Global`*

<<xAct`TexAct`

Define a manifold

DefManifold[M4,4,{a,b,c,d,f,h,l,m,p,q}]

Define a two index tensor without symmetries.

$MixedSymVerbose=False;

DefTensor[T[-a,-b,-c,-d],M4]

DefTensor[S[-a,-b,-c,-d],M4]

ImposeSym[T[-a,-b,-c,-d]T[-f,-h,-l,-m],IndexList[-d,-b,-c],Symmetric[{1,2,3}]]

ToCanonical@%

DefTensor[U[-a,-b,-c,-d],M4, GenSet[-Cycles[{1,3},{2,4}]]]

S[-d,-f,-a,-q]T[a,q,f,d]+S[-f,-a,-d,-q]T[a,d,f,q]+U[-f,-a,-b,-q]T[a,q,b,f]+
S[-a,-c,-b,c]T[a,f,b,-f]

FindTensorCoefficients[%,{T[a,b,c,d]}]

Plus@@Times@@@%

ToCanonical[%-%%%]

UndefTensor[S]

ImposeSym[U[-a,-b,-c,-d],IndexList[-a,-b,-c,-d],JoinSGS[Symmetric[{1,2}],Symmetric[{3,4}]]]

SymmetryGroupOfTensor@Head@%

UndefTensor[U]

DefTensor[V[-a,-b,-c,-d,-f,-h,-l,-m,-p],M4,GenSet[Cycles[{1,4,7},{2,5,8},{3,6,9}]]]

SymmetryGroupOfTensor@V

ImposeSym[V[-a,-b,-c,-d,-f,-h,-l,-m,-p],IndexList[-a,-b,-c,-d,-f,-h,-l,-m,-p],JoinSGS[Symmetric[{1,2,3}],Symmetric[{4,5,6}],Symmetric[{7,8,9}]]]

SymmetryGroupOfTensor@Head@%

UndefTensor[V]

DefTensor[S[-a,-b,-c,-d],M4, Symmetric[{2,3,4}]]

ImposeSym[S[-a,-b,-c,-d],IndexList[-b,-c]]

ImposeSym[S[-a,-b,-c,-d],IndexList[-b,-c],ImposeSuperfluousSym->True]

ExpandSym[%,SmartExpand->True]

ImposeSym[S[-a,-b,-c,-d]T[-f,-h,-l,-m],IndexList[-d,-m],Symmetric[{1,2}],ImposeLargerSym->True]

SymmetryGroupOfTensor@Head@%

ImposeSym[S[-a,-b,-c,-d]T[-f,-h,-l,-m],IndexList[-d,-m],Symmetric[{2,4}],ImposeLargerSym->True]

CanonicalizeGroupInSym@%

ToCanonical@ExpandSym[%-%%]

UndefTensor[T]

UndefTensor[S]

DefTensor[T[-a,-b],M4]

Impose the symmetric group on the free indices:

ImposeSym[T[-b,-a]T[-c,b]]

The total symmetry group is just the imposed symmetry.

SymmetryGroupOfTensor@%

Find a canonical representation.

CanonicalizeGroupInSym@%%

ToCanonical@%

Print as TEX:

TexPrint@%

A larger example:

GS1=GenSet[Cycles[{1,3},{2,4}],Cycles[{4,5,6,7}],Cycles[{8,9}],Cycles[{9,10}]]

DefTensor[S[-a,-b,-c,-d,-f,-h,-l,-m,-p,-q],M4,GS1]

OrderOfGroup@SymmetryGroupOfTensor@S

expr1a=ImposeSym[S[-a,-b,-c,-d,-f,-h,d,c,-p,-q],IndexList[-c,-b,-f,-p,-q]]

expr2a=CanonicalizeGroupInSym@expr1a

SymmetryOf@expr2a

expr3a=ToCanonical@expr2a

Check that the original and the canonicalized expression is the same thing.

ToCanonical@ExpandSym[expr3a-expr1a]

A few examples where the imposed group is the antisymmetric group.

ImposeSym[T[-a,-b]T[-c,-d],IndexList[-b,-c],Antisymmetric[{1,2}]]

CanonicalizeGroupInSym@%

Expand@ExpandSym[%%-%]

ImposeSym[T[-a,-b]T[-c,-d],IndexList[-b,-c],Antisymmetric[{1,2}]]

SymmetryGroupOfTensor@Head@%

ToCanonical@%%

CanonicalizeGroupInSym@%

ToCanonical@%

Expand@ExpandSym[%%-%]

ImposeSym[T[-a,-b]T[-c,-d],IndexList[-b,-c],Antisymmetric[{1,2}]]

TexPrint@%

ToCanonical@CanonicalizeGroupInSym@%%

TexPrint@%

DefTensor[U[-a,-b,-c,-d],M4,GenSet[-Cycles[{1,2},{3,4}]]]

expr1=ImposeSym[U[-a,-b,-c,-d],IndexList[-a,-b,-c,-d],Antisymmetric[{1,2}]]

Checking that the mixed symmetry group is correct:

SymmetryGroupOfTensor@Head@expr1

List@@(OrderOfGroup[%]*ImposeSymmetry[expr1,IndexList[-a,-b,-c,-d],%])

ExpandSym/@%

ToCanonical@%

All four elements in the mixed symmetry group gives the same expanded expression.

Union@%

An example where RemoveSuperfluousInnerSym is needed.

DefTensor[V[-a,-b,-c],M4,GenSet[-Cycles[{2,3}]]]

ImposeSym[V[-a,-b,-c],IndexList[-a,-b],Symmetric[{1,2}]]

ImposeSym[%,IndexList[-a,-b,-c],Symmetric[{1,2,3}]]

CanonicalizeGroupInSym@%

The canonicalizers can not do anything about this.

ToCanonical@CanonicalizeGroupInSym@%

However, the inner symmetrization is not needed, because it is a part of the outer symmetrization. Therefore it can be removed and the combination with the inner antisymmetry makes the expression zero.

RemoveSuperfluousInnerSym@%

We can canonicalize internal expressions too, but only if the SymH is the only tensor in the headlist...

ImposeSym[S[-a,-b,-c,-d,-f,-h,d,c,-p,-q],IndexList[-c,-b,-f,-p,-q]]

ImposeSym[%,IndexList[-a,-f]]

CanonicalizeGroupInSym@%

ToCanonical@ExpandSym@ExpandSym[%-%%]

Internally the ZeroTensor is used to represent 0 as a tensor head.

ZeroTensor[{TangentM4}]

xTensorQ[%]

ZeroTensor[{TangentM4}][-a]

Using the ordinary 0 can cause problems

0[-a]

SymmetryGroupOfTensor for SymH objects returns the correct total symmetry group, but the general algorithm is not very efficeint. The special algorithms for for the cases of imposed symmetric or the antisymmetric group are more efficient though.

Canonicalization is now implemented for cases where the imposed symmetry is not the symmetric or the antisymmetric group.
But this code is not very efficient.

expr1=ImposeSym[T[-a,-b]T[-c,-d],IndexList[-a,-b,-c,-d],GenSet[-Cycles[{1,2},{3,4}]]]

CanonicalizeGroupInSym@expr1

SymmetryGroupOfTensor@Head@expr1

List@@(OrderOfGroup[%]*ImposeSymmetry[expr1,IndexList[-a,-b,-c,-d],%])

ExpandSym/@%

ToCanonical@%

All four elements in the mixed symmetry group gives the same expanded expression.

Union@%

DefTensor[DH[-a,-b,-c,-d],M4,StrongGenSet[{1,2,3,4},GenSet[Cycles[{1,2}],Cycles[{3,4}]]]]

ImposeSym[DH[-a,-b,-c,-d],IndexList[-a,-b,-c,-d],StrongGenSet[{1,2,3,4},GenSet[-Cycles[{1,3}],-Cycles[{2,4}]]]]

CanonicalizeGroupInSym@%

ExpandSym[%-%%]

ToCanonical@%

DefTensor[DH2[-a,-b,-c,-d],M4,Antisymmetric[{-b,-d}]]

SymH[{SymH[{DH2},StrongGenSet[{1,2,3,4},GenSet[Cycles[{1,2}],Cycles[{3,4}]]],"(12)(34)"]},StrongGenSet[{1,2,3,4},GenSet[-Cycles[{1,3}],-Cycles[{2,4}]]],"[13][24]"][-a,-b,-c,-d]

CanonicalizeGroupInSym@%

ToCanonical@ExpandSym@ExpandSym[%-%%]

UndefTensor[S]

UndefTensor[T]

UndefTensor[U]

UndefTensor[V]

Test of the irreducible decomposition code

<<xAct`Spinors`

Definition of a Lorentzian metric on the manifold:

SetOptions[DefTensor,DefInfo->False]

DefMetric[{1,3,0},g[-a,-b],CD,{";","▽"},DefInfo->False]

DefSpinStructure[g,Spin,{A,B,C,D,F,H,L,M,P,Q,R},ϵ,σ,CDe,{";","▽"},SpinorPrefix->SP,SpinorMark->"S"]

DefSpinor[T1[-A],M4]

DefSpinor[T2[-A,-B],M4]

DefSpinor[T4[-A,-B,-C,-D],M4,Symmetric[{-A,-B,-C}]]

DefSpinor[Ω4[-A,-B,-C,-D],M4,Symmetric[{-A,-B,-C,-D}]]

DefSpinor[Ω4b[-A,-B,-C,-D],M4,Symmetric[{-A,-B,-C,-D}]]

IrrDecompose[T1[-A]T1[A]]

IrrDecompose[T1[-A]T1[-B]]

IrrDecompose[T2[-A,-B]]

IrrDecompose[T2[-A,-B],ResultType->"Rule"]

IrrDecompose[T2[-A,-B],ResultType->"Equation"]

ToCanonical@ContractMetric[Times[ϵ[A,B],#]&/@%]

IrrDecompose[T2[A,B],ResultType->"Equation"]

ToCanonical@ContractMetric[Times[ϵ[-A,-B],#]&/@%]

IrrDecompose[T4[-A,-B,-C,-D],UseSym->False,ResultType->"Equation"]

IrrDecompose[T4[-A,-B,-C,-D],UseSym->True,ResultType->"Equation"]

Ω4IrrDec1=IrrDecompose[Ω4[-A,-B,-C,-D]Ω4[-F,-H,-L,-M],ResultType->"Equation",TimeVerbose->True]

Ω4IrrDec2=IrrDecompose[Ω4b[-A,-B,-C,-D]Ω4b[-F,-H,-L,-M],ResultType->"Equation",TimeVerbose->True,IrrDecImposeNewMethod->True]

The new method gives the same result slightly quicker.

(Ω4IrrDec2[[2]]/.Ω4b->Ω4)-Ω4IrrDec1[[2]]

DefSpinor[ϕ2[-A,-B],M4,Symmetric[{1,2}],PrintAs->"ϕ"]

SymH[{ϵ, ϵ, ϕ2}, StrongGenSet[{1, 5}, GenSet[Cycles[{1, 5}]]], 
  "(15)"][-C, A, -P, -Q, -L, -A]

ContractMetricsInsideSym@%

ExpandSym[%==%%]

SymH[{ϵ, ϵ, ϕ2}, StrongGenSet[{1, 5}, GenSet[Cycles[{1, 5}]]], 
  "(15)"][-C, -D, P, -P, -L, -F]

ContractMetricsInsideSym@%

ExpandSym[%==%%]

ToCanonical@ContractMetric@%

Test the irreducible decompositon of

ToCanonical@Antisymmetrize[ϵ[-A, -B]ϵ[-C, -D],{-A,-B,-C}]

IrrDecompose@%

Is is easy to see theoretically because we antisymmetrize over 3 indices on a 2 dimensional VBundle.

ϵ[-B, -C]*ϕ2[-D,-L]

IrrDecompose@%

ToCanonical@ContractMetric[1/2ϵ[B,C]*%]

ϵ[-B, -C]*ϕ2[-D,-L]

ToCanonical@ContractMetric[1/2ϵ[B,C]*%]

ϵ†[B†,-C†] ϕ2[-A,-C] ϕ2[-B,-D]ϕ2[-F,-H]

ImposeSym[%,IndexList[-C,-D]]

MoveTensorsOutsideSym[%,True]

MoveTensorsOutsideSym[%%]

ImposeSym[CDe[-A,A†]@PhiCDe[-B,-C,B†,C†],IndexList[-A,-B,-C,B†,C†]]

RemoveSuperfluousPartialSym@%

ToCanonical@ExpandSym[%-%%]

UndefSpinor/@{T1,T2,T4,ϕ2,Ω4};

The set stabilizer is not always the correct thing for the InternalCommutingSymmetry. The new code gives the correct answer though.

xAct`SymManipulator`Private`ClearRememberedSymHSymmetry[g];

$MixedSymSpecialCode=SetStabilizer;

SymmetryGroupOfTensor@SymH[{g, g}, Symmetric[{2, 3}], "(23)"]

InternalCommutingSymmetry@SymH[{g, g}, Symmetric[{2, 3}], "(23)"]

xAct`SymManipulator`Private`ClearRememberedSymHSymmetry[g];

$MixedSymSpecialCode=True;

SymmetryGroupOfTensor@SymH[{g, g}, Symmetric[{2, 3}], "(23)"]

InternalCommutingSymmetry@SymH[{g, g}, Symmetric[{2, 3}], "(23)"]

xAct`SymManipulator`Private`ClearRememberedSymHSymmetry[g];

$MixedSymSpecialCode=False;

SymmetryGroupOfTensor@SymH[{g, g}, Symmetric[{2, 3}], "(23)"]

InternalCommutingSymmetry@SymH[{g, g}, Symmetric[{2, 3}], "(23)"]

ImposeSym[g[a,b]g[c,d],IndexList[b,c],Symmetric[{1,2}]]

SymmetryGroupOfTensor@Head@%

Dimino@%

ImposeSym[g[a,b]g[c,d],IndexList[b,c],Antisymmetric[{1,2}]]

%-PermuteIndices[%,IndexList[a,b,c,d],Cycles[{1,2},{3,4}]]

Expand@ExpandSym@%

ToCanonical@%