-- -*- coding: utf-8 -*-
----------------------------------------------------
----------------------------------------------------
-- previous version: 0.2 30Jun08, submitted by Josephine Yu.
-- author: Mike Stillman -- 
--     	    	      	   core
-- author: Josephine Yu -- 
--     	    	      	   all remaining functions; documentation
-- Sonja Petrovic -- 
--     	    	      	   interface for windows; edited documentation; tests
--     	    	      	   
-- latest major update: 6Jul08;  
-- small revision in Documentation: 6oct08; SP.
----------------------------------------------------
----------------------------------------------------

newPackage(
	"FourTiTwo",
    	Version => "1.0", 
    	Date => "February 8, 2009",
    	Authors => {
	     {Name => "Mike Stillman", Email => "mike@math.cornell.edu"},
	     {Name => "Josephine Yu", Email => "jyu@math.mit.edu"},
	     {Name => "Sonja Petrovic", Email => "petrovic@psu.edu"}
	     },
    	Headline => "Interface to 4ti2",
	Keywords => {"Interfaces"},
	Configuration => { "path" => "",
	     "keep files" => true
	      },
    	DebuggingMode => false
    	)

export {
     "toBinomial",
     "getMatrix",
     "putMatrix",
     "toricMarkov",
     "toricGroebner",
     "toricCircuits",
     "toricGraver",
     "hilbertBasis",
     "rays",
     "InputType",
     "toricGraverDegrees"
     }

-- for backward compatibility
if not programPaths#?"4ti2" and FourTiTwo#Options#Configuration#"path" != ""
    then programPaths#"4ti2" = FourTiTwo#Options#Configuration#"path"

fourTiTwo = null

run4ti2 = (exe, args) -> (
    if fourTiTwo === null then
	fourTiTwo = findProgram("4ti2", "markov -h",
	    Prefix => {(".*", "4ti2-"), -- debian
		       (".*", "4ti2_")}, -- suse
	    AdditionalPaths =>
		{"/usr/lib/4ti2/bin", "/usr/lib64/4ti2/bin"}); -- fedora
     runProgram(fourTiTwo, exe, args)
)

getFilename = () -> (
     filename := temporaryFileName();
     while fileExists(filename) or fileExists(filename|".mat") or fileExists(filename|".lat") do filename = temporaryFileName();
     filename)

putMatrix = method()
putMatrix(File,Matrix) := (F,B) -> (
     B = entries B;
     numrows := #B;
     numcols := #B#0;
     F << numrows << " " << numcols << endl;
     for i from 0 to numrows-1 do (
	  for j from 0 to numcols-1 do (
	       F << B#i#j << " ";
	       );
	  F << endl;
	  );
     )

getMatrix = method()
getMatrix String := (filename) -> (
     L := lines get filename;
     l := toString first L;
     v := value("{" | replace(" +",",",l)|"}"); 
     dimensions := select(v, vi -> vi =!= null);
     if dimensions#0 == 0 then ( -- matrix has no rows
	  matrix{{}}
     ) else(
     	  L = drop(L,1);
     	  --L = select(L, l-> regex("^[:space:]*$",l) === null);--remove blank lines
     	  matrix select(
	       apply(L, v -> (w := value("{" | replace(" +",",",v)|"}"); w = select(w, wi -> wi =!= null))),
	       row -> row =!= {}
     	       ) 
     ))    

toBinomial = method()
toBinomial(Matrix,Ring) := (M,S) -> (
     toBinom := (b) -> (
       pos := 1_S;
       neg := 1_S;
       scan(#b, i -> if b_i > 0 then pos = pos*S_i^(b_i)
                   else if b_i < 0 then neg = neg*S_i^(-b_i));
       pos - neg);
     ideal apply(entries M, toBinom)
     )

toricMarkov = method(Options=> {InputType => null})
toricMarkov Matrix := Matrix => o -> (A) -> (
     filename := getFilename();
     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
     if o.InputType === "lattice" then
     	  F := openOut(filename|".lat")
     else 
       	  F = openOut(filename|".mat");
     putMatrix(F,A);
     close F;
     run4ti2("markov", rootPath | filename);
     getMatrix(filename|".mar")
     )
toricMarkov(Matrix,Ring) := o -> (A,S) -> toBinomial(toricMarkov(A,o), S)

toricGroebner = method(Options=>{Weights=>null})
toricGroebner Matrix := o -> (A) -> (
     filename := getFilename();
     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
     F := openOut(filename|".mat");
     putMatrix(F,A);
     close F;
     if o.Weights =!= null then (
	  cost := concatenate apply(o.Weights, x -> (x|" "));
	  (filename|".cost") << "1 " << #o.Weights << endl << cost << endl  << close;
	  );
     run4ti2("groebner", rootPath | filename);
     getMatrix(filename|".gro")
     )
toricGroebner(Matrix,Ring) := o -> (A,S) -> toBinomial(toricGroebner(A,o), S)

toricCircuits = method()
toricCircuits Matrix := Matrix => (A ->(
     filename := getFilename();
     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
     F := openOut(filename|".mat");
     putMatrix(F,A);
     close F;
     run4ti2("circuits", rootPath | filename);
     getMatrix(filename|".cir")
     ))

toricGraver = method()
toricGraver Matrix := Matrix => (A ->(
     filename := getFilename();
     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
     F := openOut(filename|".mat");
     putMatrix(F,A);
     close F;
     run4ti2("graver -q ", rootPath | filename);
     getMatrix(filename|".gra")
     ))
toricGraver (Matrix,Ring) := Ideal => ((A,S)->toBinomial(toricGraver(A),S))

hilbertBasis = method(Options=> {InputType => null})
hilbertBasis Matrix := Matrix => o -> (A ->(
     filename := getFilename();
     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
     if o.InputType === "lattice" then
     	  F := openOut(filename|".lat")
     else 
       	  F = openOut(filename|".mat");
     putMatrix(F,A);
     close F;
     run4ti2("hilbert", rootPath | filename);
     getMatrix(filename|".hil")
     ))


rays = method()
rays Matrix := Matrix => (A ->(
     filename := getFilename();
     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
     F := openOut(filename|".mat");
     putMatrix(F,A);
     close F;
     run4ti2("rays", rootPath | filename);
     getMatrix(filename|".ray")
     ))

-- SP: the output command interface
-- I would like to have a command that gives the list of degrees of Graver/Groebner/Circuit/Markov file;
-- the way 4ti2 does this is you tell it the whatever.mar or whatever.cir file and it writes the degrees
-- to the screen.
-- On the other hand, it doesn't matter because you can ask M2 for those degrees directly! 
toricGraverDegrees = method()
toricGraverDegrees Matrix := Matrix => (A ->(
     filename := getFilename();
     if debugLevel >= 1 then << "using temporary file name " << filename << endl;
     F := openOut(filename|".mat");
     putMatrix(F,A);
     close F;
     run4ti2("graver", rootPath | filename);
     ret := run4ti2("output", "--degrees " | rootPath | filename|".gra");
     print ret#"output"
     ))


beginDocumentation()

doc ///
     Key 
          FourTiTwo
     Headline
     	  Interface for 4ti2
     Description
          Text
	       Interfaces most of the functionality of the software {\tt 4ti2} available at
	       @HREF"http://www.4ti2.de/"@.
	       (The user needs to have {\tt 4ti2} installed on his/her machine.)
	        
	       A $d\times n$ integral matrix $A$ (with nonnegative entries) specifies a map from a polynomial 
	       ring in d variables to a polynomial ring with n variables by specifying exponents of the variables indexing
	       its columns. For example, if $A$ is a matrix 
	       $$\begin{pmatrix}
	       3&2&1&0\\
	       0&1&2&3
	       \end{pmatrix}$$
	       the map from $k[s,t]$ to $k[a,b,c,d]$ is given by
	       $(s,t) \mapsto \ (s^3,s^2t,st^2,t^3)$.
	       
	       The toric ideal $I_A$ is the kernel of this map. 
	       It is minimally generated by the 2-minors of the matrix
	       $$\begin{pmatrix}
	       x&y&z\\
	       y&z&w
	       \end{pmatrix}$$
	       Given the matrix $A$, one can compute its lattice basis ideal specified by the integral basis
	       of the lattice $A$, the toric ideal $I_A$, its Groebner bases, etc. In practice, however, 
	       these are nontrivial computational tasks. 
	       The software {\tt 4ti2} is very efficient in computing these objects. 	      
	       
	       For more theoretical details (and more generality), see the standard reference: 
	       B. Sturmfels, {\bf Gr\"obner bases and convex polytopes.} 
	       American Mathematical Society, University Lectures Series, No 8, Providence, Rhode Island, 1996. 
	       
               {\bf Note for cygwin users:} 
	       If a problem occurs during package installation and/or loading, it should be fixed 
	       by setting the path inside the file {\tt .Macaulay2/init-FourTiTwo.m2}  to whatever folder {\tt 4ti2} is installed.
	       For example, if  {\tt 4ti2} has been installed in {\tt C:/cygwin/4ti2/win32}, then the line 
	       inside the {\tt init-FourTiTwo.m2} file will look like this:  {\tt "path" => "C:/cygwin/4ti2/win32/"}  .
	       Alternately, the path for {\tt 4ti2} may be set when loading the package using the following command:
	         loadPackage("FourTiTwo", Configuration=>{"path"=>"C:/cygwin/4ti2/win32/"})
	       assuming that 4ti2 has been installed in C:/cygwin/4ti2/win32.
      	       
     Caveat
      	       If the package SimpleDoc is not found when installing {\tt FourTiTwo.m2}, 
	       see questions and answers 6, 7, and 8 on the Macaulay 2 web site.	       
///;

doc ///
    Key
	getMatrix
	(getMatrix, String)
    Headline
	reads a matrix from a 4ti2-formatted input file
    Usage
	getMatrix s
    Inputs
    	s:String
	     file name
    Outputs
	A:Matrix
    Description
	Text
	     The file should contain a matrix in the format such as

	     2 4\break
	     1 1 1 1\break
	     1 2 3 4\break

	     The first two numbers are the numbers of rows and columns.
    SeeAlso
	putMatrix
///;

doc ///
     Key 
	  putMatrix
     	  (putMatrix,File,Matrix)
     Headline
     	  writes a matrix into a file formatted for 4ti2
     Usage
     	  putMatrix(F,A)
     Inputs
     	  F:File
       	  A:Matrix
     Description
     	  Text
		Write the matrix {\tt A} in file {\tt F} in {\tt 4ti2} format.
	  Example
		A = matrix "1,1,1,1; 1,2,3,4"
		s = temporaryFileName()
		F = openOut(s)
		putMatrix(F,A)
		close(F)
		getMatrix(s)
     SeeAlso
     	  getMatrix
///;

doc ///
     Key
          toBinomial
     	  (toBinomial, Matrix, Ring)	  
     Headline
     	  creates a toric ideal from a given set of exponents of its generators
     Usage
     	  toBinomial(M,R)
     Inputs
     	  M: Matrix
	  R: Ring
	       ring with as least as many generators as the columns of {\tt M}
     Outputs
     	  I: Ideal
     Description
     	  Text
	       Equivalent to "output --binomials" in 4ti2.
	       Returns the ideal in the ring {\tt R} generated by the binomials corresponding to rows of {\tt M}.
	  Example
	       A = matrix "1,1,1,1; 1,2,3,4"
	       B = syz A 
	       R = QQ[a..d]	     
	       toBinomial(transpose B,R)
///;

doc ///
     Key
     	  toricGroebner
          (toricGroebner, Matrix)
     	  (toricGroebner, Matrix, Ring)
     	  [toricGroebner, Weights]
     Headline
     	  calculates a Groebner basis of the toric ideal I_A, given A; invokes "groebner" from 4ti2
     Usage
     	  toricGroebner(A) or toricGroebner(A,R)
     Inputs
      	  A:Matrix    
	       whose columns parametrize the toric variety. The toric ideal $I_A$ is the kernel of the map defined by {\tt A}.
     	  R:Ring
	       ring with as least as many generators as the columns of {\tt A}
     Outputs
     	  B:Matrix 
	       whose rows give binomials that form a Groebner basis of the toric ideal of {\tt A}
     	  I:Ideal
	       whose generators form a Groebner basis for the toric ideal
     Description
	  Example
	       A = matrix "1,1,1,1; 1,2,3,4"
	       toricGroebner(A)
	  Text
	       Note that the output of the command is a matrix whose rows are the exponents of the binomials that for a Groebner basis of the 
	       toric ideal $I_A$.
	       As a shortcut, one can ask for the output to be an ideal instead:
	  Example
	       R = QQ[a..d]
	       toricGroebner(A,R)
	  Text
	       {\tt 4ti2} offers the use of weight vectors representing term orders, as follows:
	  Example
	       toricGroebner(A,Weights=>{1,2,3,4})
     Caveat
      	       It seems that some versions of 4ti2 do not pick up on the weight vector. It may be better to run gb computation in M2 directly with specified weights.	       

 ///;
 
 doc ///
     Key
     	  toricMarkov
          (toricMarkov, Matrix)
	  (toricMarkov, Matrix, Ring)
	  [toricMarkov, InputType]
     Headline
     	  calculates a generating set of the toric ideal I_A, given A; invokes "markov" from 4ti2
     Usage
     	  toricMarkov(A) or toricMarkov(A, InputType => "lattice") or toricMarkov(A,R)
     Inputs
      	  A:Matrix
	       whose columns parametrize the toric variety; the toric ideal is the kernel of the map defined by {\tt A}.
	       Otherwise, if InputType is set to "lattice", the rows of {\tt A} are a lattice basis and the toric ideal is the 
	       saturation of the lattice basis ideal.	       
	  InputType=>String
	       which is the string "lattice" if rows of {\tt A} specify a lattice basis
	  R:Ring
	       polynomial ring in which the toric ideal $I_A$ should live
     Outputs
     	  B:Matrix 
	       whose rows form a Markov Basis of the lattice $\{z {\rm integral} : A z = 0\}$
	       or the lattice spanned by the rows of {\tt A} if the option {\tt InputType => "lattice"} is used
     Description
     	  Text
	       Suppose we would like to comput the toric ideal defining the variety parametrized by the following matrix:
	  Example
	       A = matrix"1,1,1,1;0,1,2,3"
	  Text
	       Since there are 4 columns, the ideal will live in the polynomial ring with 4 variables.
	  Example
	       R = QQ[a..d]
	       M = toricMarkov(A)
	  Text
	       Note that rows of M are the exponents of minimal generators of $I_A$.  To get the ideal, we can do the following:
	  Example
      	       I = toBinomial(M,R)
	  Text
	       Alternately, we might wish to give a lattice basis ideal instead of the matrix A. The lattice basis will be specified 
	       by a matrix, as follows:
	  Example
	       B = syz A 
	       N = toricMarkov(transpose B, InputType => "lattice")	  
	       J = toBinomial(N,R) -- toricMarkov(transpose B, R, InputType => "lattice")	     
	  Text
	       We can see that the two ideals are equal:
	  Example
     	       I == J
	  Text
	       Also, notice that instead of the sequence of commands above, we could have used the following:
	  Example
	       toricMarkov(A,R)	       
///;

doc ///
     Key
     	  toricGraver
          (toricGraver, Matrix)
     	  (toricGraver, Matrix, Ring)
     Headline
     	  calculates the Graver basis of the toric ideal; invokes "graver" from 4ti2
     Usage
     	  toricGraver(A) or toricGraver(A,R)
     Inputs
      	  A:Matrix    
	       whose columns parametrize the toric variety. The toric ideal $I_A$ is the kernel of the map defined by {\tt A}
	  R:Ring
	       polynomial ring in which the toric ideal $I_A$ should live
     Outputs
     	  B:Matrix 
	       whose rows give binomials that form the Graver basis of the toric ideal of {\tt A}, or
     	  I:Ideal
	       whose generators form the Graver basis for the toric ideal
     Description
     	  Text
	       The Graver basis for any toric ideal $I_A$ contains (properly) the union of all reduced Groebner basis of $I_A$.  
	       Any element in the Graver basis of the ideal is called a primitive binomial.
	  Example
	       A = matrix "1,1,1,1; 1,2,3,4"
	       toricGraver(A)
	  Text
	       If we prefer to store the ideal instead, we may use:
	  Example
	       R = QQ[a..d]
	       toricGraver(A,R)
	  Text
	       Note that this last ideal equals the toric ideal $I_A$ since every Graver basis element is actually in $I_A$.
///;

doc ///
     Key
     	  toricGraverDegrees
          (toricGraverDegrees, Matrix)
     Headline
     	  displays the degrees of all Graver basis elements for the toric ideal I_A
     Usage
     	  toricGraverDegrees(A) 
     Inputs
      	  A:Matrix    
	       whose columns parametrize the toric variety. The toric ideal $I_A$ is the kernel of the map defined by {\tt A}
     Description
     	  Text
	       Equivalent to "output --degrees foo.gra" in 4ti2.
	       Very often the Graver basis consists of too many binomials, and one is only interested in their degrees. In this case,
	       instead of looking at the Graver basis of $I_A$, we may just want to look for the degrees of binomials which show up:
	  Example
	       A = matrix "1,1,1,1; 1,2,3,4"
	       toricGraver(A) -- the Graver basis
	       toricGraverDegrees(A) -- just the degrees
	  Text
	       Note that these are all 1-norms of the vectors. Since $I_A$ is homogeneous, there are 3 binomials of degree 2 (norm 4) 
	       and 2 binomials of degree 3 (norm 6).
	       
	       Here is a more complicated example, where one may not want to see the Graver basis elements explicitly. 
	       The toric ideal I_M is the ideal of the rational normal scroll S(3,2,3):
	  Example
	       M = matrix "1,1,1,1,1,1,1,1,1,1,1; 1,1,1,1,0,0,0,0,0,0,0; 0,0,0,0,1,1,1,0,0,0,0; 0,0,0,0,0,0,0,1,1,1,1; 1,2,3,4,1,2,3,1,2,3,4"
	       toricGraverDegrees(M)
	  Text
	       Here is another example where with many Graver basis elements. The following matrix is a design matrix for a particular statistical model for a 4-node p1 network; see Fienberg-Petrovic-Rinaldo.
	  Example
	        A = matrix "1,0,1,1,0,1,1,0,1,0,0,0,0,0,0,0,0,0;0,1,1,0,0,0,0,0,0,1,0,1,1,0,1,0,0,0;0,0,0,0,1,1,0,0,0,0,1,1,0,0,0,1,0,1;0,0,0,0,0,0,0,1,1,0,0,0,0,1,1,0,1,1;0,1,1,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0;1,0,1,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0;0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,0,1,1;0,0,0,0,0,0,1,0,1,0,0,0,1,0,1,1,0,1;0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1";
                toricGraverDegrees(A)
///;


doc ///
     Key
     	  hilbertBasis
          (hilbertBasis, Matrix)
	  [hilbertBasis, InputType]
     Headline
     	  calculates the Hilbert basis of the cone; invokes "hilbert" from 4ti2
     Usage
     	  hilbertBasis(A) or hilbertBasis(A, InputType => "lattice")
     Inputs
      	  A:Matrix    
	       defining the cone $\{z : Az = 0, z \ge 0 \}$
     Outputs
     	  B:Matrix 
	       whose rows form the Hilbert basis of the cone $\{z : Az = 0, z \ge 0 \}$
	       or the cone $\{z A : z {\rm {} is an integral vector and } z A \ge 0 \}$ if {\tt InputType => "lattice"} is used
     Description
	  Example
	       A = matrix "1,1,1,1; 1,2,3,4"
	       B = syz A
	       hilbertBasis(transpose B)
     	       hilbertBasis(A, InputType => "lattice")     	       
///;


doc ///
     Key
     	  rays
          (rays, Matrix)
     Headline
     	  calculates the extreme rays of the cone; invokes "rays" from 4ti2
     Usage
     	  rays(A)
     Inputs
      	  A:Matrix   
	       defining the cone $\{z : Az = 0, z \ge 0 \}$
     Outputs
     	  B:Matrix 
	       whose rows are the extreme rays of the cone $\{z : Az = 0, z \ge 0 \}$
     Description
	  Example
	       A = matrix "1,1,1,1; 1,2,3,4"
	       B = syz A
	       rays(transpose B)
///;



doc ///
     Key
     	  toricCircuits
          (toricCircuits, Matrix)
     Headline
     	  calculates the circuits of the toric ideal; invokes "circuits" from 4ti2
     Usage
     	  toricCircuits(A)
     Inputs
      	  A:Matrix    
               whose columns parametrize the toric variety. The toric ideal $I_A$ is the kernel of the map defined by {\tt A} 
     Outputs
     	  B:Matrix 
	       whose rows form the circuits of A
     Description
     	  Text
	       The circuits are contained in the Graver basis of $I_A$. In fact, they are precisely the primitive binomials in the ideal
	       with  minimal support.
	  Example
	       A = matrix "1,1,1,1; 1,2,3,4"
	       C = toricCircuits A
	  Text 
	       The ideal generated by the circuits of A in general differs from the toric ideal of A. For example:
	  Example
	       R = QQ[a..d]
	       Icircuit = toBinomial(toricCircuits(A), R) -- this is the circuit ideal of A
	       I = toBinomial(toricMarkov(A), R)
	       I==Icircuit  
	  Text
	       The two ideals are not the same. There is a minimal generator of I which is not a circuit:
	  Example
	       a*d-b*c % I -- this binomial is in I:
	       a*d-b*c % Icircuit -- but not in Icircuit:
///;

doc ///
     Key
     	  InputType
     Description
          Text
     	      Put {\tt InputType => "lattice"} as an argument in the functions toricMarkov and hilbertBasis
     SeeAlso
     	  toricMarkov
	  hilbertBasis
///;


TEST/// 
  A = matrix "1,1,1,1; 1,2,3,4"
  M = toricMarkov(A)
  R = QQ[x_0,x_1,x_2,x_3]
  I = toBinomial(M,R)
  Irnc3 = ideal(x_0*x_2-x_1^2,x_1*x_3-x_2^2,x_0*x_3-x_1*x_2)
  assert(I==Irnc3)
///
TEST ///   
  B = matrix "1,-2,1,0; 0,1,-2,1"
  M = toricMarkov(B, InputType => "lattice")
  R = QQ[x_0,x_1,x_2,x_3]
  I = toBinomial(M,R)
  Irnc3 = ideal(x_0*x_2-x_1^2,x_1*x_3-x_2^2,x_0*x_3-x_1*x_2)
  assert(I== Irnc3)  
///
TEST ///     
  R=CC[x_0,x_1,x_2,x_3]
  A = matrix "1,1,1,1; 1,2,3,4"
  C = toricCircuits(A)  --circuits returned by 4ti2
  Icir = toBinomial(C,R) -- circuit ideal returned by 4ti2
  Ctrue = matrix{{0,1,-2,1},{1,-2,1,0},{1,0,-3,2},{2,-3,0,1}} --known: all circuits
  IcirTrue = toBinomial(Ctrue,R) --known: circuit ideal
  Irnc3 = ideal(x_0*x_2-x_1^2,x_1*x_3-x_2^2,x_0*x_3-x_1*x_2)
  assert(Icir==IcirTrue)
  shouldBeFalse = (Icir==Irnc3)
  assert(shouldBeFalse==false)
  assert(target C == target Ctrue)
  assert(source C == source Ctrue)
///
TEST ///     
  B = matrix "1,-2,1,0; 0,1,-2,1"  
  R = QQ[a..d]
  I = toBinomial(B,R)
  assert(a*c-b^2 % I == 0)
  assert(a*c-d^2 %I =!= 0)
  assert(degree I == 4)
  M = hilbertBasis B
  assert(numrows M == 3)
  assert(numcols M == 4)     
  M1 = rays B
  assert(numrows M1 == 2)
  assert(numcols M1 == 4)
///
TEST///
  A = matrix "1,0,1,1,0,1,1,0,1,0,0,0,0,0,0,0,0,0;0,1,1,0,0,0,0,0,0,1,0,1,1,0,1,0,0,0;0,0,0,0,1,1,0,0,0,0,1,1,0,0,0,1,0,1;0,0,0,0,0,0,0,1,1,0,0,0,0,1,1,0,1,1;0,1,1,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0;1,0,1,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0;0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,0,1,1;0,0,0,0,0,0,1,0,1,0,0,0,1,0,1,1,0,1;0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1"
  M = toricGroebner(A); --note this matrix is the design matrix for the p1 statistical model on 4 nodes using a constant rho. (see fienberg/rinaldo/petrovic; in prep-missing reference).
  assert(numrows M == 137)
  assert(numcols M == 18)
  R = QQ[x_1..x_18]
  I = toBinomial(M,R);
  assert(degree I == 192)
  A1 = matrix "3,2,1,0;0,1,2,3" --one more example
  R=QQ[x_0..x_3]
  G4ti2=gens toricGroebner(A1,R)
  GM2 =gens gb toricMarkov(A1,R)
  Gtrue=toList flatten entries GM2  
  G = toList flatten entries G4ti2
  apply(0..#Gtrue-1, j->  (isSubset({Gtrue_j},G) )) --checking 4ti2's gb against M2's gb
  assert(numrows GM2 == numrows G4ti2)
  assert(numcols GM2 == numcols G4ti2)
  Rwt=QQ[x_0..x_3,Weights=>{3,2,4,1}] --with wt vector
  G4ti2=gens toBinomial(toricGroebner(A1,Weights=>{3,2,4,1}),Rwt) 
  GM2=gens gb toricMarkov(A1,Rwt)
  Gtrue=toList flatten entries GM2  
  G = toList flatten entries G4ti2
  assert(       numrows GM2 == numrows G4ti2       )
  assert(       numcols GM2 == numcols G4ti2       )
  apply(0..#Gtrue-1, j-> assert(isSubset({Gtrue_j},G)) ) --checking 4ti2's gb against M2's gb
///
TEST///
  needsPackage "FourTiTwo"   --testing graver  
  A1 = matrix "3,2,1,0;0,1,2,3"
  R=QQ[x_0..x_3]
  G = toricGraver(A1)
  assert( numrows G==5)
  assert(numcols G==4)
  Gtrue = toBinomial(matrix{{1,-2,1,0},{0,1,-2,1},{1,-1,-1,1},{2,-3,0,1},{1,0,-3,2}},R) --known: Graver basis
  Gtrue=toList flatten entries gens Gtrue  
  G = toList flatten entries gens toBinomial(G,R)
  apply(0..#Gtrue-1, j->  assert(isSubset({Gtrue_j},G)) ) --testing 4ti2 output against by-hand calculation!
  A = matrix "1,0,1,1,0,1,1,0,1,0,0,0,0,0,0,0,0,0;0,1,1,0,0,0,0,0,0,1,0,1,1,0,1,0,0,0;0,0,0,0,1,1,0,0,0,0,1,1,0,0,0,1,0,1;0,0,0,0,0,0,0,1,1,0,0,0,0,1,1,0,1,1;0,1,1,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0;1,0,1,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0;0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,0,1,1;0,0,0,0,0,0,1,0,1,0,0,0,1,0,1,1,0,1;0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1"
  M = toricGraver(A);   --note this matrix is the design matrix for the p1 statistical model on 4 nodes using a constant rho. (see fienberg/rinaldo/petrovic; in prep-missing reference).
  assert(numrows M == 7462)
  assert(numcols M == 18)
  AS=matrix"1,1,1,1,1,1,1,1;1,1,1,0,0,0,0,0;1,2,3,1,2,3,4,5"--scroll S(2,4)
  R=QQ[x_1..x_8]
  G4ti2 = toList flatten entries gens toBinomial(toricGraver(AS),R)
  assert(#G4ti2 == 82)
  Gtrue=toList flatten entries gens toBinomial(matrix" 1,-2,1,0,0,0,0,0;1,-1,0,-1,1,0,0,0;2,-2,0,-1,0,1,0,0;3,-3,0,-1,0,0,1,0;4,-4,0,-1,0,0,0,1;3,-3,0,0,-1,0,0,1;2,-2,0,0,0,-1,0,1;1,-1,0,0,0,0,-1,1;2,-2,0,0,-1,0,1,0;1,-1,0,0,0,-1,1,0;1,-1,0,0,-1,1,0,0;0,1,-1,-1,1,0,0,0;0,1,-1,0,0,-1,1,0;0,1,-1,0,0,0,-1,1;0,1,-1,0,-1,1,0,0;1,0,-1,-1,0,1,0,0;1,0,-1,0,-1,0,1,0;1,0,-1,0,0,-1,0,1;2,-1,-1,-1,0,0,1,0;2,-1,-1,0,-1,0,0,1;3,-2,-1,-1,0,0,0,1;0,2,-2,0,-1,0,1,0;1,1,-2,-1,0,0,1,0;0,2,-2,-1,0,1,0,0;2,0,-2,-1,0,0,0,1;0,2,-2,0,0,-1,0,1;1,1,-2,0,-1,0,0,1;0,0,0,1,-1,-1,1,0;0,0,0,1,-1,0,-1,1;0,0,0,1,-2,1,0,0;1,-1,0,1,-2,0,1,0;1,-1,0,1,-1,-1,0,1;2,-2,0,1,-2,0,0,1;0,0,0,0,0,1,-2,1;0,0,0,0,1,-2,1,0;1,-1,0,-1,0,2,-1,0;0,0,0,0,1,-1,-1,1;1,-1,0,-1,0,1,1,-1;1,-1,0,0,-1,0,2,-1;2,-2,0,-1,0,0,2,-1;1,-1,0,0,1,-2,0,1;0,0,0,1,0,-2,0,1;0,1,-1,0,1,-2,0,1;0,1,-1,-1,0,1,1,-1;0,1,-1,-1,0,2,-1,0;1,0,-1,-1,0,0,2,-1;0,1,-1,0,-1,0,2,-1;0,1,-1,1,-1,-1,0,1;0,1,-1,1,-2,0,1,0;1,0,-1,1,-2,0,0,1;0,3,-3,-1,0,0,1,0;1,2,-3,-1,0,0,0,1;0,3,-3,0,-1,0,0,1;0,2,-2,1,-2,0,0,1;0,2,-2,-1,0,0,2,-1;0,4,-4,-1,0,0,0,1;0,0,0,1,-2,0,2,-1;0,0,0,2,-2,-1,0,1;0,0,0,2,-3,0,1,0;1,-1,0,2,-3,0,0,1;0,0,0,1,0,-3,2,0;0,0,0,1,0,-1,-2,2;1,-1,0,-1,0,0,3,-2;0,0,0,0,1,0,-3,2;0,0,0,0,2,-3,0,1;0,1,-1,-1,0,0,3,-2;0,1,-1,2,-3,0,0,1;0,0,0,3,-4,0,0,1;0,0,0,1,0,0,-4,3;1,0,-1,0,-1,1,-1,1;1,0,-1,0,-2,2,0,0;1,0,-1,-1,1,-1,1,0;1,0,-1,-1,1,0,-1,1;1,0,-1,-2,2,0,0,0;1,0,-1,0,0,-2,2,0;1,0,-1,0,0,0,-2,2;2,0,-2,0,-2,1,0,1;2,0,-2,-2,1,0,1,0;2,0,-2,-1,0,-1,2,0;2,0,-2,0,-1,0,-1,2;3,0,-3,-2,0,0,2,0;3,0,-3,0,-2,0,0,2"  ,R)
   apply(0..#Gtrue-1, j->  assert(isSubset({Gtrue_j},G4ti2)) ) -- checking 4ti2 output against by hand input!!
///


end



restart
--load "4ti2.m2"
installPackage ("FourTiTwo", RemakeAllDocumentation => true, UserMode=>true)
installPackage("FourTiTwo",UserMode=>true,DebuggingMode => true)
viewHelp FourTiTwo


check FourTiTwo

debug FourTiTwo


A = matrix{{1,1,1,1},{0,1,2,3}}
A = matrix{{1,1,1,1},{0,1,3,4}}
B = syz A
time toricMarkov A
A
toricMarkov(A, InputType => "lattice")
R = QQ[a..d]
time toricGroebner(A)
toBinomial(transpose B, R)
toricCircuits(A)
H = hilbertBasis(A)
hilbertBasis(transpose B)
toBinomial(H,QQ[x,y])
toricGraver(A)
A
toricMarkov(A)

7 9
A = matrix"
1,1,1,-1,-1,-1, 0, 0, 0;
1,1,1, 0, 0, 0,-1,-1,-1;
0,1,1,-1, 0, 0,-1, 0, 0;
1,0,1, 0,-1, 0, 0,-1, 0;
1,1,0, 0, 0,-1, 0, 0,-1;
0,1,1, 0,-1, 0, 0, 0,-1;
1,1,0, 0,-1, 0,-1, 0, 0"
transpose A
toricMarkov transpose A
hilbertBasis transpose A
toricGraver transpose A
toricCircuits transpose A

27 27
A = matrix"
1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0;
0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0;
0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0;
0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0;
0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0;
0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0;
0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0;
0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0;
0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1;
1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1;
1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1"
toricMarkov A
R = QQ[x_1..x_27]
toricMarkov(A,R)
toricGroebner(A,R)
gens gb oo

I = toBinomial(matrix{{}}, QQ[x])
gens I
gens gb I


-- Notes from Mike talking with Peter Malkin, 4/21/09

in 1.3.2 (and in 1.3.1):

These routines use the structure below:

groebner, markov, 
hilbert, graver, zsolve
rays, circuits, qsolve

also: minimise, walk, normalform

a.mat: m by n
a.rel: 1 by m: symbols: >, =, <  (means: >= 0, == 0, <= 0)
a.sign: 1 by n matrix: 0,1,-1,2

a.sign: 0:  x_i unrestricted in sign
        1:  x_i >= 0
	-1: x_i <= 0.
	2:  x_i is a Graver component

a.mat, a.rel, a.sign.

groebner, markov, zsolve: also can give a.rhs (1 x n matrix).

groebner does this:
Ax >= 0, x>=0.

Ax-Iy = 0, x >= 0, y >= 0.

doc on main page of 4ti2 web site: manual, and the slides.
.rhs doesn't work for qsolve though, possibly.

for hilbert, graver, zsolve, have the following filter:
.ub, .lb can be used to provide upper and lower bounds (lower only for Graver components, or for <= vars).

install glpk first, and make sure gmp is visible.
./configure --with-gmp=.... --with-glpk=....

call the different routines as: -p32, -p64, -parb

tests are in 
4ti2-1.3.2/test

email Peter if I have more questions
4ti2 google group: joined.

glpk:open source linear programming