#########################################################################
# mhomo-unmixed.mpl                                                     #
# Maple V code under development; the author assumes no responsibility  #
# but comments are welcome to amantzaf@math.uoa.gr                      #
# created 06/2007, last modified 10/2008                                #
# Authors: Ioannis Emiris & Angelos Mantzaflaris, 2009                  #
# Implementing scaled multihomogeneous resultant matrix formulae        #
#########################################################################

with(linalg):
with(combinat):

#########################################################################
#########################################################################
#
# CONTENTS (in this order)
#
# mBezout          - computes the multihomogeneous bezout bound
# computeP         - computes the union of Pk sets
# coHzero          - true if cohomology H^q(m-pd)=0 is zero
# coHdim           - dimension of H^q(m-pd)=0
# dimKn            - dimension of term K_{n}
# dimKnp           - dimension of summand K_{n,p}
# critVect         - critical degree vector
# allDetVecs       - search for all determinantal vectors
# detboxes         - determinantal boxes
# has_determ       - true if data is determinantal
# allvars          - return all the variables of the system
# mdegree          - multidegree of a polynomial
# monbasis         - momomial basis of S(u)
# multmap          - Multiplication matrix of M(f_i)
# Sdim             - dimension of S(u)
# Sylvmat          - Sylvester matrix
# BezoutianPoly    - Bezoutian block
# Bezoutmat        - Bezoutmatrix
# makeMatrix       - matrix of K_1 -> K_0
# Makepoly         - make a multihomogeneous polynomial
# MakeExtremePoly  - make a multihomogeneous polynomial
#                    with only extreme terms
# Makesystem       - make a multihomogeneous system
# lstmonof         - list monomials of a polynomial
# MakeComplex      - compute the Weyman complex of an m-vector
# MakeKv           - make term K_{\nu} of the complex
# makeKnp          - make summand K_{\nu,p}
# CoHzero          - true if cohomology is zero
# CohDim
# allsums          - all q-sums
# printBlocks      - print complex as blocks
# printCohs        - print complex as cohomologies
# Jdiscr           - determinant of discrete jacobian
# deltaof          - compute the matrix of a polynomial
# Jtoric           - determinant of toric jacobian
# sparsematrixplot - plot tool
# visVects         - visualize det vects
# visVect          -
# map2axel         - output 3d data
# has_deter        - true if data is determinantal
# sigmaK           - sum l_i , s(k)<-s(i)
# testboxes        - run experiments

#########################################################################
#########################################################################

# Multihomogeneous Bezout Bound

mBezout:= proc(nis::vector, dis::vector, sis::vector)
	local grps,n;
	n:=convert(convert(nis,list),`+`);
	grps := vectdim(nis);
	RETURN(	convert([seq(dis[i]^nis[i],i=1..grps)],`*`)*
   		convert(convert(sis,list),`*`)*
		convert([seq(1/sis[i], i=1..n+1)],`+`)*factorial(n)/
		convert([seq(factorial(nis[i]),i=1..grps)],`*`)	);
end:

### union of all P_k sets
computeP:= proc(nis::vector, dis::vector, mis::vector)
	local grps, i, k, P;
	P:= {};
	grps := vectdim(nis);

	for k from 1 to grps do
		P:= P union
            {seq(floor(mis[k]/dis[k]+1)..floor((mis[k]+nis[k])/dis[k]))};
	od;
	RETURN(P);
end: #computeP

### true iff H^q(m-pd)=0
coHzero := proc(nis::vector,q::integer,mpd::vector,summs::array)
	local r, i,n, nonz, goodsum;
	n:=convert(convert(nis,list), `+`);
	r:=vectdim(nis):
	if r <> vectdim(mpd) then ERROR(`BAD DIMS`) fi;
	if not type(summs[0],set) or not type(summs[sum(nis[i],i=1..r)],set)
	 then ERROR(`BAD summs`) fi;  if q<0 or q>n then RETURN(true); fi;
	for nonz in summs[q] do
		goodsum:=false;
		for i from 1 to r do
			if member(i,nonz) then
				if mpd[i] >= -nis[i] then
					goodsum:=true; fi;
				else if mpd[i] < 0 then
					goodsum:=true;
				fi;
			fi;
		od:
		if not goodsum then RETURN(false) fi;
	od:
true;
end:#coHzero

### Cohomology  H^q(X,mpd)  dimension
coHdim := proc(nis::vector,q::integer,mpd::vector)
 local summs,r, i, nonz, goodsum, dimH,dimHs, zerotest;

 summs:=allsums(nis);
 r:=vectdim(nis): if r <> vectdim(mpd) then ERROR(`BAD DIMS`) fi;
 dimH:=0;
 for nonz in summs[q] do  # sum_summs
  goodsum:=false;dimHs:=1;
  i:=0;
  while dimHs>0 and i<r do
    i:=i+1;  # prod
    if member(i,nonz) then # H^ni (mpd_i)
      if mpd[i] >= -nis[i] then #print(i,nonz,mpd[i],nis[i]);
        goodsum:=true; dimHs:=0;
        else dimHs := dimHs * numbcomb(-mpd[i]-1,nis[i]); fi;
    else if mpd[i] < 0 then #print(i,mpd[i]);
        goodsum:=true; dimHs := 0;
        else dimHs := dimHs * numbcomb(mpd[i]+nis[i],nis[i]); fi;
    fi;
  od:#while (i)
  if not goodsum then
     if dimHs=0 then ERROR(`BAD computation`) fi;
  fi;
  dimH := dimH + dimHs;
 od:
 eval(dimH);
end:    # coHdim

### Dimension of Kv
dimKv:=proc(nis::vector,dis::vector,sis::vector,mis::vector,summs::array,Nu::integer)
 local n,grps, i,p, dimK;

 grps:=vectdim(nis):
 if grps <> vectdim(dis) then ERROR(`BAD DIMS`) fi;
 if grps <> vectdim(mis) then ERROR(`BAD DIMS for m-vector`) fi;
 unassign('i'); n:=sum(nis[i],i=1..grps);
 dimK := 0;
 for p from max(0,Nu) to min(n+1,Nu+n) do
     dimK := dimK + dimKvp(nis,dis,sis,mis,Nu,p);
 od:
 eval(dimK);
end:# dimKv

#### Dimension of Kvp
dimKvp:= proc(nis::vector, dis::vector, sis::vector,mis::vector, nu::integer,p::integer)
	local c,Kvp,n,grps,q,dim;
	grps:=vectdim(nis);
	n:=convert(convert(nis,list), `+`);	q:=p-nu;
	dim:=0;	for c in choose([seq(1..n+1)],p) do
		dim:=dim+ coHdim(nis,q, evalm( mis- convert([seq(sis[i],i=c)],`+`)*dis) );
	od;
	RETURN(dim);
end:# dimKvp

##### Critical degree  vector
critVect:= proc(nis::vector, dis::vector, sis::vector)
	local grps,n,r,s;
	s:=convert(convert(sis,list),`+`);
	grps := vectdim(nis);
	n:=convert(convert(nis,list),`+`);
	r:= s*dis - nis - vector(grps,1);
	evalm(r);
end:


### Bounds for determinantal degree vectors:
allDetVecs := proc(nis::vector,dis::vector,sis::vector)
 local	cand,
	tmp, misD, low,upp, grps,n,i, goodmis, summs, cnt,
	msmall,small,mbig,big;

 grps:=vectdim(nis):
 if grps <> vectdim(dis) then ERROR(`BAD DIMS`) fi;
 n:=sum(nis[i],i=1..grps);
 goodmis:={};
 summs:=allsums(nis):

 #compute the bounds
 low:=evalm(scalarmul(dis,-1));
 for i from 1 to grps do
   if nis[i]<dis[i] then low[i]:=-nis[i] fi;
 od;
 unassign('i');
 upp:=matadd(dis,vector(grps,1), sum(sis[i],i=1..n+1) ,-1);
 for i from 1 to grps do
   if nis[i]>dis[i] then upp[i]:=upp[i]+dis[i]-nis[i] fi;
 od;

 cand:=evalm(low);

 cnt:=1;
 unassign('i'): small:= sis[n]+sis[n+1]:
 unassign('i'); big:= sum(sis[i],i=1..n-1): unassign('i');

 while cand <> 0 do

  # necessary condition
  msmall := false; mbig:=false;
 for i from 1 to grps do
   if cand[i]<small*dis[i] then msmall:=true; break; fi;
 od:#i
 for i from 1 to grps do
   if cand[i] >= big*dis[i]-nis[i] then mbig:=true; break; fi;
 od:#i

  if msmall and mbig and is_determ(nis,dis,sis,eval(cand),summs) then
      misD:=[op(convert(cand,list)),dimKv(nis,dis,sis,cand,summs,0)];
      goodmis:=eval(goodmis) union {eval(misD)};
  fi;
  tmp := next_lex_rect(eval(cand),evalm(low),evalm(upp));
  cand := eval(tmp);
  cnt:=cnt+1;
 od;

 matadd(matadd(upp,low,1,-1),vector(grps,1));
 unassign('i');
 tmp :=product(%[i],i=1..grps);
 if cnt-1 <> tmp then ERROR(`bad count`,cnt-1,tmp) fi;

 sort(convert(goodmis,list), sort_dim);
end:	# allDetVecs

# find/return next vector (following mis) in lex-order
next_lex_rect:=proc(mis::vector,low::vector,upp::vector)
 local i, j, r;

 r :=vectdim(mis);if r <> vectdim(low) then ERROR(`BAD DIMS`) fi;
 if r <> vectdim(upp) then ERROR(`BAD DIMS`) fi;
 for i from r to 1 by -1 do if mis[i]<upp[i] then
   mis[i]:=mis[i]+1;
   for j from i+1 to r do mis[j]:=low[j]; od;

   RETURN(mis);
   fi;
 od;
 eval(0);
end:# next_lex_rect

detboxes := proc(nis::vector,dis::vector,sis::vector )
local p, r, n, B, box, k, left, right, ss;

r:= vectdim(dis);
n:= convert(convert(nis,list),`+`);

B := NULL;
for p in permute(r) do
	box:= NULL;
	for k to r do
		ss   := sigmaK(k,nis,p);
		left := dis[k]* convert( [seq(sis[i], i=n-ss+3..n+1)   ], `+`) - nis[k];
		right:= dis[k]* convert( [seq(sis[i], i=1..ss-nis[k]+2)], `+`) - 1; #print(p,[left,right]);
		if left<=right then box:= box, [left,right] ; else box:=NULL; break;  fi;
	od;
	if not box=NULL then B:= B, [ box ]; fi; #print(p);
od;
B;
end: #detboxes


has_deter := proc (nis::vector, dis::vector, sis::vector)
	local ss,k,r,n,p,d;
	n:= convert(convert(nis,list),`+`);
	r:= vectdim(nis);

	for k to r do
 	       if not dis[k]*convert([seq(sis[i],i=n-nis[k]+3..n+1)],`+`) < dis[k]*(sis[1]+sis[2])+nis[k] then
   	   	  RETURN(false); fi;
	od;

        for p in permute(r) do
	    d:=true;
     	    for k to r do
	       ss   := sigmaK(k,nis,p);
	       if dis[k]* convert( [seq(sis[i], i=n-ss+3..n+1)   ], `+`) - nis[k] >
	          dis[k]* convert( [seq(sis[i], i=1..ss-nis[k]+2)], `+`) - 1
	       then d:=false; fi;
	    od;
	if d then RETURN(true); fi;
	od;
false;
end: #has_deter

is_determ := proc(nis::vector,dis::vector,sis::vector,mis::vector,summs::array)
 local n,grps, Nu, i,p,c, P;
 grps:=vectdim(nis):
 if grps <> vectdim(dis) then ERROR(`BAD DIMS`) fi;
 if grps <> vectdim(mis) then ERROR(`BAD DIMS for m-vector`) fi;
 unassign('i');
 n:=sum(nis[i],i=1..grps);
 P:= computeP(nis,dis,mis);

 for Nu in {-1,2} do for p from max(0,Nu) to min(n+1,Nu+n) do
     for c in choose([seq(1..n+1)],p) do
           if not coHzero( nis , p-Nu ,
              evalm( mis- convert([seq(sis[i],i=c)],`+`)*dis) , summs)
		then RETURN(false) fi:
 od:od:od:

 evalb(true);
end:	# is_determ


### true iff l1<l2
sort_dim:=proc(l1,l2)
 local i;
 if l1[-1]<l2[-1] then RETURN(true)
 else RETURN(false)
 fi;
end:	# sort_dim


#all variables, var the vector of group names, for example bihomo: [x,y]
allvars:= proc(nis::vector,var)
local i,vars;

vars:= NULL;
for i to vectdim(nis) do
if not var[i]=0 then
vars:= vars, seq( var[i][k], k=1..nis[i] );
fi; od;
[vars];
end:

#Multidegree with respect to vector nis
mdegree:= proc(g, nis::vector, var )
local deg, i, vars;
vars:=allvars(nis, var);

deg:=NULL;
for i to nops(var) do
    deg:=deg, degree(g, [seq( var[i][k], k=1..nis[i] )] );
od;
vector([deg]);
end:

### The natural monomial basis of S(mdeg)
monbasis := proc(nis::vector,mdeg::vector,var)
  local p, g, gvar, vars;
  vars:=NULL;
  if not ( nops(var)= vectdim(nis)and nops(var)=vectdim(mdeg)) then
    ERROR(`Wrong dimensions`); fi;
  p:=1;
  for g from 1 to vectdim(nis) do
       if var[g]=0 then next fi;

       gvar:= seq( var[g][k], k=1..nis[g] );
      vars:= vars, gvar;
      if mdeg[g]>=0 then gvar:=gvar;
	p:= p*expand( (1+convert([gvar] , `+`))^mdeg[g] );
      else
        gvar:=gvar,gvar[1]; #degeneracy quick fix
        p:= p*expand( (1+convert( [seq(1/gvar[k], k=1..nis[g])] , `+`) )^(-mdeg[g]) );
      fi;
  od;
  lstmonof(p,[vars]);
end:

### Matrix of multihomogeneous multiplication map
multmap:= proc(f,sp1::vector,sp2::vector,nis::vector,var )
  local i,j,row,col, vars, mat;
  row:= [monbasis(nis,sp1,var)];
  col:= [monbasis(nis,sp2,var)];

  vars:=allvars(nis,var);

  mat:= NULL;
  for i to nops(row) do
    for j to nops(col) do
     mat := mat, coeffof(expand(col[j]/row[i]), f , vars ) ;
    od;
  od;
  matrix(nops(row),nops(col),[mat]);
end:


### Dimension of S(uis)
Sdim:= proc(nis::vector,uis::vector)
local k, grps, dim;
grps:= vectdim(nis); dim:=1;

for k in 1..grps do
    if uis[k]>=0 then dim:=dim*numbcomb(uis[k]+nis[k],nis[k]);
    else  dim:=dim*numbcomb(-uis[k]+nis[k],nis[k]); fi;
od;
dim;
end:

### Sylvester Matrix block K_{1,p}->K_{0,p-1}
Sylvmat:= proc(KK1 ,  KK0, f, nis::vector, dis::vector, sis::vector, mis::vector , var)
  local p,n,i, r,c, u ,v, mat, cols,rows, grps, k ;

  n:=convert(convert(nis,list),`+`);
  grps := vectdim(nis);
  p:= KK1[1];

  rows:=NULL;
  for r from 2 to nops(KK1) do
     cols:=NULL;
     u:= evalm(mis-convert([seq(sis[i],i=KK1[r][2])],`+` )*dis);
     for k to grps do if u[k]<0 then u[k]:=u[k]+nis[k]+1; fi;od;
     for c from 2 to nops(KK0) do
        v:= evalm( mis- convert([seq(sis[i],i=KK0[c][2])],`+`)*dis ) ;
        for k to grps do if v[k]<0 then v[k]:=v[k]+nis[k]+1; fi;od;
        if convert(KK0[c][2],set) subset convert(KK1[r][2],set) then
              i:=1; while i<p and KK1[r][2][i]=KK0[c][2][i] do i:=i+1; od; k:=KK1[r][2][i];
              mat:= evalm( ((-1)^(i+1))*multmap(f[k],u, v, nis,var ) );
        else
              mat:= matrix( Sdim(nis,u), Sdim(nis,v), 0);
        fi;
        cols:= cols, convert(mat, Matrix);
     od;
     rows:= rows,[cols] ;
  od;
Matrix([rows]);
end:

### Bezoutian block of S(sp1)->S(sp2)
BezoutianPoly:= proc(f,sp1::vector,sp2::vector,nis::vector,var, grp:={} )
  local i,j,row,col, nvar, nvars, ovar, ovars, Bpol, mat;

  if grp={} then grp:={seq(1..nops(var))} fi;

  nvar:=NULL; ovar:=NULL;
  for i to nops(var) do
         if member(i,grp ) then
	    nvar:= nvar,cat(_,var[i]);
	    ovar:= ovar, var[i];
	 else
	    nvar:= nvar, 0;
	    ovar:= ovar, 0;
	 fi;
  od;
  nvar:= [nvar]; ovar:= [ovar];

  nvars := allvars(nis,nvar);
  ovars := allvars(nis,ovar);

  Bpol:= Jdiscr(f,ovars );
  row:= [monbasis(nis,sp1, var)];
  col:= [monbasis(nis,sp2,nvar)];

  nvars:= [op(allvars(nis,var)),op(nvars)];
  mat:= NULL;
  for i to nops(row) do
    for j to nops(col) do
     mat := mat, coeffof(col[j]*row[i], Bpol , nvars );
    od;
  od;
  matrix(nops(row),nops(col),[mat]);
end:


### Bezoutian block K_{1,a}->K_{0,b}
Bezoutmat:= proc(KK1 ,  KK0, f, nis::vector, dis::vector, sis::vector, mis::vector , var)
  local p, pols, subsvar, n,grps, r,c, rows, cols, u,v, k, mat;

  n:=convert(convert(nis,list),`+`);
  grps := vectdim(nis);
  p:= KK0[1];

  rows:=NULL;
  for r from 2 to nops(KK1) do
     cols:=NULL;
     u:= evalm(mis-convert([seq(sis[i],i=KK1[r][2])],`+` )*dis);
     for k to grps do if u[k]<0 then u[k]:=-u[k]-nis[k]-1; fi;od;
     for c from 2 to nops(KK0) do
        v:= evalm( mis- convert([seq(sis[i],i=KK0[c][2])],`+`)*dis ) ;
        for k to grps do if v[k]<0 then v[k]:=-v[k]-nis[k]-1; fi;od;
	if convert(KK0[c][2],set) subset convert(KK1[r][2],set) then
	   pols:= convert(convert(KK1[r][2],set) minus convert(KK0[c][2],set), list);
           subsvar:= KK1[r][1] minus KK0[c][1];
           mat:= BezoutianPoly( [seq(f[k],k=pols)] ,u,v,nis,var, subsvar );
        else
	   mat:= matrix( Sdim(nis,u), Sdim(nis,v), 0);
	fi;
	cols:= cols, convert(mat, Matrix);
     od;
     rows:= rows,[cols] ;
  od;
Matrix([rows]);
end:

### The Matrix K_1 -> K_0
makeMatrix:= proc(nis::vector, dis::vector, sis::vector, mis::vector, sysp:=[], varp:=[] )
  local sys, var, n, KK, matr, row, col, rows, cols, d1, d2;

n:=convert(convert(nis,list), `+`);
if nargs<5 then
  var:= [seq( cat(x,i), i=1..vectdim(nis))];
  sys:= Makesystem(nis,dis,sis,[seq( cat(c,i-1), i=1..n+2)],var) ;
else
  var:=varp;
  sys:=Makesystem(nis,dis,sis,sysp,var) ;
fi;

  #compute the complex
  KK:=MakeComplex(nis,dis,sis,mis);

  #for now demand det/al complex
  if not nops(KK[1])=2 then ERROR(`Not Determinantal`) fi;

  #compute block dimensions
  d1:=0;
  for row from 2 to nops(KK[1][1]) do
	d1:= d1, CohDim(KK[1][1][row], nis, dis, sis, mis);
  od; d1:= [d1];
  d2:=0;
  for col from 2 to nops(KK[1][2]) do
	d2:= d2, CohDim(KK[1][2][col], nis, dis, sis, mis);
  od; d2:= [d2];

  rows:= NULL;
  for row from 2 to nops(KK[1][1]) do
     cols:= NULL;
     for col from 2 to nops(KK[1][2]) do
	    if KK[1][1][row][1]-1 < KK[1][2][col][1] then matr:= matrix( d1[row] , d2[col] , 0 ); #Zero map
       else if KK[1][1][row][1]-1 = KK[1][2][col][1] then
             matr:=   Sylvmat(KK[1][1][row],KK[1][2][col], sys, nis,dis,sis,mis,var) ;
       else  matr:= Bezoutmat(KK[1][1][row],KK[1][2][col], sys, nis,dis,sis,mis,var) ;
       fi;fi;
       cols:= cols, convert(matr,Matrix);
     od;
     rows:= rows, [cols];
  od;
  convert(Matrix([rows]),matrix);
end:


### Make polynomial with m-degree scalar*dis
Makepoly:= proc(nis::vector,dis::vector,scalar, c, var )
  local p, i, s;

  i:=0;p:=0;
  s:=[monbasis(nis, evalm(scalar*dis),var)];
  for i from 1 to nops(s) do
     p:= p+ c[i-1]*s[i];
  od;
  p;
end:

### Make extreme polynomial with degree scalar*dis
MakeExtremepoly:= proc(nis::vector,dis::vector,scalar, c, var )
  local v, p, i,j,dd, vars;

  vars:= allvars(nis,var);
  dd:= evalm(scalar*dis);
  p:=1;v:=1;
  for i to vectdim(nis) do
     for j to nis[i] do
       p:= p*(1+ vars[v]^dd[i] );
       v:= v+1;
  od;od;

  v:= [lstmonof(evalm(p),vars)];
  p:=0;
  for i from 1 to nops(v) do
     p:= p+ c[i-1]*v[i];
  od;
  p;
end:

### Make multihomogeneous system (nis;dis;sis)
Makesystem:= proc(nis::vector,dis::vector,sis::vector, coef, var )
  local i,n,lst;

  n:=convert(convert(nis,list),`+`);
  lst:=NULL;
  for i from 1 to n+1 do
      lst:= lst, Makepoly(nis,dis,sis[i], coef[i] ,var);
  od;
  [lst];
end:

###List the monomials of a polynomial p in the variables var:
lstmonof := proc(p,var)
local lm,r,c;
  lm := NULL;
  r := sort(p,var,tdeg);
  while r <> 0 do
     c := tcoeff(r,var,'m');
     lm := lm, m;
     r := expand(r-c*m);
  od;
  lm;
end:

###return the coeff. of p in variables var of the monomial m:
coeffof := proc(m,p,var)
local lm,lc,k;
  lc := [coeffs(p,var,'lm')];
  lm := [lm];
  if member(m,lm,'k') then lc[k] else 0 fi;
end:

### Compute the Wayman Complex
MakeComplex:=proc(nis::vector, dis::vector, sis::vector, mis::vector)
	local i,summs, Kv, c,n, K;

	summs:=allsums(nis);
	n:=convert(convert(nis,list),`+`);
	K:= NULL;
	for i in seq(-e, e=-n-1..n) do
		Kv := MakeKv(nis,dis,sis,mis,i,summs);
		if not Kv=[] then
			K:= K, Kv;
		fi;
	od;
[[K], evalm(nis), evalm(dis), evalm(sis), evalm(mis)];
end:

### Make term Kv###
MakeKv:= proc(nis::vector, dis::vector, sis::vector, mis::vector, nu::integer,summs::array)
	local b, Kv, n, p, Kvp;
	n:=convert(convert(nis,list),`+`);
	Kv:=nu;	for p in seq(0..n+1) do
		Kvp := MakeKvp(nis,dis,sis,mis,nu,p, summs );
		if not Kvp=[] then
			Kv:= Kv, Kvp;
		fi;	od;
	if Kv=nu then RETURN([]); else RETURN([Kv]); fi;end:

###### Make term Kv,p
MakeKvp:= proc(nis::vector, dis::vector, sis::vector,mis::vector, nu::integer,p::integer,summs::array)
	local c,Kvp,n,q,s;
	n:=convert(convert(nis,list), `+`);
	q:=p-nu;
	if q<0 or q>n then RETURN([]); fi;
	Kvp:= p ; 	for c in choose([seq(1..n+1)],p) do # improve ?
		for s in summs[q] do
		if not CoHzero(s,q, nis, evalm( mis- convert([seq(sis[i],i=c)],`+`)*dis)) then #improve ?
	                Kvp:=Kvp, [ s , c];
			break;
                fi;od;od;
	if Kvp=p then RETURN([]); else RETURN([Kvp]); fi;end:

### true iff H^q(mpd)=0
CoHzero := proc(grp::set, q::integer, nis::vector, mpd::vector)
  local r, i,n, nonz, goodsum;
  n:=convert(convert(nis,list),`+`);
  r:= vectdim(nis);
  for i from 1 to r do    if member(i,grp) then
         if mpd[i] >= -nis[i] then RETURN(true); fi;
    else
	if mpd[i] < 0        then RETURN(true); fi;    fi;  od;
  RETURN(false);
end:

### Dimension of Kv,p
CohDim:= proc(KKvp, nis::vector, dis::vector, sis::vector, mis::vector)
  local i, h, dim, prod, u ;

  dim:= 0;
  for h from 2 to nops(KKvp) do
	u:= evalm(mis-convert([seq(sis[i],i=KKvp[h][2])],`+` )*dis);
	prod:=1; #compute dim H^q
	for i from 1 to vectdim(nis) do
	   if u[i]<0 then
		prod:= prod* numbcomb( -u[i]-1    ,nis[i] );
	   else prod:= prod* numbcomb( u[i]+nis[i],nis[i] );
	   fi;
	od;
	dim:= dim + prod ;
  od;
  dim;
end:

### all q-sums by dynamic programming
allsums := proc(nis::vector)
	local n, grps, i, q, Subs, ss, summs;
	grps:=vectdim(nis): unassign('i');
	n:=convert(convert(nis,list), `+`);
	summs := array(0..n):
	for i from 0 to n do
	       summs[i]:={}
	od:
	unassign('i');
        Subs := subsets({seq(i,i=1..grps)}):
	       while not Subs[finished] do
	          ss:=Subs[nextvalue]();
	          q := 0:
                  for i in ss do q := q + nis[i] od;
                  summs[q] := summs[q] union {ss};
                  od:
eval(summs);
end:


### Print Kv,p 's
printBlocks:= proc( KK )
  local v, p, sum;

  for v to nops(KK[1]) do
     sum:=0;
     for p from 2 to nops(KK[1][v]) do
	sum:= sum+  K[ KK[1][v][1] , KK[1][v][p][1] ] ;
     od;
  print(K[ KK[1][v][1] ]= sum  );
  od;
end:

### Print H^q 's
printCohs:= proc( KK )
  local u, r, q, v, p, h, sum;
  r:= vectdim(evalm(KK[2]));

  for v to nops(KK[1]) do
     sum:=1;
     for p from 2 to nops(KK[1][v]) do
	 q:= KK[1][v][p][1] - KK[1][v][1];
	 for h from 2 to nops(KK[1][v][p]) do
	   u:= evalm(KK[5]-convert([seq( KK[4][i],i=KK[1][v][p][h][2])],`+` )*KK[3]);
	   sum:= sum* H[q]( seq(u[i], i=1..r) );
	 od;
     od;
  print( K[KK[1][v][1]]= sum );
  od;
end:

### OLD: true iff H^q(mpd)=0
is_coHzero := proc(nis::vector,q::integer,mpd::vector,summs::array)
    local r, i,n, nonz, goodsum;

    n:=convert(convert(nis,list), `+`);
    r:=vectdim(nis): if r <> vectdim(mpd) then ERROR(`BAD DIMS`) fi;
    if not type(summs[0],set) or
       not type(summs[sum(nis[i],i=1..r)],set)
    then ERROR(`BAD summs`)
    fi;

    if q<0 or q>n then RETURN(true); fi;
    for nonz in summs[q] do
      goodsum:=false;
      for i from 1 to r do
      if member(i,nonz) then
      if mpd[i] >= -nis[i] then
         goodsum:=true; fi;
      else if mpd[i] < 0 then
         goodsum:=true; fi;    fi;  od:
      if not goodsum then RETURN(false) fi;
    od:
true;
end:

# Determinant of the discrete Jacobian
Jdiscr := proc(lp::list, vx)
local N,i,j,k,vxi,mtr,s;

  N := nops(lp);
  if nops(vx) <> N-1 then print(nops(lp), vx );
    ERROR(`the number of polynomials must be the number of variables +1 `);
  fi;
  mtr := array(1..N,1..N);
  for i from 1 to N do mtr[i,1] := lp[i] od;

  for j from 2 to N do
     for i from 1 to N do
       mtr[i,j]:= subs(vx[j-1]= cat(_,vx[j-1]) ,mtr[i,j-1]); ## vy[j-1]
     od;
  od;

  for j from N to 2 by -1 do
     for i from 1 to N do
       mtr[i,j] := quo(mtr[i,j]-mtr[i,j-1], cat(_,vx[j-1]) - vx[j-1], cat(_,vx[j-1]));
     od;
  od;

  det(convert(mtr,matrix));
end:


#List operations
append := (x, L) -> [x,op(L)]:
merge  := (L, Q) -> [op(L),op(Q)]:

deltaof := proc(Th,lx ::list,lxi::list)
local var,lm,i,j,ll,c,vx,vxi,p;

  vx  := [lstmonof( Th, lx  )] ; print(nops(%),"--",% );
  vxi := [lstmonof( Th, lxi )] ; print(nops(%),"--",% );

  ll := NULL;
  p:= expand(Th);
  for i to nops(vx) do
    c := coeffof(vx[i],p,lx);
    for j to nops(vxi) do
     ll := ll, coeffof(vxi[j],c,lxi);
    od;
  od;
  matrix(nops(vx),nops(vxi),[ll]);
end:

#Old sylvester matrix
OSylvmat:= proc(f, p, nis::vector, dis::vector, sis::vector, mis::vector , var)
  local n,i, r,c, u ,v, mat, cols,rows, grps, k , summs;

  summs:=allsums(nis);

  n:=convert(convert(nis,list),`+`);
  grps := vectdim(nis);

  rows:= NULL;
  for r in choose(n+1,p) do
     cols:= NULL;
     u:= evalm( mis- convert([seq(sis[k],k=r)],`+`)*dis ) ;
     if is_coHzero(nis,p-1,u,summs) then next; fi;
     for k from 1 to grps do if u[k]<0 then u[k]:=u[k]+nis[k]+1; fi;od;
     for c in choose(n+1,p-1) do
	  v:= evalm( mis- convert([seq(sis[k],k=c)],`+`)*dis ) ;
	  if is_coHzero(nis,p-1,v,summs) then next; fi;
          if convert(c,set) subset convert(r,set) then
              i:=1; while i<=nops(c) and r[i]=c[i] do i:=i+1; od;
                mat:= evalm( ((-1)^(i+1))*multmap(f[r[i]],u, evalm(u+sis[r[i]]*dis), nis,var ) );
          else
                for k from 1 to grps do if v[k]<0 then v[k]:=v[k]+nis[k]+1; fi;od;
                mat:= matrix( Sdim(nis,u), Sdim(nis,v), 0);
          fi;
          cols:= cols, convert(mat, Matrix);
     od;
         rows:= rows,[cols] ;
  od;
convert(Matrix([rows]),matrix);
end:

# OLD Partial Bezoutian matrix (sp1)->(sp2)
OBezoutmat:= proc(fset,sp1::vector,sp2::vector,nis::vector,var, grp:={} )
  local i,j,row,col, nvar, nvars, ovar, ovars, Bpol, mat;

  if grp={} then grp:={$1..nops(var)} fi;

  nvar:=NULL; ovar:=NULL;
  for i to nops(var) do
         if member(i,grp) then
	    nvar:= nvar,cat(_,var[i]);
	    ovar:= ovar, var[i];
	 else
	    nvar:= nvar, 0;
	    ovar:= ovar, 0;
	 fi;
  od;
  nvar:= [nvar]; ovar:= [ovar];

  nvars := allvars(nis,nvar);
  ovars := allvars(nis,ovar);

  Bpol:= Jdiscr(fset,ovars );
  if convert(sp1,list)=convert(mdegree(Bpol, nis, var),list) then		print(case1);
     row:= [monbasis(nis,sp1, var)]; print(%);
     col:= [monbasis(nis,sp2,nvar)]; print(%);
  else if convert(sp2,list)=convert(mdegree(Bpol, nis, var),list) then		print(case2);
     row:= [monbasis(nis,sp1, nvar)]; print(%);
     col:= [monbasis(nis,sp2,  var)]; print(%);
  else				## Dirty
     Bpol:= Jtoric(fset,ovars );
     nvars:=[];
  if convert(sp1,list)=[0,0] then						print(case3);
       row:= [1];
       col:= [ lstmonof(Bpol,ovars) ]; print(%);
  else if convert(sp2,list)=[0,0] then						print(case4);
       row:= [ lstmonof(Bpol,ovars) ]; print(%);
       col:= [1];
  else print(__UNDETERMINED__);
  fi;fi;fi;fi;

  nvars:= [op(allvars(nis,var)),op(nvars)];
  mat:= NULL;
  for i to nops(row) do
    for j to nops(col) do
     mat := mat, coeffof(col[j]*row[i], Bpol , nvars );
    od;
  od;
  matrix(nops(row),nops(col),[mat]);
end:

# Determinant of the toric Jacobian
Jtoric := proc(lp::list, vx)
local N,i,j,k,vxi,mtr,s;

  N := nops(lp);
  if nops(vx) <> N-1 then print(nops(lp), vx);
    ERROR(`the number of polynomials must be the number of variables +1 `);
  fi;
  mtr := array(1..N,1..N);
  for i from 1 to N do mtr[i,1] := lp[i] od;

  for j from 2 to N do
     for i from 1 to N do
       mtr[i,j]:= diff( lp[i], vx[j-1]);   subs(vx[j-1]= cat(_,vx[j-1]) ,mtr[i,j-1]); ## vy[j-1]
     od;
  od;

  det(convert(mtr,matrix));
end: #Jtoric

# MATRIX OF DISTRIBUTION OF DETERMINANTAL m-VECTORS
allDetVecsMatrix := proc(nis::vector,dis::vector,sis::vector)
 local	OPT, MM, cand,	# candidate m-vector
	tmp, misD, low,upp, grps,n,i, goodmis, summs, cnt,
	msmall,small,mbig,big;

 grps:=vectdim(nis):
 if grps <> vectdim(dis) then ERROR(`BAD DIMS`) fi;
 if grps <> 2 then ERROR(`MATRIX WORKS ONLY FOR 2 GRPS`) fi;
 n:=sum(nis[i],i=1..grps);
 goodmis:={};
 summs:=allsums(nis):

 #compute the bounds
 low:=evalm(scalarmul(dis,-1));
 for i from 1 to grps do
   if nis[i]<dis[i] then low[i]:=-nis[i] fi;
 od;
 unassign('i');
 upp:=matadd(dis,vector(grps,1), sum(sis[i],i=1..n+1) ,-1);
 for i from 1 to grps do
   if nis[i]>dis[i] then upp[i]:=upp[i]+dis[i]-nis[i] fi;
 od;

MM:= Matrix(upp[1]-low[1]+1, upp[2]-low[2]+1); #print( low[1]..upp[1], low[2]..upp[2]  );

 cand:=evalm(low); # (watch evalM!)
 #print(`From`,low,`to`,upp,`start at`,cand);
 cnt:=1;
 unassign('i'): small:= sis[1]+sis[2]:
 unassign('i'); big:= sum(sis[i],i=1..n-1): unassign('i');

 while cand <> 0 do

  # necessary condition for mis to be determinantal
  msmall := false; mbig:=false;
 for i from 1 to grps do
   if cand[i]<small*dis[i] then msmall:=true; break; fi;
 od:#i
 for i from 1 to grps do
   if cand[i] >= big*dis[i]-nis[i] then mbig:=true; break; fi;
 od:#i

  if msmall and mbig and is_determ(nis,dis,sis,eval(cand),summs) then # printf("Determinantal %a\n",evalm(cand));
      misD:=[op(convert(cand,list)),dimKv(nis,dis,sis,cand,summs,0)];
      goodmis:=eval(goodmis) union {eval(misD)};
      MM[ upp[1]+1 -cand[1],cand[2]-low[2]+1]:= 1;  upp[1]-low[1]+1 # MM[cand[1]-low[1]+1,cand[2]-low[2]+1]:= 1;
  fi;
  tmp := next_lex_rect(eval(cand),evalm(low),evalm(upp));#print(%);
  cand := eval(tmp);
  cnt:=cnt+1;
 od;

 matadd(matadd(upp,low,1,-1),vector(grps,1));
 unassign('i');
 tmp :=product(%[i],i=1..grps);
 if cnt-1 <> tmp then ERROR(`bad count`,cnt-1,tmp) fi;

 tmp:= sort(convert(goodmis,list), sort_dim);
 if tmp=[] then RETURN( [transpose(convert((MM),array)), low[1], low[2],[] ] ); else cnt:= tmp[1][3]; fi;

 OPT := NULL;

 for i in tmp do
	if i[3]=cnt then OPT:= OPT, [ i[1]-low[1]+1,i[2]-low[2]+1];
 fi;
 od;

 [transpose(convert((MM),array)), low[1], low[2],[OPT] ]
end:	# allDetVecsMatrix

sparsematrixplot := proc (A)
local nn, mm, i, j, m, n, B, xticks, yticks, nq, mq, nticlabels, mticlabels,
      plot1, plot2, gt, indexproc;
 if nargs = 0 then ERROR(`wrong number of arguments`) fi;
 if not type(A[1],matrix) then ERROR(`matrix expected`) fi;
 m := linalg['rowdim'](A[1]);
 n := linalg['coldim'](A[1]);

 indexproc := proc(x,i,j,m)
     if x <> 0 then RETURN([j,m-i+1]) else RETURN(NULL) fi;
 end;

 B :=[seq(seq(indexproc(A[1][i,j],i,j,m),j=1..n),i=1..m)];

 xticks := ('xtickmarks') = 0;
 yticks := ('ytickmarks') = 0;
 nn := min(n,5);  nq:=iquo(n,nn);  mm := min(m,5); mq:=iquo(m,mm);
 gt := proc(x,n) if x > n then NULL else x fi end;
 nticlabels := {1,seq(gt(i*nq,n),i=1..n),n};
 nticlabels := [seq(i=A[2]-1 + i,i=nticlabels)];
 mticlabels := {1,seq(gt(i*mq,m),i=1..m),m};
 mticlabels := [seq(i=A[3]-1 +i,i=mticlabels)];
 plot1 := plots['pointplot'](B,xticks,labels=[group1,group2],yticks,args[2 .. nargs],
                ('style') = ('POINT'),view=[0..n,0..m], color=blue );
 plot1 := [op(plot1)];
 plot1 := remove(has,plot1,{AXESSTYLE,AXESTICKS});
 plot1 := PLOT(op(plot1),AXESTICKS(nticlabels,mticlabels));
 plot2:= plots['pointplot'](A[4], color=red, symbol=box,symbolsize=15);

 plot2:= plots[display]([plot1,plot2]);
 #plot1;
end: #sparsematrixplot

visVects := proc ( l::vector )
#A:= allDetVecsMatrix(nis,dis,sis);
#sparsematrixplot(A);
local A, flag, dets, bb, i,j, d, d1, s, s0,  s1, n, r, SS, DD, ss ;

r:= vectdim(l);
n:=convert(convert(l,list),`+`);

s0 := vector([seq(1, i=1..n+1)]);
s1 := vector([seq(3, i=1..n+1)]);
d  := vector([seq(1, i=1..r)]);
d1 := vector([seq(2, i=1..r)]);
DD := cartprod( [seq( [seq( i, i=d[j]..d1[j])], j=1..r  )] );

while not DD[finished] do   ################# Loop degrees

    # Next Degree
       d:= vector( DD[nextvalue]() );
    # end Next Degree

    SS:= combinat:-cartprod( [seq( [seq( i, i=s0[j]..s1[j])], j=1..n+1  )] )  ;

    while not SS[finished] do   ################# Loop scalars

    # Next scalar vector
        s:= vector( SS[nextvalue]() );
        if not 1 = igcd(seq(s[i],i=1..n+1)) then next; fi;
        flag:= false;
        for i from 1 to n do
            if s[i]>s[i+1] then flag:=true; fi;
        od;
        if flag=true then next; fi;
    # end Next Scalar
	if not has_deter(l,d,s) then next; fi;
	print(`Degree:`, d, `Scalars:`, s, `Critical:`,  critVect(l,d,s) );

	A:= allDetVecsMatrix(l,d,s);
	print(sparsematrixplot(A) );

    od;
od;
l;
end:

visVect := proc ( nis::vector,dis::vector,sis::vector )
local A;
sparsematrixplot( allDetVecsMatrix(nis,dis,sis) );
end:#visVect


map2axel := proc (nis::vector, dis::vector, sis::vector)
local dv, axl, i ;

if not vectdim(nis) = 3 then ERROR("3D data only"); fi;

dv:= allDetVecs(nis,dis,sis);
axl:= "/home/angelo/out.axl" ;

fopen(axl,WRITE,TEXT);
fprintf(axl,"<axl>\n 	<pointset size=\"%d\">\n", nops(dv) );

for i in dv do
	fprintf(axl, "%d  %d  %d\n", i[1],i[2],i[3] ) ;
od;

fprintf(axl,"</pointset> \n </axl>" );
fclose(axl);
dv;
end:#map2axl

sigmaK:= proc(k::integer, nis::vector, p)
local i, sum;
	sum:= 0;
	for i to vectdim(nis) do
		if p[i]<= p[k] then
			sum:= sum + nis[i];
		fi;
	od;
sum;
end: #sigmaK

testBoxes := proc ( l::vector )
#A:= allDetVecsMatrix(nis,dis,sis);
#sparsematrixplot(A);
local A, flag, dets, bb, i,j, d, d1, s, s0,  s1, n, r,k, SS, DD, ss ;

r:= vectdim(l);
n:=convert(convert(l,list),`+`);

s0 := vector([seq(1, i=1..n+1)]);
s1 := vector([seq(3, i=1..n+1)]);
d  := vector([seq(1, i=1..r)]);
d1 := vector([seq(2, i=1..r)]);
DD := cartprod( [seq( [seq( i, i=d[j]..d1[j])], j=1..r  )] );

while not DD[finished] do   ################# Loop degrees

    # Next Degree
       d:= vector( DD[nextvalue]() );
    # end Next Degree

    SS:= combinat:-cartprod( [seq( [seq( i, i=s0[j]..s1[j])], j=1..n+1  )] )  ;

    while not SS[finished] do   ################# Loop scalars

    # Next scalar vector
        s:= vector( SS[nextvalue]() );
        if not 1 = igcd(seq(s[i],i=1..n+1)) then next; fi;
        flag:= false;
        for i from 1 to n do
            if s[i]>s[i+1] then flag:=true; fi;
        od;
        if flag=true then next; fi;
    # end Next Scalar

        if not has_deter(l,d,s) then next; fi;
	print(`Degree:`, d, `Scalars:`, s, `Critical:`,  critVect(l,d,s) );
    od;
od;
l;
end:


printf("Scaled Multihomogeneous Resultants");