--pairwise Markov property
restart
R1 = QQ[p_(0,0,0,0)..p_(1,1,1,1)]
M1 = matrix{{p_(0,0,0,0),p_(0,0,0,1),p_(0,0,1,0),p_(0,0,1,1)},{p_(1,0,0,0),p_(1,0,0,1),p_(1,0,1,0),p_(1,0,1,1)}}
M2 = matrix{{p_(0,1,0,0),p_(0,1,0,1),p_(0,1,1,0),p_(0,1,1,1)},{p_(1,1,0,0),p_(1,1,0,1),p_(1,1,1,0),p_(1,1,1,1)}}
IP = ideal(det(M1_{0,2}),det(M1_{1,3}),det(M2_{0,2}),det(M2_{1,3}),det(M1_{0,1}),det(M1_{2,3}),det(M2_{0,1}),det(M2_{2,3}))

--global Markov property
IG = minors(2,M1) + minors(2,M2)

--factorization according to the graph
R3 = QQ[p_(0,0,0,0)..p_(1,1,1,1),a_(0,0)..a_(1,1),b_(0,0,0)..b_(1,1,1)]
IF = ideal flatten flatten flatten for i to 1 list for j to 1 list for k to 1 list for l to 1 list p_(i,j,k,l)-a_(i,j)*b_(j,k,l)
JF = eliminate(IF,join(toList(a_(0,0)..a_(1,1)), toList(b_(0,0,0)..b_(1,1,1))))
JF = sub(JF,R1)

--comparison of the ideals
IP == JF
decompose IP
(decompose IP)#0 == JF
IG == JF