-- -*- coding: utf-8-unix -*-

-*
Copyright 2020 Carlos Amendola, Luis David Garcia Puente, Roser Homs Pons,
Olga Kuznetsova, Harshit J Motwani, Elina Robeva and David Swinarski.

You may redistribute this file under the terms of the GNU General Public
License as published by the Free Software Foundation, either version 2 of
the License, or any later version.
*-
newPackage(
     "GraphicalModelsMLE",
     Version => "1.0",
     Date => "April 19, 2022",
     Authors => {
	  {Name=> "Carlos Amendola",
	   Email=> "carlos.amendola@mis.mpg.de",
	   HomePage=>"http://www.carlos-amendola.com/"},

	  {Name => "Luis David Garcia Puente",
	   Email => "lgarciapuente@coloradocollege.edu",
	   HomePage => "https://luisgarciapuente.github.io"},

          {Name=> "Roser Homs Pons",
	   Email=> "roser.homs@tum.de",
	   HomePage=>"https://personal-homepages.mis.mpg.de/homspons/index.html"},

          {Name=> "Olga Kuznetsova",
	   Email=> "kuznetsova.olga@gmail.com",
	   HomePage=>"https://okuznetsova.com"},

          {Name=> "Harshit J Motwani",
	   Email=> "harshitmotwani2015@gmail.com",
	   HomePage=>"https://sites.google.com/view/harshitjmotwani/home"},

          {Name=> "Elina Robeva",
	   Email=> "erobeva@gmail.com",
	   HomePage=>"https://www.math.ubc.ca/~erobeva/"},

          {Name=> "David Swinarski",
	   Email=> "dswinarski@fordham.edu",
	   HomePage=>"http://faculty.fordham.edu/dswinarski"}
	  },
     Headline => "maximum likelihood estimates for graphical statistical models",
     Keywords => {"Algebraic Statistics"},
     DebuggingMode => true,
     PackageExports => {"GraphicalModels","Graphs","EigenSolver","NumericalAlgebraicGeometry","StatGraphs"}
     )
export {
    "checkPD",
    "checkPSD",
    "Solver",--optional argument in solverMLE
    "ConcentrationMatrix",-- optional argument in solverMLE
    "CovarianceMatrix", -- optional argument in scoreEquations
    "Saturate",-- optional argument in scoreEquations and solverMLE
    "jacobianMatrixOfRationalFunction",
    "MLdegree",
    "OptionsEigenSolver",--optional argument in solverMLE
    "OptionsNAG4M2",--optional argument in solverMLE
    "RealPrecision",--optional argument in solverMLE and scoreEquations
    "sampleCovarianceMatrix",
    "SampleData",-- optional argument in scoreEquations and solverMLE
    "SaturateOptions", -- optional argument in scoreEquations and solverMLE
    "scoreEquations",
    "solverMLE",
    "ZeroTolerance"
     }

--**************************--
--  INTERNAL ROUTINES       --
--**************************--



--*************************************--
--  Functions (local) used throughout  --
--*************************************--

------------------------------------------------------
-- Substitutes a list of points on a list of matrices
    -- input -  list of points from sols
    --          matrix whose entries are variables
    --          (expect it to be an inverse of a covariance matrix, Sin)
    -- output - list of matrices after substituting these values
------------------------------------------------------
genListMatrix = (L,A) ->
(
		T:= for l in L list coordinates(l);
		M:= for t in T list substitute(A,matrix{t});
    return M
);

----------------------------------------------
-- Selects all argmax for log det K- trace(S*K),
-- where K is an element of the list L.
-- We assume that L is the intersection of the
-- variety of the ideal generated by the Jacobian
-- of critical equations and the cone of PD matrices.

-- input - list L of candidate Sinv matrices (Sinv is Sigma^{-1}) and
-- sample covariance matrix V. Notation in line with scoreEquationsFromCovarianceMatrix
-- output -list of argmax
----------------------------------------------
maxMLE=(L,V)->(
    if #L==0 then  error("No critical points to evaluate");
    if #L==1 then  (E:=L_0; maxPt:=log det L_0- trace (V*L_0))
    else
    	(eval:=for Sinv in L list log det Sinv- trace (V*Sinv);
	evalReal:=for pt in eval when isReal pt list pt;
	if #evalReal==0 then  error("No critical point evaluates to a real solution");
	maxPt=max evalReal;
	indexOptimal:=positions(eval, i ->i== maxPt);
		E= for i in indexOptimal list L_i;);
    return (maxPt, E)
    );

-------------------------------------------
-- scoreEquationsInternal - function that returns
-- both the ideal and the corresponding SInv matrix.
-- The user-facing scoreEquations method returns only
-- the ideal, whereas SInv is used in solverMLE
-------------------------------------------
scoreEquationsInternal={Saturate => true, SaturateOptions => options saturate, SampleData=>true, RealPrecision=>53, CovarianceMatrix => false}>>opts->(R,U)->(
    ----------------------------------------------------
    -- Extract information about the graph
    ----------------------------------------------------
    -- Lambda
    L := directedEdgesMatrix R;
    -- Psi
    P := bidirectedEdgesMatrix R;

    -- If the mixedGraph only has undirected part, call specific function for undirected.
    if L==0 and P==0 then
    return scoreEquationsInternalUndir(R,U,opts);

    -- K
    K := undirectedEdgesMatrix R;
    ----------------------------------------------------
    -- Create an auxiliary ring and its fraction field
    -- which do not have the s variables
    ----------------------------------------------------
    -- create a new ring, lpR, which does not have the s variables
    lpR:=coefficientRing(R)[gens R-set support covarianceMatrix R];
    -- create its fraction field
    FR := frac(lpR);

    -----------------------------------------------------
    -- Compute Sinv
    -----------------------------------------------------
    -- Kinv
    K=sub(K, FR);
    Kinv:=inverse K;
    P=sub(P,FR);

     --Omega
    if K==0 then W:=P else (if P==0 then W=Kinv else W = directSum(Kinv,P));

    -- move to FR, the fraction field of lpR
    L= sub(L,FR);

    -- Sigma
    d:=numcols L;
    if L==0 then S:=W else (
	IdL := inverse (id_(FR^d)-L);
    	S = (transpose IdL) * W * IdL
	);
    Sinv := inverse S;

    -----------------------------------------------------
    -- Compute score equations ideal
    ----------------------------------------------------
    -- Sample covariance matrix
    if opts.SampleData then V:= sampleCovarianceMatrix(U) else V=U;
    if ring V===RR_53 then V = matrix apply(entries V,r->r/(v->lift(numeric(opts.RealPrecision,v),QQ)));
    -- Jacobian of log-likelihood function
    C1 := trace(Sinv * V);
    C1derivative := jacobianMatrixOfRationalFunction(C1);
    LL :=jacobianMatrixOfRationalFunction (det Sinv)*matrix{{1/det(Sinv)}} - C1derivative;
    LL=flatten entries(LL);
    denoms := apply(#LL, i -> lift(denominator(LL_i), lpR));
    J:=ideal apply(#LL, i -> lift(numerator(LL_i),lpR));
    --Saturate
    if opts.Saturate then (
        argSaturate:=opts.SaturateOptions  >>newOpts-> args ->(args, newOpts);
    	for i from 0 to (#denoms-1) do (
	    if degree denoms_i =={0} then J=J else
	    J=saturate(argSaturate(J,denoms_i))
	    );
	);
    return (J,Sinv);
);

----------------------------------------------------
--scoreEquationsInternalUndir for undirected graphs
----------------------------------------------------
scoreEquationsInternalUndir={Saturate => true, SaturateOptions => options saturate, SampleData=>true, RealPrecision=> 53, CovarianceMatrix => false}>>opts->(R,U)->(
    -- Sample covariance matrix
    if opts.SampleData then V := sampleCovarianceMatrix(U) else V=U;
    if ring V===RR_53 then V = matrix apply(entries V,r->r/(v->lift(numeric(opts.RealPrecision,v),QQ)));
    -- Concentration matrix K
    K:=undirectedEdgesMatrix R;
    -- move to a new ring, lpR, which does not have the s variables
    lpR:=coefficientRing(R)[gens R - set support covarianceMatrix R];
    K=sub(K,lpR);
    J:=ideal{jacobian ideal{determinant(K)}-determinant(K)*jacobian(ideal{trace(K*V)})};
    if opts.Saturate then
    (  	argSaturate:=opts.SaturateOptions  >>newOpts-> args ->(args, newOpts);
    	J=saturate(argSaturate(J,ideal{determinant(K)}));
	);
    return (J,K);
 );


-------------------------------------------
-- Method copied from package DeterminantalRepresentations
-------------------------------------------
-----------------------------------------------

-- Method for retrieving the real part of a matrix. 
--The code of this function is directly taken from DeterminantalRepresentations package in M2.

realPartMatrix = A -> matrix apply(entries A, r -> r/realPart)

--**************************--
--  METHODS 	      	   	  --
--**************************--
sampleCovarianceMatrix = method(TypicalValue =>Matrix);
sampleCovarianceMatrix(Matrix) := (U) -> (
   n := numRows U;
   --Convert from integers to rationals if needed
   if ring U===ZZ then U=sub(U,QQ);
   --Convert matrix into list of row matrices
   U = for i to n-1 list U^{i};
   --Compute the mean vector
   Ubar := matrix{{(1/n)}} * sum(U);
   --Compute sample covariance matrix
   return ((1/n)*(sum apply(n, i -> (transpose (U#i-Ubar))*(U#i-Ubar))));
);

sampleCovarianceMatrix(List) := (U) -> (
   return sampleCovarianceMatrix(matrix U);
);

jacobianMatrixOfRationalFunction = method(TypicalValue =>Matrix);
jacobianMatrixOfRationalFunction(RingElement) := (F) -> (
    if not instance(ring F,FractionField) then error "Expected element in a field of fractions";
    f:=numerator(F);
    g:=denominator(F);
    R:=ring(f);
    answer:=diff(vars(R), f) * g - diff(vars(R), g)*f;
    answer=substitute(answer, ring(F));
    return transpose(matrix({{(1/g)^2}})*answer)
);

scoreEquations = method(TypicalValue =>Sequence, Options =>{SampleData => true, Saturate => true, SaturateOptions => options saturate, RealPrecision => 53, CovarianceMatrix => false});
scoreEquations(Ring,Matrix) := opts -> (R, U) -> (
    ----------------------------------------------------
    --Check input
    ----------------------------------------------------
    if not R.?graph then error "Expected a ring created with gaussianRing of a Graph, Bigraph, Digraph or MixedGraph";
    if not numRows U==#vertices R.graph  then error "Size of sample data does not match the graph.";
    if not opts.SampleData then (if not U==transpose U then error "The sample covariance matrix must be symmetric.");
    ---------------------------------------------------
    -- Apply appropriate scoreEquations routine
    ---------------------------------------------------
    if R.graphType===Graph
    then (J,Sinv):=scoreEquationsInternalUndir(R,U,opts)
    else (J,Sinv)=scoreEquationsInternal(R,U,opts);

    if not opts.CovarianceMatrix
    then return J
    else return (J, inverse Sinv)
    ;
);

scoreEquations(Matrix,Ring) := opts ->(U,R) -> (
    return scoreEquations(R,U,opts);
    );

scoreEquations(Ring,List) := opts ->(R, U) -> (
    ----------------------------------------------------
    --Check input
    ----------------------------------------------------
    if not opts.SampleData then error "The sample covariance matrix must be a matrix.";

    ---------------------------------------------------
    -- Call scoreEquations routine with a matrix
    ---------------------------------------------------
    return scoreEquations(R,matrix U,opts);
);

scoreEquations(List,Ring) := opts ->(U,R) -> (
    return scoreEquations(R,U,opts);
    );

--Allow for graphs as input 
scoreEquations(MixedGraph, Matrix) := opts ->(G,U) -> (
    return scoreEquations(gaussianRing(G),U,opts);
    );

scoreEquations(MixedGraph,List) := opts ->(G,U) -> (
    return scoreEquations(gaussianRing(G),U,opts);
    );

--All permutations of input for MixedGraphs and other type of graphs
scoreEquations(Graph,List) := opts -> (G, U) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );

scoreEquations(Digraph,List) := opts -> (G, U) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );
scoreEquations(Bigraph,List) := opts -> (G, U) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );
scoreEquations(Graph,Digraph,List) := opts -> (G,D,U) -> (
    return scoreEquations(mixedGraph (G,D),U, opts);
    );
scoreEquations(Digraph,Graph,List) := opts -> (D,G,U) -> (
    return scoreEquations(mixedGraph (D,G),U, opts);
    );
scoreEquations(Digraph,Bigraph,List) := opts -> (D,B,U) -> (
    return scoreEquations(mixedGraph (D,B),U, opts);
    );
scoreEquations(Bigraph,Digraph,List) := opts -> (B,D,U) -> (
    return scoreEquations(mixedGraph (B,D),U, opts);
    );
scoreEquations(Graph, Bigraph,List) := opts -> (G,B,U) -> (
    return scoreEquations(mixedGraph (G,B),U, opts);
    );
scoreEquations(Bigraph,Graph,List) := opts -> (B,G,U) -> (
    return scoreEquations(mixedGraph (B,G),U, opts);
    );
scoreEquations(Graph, Digraph, Bigraph, List) := opts -> (G,D,B,U) -> (
    return scoreEquations(mixedGraph (G,D,B),U, opts);
    );
scoreEquations(Digraph, Bigraph, Graph, List) := opts -> (D,B,G,U) -> (
    return scoreEquations(mixedGraph (D,B,G),U, opts);
    );
scoreEquations(Bigraph, Graph, Digraph, List) := opts -> (B,G,D,U) -> (
    return scoreEquations(mixedGraph (B,G,D),U, opts);
    );
scoreEquations(Graph,Bigraph, Digraph, List) := opts -> (G,B,D,U) -> (
    return scoreEquations(mixedGraph (G,B,D),U, opts);
    );
scoreEquations(Bigraph, Digraph,Graph, List) := opts -> (B,D,G,U) -> (
    return scoreEquations(mixedGraph (B,D,G),U, opts);
    );
scoreEquations(Digraph, Graph, Bigraph, List) := opts -> (D,G,B,U) -> (
    return scoreEquations(mixedGraph (D,G,B),U, opts);
    );

scoreEquations(List,MixedGraph) := opts -> (U,G) -> (
    return scoreEquations(G,U, opts);
    );

scoreEquations(List,Graph) := opts -> (U,G) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );

scoreEquations(List,Digraph) := opts -> (U,G) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );
scoreEquations(List,Bigraph) := opts -> (U,G) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );
scoreEquations(List,Graph,Digraph) := opts -> (U,G,D) -> (
    return scoreEquations(mixedGraph (G,D),U, opts);
    );
scoreEquations(List,Digraph,Graph) := opts -> (U,D,G) -> (
    return scoreEquations(mixedGraph (D,G),U, opts);
    );
scoreEquations(List,Digraph,Bigraph) := opts -> (U,D,B) -> (
    return scoreEquations(mixedGraph (D,B),U, opts);
    );
scoreEquations(List,Bigraph,Digraph) := opts -> (U,B,D) -> (
    return scoreEquations(mixedGraph (B,D),U, opts);
    );
scoreEquations(List,Graph, Bigraph) := opts -> (U,G,B) -> (
    return scoreEquations(mixedGraph (G,B),U, opts);
    );
scoreEquations(List,Bigraph,Graph) := opts -> (U,B,G) -> (
    return scoreEquations(mixedGraph (B,G),U, opts);
    );
scoreEquations(List,Graph, Digraph, Bigraph) := opts -> (U,G,D,B) -> (
    return scoreEquations(mixedGraph (G,D,B),U, opts);
    );
scoreEquations(List,Digraph, Bigraph, Graph) := opts -> (U,D,B,G) -> (
    return scoreEquations(mixedGraph (D,B,G),U, opts);
    );
scoreEquations(List,Bigraph, Graph, Digraph) := opts -> (U,B,G,D) -> (
    return scoreEquations(mixedGraph (B,G,D),U, opts);
    );
scoreEquations(List,Graph,Bigraph, Digraph) := opts -> (U,G,B,D) -> (
    return scoreEquations(mixedGraph (G,B,D),U, opts);
    );
scoreEquations(List,Bigraph, Digraph,Graph) := opts -> (U,B,D,G) -> (
    return scoreEquations(mixedGraph (B,D,G),U, opts);
    );
scoreEquations(List,Digraph, Graph, Bigraph) := opts -> (U,D,G,B) -> (
    return scoreEquations(mixedGraph (D,G,B),U, opts);
    );

scoreEquations(Graph,Matrix) := opts -> (G, U) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );

scoreEquations(Digraph,Matrix) := opts -> (G, U) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );
scoreEquations(Bigraph,Matrix) := opts -> (G, U) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );
scoreEquations(Graph,Digraph,Matrix) := opts -> (G,D,U) -> (
    return scoreEquations(mixedGraph (G,D),U, opts);
    );
scoreEquations(Digraph,Graph,Matrix) := opts -> (D,G,U) -> (
    return scoreEquations(mixedGraph (D,G),U, opts);
    );
scoreEquations(Digraph,Bigraph,Matrix) := opts -> (D,B,U) -> (
    return scoreEquations(mixedGraph (D,B),U, opts);
    );
scoreEquations(Bigraph,Digraph,Matrix) := opts -> (B,D,U) -> (
    return scoreEquations(mixedGraph (B,D),U, opts);
    );
scoreEquations(Graph, Bigraph,Matrix) := opts -> (G,B,U) -> (
    return scoreEquations(mixedGraph (G,B),U, opts);
    );
scoreEquations(Bigraph,Graph,Matrix) := opts -> (B,G,U) -> (
    return scoreEquations(mixedGraph (B,G),U, opts);
    );
scoreEquations(Graph, Digraph, Bigraph, Matrix) := opts -> (G,D,B,U) -> (
    return scoreEquations(mixedGraph (G,D,B),U, opts);
    );
scoreEquations(Digraph, Bigraph, Graph, Matrix) := opts -> (D,B,G,U) -> (
    return scoreEquations(mixedGraph (D,B,G),U, opts);
    );
scoreEquations(Bigraph, Graph, Digraph, Matrix) := opts -> (B,G,D,U) -> (
    return scoreEquations(mixedGraph (B,G,D),U, opts);
    );
scoreEquations(Graph,Bigraph, Digraph, Matrix) := opts -> (G,B,D,U) -> (
    return scoreEquations(mixedGraph (G,B,D),U, opts);
    );
scoreEquations(Bigraph, Digraph,Graph, Matrix) := opts -> (B,D,G,U) -> (
    return scoreEquations(mixedGraph (B,D,G),U, opts);
    );
scoreEquations(Digraph, Graph, Bigraph, Matrix) := opts -> (D,G,B,U) -> (
    return scoreEquations(mixedGraph (D,G,B),U, opts);
    );

scoreEquations(Matrix,MixedGraph) := opts -> (U,G) -> (
    return scoreEquations(G,U, opts);
    );

scoreEquations(Matrix,Graph) := opts -> (U,G) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );

scoreEquations(Matrix,Digraph) := opts -> (U,G) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );
scoreEquations(Matrix,Bigraph) := opts -> (U,G) -> (
    return scoreEquations(mixedGraph (G),U, opts);
    );
scoreEquations(Matrix,Graph,Digraph) := opts -> (U,G,D) -> (
    return scoreEquations(mixedGraph (G,D),U, opts);
    );
scoreEquations(Matrix,Digraph,Graph) := opts -> (U,D,G) -> (
    return scoreEquations(mixedGraph (D,G),U, opts);
    );
scoreEquations(Matrix,Digraph,Bigraph) := opts -> (U,D,B) -> (
    return scoreEquations(mixedGraph (D,B),U, opts);
    );
scoreEquations(Matrix,Bigraph,Digraph) := opts -> (U,B,D) -> (
    return scoreEquations(mixedGraph (B,D),U, opts);
    );
scoreEquations(Matrix,Graph, Bigraph) := opts -> (U,G,B) -> (
    return scoreEquations(mixedGraph (G,B),U, opts);
    );
scoreEquations(Matrix,Bigraph,Graph) := opts -> (U,B,G) -> (
    return scoreEquations(mixedGraph (B,G),U, opts);
    );
scoreEquations(Matrix,Graph, Digraph, Bigraph) := opts -> (U,G,D,B) -> (
    return scoreEquations(mixedGraph (G,D,B),U, opts);
    );
scoreEquations(Matrix,Digraph, Bigraph, Graph) := opts -> (U,D,B,G) -> (
    return scoreEquations(mixedGraph (D,B,G),U, opts);
    );
scoreEquations(Matrix,Bigraph, Graph, Digraph) := opts -> (U,B,G,D) -> (
    return scoreEquations(mixedGraph (B,G,D),U, opts);
    );
scoreEquations(Matrix,Graph,Bigraph, Digraph) := opts -> (U,G,B,D) -> (
    return scoreEquations(mixedGraph (G,B,D),U, opts);
    );
scoreEquations(Matrix,Bigraph, Digraph,Graph) := opts -> (U,B,D,G) -> (
    return scoreEquations(mixedGraph (B,D,G),U, opts);
    );
scoreEquations(Matrix,Digraph, Graph, Bigraph) := opts -> (U,D,G,B) -> (
    return scoreEquations(mixedGraph (D,G,B),U, opts);
    );


checkPD = method(TypicalValue =>List, Options =>{ZeroTolerance=>1e-10});
checkPD(List) :=  opts -> (L) -> (
		for l in L
		list (
			if not length (select(eigenvalues l, i->realPart i<= opts.ZeroTolerance
				or abs(imaginaryPart i )>opts.ZeroTolerance))==0 then continue;
			realPartMatrix l)
);
checkPD(Matrix):= opts -> (L)->{
    return checkPD({L},opts);
};

checkPSD = method(TypicalValue =>List, Options =>{ZeroTolerance=>1e-10});
checkPSD(List) :=  opts -> (L) -> (
	for l in L
		list (
			if not length (select(eigenvalues l, i->realPart i< -opts.ZeroTolerance
				or abs(imaginaryPart i )>opts.ZeroTolerance))==0 then continue;
				realPartMatrix l)
);
checkPSD(Matrix):= opts -> (L)->{
    return checkPSD({L},opts);
};

MLdegree = method(TypicalValue =>ZZ);
MLdegree(Ring):= (R) -> (
   if not R.?graph then error "Expected gaussianRing created from a graph, digraph, bigraph or mixedGraph";
   n:=# vertices R.graph;
   J:=scoreEquations(R,random(RR^n,RR^n));
   dimJ := dim J;
   if dimJ > 0 then error concatenate("the ideal of score equations has dimension ",toString dimJ, " > 0,
       so ML degree is not well-defined. The degree of this ideal is ", toString degree J,".");
   return degree J;
);



solverMLE = method(TypicalValue =>Sequence, Options =>{SampleData=>true, ConcentrationMatrix=> false, Saturate => true, SaturateOptions => options saturate, Solver=>"EigenSolver", OptionsEigenSolver => options zeroDimSolve, OptionsNAG4M2=> options solveSystem, RealPrecision => 53, ZeroTolerance=>1e-10});
solverMLE(MixedGraph,Matrix) := opts -> (G, U) -> (    
    return  solverMLE(gaussianRing G,U,opts);
);

-- Allow list instead of matrix
solverMLE(MixedGraph,List):=  opts ->(G,U) -> (
    return  solverMLE(gaussianRing G,U,opts);
);

--Allow ring instead of graph
solverMLE(Ring,Matrix):= opts ->(R,U) -> (
        -- check input
    if not numgens source U ==R.gaussianRingData#nn then error "Size of sample data does not match the graph.";
    -- sample covariance matrix
    if opts.SampleData then V := sampleCovarianceMatrix(U)
    else (V=U;
    if not V==transpose V then error "The sample covariance matrix must be symmetric.");
    -- generate the ideal of the score equations
    if opts.Saturate then (
	 argSaturate:=opts.SaturateOptions  >>newOpts-> args ->(args, SaturateOptions=>newOpts,SampleData=>false);
         (J,SInv):=scoreEquationsInternal(argSaturate(R,V));)
    else (J,SInv)= scoreEquationsInternal(R,V,Saturate=>false, SampleData=>false);
    -- check that the system has finitely many solutions
    if dim J =!= 0 then return J
	else (
    ML:=degree J;
    -- solve system
    if opts.Solver=="EigenSolver" then(
	argES:=opts.OptionsEigenSolver  >>newOpts-> args ->(args, newOpts);
        sols:=zeroDimSolve(argES(J));
	) else (
	if opts.Solver=="NAG4M2" then (
	   sys:= flatten entries gens J;
	   argNAG4M2:=opts.OptionsNAG4M2 >>newOpts-> args ->(args, newOpts);
           sols=solveSystem(argNAG4M2(sys));
	    )
	else error "Accepted solver options are EigenSolver
	(which uses function zeroDimSolve) or NAG4M2 (which uses solveSystem). Options should
	be given as strings.";
	);

    --evaluate matrices on solutions
    M:=genListMatrix(sols,SInv);
    --consider only PD matrices
    L:=checkPD (M, ZeroTolerance=>opts.ZeroTolerance);
    --find the optimal points
    (maxPt, E):=maxMLE(L,V);
    if not opts.ConcentrationMatrix then (
	if instance(E,List) then E=(for e in E list e=inverse e) else  E=inverse E
	);
    return  (maxPt,E,ML));


);

-- Allow ring instead of graph and list instead of matrix
solverMLE(Ring,List):=  opts ->(R,U) -> (
    -- check input
    if not opts.SampleData then error "The sample covariance matrix must be a matrix.";
    -- call solverMLE for a matrix
    return  solverMLE(R,matrix U,opts);
);


-- Permutations of input

solverMLE(Matrix,Ring) := opts -> (U, R) -> (
    return solverMLE(R,U, opts);
    );
solverMLE(List,Ring) := opts -> (U, R) -> (
    return solverMLE(R,U, opts);
    );

solverMLE(Graph,List) := opts -> (G, U) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );

solverMLE(Digraph,List) := opts -> (G, U) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );
solverMLE(Bigraph,List) := opts -> (G, U) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );
solverMLE(Graph,Digraph,List) := opts -> (G,D,U) -> (
    return solverMLE(mixedGraph (G,D),U, opts);
    );
solverMLE(Digraph,Graph,List) := opts -> (D,G,U) -> (
    return solverMLE(mixedGraph (D,G),U, opts);
    );
solverMLE(Digraph,Bigraph,List) := opts -> (D,B,U) -> (
    return solverMLE(mixedGraph (D,B),U, opts);
    );
solverMLE(Bigraph,Digraph,List) := opts -> (B,D,U) -> (
    return solverMLE(mixedGraph (B,D),U, opts);
    );
solverMLE(Graph, Bigraph,List) := opts -> (G,B,U) -> (
    return solverMLE(mixedGraph (G,B),U, opts);
    );
solverMLE(Bigraph,Graph,List) := opts -> (B,G,U) -> (
    return solverMLE(mixedGraph (B,G),U, opts);
    );
solverMLE(Graph, Digraph, Bigraph, List) := opts -> (G,D,B,U) -> (
    return solverMLE(mixedGraph (G,D,B),U, opts);
    );
solverMLE(Digraph, Bigraph, Graph, List) := opts -> (D,B,G,U) -> (
    return solverMLE(mixedGraph (D,B,G),U, opts);
    );
solverMLE(Bigraph, Graph, Digraph, List) := opts -> (B,G,D,U) -> (
    return solverMLE(mixedGraph (B,G,D),U, opts);
    );
solverMLE(Graph,Bigraph, Digraph, List) := opts -> (G,B,D,U) -> (
    return solverMLE(mixedGraph (G,B,D),U, opts);
    );
solverMLE(Bigraph, Digraph,Graph, List) := opts -> (B,D,G,U) -> (
    return solverMLE(mixedGraph (B,D,G),U, opts);
    );
solverMLE(Digraph, Graph, Bigraph, List) := opts -> (D,G,B,U) -> (
    return solverMLE(mixedGraph (D,G,B),U, opts);
    );

solverMLE(List,MixedGraph) := opts -> (U,G) -> (
    return solverMLE(G,U, opts);
    );

solverMLE(List,Graph) := opts -> (U,G) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );

solverMLE(List,Digraph) := opts -> (U,G) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );
solverMLE(List,Bigraph) := opts -> (U,G) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );
solverMLE(List,Graph,Digraph) := opts -> (U,G,D) -> (
    return solverMLE(mixedGraph (G,D),U, opts);
    );
solverMLE(List,Digraph,Graph) := opts -> (U,D,G) -> (
    return solverMLE(mixedGraph (D,G),U, opts);
    );
solverMLE(List,Digraph,Bigraph) := opts -> (U,D,B) -> (
    return solverMLE(mixedGraph (D,B),U, opts);
    );
solverMLE(List,Bigraph,Digraph) := opts -> (U,B,D) -> (
    return solverMLE(mixedGraph (B,D),U, opts);
    );
solverMLE(List,Graph, Bigraph) := opts -> (U,G,B) -> (
    return solverMLE(mixedGraph (G,B),U, opts);
    );
solverMLE(List,Bigraph,Graph) := opts -> (U,B,G) -> (
    return solverMLE(mixedGraph (B,G),U, opts);
    );
solverMLE(List,Graph, Digraph, Bigraph) := opts -> (U,G,D,B) -> (
    return solverMLE(mixedGraph (G,D,B),U, opts);
    );
solverMLE(List,Digraph, Bigraph, Graph) := opts -> (U,D,B,G) -> (
    return solverMLE(mixedGraph (D,B,G),U, opts);
    );
solverMLE(List,Bigraph, Graph, Digraph) := opts -> (U,B,G,D) -> (
    return solverMLE(mixedGraph (B,G,D),U, opts);
    );
solverMLE(List,Graph,Bigraph, Digraph) := opts -> (U,G,B,D) -> (
    return solverMLE(mixedGraph (G,B,D),U, opts);
    );
solverMLE(List,Bigraph, Digraph,Graph) := opts -> (U,B,D,G) -> (
    return solverMLE(mixedGraph (B,D,G),U, opts);
    );
solverMLE(List,Digraph, Graph, Bigraph) := opts -> (U,D,G,B) -> (
    return solverMLE(mixedGraph (D,G,B),U, opts);
    );

solverMLE(Graph,Matrix) := opts -> (G, U) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );

solverMLE(Digraph,Matrix) := opts -> (G, U) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );
solverMLE(Bigraph,Matrix) := opts -> (G, U) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );
solverMLE(Graph,Digraph,Matrix) := opts -> (G,D,U) -> (
    return solverMLE(mixedGraph (G,D),U, opts);
    );
solverMLE(Digraph,Graph,Matrix) := opts -> (D,G,U) -> (
    return solverMLE(mixedGraph (D,G),U, opts);
    );
solverMLE(Digraph,Bigraph,Matrix) := opts -> (D,B,U) -> (
    return solverMLE(mixedGraph (D,B),U, opts);
    );
solverMLE(Bigraph,Digraph,Matrix) := opts -> (B,D,U) -> (
    return solverMLE(mixedGraph (B,D),U, opts);
    );
solverMLE(Graph, Bigraph,Matrix) := opts -> (G,B,U) -> (
    return solverMLE(mixedGraph (G,B),U, opts);
    );
solverMLE(Bigraph,Graph,Matrix) := opts -> (B,G,U) -> (
    return solverMLE(mixedGraph (B,G),U, opts);
    );
solverMLE(Graph, Digraph, Bigraph, Matrix) := opts -> (G,D,B,U) -> (
    return solverMLE(mixedGraph (G,D,B),U, opts);
    );
solverMLE(Digraph, Bigraph, Graph, Matrix) := opts -> (D,B,G,U) -> (
    return solverMLE(mixedGraph (D,B,G),U, opts);
    );
solverMLE(Bigraph, Graph, Digraph, Matrix) := opts -> (B,G,D,U) -> (
    return solverMLE(mixedGraph (B,G,D),U, opts);
    );
solverMLE(Graph,Bigraph, Digraph, Matrix) := opts -> (G,B,D,U) -> (
    return solverMLE(mixedGraph (G,B,D),U, opts);
    );
solverMLE(Bigraph, Digraph,Graph, Matrix) := opts -> (B,D,G,U) -> (
    return solverMLE(mixedGraph (B,D,G),U, opts);
    );
solverMLE(Digraph, Graph, Bigraph, Matrix) := opts -> (D,G,B,U) -> (
    return solverMLE(mixedGraph (D,G,B),U, opts);
    );

solverMLE(Matrix,MixedGraph) := opts -> (U,G) -> (
    return solverMLE(G,U, opts);
    );

solverMLE(Matrix,Graph) := opts -> (U,G) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );

solverMLE(Matrix,Digraph) := opts -> (U,G) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );
solverMLE(Matrix,Bigraph) := opts -> (U,G) -> (
    return solverMLE(mixedGraph (G),U, opts);
    );
solverMLE(Matrix,Graph,Digraph) := opts -> (U,G,D) -> (
    return solverMLE(mixedGraph (G,D),U, opts);
    );
solverMLE(Matrix,Digraph,Graph) := opts -> (U,D,G) -> (
    return solverMLE(mixedGraph (D,G),U, opts);
    );
solverMLE(Matrix,Digraph,Bigraph) := opts -> (U,D,B) -> (
    return solverMLE(mixedGraph (D,B),U, opts);
    );
solverMLE(Matrix,Bigraph,Digraph) := opts -> (U,B,D) -> (
    return solverMLE(mixedGraph (B,D),U, opts);
    );
solverMLE(Matrix,Graph, Bigraph) := opts -> (U,G,B) -> (
    return solverMLE(mixedGraph (G,B),U, opts);
    );
solverMLE(Matrix,Bigraph,Graph) := opts -> (U,B,G) -> (
    return solverMLE(mixedGraph (B,G),U, opts);
    );
solverMLE(Matrix,Graph, Digraph, Bigraph) := opts -> (U,G,D,B) -> (
    return solverMLE(mixedGraph (G,D,B),U, opts);
    );
solverMLE(Matrix,Digraph, Bigraph, Graph) := opts -> (U,D,B,G) -> (
    return solverMLE(mixedGraph (D,B,G),U, opts);
    );
solverMLE(Matrix,Bigraph, Graph, Digraph) := opts -> (U,B,G,D) -> (
    return solverMLE(mixedGraph (B,G,D),U, opts);
    );
solverMLE(Matrix,Graph,Bigraph, Digraph) := opts -> (U,G,B,D) -> (
    return solverMLE(mixedGraph (G,B,D),U, opts);
    );
solverMLE(Matrix,Bigraph, Digraph,Graph) := opts -> (U,B,D,G) -> (
    return solverMLE(mixedGraph (B,D,G),U, opts);
    );
solverMLE(Matrix,Digraph, Graph, Bigraph) := opts -> (U,D,G,B) -> (
    return solverMLE(mixedGraph (D,G,B),U, opts);
    );


--******************************************--
-- DOCUMENTATION     	       	    	    --
--******************************************--

beginDocumentation()

doc ///
    Key
        GraphicalModelsMLE
    Headline
        a package for MLE of parameters for Gaussian graphical models
    Description
     Text
      {\bf Graphical Models MLE} is a package for algebraic statistics that broadens the functionalities of @TO GraphicalModels@.
      It computes the maximum likelihood estimates (MLE) of the covariance matrix of Gaussian graphical models associated to loopless mixed graphs(LMG).

      The main features of the package are the computation of the @TO sampleCovarianceMatrix@ of sample data,
      the ideal generated by @TO scoreEquations@ of log-likelihood functions of Gaussian graphical model,
      the @TO MLdegree@ of such models and the MLE for the covariance or concentration matrix via @TO solverMLE@.

      For more details on the type of graphical models that are accepted see @TO gaussianRing@.
      In particular, for further information about LMG with undirected, directed and bidirected edges, check @TO partitionLMG@.


      {\bf References:}

      An introduction to key notions such as MLE and ML-degree can be found in the books:

      Seth Sullivant, {\em Algebraic statistics}, American Mathematical Society, Vol 194, 2018.

      Mathias Drton, Bernd Sturmfels and Seth Sullivant, {\em Lectures on Algebraic Statistics}, Oberwolfach Seminars, Vol 40, Birkhauser, Basel, 2009.

      The definition and classification of loopless mixed graphs (LMG) can be found in the paper:

      Kayvan Sadeghi and Steffen Lauritzen, {\em Markov properties for mixed graphs}, Bernoulli, 20 (2014), no 2, 676-696.


      {\bf Examples:}

      Computation of a sample covariance matrix from sample data:
     Example
      U= matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
      sampleCovarianceMatrix(U)
     Text
      The ideal generated by the score equations of the log-likelihood function of the graphical model associated to the
      graph $1\rightarrow 2,1\rightarrow 3,2\rightarrow 3,3\rightarrow 4,3<-> 4$ is computed as follows:
     Example
      G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}});
      R = gaussianRing(G);
      U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
      scoreEquations(R,U)
     Text
      Computation of the ML-degree of the 4-cycle:
     Example
      G=graph{{1,2},{2,3},{3,4},{4,1}};
      MLdegree(gaussianRing G)
     Text
      Next compute the MLE for the covariance matrix of the graphical model associated
      to the graph $1\rightarrow 3,2\rightarrow 4,3<-> 4,1 - 2$.
      The input is the sample covariance instead of the sample data.
     Example
      G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}},graph {{1,2}});
      V =  matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}};
      solverMLE(G,V,SampleData=>false)

     Text
      As an application of @TO solverMLE@: positive definite matrix completion
     Text
      Consider the following symmetric matrix with some unknown entries:
     Example
      R=QQ[x,y];
      M=matrix{{115,-13,x,47},{-13,5,7,y},{x,7,27,-21},{47,y,-21,29}}
     Text
      Unknown entries correspond to the non-edges of the 4-cycle. A positive definite completion of this matrix
      is obtained by giving values to x and y and computing the MLE for the covariance matrix in the Gaussian graphical model
      given by the 4-cycle. Check @TO solverMLE@ for more details.
     Example
      G=graph{{1,2},{2,3},{3,4},{1,4}};
      V=matrix{{115,-13,-29,47},{-13,5,7,-11},{-29,7,27,-21},{47,-11,-21,29}};
      (mx,MLE,ML)=solverMLE(G,V,SampleData=>false)

    Caveat
     GraphicalModelsMLE requires @TO Graphs@,  @TO StatGraphs@ and  @TO GraphicalModels@.
     In order to use the default numerical solver, it also requires @TO EigenSolver@.

     @TO Graphs@ allows the user to create graphs whose vertices are labeled arbitrarily.
     However, several functions in GraphicalModels sort the vertices of the graph. Hence, graphs used as input to methods
     in GraphicalModelsMLE must have sortable vertex labels, e.g., all numbers or all letters.

     @TO StatGraphs@ allows the user to work with objects such as bigraphs and mixedGraphs.

     @TO GraphicalModels@ is used to generate @TO gaussianRing@, i.e. rings encoding graph properties.


///

--------------------------------
-- Documentation
--------------------------------

doc ///
    Key
        sampleCovarianceMatrix
        (sampleCovarianceMatrix, List)
	(sampleCovarianceMatrix, Matrix)
    Headline
        sample covariance matrix of observation vectors
    Usage
        sampleCovarianceMatrix U
    Inputs
	U:Matrix
	   or @TO List@ of sample data
    Outputs
         :Matrix
           sample covariance matrix of the sample data
    Description
        Text
	    The sample covariance matrix is $S = \frac{1}{n} \sum_{i=1}^{n} (X^{(i)}-\bar{X}) (X^{(i)}-\bar{X})^T$.
	    Note that for normally distributed random variables, $S$ is the maximum likelihood estimator (MLE) for the
	    covariance matrix. This is different from the unbiased estimator, which uses a denominator of $n-1$ instead of $n$.

	    Sample data is input as a matrix or a list.
	    The rows of the matrix or the elements of the list are observation vectors.
	Example
          L= {{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
	  sampleCovarianceMatrix(L)
	  U= matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
	  sampleCovarianceMatrix(U)

     ///

doc ///
    Key
        jacobianMatrixOfRationalFunction
        (jacobianMatrixOfRationalFunction,RingElement)
    Headline
        Jacobian matrix of a rational function
    Usage
        jacobianMatrixOfRationalFunction(F)
    Inputs
        F:RingElement
           in @TO frac@
    Outputs
         :Matrix
	   the Jacobian matrix of a rational function
    Description
        Text
	    This function computes the Jacobian matrix of a rational function.
	    The input is an element in a fraction field.
        Example
	    R=QQ[x,y];
	    FR=frac R;
	    F=1/(x^2+y^2);
            jacobianMatrixOfRationalFunction(F)
        Example
	    R=QQ[t_1,t_2,t_3];
	    FR=frac R;
	    jacobianMatrixOfRationalFunction( (t_1^2*t_2)/(t_1+t_2^2+t_3^3) )
   ///

-------------------------------------------------------
-- Documentation scoreEquations -----------------------
-------------------------------------------------------

doc ///
    Key
        scoreEquations
        (scoreEquations, Ring, List)
	(scoreEquations, Ring, Matrix)
	(scoreEquations, List, Ring)
	(scoreEquations, Matrix, Ring)	
	(scoreEquations, MixedGraph, Matrix)
	(scoreEquations, MixedGraph, List)
	(scoreEquations, Graph, List)
	(scoreEquations, Digraph, List)
	(scoreEquations, Bigraph, List)
	(scoreEquations, Graph, Digraph,List)
	(scoreEquations, Digraph, Graph,List)
	(scoreEquations, Digraph, Bigraph, List)
	(scoreEquations, Bigraph, Digraph, List)
	(scoreEquations, Graph, Bigraph, List)
	(scoreEquations, Bigraph, Graph, List)
	(scoreEquations, Graph, Digraph, Bigraph, List)
	(scoreEquations, Digraph, Bigraph, Graph, List)
	(scoreEquations, Bigraph, Graph, Digraph, List)
	(scoreEquations, Graph, Bigraph, Digraph, List)
	(scoreEquations, Bigraph, Digraph, Graph, List)
	(scoreEquations, Digraph, Graph, Bigraph, List)
	(scoreEquations, List, MixedGraph)
	(scoreEquations, List, Graph)
	(scoreEquations, List, Digraph)
	(scoreEquations, List, Bigraph)
	(scoreEquations, List, Graph, Digraph)
	(scoreEquations, List, Digraph,Graph)
	(scoreEquations, List, Digraph,Bigraph)
	(scoreEquations, List, Bigraph,Digraph)
	(scoreEquations, List, Graph, Bigraph)
	(scoreEquations, List, Bigraph, Graph)
	(scoreEquations, List, Graph, Digraph, Bigraph)
	(scoreEquations, List, Digraph, Bigraph, Graph)
	(scoreEquations, List, Bigraph, Graph, Digraph)
	(scoreEquations, List, Graph, Bigraph, Digraph)
	(scoreEquations, List, Bigraph, Digraph, Graph)
	(scoreEquations, List, Digraph, Graph, Bigraph)
	(scoreEquations, Graph, Matrix)
	(scoreEquations, Digraph, Matrix)
	(scoreEquations, Bigraph, Matrix)
	(scoreEquations, Graph, Digraph,Matrix)
	(scoreEquations, Digraph, Graph,Matrix)
	(scoreEquations, Digraph, Bigraph, Matrix)
	(scoreEquations, Bigraph, Digraph, Matrix)
	(scoreEquations, Graph, Bigraph, Matrix)
	(scoreEquations, Bigraph, Graph, Matrix)
	(scoreEquations, Graph, Digraph, Bigraph, Matrix)
	(scoreEquations, Digraph, Bigraph, Graph, Matrix)
	(scoreEquations, Bigraph, Graph, Digraph, Matrix)
	(scoreEquations, Graph, Bigraph, Digraph, Matrix)
	(scoreEquations, Bigraph, Digraph, Graph, Matrix)
	(scoreEquations, Digraph, Graph, Bigraph, Matrix)
	(scoreEquations, Matrix, MixedGraph)
	(scoreEquations, Matrix, Graph)
	(scoreEquations, Matrix, Digraph)
	(scoreEquations, Matrix, Bigraph)
	(scoreEquations, Matrix, Graph, Digraph)
	(scoreEquations, Matrix, Digraph,Graph)
	(scoreEquations, Matrix, Digraph,Bigraph)
	(scoreEquations, Matrix, Bigraph,Digraph)
	(scoreEquations, Matrix, Graph, Bigraph)
	(scoreEquations, Matrix, Bigraph, Graph)
	(scoreEquations, Matrix, Graph, Digraph, Bigraph)
	(scoreEquations, Matrix, Digraph, Bigraph, Graph)
	(scoreEquations, Matrix, Bigraph, Graph, Digraph)
	(scoreEquations, Matrix, Graph, Bigraph, Digraph)
	(scoreEquations, Matrix, Bigraph, Digraph, Graph)
	(scoreEquations, Matrix, Digraph, Graph, Bigraph)
    Headline
        score equations of the log-likelihood function of a Gaussian graphical model
    Usage
        scoreEquations(R,U)
    Inputs
        R:Ring
	   defined as a @TO gaussianRing@ of @TO Graph@, or @TO Digraph@, or @TO Bigraph@, or @TO MixedGraph@

	U:Matrix
	   or @TO List@ of sample data.
	   Alternatively, the input can be the sample covariance @TO Matrix@ by setting the optional input @TO SampleData@ to false
    Outputs
         :Sequence
	   consisting either of (Ideal) or (Ideal,Matrix)
	   where the ideal is generated by the score equations of the log-likelihood
	   function of the Gaussian model and the matrix is the symbolic covariance
	   matrix of the model
    Description
        Text
	    This function computes the score equations that arise from taking
	    partial derivatives of the log-likelihood function of the concentration matrix
	    (the inverse of the covariance matrix) of a Gaussian graphical
	    statistical model and returns the ideal generated by such equations.

	    The input of this function is a @TO gaussianRing@ and statistical data.
	    The latter can be given as a matrix or a list of observations. The rows of the matrix or the elements of the list are observation vectors given as lists.
            It is possible to input the sample covariance matrix directly by using the optional input @TO SampleData@.
        Example
	    G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}});
	    R = gaussianRing(G);
	    U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
	    JU=scoreEquations(R,U)
	    V = sampleCovarianceMatrix U
	    JV=scoreEquations(R,V,SampleData=>false)
        Text
	    @TO SaturateOptions@ allows to use all functionalities of @TO saturate@.
	    @TO Saturate@ determines whether to saturate. Note that the latter will not
	    provide the score equations of the model.

        Example
	    G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}});
	    R = gaussianRing(G);
            U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
            J=scoreEquations(R,U,SaturateOptions => {Strategy => Eliminate})
            JnoSat=scoreEquations(R,U,Saturate=>false)
	Text
	   The ML-degree of the model is the degree of the score equations ideal. The ML-degree
	   of the running example is 1:

	Example
	    G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}});
	    R = gaussianRing(G);
	    U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
	    J = scoreEquations(R,U)
	    dim J, degree J

    ///

doc ///
  Key
    CovarianceMatrix
  Headline
    optional input to output covariance matrix
  SeeAlso
     scoreEquations
   ///

doc ///
  Key
    [scoreEquations,CovarianceMatrix]
  Headline
    output covariance matrix
  Usage
    scoreEquations(R,U,CovarianceMatrix=>b)
  Inputs 
     b:Boolean
        default is false
  Description
    Text
     @TO [scoreEquations,CovarianceMatrix]@ is set to false by default. If b is true,
     @TO scoreEquations@ gives an additional output:
     the covariance matrix with rational entries in the same variables as the ideal of score equations.

    Example
     G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
     R=gaussianRing(G);
     U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
     (J,Sigma)=scoreEquations(R,U,CovarianceMatrix=>true);
     Sigma

  SeeAlso
     scoreEquations
///



doc ///
  Key
    Saturate
  Headline
    optional input whether to saturate
  SeeAlso
     scoreEquations
     solverMLE
   ///

doc ///
  Key
    [scoreEquations, Saturate]
  Headline
     whether to saturate
  Usage
    scoreEquations(R,U,Saturate=>b)
  Inputs 
     b:Boolean
        default is true
  Description
    Text
     @TO [scoreEquations,Saturate]@ is set to true by default. If b is false, saturation is not performed.

     Avoiding saturation is only intended for big computations
     when saturation cannot be computed or the computational time is very high.
     When @TO Saturate@ is set to false, @TO scoreEquations@ might not output the ideal
     generated by score equations because of the existence of vanishing denominators.

     For graphs with only undirected edges, the ideal of score equations is the saturation of the
     outputted ideal by the determinant of the concentration matrix. In the general case,
     the ideal of score equations consist of the saturation of the outputted ideal by the product of denominators
     of the Jacobian matrix.

     For example, in the following case the degree of the ideal prior to saturation already gives the right ML-degree:
    Example
     G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
     R=gaussianRing(G);
     U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
     JnoSat=scoreEquations(R,U,Saturate=>false);
     dim JnoSat
     degree JnoSat
     J=scoreEquations(R,U)
     degree JnoSat==degree J

  SeeAlso
     scoreEquations
///

doc ///
  Key
    [solverMLE, Saturate]
  Headline
    whether to saturate
  Usage
    solverMLE(G,U,Saturate=>b)  
  Inputs  
     b:Boolean
        default is true
  Description
    Text

     @TO [solverMLE,Saturate]@ is set to true by default.
     If we set @TO Saturate@ to false in @TO solverMLE@, saturation will not be
     performed when computing the score equations of the log-likelihood function,
     see @TO [scoreEquations,Saturate]@.

     If the ideal returned by @TO scoreEquations@ has positive dimension, @TO solverMLE@
     gives this ideal as output.

     On the other hand, if we obtain a zero-dimensional ideal in @TO scoreEquations@,
     @TO solverMLE@ computes the solutions of this polynomial system and returns:

     * the maximal value of the log-likelihood function attained by positive definite matrices
     corresponding to such solutions,

     * the positive definite matrices where the maximum is attained,

     * the degree of the zero-dimensional ideal.

     Be aware that this output might not correspond to the right MLE.

    Example
     G=graph{{1,2},{2,3},{3,4},{1,4}}
     U=random(ZZ^4,ZZ^4)
     solverMLE(G,U,Saturate=>false)

  SeeAlso
     solverMLE
     scoreEquations
     [scoreEquations,Saturate]
///


doc ///
  Key
    SaturateOptions
  Headline
    optional input to use options "saturate"
  SeeAlso
     scoreEquations
     solverMLE
     Saturate
     saturate
   ///
doc ///
  Key
    [scoreEquations, SaturateOptions]
  Headline
    use options from "saturate"
  Usage
    scoreEquations(R,U,SaturateOptions=>L)
  Inputs 
    L: List
       of options to set up saturation. Accepts any option from the function
       @TO saturate@
  Description
    Text
     Default @TO SaturateOptions@ in @TO scoreEquations@ are the
     default options in @TO saturate@. All optional input in @TO saturate@ is allowed.

    Example
     G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}})
     R=gaussianRing(G)
     U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
     J=scoreEquations(R,U,SaturateOptions => {DegreeLimit=>1, MinimalGenerators => false})

  SeeAlso
     scoreEquations
     Saturate
     saturate
     solverMLE
///

doc ///
  Key
    [solverMLE, SaturateOptions]
  Headline
    use options from "saturate"
  Usage
    solverMLE(G,U,SaturateOptions=>L)  
  Inputs 
    L: List
       of options to set up saturation. Accepts any option from the function
       @TO saturate@

  Description
    Text
     Default @TO SaturateOptions@ in @TO solverMLE@ are the default options in @TO saturate@.
     All optional input in @TO saturate@ is allowed.

    Example
     G=graph{{1,2},{2,3},{3,4},{1,4}}
     U=random(ZZ^4,ZZ^4)
     solverMLE(G,U,SaturateOptions => {DegreeLimit=>1, MinimalGenerators => false})

  SeeAlso
     scoreEquations
     Saturate
     saturate
     solverMLE
///

doc ///
  Key
    Solver
  Headline
    optional input to choose numerical solver
  SeeAlso
    solverMLE
    EigenSolver
    NumericalAlgebraicGeometry
    zeroDimSolve
    solveSystem
   ///
doc ///
  Key
    [solverMLE, Solver]
  Headline
    choose numerical solver
  Usage
    solverMLE(G,U,Solver=>P)  
  Inputs 
    P: String
       name of the corresponding package

  Description
    Text
      This option allows to choose which numerical solver to use to estimate the critical
      points. There are two options: "EigenSolver" or "NAG4M2" (Numerical Algebraic Geometry for Macaulay2).

      The default and strongly recommended option is "EigenSolver", in which case
      the function @TO zeroDimSolve@ is used.
      If "NAG4M2" is chosen, then @TO solveSystem@ is used.
    Example
     G=mixedGraph(graph{{a,b}},digraph{{a,d}},bigraph{{c,d}})
     U=matrix{{1, 2, 5, 1}, {5, 3, 2, 1}, {4, 3, 5, 10}, {2, 5,1, 3}};
     solverMLE (G,U,Solver=>"EigenSolver")
     solverMLE (G,U,Solver=>"NAG4M2")
  SeeAlso
     solverMLE
     EigenSolver
     NumericalAlgebraicGeometry
     zeroDimSolve
     solveSystem
///

doc ///
  Key
    OptionsEigenSolver
  Headline
    optional input to use options of "zeroDimSolve" in "EigenSolver"
  SeeAlso
    solverMLE
    EigenSolver
    zeroDimSolve
   ///
doc ///
  Key
    [solverMLE, OptionsEigenSolver]
  Headline
    use options of "zeroDimSolve" in "EigenSolver"
  Usage
    solverMLE(G,U,Solver=>"EigenSolver",OptionsEigenSolver=>L)
  Inputs 
    L: List
       of optional inputs in @TO zeroDimSolve@

  Description
    Text
     Default @TO OptionsEigenSolver@ in @TO solverMLE@ are the default options in @TO zeroDimSolve@ in @TO EigenSolver@.
     All optional input in @TO zeroDimSolve@ is allowed.
    Example
     G=mixedGraph(graph{{a,b},{b,c}})
     solverMLE(G,matrix{{1,0,0},{0,1,0},{0,0,1}},Solver=>"EigenSolver",OptionsEigenSolver=>{Multiplier =>1, Strategy=>"Stickelberger"})
    Text
     In fact, since @TO zeroDimSolve@ is the current default numerical solver.
     the sintaxis below is enough:
    Example
     G=mixedGraph(graph{{a,b},{b,c}})
     solverMLE(G,matrix{{1,0,0},{0,1,0},{0,0,1}},OptionsEigenSolver=>{Multiplier =>1, Strategy=>"Stickelberger"})

  SeeAlso
     solverMLE
     EigenSolver
     zeroDimSolve
///

doc ///
  Key
   OptionsNAG4M2
  Headline
   optional parameter to set the parameters of solveSystem in NumericalAlgebraicGeometry
  SeeAlso
   solverMLE
   NumericalAlgebraicGeometry
   solveSystem
   ///

doc ///
  Key
    [solverMLE, OptionsNAG4M2]
  Headline
    use options of "solveSystem" in "NumericalAlgebraicGeometry"
  Usage
    solverMLE(G,U,Solver=>"NAG4M2",OptionsNAG4M2=>L)
  Inputs 
    L: List
       of optional inputs to @TO solveSystem@

  Description
    Text
     Default @TO OptionsNAG4M2@ in @TO solverMLE@ when setting @TO [solverMLE,Solver]@
     to "NAG4M2" are the default options of @TO solveSystem@ in @TO NumericalAlgebraicGeometry@.
     All optional input in @TO solveSystem@ is allowed.

    Example
     G=mixedGraph(graph{{a,b},{b,c}})
     solverMLE(G,matrix{{1,0,0},{0,1,0},{0,0,1}},Solver=>"NAG4M2",OptionsNAG4M2=>{tStep =>.01,numberSuccessesBeforeIncrease => 5})

  SeeAlso
     solverMLE
     NumericalAlgebraicGeometry
     solveSystem
///

doc ///
  Key
    RealPrecision
  Headline
    optional input to choose the number of decimals used to round to QQ when inputing data in RR
  SeeAlso
     scoreEquations
     solverMLE
   ///
doc ///
  Key
    [scoreEquations, RealPrecision]
  Headline
    the number of bits of precision to use in the computation
  Usage
    scoreEquations(R,U,RealPrecision=>n)
  Inputs 
    n: ZZ
        default is 53
  Description
    Text
     This optional input only applies when the sample data or the sample covariance matrix has real entries.
     By default, the precision is 53 (the default precision for real numbers in M2).
    Example
     G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
     R=gaussianRing(G);
     U = matrix{{6.2849049, 10.292875, 1.038475, 1.1845757}, {3.1938475, 3.2573, 1.13847, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
     J=scoreEquations(R,U,RealPrecision=>3)

  SeeAlso
     scoreEquations
///

doc ///
  Key
    [solverMLE, RealPrecision]
  Headline
    the number of bits of precision to use in the computation
  Usage
    solverMLE(G,U,RealPrecision=>n)
  Inputs 
    n: ZZ
        default is 53
  Description
    Text
     This optional input only applies when the sample data or the sample covariance matrix has real entries.
     By default, the precision is 53 (the default precision for real numbers in M2).
    Example
     G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
     U = matrix{{6.2849049, 10.292875, 1.038475, 1.1845757}, {3.1938475, 3.2573, 1.13847, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
     solverMLE(G,U,RealPrecision=>10)

  SeeAlso
     scoreEquations
     solverMLE
///

doc ///
  Key
    SampleData
  Headline
    optional input to allow to input the sample covariance matrix instead of sample data
  SeeAlso
     scoreEquations
     solverMLE
   ///
doc ///
  Key
    [scoreEquations, SampleData]
  Headline
    input sample covariance matrix instead of sample data
  Usage
    scoreEquations(R,U,SampleData=>b)
  Inputs 
    b: Boolean
        default is true
  Description
    Text
      @TO scoreEquations@ requires a matrix or a list of sample data as part of the default input.
      Setting @TO [scoreEquations,SampleData]@ to false allows the user to enter a sample
      covariance matrix instead of sample data. It must be a symmetric matrix.
    Example
     G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
     R=gaussianRing(G);
     U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}}
     J=scoreEquations(R,U,SampleData=>true)
     V=sampleCovarianceMatrix(U)
     I=scoreEquations(R,V,SampleData=>false)

  SeeAlso
     scoreEquations
///

doc ///
  Key
    [solverMLE, SampleData]
  Headline
    input sample covariance matrix instead of sample data
  Usage
    solverMLE(G,U,SampleData=>b)  
  Inputs 
    b: Boolean
        default is true
  Description
    Text
      @TO solverMLE@ requires a matrix or a list of sample data as part of the default input.
      Setting @TO [solverMLE,SampleData]@ to false allows the user to enter a sample
      covariance matrix as input. It must be a symmetric matrix.
    Example
     G=graph{{1,2},{2,3},{3,4},{1,4}}
     U=random(ZZ^4,ZZ^4)
     solverMLE(G,U,SampleData=>true)
     V=sampleCovarianceMatrix(U)
     solverMLE(G,V,SampleData=>false)
  SeeAlso
     scoreEquations
     solverMLE
///

doc ///
  Key
    ConcentrationMatrix
  Headline
    optional input to output MLE for concentration matrix instead of MLE for covariance matrix
  SeeAlso
    solverMLE
   ///
doc ///
  Key
    [solverMLE, ConcentrationMatrix]
  Headline
    output MLE for concentration matrix instead of MLE for covariance matrix
  Usage
    solverMLE(G,U,ConcentrationMatrix=>b)  
  Inputs 
    b: Boolean
        false by default
  Description
    Text
     By default, @TO solverMLE@ outputs the MLE for the covariance matrix. By providing true as the value for the option
     @TO [solverMLE, ConcentrationMatrix]@, @TO solverMLE@ provides the MLE for
     the concentration matrix.

     Note that both the maximum value attained in the log-likelihood function and
     the ML-degree remain the same.
    Example
      G= mixedGraph(graph{{a,b},{b,c}},digraph {{a,d},{c,e},{f,g}},bigraph {{d,e}})
      solverMLE (G, random(QQ^7,QQ^7))
      solverMLE (G, random(QQ^7,QQ^7), ConcentrationMatrix=>true)
  SeeAlso
     solverMLE
///


doc   ///
    Key
    	checkPD
	(checkPD,List)
	(checkPD,Matrix)
    Headline
    	returns positive definite matrices from a list of symmetric  matrices
    Usage
    	checkPD(L)
    Inputs
    	L: List  
	   of symmetric matrices, or a single symmetric @TO Matrix@
    Outputs
    	 : List
	   of positive definite matrices
    Description
    	Text
	   This function takes a list of symmetric matrices (or a single symmetric matrix) and returns the sublist of its positive definite matrices. 
	   If there are no positive definite matrices in the list, it returns an empty list.
	   
	   Note that the function does not check whether the matrices in the original list are symmetric.
	   
	   If a matrix contains an imaginary part below the tolerance level, then only
	   the real part is reported in the output. (See  @TO [checkPD, ZeroTolerance]@)

      	Example
	    L={matrix{{1,0},{0,1}},matrix{{-2,0},{0,1}},matrix{{sqrt(-1),0},{0,sqrt (-1)}}}
    	    checkPD(L)
     	 ///

doc   ///
    Key
    	checkPSD
	(checkPSD,List)
	(checkPSD,Matrix)
    Headline
    	returns positive semidefinite matrices from a list of symmetric matrices
    Usage
    	checkPSD(L)
    Inputs
    	L: List  
	    of symmetric matrices, or a single symmetric @TO Matrix@
    Outputs
    	 : List
	   of positive semidefinite matrices
    Description
    	Text
	   This function takes a list of symmetric matrices (or a single symmetric matrix) and returns the sublist of its
	   positive semidefinite matrices.
          
	   Note that the function does not check whether the matrices in the original list are symmetric.
	  
	   If there are no positive semidefinite matrices in the list, it returns an empty list.
      	Example
	    L={matrix{{1,0},{0,1}},matrix{{-2,0},{0,1}},matrix{{sqrt(-1),0},{0,sqrt (-1)}},matrix{{0,0},{0,0}}}
    	    checkPSD(L)
     	 ///

doc ///
  Key
    ZeroTolerance
  Headline
    optional input to set the largest absolute value that should be treated as zero
  SeeAlso
     checkPD
     checkPSD
     solverMLE
   ///

doc ///
  Key
    [solverMLE, ZeroTolerance]
  Headline
    optional input to set the largest absolute value that should be treated as zero
  Usage
    solverMLE(G,U,ZeroTolerance=>n)  
  Inputs 
    n: RR
        default is 1e-10
  Description
    Text
     After computing the variety of the zero-dimensional ideal of the score equations,
     {\tt solverMLE} needs to determine which points lie in the PD cone. A matrix is
     assumed to be positive definite if for all eigenvalues e:

     -  @TO realPart@ e > ZeroTolerance

     -  @TO abs @ @TO imaginaryPart@ e <= ZeroTolerance

     If an MLE matrix contains an imaginary part below the tolerance level, then only
     the real part is reported in the output. (See  @TO [checkPD, ZeroTolerance]@)


  SeeAlso
     checkPD
///


doc ///
  Key
    [checkPD, ZeroTolerance]
  Headline
    optional input to set the largest absolute value that should be treated as zero
  Usage
    checkPD(L,ZeroTolerance=>n)  
  Inputs 
    n: RR
        default is 1e-10
  Description
    Text
     A matrix is assumed to be positive definite if for all eigenvalues e:

     -  @TO realPart@ e > ZeroTolerance

     -  @TO abs @ @TO imaginaryPart@ e <= ZeroTolerance

     If a matrix contains an imaginary part below the tolerance level, then only
     the real part is reported in the output.

    Example
     L={matrix{{10^(-9)+10^(-10)*sqrt(-1),0},{0,10^(-9)+10^(-10)*sqrt (-1)}},
	 matrix{{10^(-10)+10^(-10)*sqrt(-1),0},{0,10^(-10)+10^(-10)*sqrt (-1)}},
	 matrix{{1+10^(-10)*sqrt(-1),0},{0,1+10^(-10)*sqrt (-1)}},
	 matrix{{1-10^(-9)*sqrt(-1),0},{0,1+10^(-9)*sqrt (-1)}}
	 }
     checkPD L

  SeeAlso
     [checkPSD, ZeroTolerance]
     solverMLE
///

doc ///
  Key
    [checkPSD, ZeroTolerance]
  Headline
    optional input to set the largest absolute value that should be treated as zero
  Usage
    checkPD(L,ZeroTolerance=>n)  
  Inputs 
    n: RR
        default is 1e-10
  Description
    Text
     A matrix is assumed to be positive semidefinite if for all eigenvalues e:

     -  @TO realPart@ e >= -ZeroTolerance

     -  @TO abs @ @TO imaginaryPart@ e <= ZeroTolerance

     If a matrix contains an imaginary part below the tolerance level, then only
     the real part is reported in the output.

    Example
     L={matrix{{10^(-9)+10^(-10)*sqrt(-1),0},{0,10^(-9)+10^(-10)*sqrt (-1)}},
	 matrix{{10^(-10)+10^(-10)*sqrt(-1),0},{0,10^(-10)+10^(-10)*sqrt (-1)}},
	 matrix{{1+10^(-10)*sqrt(-1),0},{0,1+10^(-10)*sqrt (-1)}},
	 matrix{{1-10^(-9)*sqrt(-1),0},{0,1+10^(-9)*sqrt (-1)}}
	 }
     checkPD L

  SeeAlso
     [checkPD, ZeroTolerance]
///

doc   ///
    Key
    	MLdegree
	(MLdegree,Ring)
    Headline
    	ML-degree of a graphical model
    Usage
    	MLdegree(R)
    Inputs
    	R: Ring
	    defined as a @TO gaussianRing@ of @TO Graph@, or @TO Digraph@, or @TO Bigraph@, or @TO MixedGraph@
    Outputs
    	 : ZZ
	   the ML-degree of the model
    Description
    	Text
	   This function computes the ML-degree of a graphical model. It takes as input
	   a @TO gaussianRing@ of a @TO Graph@, or a @TO Digraph@, or a @TO Bigraph@, or a  @TO MixedGraph@.
	   It computes the degree of the score equation ideal given by @TO scoreEquations@
	   with a random sample data matrix.

	   We compute the ML-degree of the 4-cycle:
	Example
	    G=graph{{1,2},{2,3},{3,4},{4,1}}
            MLdegree(gaussianRing G)

     	 ///


doc ///
    Key
        solverMLE
	(solverMLE, Ring, List)
	(solverMLE, Ring, Matrix)
	(solverMLE, List, Ring)
	(solverMLE, Matrix, Ring)	
	(solverMLE,MixedGraph,List)
	(solverMLE, MixedGraph, Matrix)
	(solverMLE, Graph, List)
	(solverMLE, Digraph, List)
	(solverMLE, Bigraph, List)
	(solverMLE, Graph, Digraph,List)
	(solverMLE, Digraph, Graph,List)
	(solverMLE, Digraph, Bigraph, List)
	(solverMLE, Bigraph, Digraph, List)
	(solverMLE, Graph, Bigraph, List)
	(solverMLE, Bigraph, Graph, List)
	(solverMLE, Graph, Digraph, Bigraph, List)
	(solverMLE, Digraph, Bigraph, Graph, List)
	(solverMLE, Bigraph, Graph, Digraph, List)
	(solverMLE, Graph, Bigraph, Digraph, List)
	(solverMLE, Bigraph, Digraph, Graph, List)
	(solverMLE, Digraph, Graph, Bigraph, List)
	(solverMLE, List, MixedGraph)
	(solverMLE, List, Graph)
	(solverMLE, List, Digraph)
	(solverMLE, List, Bigraph)
	(solverMLE, List, Graph, Digraph)
	(solverMLE, List, Digraph,Graph)
	(solverMLE, List, Digraph,Bigraph)
	(solverMLE, List, Bigraph,Digraph)
	(solverMLE, List, Graph, Bigraph)
	(solverMLE, List, Bigraph, Graph)
	(solverMLE, List, Graph, Digraph, Bigraph)
	(solverMLE, List, Digraph, Bigraph, Graph)
	(solverMLE, List, Bigraph, Graph, Digraph)
	(solverMLE, List, Graph, Bigraph, Digraph)
	(solverMLE, List, Bigraph, Digraph, Graph)
	(solverMLE, List, Digraph, Graph, Bigraph)
	(solverMLE, Graph, Matrix)
	(solverMLE, Digraph, Matrix)
	(solverMLE, Bigraph, Matrix)
	(solverMLE, Graph, Digraph,Matrix)
	(solverMLE, Digraph, Graph,Matrix)
	(solverMLE, Digraph, Bigraph, Matrix)
	(solverMLE, Bigraph, Digraph, Matrix)
	(solverMLE, Graph, Bigraph, Matrix)
	(solverMLE, Bigraph, Graph, Matrix)
	(solverMLE, Graph, Digraph, Bigraph, Matrix)
	(solverMLE, Digraph, Bigraph, Graph, Matrix)
	(solverMLE, Bigraph, Graph, Digraph, Matrix)
	(solverMLE, Graph, Bigraph, Digraph, Matrix)
	(solverMLE, Bigraph, Digraph, Graph, Matrix)
	(solverMLE, Digraph, Graph, Bigraph, Matrix)
	(solverMLE, Matrix, MixedGraph)
	(solverMLE, Matrix, Graph)
	(solverMLE, Matrix, Digraph)
	(solverMLE, Matrix, Bigraph)
	(solverMLE, Matrix, Graph, Digraph)
	(solverMLE, Matrix, Digraph,Graph)
	(solverMLE, Matrix, Digraph,Bigraph)
	(solverMLE, Matrix, Bigraph,Digraph)
	(solverMLE, Matrix, Graph, Bigraph)
	(solverMLE, Matrix, Bigraph, Graph)
	(solverMLE, Matrix, Graph, Digraph, Bigraph)
	(solverMLE, Matrix, Digraph, Bigraph, Graph)
	(solverMLE, Matrix, Bigraph, Graph, Digraph)
	(solverMLE, Matrix, Graph, Bigraph, Digraph)
	(solverMLE, Matrix, Bigraph, Digraph, Graph)
	(solverMLE, Matrix, Digraph, Graph, Bigraph)
    Headline
    	Maximum likelihood estimate of a loopless mixed graph
    Usage
    	solverMLE(G,U)
    Inputs
    	G:Graph
 	  or @ofClass Digraph@, or @ofClass Bigraph@, or @ofClass MixedGraph@
	U:Matrix
	  or @ofClass List@ of sample data.
	  Alternatively, the sample covariance matrix can be given as input
	  by setting @TO SampleData@ to false

    Outputs
        : Sequence
	   of length 3, whose first element is the maximum value attained in the log-likelihood function (of type @TO RR@), 
	   its second element is the MLE (or MLEs) of the covariance matrix (of types @TO Matrix@ or @TO List@), 
	   and its third element is the ML-degree of the model (of type @TO ZZ@).
	   
	   By providing true as the value for the option @TO ConcentrationMatrix@, the MLE for the concentration matrix can be given as output.

    Description

	Text
	    This function takes as input a @TO2{Graph,"graph"}@, or a @TO2{Digraph,"digraph"}@, or a @TO2{Bigraph,"bigraph"}@ or a @TO2{MixedGraph,"mixed graph"}@ and a list or matrix that encodes, by default, the sample data.
	    It computes the critical points of the score equations and 
	    selects the maximum value achieved among those that lie in the cone of positive-definite matrices.
	    The default output is the maximum value in the log-likelihood function, maximum likelihood estimators (MLE) for the covariance matrix
	    and the ML-degree of the model.
	    MLE for the concentration matrix can be obtained by setting the optional input @TO ConcentrationMatrix@ to false.

	    The same optional inputs as in @TO scoreEquations@ can be used, plus extra optional inputs related to
	    the numerical solver (EigenSolver by default, NAG4M2 alternatively) and its functionalities.

	    Below we reproduce Example 2.1.13 for the 4-cycle in the book: Mathias Drton, Bernd Sturmfels and Seth Sullivant,
	    {\em Lectures on Algebraic Statistics}, Oberwolfach Seminars, Vol 40, Birkhauser, Basel, 2009.

	Example
	    G=graph{{1,2},{2,3},{3,4},{1,4}};
	    U =matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
	    solverMLE(G,U)

	Text
	    The data sample can also be given as a list:

	Example
            G=graph{{1,2},{2,3},{3,4},{1,4}};
	    U = {{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
	    solverMLE(G,U)

	Text
	    In the following example we compute the MLE for the covariance matrix of
	    the graphical model associated to the graph $1\rightarrow 2,1\rightarrow 3,2\rightarrow 3,3\rightarrow 4,3<-> 4$
	    In this case we give as input the sample covariance matrix:

	Example
	    G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}});
            S =  matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}};
            solverMLE(G,S,SampleData=>false)

	Text
	   Next we provide the MLE for the concentration matrix of the graphical model
	   associated to the graph $1\rightarrow 3,2\rightarrow 4,3<->4,1 - 2$.
           Again the sample covariance matrix is given as input.

	Example
	   G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}},graph {{1,2}});
           S =  matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}};
           solverMLE(G,S,SampleData=>false,ConcentrationMatrix=>true)

	Text
          {\bf Application to positive definite matrix completion problems}
        Text
          Consider the following symmetric matrix with some unknown entries:
        Example
          R=QQ[x,y];
          M=matrix{{115,-13,x,47},{-13,5,7,y},{x,7,27,-21},{47,y,-21,29}}
        Text
          Unknown entries correspond to non-edges of the 4-cycle. A positive definite completion of this matrix
          is obtained by giving values to x and y and computing the MLE for the covariance matrix in the Gaussian graphical model
          given by the 4-cycle. To understand which values of x and y will result in a maximum likelihood estimate,
	  see Example 12.16 in the book: Mateusz Michalek and Bernd Sturmfels,
	  {\em Invitation to Nonlinear Algebra}, Graduate Studies in Mathematics, Vol ???, American Mathematical Society, 2021.

        Example
          G=graph{{1,2},{2,3},{3,4},{1,4}};
          V=matrix{{115,-13,-29,47},{-13,5,7,-11},{-29,7,27,-21},{47,-11,-21,29}}
          (mx,MLE,ML)=solverMLE(G,V,SampleData=>false)
        Text
	  The MLE of the covariance matrix is the unique positive definite completion of the matrix M such that its inverse, namely the concentration matrix,
	  has zero's in the entries corresponding to non-edges of the graph.
	  Observe that all entries of V remain the same in the MLE except for those that correspond to non-edges of the graph.

    SeeAlso
       checkPD
       checkPSD
       scoreEquations
       jacobianMatrixOfRationalFunction
       sampleCovarianceMatrix
///

--******************************************--
-- TESTS     	       	    	      	    --
--******************************************--

TEST /// --test jacobianMatrixOfRationalFunction
R=QQ[x,y];
FR=frac R;
F=1/(x^2+y^2);
M=entries jacobianMatrixOfRationalFunction(F);
N=transpose {{-2*x/(x^2 + y^2)^2,-2*y/(x^2 + y^2)^2 }};
assert(M === N)
///

TEST /// --test jacobianMatrixOfRationalFunction
R=QQ[x_1,x_2,x_3];
FR=frac R;
M=entries jacobianMatrixOfRationalFunction( (x_1^2*x_2)/(x_1+x_2^2+x_3^3) );
N=transpose {{2*x_1*x_2/(x_2^2 + x_3^3 + x_1) - x_1^2*x_2/(x_2^2 + x_3^3 + x_1)^2, -2*x_1^2*x_2^2/(x_2^2 + x_3^3 + x_1)^2 + x_1^2/(x_2^2 + x_3^3 + x_1) , -3*x_1^2*x_2*x_3^2/(x_2^2 + x_3^3 + x_1)^2 }};
assert(M === N)
///

TEST /// --test sampleCovarianceMatrix
M = matrix{{1, 2, 0},{-1, 0, 5/1},{3, 5, 2/1},{-1, -4, 1/1}};
N = sampleCovarianceMatrix(M);
A = matrix {{11/4, 39/8, -1}, {39/8, 171/16, 0}, {-1, 0, 7/2}};
assert(N===A)
///

TEST /// --test sampleCovarianceMatrix with sample of bigger size than the covariance matrix
X = matrix{{36, -3, -25, -36},{-10, 11, -29, -20},{-18, 33, -15, -11},{-42, 0, 20, 43}, {-30, -26, 32, 2},{2, -38, -24, -43}};
Y = sampleCovarianceMatrix(X);
B = matrix {{5621/9, -1037/18, -7835/18, -10565/18}, {-1037/18, 19505/36, -4897/36, 5147/36}, {-7835/18, -4897/36, 20465/36, 18941/36}, {-10565/18, 5147/36, 18941/36, 28889/36}};
assert(Y===B)
///

TEST /// --test sampleCovarianceMatrix with list input
X = {{48,89,27,28},{23,19,29,94},{135,23,44,71},{91,75,24,98}};
Y = sampleCovarianceMatrix(X);
B = matrix {{29147/16, -1313/8, 220, 1609/16}, {-1313/8, 3827/4, -155, -3451/8}, {220, -155, 119/2, -63/4}, {1609/16, -3451/8, -63/4, 12379/16}};
assert(Y===B)
///

TEST /// --test scoreEquations for a mixedGraph with directed and bidirected edges
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
R=gaussianRing(G);
U = matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
J=scoreEquations(R,U);
I=ideal(20*p_(3,4)+39,50*p_(4,4)-271,440104*p_(3,3)-742363,230*p_(2,2)-203,16*p_(1,1)-115,5*l_(3,4)+2,110026*l_(2,3)-2575,55013*l_(1,3)-600,115*l_(1,2)+26);
assert(J===I)
///

TEST ///  --test scoreEquations for a graph (and random input data)
G=graph{{1,2},{2,3},{3,4},{1,4}};
R=gaussianRing(G);
U=random(ZZ^4,ZZ^4);
J=scoreEquations(R,U);
assert(dim J===0)
assert(degree J===5)
///

TEST /// --score equations of sample data equals score equations of its sample covariance data with SampleData=>false
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}})
R = gaussianRing(G)
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}}
JU=scoreEquations(R,U)
RU=ring(JU)
V = sampleCovarianceMatrix U
JV=scoreEquations(R,V,SampleData=>false)
JV=sub(JV,RU)
assert(JU==JV)
///

TEST /// --score equations with elimination strategy equals default saturation strategy
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}})
R = gaussianRing(G)
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}}
JU=scoreEquations(R,U)
RU=ring(JU)
J=scoreEquations(R,U,SaturateOptions => {Strategy => Eliminate})
J=sub(J,RU)
assert(J==JU)
///

TEST /// --test checkPD
L={matrix{{1,0},{0,1}},matrix{{-2,0},{0,1}},matrix{{sqrt(-1),0},{0,sqrt (-1)}},matrix{{0.0001*sqrt(-1),0},{0,0.0000001*sqrt (-1)}}};
Y = checkPD(L);
B = {matrix{{1, 0}, {0, 1}}};
assert(Y===B)
///

TEST /// --test ZeroTolerance checkPD
L={matrix{{10^(-9)+10^(-10)*sqrt(-1),0},{0,10^(-9)+10^(-10)*sqrt (-1)}}, matrix{{10^(-10)+10^(-10)*sqrt(-1),0},{0,10^(-10)+10^(-10)*sqrt (-1)}},matrix{{1+10^(-10)*sqrt(-1),0},{0,1+10^(-10)*sqrt (-1)}},matrix{{1-10^(-9)*sqrt(-1),0},{0,1+10^(-9)*sqrt (-1)}}}
assert(checkPD L ==={matrix {{1e-9, 0}, {0, 1e-9}}, sub(matrix {{1, 0}, {0, 1}},RR)})
///

TEST /// --test ZeroTolerance checkPSD
L={matrix{{10^(-9)+10^(-10)*sqrt(-1),0},{0,10^(-9)+10^(-10)*sqrt (-1)}}, matrix{{-10^(-10)+10^(-10)*sqrt(-1),0},{0,-10^(-10)+10^(-10)*sqrt (-1)}},matrix{{1+10^(-10)*sqrt(-1),0},{0,1+10^(-10)*sqrt (-1)}},matrix{{1-10^(-9)*sqrt(-1),0},{0,1+10^(-9)*sqrt (-1)}}}
assert(checkPSD L ==={matrix {{1e-9, 0}, {0, 1e-9}},matrix {{-1e-10, 0}, {0, -1e-10}}, sub(matrix {{1, 0}, {0, 1}},RR)})

///
TEST /// --test checkPSD
L={matrix{{1,0},{0,1}},matrix{{-2,0},{0,1}},matrix{{sqrt(-1),0},{0,sqrt (-1)}},matrix{{0.0001*sqrt(-1),0},{0,0.0000001*sqrt (-1)}},matrix{{0,0},{0,0}}};
Y = checkPSD(L);
B = {matrix{{1, 0}, {0, 1}},matrix{{0,0},{0,0}}};
assert(Y===B)
///

TEST /// --test MLdegree
G=graph{{1,2},{2,3},{3,4},{4,1}}
assert(MLdegree(gaussianRing G)==5)
///

TEST /// --test solverMLE graph
G=graph{{1,2},{2,3},{3,4},{1,4}}
U =matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}}
(mx,MLE,ML)=solverMLE(G,U)
assert(round(4,mx)==-6.2615)
assert(ML==5)
///

TEST /// --test solverMLE for mixedGraph with directed and bidirected edges
G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}})
S =  matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}}
(mx,MLE,ML)=solverMLE(G,S,SampleData=>false)
assert(ML==5)
assert(round(5, mx)==-1.14351)
///

TEST /// --test solverMLE for mixedGraph with all kind of edges and concentration matrix
G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}},graph {{1,2}})
S =  matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}}
(mx,MLE,ML)=solverMLE(G,S,SampleData=>false,ConcentrationMatrix=>true)
assert(ML==5)
assert(round(4,mx)==-.8362)
///

TEST /// --test solverMLE graph with complete test
G=graph{{1,2},{2,3},{3,4},{1,4}}
V =matrix{{5,1,3,2},{1,5,1,6},{3,1,5,1},{2,6,1,5}}
(mx,MLE,ML)=solverMLE(G,V,SampleData=>false,ConcentrationMatrix=>true)
assert(round(4,mx)==-10.1467)
assert(ML==5)
assert(round(6,MLE_(0,1))== -.038900)
assert(round(6,MLE_(3,1))== 0)
///


TEST /// -- test solverMLE with NAG4M2
G=mixedGraph(graph{{a,b}},digraph{{a,d}},bigraph{{c,d}})
U=matrix{{1, 2, 5, 1}, {5, 3, 2, 1}, {4, 3, 5, 10}, {2, 5,1, 3}}
(mx,MLE,ML)= solverMLE (G,U,Solver=>"NAG4M2")
assert(round(5,mx)==-8.4691)
assert(ML==1)
assert(round(6,MLE_(2,0))==0)
assert(round(6,MLE_(1,1))== 1.1875)
assert(round(6,MLE_(3,2))==   3.264929)
///
--------------------------------------
--------------------------------------
end--
--------------------------------------
--------------------------------------