-----------------------------------------------------------------
-- Preamble
-----------------------------------------------------------------
newPackage(
	"CharToric",
	Version => "1.1", 
    	Date => "April 20, 2015",
    	Authors => {{Name => "Martin Helmer", 
		  Email => "martin.helmer2@gmail.com", 
		  HomePage => "http://publish.uwo.ca/~mhelmer2/"}},
    	Headline => "Computes CSM classes, Segre classes and the Euler Char.",
    	DebuggingMode => false,
	Reload => true
    	)

needsPackage "NormalToricVarieties";

--Exported functions/variables
export{
   EulerToric,
   CSMToric,
   ChowRing,
   CheckSmooth,
   mult
}
mult=method(TypicalValue=>ZZ);
mult:=RayMatrix->(
    r:=rank RayMatrix;
    return multr(RayMatrix,r);
    )

ChowRing=method(TypicalValue=>QuotientRing);
ChowRing NormalToricVariety:=TorVar->(
    --assert isSimplicial TorVar;
    --assert isComplete TorVar;
    --First build Chow ring, need Stanley-Reisner Ideal (SR) and the ideal
    -- generated by the linear relations of the rays (J)
    --See Cox, Little, Schenck Th. 12.5.3 and comments after proof
    R:=ring(TorVar);
    A:=0;
    --For simplical toric var. Lemma 3.5 of Euler characteristic of coherent sheaves on simplicial torics via the Stanley-Reisner ring
    -- (and probally other sources) tell us that the SR ideal is the Alexander
    --dual of the toric irrelevant ideal 
    SR:=dual monomialIdeal TorVar;
    F:=fan TorVar;
    Fd:=dim(F);
    --Build ideal generated by linear relations of the rays
    Jl:={};
    for j from 0 to dim(F)-1 do(
	Jl=append (Jl,sum(0..((length rays(TorVar))-1), i->(((rays TorVar)_i)_j)*(gens R)_i ));
    	);
    J:=ideal(Jl);
    --Chow ring
    if isSmooth(TorVar) then(
	--if smooth our Chow ring should be over ZZ
	 C:=ZZ[gens R, Degrees=>degrees R, Heft=>heft R];
         A=C/substitute(SR+J,C);
	 )
     else (A=R/(SR+J););
    --Generators (as a ring) of the quotient ring representation of the chow ring corespond to 
    --the divisors associated to the rays in the fan Theorem 12.5.3. Cox, Little, Schenck and 
    --comments above
    return A;
    
    
    );
EulerToric=method(TypicalValue => RingElement,Options => {CheckSmooth=>true});


EulerToric NormalToricVariety :=opts->TorVar->(
    A:=ChowRing(TorVar);
    return EulerToric(A,TorVar,CheckSmooth=>opts.CheckSmooth);
    );
EulerToric (QuotientRing, NormalToricVariety) :=opts->(A,TorVar)->(
    --For Euler we need only the zero dimensional part of csm, so just calculate that part
    L:=gens(A);
    Trays:=rays TorVar;
    csm0:=0_A;
    Rmat:=0;
    prodj:=0;
    Ssets:=0;
    indsubsets:=0;
    dimTorVar:=dim(TorVar);
    TorVarIsSmooth:=false;
    if opts.CheckSmooth==true then TorVarIsSmooth=isSmooth(TorVar);
    if  TorVarIsSmooth then(
	indsubsets=subsets((0..numgens(A)-1),dimTorVar);
	Ssets=for l in indsubsets list L_l;
	csm0=sum(0..(length(Ssets)-1),j-> product(Ssets_j));	
	)
    else(
	indsubsets=subsets((0..numgens(A)-1),dimTorVar);
	Ssets=for l in indsubsets list L_l;
	--csm0=sum(0..(length(Ssets)-1),j-> mult(transpose matrix Trays_(indsubsets_j))*product(Ssets_j));
	for j from 0 to length(Ssets)-1 do(
	    Rmat=transpose matrix Trays_(indsubsets_j);
	    prodj=product(Ssets_j);
	    if prodj!=0 then(
		csm0=csm0+ mult(Rmat)*prodj;
		);
	    );	   
      	);
    (m,c):=coefficients(csm0);
    Echar:=(sum entries c)_0; 
    <<"Euler Charcteristic = "<<Echar<<endl;
    return Echar;
    
    );
CSMToric=method(TypicalValue => RingElement,Options => {CheckSmooth=>true});

CSMToric NormalToricVariety :=opts->TorVar->(
    A:=ChowRing(TorVar);
    return CSMToric(A,TorVar,CheckSmooth=>opts.CheckSmooth);
    );

CSMToric (QuotientRing, NormalToricVariety) :=opts->(A,TorVar)->(
    --Generators (as a ring) of the quotient ring representation of the chow ring corespond to 
    --the divisors associated to the rays in the fan Theorem 12.5.3. Cox, Little, Schenck and 
    --comments above
    L:=gens(A);
    Trays:=rays TorVar;
    csm:=1_A;
    Rmat:=0;
    prodj:=0;
    --The following implements the method described 
    --in Barthel, Brasselet, and Fieseler.
    --Lemma 12.5.2 of Cox, Little, Schenck is used to find the chow ring class of the 
    --orbit closure from divisors
    --if the toric variety is smooth the multiplicity is 1. 
    Ssets:=0;
    indsubsets:=0;
    TorVarIsSmooth:=false;
    if opts.CheckSmooth==true then TorVarIsSmooth=isSmooth(TorVar);
    if  TorVarIsSmooth then(
     	for i from 1 to dim(TorVar) do(
	    indsubsets=subsets((0..numgens(A)-1),i);
	    Ssets=for l in indsubsets list L_l;
	    csm=csm+sum(0..(length(Ssets)-1),j-> product(Ssets_j));	
	    );
	)
    else(
	for i from 1 to dim(TorVar) do(
	    indsubsets=subsets((0..numgens(A)-1),i);
	    Ssets=for l in indsubsets list L_l;
	    --csm=csm+sum(0..(length(Ssets)-1),j-> mult(transpose matrix Trays_(indsubsets_j))*product(Ssets_j));
	    for j from 0 to length(Ssets)-1 do(
	    	Rmat=transpose matrix Trays_(indsubsets_j);
	    	prodj=product(Ssets_j);
		if prodj!=0 then(
		    csm=csm+ mult(Rmat)*prodj;
		    );
	    	);		
	    );
     	);
    <<"csm = "<<csm<<endl;
    return csm;
)

---------------------------
--Internal functions 
---------------------------

--Find the multiplicity (or index) of a simplicial cone defined by the
-- rays given by the columns of the input RayMatrix pg. 300
-- Cox, Little, Schenck.
--pg 66-68, and others.... 
--Find the multiplicity (or index) of a simplicial cone defined by the
-- rays given by the columns of the input RayMatrix pg. 300
-- Cox, Little, Schenck. 


multr=(RayMatrix,r)->(
    m:=RayMatrix;
    if m==0 then return 0;
    if (numRows(m)==1 or numColumns(m)==1) then( return 1);
    --r:=rank m;
    if r<numRows(m) then(
	m=transpose groebnerBasis(transpose(m), Strategy=>"MGB");
	);
    if r<numColumns(m) then(
	--this shouldn't be reached when using mult in csm/euler calc
	m= groebnerBasis(m, Strategy=>"MGB");
	);
    if numRows(m)==numColumns(m) then(	
	mymult:=abs(determinant(m,Strategy =>Cofactor));
	--<<mymult<<endl;
	return mymult;
	)
    else (
	error "multiplicity computation error";
	return 0;
	);
    
    )

-------------------------------
--Examples for thesis
--
------------------------------
TEST ///
{*
    restart
    needsPackage "CharToric"
    needsPackage "NormalToricVarieties"
*}
    Rho = {{-1,-1,1},{3,-1,1},{0,0,1},{1,0,1},{0,1,1},{-1,3,1},{0,0,-1}};
    Sigma = {{0,1,3},{0,1,6},{0,2,3},{0,2,5},{0,5,6},{1,3,4},{1,4,5},{1,5,6},{2,3,4},{2,4,5}};
    X = normalToricVariety(Rho,Sigma);
    isSmooth X
    isComplete X
    isSimplicial X
    
    time A=ChowRing(X)
   
    time csm=CSMToric(A,X)
    time Ec=EulerToric(A,X)
    
///

TEST ///
{*
    restart
    needsPackage "CharToric"
    needsPackage "NormalToricVarieties"
*}    
    PP=projectiveSpace(6);
    time A=ChowRing(PP)
    
    time CSMToric(A,PP)
    time EulerToric(A,PP)
    
    time CSMToric(A,PP,CheckSmooth=>false);
    time EulerToric(A,PP,CheckSmooth=>false);
///
TEST ///
{*
    restart
    needsPackage "CharToric"
    needsPackage "NormalToricVarieties"
*}  
    PP8=projectiveSpace(8);
    PP6=projectiveSpace(6);
    PP5=projectiveSpace(5);
    PP2=projectiveSpace(2);
    
      
    --PP5xPP6  
    time A=ChowRing(PP5**PP6)
    
    time CSMToric (A,PP5**PP6)
    time EulerToric(A,PP5**PP6)
    
    time CSMToric (A,PP5**PP6,CheckSmooth=>false)
    time EulerToric(A,PP5**PP6,CheckSmooth=>false)
    
    --PP5xPP8 
    time A=ChowRing(PP5**PP8)
    
    time CSMToric (A,PP5**PP8)
    time EulerToric (A,PP5**PP8)
    
    time CSMToric (A,PP5**PP8,CheckSmooth=>false)
    time EulerToric (A,PP5**PP8,CheckSmooth=>false)
    --PP8xPP8
    time A=ChowRing(PP8**PP8)
    
    time CSMToric (A,PP8**PP8)
    time EulerToric (A,PP8**PP8)
    
    time CSMToric (A,PP8**PP8,CheckSmooth=>false)
    time EulerToric (A,PP8**PP8,CheckSmooth=>false)
    
    --PP5xPP5xPP5
    time A=ChowRing(PP5**PP5**PP5)
    
    time CSMToric (A,PP5**PP5**PP5);
    time EulerToric (A,PP5**PP5**PP5)
    
    time CSMToric (A,PP5**PP5**PP5,CheckSmooth=>false);
    time EulerToric (A,PP5**PP5**PP5,CheckSmooth=>false)
    
    --PP5xPP5xPP6 
    time A=ChowRing(PP5**PP5**PP6)
    
    time CSMToric (A,PP5**PP5**PP6);
    time EulerToric (A,PP5**PP5**PP6)
    
    time CSMToric (A,PP5**PP5**PP6,CheckSmooth=>false);
    time EulerToric (A,PP5**PP5**PP6,CheckSmooth=>false)
///   
    
TEST ///
{*
    restart
    needsPackage "CharToric"
    needsPackage "NormalToricVarieties"
*}    
    sF1=smoothFanoToricVariety(6,123);
    sF2=smoothFanoToricVariety(6,423);
    sF3=smoothFanoToricVariety(6,1007);
    
    --Smooth Fano 1
    time A=ChowRing(sF1)
    
    time CSMToric(A,sF1);
    time EulerToric(A,sF1);

    time CSMToric(A,sF1,CheckSmooth=>false);
    time EulerToric(A,sF1,CheckSmooth=>false);    
    
    
    --Smooth Fano 2
    time A=ChowRing(sF2)
    
    time CSMToric(A,sF2);
    time EulerToric(A,sF2);
    
    time CSMToric(A,sF2,CheckSmooth=>false);
    time EulerToric(A,sF2,CheckSmooth=>false); 
    
    --Smooth fano 3
    time A=ChowRing(sF3)
    
    time CSMToric(A,sF3);
    time EulerToric(A,sF3);
    
    time CSMToric(A,sF3,CheckSmooth=>false);
    time EulerToric(A,sF3,CheckSmooth=>false);
    
///     
TEST ///
{*
    restart
    needsPackage "CharToric"
    needsPackage "NormalToricVarieties"
*}    
    hiz=hirzebruchSurface 23
    time csm=CSMToric hiz
    
    
///

TEST ///
{*
    restart
    needsPackage "CharToric"
    needsPackage "NormalToricVarieties"
*}    
    PP=projectiveSpace(12);
    rP=rays(PP)    
    m=transpose matrix(rP_{0..2,5..9}) 
    j= groebnerBasis(transpose groebnerBasis(m));
    time(
	for i from 0 to 2^13 do(
	    --abs(determinant(j,Strategy =>Cofactor));
	    j= groebnerBasis(transpose groebnerBasis(m));
	    abs(determinant(j,Strategy =>Cofactor));
	    );
	)
	 
    
///