(* ::Package:: *)

(*===========================================================================================*========
HnfDecomposition package
Last Change: Friday February 8 th 2013
Mathematica package for decomposition of ideals of cancellative semigroups
Authors: J. I. Garc\[IAcute]a Garc\[IAcute]a (ignacio.garcia@uca.es) , A. Vigneron Tenorio (alberto.vigneron@uca.es)
=====================================================================================================*)

BeginPackage["HnfDecomposition`"];

Authors="J. I. Garc\[IAcute]a-Garc\[IAcute]a and A. Vigneron-Tenorio,\nDpto. Matem\[AAcute]ticas, Universidad de C\[AAcute]diz\n
ignacio.garcia@uca.es, alberto.vigneron@uca.es\n";
Commands="The package defines the following commands:\n
hnfDecomposition, exec4ti2, exec4ti2Blocks\n"


Print["\nHnfDecomposition package by\n\n"<>Authors<>"\n"
<>Commands<>
"Use the Information command to know how to use them\n"];

Needs["Combinatorica`"];

hnfDecomposition::usage="{permutations,blocks} = hnfDecomposition[M]\n
Computes the HNF-decomposition of a matrix \"M\".\n 
The returning values are the permutations of the columns of the matrix and the block matrices";

exec4ti2::usage="exec4ti2[comm,path,M]\n
Execute the 4ti2 command \"comm\" which can be equal to \"4ti2int64 groebner\", \"4ti2int32 groebner\", \"4ti2gmp groebner\" or \"zsolve -G\". The entries are the command \"comm\", the path \"path\" where 4ti2 will store the output matrices and the matrix \"M\".\n 
The returning value is the time of the execution and the returning code of the 4ti2 command.\n
The directory of the program 4ti2 must be in the PATH enviroment variable of your computer.
";

exec4ti2Blocks::usage="exec4ti2List[cmm,path,listOfMatrices]\n
Execute the 4ti2 command \"comm\" which can be equal to 
\"4ti2int64 groebner\", \"4ti2int32 groebner\", \"4ti2gmp groebner\" or \"zsolve -G\". The entries are the command \"comm\", the path \"path\" where 4ti2 will store the output matrices and the list of matrices listOfMatrices.\n 
The returning value is the time of the execution and the returning code of the 4ti2 command.\n
The execution is done in parallel and the time of execution depends mainly on the number of processor of the computer.\n
The directory of the program 4ti2 must be in the PATH enviroment variable of your computer.
";



Begin["Private`"];


matrizAFichero[matriz_,workingDir_,nameFile_]:=Module[{i=0,j=0,cad="",k,fich},
k=Dimensions[matriz];
cad=ToString[k[[1]]]<>" "<>ToString[k[[2]]]<>"\n";
For[i=1,i<=Length[matriz],i++,
For[j=1,j<=Length[matriz[[i]]],j++,cad=cad<>ToString[matriz[[i]][[j]]] <>" "];
cad=cad<>"\n";
];
fich=OpenWrite[workingDir<>"\\"<>nameFile<>".mat"];
WriteString[fich,cad];
Close[fich]
];

listaMAFicheros[lista_,wd_,prefijo_]:=Module[{l=0,i},
l=Length[lista];
For[i=1,i<=l,i++,
matrizAFichero[lista[[i]],wd,prefijo<>ToString[i]];
];
];

ejecuta4ti2[fti2_,p_,nFichero_]:=Run[fti2<>" "<>"\""<>p<>"\\"<>nFichero<>"\""];

exec4ti2[fti2_,path_,m_]:=Module[{},
matrizAFichero[m,path,"aa"];
AbsoluteTiming[ejecuta4ti2[fti2,path,"aa"]]
];

ejecuta4ti2Paralelo[fti2_,p_,prefijo_,listaMatrices_]:=AbsoluteTiming [ParallelMap[  ejecuta4ti2Gro[fti2<>" ",p,#]&,listaNombres[prefijo,Length[ listaMatrices ] ] ]  ];

listaNombres[prefijo_,longlista_]:=Table[prefijo<>ToString[i],{i,longlista}];

exec4ti2Blocks[fti2_,path_,lm_]:=Module[{longlista=Length[lm]},
listaMAFicheros[lm,path,"aaa"];
AbsoluteTiming [ParallelMap[  ejecuta4ti2[fti2,path,#]&,listaNombres["aaa",longlista] ]  ]
];


listarel2listarel[lr_,{}]:=lr;
listarel2listarel[lr_,{a_}]:=lr;
listarel2listarel[lr_,rel_]:=Module[{v,i,laux=lr},
If[lr[[ First[rel] ]]<First[rel],
v=lr[[ First[rel] ]],
v=First[rel];
];
For[i=1,i<=Length[rel],i++,laux=laux/.(laux[[ rel[[i]] ]])-> v];
laux
];
listarelaux[l_,{}]:=l;
listarelaux[l_,ll_]:=listarelaux[listarel2listarel[l,First[ll]],Rest[ll]];
listarel[m_,nc_]:=listarelaux[Table[i,{i,1,nc }],m];
agrupaciones[lista_]:=Module[{lon=Length[lista],i,laux={}},
For[i=1,i<=Length[lista],i++,
laux=Append[laux, Select[Table[i,{i,lon}],lista[[#]]==lista[[i]]&]];
];
laux=DeleteDuplicates[laux];
laux
];

divensubmatrices[A_,{}]:={};
divensubmatrices[A_,{{}}]:={};
divensubmatrices[A_,{pult_}]:={tomaCols[A,pult]};
divensubmatrices[A_,lp_]:=
If[First[lp]!={},
Join[{tomaCols[A,First[lp]]},
(divensubmatrices[A,Rest[lp]])],
divensubmatrices[A,Rest[lp]]
];




noeslistacero[{}]:=False;

noeslistacero[lista_]:=If[
First[lista]!=0,
True,
noeslistacero[Rest[lista]]
];

quitaceros[A_]:=Select[A,noeslistacero[#] & ];


tomaCols[A_,{}]:={};
tomaCols[A_,lista_]:=quitaceros[


\!\(\*SuperscriptBox[
RowBox[{"(", " ", 
RowBox[{"Part", "[", 
RowBox[{
SuperscriptBox["A", "\[Transpose]"], ",", "lista"}], "]"}], " ", ")"}], "\[Transpose]"]\)];


hnfDecomposition[m_]:=Module[{mi,ts,perm,maux,nbl,m1,t1,matriz2listapNN,tieneComunes,agrupaAux,agrupa,listapNN1,numcolsm,lrels,orden,numfilasm,
agr},
<<Combinatorica`;
$RecursionLimit=400;


listapNN1[l_]:=Select[Table[i,{i,Length[l]}],l[[#]]!=0 &];
matriz2listapNN[mm_]:=listapNN1/@ mm;

{mi,ts}=HermiteDecomposition[m];
maux=matriz2listapNN[ts];
numfilasm=Dimensions[m][[1]];
numcolsm=Dimensions[m][[2]];
lrels=listarel[maux,numcolsm];
agr=agrupaciones[lrels];
orden=Flatten[agr];
{mi,m,mi.m,
agr,
mi.m.



\!\(\*SuperscriptBox[
RowBox[{"(", 
RowBox[{"Permute", "[", 
RowBox[{
RowBox[{"IdentityMatrix", "[", "numcolsm", "]"}], ",", "orden"}], "]"}], ")"}], "\[Transpose]"]\),
divensubmatrices[mi.m,agr]};
{agr,divensubmatrices[mi.m,agr]}
];

End[];

EndPackage[];



