-----------------------------------------------------
-- A translation of Luis' Markov package to M2   -----
------------------------------------------------------

-- Give a graph as a hash table i => descendents

-- Make a graph
-- Input: a directed acyclic graph in the form of a 
--        list of lists of children.
--        the vertices must be named 1..n, some n.
--        ASSUMPTION: we assume that the descendents of vertex
--        i are all less than i.  This only represents DAGS.
-- Output: A hashtable G with keys 1..n, and G#i is the
--        the set of all children of the vertex i.
-- This routine produces a useful version of a 'graph'
-- which we use in routines throughout this package.
makeGraph = (g) -> (
     h := new MutableHashTable;
     scan(#g, i -> h#(i+1) = set g#i);
     new HashTable from h)

descendents = (G,v) -> (
     -- returns a set of vertices
     result := G#v;
     scan(reverse(1..v-1), i -> (
	  if member(i,result) then result = result + G#i;
     ));
     result)

nondescendents = (G,v) -> set(1..#G) - descendents(G,v) - set {v}

parents = (G,v) -> set select(1..#G, i -> member(v, G#i))

children = (G,v) -> G#v

  --------------------
  -- Graph features --
  --------------------
  subgraph = (G,vertices) -> (
       -- make a graph on 1..#vertices
       H := new MutableHashTable;
       scan(keys G, i -> H#i = set{});
       scan(#vertices, i -> H#(vertices#i) = set{i+1});
       tr := (s) -> if #s === 0 then {} else toList sum apply(toList s, i -> H#i);
       G' = apply(vertices, i -> tr(G#i));
       makeGraph G')

  H4a = makeGraph {{},{1},{2},{1,3}};
  H4b = makeGraph {{},{1},{1},{2,3}};

  is4cycle = (G) -> G === H4a or G === H4b
       -- G should be a graph on 4 elements


  n4cycles = (G) -> (
       x := subsets(toList(1..#G), 4);
       # select(x, s -> is4cycle subgraph(G,s)))
  ///
      G = makeGraph D5s_65
      G1 = subgraph(G,{1,2,3,4})
      n4cycles G
      G1 === H4b
      subgraph(G,{1,2,3,5})
      sort apply(#D5s, i -> (G := makeGraph D5s#i; {n4cycles G, i}))
      select(0..#D5s-1, i -> (n4cycles makeGraph D5s#i > 0))
      select(0..#D5s-1, i -> (n4cycles makeGraph D5s#i == 3))
      select(0..#D5s-1, i -> #D5s#1>=1 and  # parents(makeGraph D5s#i, 2) >= 3)
  ///
--------------------------
-- Statement calculus ----
--------------------------
-- A dependency is a list {A,B,C}
--  where A,B,C are (disjoint) subsets of positive integers.
--  The meaning is: A is independent of B given C.
-- A dependency list is a list of dependencies

-- No serious attempt is made to remove redundant dependencies.
-- However, we have several very simple routines to remove
-- the most obvious redundant elements

-- If S and T represent exactly the same dependency, return true.
equivStmts = (S,T) -> S#2 === T#2 and set{S#0,S#1} === set{T#0,T#1}

removeRedundants = (Ds) -> (
     -- Ds is a list of triples of sets {A,B,C}
     -- test1: returns true if D1 can be removed
     -- Return a sublist of Ds which removes any 
     --  that test1 declares not necessary.
     test1 := (D1,D2) -> (D1_2 === D2_2 and 
                          ((isSubset(D1_0, D2_0) and isSubset(D1_1, D2_1))
	               or (isSubset(D1_1, D2_0) and isSubset(D1_0, D2_1))));
     -- first remove non-unique elements, if any
     Ds = apply(Ds, d -> {set{d#0,d#1}, d#2});
     Ds = unique Ds;
     Ds = apply(Ds, d -> append(toList(d#0), d#1));
     c := toList select(0..#Ds-1, i -> (
	       a := Ds_i;
	       D0 := drop(Ds,{i,i});
	       all(D0, b -> not test1(a,b))));
     Ds_c)

-- return the set of variables which occur in some
-- dependency.  The input 'Ds' is a list of {A,B,C} 
-- dependencies
occuring = (Ds) -> sum apply(Ds, c -> sum c) 

cleardep = (Ds) -> (
     -- if there is any element which is in every b_i_2
     -- remove it.
     -- This is used to remove a variable if possible.
     -- USE IT WITH CAUTION: it will change the
     -- resulting ideal
     commons := product apply(Ds, c -> c_2);
     apply(Ds, c -> {c_0,c_1,c_2 - commons}))

--------------------------
-- Bayes ball algorithm --
--------------------------
bayesBall = (A,C,G) -> (
     -- A is a set in 1..n (n = #G)
     -- C is a set in 1..n (the "blocking set")
     -- G is a DAG
     -- Returns the subset B of 1..n which is
     --   independent of A given C.
     -- The algorithm is the Bayes Ball algorithm,
     -- as implemented by Luis Garcia, after
     -- the paper of Ross Schlacter
     n := #G;
     zeros := toList((n+1):false);
     visited := new MutableList from zeros;
     blocked := new MutableList from zeros;
     up := new MutableList from zeros;
     down := new MutableList from zeros;
     top := new MutableList from zeros;
     bottom := new MutableList from zeros;
     vqueue := sort toList A;
     -- Now initialize vqueue, set blocked
     scan(vqueue, a -> up#a = true);
     scan(toList C, c -> blocked#c = true);
     local pa;
     local ch;
     while #vqueue > 0 do (
	  v := vqueue#-1;
	  vqueue = drop(vqueue,-1);
	  visited#v = true;
	  if not blocked#v and up#v
	  then (
	       if not top#v then (
		    top#v = true;
		    pa = toList parents(G,v);
		    scan(pa, i -> up#i = true);
		    vqueue = join(vqueue,pa);
		    );
	       if not bottom#v then (
		    bottom#v = true;
		    ch = toList children(G,v);
		    scan(ch, i -> down#i = true);
		    vqueue = join(vqueue,ch);
		    );
	       );
	  if down#v
	  then (
	       if blocked#v and not top#v then (
		    top#v = true;
		    pa = toList parents(G,v);
		    scan(pa, i -> up#i = true);
		    vqueue = join(vqueue,pa);
		    );
	       if not blocked#v and not bottom#v then (
		    bottom#v = true;
		    ch = toList children(G,v);
		    scan(ch, i -> down#i = true);
		    vqueue = join(vqueue,ch);
		    );
	       );
	  ); -- while loop
     set toList select(1..n, i -> not blocked#i and not bottom#i)
     )
--------------------------
-- Markov relationships --
--------------------------
pairMarkov = (G) -> (
     -- given a graph G, returns a list of triples {A,B,C}
     -- where A,B,C are disjoint sets, and for every vertex v
     -- and non-descendent w of v,
     -- {v, w, nondescendents(G,v) - w}
     flatten apply(1..#G, v -> (
	       ND := nondescendents(G,v);
	       W := ND - parents(G,v);
	       apply(toList W, w -> {set {v}, set{w}, ND - set{w}}))))
			 
localMarkov = (G) -> (
     -- Given a graph G, return a list of triples {A,B,C}
     -- of the form {v, nondescendents - parents, parents}
     result := {};
     scan(1..#G, v -> (
	       ND := nondescendents(G,v);
	       P := parents(G,v);
	       if #(ND - P) > 0 then
	         result = append(result,{set{v}, ND - P, P})));
     result)

globalMarkov = (G) -> (
     -- Given a graph G, return a complete list of triples {A,B,C}
     -- so that A and B are d-separated by C (in the graph G).
     -- If G is large, this should maybe be rewritten so that
     -- one huge list of subsets is not made all at once
     n := #G;
     vertices := toList(1..n);
     result := {};
     AX := subsets vertices;
     AX = drop(AX,1); -- drop the empty set
     AX = drop(AX,-1); -- drop the entire set
     scan(AX, A -> (
	       A = set A;
	       Acomplement := toList(set vertices - A);
	       CX := subsets Acomplement;
	       CX = drop(CX,-1); -- we don't want C to be the entire complement
	       scan(CX, C -> (
			 C = set C;
			 B := bayesBall(A,C,G);
			 if #B > 0 then (
			      B1 := {A,B,C};
			      if all(result, B2 -> not equivStmts(B1,B2))
			      then 
			          result = append(result, {A,B,C});
	       )))));
     result
     )
---------------------------
-- Constructing the ring --
---------------------------
probring = d -> (
     -- Returns a ring in d1*...*dn variables
     -- d should be a list of integers (d1, d2, ..., dn)
     -- where each di >= 1
     start := (#d):1;
     QQ[(symbol p)_start .. (symbol p)_d])

probring1 = (parents, d) -> (
     -- Returns a ring with variables
     --   a_(i1,...,ir) and q_(j2,...,jn), where
     -- r-1 is the number of parents of 1
     -- each j_p ranges over 1..d_p
     -- each i_p ranges over 1..(d_(parents_p))
     -- parents should be a list of indices 1..#d which 
     -- are the parents of vertex 1.
     -- d should be a sequence of integers (d1, d2, ..., dn)
     -- where each di >= 1
     start := (#d-1):1;
     lastd := drop(d,1);
     startpa := (#parents+1):1;
     lastpa := toSequence(prepend(d#0, apply(parents, i -> d#(i-1))));
     QQ[(symbol a)_startpa .. (symbol a)_lastpa, (symbol q)_start .. (symbol q)_lastd])

toring1 = (parents,d,R,R1) -> (
     -- R should be probring d
     -- R1 should be probring(parents,d)
     -- Constructs a ring map F : R --> R1
     -- which is the partial factorization map
     F := toList apply(((#d):1) .. d, i -> (
	       pa := toSequence(prepend(i#0, apply(parents, j -> i#(j-1))));
	       a_pa * q_(drop(i,1))));
     map(R1,R,F)
     )

margin1 = (d,R) -> (
     -- R should be probring d.
     -- construct the ring map which replaces p_(1,i2,...,in)
     --   with the sum of all p_(j,i2,...,in), all j.
     F := toList apply(((#d):1) .. d, i -> (
	       if i#0 > 1 then p_i
	       else (
		    i0 := drop(i,1);
		    p_i - sum(apply(toList(2..d#0), j -> p_(prepend(j,i0)))))));
     map(R,R,F))

margin = (v) -> (d,R) -> (
     -- R should be probring d.
     -- construct the ring map which replaces p_(1,i2,...,in)
     --   with the sum of all p_(j,i2,...,in), all j.
     v = v-1;
     F := toList apply(((#d):1) .. d, i -> (
	       if i#v > 1 then p_i
	       else (
		    i0 := drop(i,1);
		    p_i - sum(apply(toList(2..d#v), j -> (
			      newi := join(take(i,v), {j}, take(i,v-#d+1));
			      print p_newi;
			      p_newi))))));
     map(R,R,F))
-------------------------------------------------------
-- Constructing the ideal of a independence relation --
-------------------------------------------------------
cartesian = (L) -> (
     if #L == 1 then 
	return toList apply (L#0, e -> singleton e);
     L0 := L#0;
     Lrest := drop (L,1);
     C := cartesian Lrest;
     flatten apply (L0, s -> apply (C, c -> prepend (s,c))))

possibleValues = (d,A) ->
     cartesian (toList apply(1..#d, i -> 
	       if member(i,A) 
	       then toList(1..d#(i-1)) 
	       else {0}))

prob = (d,s) -> (
     L := cartesian toList apply (#d, i -> 
	   if s#i === 0 
	   then toList(1..d#i) 
	   else {s#i});
     sum apply (L, v -> p_v))

indepMatrices = (d,A,B,C) -> (
     -- d is a sequence of integers
     -- A,B,C should be disjoint sets of integers in 1..#d
     Avals = possibleValues(d,A);
     Bvals = possibleValues(d,B);
     Cvals = possibleValues(d,C);
     apply(Cvals, c -> (
       matrix apply(Avals, a -> apply(Bvals, b -> (
		      e := toSequence(toList a + toList b + toList c);
		      prob(d,e))))))
     )

indepIdeal = (d,A,B,C) -> 
     sum apply(indepMatrices(d,A,B,C), m -> minors(2,m))

Ideps = (d, Ds) -> (
     R = probring d;
     sum apply(Ds, D -> indepIdeal(d,D_0,D_1,D_2)))

-----------------------------------------------------
-- local Markov ideal
localMarkovIdeal = (d,gr) -> (
     R = probring d;
     G := makeGraph gr;
     D := localMarkov G;
     D = removeRedundants D;
     sum apply(D, D0 -> indepIdeal(d, D0_0,D0_1,D0_2)))

----------------------
-- Binary variables --
----------------------
  -- binaryR#0, binaryR#1 are not very useful...
binaryD = apply(0..7, d -> toSequence(d:2));
binaryR = apply(binaryD, probring);

simplifyDeps = (b) -> (
     b = removeRedundants b;
     b = cleardep b;
     s := sort toList occuring b;
     d := #s;
     h := new MutableHashTable;
     scan(d, i -> h#(s#i) = i+1);
     g := b/(c -> (
	       -- c is a list of three sets
	       {set apply(toList c#0, j -> h#j),
		set apply(toList c#1, j -> h#j),
		set apply(toList c#2, j -> h#j)}));
     R = binaryR#d;
     use R;
     g)     

IBinaryDeps = (b) -> (
     -- Given a set of dependencies, compute the
     -- ideal in the smallest possible ring.
     -- IE, if any variables don't occur, use a smaller
     -- size ring.
     s := sort toList occuring b;
     d := #s;
     h := new MutableHashTable;
     scan(d, i -> h#(s#i) = i+1);
     g = b/(c -> (
	       -- c is a list of three sets
	       {set apply(toList c#0, j -> h#j),
		set apply(toList c#1, j -> h#j),
		set apply(toList c#2, j -> h#j)}));
     R = binaryR#d;
     use R;
     << "using" << endl;
     g/print;
     sum apply(g, D -> indepIdeal(binaryD#d,D_0,D_1,D_2)))

localMarkovConditions = (gr) -> (
     D := localMarkov makeGraph gr;
     --simplifyDeps D
     removeRedundants D
     )

globalMarkovConditions = (gr) -> (
     D := globalMarkov makeGraph gr;
     --simplifyDeps D
     removeRedundants D
     )
     
localIdeal = (gr) -> (
     D := localMarkov makeGraph gr;
     D = simplifyDeps D;
     d := #(occuring D);
     J := sum apply(D, c -> indepIdeal(binaryD#d, c_0,c_1,c_2));
     F := margin1(binaryD#d, binaryR#d);
     F J)

globalIdeal = (gr) -> (
     D := globalMarkov makeGraph gr;
     D = simplifyDeps D;
     d := #(occuring D);
     J := sum apply(D, c -> indepIdeal(binaryD#d, c_0,c_1,c_2));
     F := margin1(binaryD#d, binaryR#d);
     F J)

Iglobal = (D) -> (
     -- D is a list of dependencies
     -- d is the largest variable which occurs
     d := max apply(D, b -> max join(apply(b, c -> max toList c)));
     J := sum apply(D, c -> indepIdeal(binaryD#d, c_0,c_1,c_2));
     F := margin1(binaryD#d, binaryR#d);
     F J)
///
-- tests and examples of using these routines
restart
load "markov.m2" -- this file
G = makeGraph{{}, {1}, {2}, {1, 2, 3}, {4}}
parents(G,1)
parents(G,2)
parents(G,3)
parents(G,4)
parents(G,5)
descendents(G,5)
descendents(G,4)
descendents(G,3)
nondescendents(G,3)
G

scan(pairMarkov G, D -> print toString D)
scan(localMarkov G, D -> print toString D)
scan(globalMarkov G, D -> print toString D)
scan(removeRedundants globalMarkov G, D -> print toString D)
Ds = removeRedundants globalMarkov G;
occuring Ds

R = probring(2,3,3,3)
numgens R
gens R

-- Here are useful already created values for 
-- up through 7 random variables
binaryD_5
binaryR_5

R = binaryR_4

-- Marginalizing over vertex 1.
-- This often simplifies the Markov ideals.
F = margin1((2,2,2,2),R)
transpose F.matrix
F p_(1,1,1,1)

-- The partial factorization map
R1 = probring1({2,3},(2,2,2,2))
  -- {2,3} is the list of parents of 1.
F = toring1({2,3},(2,2,2,2),R,R1)
transpose F.matrix

-- prob
use R
gens R
prob((2,2,2,2),{0,1,1,0})

-- making the ideals/matrices
-- Use the binary 4 variable ring R from above
indepMatrices((2,2,2,2),set{1,2},set{3},set{4})
J = indepIdeal((2,2,2,2),set{1,2},set{3},set{4})
(codim J, degree J)

-- local Markov ideal, starting with an externally
-- constructed graph (ie: a list of children.)
-- This defines the ring too, removes redundant
-- variables.  WARNING: it might remove random
-- variables which appear as the dependent in every dependency.
-- (This simplifies the ideal).

J = localMarkovIdeal((2,2,2,2),{{},{1},{2},{3}})
transpose gens J

G = makeGraph{{},{1},{2},{3}}
D = localMarkov G
J1 = Ideps D


R = probring(2,2,3)
R = probring(4:2)
prob((2,2,2,2),{0,1,0,2})
gens R
#oo
cartesian{{1,2},{1,2},{1,2,3}}
possibleValues((2,2,2,2),set(1,3))
indepMatrices((2,2,2,2),set{1},set{2},set{3,4})
J = indepIdeal((2,2,2,2),set{1},set{2},set{3,4})

-- Test of Bayes Ball
G = makeGraph{{}, {1}, {2}, {1, 2, 3}, {4}}
Ds = globalMarkov G
oo/print;
removeRedundants Ds
oo/print;
bayesBall(set{1},set{2},G)

restart
load "luis.m2"  -- this is THIS file
J = localIdeal{{}, {1}, {2}, {1, 2, 3}, {4}};
JG = globalIdeal{{}, {1}, {2}, {1, 2, 3}, {4}};
F = margin1(binaryD#4, R)
J = F J
JG = F JG;
apply(flatten entries gens JG, f -> (<< "." << flush; J : f));
intersect oo;
intersect(o35,JG) == J -- yes
C = primaryDecomposition o35;
transpose gens J
transpose gens JG
load "facGB.m2"
load "primdecomp-tools.m2"
JG1 = ideal compress((gens JG) % J);
JG1 = flatten entries gens JG1;
JG1 = ideal select(JG1, f -> size f === 2);
transpose gens JG1
J1 = J : JG1;
J1 = J : JG_0;
J2 = J1 : JG_1;
L = apply(flatten entries gens JG, f -> J : f);
info J
info JG
J1 = J : JG
J = localMarkovIdeal((2,2,2,2,2), {{}, {1}, {2}, {1, 2, 3}, {4}});


info J
(margin1((2,2,2,2),R)) J
transpose gens oo

R = probring(2,2,2,2,2)
R1 = probring1({2,3},(2,2,2,2,2))
toring1({2,3},(2,2,2,2,2),R,R1)
margin1((2,2,2,2,2),R)
transpose oo.matrix

D1 = {set{1},set{2},set{3,4}}
D2 = {set{1},set{2},set{3}}
R = probring(2,2,2,2)
J1 = indepIdeal((2,2,2,2),D1_0,D1_1,D1_2)
J2 = indepIdeal((2,2,2,2),D2_0,D2_1,D2_2)
codim J1
codim J2
(gens J2) % J1
(gens J1) % J2
///