---------------------------------------------------------------------
-- gradient
grad = method()
grad RingElement := f -> (
     R := ring f; -- assume R = QQ[x_1,...,x_n,a_1...,a_n]
     n := numgens R // 2;
     transpose (jacobian ideal f)^{0..n-1}
     ) 
----------------------------------------------------------------------
-- processComponent(I,f) --
-- computes the cotangent bundle of a component Y=V(I) relative to f
-- (optional: coDim = codimension of Y, default = the number of gens of I)
----------------------------------------------------------------------
-- output: the defining ideal J_f of the bundle
-- caveat: the correctness is guaranteed ONLY for top-dimensional components 
--         of the projection of V(J_f) 
processComponent = method(Options => {coDim => 0})
processComponent(Ideal, RingElement)  := o -> (I,f) -> (
     R := ring f; -- assume R = QQ[x_1,...,x_n,a_1...,a_n]
     n := numgens R // 2;
     GRAD := grad f;
     Nspace := GRAD;
     scan(flatten entries gens I, g-> Nspace = Nspace||grad g);
     nRows := numgens target Nspace;
     d := codim I;
     -- d := o.coDim;     
     if d == 0 then d = nRows - 1; -- should be the codim of V(I) 
     I0 := minors(d+1, Nspace);
     K := ker map(R^nRows/I, R^n, Nspace);
     J := I + ideal (matrix{{R_n..R_(2*n-1)}}*(gens K));
     satJf := saturate(J, I0);
     out := satJf + ideal f;
     out
     );
----------------------------------------------------------------------
-- correctBundle(B,I,f,J) --
-- takes the bundle B computed for V(I) relative to f and 
-- "corrects" the component V(J) of B
-- output: ideal J' such that V(J')=V(J), but the multiplicity of J' is correct
correctBundle = method()
correctBundle List  := input -> ( --(B,I,f,J) 
     B := input#0; I := input#1; f := input#2; J := input#3;
     R := ring f; -- assume R = QQ[x_1,...,x_n,a_1...,a_n]
     n := numgens R // 2;
     GRAD := grad f;
     Nspace := GRAD;
     scan(flatten entries gens I, g-> Nspace = Nspace||grad g);
     nRows := numgens target Nspace;

     Kf := ker map( R^nRows/(J + ideal f), R^n, Nspace);
     out  := B + ideal (matrix{{R_n..R_(2*n-1)}}*(gens Kf))
     );

---------------------------------------------------------------------
-- old cotanbun
cotanbun = method()
cotanbun RingElement := f -> (
     R := ring f; -- assume R = QQ[x_1,...,x_n,a_1...,a_n]
     processComponent(ideal map(R^0,R^0),f)
     )

---------------------------------------------------------------------
-- multiplicity -- 
-- given a Poincare polynomial for a module outputs its multiplicity
multiplicity RingElement := f->(
     R := ring f;
     while (f % (1-R_0) == 0) do f = f//(1-R_0);
     sub(f, {R_0=>1})
     );

---------------------------------------------------------------------
-- projection onto x_1..x_6
project = method()
project Ideal := I -> (
     R := ring I; -- assume R = QQ[x_1,...,x_n,a_1...,a_n]
     n := numgens R // 2;
     PW :=  QQ[a_1..a_n, x_1..x_n, MonomialOrder=>Eliminate n];
     out := sub(ideal selectInSubring(1, gens gb sub(I,PW)), R);
     use R;
     out
     )

----------------------------------------------------------------------
-- get all the components
BMM = method(Options => {coDim => 0})
BMM(Ideal, RingElement)  := o -> (I,f) -> (
     if BMMstashedCC#?(I,f) then return BMMstashedCC#(I,f);
     if f % I == 0 then return {};
     comps := {}; toCorrect := {};
     B := processComponent(I,f,coDim=>o.coDim);
     dec := ass B; 
     prd := primaryDecomposition B;
     scan(#dec, i->(
	 compI := project dec#i;      
	 projdec := drop(dec,{i,i})/project;
	 if any(projdec, c->isSubset(c,compI)) 
	 then toCorrect = toCorrect | {i} 
	 else (  
	      comps = comps|{compI => (multiplicity@@poincare prd#i //  
	      		multiplicity@@poincare dec#i)};
	      )
	 ));
     while #toCorrect > 0 do (
	  i := first toCorrect; toCorrect = drop(toCorrect,1);
     	  compI := project dec#i;
	  << "-- " << compI << " is a projection of " << dec#i << endl;
	  corB := correctBundle{B,I,f,compI};
     	  cordec = ass corB; 
     	  corprd = primaryDecomposition corB;
	  -- here assume that the largest component comes first
	  comps = comps|{compI => (multiplicity@@poincare corprd#0 //  
		    multiplicity@@poincare cordec#0)}; 
	  );
     comps = {I=>1} | comps ;
     if BMMidealREMEMBER then BMMstashedCC#(I,f) = comps;
     return comps	  
     )
BMM(List,RingElement)  := o -> (cc,f) -> (
     out := {};
     scan(cc, c->if f%c#0!=0 then (
	       comps := BMM(c#0,f,coDim=>codim c#0);
	       scan(comps, comp->(
			 X := select(1,toList(0..#out-1),i->out#i#0==comp#0);
			 if #X==0 then out = out | {comp#0=>c#1*comp#1}
			 else (
			      i := first X;
			      out = take(out,i)|{comp#0=>c#1*comp#1+out#i#1}|drop(out,i+1);
			      )
			 ))  
	       ));
     out
     )

------------------------------------------------------------------------------
-- By default, BMM procedure stashed the result
BMMstashedCC = new MutableHashTable from {}
BMMidealREMEMBER = true;

--------------------------------------------------------------------------
-- "Prunes" a pair of CCs, 
-- i.e: removes identical components (accounting for multiplicity)
pruneCC = method()
pruneCC(List,List) := (A,B) -> (
     nA := {};
     nB := B; 
     scan(A, c->(
	       icc := position(nB, cc -> cc#0==c#0);
	       if icc === null then nA = nA | {c}
	       else (
		    cc :=  nB#icc;
		    m := min(c#1,cc#1); -- take minimum of the multiplicities
		    if c#1-m > 0 then nA = nA | {c#0=>c#1-m};
		    if cc#1-m > 0 then nB = take(nB,icc)|{c#0=>cc#1-m}|take(nB,icc-#nB+1)
		    else nB = take(nB,icc)|take(nB,icc-#nB+1); -- take cc out
		    )
	       )); 
     return (nA,nB)
     ) 

--------------------------------------------------------------------------------
-- Local cohomology computation by "pruning" the CCs of Cech complex 
-- given by a mutable hashtable B
---------------------------------------
-- output: none (B is modified)
PruneCechComplexCC = method() 
PruneCechComplexCC(MutableHashTable) := B -> (
     ind := keys B;
     n := max(ind/(i->#i));
     scan(n, j->scan(ind, a->(
		    if not (set a)#?j then (
			 aj := sort (a|{j});
			 NEWaj := pruneCC(B#a, B#aj);
			 B#a = NEWaj#0; 
			 B#aj = NEWaj#1;
			 )
		    )) ) 
     )