(*
   Eliminate.nb - M. Q. Rieck 

 Eliminate some angles and reduce system of inequalities 
 associated with the P3P triangle.
*)

(* 
   Displays the inequalities represented by the matrix. 
 Each row holds the coefficients of π, ∡A, ∡B, α, β, 
   γ, α', β', γ' for a strictly positive quantity. 
*)
showIneqs[mat_] := Module[ {i,j,len,abs,s}, 
len = Length[mat];
For[i=1,i<=len,i++,
s = ""; 
For[j=1,j<=9,j++,
If[mat[[i]][[j]]>0,s = s <> "+"];
If[mat[[i]][[j]]<0,s = s <> "-"];
abs=Abs[mat[[i]][[j]]];
If[Not[IntegerQ[abs]],abs=N[abs]];
If[abs!=0 && abs!=1, s=s <>ToString[abs]]; 
If[abs!=0, 
If[j==1, s = s <> "π"];
If[j==2, s = s <> "A"];
If[j==3, s = s <> "B"];
If[j==4, s = s <> "α"];
If[j==5, s = s <> "β"];
If[j==6, s = s <> "γ"];
If[j==7, s = s <> "α'"];
If[j==8, s = s <> "β'"];
If[j==9, s = s <> "γ'"]
]
];
Print[ToString[i]<>": "<>s]
]
]

(* example *)
showIneqs[{{3, 0, 1, -2, 1/4, 0, 0, -5, -1/2}}]

(* 
  Assuming all entries in specified column are -1, 0 or 1, 
  eliminate the corresponding variable (Fourier-Motzkin)
 to produce a new matrix.
*)
elim[mat_,index_]:= Module[ {i,j,len, newmat},
newmat={};
len = Length[mat] ;
For[i=1, i<=len, i++, 
If[ mat[[i]][[index]] == 0, AppendTo[newmat,mat[[i]]]]];
For[i=1, i<=len, i++, 
If[ mat[[i]][[index]] == 1, 
For[j=1, j<=len, j++, 
If[ mat[[j]][[index]] == -1, 
AppendTo[newmat,mat[[i]]+mat[[j]]]]]]];
newmat
]

(*
  Scale each row of the matrix so that the entry in the 
specified column is either -1, 0 or 1.
*)
normal[ mat_, index_] := Module[ {i,j,len,newmat}, 
newmat={};
len = Length[mat]; 
For[i=1,i⩽len,i++, 
If[ mat[[i]][[index]] == 0, 
AppendTo[newmat, mat[[i]]], 
AppendTo[newmat, mat[[i]]/Abs[mat[[i]][[index]]]]
]
];
newmat
]

(* 
  Scale each row of integers so that its GCD is 1;  
  Also, remove duplicate rows. 
*)
clean[mat_] := Module[ {i,j,len,len2,newmat} , 
newmat = {}; 
len = Length[mat]; 
For[i=1, i<=len, i++,
If[ And@@IntegerQ /@ mat[[i]], 
AppendTo[newmat, mat[[i]] / GCD @@ mat[[i]]], 
AppendTo[newmat, mat[[i]] ]
]; 
];
newmat = DeleteDuplicates[newmat]; 
newmat
]

(* 
  Used to convert final matrix to integers only,
  and to remove repeated rows.
*)
tweak[mat_] := Module[ {i,j,len,len2,newmat} , 
newmat = {}; 
len = Length[mat]; 
For[i=1, i<=len, i++,
If[ And@@IntegerQ /@ mat[[i]], 
AppendTo[newmat, mat[[i]]], 
AppendTo[newmat, 2 mat[[i]] ]
]; 
];
newmat = DeleteDuplicates[newmat]; 
newmat
]

(*
  Find and remove redundant inequality rows from a 
 matrix, using linear programming.
*)
reduce[mat_, startindex_] := Module[ 
{i,len,len0,newmat,lower,vector}, 
newmat = mat;
len0 = Length[mat];   
For[ i = len0, i ⩾ startindex, i--, 
len = Length[newmat]; 
lower = ConstantArray[0,len];
lower[[i]]= -1;
vector = LinearProgramming[newmat[[i]],newmat,lower, {{1,1},{0,1},{0,1},{0,1},{0,1},{0,1},{0,1},{0,1},{0,1}}];
If[ vector. newmat[[i]] == 0, newmat = Delete[newmat,i]]
];
newmat
]

(* 
  Matrix describing tetrahedron angles inequalities 
*)
M1 = {
{  1,  0,  0,  0,  0,  0,  0,  0,  0 }, 
{  0,  1,  0,  0,  0,  0,  0,  0,  0 }, 
{  0,  0,  1,  0,  0,  0,  0,  0,  0 }, 
{  0,  0,  0,  1,  0,  0,  0,  0,  0 }, 
{  0,  0,  0,  0,  1,  0,  0,  0,  0 }, 
{  0,  0,  0,  0,  0,  1,  0,  0,  0 }, 
{  0,  0,  0,  0,  0,  0,  1,  0,  0 }, 
{  0,  0,  0,  0,  0,  0,  0,  1,  0 }, 
{  0,  0,  0,  0,  0,  0,  0,  0,  1 }, 
{  1,-1,-1, 0,  0,  0, 0,  0,  0 }, 
{  1,  0,  0,-1,  0,  0,-1,  0,  0 }, 
{  1,  0,  0,  0,-1,  0,  0,-1,  0 }, 
{  1,  0,  0,  0,  0,-1,  0,  0,-1 }, 
{  2,  0,  0,-1,-1,-1,  0,  0,  0 }, 
{  0,  0,  0,-1,  1,  1,  0,  0,  0 }, 
{  0,  0,  0,  1,-1,  1,  0,  0,  0 }, 
{  0,  0,  0,  1,  1,-1,  0,  0,  0 }, 
{  1,-1,  0,  0,  0,  1,  0,-1,  1 }, 
{  1,-1,  0,  0,  0,-1,  0,  1,-1 }, 
{  1,  1,  0,  0,  0,-1,  0,-1,-1 }, 
{-1,  1,  0,  0,  0,  1,  0,  1,  1 }, 
{  1,  0,-1,  1,  0,  0,  1,  0,-1 }, 
{  1,  0,-1,-1,  0,  0,-1,  0,  1 },
{  1,  0,  1,-1,  0,  0,-1,  0,-1 }, 
{-1,  0,  1,  1,  0,  0,  1,  0,  1 }, 
{  0,  1,  1,  0,  1,  0,-1,  1,  0 }, 
{  0,  1,  1,  0,-1,  0,  1,-1,  0 }, 
{  2,-1,-1,  0,-1,  0,-1,-1,  0 }, 
{  0,-1,-1,  0,  1,  0,  1,  1,  0 }
} ; 
showIneqs[M1]

(* Eliminate γ' and normalize for β' *)
showIneqs[ M2 = clean[normal[elim[M1,9],8]]]

showIneqs[ M3 = reduce[M2,9]]

(* Eliminate β' and normalize for α' *)
showIneqs[ M4 = clean[normal[elim[M3,8],7]]]

showIneqs[ M5 = reduce[M4,8]]

(* Eliminate α' *)
showIneqs[ M6 = tweak[clean[elim[M5,7]]]]

showIneqs[ M7 = reduce[M6,7]]