(* ::Package:: *)

(* ::Input:: *)
(*Usage*)


(* ::Input:: *)
(*(*SackinMinimizer considers only rooted binary trees, ordered rooted binary trees and rooted binary phylogenetic trees -- all with minimal Sackin index.*)
(**)
(*Evaluate the entire noteboook. Keep it open in the background. Open a new (blank) notebook to work in.*)
(**)
(*Then decide what you want to do.*)
(**)
(*ROOTED BINARY TREES: *)
(**)
(*You can get the value of the minimal Sackin index for n leaves: type getminSackinvalue[n] in an empty cell of your file and evaluate the cell.*)
(**)
(*You can get the number of rooted binary trees with n leaves that have minimal Sackin index:*)
(*type numberoftreeswithMinSackinIndex[n] in an empty cell of your file and evaluate the cell.*)
(**)
(*You can get a complete list of Sackin minimal trees with n leaves in a Newick-like format, i.e. as nested list {{x,x},x}, where x is the placeholder for leaves:*)
(*type getSackinMinima[n] in an empty cell of your file and evaluate the cell.*)
(**)
(*You can get a complete list of Sackin minimal trees with n leaves plotted in the TreeGraph format:*)
(*type plotSackinMinima[n] in an empty cell of your file and evaluate the cell.*)
(**)
(**)
(*ORDERED ROOTED BINARY TREES:*)
(**)
(*You can get the number of ordered rooted binary trees with n leaves that have minimal Sackin index:*)
(*type numberoforderedtreeswithMinSackinIndex[n] in an empty cell of your file and evaluate the cell.*)
(**)
(*You can get a complete list of Sackin minimal ordered trees with n leaves in a Newick-like format, i.e. as nested list {{x,x},x}, where x is the placeholder for leaves:*)
(*type getOrderedSackinMinima[n] in an empty cell of your file and evaluate the cell.*)
(**)
(*You can get a complete list of Sackin minimal ordered trees with n leaves plotted in the TreeGraph format:*)
(*type plotOrderedSackinMinima[n] in an empty cell of your file and evaluate the cell.*)
(**)
(**)
(*ROOTED BINARY PHYLOGENETIC TREES:*)
(**)
(*You can get the number of ordered rooted binary trees with n leaves that have minimal Sackin index:*)
(*type numberofphylotreeswithMinSackinIndex[n] in an empty cell of your file and evaluate the cell.*)
(**)
(**)*)


(* ::Input:: *)
(*Main functions *)


(*Generate all ORDERED rooted binary trees with n leaves and with minimal Sackin index. The output is given in the anonymous Newick format (cf. Usage).*)
getOrderedSackinMinima[n_]:=Module[{treeset,rule},treeset=generatePreliminaryOrderedSackinMinima[n];
rule=Join[Table[ToString[i]->"x",{i,1,n}],Table[StringJoin[ToString[i],"a"]->"x",{i,1,n}],Table[StringJoin[ToString[i],"b"]->"x",{i,1,n}]];
Return[treeset/.rule];


];




(*Generate all (non-ordered) rooted binary trees with n leaves and with minimal Sackin index. The output is given in the anonymous Newick format (cf. Usage).*)
getSackinMinima[n_]:=DeleteDuplicates[completeSort[getOrderedSackinMinima[n]]];



(*Plot all ORDERED rooted binary trees with n leaves and with minimal Sackin index. The output is depicted as a graph and is a TreeGraph object.*)
plotOrderedSackinMinima[n_]:=Module[{treelist,i},treelist=generatePreliminaryOrderedSackinMinima[n];For[i=1,i<=Length[treelist],i++,
Print[TreeGraph[getEdges[treelist[[i]]],GraphLayout->{"LayeredEmbedding","RootVertex"->treelist[[i]]}]]


]];




(*Plot all (non-ordered) rooted binary trees with n leaves and with minimal Sackin index. The output is depicted as a graph and is a TreeGraph object.*)
plotSackinMinima[n_]:=Module[{treelist,i},treelist=generatePreliminarySackinMinima[n];For[i=1,i<=Length[treelist],i++,
Print[TreeGraph[getEdges[treelist[[i]]],GraphLayout->{"LayeredEmbedding","RootVertex"->treelist[[i]]}]]


]];


(*Calculate minimal Sackin index for trees with n leaves*)
getminSackinvalue[n_]:=Module[{k=Ceiling[Log2[n]]},Return[(-2^k)+n*(k+1)]];

(*Calculate number of rooted binary trees with minimal Sackin index*)
numberoftreeswithMinSackinIndex[n_]:=Module[{parts=IntegerPartitions[n,{2}],s,na,nb,total=0,ka,kb,k},
If[n==1,Return[1]];
If[n==2,Return[1]];
If[n==3,Return[1]];

For[s=1,s<=Length[parts],s++,
na=parts[[s]][[1]];nb=parts[[s]][[2]];k=Ceiling[Log2[n]];
ka=Ceiling[Log2[na]];kb=Ceiling[Log2[nb]];

If[na>=n/2&&ka==k-1&&kb==k-1&&na!=nb,
total=total+numberoftreeswithMinSackinIndex[na]*numberoftreeswithMinSackinIndex[nb],

If[na>=n/2&&ka==k-1&&kb==k-1&&na==nb,total=total+Binomial[numberoftreeswithMinSackinIndex[na]-1+2,2],If[na>=n/2&&nb==2^(k-2)&&ka==k-1,total=total+numberoftreeswithMinSackinIndex[na]]]];


];


Return[total];


]




(*Calculate number of ordered trees with minimal Sackin index*)
numberoforderedtreeswithMinSackinIndex[n_]:=Module[{parts,s,na,nb,total=0,ka,kb,k,length},
If[n==1||n==2,Return[1],parts=IntegerPartitions[n,{2}];length=Length[parts];
For[s=1,s<=length,s++,

na=parts[[s]][[1]];nb=parts[[s]][[2]];k=Ceiling[Log2[n]];
ka=Ceiling[Log2[na]];kb=Ceiling[Log2[nb]];

If[n==2^k,Return[1],
If[na>=n/2&&ka==k-1&&kb==k-1&&na!=nb,
total=total+2*numberoforderedtreeswithMinSackinIndex[na]*numberoforderedtreeswithMinSackinIndex[nb],
If[na>=n/2&&ka==k-1&&nb==2^(k-2),total=total+2*numberoforderedtreeswithMinSackinIndex[na],

If[na>=n/2&&ka==k-1&&kb==k-1&&na==nb&&n!=2^k,

total=total+numberoforderedtreeswithMinSackinIndex[na]*numberoforderedtreeswithMinSackinIndex[nb];];

]
]
]]
];
Return[total];

]


(*Calculate number of phylogenetic trees with minimal Sackin index*)
numberofphylotreeswithMinSackinIndex[n_]:=Module[{parts,s,na,nb,total=0,ka,kb,k,length},
If[n==1||n==2,Return[1],parts=IntegerPartitions[n,{2}];length=Length[parts];
For[s=1,s<=length,s++,

na=parts[[s]][[1]];nb=parts[[s]][[2]];k=Ceiling[Log2[n]];
ka=Ceiling[Log2[na]];kb=Ceiling[Log2[nb]];

If[n==2^k,Return[n!/(2^(n-1))],
If[

na>=n/2&&ka==k-1&&kb==k-1&&na!=nb,
total=total+Binomial[n,na]*numberofphylotreeswithMinSackinIndex[na]*numberofphylotreeswithMinSackinIndex[nb],
If[na>=n/2&&ka==k-1&&nb==2^(k-2),total=total+(nb!/(2^(nb-1)))*numberofphylotreeswithMinSackinIndex[na]*Binomial[n,nb],

If[na>=n/2&&ka==k-1&&kb==k-1&&na==nb&&n!=2^k,

total=total+1/2*Binomial[n,na]*(numberofphylotreeswithMinSackinIndex[na]*numberofphylotreeswithMinSackinIndex[na]);];

]
]
]]
];
Return[total];

]


(* ::Input:: *)
(*Supporting internal functions *)


(*Sort all levels of a nested list*)
completeSort[expr_]:=Map[Sort,expr,{0,-2}];
(*Subdivide a list of an even number n of elements into two equally long sublists, such that the first list contains elements 1 to n/2 and the second list contains elements n/2+1 to n.*)bipartiteSetpower2[list_]:=Module[{n=Length[list]},If[n>2&&Mod[n,2]==0,Return[{bipartiteSetpower2[list[[1;;n/2]]],bipartiteSetpower2[list[[n/2+1;;n]]]}],Return[list]]];



(*For a given number n, return the smallest possible fully balanced binary tree that contains at most n leaves.*)
getStartTree[n_]:=Module[{tree,k=Ceiling[Log2[n]],m},m=2^(k-1);
If[n==2^k,Return[bipartiteSetpower2[Range[n]]],Return[bipartiteSetpower2[Range[m]]]];

];



(*Internal function that is required to generate all binary ordered trees with minimal Sackin index*)generatePreliminaryOrderedSackinMinima[n_]:=Module[{tree,m,d,i,list,set,treeset,rule,k,table,subsets,numberofminima,j},tree=getStartTree[n];list=Flatten[tree];m=Length[list];d=n-m;If[d==0,rule=Table[i->ToString[i],{i,1,m}];
Return[{tree}/.rule]
Return[{tree}],
subsets=Subsets[list,{d}];
treeset={};
For[i=1,i<=Length[subsets],i++,
set=subsets[[i]];table=Table[{"x","x"},{k,1,d}];
(*rule=Table[set[[j]]\[Rule]table[[j]],{j,1,d}];*)
rule=Table[set[[j]]->{StringJoin[ToString[set[[j]]],"a"],StringJoin[ToString[set[[j]]],"b"]},{j,1,d}];
AppendTo[treeset,tree/.rule];


];

rule=Table[i->ToString[i],{i,1,m}];
Return[treeset/.rule];

]
]



(*Internal function that is required to generate all (non-ordered)binary trees with minimal Sackin index*)generatePreliminarySackinMinima[n_]:=Module[{treeset,length,i,j,rule,delete},treeset=generatePreliminaryOrderedSackinMinima[n];
length=Length[treeset];rule=Join[Table[ToString[i]->"x",{i,1,n}],Table[StringJoin[ToString[i],"a"]->"x",{i,1,n}],Table[StringJoin[ToString[i],"b"]->"x",{i,1,n}]];delete={};
For[i=1,i<=length-1,i++,
For[j=i+1,j<=length,j++,
If[completeSort[treeset[[i]]/.rule]==completeSort[treeset[[j]]/.rule],AppendTo[delete,{j}]];
]


]; treeset=Delete[treeset,delete];
Return[treeset];
];

(*Internal function to get all edges from a nested (Newick) tree representation*)
getEdges[rootedbinarytree_]:=Module[{i,j,list,length,edges={}},list=Append[Sort[Level[rootedbinarytree,Depth[rootedbinarytree]]],rootedbinarytree];
length=Length[list];
For[i=1,i<=length-1,i++,
For[j=i+1,j<=length,j++,

If[MemberQ[list[[i]],list[[j]]],AppendTo[edges,list[[i]]\[DirectedEdge]list[[j]]],
If[MemberQ[list[[j]],list[[i]]],AppendTo[edges,list[[j]]\[DirectedEdge]list[[i]]]]];
];
];Return[edges];
]

