newPackage(
	"EDPolytope",
	Version => "2.2", 
    	Date => "April 20, 2016",
    	Authors => {{Name => "Martin Helmer", 
		  Email => "martin.helmer@berkeley.edu", 
		  HomePage => "https://math.berkeley.edu/~mhelmer/"}},
    	Headline => "Computes ED degrees of toric varieites X_A from the polytope conv(A)",
    	DebuggingMode => false,
	Reload => true
    	);

needsPackage "Polyhedra"
export{"PolarDegrees",
       "EDdeg",
       "CMVolumes",
       "CMClass",
       "normalizedVolume",
       "hyperSimplexVertices",
       "Output",
       "TextOutput"
    }

hyperSimplexVertices=method(TypicalValue=>Matrix)
hyperSimplexVertices (ZZ,ZZ):=(d,k)->(
    hsl:=join(toList(k:1),toList((d-k):0));
    return transpose matrix unique permutations(hsl);
    );
CMClass=method(TypicalValue=>RingElement,Options => {TextOutput=>"Quiet",Output=>RingElement})
CMClass (Matrix,QuotientRing):=opts->(A,Chring)->(
    if opts.Output===HashTable then(
	return PolarDegrees(A,Chring,Output=>HashTable,TextOutput=>opts.TextOutput);
	);
    return (PolarDegrees(A,Chring,Output=>HashTable,TextOutput=>opts.TextOutput))#"CM class";
    );
CMClass (Matrix):=opts->(A)->(
        if opts.Output===HashTable then(
	    return PolarDegrees(A,Output=>HashTable,TextOutput=>opts.TextOutput);
	);
    return (PolarDegrees(A,Output=>HashTable,TextOutput=>opts.TextOutput))#"CM class"
    );

PolarDegrees=method(TypicalValue=>List,Options => {TextOutput=>"Quiet",Output=>List})
PolarDegrees (Matrix):=opts->(A)->(
    n:=numColumns(A);
    h:=symbol h;
    Chring:=ZZ[h]/(h^n);
    return PolarDegrees(A,Chring,Output=>opts.Output,TextOutput=>opts.TextOutput);
    );
PolarDegrees (Matrix,QuotientRing):=opts->(A,Chring)->(
    EV:=CMVolumes(A,TextOutput=>opts.TextOutput);
    m:=#EV-1;
    delta:={};
    temp:=0;
    for i from 0 to m do(
	temp=(-1)^(m-i)*sum(i..m,j->(-1)^(j-i)*binomial(j+1,i+1)*(EV_(m-j)));
	delta=append(delta, temp);
	);
    EVList:=reverse EV;
    A=transpose matrix unique entries transpose A;
    n:=numColumns(A);
    h:=Chring_0;
    CMClass:=sum(#(EVList),i->EVList_i*h^(n-1-i));
    FirstNonZero:=position (delta,i->(not i==0) );
    if opts.Output===List then(
	<<"The toric variety has degree = "<<last(EVList)<<endl;
	<<"The dual variety has degree = "<<delta_FirstNonZero<<endl;
    	<<"Chern-Mather Volumes: (V_0,..,V_(d-1)) = "<<EVList<<endl;
    	<<"Polar Degrees: "<<delta<<endl;
    	<<"ED Degree = "<<sum(delta)<<endl;
    	<<endl;
    	<<"Chern-Mather Class: "<<CMClass<<endl;
	);
    if opts.Output===HashTable then(
	delHash:=hashTable{"degree"=>last(EVList),"dual degree"=>first(delta),"hyperplane class"=>Chring_0,"Chow ring"=>Chring,"ED"=>sum(delta),"polar degrees"=>delta,"CM class"=>CMClass,"CM volumes"=>EVList};
	return delHash;
	);
    return delta;
    );
EDdeg=method(TypicalValue=>ZZ,Options => {TextOutput=>"Quiet",Output=>ZZ})
EDdeg (Matrix):=opts->(A)->(
    if opts.Output===HashTable then(
	return PolarDegrees(A,Output=>opts.Output,TextOutput=>opts.TextOutput);
	);
    return sum PolarDegrees(A,TextOutput=>opts.TextOutput);    
    );
CMVolumes=method(TypicalValue=>List,Options => {TextOutput=>"Quiet"})
CMVolumes (Matrix):=opts->(A)->(
    AinNum:=numColumns(A);
    A=transpose matrix unique entries transpose A;
    if opts.TextOutput!="Quiet" and AinNum!=numColumns(A) then (
	print "Warning: Input matrix contains duplicate columns, duplicate columns have been removed";
	);
    A=transpose hermite transpose A;
    if minors(numRows(A),A)!=ideal(1) then (
	if opts.TextOutput!="Quiet" then print "Input matrix does not generate the integer lattice, attempting to build a new matrix which generates the lattice and defines the same toric ideal";
	betterA:= transpose gens kernel transpose gens kernel A;
	if minors(numRows(betterA),betterA)==ideal(1) then(
	    if opts.TextOutput!="Quiet" then print "New matrix generated";
	    A=transpose hermite transpose betterA;
	    )
	else(
	    print "Input matrix does not generate the integer lattice, attempting to build a new matrix which generates the lattice and defines the same toric ideal";
	    error"Matrix generation failed, please enter a matrix whose columns span the full integer lattice";
	    return 0; 
	    );
	);
    if rank(A)!=min(numRows(A),numColumns(A)) then (
	error "Input matrix expected to have maximal rank";
	return 0;
	);
    if rank(A)<rank(A||matrix{toList(numColumns(A):1)}) then(
	if opts.TextOutput!="Quiet" then print "Adding a row of ones";
	A=(A||matrix{toList(numColumns(A):1)});
	)
    else(
	if not any(entries A,l->l==toList(numColumns(A):1)) then(
	    if opts.TextOutput!="Quiet" then print "Ones vector (1,1,...,1) in row span. Replacing last row with row of ones";
	    A=(A^{0..(numRows(A)-2)}||matrix{toList(numColumns(A):1)});
	    )
	else(
	    ne:=position(entries A,l->l==toList(numColumns(A):1));
	    A=submatrix'(A,{ne},)||A^{ne};
	    );
	);

    if opts.TextOutput!="Quiet" and ((numRows(A)+numColumns(A))<35) then <<"Proceeding with A-matrix="<<A<<endl;
    P:=convexHull(A);
    m:=dim(P);
    FacePoset:={};
    ifaces:={};
    tempFace:={};
    curFace:=0;
    curVol:=0;
    Acols:=entries transpose(A);
    notFcols:=0;
    interiorCols:={};
    faceCounts:={};
    for i from 0 to m do(
	ifaces={};
	tempFace=faces(i,P);
	faceCounts=append(faceCounts, #tempFace);
	for f in tempFace do(
	    curFace=new MutableHashTable;
	    curFace#"dim"=m-i;
	    curFace#"verticesMat"=lift(vertices(f),ZZ);
	    curFace#"verticesList"=entries transpose curFace#"verticesMat";
	    if i==m then(
		  curVol=1;
		 )
	     else(
		 curVol=((m-i)!)*volume(f);
		 );
	    curFace#"volume"=curVol;
	    notFcols=toList(set(Acols)-set(curFace#"verticesList"));
	    interiorCols={};
	    for c in notFcols do(
		if contains(f,transpose matrix{c}) then(
		    interiorCols=append(interiorCols,c);
		    );
		);
	    curFace#"interiorCol"=interiorCols;
	    curFace#"Aalpha"=join(interiorCols,curFace#"verticesList");
	    ifaces=append(ifaces,curFace);	    
	    );
	FacePoset=append(FacePoset,ifaces);
	);
    V:={};
    tempEL:=0;
    tempVL:=0;
    tempAlpha:=0;
    tempBeta:=0;
    tempsum:=0;
    volBeta:=0;
    tempEV:=0;
    ((FacePoset_0)_0)#"EulerOb"=1;
    V=append(V,((FacePoset_0)_0)#"EulerOb"*((FacePoset_0)_0)#"volume");
    for i from 1 to m do(
	 tempBeta=FacePoset_i;
	 tempEL={};
	 tempVL={};
	 tempEV=0;
         for j from 0 to #(tempBeta)-1 do(
	     tempsum=0;
		 for l from 0 to i-1 do(
		     tempAlpha=FacePoset_l;
		     for b from 0 to (#(tempAlpha)-1) do(
			 if isSubset((tempBeta_j)#"verticesList",(tempAlpha_b)#"verticesList") then(
			     tempsum=tempsum+(-1)^(i-l-1)*mu(A,tempAlpha_b,tempBeta_j)*((tempAlpha_b)#"EulerOb");
			     );
			 );
		     );
		 ((FacePoset_i)_j)#"EulerOb"=tempsum;
		 tempEL=append(tempEL,tempsum);
		 tempEV=tempEV+tempsum*((tempBeta_j)#"volume");
	     );
	 V=append(V,tempEV);
	);
    return V;
    );
normalizedVolume=method(TypicalValue=>ZZ)
normalizedVolume (Polyhedron,ZZ):=(P,n)->(
    if n>dim(P) then return 0;
    if n==0 then (
	    return 1;
	    ) 
	else return ((dim(P))!)*volume(P);
        );    
normalizedVolume (Polyhedron):=(P)->(
    if dim(P)==0 then return 1 else return ((dim(P))!)*volume(P);
        );     
---------------------------------------------------
--Internal Functions
--      
--mu=method(TypicalValue=>ZZ)
--------------------------------------------------
mu =(A,alpha,beta)->(
    dBeta:=beta#"dim";
    dAlpha:=alpha#"dim";
    r:=dAlpha-dBeta;
    d:=numRows(A);
    n:=numColumns(A);
    vbeta:=beta#"verticesMat";
    mbeta:=beta#"verticesList";
    valpha:=alpha#"verticesMat";
    malpha:=alpha#"verticesList";
    Atemp:=0;
    Asort3:=beta#"Aalpha";
    Asort2:=toList(set(alpha#"Aalpha")-set(beta#"Aalpha"));
    Asort1:=toList(set(entries transpose(A))-set(alpha#"Aalpha")); 
    M:=transpose(matrix(join(Asort1,Asort2,Asort3)));
    if numColumns(A)!=numColumns(M) then (
	print "There seems to be an error somewhere";
	);
    W:=(M);
    Anew:=transpose hermite(transpose(W));
    cs:=n-(#Asort2+#Asort3);
    ce:=n-#Asort3-1;
    rowInds:={};
    inc:=0;
    dind:=d-1;
    for w from 0 to dind do(
	if flatten(entries((Anew^{dind-w})_(toList((ce+1)..(n-1)))))==toList((#Asort3):0) then(
	    inc=inc+1;
	    rowInds=append(rowInds,dind-w);
	    );
	if inc==r then break;    
	);
    C:=(Anew_{cs..ce})^rowInds;
    if rank(C)==0 then return 1;
    big:=convexHull((transpose(matrix{toList(numRows(C):0)})|C));
    C1:=transpose matrix delete(toList(numRows(C):0), entries transpose C);
    little:=convexHull(C1);
    vol1:=normalizedVolume(big,r);
    vol2:=normalizedVolume(little,r);
    vol:=(vol1-vol2);
    if numRows(C)!=r then(
	print "Something may be wrong...we seem to have picked a C matrix with the wrong number of rows";
    );
    return vol;
    );

TEST ///
{*  
    restart
    needsPackage "EDPolytope"
*}  

    A=transpose(matrix{{0,0,1},{0,7,1},{3,0,1},{5,0,1}});
    time pd=EDdeg(A);
    time pd=PolarDegrees(A,TextOutput=>"All");
    A=hyperSimplexVertices(4,2);
    time pd=PolarDegrees(A);
    assert(pd=={4,12,8,4});
    time AHash=PolarDegrees(A,Output=>HashTable);
    CMVolumes(A);
    --peek AHash;
    assert(AHash#"CM volumes"=={12, 12, 8, 4});
    assert(AHash#"polar degrees"==pd);
    assert(AHash#"ED"==28);  
///

TEST ///
{*  
    restart
    needsPackage "EDPolytope"
*}
    A=transpose matrix{{1,0},{0,2},{3,2},{0,1}};
    c1=CMClass(A);
    assert(c1==CMClass(A,ring c1));
    assert(EDdeg(A)==26);
    Avals=(EDdeg(A,Output=>HashTable));
    --peek Avals
    assert(Avals#"ED"==26);
    assert(Avals#"dual degree"==7);
    assert(Avals#"degree"==7);
    assert(Avals#"polar degrees"=={7,12,7});
    assert(Avals#"CM volumes"=={4,9,7});
///