restart
load "it.m2"
W =  QQ[x_1..x_6, a_1..a_6];
I = minors(2, matrix{{x_1, x_2, x_3}, {x_4, x_5, x_6}});
inds = drop(subsets{0,1,2},1) -- nonempty subsets of {0,1,2} 
B = new MutableHashTable from {}
scan(inds, i -> time B#i = cotanbun(product(i,j->I_j))); -- cotanbun for R_{a product of 1,2,or 3 minors} 
ASS = new MutableHashTable from apply(inds, i -> time i=>ass B#i);
PRD = new MutableHashTable from apply(inds, i -> time i=>primaryDecomposition B#i); 
HASS = new MutableHashTable from apply(inds, i -> time i=>ASS#i/poincare);
HPRD = new MutableHashTable from apply(inds, i -> time i=>PRD#i/poincare);
MPRD = new MutableHashTable from apply(inds, i -> i=>HPRD#i/multiplicity);
MASS = new MutableHashTable from apply(inds, i -> i=>HASS#i/multiplicity);

-- projections of all of the components 
PW =  QQ[a_1..a_6, x_1..x_6, MonomialOrder=>Eliminate 6];
--project := I -> ideal selectInSubring(1, gens gb sub(I,PW));
PROJ = new MutableHashTable from apply(inds, i -> time i=>ASS#i/project);

-- computes the multiplicities with which the given component enters every block of the Cech complex
statsForComponent = I -> (
     new HashTable from 
     apply(inds, i->i=>(sum apply(#(PROJ#i), j->(
			 if sub(PROJ#i#j,W) != sub(I,W)then 0
			 else MPRD#i#j//MASS#i#j
		    ))))
     )
-- create the list of all the components (no repetitions)
AC = flatten apply(inds, i -> PROJ#i/(I->sub(I,W)));
allASS = apply( select(toList(0..#AC-1), i->all(take(AC,i-1), I->I!=AC#i)), j->AC#j) 
-- take the alternating sum of multiplicities
sumStats = s -> sum(keys s, i->(-1)^(#i)*s#i);
ASS#{0,1,2}/sumStats@@statsForComponent
scan(allASS, I ->(
	  s := statsForComponent I;
	  ss := sumStats s;
	  <<"Component = V(" << I << ")" << endl 
	  << "entries-> " << s << " alternated sum of multiplicities = " << ss << endl;
	  ))

-- a component is present in the local cohomology modules iff 
-- the alternating sum of its multiplicities in the Cech complex is nonzero
-- (this is true for this particular example, but false in general)
scan(allASS, I ->(
	  s := statsForComponent I;
	  ss := sumStats s;
	  if ss!=0 then 
	  <<"Component = V(" << I << ")" << endl 
	  << "entries-> " << s << " alternated sum of multiplicities = " << ss << endl;
	  ))

-- there are two components above: one accounts for H^3 the other for H^2.