(* ------------------------------------------------------------------------- *)
(*
Unitarize

Mathematica package for testing spherical polygon inequalities
and simultaneously unitarizing sets of matrices in SL_2(C).

Provides the following functions:

  PolygonInequalities[] -- check the spherical polygon inequalities
  Unitarizable[] -- check simultaneous unitarizablility of matrices
  Unitarizer[] -- compute simultaneous unitarizer of matrices
  SL2[] -- check whether a matrix is in SL_2(C)
  SU2[] -- check whether a matrix is in SU_2

See below for complete documentation.

Nick Schmitt
July 29 2002
*)
(* ------------------------------------------------------------------------- *)

BeginPackage[
  "Unitarize`",
    {"DiscreteMath`Combinatorica`", "LinearAlgebra`MatrixManipulation`" }
]

PolygonInequalities::usage=
"PolygonInequalities[value,...] returns True if the values satisfy the
spherical polygon inequalities and False otherwise.";

SL2::usage=
"SL2[M] returns True if the matrix M is in SL_2(C) and False otherwise."

SU2::usage=
"SU2[M] returns True if the matrix M is in SU_2 and False otherwise."

Unitarizable::usage=
"Unitarizable[M1,...] returns True if the matrices M1,... in SL_2(C)
are simultaneously unitarizable and False otherwise."

Unitarizer::usage=
"Unitarizable[M1,...] returns a simultaneous unitarizer T of the
matrices M1,... in SL_2(C). T is upper triangular and has positive diagonal.
If the matrices are not simultaneously unitarizable, False is returned."

Tolerance::usage=
"An option to functions in the Unitarize package. Tolerance->ep causes checks
for equality and inequality to be performed within tolerance ep."

Verbose::usage=
"Option to PolygonInequalities[]. The option Verbose->True causes
PolygonInequalities[] to print which spherical polygon inequality first
failed.";

CheckStrict::usage=
"Option to PolygonInequalities[]. The option CheckStrict->True causes
PolygonInequalities[] to return a symbol from the list
{Fail, Strict, Static} rather than returning {False, True}.";

Fail::usage=
"Return value of PolygonInequality[] when the option CheckStrict->True is
used. Fail is returned if at least one of the spherical polygon inequalities
fails."

Strict::usage=
"Return value of PolygonInequality[] when the option CheckStrict->True is
used. Strict is returned if all the spherical polygon inequalities are
strict inequalities."

Static::usage=
"Return value of PolygonInequality[] when the option CheckStrict->True is
used. Static is returned if at least one of the spherical polygon inequalities
is a strict equality."

(* 
The DiscreteMath`Combinatorica` package defines functions M[] and V[].
These are Remove[]'d so they can be used in other contexts.
*)
Unprotect[M];
Remove[M];
Unprotect[V];
Remove[V];

(* ------------------------------------------------------------------------- *)
Begin["`Private`"]

Unitarize::tolerance=
"Tolerance is negative."

(* ------------------------------------------------------------------------- *)
(* PolygonInequalities[] *)

Options[PolygonInequalities]=
{
  Verbose->False,
  CheckStrict->False,
  Tolerance->0
};

(* list version *)
PolygonInequalities[{v__?NumberQ}, opts___?OptionQ] :=
Module[{verbose, checkstrict, tolerance, retval},

  (* parse options *)
  {verbose, checkstrict, tolerance} = {Verbose, CheckStrict, Tolerance}/.
    Flatten[{opts, Options[PolygonInequalities]}];

  If[checkstrict && tolerance != 0, Message[PolygonInequalities::ep];Return[]];

  (* turn off messages if !verbose *)
  If[!verbose,
    Off[PolygonInequalities::count, PolygonInequalities::range,
      PolygonInequalities::fail]];

  (* compute polygon inequalities *)
  retval = Catch[polygonInequalities[{v},opts]];

  (* turn on messages if !verbose *)
  If[!verbose,
    On[PolygonInequalities::count, PolygonInequalities::range,
      PolygonInequalities::fail]];

  (* return *)
  retval = If[retval === Exception, Fail, retval];
  retval = If[checkstrict,retval,If[retval===Fail,False,True]];

  Return[retval];
];

(* sequence version *)
PolygonInequalities[v__?NumberQ, opts___?OptionQ] :=
  PolygonInequalities[{v},opts];


PolygonInequalities::ep =
  "CheckStrict->True and Tolerance is not 0.";
PolygonInequalities::count =
  "Not enough input values (must be at least two).";
PolygonInequalities::range =
  "Value `1` is not in the interval [0, 1/2].";
PolygonInequalities::fail =
  "Values failed the spherical polygon inequality `1`.";

(* private list version *)
polygonInequalities[{v0__?NumberQ},opts___?OptionQ] :=
Module[{v = {v0}, tolerance, n, S, P, A, B, s, static=False},

  n = Length[v];

  (* parse options *)
  {tolerance} = {Tolerance}/.
    Flatten[{opts, Options[PolygonInequalities]}];
  If[tolerance < 0, Message[Unitarize::tolerance];Throw[Exception]];

  (* check the number of values *)
  If[n < 2, Message[PolygonInequalities::count, n]; Throw[Exception]];
  
  (* check that the values lie in the interval [0,1/2] *)
  Do[
    If[! NumberQ[v[[i]]], Message[PolygonInequalities::num, i];
      Throw[Exception]];
    If[v[[i]] > 1/2 || v[[i]] < 0, Message[PolygonInequalities::range, i];
      Throw[Exception]];
    , {i, n}];
  
  (* check the spherical polygon inequalities *)
  S = Range[n];
  Do[
    P = KSubsets[S, k];
    Do[
      A = P[[i]];
      B = Complement[S, A];
      s = Sum[v[[A[[j]]]], {j, Length[A]}] - Sum[v[[B[[j]]]], {j, Length[B]}];
      s = s - (k-1)/2;
      If[ s > tolerance,
        s = StringJoin[
            Join[
              Table["+\[Nu]" <> ToString[A[[j]]], {j, Length[A]}],
              Table["-\[Nu]" <> ToString[B[[j]]], {j, Length[B]}]
              ], " \[LessSlantEqual] ", ToString[(k - 1)/2]];
        Message[PolygonInequalities::fail, s]; Throw[Exception]];
      If[tolerance==0, static = static || (s == 0)];
    , {i, Length[P]}];
  , {k, 1, n, 2}];

  Return[ If[static, Static, Strict] ];
];

(* ------------------------------------------------------------------------- *)
(* Unitarizable[] *)

Options[Unitarizable]=
{
  Tolerance->0
};

Unitarizable::sl2=
"Matrix is not in SL_2(C).";

(* single matrix version *)
Unitarizable[{M_?MatrixQ}, opts___?OptionQ] :=
Module[{r},
  r = Catch[unitarizable[{M},opts]];
  Return[If[r === Exception,Null,r]];
];

(* list of matrices version *)
Unitarizable[{M__?MatrixQ}, opts___?OptionQ] := MatrixQ[Unitarizer[{M},opts]]

(* sequence of matrices version *)
Unitarizable[M__?MatrixQ, opts___?OptionQ] := Unitarizable[{M},opts];

(* private single matrix version *)
unitarizable[{M_?MatrixQ}, opts___?OptionQ] :=
Module[{tolerance, t},

  (* parse options *)
  {tolerance} = {Tolerance}/.
    Flatten[{opts, Options[Unitarizable]}];
  If[tolerance < 0, Message[Unitarize::tolerance];Throw[Exception]];

  (* check that the matrix are in SL_2(C) *)
  If[!SL2[M,opts],Message[Unitarizable::sl2];Throw[Exception]];

  t = Chop[Simplify[Tr[M]/2],tolerance];
  If[Im[t] != 0, Return[False]];
  If[Abs[t] == 1, Return[M==IdentityMatrix[2] || M==-IdentityMatrix[2]]];
  Return[t > -1 && t < 1];
];

(* ------------------------------------------------------------------------- *)
(* Unitarizer[] *)

Options[Unitarizer]=
{
  Tolerance->0
};

Unitarizer::sl2=
"Matrix `1` is not in SL_2(C).";
Unitarizer::fail1=
"The dimension of the kernel is >1: not implemented";

(* list of matrices version *)
Unitarizer[{M__?MatrixQ}, opts___?OptionQ] :=
Module[{r},
  r = Catch[unitarizer[{M},opts]];
  Return[If[r === Exception,Null,r]];
];

(* sequence of matrices version *)
Unitarizer[M__?MatrixQ, opts___?OptionQ] := Unitarizer[{M},opts];

(* private list of matrices version *)
unitarizer[{M0__?MatrixQ}, opts___?OptionQ] :=
Module[{M={M0}, tolerance, n, f, L, nullity, X, det},

  (* parse options *)
  {tolerance} = {Tolerance}/.
    Flatten[{opts, Options[Unitarizer]}];
  If[tolerance < 0, Message[Unitarize::tolerance];Throw[Exception]];

  n = Length[M];
  
  (* check that the matrices are in SL_2(C) *)
  Do[If[!SL2[M[[i]],opts],Message[Unitarizer::sl2,i];Throw[Exception]],{i,n}];

  (* check that each matrix is unitarizable *)
  Do[If[!Unitarizable[M[[i]],opts],f=False],{i,n}];
  If[!f,Return[False]];

  L = linearMap[M];
  L = NullSpace[L];

  nullity = Length[L];
  If[nullity==0,Return[False]];

  (* the case when the kernel dimension is >1: NOT YET IMPLEMENTED *)
  If[nullity>1,
     Message[Unitarizer::fail1];Throw[Exception];
    ];

  (* the case when the kernel dimension == 1 *)
  X = unflatten[L[[1]]];
  X = Expand[(X + star[X])/2];

  If[X === ZeroMatrix[2],
     X = unflatten[L[[1]]];
     X = Expand[(I X + star[I X])/2];
    ];
  X = Chop[X,tolerance];
  If[X[[1,1]]<0,X=-X];

  (* compute unitarizer A such that X = A* A *)
  det = Chop[Simplify[Det[X]],tolerance];
  If[det<=0,Return[False]];
  X = {X[[1]],{0,Sqrt[Det[X]]}};

  (* make Det[X]==1 *)
  det = Chop[Simplify[Det[X]],tolerance];
  X = X/Sqrt[det];

  X = Chop[X,tolerance];
  Return[X];
];

(* ------------------------------------------------------------------------- *)

Options[SL2]=
{
  Tolerance->0
};

SL2[M_?MatrixQ, opts___?OptionQ] :=
Module[{tolerance},

  (* parse options *)
  {tolerance} = {Tolerance}/.
    Flatten[{opts, Options[SL2]}];
  If[tolerance < 0, Message[Unitarize::tolerance];Throw[Exception]];

  (* check whether M is in SL_2(C) *)
  Return[ Dimensions[M]==={2,2} && Chop[Det[M]-1,tolerance]==0 ];
]

(* ------------------------------------------------------------------------- *)

Options[SU2]=
{
  Tolerance->0
};

SU2[M_?MatrixQ, opts___?OptionQ] :=
Module[{tolerance},

  (* parse options *)
  {tolerance} = {Tolerance}/.
    Flatten[{opts, Options[SU2]}];

  (* check whether M is in SU_2(C) *)
  Return[ SL2[M,opts] &&
          Chop[Expand[star[M]-Inverse[M]],tolerance] == ZeroMatrix[2] ];
]

(* ------------------------------------------------------------------------- *)
(* private utility functions *)

toLinear[M_,V_]:=
  Transpose[Table[Coefficient[M,V[[i]]],{i,1,Length[V]}]];

linearMap[M_?MatrixQ] :=
Module[{X, X11, X12, X21, X22, Mx, Q},
  X = {{X11, X12}, {X21, X22}};
  Mx = star[Inverse[M]];
  Q = X.M - Mx.X;
  Return[toLinear[Flatten[Q], Flatten[X]]];
];

linearMap[M_List]:= Flatten[Map[linearMap[#]&,M],1];

star[A_?MatrixQ]:=Transpose[ComplexExpand[Conjugate[A]]];

unflatten[X_List]:={{X[[1]],X[[2]]},{X[[3]],X[[4]]}};


End[]
EndPackage[]