New Binary Tree Operations (Prune, Contract)

src/Kernel/BTrees/Algebra.wl

@@ -0,0 +1,100 @@
+(* ::Package:: *)
+
+(* ::Section:: *)
+(*Private Members*)
+
+
+(*Inner sub tree function for only formal expansions*)
+innerBSubtrees[t_]:=Complement[DeleteDuplicates[subTrees[t]], {t}];
+
+ (*Expand out power terms in a list to individual terms*)
+listPowerRemove[x_]:=Flatten[Replace[x,{Power[y_,p_] :>ConstantArray[y,p]},{0,1}]];
+
+(*Flatten a Forest Space into a flat List *)
+linearCombToList[x_Plus]:=Flatten[Replace[List@@x,{Times[c_Integer,y_] :>ConstantArray[y,c]},{0,1}]];
+linearCombToList[Times[y_Integer,x_]]:=ConstantArray[x,y];
+linearCombToList[x_]:={x};
+
+(*Generalized Outer Product on two list with a Cofactor Expansion applied where there are no alternating signs*)
+outerACE[op_,t_List /;Length[t] ==2,k_List /;Length[k] ==2]:=op[t[[1]],k[[1]]]*op[t[[2]],k[[2]]] + op[t[[1]],k[[2]]]*op[t[[2]],k[[1]]];
+outerACE[op_,t_List,k_List] := Sum[op[First[t],k[[i]]]*outerACE[op,Drop[t,1],Delete[k,i]],{i,1,Length[k]}];
+
+(*Given a list of trees with powers in entries get the factorial of each power and product*)
+factorialProd[t_]:=Times@@(#!&/@Replace[Cases[t,_Power,{0,1}],{Power[y_,p_] :>p},{0,1}]);
+
+(*Given a Polynomial it will map a function to each varaible*)
+PolynomialMap[f_,poly_]:=FromCoefficientRules[CoefficientRules[poly],f/@Variables[poly]];
+
+
+getRoot[t_Symbol]:=t;
+getRoot[t_Subscript]:=t;
+getRoot[t_Power]:=0;
+getRoot[t_Times]:=0;
+getRoot[t_Plus]:=0;
+getRoot[t_]:=Head[t];
+
+
+getChildren[_Symbol | _Subscript | _Power | _Times | _Plus]:=0;
+getChildren[t_[h_]]:=h;
+
+
+subTrees[k_]:={k};
+subTrees[k_[t_]]:=Append[Map[k,subTrees[t]],k];
+subTrees[t_Times]:=Times@@@Most[Tuples[Map[Append[subTrees[#],1]&,List@@t]]];
+subTrees[t_^p_]:=Times@@@Most[Tuples[Append[subTrees[t],1],p]];
+
+
+prune[t_,t_]:=1;
+prune[t_,k_]:=0;
+prune[t_,1]:=t;
+prune[t_[u_],t_]:=u;
+prune[t_[u_],t_[g_]]:=prune[u,g];
+prune[t_Times, k_]:=Sum[prune[t[[i]],k]*(Length[t] - 1)!*Delete[t,i],{i,1,Length[t]}];
+prune[t_,k_Times]:=0;
+prune[t_Times, k_Times]:=With[{tpad = listPowerRemove[List@@t], kpad =listPowerRemove[List@@k] },If[Length[tpad]>=Length[kpad],(1/factorialProd[List@@k])*outerACE[prune,tpad,PadRight[kpad,Length[tpad],1]],0]];
+prune[t_^p_,k_Times]:=With[{kpad = listPowerRemove[List@@k]},If[p>=Length[kpad],(p! /factorialProd[List@@k] )*Product[prune[t,kpad[[i]]],{i,1,Length[kpad]}]*prune[t,1]^(p-Length[kpad]),0]];
+prune[t_Times,k_^q_]:=With[{tpad=listPowerRemove[List@@t]},If[Length[tpad]>=q,(1/q!)*outerACE[prune,tpad,PadRight[ConstantArray[k,q],Length[tpad],1]],0]];
+prune[t_^p_,k_]:=p*prune[t,k]*prune[t,1]^(p-1);
+prune[t_^p_,k_^q_]:=If[p>=q,Binomial[p,q]*prune[t,k]^(q)*prune[t,1]^(p-q),0];
+
+
+contract[t_,k_]:=If[getRoot[t]===getRoot[k],t,0];
+contract[t_,k_[u_]]:=If[getRoot[t]===k,Total[#*innerContract[ Flatten[linearCombToList[Expand[prune[t,#]]]],u]&/@innerBSubtrees[t],2],0];
+innerContract[tl_,u_]:=contract[#,u]&/@tl;
+contract[t_Times,k_]:=0;
+contract[t_,k_Times]:=0;
+contract[t_Times,k_[u_]]:=0;
+contract[t_[w_],k_Times]:=0;
+contract[t_^p_,k_]:=0;
+contract[t_,k_^q_]:=0;
+contract[t_^p_,k_[u_]]:=0;
+contract[t_[w_],k_^q_]:=0;
+contract[t_Times,k_Times]:=With[{tpad = listPowerRemove[List@@t],kpad=listPowerRemove[List@@k]},If[Length[tpad]==Length[kpad],(1/factorialProd[List@@k])*outerACE[contract,tpad,kpad],0]];
+contract[t_Times,k_^q_]:=With[{tpad=listPowerRemove[List@@t]},If[Length[tpad]==q,Product[contract[tpad[[i]],k],{i,1,Length[tpad]}],0]];
+contract[t_^p_,k_Times]:=With[{kpad = listPowerRemove[List @@k]},If[p==Length[kpad],(p!/factorialProd[List@@k])*Product[contract[t,kpad[[i]]],{i,1,Length[kpad]}],0]];
+contract[t_^p_,k_^q_]:=If[p==q,contract[t,k]^p,0];
+
+
+(* ::Section:: *)
+(*Package Definitions*)
+
+
+BTreePrune[h_[t_,p___],h_[t_,p___]]:=1;
+BTreePrune[h_[t_,p___],h_[\[FormalY],p___]]:=h[t,p];
+BTreePrune[h_[t_,p___],h_[k_,p___]]:=With[{poly = Expand[prune[t,k]]},If[NumericQ[poly],poly,PolynomialMap[h[#,p]&,poly]]];
+
+
+BTreeContract[h_[\[FormalY],p___],h_[\[FormalY],p___]]:=h[\[FormalY],p];
+BTreeContract[h_[t_,p___],h_[k_,p___]]:=With[{poly = Expand[contract[t,k]]},If[NumericQ[poly],poly,PolynomialMap[h[#,p]&,poly]]];
+
+
+BTreeSubTrees[h_[t_,p___]]:=Map[h[#,p]&,Complement[DeleteDuplicates[subTrees[t]], {t}]];
+
+
+BTreeRoot[h_[t_,p___]]:=h[getRoot[t],p];
+SetAttributes[BTreeRoot, Listable];
+
+
+BTreeChildren[h_[_Symbol | _Subscript,p___]]:=h[\[FormalY],p];
+BTreeChildren[h_[t_[k_],p___]]:=PolynomialMap[h[#,p]&,Expand[getChildren[t]]];
+SetAttributes[BTreeChildren, Listable];

src/Kernel/BTrees/BTrees.wl

@@ -0,0 +1,41 @@
+(* ::Package:: *)
+
+BeginPackage["Integreat`BTrees`"];
+
+
+BTree::usage =
+	"BTree[p] generates all B-Trees up to order p grouped by order.\n" <>
+	"BTree[{p}] generates a list of B-Trees of order p.";
+BTreeN::usage =
+	"BTreeN[p] generates all N-Trees up to order p grouped by order.\n" <>
+	"BTreeN[{p}] generates a list of N-Trees of order p.";
+BTreeDAE::usage =
+	"BTreeDAE[p] generates all DAE-Trees up to order p grouped by order.\n" <>
+	"BTreeDAE[{p}] generates a list of DAE-Trees of order p.";
+BTreeFormalForm::usage = "Get the Formal Representation of the Tree";
+BTreeQ::usage = "Check if the argument is a BTree"
+
+
+BTreeOrder::usage = "BTreeOrder[t] computes the order of the tree t.";
+BTreeAlpha::usage = "BTreeAlpha[t] computes the number of monotonic labelings of the tree t.";
+BTreeGamma::usage = "BTreeGamma[t] computes the density of the tree t.";
+BTreeSigma::usage = "BTreeSigma[t] computes the number of symmetries of the tree t.";
+
+
+BTreePrune::usage = "Get all prunings of tree T with tree S and get its Forest Space";
+BTreeContract::usage = "Get all cuts of tree T which contract to tree S and get its Forest Space";
+BTreeSubTrees::usage = "Get All Subtrees of a Tree as a List";
+BTreeRoot::usage = "Get the root node of a tree";
+BTreeChildren::usage = "Get children of a tree and return as a Forest Space";
+
+
+Begin["`Private`"];
+
+<<Integreat`BTrees`Core`;
+<<Integreat`BTrees`IntegerProperties`;
+<<Integreat`BTrees`Algebra`;
+
+End[];
+
+
+EndPackage[];

src/Kernel/BTrees.wl → src/Kernel/BTrees/Core.wl

@@ -1,40 +1,18 @@
 (* ::Package:: *)
 
-(* ::Section:: *)
-(*Usage*)
-
-
-BeginPackage["Integreat`BTrees`"];
-
-BTree::usage =
-	"BTree[p] generates all B-Trees up to order p grouped by order.\n" <>
-	"BTree[{p}] generates a list of B-Trees of order p.";
-BTreeN::usage =
-	"BTreeN[p] generates all N-Trees up to order p grouped by order.\n" <>
-	"BTreeN[{p}] generates a list of N-Trees of order p.";
-BTreeDAE::usage =
-	"BTreeDAE[p] generates all DAE-Trees up to order p grouped by order.\n" <>
-	"BTreeDAE[{p}] generates a list of DAE-Trees of order p.";
-
-BTreeOrder::usage = "BTreeOrder[t] computes the order of the tree t.";
-BTreeAlpha::usage = "BTreeAlpha[t] computes the number of monotonic labelings of the tree t.";
-BTreeGamma::usage = "BTreeGamma[t] computes the density of the tree t.";
-BTreeSigma::usage = "BTreeSigma[t] computes the number of symmetries of the tree t.";
-
-
 (* ::Section:: *)
 (*Private Members*)
 
 
-Begin["`Private`"];
-
 treeMap[head_, treeFun_, {p_}] := Map[head, treeFun[p]];
 treeMap[head_, treeFun_, p_] := Map[head, Table[treeFun[i], {i, p}], {2}];
 
+
 (* Allow KroneckerProduct to accept 1 argument. For lists, this works like TensorProduct but has less List nesting and is faster. *)
 kron[x_] := x;
 kron[x__] := KroneckerProduct[x];
 
+
 (* Generates all sets of subtrees with at least minSubtrees elements and p total vertices *)
 subtrees[head_, p_, minSubtrees_:1] := DeleteDuplicates[Flatten[kron @@@ Map[head, IntegerPartitions[p, {minSubtrees, p}], {2}]]];
 
@@ -63,50 +41,22 @@ treeN[n_][p_] := treeN[n][p] = Flatten[Outer[Construct, treeN[n][1], subtrees[tr
 
 treeDiffAlg[idx_][p_] := Join[treeDiff[idx][p], treeAlg[idx][p]];
 
+
 treeDiff[idx_, start_:1, q_:0][p_] := Catenate[Table[treeDiff[idx, {i}, q][p], {i, start, Max[1, idx - 1]}]];
 treeDiff[_, {part_}, _][0] := {Subscript[\[FormalY], part]};
 treeDiff[_, {part_}, _][1] := {Subscript[\[FormalF], part]};
 treeDiff[idx_, {1}, _][p_] := treeDiff[idx, {1}, _][p] = Map[Subscript[\[FormalF], 1], subtrees[treeDiffAlg[idx], p - 1]];
 treeDiff[idx_, {1}, p_][p_] := treeDiff[idx, {1}, p][p] = Map[Subscript[\[FormalF], 1], Join[subtrees[treeDiffAlg[idx], p - 1, 2], treeDiff[idx][p - 1]]];
 treeDiff[idx_, {part_}, q_][p_] := treeDiff[idx, {part}, q][p] = Map[Subscript[\[FormalF], part], subtrees[treeDiff[idx, part - 1, q], p - 1]];
 
+
 treeAlg[idx_][0] := {Subscript[\[FormalY], Max[2, idx]]};
 treeAlg[0][1] := {Subscript[\[FormalG], 0]};
 treeAlg[0][p_] := treeAlg[0][p] = Map[Subscript[\[FormalG], 0], subtrees[treeDiffAlg[0], p - 1]];
 treeAlg[1][p_] := treeAlg[1][p] = Map[Subscript[\[FormalG], 1], Join[subtrees[treeDiffAlg[1], p, 2], treeDiff[1][p]]];
 treeAlg[idx_][p_] := treeAlg[idx][p] = Map[Subscript[\[FormalG], idx], subtrees[treeDiff[idx, {idx - 1}, p + 1], p + idx - 1]];
 
 
-(* ::Subsection:: *)
-(*Tree Functions*)
-
-
-treeOrder[\[FormalY] | Subscript[\[FormalY], _]] := 0;
-treeOrder[_Symbol | _Subscript] := 1;
-treeOrder[Power[t_, p_]] := p * treeOrder[t];
-treeOrder[t_Times] := Map[treeOrder, Plus @@ t];
-treeOrder[r:(\[FormalF] | Subscript[\[FormalF], _])[t_]] := treeOrder[r] = treeOrder[t] + 1;
-treeOrder[r:Subscript[\[FormalG], idx_][t_]] := treeOrder[r] = treeOrder[t] + 1 - idx;
-
-treeAlpha[_Symbol | _Subscript] := 1;
-treeAlpha[Power[t_, p_]] := With[{o = treeOrder[t]},
-	Quotient[Pochhammer[p + 1,  p * (o - 1)], (o!)^p] * treeAlpha[t]^p
-];
-treeAlpha[t_Times] := Apply[Multinomial, Map[treeOrder, List @@ t]] * Map[treeAlpha, t];
-treeAlpha[r:_[t_]] := treeAlpha[r] = treeAlpha[t];
-
-treeGamma[_Symbol | _Subscript] := 1;
-treeGamma[Power[t_, p_]] := treeGamma[t]^p;
-treeGamma[t_Times] := Map[treeGamma, t];
-treeGamma[r:(\[FormalF] | _Subscript)[t_]] := treeGamma[r] = treeOrder[r] * treeGamma[t];
-treeGamma[r:Subscript[\[FormalG], idx_][t_]] := treeGamma[r] = treeGamma[t] / Pochhammer[treeOrder[r] + 1, idx - 1];
-
-treeSigma[_Symbol | _Subscript] := 1;
-treeSigma[Power[t_, p_]] := Factorial[p] * treeSigma[t]^p;
-treeSigma[t_Times] := Map[treeSigma, t];
-treeSigma[r:_[t_]] := treeSigma[r] = treeSigma[t];
-
-
 (* ::Subsection:: *)
 (*Conversions*)
 
@@ -118,6 +68,7 @@ toTree[Power[t_, p_]] := Splice[ConstantArray[toTree[t], p]];
 toTree[t_Times] := Splice[Map[toTree, List @@ t]];
 toTree[r:h_[t_]] := toTree[r] = Tree[h, {toTree[t]}];
 
+
 toGraph[t_] := With[{
 		tree = toTree[t]
 	},
@@ -143,31 +94,28 @@ toGraph[t_] := With[{
 	]
 ];
 
+
 toBoxes[\[FormalY], format_] := ToBoxes["\[EmptySet]", format];
 toBoxes[Subscript[\[FormalY], i_], format_] := ToBoxes[Subscript["\[EmptySet]", i], format];
-If [
-	TrueQ[$VersionNumber >= 12.3],
-	toBoxes[t_, format_] := ToBoxes[toGraph[t], format];,
-	toBoxes[t_, format_] := ToBoxes[t, format];
-];
+toBoxes[t_, format_] := ToBoxes[toGraph[t], format];
 
 
 (* ::Section:: *)
 (*Package Definitions*)
 
 
 BTree[p:_Integer?Positive | {_Integer?NonNegative}] := treeMap[BTree, tree, p];
-BTree /: Tree[HoldPattern[BTree[t_]]] := toTree[t];
-BTree /: MakeBoxes[HoldPattern[BTree[t_]], format_] := toBoxes[t, format];
+BTree /: MakeBoxes[HoldPattern[r:BTree[t_]], format_] := With[{b = toBoxes[t, format]}, InterpretationBox[b, r]];
+
 
 Options[BTreeN] = {"Partitions" -> 2};
 BTreeN[p:_Integer?Positive | {_Integer?NonNegative}, OptionsPattern[]] := With[{
 		n = OptionValue[Partitions]
 	},
 	treeMap[BTreeN[#, n] &, treeN[n], p] /; IntegerQ[n] && Positive[n]
 ];
-BTreeN /: Tree[HoldPattern[BTreeN[t_, _Integer]]] := toTree[t];
-BTreeN /: MakeBoxes[HoldPattern[BTreeN[t_, _Integer]], format_] := toBoxes[t, format];
+BTreeN /: MakeBoxes[HoldPattern[r:BTreeN[t_, _Integer]], format_] := With[{b = toBoxes[t, format]}, InterpretationBox[b, r]];
+
 
 Options[BTreeDAE] = {"Index" -> 1, "Partition" -> All};
 BTreeDAE[p:_Integer?Positive | {_Integer?NonNegative}, OptionsPattern[]] := With[{
@@ -181,25 +129,14 @@ BTreeDAE[p:_Integer?Positive | {_Integer?NonNegative}, OptionsPattern[]] := With
 	},
 	treeMap[BTreeDAE[#, idx] &, treeFun[idx], p] /; IntegerQ[idx] && NonNegative[idx] && !TrueQ[treeFun] 
 ];
-BTreeDAE /: Tree[HoldPattern[BTreeDAE[t_, _Integer]]] := toTree[t];
-BTreeDAE /: MakeBoxes[HoldPattern[BTreeDAE[t_, _Integer]], format_] := toBoxes[t, format];
-
-BTreeOrder[(BTree | BTreeN | BTreeDAE)[t_, ___]] := treeOrder[t];
-SetAttributes[BTreeOrder, Listable];
-
-BTreeAlpha[(BTree | BTreeN | BTreeDAE)[t_, ___]] := treeAlpha[t];
-SetAttributes[BTreeAlpha, Listable];
+BTreeDAE /: MakeBoxes[HoldPattern[r:BTreeDAE[t_, _Integer]], format_] := With[{b = toBoxes[t, format]}, InterpretationBox[b, r]];
 
-BTreeGamma[(BTree | BTreeN | BTreeDAE)[t_, ___]] := treeGamma[t];
-SetAttributes[BTreeGamma, Listable];
 
-BTreeSigma[(BTree | BTreeN | BTreeDAE)[t_, ___]] := treeSigma[t];
-SetAttributes[BTreeSigma, Listable];
-
-
-(* ::Section:: *)
-(*End Package*)
+BTreeQ[BTree[_] | BTreeN[_, _] | BTreeDAE[_, _]] := True;
+BTreeQ[BTree | BTreeN | BTreeDAE] := True;
+BTreeQ[_] := False;
 
 
-End[];
-EndPackage[];
+BTreeFormalForm[h_[t_,___]]/;BTreeQ[h]:=t;
+BTreeFormalForm[h_]/;NumericQ[h]:=h;
+SetAttributes[BTreeFormalForm, Listable];

src/Kernel/BTrees/IntegerProperties.wl

@@ -0,0 +1,53 @@
+(* ::Package:: *)
+
+(* ::Section:: *)
+(*Private Members*)
+
+
+treeOrder[\[FormalY] | Subscript[\[FormalY], _]] := 0;
+treeOrder[_Symbol | _Subscript] := 1;
+treeOrder[Power[t_, p_]] := p * treeOrder[t];
+treeOrder[t_Times] := Map[treeOrder, Plus @@ t];
+treeOrder[r:(\[FormalF] | Subscript[\[FormalF], _])[t_]] := treeOrder[r] = treeOrder[t] + 1;
+treeOrder[r:Subscript[\[FormalG], idx_][t_]] := treeOrder[r] = treeOrder[t] + 1 - idx;
+
+
+treeAlpha[_Symbol | _Subscript] := 1;
+treeAlpha[Power[t_, p_]] := With[{o = treeOrder[t]},
+	Quotient[Pochhammer[p + 1,  p * (o - 1)], (o!)^p] * treeAlpha[t]^p
+];
+treeAlpha[t_Times] := Apply[Multinomial, Map[treeOrder, List @@ t]] * Map[treeAlpha, t];
+treeAlpha[r:_[t_]] := treeAlpha[r] = treeAlpha[t];
+
+
+treeGamma[_Symbol | _Subscript] := 1;
+treeGamma[Power[t_, p_]] := treeGamma[t]^p;
+treeGamma[t_Times] := Map[treeGamma, t];
+treeGamma[r:(\[FormalF] | _Subscript)[t_]] := treeGamma[r] = treeOrder[r] * treeGamma[t];
+treeGamma[r:Subscript[\[FormalG], idx_][t_]] := treeGamma[r] = treeGamma[t] / Pochhammer[treeOrder[r] + 1, idx - 1];
+
+
+treeSigma[_Symbol | _Subscript] := 1;
+treeSigma[Power[t_, p_]] := Factorial[p] * treeSigma[t]^p;
+treeSigma[t_Times] := Map[treeSigma, t];
+treeSigma[r:_[t_]] := treeSigma[r] = treeSigma[t];
+
+
+(* ::Section:: *)
+(*Package Definitions*)
+
+
+BTreeOrder[t_?BTreeQ] := treeOrder[First[t]];
+SetAttributes[BTreeOrder, Listable];
+
+
+BTreeAlpha[t_?BTreeQ] := treeAlpha[First[t]];
+SetAttributes[BTreeAlpha, Listable];
+
+
+BTreeGamma[t_?BTreeQ] := treeGamma[First[t]];
+SetAttributes[BTreeGamma, Listable];
+
+
+BTreeSigma[t_?BTreeQ] := treeSigma[First[t]];
+SetAttributes[BTreeSigma, Listable];