-- -*- utf-8 -*-
newPackage(
"PhylogeneticTrees",
Version => "2.0",
Date => "November 15, 2019",
Headline => "invariants for group-based phylogenetic models",
--HomePage => "",
Authors => {
{Name => "Hector Baños", Email => "hbanos@gatech.edu"},
{Name => "Nathaniel Bushek", Email => "nbushek@css.edu"},
{Name => "Ruth Davidson", Email => "ruth.davidson.math@gmail.com"},
{Name => "Elizabeth Gross", Email => "egross@hawaii.edu"},
{Name => "Pamela Harris", Email => "peh2@williams.edu"},
{Name => "Robert Krone", Email => "rckrone@gmail.com"},
{Name => "Colby Long", Email => "clong@wooster.edu"},
{Name => "AJ Stewart", Email => "stewaral@seattleu.edu"},
{Name => "Robert Walker", Email => "robmarsw@umich.edu"}
},
PackageImports => {
"FourTiTwo"
},
PackageExports => {
"Graphs",
"Posets"
},
DebuggingMode => false
)
export {
"qRing",
"pRing",
"secant",
"joinIdeal",
"phyloToricFP",
"phyloToric42",
"phyloToricLinears",
"phyloToricQuads",
"phyloToricRandom",
"phyloToricAMatrix",
"toricSecantDim",
"toricJoinDim",
"Model",
"CFNmodel", "JCmodel", "K2Pmodel", "K3Pmodel",
"leafColorings",
"model",
"buckets",
"group",
"LeafTree",
"leafTree",
"internalEdges",
"internalVertices",
"edgeCut",
"vertexCut",
"edgeContract",
"QRing",
"fourierToProbability",
"labeledTrees",
"labeledBinaryTrees",
"rootedTrees",
"rootedBinaryTrees",
"unlabeledTrees",
"isIsomorphic"
}
protect \ {Group, Automorphisms, AList, Buckets}
--------------------------------------------------------------------
Model = new Type of HashTable
LeafTree = new Type of List
group = method()
group(Model) := M -> M.Group
buckets = method()
buckets(Model) := M -> M.Buckets
aList = (M,g) -> M.AList#g
model = method()
model(List,List,List) := (G,buckets,auts) -> (
modelAuts := hashTable for l in auts list (
l#0 => (hashTable for i to #G-1 list G#i => G#(l#1#i)));
AL := hashTable for g in G list (
g => apply(buckets, b->if member(g,b) then 1 else 0));
new Model from hashTable {
Group => G,
Buckets => buckets,
Automorphisms => modelAuts,
AList => AL
}
)
--ZZ/2 models
F0 = 0_(ZZ/2)
F1 = 1_(ZZ/2)
ZZ2 = {F0,F1}
--CFN
CFNmodel = model(ZZ2, {{F0}, {F1}}, {})
--ZZ/2 x ZZ/2 models
F00 = {0_(ZZ/2),0_(ZZ/2)}
F01 = {0_(ZZ/2),1_(ZZ/2)}
F10 = {1_(ZZ/2),0_(ZZ/2)}
F11 = {1_(ZZ/2),1_(ZZ/2)}
ZZ2ZZ2 = {F00,F01,F10,F11}
--Jukes-Cantor
JCmodel = model(ZZ2ZZ2, {{F00}, {F01,F10,F11}}, {
({1,2},{0,2,1,3}),
({1,3},{0,3,2,1}),
({2,1},{0,2,1,3}),
({2,3},{0,1,3,2}),
({3,1},{0,3,2,1}),
({3,2},{0,1,3,2})})
--Kimura 2-parameter
K2Pmodel = model(ZZ2ZZ2, {{F00}, {F01}, {F10,F11}}, {
({2,3},{0,1,3,2}),
({3,2},{0,1,3,2})})
--Kimura 3-parameter
K3Pmodel = model(ZZ2ZZ2, {{F00}, {F01}, {F10}, {F11}}, {})
qRing = method(Options=>{Variable=>null})
qRing(LeafTree,Model) := opts -> (T,M) -> qRing(#(leaves T),M,opts)
qRing(ZZ,Model) := opts -> (n,M) -> (
qList := leafColorings(n,M);
q := opts.Variable;
if q === null then q = getSymbol "q";
qRingFromList(qList,M,q)
)
qRingFromList = (qList,M,q) -> (
G := group M;
Ghash := hashTable apply(#G,i->(G#i=>i));
QQ(monoid[apply(qList, qcolors -> (qindex := apply(qcolors, c->Ghash#c); q_qindex))])
)
pRing = method(Options=>{Variable=>null})
pRing(LeafTree,Model) := opts -> (T,M) -> pRing(#(leaves T),M,opts)
pRing(ZZ,Model) := opts -> (n,M) -> (
G := group M;
pList := (n:0)..(n:#G-1);
p := opts.Variable;
if p === null then p = getSymbol "p";
QQ(monoid[apply(pList, pindex->p_pindex)])
)
fourierToProbability = method()
fourierToProbability(Ring,Ring,ZZ,Model) := (S,R,n,M) -> (
if not member(M, set{CFNmodel, JCmodel, K2Pmodel, K3Pmodel}) then
error "model must be CFNmodel, JCmodel, K2Pmodel, or K3Pmodel";
K := keys M;
G := M#(K_1);
qList := leafColorings(n,M);
Ghash := hashTable apply(#G,i->(G#i=>i));
varIndex := apply(qList, qcolors -> apply(qcolors, c->Ghash#c));
L := (n:0)..(n:#G-1);
Char := matrix{
{1,1,1,1},
{1,-1,1,-1},
{1,1,-1,-1},
{1,-1,-1,1}};
SubVars := for vi in varIndex list (
sum for i to #L-1 list (
s := product(n, j->(Char_(vi#j, (L#i)#j)));
if s>0 then S_i else -S_i
)
);
map(S,R,matrix{SubVars})
)
phyloToric42 = method(Options=>{QRing=>null})
phyloToric42(ZZ,List,Model) := opts -> (n,E,M) -> phyloToric42(leafTree(n,E),M,opts)
phyloToric42(Graph,Model) := opts -> (G,M) -> phyloToric42(leafTree G,M,opts)
phyloToric42(LeafTree,Model) := opts -> (T,M) -> (
A := phyloToricAMatrix(T,M);
S := if opts.QRing =!= null then opts.QRing else qRing(T,M);
toricMarkov(A,S)
)
phyloToricAMatrix = method()
phyloToricAMatrix(ZZ,List,Model) := (n,E,M) -> phyloToricAMatrix(leafTree(n,E),M)
phyloToricAMatrix(Graph,Model) := (G,M) -> phyloToricAMatrix(leafTree G,M)
phyloToricAMatrix(LeafTree,Model) := (T,M) -> (
ECs := edgeColorings(T,M);
A := for ec in ECs list flatten for g in ec list aList(M,g);
transpose matrix A
)
leafColorings = method()
leafColorings(LeafTree,Model) := (T,M) -> leafColorings(#(leaves T),M)
leafColorings(ZZ,Model) := (n,M) -> (
G := group M;
qList := toList(((n-1):0)..((n-1):(#G-1))); -- list of q variable indices
for qindex in qList list (
qcolors := apply(qindex,j->G#j);
append(qcolors,sum toList qcolors)
)
)
--List all consistent edge colorings of a tree
edgeColorings = (T,M) -> (
L := leavesList T;
Lhash := hashTable apply(#L, i->L#i=>i);
qList := leafColorings(T,M);
for qcolors in qList list for e in edges T list (
sum apply(toList e, l->qcolors#(Lhash#l))
)
)
-- Uses Local Structure of Invariants (Theorem 24 in Stumfels and Sullivant) to inductively determine the ideal of
-- phylogenetic invariants for any k-valent tree
-- with three types of generators: linear, quadratic and claw tree generators
phyloToricFP = method(Options=>{QRing=>null})
phyloToricFP(ZZ,List,Model) := opts -> (n,E,M) -> phyloToricFP(leafTree(n,E),M,opts)
phyloToricFP(LeafTree,Model) := opts -> (T,M) -> (
S := if opts.QRing =!= null then opts.QRing else qRing(#(leaves T),M);
Inv1 := phyloToricLinears(T,M,QRing=>S);
Inv2 := phyloToricQuads(T,M,QRing=>S);
Inv3 := phyloToricClaw(T,M,QRing=>S);
gensList := Inv1|Inv2|Inv3;
ideal gensList
)
phyloToricLinears = method(Options=>{QRing=>null,Random=>false})
phyloToricLinears(ZZ,List,Model) := opts -> (n,E,M) -> phyloToricLinears(leafTree(n,E),M,opts)
phyloToricLinears(LeafTree,Model) := opts -> (T,M) -> (
S := if opts.QRing =!= null then opts.QRing else qRing(T,M);
ECs := edgeColorings(T,M);
P := partition(i -> for g in ECs#i list aList(M,g), toList(0..#ECs-1));
gensList := flatten for p in values P list (
if #p < 2 then continue;
for j to #p-2 list sub(S_(p#j)-S_(p#(j+1)),S)
);
if not opts.Random then gensList else randomElement gensList
)
-- Produce the "edge invariants", quadratic invariants for each internal edge of T
phyloToricQuads = method(Options=>{QRing=>null,Random=>false})
phyloToricQuads(ZZ,List,Model) := opts -> (n,E,M) -> phyloToricQuads(leafTree(n,E),M,opts)
phyloToricQuads(LeafTree,Model) := opts -> (T,M) -> (
S := if opts.QRing =!= null then opts.QRing else qRing(T,M);
quadTemplates := apply(#(group M), g->({{g,g},{g,g}},{{g,g},{g,g}},{0,3,2,1}));
if opts.Random then quadTemplates = randomElement quadTemplates;
newl := symbol newl;
intEdges := internalEdges T;
if opts.Random then intEdges = randomElement intEdges;
gensList := flatten for e in intEdges list (
P := edgeCut(T,e,newl);
fillTemplates(T,M,S,P,quadTemplates,newl,opts.Random)
);
select(gensList, f -> f != 0_(ring f))
)
phyloToricClaw = method(Options=>{QRing=>null,Random=>false})
phyloToricClaw(ZZ,List,Model) := opts -> (n,E,M) -> phyloToricClaw(leafTree(n,E),M,opts)
phyloToricClaw(LeafTree,Model) := opts -> (T,M) -> (
S := if opts.QRing =!= null then opts.QRing else qRing(T,M);
l := first leaves T;
newl := symbol newl;
clawHash := new MutableHashTable; --store claw invariants to avoid recomputing
intVerts := (internalEdges T)|{set{l}};
if opts.Random then intVerts = randomElement intVerts;
gensList := flatten for e in intVerts list (
P := vertexCut(T,e,l,newl);
if not clawHash#?(#P,M) then clawHash#(#P,M) = clawInvariants(#P,M);
fillTemplates(T,M,S,P,clawHash#(#P,M),newl,opts.Random)
);
select(gensList, f -> f != 0_(ring f))
)
phyloToricRandom = method(Options=>{QRing=>null})
phyloToricRandom(ZZ,List,Model) := opts -> (n,E,M) -> phyloToricRandom(leafTree(n,E),M,opts)
phyloToricRandom(LeafTree,Model) := opts -> (T,M) -> (
S := if opts.QRing =!= null then opts.QRing else qRing(#(leaves T),M);
n := random 3;
gensList := if n == 0 then phyloToricLinears(T,M,QRing=>S,Random=>true)
else if n == 1 then phyloToricQuads(T,M,QRing=>S,Random=>true)
else phyloToricClaw(T,M,QRing=>S,Random=>true);
if #gensList > 0 then first gensList else phyloToricRandom(T,M,opts)
)
randomElement = L -> (
if #L == 0 then return {};
n := random(#L);
{L#n}
)
----------------------
--Auxilary functions for phyloToricFP
----------------------
--A function that takes an invariant on a small tree and extends it in all possible ways to a big tree.
--P is a graph splitting: the list of connected components after deleting the vertex or edge we are focused on.
--Temps is list of "templates" meaning invariants on small trees
fillTemplates = (T,M,S,P,temps,newl,rand) -> (
G := group M;
FCs := edgeColorings(T,M); --consistent colorings
qhash := hashTable apply(#FCs, i->FCs#i => S_i); --maps from consistent edge colorings to variables in the ring
n := #P;
cem := compositeEdgeMap(T,P,newl);
PFCs := apply(P, U->partitionedFCs(U,M,set{newl})); --a List of HashTables of colorings of the graph pieces
--print PFCs;
flatten for binom in temps list (
fbinom0 := flatten binom#0; --flat list of color indices for first monomial of a claw tree invariant
fbinom1 := flatten binom#1; --same for second monomial
PFCLists := apply(#fbinom0, j->(PFCs#(j%n))#(G#(fbinom0#j))); --for each entry of fbinom0(itself a list), the list of all coloring extensions
PFCi := (#(binom#2):0)..(toSequence apply(PFCLists,l->(#l-1))); --sequence of list indices for all combinations of coloring extensions
if rand then PFCi = toSequence randomElement PFCi;
newGens := for iList in PFCi list (
CList0 := apply(#iList, j->(PFCLists#j)#(iList#j)); --a list of color extensions for first monomial
CList1 := apply(binom#2, k->CList0#k); --permutations of the color extensions for second monomial
CList1 = apply(#fbinom1, j->(
c1 := fbinom1#j;
c0 := fbinom0#(binom#2#j);
if c1 != c0 then permuteColoring(CList1#j,{c1,c0},M) else CList1#j
));
CLists := (CList0,CList1);
--print CLists;
monom := apply(CLists, CList->product for i from 0 to #(binom#0)-1 list (
comcol := compositeColoring(cem, take(CList,{i*n,i*n+n-1}));
qhash#comcol
));
--print(monom#0,monom#1);
monom#0 - monom#1
);
ultimate(flatten,newGens)
)
)
--computes the invariants for a claw tree on k leaves.
--converts each binomial to a form used by the fillTemplates function.
--the values for k==3 and the four built-in models are hard-coded
clawInvariants = (k,M) -> (
if k == 3 and M === CFNmodel then return {};
if k == 3 and M === JCmodel then return {
{{{0,1,1},{1,0,1},{1,1,0}},{{0,0,0},{1,2,3},{1,2,3}},{0,4,8,3,1,2,6,7,5}}};
if k == 3 and M === K2Pmodel then return {
{{{1,0,1},{2,1,3},{2,2,0}},{{1,1,0},{2,0,2},{2,3,1}},{0,4,8,3,1,5,6,7,2}},
{{{0,1,1},{1,2,3},{2,0,2}},{{0,2,2},{1,0,1},{2,1,3}},{0,4,5,3,7,2,6,1,8}},
{{{0,1,1},{1,2,3},{2,2,0}},{{0,2,2},{1,1,0},{2,3,1}},{0,4,5,3,1,8,6,7,2}},
{{{0,1,1},{1,0,1},{2,2,0},{2,2,0}},{{0,0,0},{1,1,0},{2,3,1},{2,3,1}},{0,4,8,3,1,11,6,7,2,9,10,5}},
{{{0,1,1},{2,0,2},{2,2,0}},{{0,0,0},{2,1,3},{2,3,1}},{0,4,8,3,1,5,6,7,2}},
{{{0,1,1},{1,1,0},{2,0,2},{2,0,2}},{{0,0,0},{1,0,1},{2,1,3},{2,1,3}},{0,7,5,3,10,2,6,1,8,9,4,11}},
{{{0,2,2},{1,0,1},{2,2,0}},{{0,0,0},{1,2,3},{2,3,1}},{0,4,8,3,1,2,6,7,5}},
{{{0,2,2},{1,1,0},{2,0,2}},{{0,0,0},{1,2,3},{2,1,3}},{0,7,5,3,1,2,6,4,8}},
{{{0,2,2},{0,2,2},{1,0,1},{1,1,0}},{{0,0,0},{0,1,1},{1,2,3},{1,2,3}},{0,7,11,3,10,8,6,1,2,9,4,5}}};
if k == 3 and M === K3Pmodel then return {
{{{2,2,0},{2,3,1},{3,0,3},{3,1,2}},{{2,0,2},{2,1,3},{3,2,1},{3,3,0}},{0,7,11,3,10,8,6,1,5,9,4,2}},
{{{1,2,3},{2,3,1},{3,1,2}},{{1,3,2},{2,1,3},{3,2,1}},{0,4,8,3,7,2,6,1,5}},
{{{1,3,2},{2,2,0},{3,0,3}},{{1,2,3},{2,0,2},{3,3,0}},{0,4,8,3,7,2,6,1,5}},
{{{1,0,1},{1,3,2},{2,1,3},{2,2,0}},{{1,1,0},{1,2,3},{2,0,2},{2,3,1}},{0,7,11,3,10,8,6,1,5,9,4,2}},
{{{1,0,1},{2,2,0},{3,1,2}},{{1,1,0},{2,0,2},{3,2,1}},{0,7,5,3,1,8,6,4,2}},
{{{1,1,0},{2,3,1},{3,0,3}},{{1,0,1},{2,1,3},{3,3,0}},{0,7,5,3,1,8,6,4,2}},
{{{1,1,0},{1,3,2},{3,0,3},{3,2,1}},{{1,0,1},{1,2,3},{3,1,2},{3,3,0}},{0,7,11,3,10,8,6,1,5,9,4,2}},
{{{0,2,2},{1,3,2},{2,1,3},{3,0,3}},{{0,3,3},{1,2,3},{2,0,2},{3,1,2}},{0,4,8,3,1,11,6,10,2,9,7,5}},
{{{0,2,2},{2,3,1},{3,0,3}},{{0,3,3},{2,0,2},{3,2,1}},{0,4,8,3,7,2,6,1,5}},
{{{0,3,3},{2,2,0},{3,1,2}},{{0,2,2},{2,1,3},{3,3,0}},{0,4,8,3,7,2,6,1,5}},
{{{0,3,3},{1,3,2},{2,2,0},{3,2,1}},{{0,2,2},{1,2,3},{2,3,1},{3,3,0}},{0,7,5,3,10,2,6,1,11,9,4,8}},
{{{0,2,2},{1,0,1},{2,3,1},{3,1,2}},{{0,1,1},{1,3,2},{2,0,2},{3,2,1}},{0,10,5,3,7,2,6,4,11,9,1,8}},
{{{0,1,1},{1,2,3},{2,0,2}},{{0,2,2},{1,0,1},{2,1,3}},{0,4,8,3,7,2,6,1,5}},
{{{0,1,1},{1,3,2},{2,2,0}},{{0,2,2},{1,1,0},{2,3,1}},{0,7,5,3,1,8,6,4,2}},
{{{0,1,1},{1,2,3},{2,2,0},{3,1,2}},{{0,2,2},{1,1,0},{2,1,3},{3,2,1}},{0,4,11,3,1,8,6,10,5,9,7,2}},
{{{0,1,1},{1,3,2},{3,0,3}},{{0,3,3},{1,0,1},{3,1,2}},{0,4,8,3,7,2,6,1,5}},
{{{0,1,1},{1,2,3},{2,3,1},{3,0,3}},{{0,3,3},{1,0,1},{2,1,3},{3,2,1}},{0,7,5,3,10,2,6,1,11,9,4,8}},
{{{0,3,3},{1,1,0},{2,3,1},{3,1,2}},{{0,1,1},{1,3,2},{2,1,3},{3,3,0}},{0,4,8,3,1,11,6,10,2,9,7,5}},
{{{0,3,3},{1,1,0},{3,2,1}},{{0,1,1},{1,2,3},{3,3,0}},{0,4,8,3,7,2,6,1,5}},
{{{0,0,0},{0,3,3},{3,1,2},{3,2,1}},{{0,1,1},{0,2,2},{3,0,3},{3,3,0}},{0,7,11,3,10,8,6,1,5,9,4,2}},
{{{0,0,0},{1,1,0},{2,3,1},{3,2,1}},{{0,1,1},{1,0,1},{2,2,0},{3,3,0}},{0,4,8,3,1,11,6,10,2,9,7,5}},
{{{0,0,0},{2,3,1},{3,1,2}},{{0,1,1},{2,0,2},{3,3,0}},{0,7,5,3,1,8,6,4,2}},
{{{0,1,1},{2,2,0},{3,0,3}},{{0,0,0},{2,1,3},{3,2,1}},{0,7,5,3,1,8,6,4,2}},
{{{0,1,1},{1,1,0},{2,0,2},{3,0,3}},{{0,0,0},{1,0,1},{2,1,3},{3,1,2}},{0,7,5,3,10,2,6,1,11,9,4,8}},
{{{0,1,1},{0,3,3},{2,0,2},{2,2,0}},{{0,0,0},{0,2,2},{2,1,3},{2,3,1}},{0,7,11,3,10,8,6,1,5,9,4,2}},
{{{0,0,0},{1,3,2},{3,2,1}},{{0,2,2},{1,0,1},{3,3,0}},{0,7,5,3,1,8,6,4,2}},
{{{0,2,2},{1,0,1},{2,2,0},{3,0,3}},{{0,0,0},{1,2,3},{2,0,2},{3,2,1}},{0,4,8,3,1,11,6,10,2,9,7,5}},
{{{0,0,0},{1,3,2},{2,2,0},{3,1,2}},{{0,2,2},{1,1,0},{2,0,2},{3,3,0}},{0,7,5,3,10,2,6,1,11,9,4,8}},
{{{0,2,2},{1,1,0},{3,0,3}},{{0,0,0},{1,2,3},{3,1,2}},{0,7,5,3,1,8,6,4,2}},
{{{0,0,0},{1,3,2},{2,3,1},{3,0,3}},{{0,3,3},{1,0,1},{2,0,2},{3,3,0}},{0,4,11,3,1,8,6,10,5,9,7,2}},
{{{0,3,3},{1,0,1},{2,2,0}},{{0,0,0},{1,2,3},{2,3,1}},{0,4,8,3,7,2,6,1,5}},
{{{0,3,3},{1,1,0},{2,0,2}},{{0,0,0},{1,3,2},{2,1,3}},{0,7,5,3,1,8,6,4,2}},
{{{0,3,3},{1,1,0},{2,2,0},{3,0,3}},{{0,0,0},{1,2,3},{2,1,3},{3,3,0}},{0,10,5,3,7,2,6,4,11,9,1,8}},
{{{0,2,2},{0,3,3},{1,0,1},{1,1,0}},{{0,0,0},{0,1,1},{1,2,3},{1,3,2}},{0,7,11,3,10,8,6,1,5,9,4,2}}};
qList := leafColorings(k,M);
q := getSymbol "q";
R := qRingFromList(qList,M,q);
Igens := flatten entries gens phyloToric42(k,{},M,QRing=>R);
Igens = select(Igens, f-> 1 < first degree f);
for f in Igens list binomialTemplate(f,k,M)
)
--converts a binomial claw tree invariant into "template" form
binomialTemplate = (f,k,M) -> (
G := group M;
Ghash := hashTable apply(#G,i->(G#i=>i));
qList := leafColorings(k,M);
termsList := for t in exponents f list join toSequence apply(#t, i->toList (t#i:toList apply(qList#i,j->Ghash#j)));
d := #(termsList#0);
h := new MutableList;
for i to d*k-1 do (
c := position((0..d-1), j->(
n := j*k + i%k;
(not h#?n or h#n === null) and
aList(M,G#((termsList#1)#j#(i%k))) == aList(M,G#((termsList#0)#(i//k)#(i%k)))));
h#(c*k + i%k) = i;
);
h = toList h;
termsList|{h}
)
--e is distinguished in T (edges outputs the list of edges)
--eindex just tells you the first time a bool is true
--"or" is necessary because each edge has two ways to be named depending on which side of the partition the edge is located
--partition breaks up list based on value of function and outputs hash table with
--key color values list of edge coloring with that color e
partitionedFCs = (T,M,e) -> (
eindex := position(edges T, f->f==e or (leaves T)-f==e);
FCs := edgeColorings(T,M);
partition(fc->fc#eindex, FCs)
)
--produces a map from the set of edges of T to the edges of the decomposition at vertex l.
--L is a list of trees in the decomposition
compositeEdgeMap = (T,L,l) -> (
EL := apply(L, U->apply(edges U, e->orientEdge'(U,e,l)));
for e in edges T list (
(i,j) := (0,0);
while i < #EL do (
j = position(EL#i, f -> e==f or (leaves T)-e==f);
if j =!= null then break else i = i+1;
);
(i,j)
)
)
compositeColoring = (cem,L) -> apply(cem, ij->L#(ij#0)#(ij#1))
permuteColoring = (col,P,M) -> (
G := group M;
aut := M.Automorphisms#P;
apply(col, c->aut#c)
)
--------------------------
--LeafTree
--------------------------
leafTree = method()
leafTree(List,List) := (L,E) -> (
E = select(E, e->#e > 1 and #e < #L-1);
leafEdges := if #L > 2 then apply(L, i->{i}) else if #L == 2 then {{L#0}} else {};
E = E|leafEdges;
E = apply(E, e->if class e === Set then e else set e);
new LeafTree from {L, E}
)
leafTree(ZZ,List) := (n,E) -> leafTree(toList(0..n-1),E)
leafTree(Graph) := G -> (
if not isTree G then error "graph must be a tree";
L := select(vertexSet G, v->isLeaf(G,v));
E := for e in edges G list (
G' := deleteEdges(G,{toList e});
side := first connectedComponents G';
select(side, v->isLeaf(G,v))
);
leafTree(L,E)
)
edges(LeafTree) := T -> T#1
internalEdges = method()
internalEdges(LeafTree) := T -> select(T#1, e->#e > 1 and #e < #(T#0)-1)
leaves(LeafTree) := T -> set T#0
leavesList = method()
leavesList(LeafTree) := T -> T#0
vertices(LeafTree) := T -> vertices graph T
internalVertices = method()
internalVertices(LeafTree) := T -> (
select(vertices graph T, v->#v > 1)
)
graph(LeafTree) := opts -> T -> (
l := first elements leaves T;
E := apply(edges T, f->orientEdge(T,f,l));
E = E|{set{}};
P := poset(E,isSubset);
G := graph coveringRelations P;
newLabels := for v in vertexSet G list (
children := select(elements neighbors(G,v), w -> #w > #v);
children = apply(children, w -> (leaves T) - w);
if #v > 0 then set({v}|children) else set children
);
graph(newLabels, adjacencyMatrix G)
)
digraph(LeafTree,List) := opts -> (T,L) -> digraph(T,set L)
digraph(LeafTree,Set) := opts -> (T,u) -> (
E := apply(edges T, e->if any(elements u, f->isSubset(e,f)) then (leaves T) - e else e);
E = E|{set{}};
P := poset(E,isSubset);
G := graph coveringRelations P;
dirE := {};
newLabels := for v in vertexSet G list (
children := select(elements neighbors(G,v), w -> #w > #v);
dirE = dirE|apply(children, w->{v,w});
children = apply(children, w -> (leaves T) - w);
if #v > 0 then set({v}|children) else set children
);
D := digraph(vertices G, dirE);
digraph(newLabels, adjacencyMatrix D)
)
Set == Set := (s,t) -> s === t
LeafTree == LeafTree := (S,T) -> (
if leaves S != leaves T then return false;
l := first leavesList S;
ES := set apply(edges S, e->orientEdge(S,e,l));
ET := set apply(edges T, e->orientEdge(T,e,l));
ES == ET
)
AHU := (G,v) -> (
chil := children(G,v);
chilAHU := flatten sort apply(elements chil, w->AHU(G,w));
{1}|chilAHU|{0}
)
isIsomorphicRooted = method()
isIsomorphicRooted(LeafTree,List,LeafTree,List) := (T1,v1,T2,v2) -> (
isIsomorphicRooted(T1,set v1,T2,set v2)
)
isIsomorphicRooted(LeafTree,Set,LeafTree,Set) := (T1,v1,T2,v2) -> (
G1 := digraph(T1,v1);
G2 := digraph(T2,v2);
AHU(G1,v1) == AHU(G2,v2)
)
isIsomorphic = method()
isIsomorphic(LeafTree,LeafTree) := (T1,T2) -> (
if #(leaves T1) != #(leaves T2) or #(edges T1) != #(edges T2) then return false;
C1 := center graph T1;
C2 := center graph T2;
if #C1 != #C2 then return false;
for v1 in C1 do for v2 in C2 do (
if isIsomorphicRooted(T1,v1,T2,v2) then return true;
);
false
)
--splits tree T at edge e into a list of two trees
edgeCut = method()
edgeCut(LeafTree,List,Thing) := (T,e,newl) -> edgeCut(T,set e,newl)
edgeCut(LeafTree,Set,Thing) := (T,e,newl) -> (
Lpart := {e, leaves(T) - e};
apply(Lpart, P->leafTree((toList P)|{newl}, edgeSelect(T,P)))
)
edgeSelect = (T,e) -> (
for f in internalEdges T list (
if f==e then continue
else if isSubset(f,e) then f
else if isSubset((leaves T) - f,e) then (leaves T) - f
else continue
)
)
--gives the side of the partition representing edge e that contains leaf l
orientEdge = (T,e,l) -> if member(l,e) then e else (leaves T) - e
--gives the side of the partition representing edge e that does not contain leaf l
orientEdge' = (T,e,l) -> if member(l,e) then (leaves T) - e else e
--splits tree T at vertex v into a list of trees
--v is specified as the vertex incident to edge e away from leaf l
vertexCut = method()
vertexCut(LeafTree,List,Thing,Thing) := (T,e,l,newl) -> vertexCut(T,set e,l,newl)
vertexCut(LeafTree,Set,Thing,Thing) := (T,e,l,newl) -> (
e = orientEdge(T,e,l);
E := apply(edges T, f->orientEdge(T,f,l));
E = E|{set{}};
P := poset(E,isSubset);
G := graph coveringRelations P;
Lpart := select(elements neighbors(G,e), w -> #w > #e);
Lpart = apply(Lpart, w -> (leaves T) - w)|{e};
apply(Lpart, P->leafTree((toList P)|{newl}, edgeSelect(T,P)))
)
edgeContract = method()
edgeContract(LeafTree,List) := (T,e) -> edgeContract(T,set e)
edgeContract(LeafTree,Set) := (T,e) -> (
L := leaves T;
E := select(edges T, f->(f != e and L - f != e));
if #e == #L-1 then L = e;
if #e == 1 then L = L - e;
leafTree(toList L, E)
)
labeledTrees = method()
labeledTrees(ZZ) := n -> (
f := L -> (
P := setPartitions L;
select(P, p -> #p > 1)
);
L := toList (1..n-1);
apply(buildBranches(L,f), T -> leafTree(n,T))
)
labeledBinaryTrees = method()
labeledBinaryTrees(ZZ) := n -> (
f := L -> for s in subsets drop(L,1) list (
if #s == 0 then continue;
{s, toList (set L - set s)}
);
L := toList (1..n-1);
apply(buildBranches(L,f), T -> leafTree(n,T))
)
rootedTrees = method()
rootedTrees(ZZ) := n -> (
f := L -> for p in partitions(#L) list (
if #p == 1 then continue;
k := 0;
for s in p list (
k = k+s;
take(L,{k-s,k-1})
)
);
L := toList (1..n-1);
apply(buildBranches(L,f), T -> leafTree(n,T))
)
rootedBinaryTrees = method()
rootedBinaryTrees(ZZ) := n -> (
f := L -> for i from 1 to (#L)//2 list {take(L,i), take(L,i-#L)};
L := toList (1..n-1);
apply(buildBranches(L,f), T -> leafTree(n,T))
)
unlabeledTrees = method()
unlabeledTrees(ZZ) := n -> (
rooted := rootedTrees n;
trees := new MutableList;
for T in rooted do (
if not any(trees, S->isIsomorphic(S,T)) then trees#(#trees) = T
);
toList trees
)
--lists all partitions of a set or list of distinct elements
setPartitions = method()
setPartitions(Set) := S -> setPartitions(toList S)
setPartitions(List) := L -> (
Lhash := new HashTable from apply(#L, i->(L#i => i));
pList := toList (#L : 0);
sps := {{L}};
i := #L-1;
while i > 0 do (
if any(i, j -> pList#j >= pList#i) then (
pList = take(pList, i)|{pList#i + 1}|toList(#L-i-1:0);
part := values partition(l->pList#(Lhash#l), L);
sps = append(sps, part);
i = #L-1;
)
else i = i-1;
);
sps
)
--recursive function for building rooted trees.
--takes a leaf set L and function f that lists how leaves can be
--partitioned at a node, and lists all possible edge sets.
buildBranches = (L,f) -> (
Trees := for p in f(L) list (
p = select(p, s -> #s > 1);
newTrees := {p};
for E in p do (
branches := buildBranches(E,f);
newTrees = flatten apply(newTrees, T->apply(branches, B->T|B));
);
newTrees
);
flatten Trees
)
--------------------------
--Secants and Joins
--------------------------
--list the monomials in ring R corresponding to the columns of matrix A
imageMonomials = method()
imageMonomials(Ring,Matrix) := (R,A) -> (
M := for i from 0 to (numcols A) - 1 list (
vect := flatten entries A_i;
product apply(#vect, j->R_j^(vect#j))
);
matrix {M}
)
--computes the nth secant of ideal I using elimination
--if degree d is specified, then the generators up to degree d will be computed (this is much faster)
secant = method(Options=>{DegreeLimit => {}})
secant(Ideal,ZZ) := opts -> (I,n) -> joinIdeal(toList (n:I), opts)
joinIdeal = method(Options=>{DegreeLimit => {}})
joinIdeal(Ideal,Ideal) := opts -> (I,J) -> joinIdeal({I,J},opts)
joinIdeal(List) := opts -> L -> (
R := ring first L;
k := numgens R;
n := #L;
T := R;
for i from 0 to n-1 do T = T**R;
T = newRing(T, MonomialOrder=>Eliminate(n*k));
Jlinears := apply(k,j->T_(n*k+j) - sum(n,i->T_(i*k+j)));
Js := apply(n, i->sub(L#i,(vars T)_(toList(i*k..(i+1)*k-1))));
J := sum(Js) + ideal Jlinears;
d := opts.DegreeLimit;
GB := gb(J, DegreeLimit=>join((n+1):d));
J = selectInSubring(1,gens GB);
ideal sub(J,matrix{toList (n*k:0_R)}|(vars R))
)
--Randomized algorithm for affine dimension of kth secant of variety defined by matrix A
toricSecantDim = method()
toricSecantDim(Matrix,ZZ) := (A,k) -> (
kk := ZZ/32003;
n := numrows A;
A = homogenizeMatrix A;
randPoints := apply(k,i->random(kk^1,kk^(n+1)));
x := symbol x;
R := kk[x_0..x_n];
J := jacobian imageMonomials(R,A);
tSpace := matrix apply(randPoints, p->{sub(J,p)});
rank tSpace
)
toricJoinDim = method()
toricJoinDim(Matrix,Matrix) := (A,B) -> toricJoinDim({A,B})
toricJoinDim(List) := L -> (
kk := ZZ/32003;
k := #L;
n := L/numrows;
L = L/homogenizeMatrix;
randPoints := apply(#L, i->random(kk^1,kk^(n#i+1)));
x := symbol x;
R := apply(#L, i->kk[x_0..x_(n#i)]);
J := apply(#L, i->jacobian imageMonomials(R#i,L#i));
tSpace := matrix apply(#L, i->{sub(J#i,randPoints#i)});
rank tSpace
)
homogenizeMatrix = A -> (
n := numrows A;
N := numcols A;
colSums := matrix{toList(n:1)}*A;
d := max flatten entries colSums;
A||(matrix{toList(N:d)}-colSums)
)
------------------------------------------
-- Documentation
------------------------------------------
beginDocumentation()
-------------------------------
-- PhylogeneticTrees
-------------------------------
doc///
Key
PhylogeneticTrees
Headline
a package to compute phylogenetic invariants associated to group-based models
Description
Text
{\em PhylogeneticTrees} is a package for phylogenetic algebraic geometry. This
package calculates generating sets for phylogenetic ideals and their joins and
secants. Additionally, the package computes lower bounds for the dimensions
of secants and joins of phylogenetic ideals.
This package handles a class of commonly used tree-based Markov models called
group-based models. These models are subject to the Fourier-Hadamard coordinate
transformation, which make the parametrizations monomial and the ideals toric.
See the following for more details: [1] and [2].
For these models, the PhylogeneticTrees package includes two methods for computing
a generating set for ideals of phylogenetic invariants. The first method calls
@TO FourTiTwo@ to compute the generating set of the toric ideal. The second
implements a theoretical construction for inductively determining the ideal of
phylogenetic invariants for any $k$-valent tree from the $k$-leaf claw tree as
described in Theorem 24 of [3].
This package also handles the joins and secants of these ideals by implementing
the elimination method described in [4].
In cases where computing a generating set for a join or secant ideal is infeasible,
the package provides a probabilistic method, based on Terracini’s Lemma, to compute
a lower bound on the dimension of a join or secant ideal.
{\em References:}
[1] S.N. Evans and T.P. Speed, it Invariants of some probability models used in phylogenetic inference, {\em Ann. Statist.} 21 (1993), no. 1, 355–377, and
[2] L. Székely, P.L. Erdös, M.A. Steel, and D. Penny, A Fourier inversion formula for evolutionary trees, {\em Applied Mathematics Letters} 6 (1993), no. 2, 13–17.
[3] Bernd Sturmfels and Seth Sullivant, Toric ideals of phylogenetic invariants, {\em J. Comp. Biol.} 12 (2005), no. 2, 204–228.
[4] Bernd Sturmfels and Seth Sullivant, Combinatorial secant varieties, {\em Quarterly Journal of Pure and Applied Mathematics} 2 (2006), 285–309.
///
-------------------------------
-- Phylogenetic invariants
-------------------------------
--phyloToricFP
doc///
Key
phyloToricFP
(phyloToricFP,ZZ,List,Model)
(phyloToricFP,LeafTree,Model)
Headline
compute the invariants of a group-based phylogenetic model with toric fiber products
Usage
phyloToricFP(T,M)
phyloToricFP(n,E,M)
Inputs
T:LeafTree
n:ZZ
the number of leaves
E:LeafTree
the internal edges of the tree, given by one part of the bipartition on leaves
M:Model
Outputs
:Ideal
Description
Text
This function computes the invariants of a group-based phylogenetic
tree model based on Theorem 24 of the paper Toric Ideals of
Phylogenetic Invariants by Sturmfels and Sullivant.
Invariants are formed in three different ways. The linear and
quadratic invariants are computed as in @TO phyloToricLinears@ and
@TO phyloToricQuads@ respectively. Finally higher degree invariants
are built using a toric fiber product construction from the invariants of
claw trees.
Example
T = leafTree(4, {{0,1}})
phyloToricFP(T, CFNmodel)
SeeAlso
phyloToric42
///
-------------------------------
--phyloToric42
doc///
Key
phyloToric42
(phyloToric42,ZZ,List,Model)
(phyloToric42,Graph,Model)
(phyloToric42,LeafTree,Model)
Headline
compute the invariants of a group-based phylogenetic model with 4ti2
Usage
phyloToric42(n,E,M)
phyloToric42(G,M)
phyloToric42(T,M)
Inputs
T:LeafTree
n:ZZ
the number of leaves
E:LeafTree
the internal edges of the tree, given by one part of the bipartition on leaves
G:Graph
a tree
M:Model
Outputs
:Ideal
Description
Text
This function computes the invariants of a group-based phylogenetic
tree model by computing the transpose of the matrix
that encodes the defining monomial map and then using the function toricMarkov of the
@TO FourTiTwo@ package.
Example
T = leafTree(4, {{0,1}})
phyloToric42(T, CFNmodel)
SeeAlso
phyloToricFP
///
-------------------------------
-- phyloToricLinears
doc///
Key
phyloToricLinears
(phyloToricLinears,LeafTree,Model)
(phyloToricLinears,ZZ,List,Model)
[phyloToricLinears,Random]
Headline
compute the linear invariants of a group-based phylogenetic model
Usage
phyloToricLinears(T,M)
phyloToricLinears(n,E,M)
Inputs
T:LeafTree
n:ZZ
the number of leaves
E:List
the internal edges of the tree, given by one part of the bipartition on leaves
M:Model
Outputs
:List
a generating set of the linear invariants
Description
Text
For models such as Jukes-Cantor (@TO "JCmodel"@) and Kimura 2-parameter (@TO "K2Pmodel"@),
multiple variables in the Fourier coordinates may map to the same monomial under the
monomial map that defines the toric variety of the model. These equivalencies give rise
to linear relations in the space of Fourier coordinates.
The number of linear invariants is the codimension of the smallest linear subspace in
which the toric variety of the model is contained.
The optional argument @TO QRing@ can be passed the ring of Fourier coordinates. Otherwise
the function will create a new ring.
Example
T = leafTree(3,{})
S = qRing(T, K2Pmodel)
phyloToricLinears(T, K2Pmodel, QRing=>S)
SeeAlso
phyloToricFP
phyloToric42
phyloToricQuads
///
-------------------------------
-- phyloToricQuads
doc///
Key
phyloToricQuads
(phyloToricQuads,LeafTree,Model)
(phyloToricQuads,ZZ,List,Model)
[phyloToricQuads,Random]
Headline
compute the quadratic invariants of a group-based phylogenetic model
Usage
phyloToricQuads(T,M)
phyloToricQuads(n,E,M)
Inputs
T:LeafTree
n:ZZ
the number of leaves
E:List
the internal edges of the tree, given by one part of the bipartition on leaves
M:Model
Outputs
:List
a generating set of the quadratic invariants
Description
Text
The quadratic invariants are also referred to as the edge invariants of the model.
Each Fourier coordinate corresponds to a consistent coloring of the edges of tree {\tt T}.
For any given internal edge {\tt e} of {\tt T}, the consistent colorings can be obtained by
coloring two smaller graphs and gluing them along {\tt e}. This corresponds to a fiber
product on the corresponding toric varieties. The quadratic invariants naturally
arise from this process by gluing a pair of colorings of one small graph to a pair of
colorings of the other small graph in two different ways.
The optional argument @TO QRing@ can be passed the ring of Fourier coordinates. Otherwise
the function will create a new ring.
Example
T = leafTree(4,{{0,1}})
S = qRing(T, CFNmodel)
phyloToricQuads(T, CFNmodel, QRing=>S)
SeeAlso
phyloToricFP
phyloToric42
phyloToricLinears
///
-------------------------------
--phyloToricRandom
doc///
Key
phyloToricRandom
(phyloToricRandom,ZZ,List,Model)
(phyloToricRandom,LeafTree,Model)
Headline
compute a random invariant of a group-based phylogenetic model
Usage
phyloToricRandom(T,M)
phyloToricRandom(n,E,M)
Inputs
T:LeafTree
n:ZZ
the number of leaves
E:LeafTree
the internal edges of the tree, given by one part of the bipartition on leaves
M:Model
Outputs
:RingElement
a randomly selected binomial invariant
Description
Text
This function computes a random invariant of a group-based phylogenetic
tree model using the toric fiber product structure.
With equal probability the algorithm decides to return a linear,
quadratic, or higher degree binomial. It then selects one of these
at random (but uniformity is not guaranteed).
This is a much more efficient way to produce a single generator than listing all of them,
which is useful for Monte Carlo random walk algorithms.
Example
phyloToricRandom(4,{{0,1}},CFNmodel)
Caveat
We currently do not guarantee a uniform distribution on the generators, even
after the choice of linear, quadratic or higher degree.
SeeAlso
phyloToricFP
///
-------------------------------
-- phyloToricAMatrix
doc///
Key
phyloToricAMatrix
(phyloToricAMatrix,LeafTree,Model)
(phyloToricAMatrix,Graph,Model)
(phyloToricAMatrix,ZZ,List,Model)
Headline
construct the design matrix of a group-based phylogenetic model
Usage
phyloToricAMatrix(T,M)
phyloToricAmatrix(G,M)
phyloToricAMatrix(n,E,M)
Inputs
T:LeafTree
G:Graph
a tree
n:ZZ
the number of leaves
E:List
the internal edges of the tree, given by the half the partition on leaves
M:Model
Outputs
:Matrix
whose columns parametrize the toric variety
Description
Example
phyloToricAMatrix(4, {{1, 2}},CFNmodel)
SeeAlso
///
-------------------------------
-- qRing
doc///
Key
qRing
(qRing,ZZ,Model)
(qRing,LeafTree,Model)
[qRing,Variable]
Headline
construct the ring of Fourier coordinates
Usage
qRing(T,M)
qRing(n,M)
Inputs
T:LeafTree
n:ZZ
the number of leaves
M:Model
Outputs
:Ring
of Fourier coordinates
Description
Text
The Fourier coordinates for a phylogenetic tree model have one coordinate for each consistent coloring
of the tree {\tt T}. A consistent coloring is an assignment of one of the group elements of the model {\tt M} to each of
the leaves of {\tt T} such that the sum of all the group elements assigned is $0$.
Each variable of the ring is indexed by a sequence representing a consistent coloring with each element of the group
represented by an integer between $0$ and $m-1$ where $m$ is the order of the group.
A variable name for the ring can be passed using the optional argument {\tt Variable}.
Otherwise the symbol {\tt q} is used.
Example
qRing(4,CFNmodel)
qRing(3,JCmodel)
SeeAlso
pRing
leafColorings
///
-------------------------------
-- leafColorings
doc///
Key
leafColorings
(leafColorings,ZZ,Model)
(leafColorings,LeafTree,Model)
Headline
list the consistent colorings of a tree
Usage
leafColorings(T,M)
leafColorings(n,M)
Inputs
T:LeafTree
n:ZZ
the number of leaves
M:Model
Outputs
:List
the consistent colorings of the tree
Description
Text
This function outputs a list of all consistent colorings of the leaves of tree {\tt T}.
That is all sequences $(g_1,\ldots,g_n)$ such that $g_1+\cdots +g_n = 0$ where each $g_i$ is an
element of the group associated to the model {\tt M}, and {\tt n} is the number of leaves of the tree.
These correspond the set of subscripts of the variables in the ring output by @TO qRing@,
and appear in the same order.
Example
leafColorings(4,CFNmodel)
leafColorings(3,JCmodel)
SeeAlso
qRing
///
-------------------------------
-- pRing
doc///
Key
pRing
(pRing,ZZ,Model)
(pRing,LeafTree,Model)
[pRing,Variable]
Headline
construct the ring of probability coordinates
Usage
pRing(T,M)
pRing(n,M)
Inputs
T:LeafTree
n:ZZ
the number of leaves
M:Model
Outputs
:Ring
of probability coordinates
Description
Text
The probability coordinates for a phylogenetic tree model have one coordinate for each possible outcome of
the model. A possible outcome is any labeling of the leaves of the tree by elements of the group $G$ of the
model. Thus the number of coordinates is $|G|^n$ where $n$ is the number of leaves.
A variable name for the ring can be passed using the optional argument {\tt Variable}.
Otherwise the symbol {\tt p} is used.
Example
pRing(4,CFNmodel)
pRing(3,JCmodel)
SeeAlso
qRing
///
-------------------------------
-- QRing
doc///
Key
QRing
[phyloToric42,QRing]
[phyloToricFP,QRing]
[phyloToricLinears,QRing]
[phyloToricQuads,QRing]
[phyloToricRandom,QRing]
Headline
optional argument to specify Fourier coordinate ring
Description
Text
For any of the functions that produce phylogenetic invariants in the ring of Fourier coordinates,
the ring can be specified with this optional argument. If {\tt null} is passed then a new ring
of Fourier coordinates will be created.
The ring passed can be any polynomial ring with sufficiently many variables. The sufficient number
is $k = |G|^{n-1}$ where $G$ is the group of labels used by the model, and $n$ is the number of leaves of
the phylogenetic tree. The ring may have more than $k$ variables, in which case only the first $k$ will be used.
Example
T = leafTree(4,{{0,1}})
phyloToricFP(T,CFNmodel)
S = QQ[a..h]
phyloToricFP(T,CFNmodel,QRing=>S)
///
-------------------------------
-- fourierToProbability
doc///
Key
fourierToProbability
(fourierToProbability,Ring,Ring,ZZ,Model)
Headline
map from Fourier coordinates to probablity coordinates
Usage
fourierToProbability(S,R,n,M)
Inputs
S:Ring
of probability coordinates
R:Ring
of Fourier coordinates
n:ZZ
the number of leaves
M:Model
Outputs
:RingMap
from Fourier coordinates to probablity coordinates
Description
Text
This function creates a ring map from the ring of Fourier coordinates to the ring of probability coordinates,
for the four predefined models, @TO "CFNmodel"@, @TO "JCmodel"@, @TO "K2Pmodel"@ or @TO "K3Pmodel"@. It will not work with
user-defined models.
The ring of probability coordinates must have at least $|G|^n$ variables where $G$ is the group
associated to the model. The ring of Fourier coordinates must have at least $|G|^{(n-1)}$ variables.
Example
M = CFNmodel;
S = pRing(3,M)
R = qRing(3,M)
m = fourierToProbability(S,R,3,M)
SeeAlso
pRing
qRing
///
------------------------------
--Models
------------------------------
--CFNmodel
doc///
Key
"CFNmodel"
Headline
the model corresponging to the Cavender-Farris-Neyman model or binary Jukes Cantor
Description
Text
The Cavender-Farris-Neyman (CFN) Model is a Markov model of base substitution. It also known as the binary Jukes-Cantor model.
It assumes the root distribution vectors describe all bases occurring uniformly in the ancestral sequence.
It also assumes that the rate of all specific base changes is the same.
The transition matrix has the form
$$\begin{pmatrix} \alpha&\beta\\
\beta&\alpha \end{pmatrix}$$
SeeAlso
Model
"JCmodel"
"K2Pmodel"
"K3Pmodel"
///
--------------------------
--JCmodel
doc///
Key
"JCmodel"
Headline
the model corresponding to the Jukes Cantor model
Description
Text
The Jukes-Cantor (JK) Model is a Markov model of base substitution.
It assumes the root distribution vectors describe all bases occurring uniformly in the ancestral sequence.
It also assumes that the rate of all specific base changes is the same.
Thus the rates of bases changes A-G, A-T and A-C are the same.
The transition matrix has the form
$$\begin{pmatrix} \alpha&\beta&\beta&\beta\\
\beta&\alpha&\beta&\beta\\
\beta&\beta&\alpha&\beta\\
\beta&\beta&\beta&\alpha \end{pmatrix}$$
SeeAlso
Model
"CFNmodel"
"K2Pmodel"
"K3Pmodel"
///
---------------------------
--K2Pmodel
doc///
Key
"K2Pmodel"
Headline
the model corresponging to the Kimura 2-parameter model
Description
Text
The Kimura 2-parameter (K2P) Model is a Markov model of base substitution. It assumes the root distribution vectors describe
all bases occurring uniformly in the ancestral sequence. It allows different probabilities of transitions and transversions.
This means that the rate of base changes A-C and A-T are the same (transversions), and the rate of
base change A-G can differ from the other two (transitions).
The transition matrix has the form
$$\begin{pmatrix} \alpha&\gamma&\beta&\beta\\
\gamma&\alpha&\beta&\beta\\
\beta&\beta&\alpha&\gamma\\
\beta&\beta&\gamma&\alpha \end{pmatrix}$$
SeeAlso
Model
"CFNmodel"
"JCmodel"
"K3Pmodel"
///
-------------------------------
--K3Pmodel
doc///
Key
"K3Pmodel"
Headline
the model corresponging to the Kimura 3-parameter model
Description
Text
The Kimura 3-parameter (K3P) Model is a Markov model of base substitution.
It assumes the root distribution vectors describe all bases occurring uniformly in the ancestral sequence.
It allows different probabilities of the base changes A-G, A-C and A-T.
This is the most general group based model on group $(\mathbb{Z}/2\mathbb{Z})^2$.
The transition matrix has the form
$$\begin{pmatrix} \alpha&\gamma&\beta&\delta\\
\gamma&\alpha&\delta&\beta\\
\beta&\delta&\alpha&\gamma\\
\delta&\beta&\gamma&\alpha \end{pmatrix}$$
SeeAlso
Model
"CFNmodel"
"JCmodel"
"K2Pmodel"
///
-------------------------------
-- Secants and Joins
-------------------------------
-- secant
doc///
Key
secant
(secant,Ideal,ZZ)
[secant,DegreeLimit]
Headline
compute the secant of an ideal
Usage
secant(I,n)
Inputs
I:Ideal
k:ZZ
order of the secant
Outputs
:Ideal
the {\tt k}th secant of {\tt I}
Description
Text
This function computes the {\tt k}th secant of {\tt I} by constructing the abstract secant and then projecting with elimination.
Here the {\tt k}th secant means the join of {\tt k} copies of {\tt I}.
Setting {\tt k} to 1 gives the dimension of the ideal, while 2 is the usual secant, and higher
values correspond to higher order secants.
Setting the optional argument @TO DegreeLimit@ to {\tt \{d\} } will produce only the generators
of the secant ideal up to degree {\tt d}.
This method is general and will work for arbitrary polynomial ideals, not just phylogenetic ideals.
Example
R = QQ[a..d]
I = ideal {a^2-b,a^3-c,a^4-d}
secant(I,2)
SeeAlso
joinIdeal
///
-------------------------------
-- joinIdeal
doc///
Key
joinIdeal
(joinIdeal,Ideal,Ideal)
(joinIdeal,List)
[joinIdeal,DegreeLimit]
Headline
compute the join of several ideals
Usage
joinIdeal(I,J)
joinIdeal L
Inputs
I:Ideal
J:Ideal
L:List
of ideals in the same ring
Outputs
:Ideal
the join of the input ideals
Description
Text
This function computes the ideal of the join by constructing the abstract join and then projecting with elimination.
Setting the optional argument @TO DegreeLimit@ to {\tt \{d\} } will produce only the generators
of the join ideal up to degree {\tt d}.
This method is general and will work for arbitrary polynomial ideals, not just phylogenetic ideals.
Example
R = QQ[a,b,c,d]
I = ideal {a-d,b^2-c*d}
J = ideal {a,b,c}
joinIdeal(I,J)
SeeAlso
secant
///
-------------------------------
-- toricSecantDim
doc///
Key
toricSecantDim
(toricSecantDim,Matrix,ZZ)
Headline
dimension of a secant of a toric variety
Usage
toricSecantDim(A,k)
Inputs
A:Matrix
the A-matrix of a toric variety
k:ZZ
order of the secant
Outputs
:ZZ
the dimension of the {\tt k}th secant of variety defined by matrix {\tt A}
Description
Text
A randomized algorithm for computing the affine dimension of a secant of a toric variety
using Terracini's Lemma.
Here the {\tt k}th secant means the join of {\tt k} copies of {\tt I}.
Setting {\tt k} to 1 gives the dimension of the ideal, while 2 is the usual secant, and higher
values correspond to higher order secants.
The matrix {\tt A} defines a parameterization of the variety. The algorithm chooses {\tt k}
vectors of parameter values at random from a large finite field. The dimension of the sum of the tangent spaces
at those points is computed.
This algorithm is much much faster than computing the secant variety.
Example
A = matrix{{4,3,2,1,0},{0,1,2,3,4}}
toricSecantDim(A,1)
toricSecantDim(A,2)
toricSecantDim(A,3)
toricSecantDim(A,4)
SeeAlso
toricJoinDim
secant
///
-------------------------------
-- toricJoinDim
doc///
Key
toricJoinDim
(toricJoinDim,Matrix,Matrix)
(toricJoinDim,List)
Headline
dimension of a join of toric varieties
Usage
toricJoinDim(A,B)
toricJoinDim L
Inputs
A:Matrix
the A-matrix of a toric variety
B:Matrix
the A-matrix of a toric variety
L:List
of A-matrices of toric varieties
Outputs
:ZZ
the dimension of the join of the toric varieties defined by the matrices
Description
Text
A randomized algorithm for computing the affine dimension of a join of toric varieties
using Terracini's Lemma.
Each input matrix defines a parameterization of the variety. For each variety, a vector of parameter values
is chosen at random from a large finite field. The dimension of the sum of the tangent spaces
at those points is computed.
This algorithm is much much faster than computing the join variety.
Example
A = matrix{{4,3,2,1,0},{0,1,2,3,4}}
B = matrix{{1,1,1,1,1}}
toricJoinDim(A,B)
toricJoinDim(B,B)
Caveat
All input matrices must have the same number of columns.
SeeAlso
toricSecantDim
joinIdeal
///
-------------------------------
-- Model functionality
-------------------------------
-- Model
doc///
Key
Model
Headline
a group-based model
Description
Text
A phylogenetic tree model on tree $T$ has outcomes that are described by assigning each leaf of the tree
any label from a particular set (typically the label set is the set of DNA bases, \{A,T,C,G\}).
The probability of a certain assignment of labels depends on transition probabilities between each ordered pair of labels.
These transition probabilities are the parameters of the model.
In a group based model, the label set is a group $G$ (typically $\mathbb{Z}/2$ or $(\mathbb{Z}/2)^2$), and the transition
probability for a pair $(g,h)$ depends only on $h-g$. This reduces the number of parameters from $|G|^2$ to $|G|$.
Depending on the model, further identifications of parameters are imposed.
An object of class @TO Model@ stores the information about a group-based model required to
compute phylogenetic invariants.
This information includes the elements of the group, how those elements are partitioned, and a set of
automorphisms of the group that preserve the partitions.
There are four built-in models: Cavender-Farris-Neyman or binary model (@TO "CFNmodel"@);
Jukes-Cantor model (@TO "JCmodel"@); Kimura 2-parameter model (@TO "K2Pmodel"@);
and Kimura 3-parameter model (@TO "K3Pmodel"@). Other models can be constructed with @TO model@.
Example
M = CFNmodel
T = leafTree(3,{})
phyloToricAMatrix(T,M)
SeeAlso
model
///
-------------------------------
-- model
doc///
Key
model
(model,List,List,List)
Headline
construct a Model
Usage
model(G,B,aut)
Inputs
G:List
the group elements
B:List
of lists of which group elements have identified parameters
aut:List
of pairs, assigning pairs of identified group elements to automorphisms of the group that switch the pair
Outputs
:Model
Description
Text
The elements of {\tt G} must have an addition operation meaning that if two elements $g, h \in {\tt G$}, then $g+h$ must work. The usual choices for {\tt G} are the list of elements of
$\mathbb{Z}/2$ or $(\mathbb{Z}/2)^2$.
Example
(a,b) = (0_(ZZ/2),1_(ZZ/2))
G = {{a,a}, {a,b}, {b,a}, {b,b}}
Text
The elements of {\tt B} are lists of the elements of {\tt G} with the same parameter value.
In the following example, the first two elements of {\tt G} receive distinct parameters, while the last two share a parameter.
This is precisely the Kimura 2-parameter model.
Example
B = {{G#0}, {G#1}, {G#2,G#3}}
Text
Finally, for every ordered pair of group elements sharing a parameter, {\tt aut} must provide an automorphism of the group
that switches those two group elements. In {\tt aut} all of the group elements are identified by their index in $G$,
and an automorphism is given by a list of permuted index values.
In our example, the pairs requiring an automorphism are {\tt \{2,3\}} and {\tt \{3,2\}}.
Example
aut = {({2,3}, {0,1,3,2}),
({3,2}, {0,1,3,2})}
model(G,B,aut)
SeeAlso
Model
///
-------------------------------
-- group
doc///
Key
group
(group,Model)
Headline
the group of a Model
Usage
group M
Inputs
M:Model
Outputs
:List
of group elements
Description
Text
Every group-based phyogenetic model has a finite group associated to it. This function
returns the group, represented as a list of elements.
Example
M = K3Pmodel
G = group M
SeeAlso
Model
///
-------------------------------
-- buckets
doc///
Key
buckets
(buckets,Model)
Headline
the equivalence classes of group elements of a Model
Usage
buckets M
Inputs
M:Model
Outputs
:List
of lists of group elements
Description
Text
Every group-based phyogenetic model has a finite group {\tt G} associated to it. Parameters
for the model are assigned to equivalence classes of group elements, which are orbits
of some subgroup of the automorphism group of {\tt G}. This function returns the equivalence
classes as a list of list of group elements.
Example
M = K2Pmodel
B = buckets M
SeeAlso
Model
///
-------------------------------
-- LeafTree functionality
-------------------------------
-- leafTree
doc///
Key
leafTree
(leafTree,ZZ,List)
(leafTree,List,List)
(leafTree,Graph)
Headline
construct a LeafTree
Usage
leafTree(n,E)
leafTree(L,E)
leafTree(G)
Inputs
n:ZZ
the number of leaves
L:List
of leaves
E:List
of lists or sets specifying the internal edges
G:Graph
a tree
Outputs
:LeafTree
Description
Text
An object of class @TO LeafTree@ is specified by listing its leaves, and for each internal edge,
the partition the edge induces on the set of leaves.
{\tt L} is the set of leaves, or if an integer {\tt n} is input then the leaves will be be named $0,\ldots,n-1$.
{\tt E} is a list with one entry for each internal edge.
Each entry is a partition specified as a list or set of the leaves in one side of the partition.
Thus each edge can be specified in two possible ways.
An object of class @TO LeafTree@ can also be constructed from a @TO Graph@ provided the graph has no cycles.
Here we construct the quartet tree, which is the tree with 4 leaves and one internal edge.
Example
T = leafTree({a,b,c,d},{{a,b}})
leaves T
edges T
Text
Here is a tree with 5 leaves given as a @TO Graph@.
Example
G = graph{{a,b},{c,b},{b,d},{d,e},{d,f},{f,g},{f,h}}
T = leafTree G
leaves T
internalEdges T
///
-------------------------------
-- LeafTree
doc///
Key
LeafTree
(symbol ==,LeafTree,LeafTree)
Headline
a tree described in terms of its leaves
Description
Text
A tree can be described in terms of its leaves by specifying a leaf set
and specifying the edges as partitions of the leaf set.
This leaf centric description is particularly useful for phylogenetic trees.
The main constructor method is @TO leafTree@.
Example
T = leafTree({a,b,c,d},{{a,b}})
leaves T
edges T
G = graph{{a,e},{b,e},{e,f},{c,f},{d,f}}
leafTree G
SeeAlso
leafTree
///
-------------------------------
-- edges
doc///
Key
(edges,LeafTree)
Headline
list the edges of a tree
Usage
edges T
Inputs
T:LeafTree
Outputs
:List
the edges of {\tt T}
Description
Text
This function lists all edges of a tree. Each entry of the list is a @TO Set@ of the leaves on one side of the edge.
Example
T = leafTree(5,{{0,1}});
leaves T
edges T
SeeAlso
internalEdges
///
-------------------------------
-- internalEdges
doc///
Key
internalEdges
(internalEdges,LeafTree)
Headline
list the internal edges of a tree
Usage
internalEdges T
Inputs
T:LeafTree
Outputs
:List
the internal edges of {\tt T}
Description
Text
An internal edge of a tree is an edge that is not incident to a leaf.
This function lists such edges. Each entry of the list is @ofClass Set@ of the leaves on one side of the edge.
Example
G = graph {{0,4},{1,4},{4,5},{5,2},{5,3}};
T = leafTree G;
internalEdges T
SeeAlso
(edges,LeafTree)
///
-------------------------------
-- vertices
doc///
Key
(vertices,LeafTree)
Headline
list the vertices of a tree
Usage
vertices T
Inputs
T:LeafTree
Outputs
:List
the vertices of {\tt T}
Description
Text
This function lists all vertices of a tree. Each vertex is specified by the partition of the set of leaves
formed by removing the vertex. Each partition is given as a list of sets.
Example
T = leafTree(4,{{0,1}});
vertices T
#(vertices T)
Caveat
The leaves of {\tt T} in the output of {\tt vertices} have a different representation from the one in the output of @TO (leaves,LeafTree)@.
SeeAlso
internalVertices
(leaves,LeafTree)
///
-------------------------------
-- internalVertices
doc///
Key
internalVertices
(internalVertices,LeafTree)
Headline
list the internal vertices of a tree
Usage
internalVertices T
Inputs
T:LeafTree
Outputs
:List
the internal vertices of {\tt T}
Description
Text
An internal vertex of a tree is a vertex that is not a leaf, meaning it has degree at least 2.
This function lists such vertices. Each vertex is specified by the partition of the set of leaves
formed by removing the vertex. Each partition is given as a list of sets.
Example
T = leafTree(4,{{0,1}});
internalVertices T
#(internalVertices T)
SeeAlso
(vertices,LeafTree)
///
-------------------------------
-- leaves
doc///
Key
(leaves,LeafTree)
Headline
list the leaves of a tree
Usage
leaves T
Inputs
T:LeafTree
Outputs
:Set
the leaves of {\tt T}
Description
Text
This function outputs the leaves of the tree as an object of class @TO Set@.
Example
T = leafTree(4,{{0,1}});
vertices T
#(vertices T)
Caveat
The leaves have a different representation from the one in the output of @TO (vertices,LeafTree)@.
SeeAlso
(vertices,LeafTree)
///
-------------------------------
-- isIsomorphic
doc///
Key
isIsomorphic
(isIsomorphic,LeafTree,LeafTree)
Headline
check isomorphism of two tree
Usage
isIsomorphic(T,U)
Inputs
T:LeafTree
U:LeafTree
Outputs
:Boolean
if U and T are isomorphic
Description
Text
This function checks if two objects of class @TO LeafTree@ are isomorphic to each other as unlabeled graphs.
This is in contrast to equality of two objects of class @TO LeafTree@, which also checks whether they have the same leaf labeling.
Example
T = leafTree(4,{{0,1}});
U = leafTree(4,{{1,2}});
isIsomorphic(T,U)
///
-------------------------------
-- edgeCut
doc///
Key
edgeCut
(edgeCut,LeafTree,List,Thing)
(edgeCut,LeafTree,Set,Thing)
Headline
break up a tree at an edge
Usage
edgeCut(T,e,newl)
edgeCut(T,E,newl)
Inputs
T:LeafTree
e:Set
an edge specified by the set of leaves on one side of it
E:List
an edge specified by a list of the leaves on one side of it
newl:Thing
the label for a new leaf
Outputs
:List
of two @TO LeafTree@s that are subtrees of {\tt T}
Description
Text
This funtion outputs the two subtrees of {\tt T} obtained by deleting edge {\tt e} from {\tt T} and then re-adding the edge
to each of the two resulting subtrees. Both subtrees share a copy of the edge {\tt e}
and the newly labeled leaf adjacent to {\tt e}. Other than this overlap, they are disjoint.
Each subtree in {\tt P} may have at most one leaf that was not a leaf of {\tt T}, and therefore previously unlabeled.
The label for this new leaf is input as {\tt newl}.
Example
T = leafTree(4,{{0,1}})
P = edgeCut(T, set {0,1}, 4);
P#0
P#1
SeeAlso
vertexCut
///
-------------------------------
-- vertexCut
doc///
Key
vertexCut
(vertexCut,LeafTree,List,Thing,Thing)
(vertexCut,LeafTree,Set,Thing,Thing)
Headline
break up a tree at a vertex
Usage
vertexCut(T,e,l,newl)
vertexCut(T,E,l,newl)
Inputs
T:LeafTree
e:Set
an edge specified by the set of leaves on one side of it
E:List
an edge specified by a list of the leaves on one side of it
l:Thing
a leaf of the tree
newl:Thing
the label for a new leaf
Outputs
:List
of @TO LeafTree@s that are subtrees of {\tt T}
Description
Text
Vertices of a tree of class @TO LeafTree@ do not have explicit names. Therefore a vertex {\tt v} is specified by naming an edge {\tt e}
incident to {\tt v}, and leaf {\tt l} on the opposite side of the edge as {\tt v}.
The function outputs the subtrees of {\tt T} obtained by deleting the vertex {\tt v} from {\tt T}
and then re-adding {\tt v} to each of the resulting subtrees as a new leaf.
The new leaf on each subtree is adjacent to the edge previously adjacent
to {\tt v} on {\tt T}. Each subtree has a copy of the vertex labeled {\tt newl}, but their edge sets form a partition
of the edge set of {\tt T}.
Each subtree in {\tt P} has one leaf that was not a leaf of {\tt T}, and therefore previously unlabeled.
The label for this new leaf is input as {\tt newl}.
Example
T = leafTree(4,{{0,1}})
P = vertexCut(T, set {0,1}, 0, 4);
P#0
P#1
P#2
SeeAlso
edgeCut
///
-------------------------------
-- edgeContract
doc///
Key
edgeContract
(edgeContract,LeafTree,List)
(edgeContract,LeafTree,Set)
Headline
contract an edge of a tree
Usage
edgeContract(T,e)
edgeContract(T,E)
Inputs
T:LeafTree
e:Set
an edge specified by the set of leaves on one side of it
E:List
an edge specified by a list of the leaves on one side of it
Outputs
:LeafTree
obtained from {\tt T} by contracting the specified edge
Description
Text
This function produces a new object of class @TO LeafTree@ obtained by contracting the edge {\tt e} of tree {\tt T}.
Example
T = leafTree(4,{{0,1}})
edgeContract(T, set {0,1})
///
-------------------------------
-- labeledTrees
doc///
Key
labeledTrees
(labeledTrees,ZZ)
Headline
enumerate all labeled trees
Usage
labeledTrees n
Inputs
n:ZZ
the number of leaves
Outputs
:List
of all trees with {\tt n} leaves
Description
Text
This function enumerates all possible homeomorphically-reduced trees
(no degree-2 vertices) with {\tt n} leaves labeled by $0,\ldots, n-1$,
including all possible labelings. The trees are represented as objects of class @TO LeafTree@.
Example
L = labeledTrees 4
SeeAlso
labeledBinaryTrees
rootedTrees
rootedBinaryTrees
unlabeledTrees
///
-------------------------------
-- labeledBinaryTrees
doc///
Key
labeledBinaryTrees
(labeledBinaryTrees,ZZ)
Headline
enumerate all binary labeled trees
Usage
labeledTrees n
Inputs
n:ZZ
the number of leaves
Outputs
:List
of all binary trees with {\tt n} leaves
Description
Text
This function enumerates all possible binary trees with {\tt n} leaves
labeled by $0,\ldots, n-1$, including all possible labelings.
The trees are represented as an object of class @TO LeafTree@.
Example
L = labeledBinaryTrees 4
SeeAlso
labeledTrees
rootedTrees
rootedBinaryTrees
unlabeledTrees
///
-------------------------------
-- rootedTrees
doc///
Key
rootedTrees
(rootedTrees,ZZ)
Headline
enumerate all rooted trees
Usage
rootedTrees n
Inputs
n:ZZ
the number of leaves
Outputs
:List
of all rooted trees with $n$ leaves
Description
Text
This function enumerates all possible homeomorphically-reduced trees
(no degree-2 vertices) with a distinguished root and {\tt n-1} unlabeled leaves.
Each tree is an object of class @TO LeafTree@. For the purposes of representation,
the root is named {\tt 0} and the unlabeled leaves are named $1,\ldots,n-1$.
In other words each class of unlabeled rooted tree is represented once by a particular
labeling of that tree.
Example
L = rootedTrees 4
SeeAlso
labeledTrees
labeledBinaryTrees
rootedBinaryTrees
unlabeledTrees
///
-------------------------------
-- rootedBinaryTrees
doc///
Key
rootedBinaryTrees
(rootedBinaryTrees,ZZ)
Headline
enumerate all rooted binary trees
Usage
rootedBinaryTrees n
Inputs
n:ZZ
the number of leaves
Outputs
:List
of all rooted binary @TO LeafTree@s with {\tt n} leaves
Description
Text
This function enumerates all possible binary trees with a distinguished root
and {\tt n-1} unlabeled leaves. Each tree is an object of class @TO LeafTree@.
For the purposes of representation, the root is named $0$ and the unlabeled leaves
are named $1,\ldots,n-1$. In other words each class of unlabeled rooted tree is
represented once by a particular labeling of that tree.
Example
L = rootedBinaryTrees 5
SeeAlso
labeledTrees
labeledBinaryTrees
rootedTrees
unlabeledTrees
///
-------------------------------
-- unlabeledTrees
doc///
Key
unlabeledTrees
(unlabeledTrees,ZZ)
Headline
enumerate all unlabeled trees
Usage
unlabeledTrees n
Inputs
n:ZZ
the number of leaves
Outputs
:List
of all binary unlabeled trees with {\tt n} leaves
Description
Text
This function enumerates all possible binary trees with {\tt n} unlabeled leaves.
Each tree is an object of class @TO LeafTree@.
Each class of unlabeled tree is represented by a particular labeling of that tree.
Some duplicates may appear in the list, but each equivalence class is guaranteed to
appear at least once.
Example
L = unlabeledTrees 5
Caveat
For {\tt n} larger than 5, some equivalence classes of trees may appear more than once.
SeeAlso
labeledTrees
labeledBinaryTrees
rootedTrees
rootedBinaryTrees
///
-------------------------------
-- graph
doc///
Key
(graph,LeafTree)
Headline
convert a LeafTree to Graph
Usage
graph T
Inputs
T:LeafTree
Outputs
:Graph
Description
Text
This converts a @TO LeafTree@ representation of a tree into a @TO Graph@.
The internal vertices of a LeafTree are not named, so each vertex is specified by the partition of the set of leaves
formed by removing the vertex. Each partition is given as a @TO List@ of @TO Set@s.
Example
T = leafTree(4,{{0,1}})
G = graph T
adjacencyMatrix G
///
-------------------------------
-- digraph
doc///
Key
(digraph,LeafTree,List)
(digraph,LeafTree,Set)
Headline
convert a LeafTree to a Digraph
Usage
digraph(T,r)
Inputs
T:LeafTree
r:List
representing a vertex
Outputs
:Digraph
Description
Text
A rooted tree can be represented by an object of class @TO LeafTree@ and a choice of vertex to be the root.
This function converts such a representation of a rooted tree into an object of class @TO Digraph@ with
edges oriented away from the root.
The internal vertices of an object of class @TO LeafTree@ are not named, so each vertex is specified by the partition of the set of leaves
formed by removing the vertex. Each partition is given as a list of sets. This is also how the root vertex should
be passed to the function.
Example
T = leafTree(4,{{0,1}})
r = {set{0,1}, set{2}, set{3}}
D = digraph(T,r)
adjacencyMatrix D
///
-----------------------------------------------------------
----- TESTS -----
-----------------------------------------------------------
--Here is a test for the leafTree function. The tests depend on the
--internalEdges and leaves methods, to be tested elsewhere.
--We test a six leaf tree by using the internalEdges and leaves methods.
--The internal edges and leaves fully determine a tree.
--First we test the (L,E) input, then we check that
--using the second (ZZ,E) and (Graph) input gives the same result.
TEST ///
S = leafTree(6, {{0,1}, {0,1,2},{0,1, 2, 3}})
D = set internalEdges(S)
L = leaves(S)
d = set {set {0, 1, 2, 3}, set {0, 1, 2}, set {0,1 }}
l = set {0, 1, 2, 3, 4, 5}
assert( D == d)
assert( L == l)
--This test is very simple and potential problems would mostly come from dependencies on graph package
--if this test is not acceptable, go back to leafTree and internalEdges and think of
--how they depend on Graphs package (unrefereed, no low level functionality tests)
///
--The following is a test for leafColorings. We check that leafColorings gives
--the correct SETS. We check that leafColorings(4,CFNmodel) gives the correct
--output set. We also check that leafColorings gives the same output
--for the JCmodel, the K2Pmodel, and the K3Pmodel on the tree with 4 leaves,
--as this method should only depend on the group, and not the acutal model.
TEST ///
A =set leafColorings(4, CFNmodel)
B =set leafColorings(4, JCmodel)
C =set leafColorings(4, K2Pmodel)
D =set leafColorings(4, K3Pmodel)
L =set {(0_(ZZ/2), 0_(ZZ/2), 0_(ZZ/2), 0_(ZZ/2)),
(0_(ZZ/2), 0_(ZZ/2), 1_(ZZ/2), 1_(ZZ/2)),
(0_(ZZ/2), 1_(ZZ/2), 0_(ZZ/2), 1_(ZZ/2)),
(0_(ZZ/2), 1_(ZZ/2), 1_(ZZ/2), 0_(ZZ/2)),
(1_(ZZ/2), 0_(ZZ/2), 0_(ZZ/2), 1_(ZZ/2)),
(1_(ZZ/2), 0_(ZZ/2), 1_(ZZ/2), 0_(ZZ/2)),
(1_(ZZ/2), 1_(ZZ/2), 0_(ZZ/2), 0_(ZZ/2)),
(1_(ZZ/2), 1_(ZZ/2), 1_(ZZ/2), 1_(ZZ/2))}
assert(A == L)
assert(B == C)
assert(C == D)
///
--Here we give a small test that qRing produces a polynomial ring
--in the correct number of variables and that the elements are as expected
TEST ///
S = leafTree(4, {{0,1}})
R = qRing(S, JCmodel)
P = qRing(4, JCmodel)
assert(dim R == 64)
assert((vars R)_(0,0) == q_(0,0,0,0))
assert(dim P == 64)
assert((vars P)_(0,0) == q_(0,0,0,0))
///
--The following gives tests for phyloToricFP.
--We test the 4-claw, which in
--Sturmfels/Sullivant is the 3-claw. We manually construct the ideal of invariants
--as the kernel of the ring homomorphism determined by the parameterization.
--Similarly, one can check that this is the same as the ideal in Sturmfels/Sullivant
--example 3. To do so, you must include the fourth index on their parameters
--to be the sum of the first three.
TEST ///
T = QQ[q_(0,0,0,0),q_(0,0,1,1), q_(0,1,0,1), q_(0,1,1,0),
q_(1,0,0,1), q_(1,0,1,0), q_(1,1,0,0), q_(1,1,1,1)]
J = phyloToricFP(4, {}, CFNmodel,QRing=>T)
R = QQ[a_0, a_1, b_0, b_1, c_0, c_1, d_0, d_1]
f = map(R, T, {a_0*b_0*c_0*d_0, a_0*b_0*c_1*d_1, a_0*b_1*c_0*d_1,
a_0*b_1*c_1*d_0, a_1*b_0*c_0*d_1, a_1*b_0*c_1*d_0,
a_1*b_1*c_0*d_0, a_1*b_1*c_1*d_1})
I = kernel f
assert(I == J)
///
--The following gives tests for phyloToricLinears and phyloToric42.
--We also include another test for phyloToricFP.
--We construct the toric ideal in the quartet tree with single
--non-trivial split using the Jukes-Cantor model. We construct the ideal as the
--kernel of the homomorphism defined by the standard parameterization. We check that
--this ideal matches the ideal defined by phyloToricFP, I == J. We extract the
--linear generators for the kernel I, and check that these generate the same
--ideal as the generatros given as output for phyloToricLinears, M == Q.
--We also check the number of linear generators defined by N and those defined
--by phyloToricLinears, P. Although these sets are not minimal, we check that
--each list is 51. This coincides with the fact that there are 13 distinct
--Fourier coordinates for this tree with the JCmodel, and there are 64 total
--parameters.
TEST ///
T = QQ[q_(0,0,0,0),
q_(0,0,1,1),q_(0,0,2,2),q_(0,0,3,3),
q_(0,1,0,1),q_(0,1,1,0),q_(0,1,2,3),
q_(0,1,3,2),q_(0,2,0,2),q_(0,2,1,3),
q_(0,2,2,0),q_(0,2,3,1),q_(0,3,0,3),
q_(0,3,1,2),q_(0,3,2,1),q_(0,3,3,0),
q_(1,0,0,1),q_(1,0,1,0),q_(1,0,2,3),
q_(1,0,3,2),q_(1,1,0,0),q_(1,1,1,1),
q_(1,1,2,2),q_(1,1,3,3),q_(1,2,0,3),
q_(1,2,1,2),q_(1,2,2,1),q_(1,2,3,0),
q_(1,3,0,2),q_(1,3,1,3),q_(1,3,2,0),
q_(1,3,3,1),q_(2,0,0,2),q_(2,0,1,3),
q_(2,0,2,0),q_(2,0,3,1),q_(2,1,0,3),
q_(2,1,1,2),q_(2,1,2,1),q_(2,1,3,0),
q_(2,2,0,0),q_(2,2,1,1),q_(2,2,2,2),
q_(2,2,3,3),q_(2,3,0,1),q_(2,3,1,0),
q_(2,3,2,3),q_(2,3,3,2),q_(3,0,0,3),
q_(3,0,1,2),q_(3,0,2,1),q_(3,0,3,0),
q_(3,1,0,2),q_(3,1,1,3),q_(3,1,2,0),
q_(3,1,3,1),q_(3,2,0,1),q_(3,2,1,0),
q_(3,2,2,3),q_(3,2,3,2),q_(3,3,0,0),
q_(3,3,1,1),q_(3,3,2,2),q_(3,3,3,3)]
R = QQ[a0, a1, b0, b1, c0, c1, d0, d1, e0, e1]
f = map(R,T, {a0*b0*c0*d0*e0,
a0*b0*c1*d1*e0, a0*b0*c1*d1*e0, a0*b0*c1*d1*e0,
a0*b1*c0*d1*e1, a0*b1*c1*d0*e1, a0*b1*c1*d1*e1,
a0*b1*c1*d1*e1, a0*b1*c0*d1*e1, a0*b1*c1*d1*e1,
a0*b1*c1*d0*e1, a0*b1*c1*d1*e1, a0*b1*c0*d1*e1,
a0*b1*c1*d1*e1, a0*b1*c1*d1*e1, a0*b1*c1*d0*e1,
a1*b0*c0*d1*e1, a1*b0*c1*d0*e1, a1*b0*c1*d1*e1,
a1*b0*c1*d1*e1, a1*b1*c0*d0*e0, a1*b1*c1*d1*e0,
a1*b1*c1*d1*e0, a1*b1*c1*d1*e0, a1*b1*c0*d1*e1,
a1*b1*c1*d1*e1, a1*b1*c1*d1*e1, a1*b1*c1*d0*e1,
a1*b1*c0*d1*e1, a1*b1*c1*d1*e1, a1*b1*c1*d0*e1,
a1*b1*c1*d1*e1, a1*b0*c0*d1*e1, a1*b0*c1*d1*e1,
a1*b0*c1*d0*e1, a1*b0*c1*d1*e1, a1*b1*c0*d1*e1,
a1*b1*c1*d1*e1, a1*b1*c1*d1*e1, a1*b1*c1*d0*e1,
a1*b1*c0*d0*e0, a1*b1*c1*d1*e0, a1*b1*c1*d1*e0,
a1*b1*c1*d1*e0, a1*b1*c0*d1*e1, a1*b1*c1*d0*e1,
a1*b1*c1*d1*e1, a1*b1*c1*d1*e1, a1*b0*c0*d1*e1,
a1*b0*c1*d1*e1, a1*b0*c1*d1*e1, a1*b0*c1*d0*e1,
a1*b1*c0*d1*e1, a1*b1*c1*d1*e1, a1*b1*c1*d0*e1,
a1*b1*c1*d1*e1, a1*b1*c0*d1*e1, a1*b1*c1*d0*e1,
a1*b1*c1*d1*e1, a1*b1*c1*d1*e1, a1*b1*c0*d0*e0,
a1*b1*c1*d1*e0, a1*b1*c1*d1*e0, a1*b1*c1*d1*e0})
I = ker(f)
--Here's the test for phyloToric42
S = leafTree(4, {{0,1}})
L = phyloToric42(S, JCmodel, QRing=>T)
assert(I ==L)
--Here's the test for phyloToricFP
J = phyloToricFP(4, {{0,1}}, JCmodel, QRing=>T)
assert(I == J)
--Here's the test for phyloToricLinears
K={1}
N = for i to 83 when degree(I_i) == K list I_i
M = ideal N
#N
P = phyloToricLinears(4, {{0,1}}, JCmodel, QRing=>T)
Q = ideal P
#P
assert(#N === 51)
assert(#P === 51)
assert(M == Q)
///
--Here's a test for phyloToricAMatrix. We only test the set of columns, since
--this only tests the parameterization up to permutation of coordinates.
TEST ///
A = phyloToricAMatrix(4, {}, CFNmodel)
B = matrix{{1,0,1,0,1,0,1,0},
{1,0,1,0,0,1,0,1},
{1,0,0,1,1,0,0,1},
{1,0,0,1,0,1,1,0},
{0,1,1,0,1,0,0,1},
{0,1,1,0,0,1,1,0},
{0,1,0,1,1,0,1,0},
{0,1,0,1,0,1,0,1}}
C = transpose B
D = set { submatrix(C, {0}), submatrix(C,{1}), submatrix(C, {2}),
submatrix(C, {3}), submatrix(C, {4}), submatrix(C, {5}),
submatrix(C, {6}), submatrix(A, {7}) }
E = set { submatrix(A, {0}), submatrix(A,{1}), submatrix(A, {2}),
submatrix(A, {3}), submatrix(A, {4}), submatrix(A, {5}),
submatrix(A, {6}), submatrix(A, {7})}
assert(D == E)
///
--Here is a test for internalVertices. We test a binary tree and also
--a claw tree. We convert all lists to sets to allow for different ordering.
TEST ///
S = leafTree(6, {{0,1}, {0,1,2},{0,1, 2, 3}})
T = leafTree(6, {{}})
internalEdges(S)
internalEdges(T)
IVS= internalVertices(S)
IVSs = set IVS
A = set {set {set {0, 1, 2, 3}, set {4}, set {5}}, set {set {0, 1, 2}, set {3},
set {4, 5}}, set {set {0, 1}, set {2}, set {3, 4, 5}}, set {set {0},
set {1}, set {2, 3, 4, 5}}}
assert( IVSs == A)
IVT = internalVertices(T)
IVTs = set IVT
B = set{set{ set{0}, set {1}, set {2}, set {3}, set {4}, set {5}}}
assert( IVTs == B)
///
--Here is a test for internalEdges. We test a quartet and a 4-claw.
TEST ///
S = leafTree(4, {{0,1}})
A = set internalEdges(S)
B = set{set{0,1}, set{2,3}}
assert( #(A * B) == 1)
T = leafTree(4, {})
C = set internalEdges(T)
assert( C == set{})
///
--Note for the vertexCut test: This test was written before user declared label
--of new leaf. Now that this label is declared, a much simpler and perhaps more
--rigourous test can be written. We will leave the old test for now.
TEST ///
--Here is a test for vertexCut.
S = leafTree(6, {{0,1}, {0,1,2},{0,1, 2, 3}})
VC=vertexCut(S,{0,1,2}, 0, 6)
--First we test that this vertex-cut gives us three connected components, of
--2,3, and 4 leaves each. We check that each tree has the appropriately labeled
--leaves. We have an additional test in which the added vertex is labeled the same on
--each connected component.
M = set{0,1,2, 3, 4, 5}
N= set{0,1,2}
NC = set {3, 4, 5}
L1= for x in VC list leaves(x)
N4 = for x in L1 list(if #x <4 then continue; x)
O4 = for x in N4 list #x
P4 =set flatten for x in N4 list elements x
assert(# set N4 == 1)
assert(O4 == {4})
assert(P4 * NC == set {})
N3 = for x in L1 list(if #x <3 or #x > 3 then continue; x)
P3 =set flatten for x in N3 list elements x
assert(# set N3 == 1)
assert(P3 * N == set{})
N2 = for x in L1 list(if #x <2 or #x > 2 then continue; x)
P2 =set flatten for x in N2 list elements x
assert(# set N2 == 1)
assert(P2 * N == set{})
N1 = for x in L1 list(if #x > 1 then continue; x)
P1 =set flatten for x in N1 list elements x
assert(P1 == set{})
L = (P4 * N ) + (P3 * NC) + (P2 * NC)
assert( L == M)
assert((P4 - N) == (P3 - NC))
assert((P3 - NC) == (P2 - NC))
--Second, we test that the unique 4 leaf component is a quartet tree (and not
--a claw).
A=flatten for x in VC list edges(x)
C =set for x in A list(if #x == 1 then continue; x)
assert(C =!= set{})
///
--Here's a test for edgeCut.
--First we test that this edge-cut gives us two connected components, of
--4 leaves each. We check that each tree has the appropriately labeled
--leaves.
--Second, we test that the 4 leaf components are the correct quartets.
TEST///
S = leafTree(6, {{0,1}, {0,1,2},{0,1, 2, 3}})
EC = edgeCut(S, {0,1,2}, 6)
N= set{0,1,2,6}
NC = set {3, 4, 5,6}
L1= for x in EC list leaves(x)
T1 = L1#0
T2 = L1#1
assert( T1 == N or T1 == NC)
assert( T2 == N or T2 == NC)
A=flatten for x in EC list edges(x)
C =set for x in A list(if #x == 1 then continue; x)
c1 = set{0,1}
d1 = set{2,6}
c2 = set{4,5}
d2 = set{3,6}
assert(C#?c1 or C#?d1)
assert(C#?c2 or C#?d2)
///
--Here is a simple test for edgeContract.
TEST ///
S = leafTree(4, {{0,1}})
T = edgeContract(S, set{0,1})
A = leaves(T)
B = set internalEdges(T)
assert(A == set {0,1,2,3})
assert(B == set {})
///
TEST ///
-- We test the function joinIdeal by computing the join of the ideal of
-- the Veronese map with n=1 and d=7 and the ideal of the Segre
-- embedding of P1 x P3. The ideal is computed directly from the
-- parameterization and using the function joinIdeal.
R = QQ[x1,x2,x3,x4,x5,x6,x7,x8]
S = QQ[a0,a1,a2,a3,b0,b1,b2,s,t]
f1 = map(S,R,{a0*b0, a1*b0, a2*b0, a3*b0,
a0*b1, a1*b1, a2*b1, a3*b1});
f2 = map(S,R,{s^7*t^0, s^6*t^1 ,s^5*t^2, s^4*t^3,s^3*t^4,s^2*t^5,s^1*t^6,s^0*t^7 });
g = map(S,R,{a0*b0 + s^7*t^0, a1*b0 + s^6*t^1, a2*b0 + s^5*t^2, a3*b0 + s^4*t^3,
a0*b1 + s^3*t^4, a1*b1 + s^2*t^5, a2*b1 + s^1*t^6, a3*b1 + s^0*t^7});
I = ker(f1);
J = ker(f2);
assert(joinIdeal(I,J) == ker(g))
///
TEST ///
-- We test the function toricSecantDim by computing the
-- dimension of a second secant of the CFN model
-- which is known to be non-defective.
-- We also verify that the dimenson of the secant for the
-- CFN model for a 4-leaf tree is no larger than the ambient dimension.
A = phyloToricAMatrix(6, {{0,1},{2,3},{4,5}},CFNmodel);
assert(toricSecantDim(A,1) == dim(phyloToric42(6, {{0,1},{2,3},{4,5}},CFNmodel)))
assert(toricSecantDim(A,2) == 20)
assert(toricSecantDim(phyloToricAMatrix(4, {{0,1}},CFNmodel),2) == 8)
///
TEST ///
-- We test the function toricJoinDim using
-- joins of 2 and 3 6-leaf trees.
-- It is known that for the JCmodel, joins
-- of 2 or 3 arbitrary trees with 6 or more leaves are
-- non-defective.
A = phyloToricAMatrix(6, {{0,1},{2,3},{4,5}},JCmodel);
B = phyloToricAMatrix(6, {{0,1},{0,1,2},{4,5}},JCmodel);
C = phyloToricAMatrix(6, {{1,2},{3,4},{0,5}},JCmodel);
assert(toricSecantDim(A,1) == 10)
assert(toricSecantDim(B,1) == 10)
assert(toricSecantDim(C,1) == 10)
assert(toricJoinDim(A,B) == 20)
assert(toricJoinDim({A,B,C}) == 30)
///
TEST ///
-- The function phyloToricQuads is tested
-- by verifying that the ideal generated by
-- the polynomials returned
-- modulo the linear invariants
-- is equal to the degree 2
-- generators of the toric ideal.
Tree = {{0,1},{0,1,2}};
n = 5;
M = JCmodel;
S = qRing(n,M);
L = phyloToricLinears(n,Tree,M,QRing=>S);
T = S/L;
I = ideal phyloToricQuads(n,Tree,M,QRing=>T);
J = phyloToric42(n,Tree,M,QRing=>T);
K = ideal(for i in flatten entries mingens J list (if (degree i)#0 > 2 then continue; i));
assert(I == K)
///
TEST ///
-- The function phyloToricRandom is tested by
-- verifying that it produces a polynomial in the
-- ideal of invariants for the appropriate model.
Tree = {{0,1}};
n = 4;
M = K2Pmodel;
S = qRing(n,M)
I = phyloToric42(n,Tree, M, QRing=> S);
f = phyloToricRandom(n,Tree,M, QRing=>S);
assert(f % I == 0)
///
TEST ///
Tree = {{0,1}};
n = 4;
M = K2Pmodel;
S = qRing(n,M)
I = phyloToric42(n,Tree, M, QRing=> S);
f = phyloToricLinears(n,Tree,M, QRing=>S,Random=>true);
g = phyloToricQuads(n,Tree,M, QRing=>S,Random=>true);
assert(f_0 % I == 0)
assert(g_0 % I == 0)
///
TEST ///
--The function fourierToProbability is tested by computing
--the ideal for the same tree in two different ways. The first is
--by computing directly from the parameterization in probability
--coordinates and the second is by using phyloToric42 to compute
--the ideal in Fourier coordinates and then forming an ideal by
--converting each of the generators into probablity coordinates.
--We assert the two ideals are equal modulo the certain linear invariants
--that are supressed when computing in the ring of Fourier coordinates.
S = pRing(4,CFNmodel);
L = ideal apply(8,i->(S_i - S_(15-i)))
R1 = S/L;
R2 = QQ[a0,a1,b0,b1,c0,c1,d0,d1,e0,e1];
f = map(R2,R1,
{a0*b0*c0*d0*e0+a0*b1*c1*d0*e1+a1*b0*c1*d1*e0+a1*b1*c0*d1*e1,
a0*b0*c1*d0*e1+a0*b1*c0*d0*e0+a1*b0*c0*d1*e1+a1*b1*c1*d1*e0,
a0*b0*c0*d0*e1+a0*b1*c1*d0*e0+a1*b0*c1*d1*e1+a1*b1*c0*d1*e0,
a0*b0*c1*d0*e0+a0*b1*c0*d0*e1+a1*b0*c0*d1*e0+a1*b1*c1*d1*e1,
a0*b0*c0*d1*e0+a0*b1*c1*d1*e1+a1*b0*c1*d0*e0+a1*b1*c0*d0*e1,
a0*b0*c1*d1*e1+a0*b1*c0*d1*e0+a1*b0*c0*d0*e1+a1*b1*c1*d0*e0,
a0*b0*c0*d1*e1+a0*b1*c1*d1*e0+a1*b0*c1*d0*e1+a1*b1*c0*d0*e0,
a0*b0*c1*d1*e0+a0*b1*c0*d1*e1+a1*b0*c0*d0*e0+a1*b1*c1*d0*e1,
a0*b0*c1*d1*e0+a0*b1*c0*d1*e1+a1*b0*c0*d0*e0+a1*b1*c1*d0*e1,
a0*b0*c0*d1*e1+a0*b1*c1*d1*e0+a1*b0*c1*d0*e1+a1*b1*c0*d0*e0,
a0*b0*c1*d1*e1+a0*b1*c0*d1*e0+a1*b0*c0*d0*e1+a1*b1*c1*d0*e0,
a0*b0*c0*d1*e0+a0*b1*c1*d1*e1+a1*b0*c1*d0*e0+a1*b1*c0*d0*e1,
a0*b0*c1*d0*e0+a0*b1*c0*d0*e1+a1*b0*c0*d1*e0+a1*b1*c1*d1*e1,
a0*b0*c0*d0*e1+a0*b1*c1*d0*e0+a1*b0*c1*d1*e1+a1*b1*c0*d1*e0,
a0*b0*c1*d0*e1+a0*b1*c0*d0*e0+a1*b0*c0*d1*e1+a1*b1*c1*d1*e0,
a0*b0*c0*d0*e0+a0*b1*c1*d0*e1+a1*b0*c1*d1*e0+a1*b1*c0*d1*e1})
I = ker(f)
T = leafTree(4,{{0,1}});
M = CFNmodel;
J = phyloToric42(T,M);
MJ = mingens J;
FToP = fourierToProbability(R1,ring J,4,M)
J = ideal(FToP MJ_(0,0), FToP MJ_(0,1))
assert(J == I)
///
end
------------------------------------------------------------
restart
installPackage "PhylogeneticTrees"
viewHelp PhylogeneticTrees