-- A moment graph is a hash table G with three kinds of entries:
--   1. G.HTpt : a ring, the symmetric algebra over a vector space/
--      lattice.  The notation means "equivariant cohomology of a point"
--   2. G.vert : a list of the names of vertices of the graph.  
--        There is a partial order on the set of vertices, and
--        it is assumed that if v > w, then v appears before w
--        in G.vert.
--   3. G.ind.v = the index of a vertex v, to speed up finding
--        vertices.  I.e., G.vert_(G.ind.v) = v.
--   4. G.edges : a list of edges.  The edges are of the form
--        {i, j, l}, where i and j are the indexes of the endpoints of 
--        the edge with i < j, and l is an element of A of degree 1, 
--        representing the direction of the edge.  

MomentGraph = new Type of MutableHashTable

--     
-- A MomentGraph Sheaf is a hash table S with entries:
--   1. S.graph = G : a MomentGraph on which the sheaf is defined
--   2. S.gs : The module of global sections, a (free) module over G.HTpt
--   3. S.restr : table of restriction homomorphisms; for the i-th
--        vertex v, S.restr.i is a (surjective) homomorphism 
--        from S.gs onto a free G.HTpt - module which is the stalk
--        of S at v. 

MomentGraphSheaf = new Type of MutableHashTable;


makeMomentGraph = {Edges => {}} >> o -> (R, V) -> (
     G := new MomentGraph;
     
     if (class V) === ZZ then
        V = toList(0..(V-1));	
     
     G.vert = V;
     G.HTpt = R;
     G.edges = {};
     
     G = addEdges(G, o.Edges);
      
     G.ind = new MutableHashTable;
     scan(#V, i -> (G.ind#(V_i) = i));     
     
     G
     );
    

findVertex = (G, v) -> (
     if (class v) === ZZ then return v
     else (
	  P = positions(G.vert, x -> (x == v));
	  if P == {} then
	     error("vertex not found.")
	  else 
	     return P_0;)
     );

addVertex = (G, v) -> (
     if (G.ind#?v) then
        error("vertex already exists")
     else (
	G2 := new MomentGraph from G;  
	G2.ind#v = #G2.vert;  
	G2.vert = append(G2.vert, v)   
        );
     G2
     );   

addEdge = (G, v, w, l) -> (
     v = findVertex(G, v);
     w = findVertex(G, w);
     
     l = substitute(l, G.HTpt);
     
     G2 := new MomentGraph from G; 
     
     if v < w then
        G2.edges = append(G2.edges, {v, w, l})
     else
        G2.edges = append(G2.edges, {w, v, l});
     G2
     );

addEdges = (G, L) -> (
    G2 := new MomentGraph from G;
    scan(L, x -> (G2 = addEdge(G2, x_0, x_1, x_2)));
    G2
    );
   
   
show(MomentGraph) := G -> (
     << "Coefficient Ring: " << G.HTpt << endl;
     << "Vertices:" << endl;
     scan(#G.vert, i -> (<< i << ": "<< G.vert#i << endl));
     << endl << "Edges" << endl;
     
     scan(G.edges, e -> (
	  << G.vert_(e_0) << " -- " << G.vert_(e_1);
	  << ":  " << e_2 << endl << endl;));
     ); 
     
show(MomentGraphSheaf) := S -> (
     << "Global Sections: " << S.gs << endl;
     << "Poincare' polynomial: " << poincare S.gs << endl;
     << "Stalks:" << endl;
     scan(#S.graph.vert, i -> (
	       L := flatten degrees target S.restr#i; 
	       << S.graph.vert#i;
	       << ": rank " << #L;
	       << ", degrees " << L;
	       << endl;));
     ); 

-- BruhatGraph(w): generate the moment graph corresponding to the
-- Bruhat interval [id,w].  An optional argument lowerVertex => x
-- allows you to calculate just [x,w].

BruhatGraph = {CoxeterType => "A", lowerVertex => {}, groundRing => QQ} >>
    opts -> upperVertex -> (
	 if opts.CoxeterType != "A" then 
	    error "Only type A Bruhat graphs implemented so far";
	    
	 P := upperVertex;
	 	 
	 G := new MomentGraph;
	 
	 i := 0;
	 V := {};
	 
	 A := G#HTpt = 
	    opts.groundRing[(symbol x)_0 .. (symbol x)_(#P-1)];  -- the coefficient ring H_T(pt)

	 L := {P};            -- list of permutations in the Bruhat interval.
         G#ind = new MutableHashTable;
	 G#ind#P = 0;          -- we step through L, adding all codimension
			       -- one "children" before moving to the next
			       -- vertex.  This ensures that L is ordered
			       -- without having to explicitly sort them.
	 G#edges = {};	 
	  
	 aboveBottom := true;    -- a flag 
	  
         k := 0; j := 0;
         while k < #L do (
   
            i = 0;
   	    while i < #P - 1 do ( 
      	       j = i + 1;
       	       while j < #P do (
             	        
		  T = ApplyReflection(L#k, i, j);
		  
     	       	  aboveBottom = true;		    
		  if opts.lowerVertex != {} then                       -- fixed by Friederike Steglich
                      aboveBottom =
                          all( 0..#P-2, l->
                             all(sort(T_{0..l}) - sort((opts.lowerVertex)_{0..l}),
                                  x -> (x >= 0))==true);

                  if L#k#i > L#k#j and aboveBottom then (  
		                            -- transposing i and j takes 
		                            -- us down, but above bottom
     	       	     G#edges =                         -- of interval
		        append(G#edges, {L#k, T, A_i - A_j});
			
                     betw = L#k_{i+1 .. j-1};
                     if all(betw, x -> (x > L#k#i) or (x < L#k#j))
                        and not member(T,set L) then (
			   G#ind#T = #L;          
                           L = append(L, T);
			   );                     
                                            -- add the perm if it's not 
					    -- already there and is codim 1 
                     );
                  
		  j = j + 1; 
                  );
	     
               i = i + 1;
               );
            k = k + 1;
            );
         G#vert = L;
	 
	 G#edges = apply(G#edges,
	      E -> ({G#ind#(E_0), G#ind#(E_1), E_2}));
       
         G
         );       
        
ApplyReflection = (x,i,j) -> (
       xs := new MutableList from x;
       temp := xs#i;
       xs#i = xs#j;
       xs#j = temp;   xs = toList xs;
       xs
       );
  
ApplySimpleReflection = (x,s) -> 
     ApplyReflection(x, s, s+1);



-- ICSheaf(G): generate the IC sheaf for a MomentGraph G, starting
-- with the first vertex, unless the optional argument topVertex is
-- specified.
     
ICSheaf = {topVertex => 0, maxDegree => infinity} >> 
     opts -> (G) -> (
	  assert(class G === MomentGraph);
	  
	  S := new MomentGraphSheaf;
	  S.graph = G;

	  S.restr = new MutableHashTable;
	  
	  A := G.HTpt;

     	  i := findVertex(G, opts.topVertex);

	  j := 0;	     
	  while j < i do (
	      S.restr#j = map(A^0,A^1,0);
	      j = j + 1;
	      ); 	
		
          S.gs = A^1;                  -- global sections module
          S.restr#i = id_(A^1);  

          while i < #(G.vert)-1 do (
   	       i = i + 1;
	       calculateICVertex(S, i, opts.maxDegree);
       	       );
	  S
	  );
     
calculateICVertex = (S, i, maxDegree) -> (     
   E := (S.graph).edges;

   edgesUp := E_(positions(E, e -> (e_1 == i)));
   
   MapsToEdges := apply(edgesUp, e -> (
	r := S.restr#(e_0);     
	inducedMap((target r)/ e_2, (target r)) * r)
        );   

   g := fold( (x,y) -> (x || y), MapsToEdges);

   Imu := image g;
   Imp := minimalPresentation Imu;
   pm := Imp.cache.pruningMap;

   Mx := cover Imp;

   f := -- inducedMap(target g, Imu) * 
        inducedMap(target g, target pm) * 
	pm * inducedMap(source pm, Mx);
	

--   f :=  inducedMap(target g, Imu) * 
--	 inducedMap(Imu, Mx, Imp.cache.pruningMap);

   Kerf := kernel f;            
   MxKernel := prune Kerf;          -- start new sections
   P := (positions(degrees MxKernel,
		      x -> (first x <= maxDegree)));

   S.restr#i = (g // f) | 
		(inducedMap(Mx, Kerf) * MxKernel.cache.pruningMap * MxKernel_P);
   MxKernel = source MxKernel_P;

   S#gs = source S.restr#i;

   if #P > 0 then (
	j = 0;
	while j < i do (
	     S.restr#j = S.restr#j | 
		     map(target S.restr#j, MxKernel, 0);
	     j = j + 1;
	     );
	);
   );

constSheaf = G -> (
     assert(class G === MomentGraph);
     n := #G.vert;

     S := new MomentGraphSheaf;
     S.graph = G;
     
     S.restr = new MutableHashTable;
     
     A := G.HTpt;
     
     M := A^n;
     f := map(A^0,M, {});
     
     i := 0;

     scan(G.edges, e -> (
	       f = f || map(A^1/e_2, M, M^{e_0} - M^{e_1})));

     S.gs = ker f;
     K = source gens S.gs;	       
     
     apply(n, i -> (S.restr#i = K^{i}));
     
     S
     );
     
restrictionMap = S -> (
     v := null;
     
     if (class S === MomentGraphSheaf) then
        v = toList(0..(#S.graph.vert-1))
     else (
	v = {findVertex((S_0).graph, S_1)};
	S = S_0
	);
     
     f:= fold((A,B) -> A || B, apply(v,
	       i -> S.restr#i));
     map(directSum(apply(v, i -> target S.restr#i)), target f, 1) * f       
     );

stalk = (S,v) -> (target restrictionMap(S,v));

globalSections = method()
globalSections(MomentGraphSheaf) := S -> S.gs;

mult = (A, B) -> (i := rank target A;      
     assert( i == rank target B);
     assert( (rank source A) == 1);
     assert((rank source B) == 1);
     
     fold((A, B) -> A || B, 
            apply(0..(i-1), j -> A_{0}^{j} * B_{0}^{j}))
     );

HT = G -> (
     assert( class G === MomentGraph);
     
     S := constSheaf(G);
     M := gens S.gs;
     D := degrees source M;
     n := #D;
     m := #degrees S.graph.HTpt;
     
     
     R := QQ[gens S.graph.HTpt, (symbol a)_1 .. (symbol a)_(n-1), 
	  Degrees => {m:{1}} | D_{1..(n-1)}];
     inc := map(R, S.graph.HTpt);
     
     
     L := flatten apply(toList(1..(n-1)), i-> 
	        apply(toList(i..(n-1)),j ->
	            a_i * a_j - 
	            (vars R)_{m..(m+n-2)} * inc(mult(M_{i}, M_{j})//M)^{1..(n-1)}));
     R / (trim (ideal L))   
     );

cohomologyRing = G -> (
          assert( class G === MomentGraph);
     
     S := constSheaf(G);
     M := gens S.gs;
     D := degrees source M;
     n := #D;
     m := #degrees S.graph.HTpt;
     
     
     R := QQ[(symbol a)_1 .. (symbol a)_(n-1), 
	  Degrees => D_{1..(n-1)}];
     inc := map(R, S.graph.HTpt);
     
     
     L := flatten apply(toList(1..(n-1)), i-> 
	        apply(toList(i..(n-1)),j ->
	            a_i * a_j - 
	            (vars R) * inc(mult(M_{i}, M_{j})//M)^{1..(n-1)}));
     R / (trim (ideal L))   
     
     );

multMatrix = (R, lambda, i) -> (
     Bi := basis(i, R);
     Biplus1 := basis(i+1, R);
     (Bi * lambda) // Biplus1
     );

isWeakLef = (R, lambda) -> (
     assert(dim R == 0);
     i := 0;
     done := false;
     
     while not done do (
	  M := multMatrix(R, lambda, i);
	  if (rank M) < min(rank source M, rank target M)
	       then return false;
	  i = i + 1;
	  if (rank target M == 0) then done = true;
	  );
     return true;
     );

isLefElt = (R, lambda) -> (
     assert(dim R == 0);
      
     L := {};
     front := {};
     back := {};     
     i := 0;
     done := false;
     M := {};
          
     while not done do (
	  M = multMatrix(R, lambda, i);
	
     	  if (rank target M == 0) then done = true
	  else (
	       L = append(L, M);
	       i = i + 1;
	       );
	  );
     
     if i == 0 then return true;
     if i == 1 then return isIsomorphism L_(0);    
     
     if even i then (
	  M = L_(i//2) * L_(i//2 - 1);
	  front = reverse L_{0..(i//2 - 2)};
	  back = L_{(i//2 + 1) .. (i-1)}
	  )
        else (
	     M = L_(i//2);
	     front = reverse L_{0..(i//2 -1)};
	     back = L_{(i//2 + 1) .. (i-1)} 	     
	     );
 
     done = false;
     
     while not done do (
	  if not (isIsomorphism M) then
	     return false;
	  if (back == {}) then
	     return true;    
	  M = (first back) * M * (first front);   
          back = drop(back, 1);
	  front = drop(front, 1);
	  ); 
     );

hasLefschetzElt = (R) -> (
     L := multMaps(R);
     i := #L;
     M := {};
     
     if i == 0 then return true;
     if i == 1 then return (rank target L_0 == 1);
     
     if even i then (
	  M = L_(i//2) * L_(i//2 - 1);
	  front = reverse L_{0..(i//2 - 2)};
	  back = L_{(i//2 + 1) .. (i-1)}
	  )
        else (
	     M = L_(i//2);
	     front = reverse L_{0..(i//2 -1)};
	     back = L_{(i//2 + 1) .. (i-1)} 	     
	     );
 
     done = false;
     
     d := {};
     
     while not done do (
	  if (rank source M != rank target M) then
	     return false;
	  d = det(M);   
	  if (d == 0_(class d)) then
	     return false;   
	  if (back == {}) then
	     return true;    
	  M = (first back) * M * (first front);   
          back = drop(back, 1);
	  front = drop(front, 1);
	  ); 
     );

LefschetzIdeal = (R) -> (
     L := multMaps(R);
     i := #L;
     
     M:= {};
     
     if i == 0 then return ideal 1_QQ;
     if i == 1 then return ideal det(L_0);
     
     if even i then (
	  M = L_(i//2) * L_(i//2 - 1);
	  front = reverse L_{0..(i//2 - 2)};
	  back = L_{(i//2 + 1) .. (i-1)}
	  )
        else (
	     M = L_(i//2);
	     front = reverse L_{0..(i//2 -1)};
	     back = L_{(i//2 + 1) .. (i-1)} 	     
	     );
 
     done = false;
     
     I := ideal(1_(ring L_0));     
     
     while not done do (
	  if (rank source M != rank target M) then (
	     << "Error: ring must be symmetric" << endl;
	     return false;
	     );
	  I = intersect(I, ideal det(M));   
	  if (back == {}) then
	     return trim I;    
	  M = (first back) * M * (first front);   
          back = drop(back, 1);
	  front = drop(front, 1);
	  ); 
     );


multMaps = (R) -> (
            assert(dim R == 0);
	    
	    
	    i := 0;
	    
	    done := false;
     	    L := {};
	    M := {};
	    
	    B1 := basis(1, R); 
	    r1 := rank source B1;	    

	    S := QQ[(symbol a)_1 .. (symbol a)_(r1)];

	    Bi := {};
	    ri := 0;
	    	    
	    Biplus1 := basis(0,R);
	    riplus1 := rank source Biplus1;
	              
     	    while not done do (
	        Bi = Biplus1;
                ri = riplus1;
	        
	        Biplus1 = basis(i+1, R);
                riplus1 = rank source Biplus1;	    
	    
	        if riplus1 == 0 then
		   done = true
		else (   
                   M = substitute((B1 ** Bi)//Biplus1, QQ);
	           L = append(L,
			substitute(M, S) * ((transpose basis(1, S)) ** map(S^ri, S^ri, 1)));
		   i = i + 1;
		   );
	      );
	 
	    L
	    );


edgeMap = (S, v1, v2) -> (
     v1 = findVertex(S.graph, v1);
     v2 = findVertex(S.graph, v2);
     
     P = positions(S.graph.edges, e -> 
	              ((e_0 == v1 and e_1 == v2) or
		       (e_1 == v1 and e_0 == v2)));
     if (P == {}) then error("edge not found.");
     
     e := S.graph.edges#(P_0);
     l := e_2;
     f0 := S.restr#(e_0);
     f1 := S.restr#(e_1);
     
     map((target f0)/l, target f1, f0 * (id_(target f1) // f1))
     );
     
     
Hom(MomentGraphSheaf, MomentGraphSheaf) := (S1, S2) -> (
     if S1.graph =!= S2.graph then
	error("cannot take Hom of sheaves on different graphs.");
     
     r1 := restrictionMap(S1);
     r2 := restrictionMap(S2);
     
     GS1 := source r1;
     GS2 := source r2;
     VertexSum1 := target r1;
     VertexSum2 := target r2; 
     v := #S1.graph.vert;
     
     f := Hom(GS1, r2);
     g := Hom(r1, VertexSum2);
     
     h := fold((x, y) -> (x | y), 
             apply(v, 
                   i -> ( Hom(VertexSum1^[i], VertexSum2) *
                         Hom((components(VertexSum1))#i, VertexSum2_[i]))));
     
     j :=  modulo(f, g * h);
     map(Hom(GS1, GS2), source j, j) 
     );
     

Hom(MomentGraphSheaf, MomentGraphSheaf, ZZ) := (S1, S2, deg) -> (
     if S1.graph =!= S2.graph then
	error("cannot take Hom of sheaves on different graphs.");
     
     r1 := restrictionMap(S1);
     r2 := restrictionMap(S2);
     
     GS1 := source r1;
     GS2 := source r2;
     VertexSum1 := target r1;
     VertexSum2 := target r2; 
     v := #S1.graph.vert;
     
     f := Hom(GS1, r2);
     f = f_(positions(degrees source f, a -> a_0 <= deg));
     g := Hom(r1, VertexSum2);
     
     h := fold((x, y) -> (x | y), 
             apply(v, 
                   i -> ( Hom(VertexSum1^[i], VertexSum2) *
                         Hom((components(VertexSum1))#i, VertexSum2_[i]))));
     
     j :=  modulo(f, g * h);
     map(Hom(GS1, GS2), source j, j)      
     );
              


getPoly = method()
getPoly(String) := filename -> (
     S := lines get("!polymake " | filename | " VERTICES");
      
     n := #separate(" ",S_1) - 1;
     
     S = S_{1..(#S-2)};
     
     
     A := QQ[(symbol x)_0..(symbol x)_(n-1)];
     
     
     V := apply(S, s -> (
	     (apply(separate(" ",s), value))_{1..n}));
     
     G := makeMomentGraph(A, V);     

     vert := apply(V, v -> sum(v, 0..(n-1), (i,j)->(i*A_j)));
     
     S = lines get("!polymake " | filename | " HASSE_DIAGRAM"); 
     I := value("{" | addCommas S_1 | "}");
     
     S = S_{3..(#S-3)};

     S = apply(S, s -> (value addCommas s));
     
     i := position(I, j -> #(S_(j-1))_0 == 2);
     
     start := I_i - 1;
     finish := I_(i+1) - 2; 
         
     scan(S_{start..finish}, s -> 
	  G = addEdge(G, s_0_0, s_0_1, vert_(s_0_0) - vert_(s_0_1)));
   
     G
     );

addCommas = s -> (
     ss := characters s;
     L := positions(ss, c -> (c == " "));
     
     ss = new MutableList from ss;
     
     scan(L, i -> (ss#i = ","));
     concatenate ss
     );

--
-- Routines to handle general reflection graphs.
--


dot = (x,y) ->                               -- dot product
   sum(x,y, (a,b)-> a*b);
   
reflect = (r,v) ->                        -- reflect v in hyperplane
    (v - 2 * (dot(r,v)/dot(r,r)) * r);    -- perpendicular to r
    
    
    
reflectionArrangement = {maxN => 100} >> 
   opts -> (R) -> (	
     V := set {};
     newV := set apply(R, v -> toVector v);
     
     while (V != newV) do (
	  V = newV;
	  newV = V + set flatten(
	       apply(toList V, w-> apply(toList newV, r-> reflect(r,w))));
	  if #newV > opts.maxN then
	     error("Too many vectors -- maybe an infinite arangement?");
	  );
     toList V
     );   
    
    
reflectionOrbit = (R,v) -> (
     
     V := set {};
     newV := set {toVector v};
     
     
     while (V != newV) do (
	  V = newV;
          newV = V + set flatten( 
	             apply(toList V, w -> apply(R, r-> reflect(r,w)))); 
	  );
     V
     );

reflectionGraph = (R,v) -> (
     G := new MomentGraph;
     
     v = toVector v;
     G#HTpt = (ring v_0)[x_0..x_(#v-1)];
     G#vert = rsort toList reflectionOrbit(R,v);
     
     G#ind = new MutableHashTable;
     
     w := {};
     apply(#(G#vert), i -> (
	       w = G#vert_i; 
	       G#ind#w = i));

     G#edges = toList set flatten apply(G#vert, w -> (
	       apply(R, r -> (rw := reflect(r,w);
			 {G#ind#w, G#ind#rw, sum(#w, i-> (w_i - rw_i) * G#HTpt_i)}))));    

     G#edges = select(G#edges, e -> (e_0 < e_1));
     G
     );

toVector = v -> (apply(v, a -> a/1));   -- make sure vector entries are (at least) in QQ, not in ZZ.

Set == Set := (S,T) -> (isSubset(S,T) and isSubset(T,S));