Skip to content

Commit 07f2396

Browse files
ollieclarke8787d-torrance
ollieclarke8787
authored and
d-torrance
committed
Update from develop repo
Adding the changes from ollieclarke8787/WorkshopDurham2024@8f94f39 and ollieclarke8787/WorkshopDurham2024@a45e5e9
1 parent 85aba8f commit 07f2396

File tree

1 file changed

+121
-21
lines changed

1 file changed

+121
-21
lines changed

‎M2/Macaulay2/packages/AllMarkovBases.m2‎

Lines changed: 121 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ newPackage(
22
"AllMarkovBases",
33
Version => "1.0",
44
Date => "May 08, 2025",
5-
Headline => "computing all minimal Markov bases of a configuration matrix",
5+
Headline => "package for computing all minimal Markov bases of a configuration matrix",
66
Authors => {
77
{Name => "Alexander Milner",
88
Email => "[email protected]",
@@ -11,7 +11,6 @@ newPackage(
1111
Email => "[email protected]",
1212
HomePage => "https://www.oliverclarkemath.com/"}
1313
},
14-
Keywords => {"Algebraic Statistics"},
1514
AuxiliaryFiles => false,
1615
DebuggingMode => false,
1716
PackageExports => {"FourTiTwo","Graphs","Normaliz"}
@@ -70,8 +69,8 @@ computeFiberInternal = method(
7069
);
7170
computeFiberInternal (Matrix,List) := opts -> (A, val) -> (
7271
if not (A.cache#"fibers")#?val or (opts.ReturnConnectedComponents and (A.cache#"fiberStarters")#?val and not (A.cache#"fiberComponents")#?val) then(
73-
if opts.FiberAlgorithm == "lattice" then computeFiberInternalLattice(A,val,ReturnConnectedComponents=>opts.ReturnConnectedComponents)
74-
else if opts.FiberAlgorithm == "fast" then computeFiberInternalFast(A,val)
72+
if opts.FiberAlgorithm == "fast" then computeFiberInternalFast(A,val)
73+
else if opts.FiberAlgorithm == "lattice" then computeFiberInternalLattice(A,val,ReturnConnectedComponents=>opts.ReturnConnectedComponents)
7574
else if opts.FiberAlgorithm == "decompose" or opts.FiberAlgorithm == "markov" then computeFiberInternalDecompose(A,val,ReturnConnectedComponents=>opts.ReturnConnectedComponents,FiberAlgorithm=>opts.FiberAlgorithm)
7675
else error("unknown input for FiberAlgorithm option");
7776
);
@@ -89,7 +88,11 @@ computeFiberInternalFast (Matrix,List) := (A,val) -> (
8988
);
9089
u
9190
);
92-
validMoves := for move in entries A.cache#"MarkovBasis" list if ((A.cache#"fiberValues")#move << val) and ((A.cache#"fiberValues")#move != val) then move else continue;
91+
validMoves := if A.cache#"NN" then(
92+
for move in entries A.cache#"MarkovBasis" list if ((A.cache#"fiberValues")#move << val) and ((A.cache#"fiberValues")#move != val) then move else continue
93+
)else(
94+
for move in entries A.cache#"MarkovBasis" list if (A.cache#"fiberValues")#move != val then move else continue
95+
);
9396
validMoves = new MutableHashTable from ((v -> {v,true}) \ validMoves);
9497
out := for i from 0 to #buildFiber - 1 list(
9598
cc := set {buildFiber#i};
@@ -121,10 +124,11 @@ computeFiberInternalLattice = method(
121124
}
122125
);
123126
computeFiberInternalLattice (Matrix,List) := opts -> (A,val) -> (
124-
M := transpose matrix for basisElement in entries A.cache#"MarkovBasis" list if (A.cache#"fiberValues")#basisElement<<val then basisElement else continue;
127+
M := transpose A.cache#"MarkovBasis";
125128
if (numColumns M)==0 then M = transpose matrix{toList((numColumns A):0)};
126-
latticeOut := if (A.cache#"fiberStarters")#?val then set latticePointsFromMoves(M,transpose matrix{first toList (A.cache#"fiberStarters")#val})
127-
else(
129+
latticeOut := if (A.cache#"fiberStarters")#?val then (
130+
set latticePointsFromMoves(M,transpose matrix{first toList (A.cache#"fiberStarters")#val})
131+
)else(
128132
u:={};
129133
for row in entries toricMarkov ((matrix vector val) | A) do(
130134
r := drop(row,1);
@@ -161,19 +165,30 @@ computeCCs List := L -> (
161165
computeFiberInternalDecompose = method(
162166
Options => {
163167
ReturnConnectedComponents => false,
164-
FiberAlgorithm => "decompose"
168+
FiberAlgorithm => "fast"
165169
}
166170
);
167171
computeFiberInternalDecompose (Matrix,List) := opts -> (A,val) -> (
168-
if opts.FiberAlgorithm == "decompose" then(
169-
A.cache#"Ring" = ZZ(monoid[Variables => numRows A + numColumns A,MonomialOrder=>Eliminate numRows A]);
170-
A.cache#"RingGenerators" = gens A.cache#"Ring";
171-
A.cache#"toricIdeal"=toricGroebner(id_(ZZ^(numRows A)) | A, A.cache#"Ring");
172+
if opts.FiberAlgorithm == "decompose" and not A.cache#?"Ring" then(
173+
if A.cache#"NN" then(
174+
A.cache#"Ring" = ZZ(monoid[Variables => numRows A + numColumns A,MonomialOrder=>Eliminate numRows A]);
175+
A.cache#"RingGenerators" = gens A.cache#"Ring";
176+
A.cache#"toricIdeal"=toricGroebner(id_(ZZ^(numRows A)) | A, A.cache#"Ring");
177+
)else(
178+
A.cache#"Ring" = ZZ(monoid[Variables => 2*numRows A + numColumns A,MonomialOrder=>Eliminate (2*numRows A)]);
179+
A.cache#"RingGenerators" = gens A.cache#"Ring";
180+
A.cache#"toricIdeal"=toricGroebner(id_(ZZ^(numRows A)) | -id_(ZZ^(numRows A)) | A, A.cache#"Ring");
181+
);
172182
);
173183
if (A.cache#"fiberStarters")#?val and opts.ReturnConnectedComponents then (
174184
out := for z in pairs A.cache#"fiberValues" list(
175-
if z#1==val or not (z#1 << val) then continue;
176-
resid := val-z#1;
185+
if z#1 == val then continue;
186+
resid := val-z#1;
187+
if A.cache#"NN" then(
188+
if not (z#1 << val) then continue;
189+
)else(
190+
if any(A.cache#"supportingHyperplanes",z->(matrix{z}*(vector resid))_0<0) then continue;
191+
);
177192
if not (A.cache#"fibers")#?resid then fibRecursion(A,resid,FiberAlgorithm=>opts.FiberAlgorithm);
178193
if #((A.cache#"fibers")#resid)==0 then continue;
179194
fibAdd((A.cache#"posNeg")#(z#0),(A.cache#"fibers")#resid)
@@ -186,6 +201,7 @@ computeFiberInternalDecompose (Matrix,List) := opts -> (A,val) -> (
186201
)else fibRecursion(A,val,FiberAlgorithm=>opts.FiberAlgorithm);
187202
);
188203

204+
189205
-- recursive method using decompose algorithm (unexported)
190206
fibRecursion = method(
191207
Options => {
@@ -195,9 +211,13 @@ fibRecursion = method(
195211
);
196212
fibRecursion (Matrix,List) := opts -> (A, val) -> (
197213
out := union for z in pairs A.cache#"fiberValues" list(
198-
if not (z#1 << val) then continue;
199-
if z#1==val and (A.cache#"fiberStarters")#?val then continue (A.cache#"posNeg")#(z#0);
200-
resid := val-z#1;
214+
if z#1==val and (A.cache#"fiberStarters")#?val then continue (A.cache#"posNeg")#(z#0);
215+
resid := val-z#1;
216+
if A.cache#"NN" then(
217+
if not (z#1 << val) then continue;
218+
)else(
219+
if any(A.cache#"supportingHyperplanes",z->(matrix{z}*(vector resid))_0<0) then continue;
220+
);
201221
if not (A.cache#"fibers")#?resid then fibRecursion(A,resid,FiberAlgorithm=>opts.FiberAlgorithm);
202222
if #((A.cache#"fibers")#resid) == 0 then continue;
203223
fibAdd((A.cache#"posNeg")#(z#0),(A.cache#"fibers")#resid)
@@ -209,8 +229,13 @@ fibRecursion (Matrix,List) := opts -> (A, val) -> (
209229
if row#0 == 1 and all(r,z->z<=0) then (out = set{-r}; break;);
210230
);
211231
)else (
212-
e:=first exponents (product(for i from 0 to #val-1 list ((A.cache#"RingGenerators")#i)^(val#i)) % A.cache#"toricIdeal");
213-
if take(e,numRows A) == toList((numRows A):0) then out = set{take(e,-numColumns A)};
232+
if A.cache#"NN" then(
233+
e1:=first exponents (product(for i from 0 to #val-1 list ((A.cache#"RingGenerators")#i)^(val#i)) % A.cache#"toricIdeal");
234+
if take(e1,numRows A) == toList((numRows A):0) then out = set{take(e1,-numColumns A)};
235+
)else(
236+
e2:=first exponents (product(for i from 0 to #val-1 list if val#i>0 then ((A.cache#"RingGenerators")#i)^(val#i) else ((A.cache#"RingGenerators")#(i+numRows A))^(-val#i)) % A.cache#"toricIdeal");
237+
if take(e2,2*numRows A) == toList((2*numRows A):0) then out = set{take(e2,-numColumns A)};
238+
);
214239
);
215240
);
216241
(A.cache#"fibers")#val = out;
@@ -232,7 +257,7 @@ fiberGraph = method(
232257
}
233258
);
234259
fiberGraph Matrix := opts -> A -> (
235-
if not A.cache#?"MarkovBasis" then setupFibers A;
260+
if not A.cache#?"MarkovBasis" then elapsedTime setupFibers A;
236261
if opts.ReturnConnectedComponents then(
237262
if not A.cache#"componentsComputed" then (
238263
for val in rsort keys A.cache#"fiberStarters" do computeFiberInternal(A,val,ReturnConnectedComponents=>true,FiberAlgorithm=>opts.FiberAlgorithm);
@@ -273,6 +298,12 @@ setupFibers Matrix := A -> (
273298
A.cache#"fiberComponents" = new MutableHashTable;
274299
A.cache#"componentsComputed" = false;
275300
A.cache#"fiberGraphs" = new MutableHashTable;
301+
A.cache#"NN"= all(for row in entries A list all(for el in row list el>=0,z->z),z->z);
302+
N:=normaliz(A,"normalization",allComputations=>true);
303+
S:=entries N#"sup";
304+
G:=entries N#"gen";
305+
tem := transpose matrix for row in entries A list S#(position(G,z->z==row));
306+
A.cache#"supportingHyperplanes" = for el in entries tem list if el == toList(#S:0) then continue else el;
276307
);
277308

278309

@@ -1400,3 +1431,72 @@ installPackage "AllMarkovBases"
14001431
check AllMarkovBases
14011432

14021433
viewHelp AllMarkovBases
1434+
1435+
1436+
1437+
1438+
1439+
-- computeFiberInternalDecompose = method(
1440+
-- Options => {
1441+
-- ReturnConnectedComponents => false,
1442+
-- FiberAlgorithm => "fast"
1443+
-- }
1444+
-- );
1445+
-- computeFiberInternalDecompose (Matrix,List) := opts -> (A,val) -> (
1446+
-- if opts.FiberAlgorithm == "decompose" and not A.cache#?"Ring" then(
1447+
-- A.cache#"Ring" = ZZ(monoid[Variables => 2*numRows A + numColumns A,MonomialOrder=>Eliminate (2*numRows A)]);
1448+
-- A.cache#"RingGenerators" = gens A.cache#"Ring";
1449+
-- A.cache#"toricIdeal"=toricGroebner(id_(ZZ^(numRows A)) | -id_(ZZ^(numRows A)) | A, A.cache#"Ring");
1450+
-- );
1451+
-- if (A.cache#"fiberStarters")#?val and opts.ReturnConnectedComponents then (
1452+
-- out := for z in pairs A.cache#"fiberValues" list(
1453+
-- print z#1;
1454+
-- if z#1==val or not (z#1 << val) then continue;
1455+
-- resid := val-z#1;
1456+
-- if not (A.cache#"fibers")#?resid then fibRecursion(A,resid,FiberAlgorithm=>opts.FiberAlgorithm);
1457+
-- if #((A.cache#"fibers")#resid)==0 then continue;
1458+
-- fibAdd((A.cache#"posNeg")#(z#0),(A.cache#"fibers")#resid)
1459+
-- );
1460+
-- G := new HashTable from for i from 0 to #out-1 list (out#i,for j from i+1 to #out-1 list if not #intersect(out#i,out#j)==0 then out#j else continue);
1461+
-- out = apply(connectedComponents graph(out,pairs G),z->toList union z);
1462+
-- out = out | apply(toList ((A.cache#"fiberStarters")#val - (flatten out)),v->{v});
1463+
-- (A.cache#"fiberComponents")#val = out;
1464+
-- (A.cache#"fibers")#val = flatten out;
1465+
-- )else fibRecursion(A,val,FiberAlgorithm=>opts.FiberAlgorithm);
1466+
-- );
1467+
1468+
-- -- recursive method using decompose algorithm (unexported)
1469+
-- fibRecursion = method(
1470+
-- Options => {
1471+
-- FiberAlgorithm => "decompose"
1472+
-- }
1473+
1474+
-- );
1475+
-- fibRecursion (Matrix,List) := opts -> (A, val) -> (
1476+
-- print val;
1477+
-- out := union for z in pairs A.cache#"fiberValues" list(
1478+
-- if not (z#1 << val) then continue;
1479+
-- if z#1==val and (A.cache#"fiberStarters")#?val then continue (A.cache#"posNeg")#(z#0);
1480+
-- resid := val-z#1;
1481+
-- if not (A.cache#"fibers")#?resid then fibRecursion(A,resid,FiberAlgorithm=>opts.FiberAlgorithm);
1482+
-- if #((A.cache#"fibers")#resid) == 0 then continue;
1483+
-- fibAdd((A.cache#"posNeg")#(z#0),(A.cache#"fibers")#resid)
1484+
-- );
1485+
-- if #out==0 then(
1486+
-- if opts.FiberAlgorithm == "markov" then (
1487+
-- for row in entries toricMarkov ((matrix vector val) | A) do(
1488+
-- r := drop(row,1);
1489+
-- if row#0 == 1 and all(r,z->z<=0) then (out = set{-r}; break;);
1490+
-- );
1491+
-- )else (
1492+
-- e:=first exponents (product(for i from 0 to #val-1 list if val#i>0 then ((A.cache#"RingGenerators")#i)^(val#i) else ((A.cache#"RingGenerators")#(i+numRows A))^(-val#i)) % A.cache#"toricIdeal");
1493+
-- if take(e,2*numRows A) == toList((2*numRows A):0) then out = set{take(e,-numColumns A)};
1494+
-- );
1495+
-- );
1496+
-- (A.cache#"fibers")#val = out;
1497+
-- );
1498+
1499+
-- fibAdd = method();
1500+
-- fibAdd (Set,Set) := (L1,L2) -> (
1501+
-- set flatten for l1 in keys L1 list for l2 in keys L2 list l1+l2
1502+
-- );

0 commit comments

Comments
 (0)