(* ::Package:: *)

(* Mathematica Package *)

BeginPackage["injectivity`"]
Unprotect @@ Names["injectivity`*"];
ClearAll @@ Names["injectivity`*"];
  
determinantCriterion::usage=" determinantCriterion[A,V,S] applies the determinant criterion to the generalized polynomial map represented
matrices A and V, with respect to S. Returns True if the family is S-injective, False otherwise"

primeMatrices::usage="prime[A,V] calculates the matrices A' and V' corresponding to the matrices A and V."

steps123::usage=" steps123[A,V] performs steps 1,2 and 3 of the method. It calculates matrices A',V' and uses 
the determinant criterion to check if the new family is S-injective. Returns True if the family
is S-injective, False otherwise"

proposition21::usage=" proposition21[A] checks if the assumptions of Proposition 2.1 are fulfilled. Returns True if the 
assumptions hold, False otherwise"

step5::usage=" step5[A,V] checks step 5 of the method. If it holds returns True (conclusion: family has multiple
S-zeros), othrwise returns False"

method::usage=" method[A,V] combines the functions above, and returns True if the family is S-injective, False if it
has multiple S-zeros, Inconlusive if no conclusion was drawn"

jacobianInjectivity::usage="jacobianInjectivity[A,V] applies the jacobian injectivity method to the generalized polynomial 
represented by matrices A and V. It returns Inconclusive if no conclusion was drawn, Injective in the situation in which it can be asserted,
returns the conditions that define a region for the parameters where we can conclude injectivity."

jacobianFunction::usage="jacobianFunction[f,x,kappa] gets the polynomial f, the variables x and the rate constants kappa, and applies the
Jacobian Injectivity method. It returns Injective, Inconclusive or the conditions that define a region for the parameters where injectivity
occurs."

(* private context *)
Begin["`Private`"] 

(* function definitions *)
determinantCriterion[A_?MatrixQ,V_?MatrixQ,S_?MatrixQ]/;Equal[Length[Transpose[A]],Length[Transpose[V]]]:=
	Module[{Z,Atilde,kappa,lambda,n=Length[V],r=Length[Transpose[A]],s=MatrixRank[A], M, p,k,l},
	
	(* if the rank of A is less than the rank of S, then we conclude that the family is not S-injective *)
	If[MatrixRank[A]<MatrixRank[S], Return[False]];
	
	
	(* Z is a matrix whose rows are a basis for the orthogonal space of S *)
	Z = NullSpace[Transpose[S]];
	
	(* Atilde has s linearly independent rows of A *)
	Atilde = Take[RowReduce[A],s];
	
	(* kappa and lambda are parameters *)
	kappa = Table[k[i],{i,r}];
	lambda = Table[l[i],{i,n}];
	
	M = Join[Z, Atilde.DiagonalMatrix[kappa].Transpose[V].DiagonalMatrix[lambda]];
	
	(* if the determinant is identically zero, return False.
	   if the determinant is not identically zero, return True,
	   the family is S-injective *)
	p = Det[M];
	
	Return[!PossibleZeroQ[p] && Length[DeleteDuplicates[Sign[CoefficientRules[p,Join[kappa,lambda]][[All,2]]]]]==1];
	];

primeMatrices[A_?MatrixQ,V_?MatrixQ]/;Equal[Length[Transpose[A]],Length[Transpose[V]]]:=
	Module[{A1,n=Length[V],r=Length[Transpose[A]],s=MatrixRank[A],
	i,j,g,x,VPrime,APrime,monomials,variables,aux,rPrime, degree, Atilde,auxA,auxV},
	ClearAll[A1,i,j,g,x,VPrime,APrime,monomials,variables,aux,rPrime, degree, Atilde,auxA,auxV];
	auxA=A;
	auxV=V;
	(* A1 is the minor of the Gauss reduction of A *)
	Atilde = Take[RowReduce[auxA],{1,s},{1,r}];
	
	(* Swap the columns of Atilde in order to get the identity matrix in the first sxs minor, and the corresponding ones in A and V*)
	For[i=2,i<=s,i++,
		j=i;
		(* increment j until a nonzero entry is reached *)
		While[Atilde[[i,j]]==0,j++];
		
		(* Swap columns i and j of Atilde, A and V *)
		aux = Atilde[[All,i]];
		Atilde[[All,i]] = Atilde[[All,j]];
		Atilde[[All,j]] = aux;
		
		aux = auxA[[All,i]];
		auxA[[All,i]] = auxA[[All,j]];
		auxA[[All,j]] = aux;
		
		aux = auxV[[All,i]];
		auxV[[All,i]] = auxV[[All,j]];
		auxV[[All,j]] = aux;		
	];	
	
	(* A1 is the minor of a matrix similar to A *)
	A1 = Take[Atilde,{1,s},{s+1,r}];
	(* assume k 1 *)
	aux= A1.Table[Product[x[i]^{auxV[[i,s+j]]}, {i,n}], {j,r-s}];
	
	(* get smaller element of V, in order to avoid non-classic polynomials *)
	degree = Abs[Min[auxV,-auxV]];
	
	(* calculate polynomial map g kappa *)
	(* multiply g by every x, so that it can be manipulated as a polynomial *)
	g= Table[-Product[x[i]^degree,{i,n}]*Product[x[i]^{-auxV[[i,j]]}, {i,n}]*aux[[j]], {j,s}]; 
	
	(* get matrices A' and V' from the expression of g. *)
	(* used Abs to avoid saving z and -z, for example. *)
	variables=Table[x[i],{i,n}];
	monomials = DeleteDuplicates[Flatten[CoefficientRules[g,variables]][[All,1]]];
	
	(* get exponents from different monomials, in order to define V' *)
	VPrime= Transpose[monomials]-degree
	
	(* number of distinct monomials *)
	rPrime = Length[monomials];
	
	(* for all g_i, check which monomials includes, and get its coefficient. *)
	APrime = 
	Table[
		If[ MemberQ[CoefficientRules[g[[i,1]],variables][[All,1]], monomials[[j]]], 
		Extract[CoefficientRules[g[[i,1]],variables],
		FirstPosition[CoefficientRules[g[[i,1]],variables][[All,1]],monomials[[j]],0]] [[2]],
		0],
		{i,s},{j,rPrime}
	]
	
	Return[Join[APrime,VPrime]];
	];
	
steps123[A_?MatrixQ,V_?MatrixQ]/;Equal[Length[Transpose[A]],Length[Transpose[V]]]:=
	Module[{APrime,VPrime,aux},
	ClearAll[APrime,VPrime,aux];
	
	(* get A' and V' *)
	aux=primeMatrices[A,V];
    APrime = aux[[1;;MatrixRank[A],All]];
    VPrime = aux[[MatrixRank[A]+1;;,All]];

	(* apply the Determinant Criterion to A',V', with respect to Im(A) *)
	Return[determinantCriterion[APrime,VPrime,A]];
	];

proposition21[A_?MatrixQ]:=
	Module[{A1,r=Length[Transpose[A]],s=MatrixRank[A],
	i,j,x,nonZero,Atilde,aux},
	ClearAll[A1,i,j,x,nonZero,Atilde,aux];

	(* A1 is the minor of the Gauss reduction of A. We are interested in its columns, 
	so save the transpose of A1 as A1. Also consider sign and absolute value. *)
	(* A1 is the minor of the Gauss reduction of A *)
	Atilde = Take[RowReduce[A],{1,s},{1,r}];
	
	(* Swap the columns of Atilde in order to get the identity matrix in the first sxs minor *)
	For[i=2,i<=s,i++,
		j=i;
		(* increment j until a nonzero entry is reached *)
		While[Atilde[[i,j]]==0,j++];
		
		(* Swap columns i and j of Atilde, A and V *)
		aux = Atilde[[All,i]];
		Atilde[[All,i]] = Atilde[[All,j]];
		Atilde[[All,j]] = aux;
	];	
	
	(* A1 is the minor of a matrix similar to A *)
	A1 = Transpose[Take[Atilde,{1,s},{s+1,r}]];
	
	(* check the assumptions. *)
	(* table whose elements give the number of nonzero entries of each column *)
	nonZero=Table[Count[A1[[i]], x_/;x!=0], {i,r-s}];
	
	(* loop to get the number of common elements in columns of A1 *)
	Return[Catch[
		For[i=1,i<=r-s,i++,
			If[nonZero[[i]]>1,
				For[j=i+1,j<=r-s,j++,
				If[And[nonZero[[j]]>1,Count[A1[[i]]+A1[[j]], x_/;x!=0]!=nonZero[[i]]+nonZero[[j]]-1], Throw[False]]
				]
			]
		]
		Throw[True];
	]];
	];

step5[A_?MatrixQ,V_?MatrixQ]/;Equal[Length[Transpose[A]],Length[Transpose[V]]]:=
	Module[{A1,n=Length[V],r=Length[Transpose[A]],s=MatrixRank[A],
	i,j,g,x,aux,Atilde},
	ClearAll[A1,i,j,g,x,aux,Atilde];
	
	(* A1 is the minor of the Gauss reduction of A *)
	Atilde = Take[RowReduce[A],{1,s},{1,r}];
	
	(* Swap the columns of Atilde in order to get the identity matrix in the first sxs minor *)
	For[i=2,i<=s,i++,
		j=i;
		(* increment j until a nonzero entry is reached *)
		While[Atilde[[i,j]]==0,j++];
		
		(* Swap columns i and j of Atilde, A and V *)
		aux = Atilde[[All,i]];
		Atilde[[All,i]] = Atilde[[All,j]];
		Atilde[[All,j]] = aux;
	];	
	
	A1 = Take[Atilde,{1,s},{s+1,r}];
	
	(* Check if all entries are non-positive and if every row has at least one nonzero entry *)
	If[And[Count[A1,x_/;x>0,2]==0, Count[Table[Count[A1[[i]],x_/;x<0],{i,s}],y_/;y>=1]==s], Return[True],Return[False]];
	
	(* If the above condition isn't met, it is necessary to check the assumptions of Prop. 2.2 *)
	
	];
	

method[A_?MatrixQ,V_?MatrixQ]/;Equal[Length[Transpose[A]],Length[Transpose[V]]]:=
	Module[{A1,n=Length[V],r=Length[Transpose[V]],s=MatrixRank[A],
	m,it,k,auxA,auxV,indices},
	
	ClearAll[A1,m,it,k,auxA,auxV,indices]; 
	
	(* loop to account for permutations of columns of A and V. *)
	(* all sets of s indices out of r *)
	indices=Subsets[Table[m,{m,r}],{s}];
	For[m=1,m<Length[indices],m++,
		ClearAll[auxA,auxV];
		(* check if the first s columns of A are linearly independent *)
		If[MatrixRank[Transpose[Transpose[A][[ indices[[m]],All ]] ]]==s,
			auxA= Transpose[Join[Transpose[A][[indices[[m]],All]], 
				                                                  Transpose[A][[Complement[Table[it,{it,r}],indices[[m]] ],All ]]]];
			auxV= Transpose[Join[Transpose[V][[indices[[m]],All]], 
				                                                  Transpose[V][[Complement[Table[it,{it,r}],indices[[m]] ],All ]]]];
		
			(* if the first sxs minor of A is nonzero and if the first s columns are linearly independent *)
			If[Total[Total[ Abs[auxA[[1;;s,1;;s]] ]]]!=0,
				(* apply the determinant criterion with A and V, with respect to Im(A) *)
				Which[determinantCriterion[auxA,auxV,A], Return[StringForm["True 0:``",m]],
	
			      (* determinant criterion failed for A,V, couldnt prove the injectivity. Then, apply steps 1,2,3 *)
	               steps123[auxA,auxV], Return[StringForm["True 123:``",m]],
	
			      (* determinant criterion failed for A', V'. Check assumptions for Prop.2.1 *)
			      And[proposition21[auxA], step5[auxA,auxV]], Return[StringForm["False 5:``",m]] 

				]		
				(* no conclusion was reached. It is necessary to permute the columns of A and V and restart. *)
	
			];
		];
	];
	Return[StringForm["Inconclusive:``",m]];
	];

jacobianInjectivity[A_?MatrixQ,V_?MatrixQ]/;Equal[Length[Transpose[A]],Length[Transpose[V]]]:=
	Module[{n=Length[V],r=Length[Transpose[A]],s=MatrixRank[A],i,x,f,y, jacobian,determinant,k,coeffDet,coefficients},

	(* check if the polynomial map has degree less or equal to 2 *)
	For[i=1,i<=r,i++,
	If[Total[Transpose[V][[i]]]>2, Return["Inconclusive"]]];

	(* calculate polynomial map from matrices *)
	f= A.Table[k[j]*Product[x[i]^{V[[i,j]]}, {i,n}], {j,r}];

	(* calculate jacobian *)
	jacobian = Flatten[Grad[f,Table[x[i],{i,n}]],1];

	(* calculate the jacobian determinant *)
	determinant = Det[jacobian];

	(* check if determinant is zero *)
	If[PossibleZeroQ[determinant], Return["Inconclusive"]];

	(* get coefficients of the determinant *)
	coeffDet = CoefficientRules[determinant,Table[x[i],{i,n}]][[All,2]];

	(* if all the coefficients are positive, conclude injectivity *)
	coefficients=CoefficientRules[coeffDet,Table[k[i],{i,r}]];
	If[Count[Flatten[coefficients][[All,2]],y_/;y>0]==Length[Flatten[coefficients][[All,2]]], Return["Injective"]];

	(* if there are coefficients with negative summands, give the expression for the region where the map is injective *)
	Print["The following conditions are necessary for injectivity: "];
	For[i=1,i<=Length[coefficients],i++,
		If[Count[coefficients[[i]][[All,2]],y_/;y<0]>0, Print[Simplify[coeffDet[[i]]>0]]];
	];	

];

jacobianFunction[f_,var_,kappa_]:=
	Module[{x,y,i,j,jacobian,determinant,coeffDet,coefficients,k},

	(* check if the polynomial map has degree less or equal to 2 *)
	For[i=1,i<=Length[f],i++,
	If[Max[Table[Total[Flatten[CoefficientRules[f,var][[i,All,All,1]],1][[j]]],{j,2}]]>2,Return["Inconclusive"]]];

	(* calculate jacobian *)
	jacobian = Flatten[Grad[f,var],1];

	(* calculate the jacobian determinant *)
	determinant = Det[jacobian];

	(* check if determinant is zero *)
	If[PossibleZeroQ[determinant], Return["Inconclusive"]];

	(* get coefficients of the determinant *)
	coeffDet = CoefficientRules[determinant,var][[All,2]];

	(* if all the coefficients are positive, conclude injectivity *)
	coefficients=CoefficientRules[coeffDet,Table[k[i],{i,Length[kappa]}]];
	If[Count[Flatten[coefficients][[All,2]],y_/;y>0]==Length[Flatten[coefficients][[All,2]]], Return["Injective"]];

	(* if there are coefficients with negative summands, give the expression for the region where the map is injective *)
	Print["The following conditions are necessary for injectivity: "];
	For[i=1,i<=Length[coefficients],i++,
		If[Count[coefficients[[i]][[All,2]],y_/;y<0]>0, Print[Simplify[coeffDet[[i]]>0]]];
	];	

];	
	
(* end private context *)
End[]
Protect @@ Names["injectivity`*"];
EndPackage[]



