(*******************************************************************************)
(* QHull[pts,dim]: simplistic interface to QHull package

   Parameters:
     dim -- dimension of space (e.g., "2"),
     pts -- list of dim-dimensional points
       (e.g., "{{1,0},{0,1},{1,1},{0,0}}").
   Returns:
     list of convex hull faces, each represented by a
     list of point indices that are incident to the face
     (e.g., "{{3,2},{1,3},{2,4},{4,1}}").
   Note: necessary options from qhull manual:
   >
   > Fv - print vertices for each facet
   >
   > The first line is the number of facets. Then each facet is printed,
   > one per line. Each line is the number of vertices followed by the
   > corresponding point ids. Vertices are listed in the order they were
   > added to the hull (the last one added is the first listed).
*)

ClearAll[QHull];

Options[QHull] = {
  Verbose    -> False,
  Executable -> "./qhull",
  Mode       -> "Fv",
  InFile     -> "./_qhi.tmp",
  OutFile    -> "./_qho.tmp" };

QHull[pts_List, dim_Integer, ops___Rule] :=
  Module[{vb,ex,md,fi,fo,str,np,nf,fs,tm,i,s},

  vb = Verbose    /. {ops} /. Options[QHull];
  ex = Executable /. {ops} /. Options[QHull];
  md = Mode       /. {ops} /. Options[QHull];
  fi = InFile     /. {ops} /. Options[QHull];
  fo = OutFile    /. {ops} /. Options[QHull];

  np = Length[pts]; (* number of points *)

  (* write input file *)
  str = OpenWrite[fi];
  WriteString[str,ToString[dim] <> " # dimension\n"];
  WriteString[str,ToString[np] <> " # number of points\n"];
  For[i = 1, i <= np, i++,
    s = StringJoin[Riffle[ToString /@ pts[[i]],{" "}]];
    WriteString[str, s <> "\n"]];
  Close[str];

  (* run qhull program *)
  If[vb, Print["running qhull..."]];
  tm = Timing[Run[ex <> " " <> md <> " < " <> fi <> " > " <> fo]];
  If[vb, Print["done: ",tm[[1]]," seconds"]];

  (* read qhull output *)
  str = OpenRead[fo];
  nf = Read[str,Number]; (* number of facets *)
  fs = ReadList[str,Number,RecordLists->True];
  Close[str];

  (* adjust facet notations for Mathematica: *)
  (* - remove vertex counter (first element) from every entry *)
  (* - increment vertex indices (qhull assumes first index 0) *)
  Return[Drop[# + 1,1]& /@ fs]]

(*******************************************************************************)
(* UF[xx,yy,z]: routine by A.V.Smirnov to find polynomials U and F

   Source: http://www-ttp.particle.uni-karlsruhe.de/~asmirnov/Tools-UF.htm
   Parameters:
     xx -- list of loop momenta (e.g., "{v,q}")
     yy -- list of denominator factors
       (e.g., "{v^2 - mm, q^2 - mm, (v - q)^2}")
     z -- list of constraints (e.g., "{mm ->  1}").
   Returns:
     {polynomial U, polynomial F, number of loops}
     (e.g., "{x[1]*x[2]+x[1]*x[3]+x[2]*x[3],-x[1]^2*x[2]-x[1]*x[2]^2-
              x[1]^2*x[3]-2*x[1]*x[2]*x[3]-x[2]^2*x[3], 2}").
   Note: uses global symbol "x", as in: "x[1]", "x[2]", etc.
*)

ClearAll[UF,x];

UF[xx_, yy_, z_] := 
  Module[{degree, coeff, i, t2, t1, t0, zz},
  zz = Map[Rationalize[##, 0] &, z, {0, Infinity}];
  degree = -Sum[yy[[i]]*x[i], {i, 1, Length[yy]}];
  coeff = 1;
  For[i = 1, i <= Length[xx], i++,
    t2 = Coefficient[degree, xx[[i]], 2];
    t1 = Coefficient[degree, xx[[i]], 1];
    t0 = Coefficient[degree, xx[[i]], 0];
    coeff = coeff*t2;
    degree = Together[t0 - ((t1^2)/(4*t2))]];
  degree = Together[-coeff*degree] //. zz;
  coeff = Together[coeff] //. zz;
  {coeff, Expand[degree], Length[xx]}]

(*******************************************************************************)
(* PExpand[fp,xs,sm]: main entry, expansion regions of a polynomial fp
     of variables xs, with small parameter x
   Parameters:
     fp -- polynomial (e.g. "x[1]*x[2]*m^2*y^2 + x[2]^2*M^2")
     xs -- variables (e.g. "{x[1],x[2]}")
     sm -- small parameter (e.g., "y")
   Returns:
     list of n-dimensional vectors (n -- number of ds),
     each entry corresponding to one region.
*)

ClearAll[PExpand];

Options[PExpand] = {
  Verbose -> False,
  IntegralDim -> Default
};

PExpand[fp_,xs_,sm_,ops___Rule] :=
  Module[{vb,id,nx,fm,cc,rm,rs,pr,re,qb,vs,vv,sy,sl,x,i,v},
  vb = Verbose /. {ops} /. Options[PExpand];
  id = IntegralDim /. {ops} /. Options[PExpand];
  nx = Length[xs];
  If[id === Default, id = nx];

  If[vb,Print["integral dimension: ",id]];

  (* polynomial -> list of rules *)
  fm = CoefficientRules[fp,xs];
  (* find GCD of coefficients *)
  cc = PolynomialGCD @@ (#[[2]] & /@ fm);
  (* remove common factors, normalize, sort *)
  rm = Sort[Prepend[#[[1]], Simplify[#[[2]]/cc]]& /@ fm];
  (* coefficients -> scaling powers of small parameter *)
  rs = Prepend[Drop[#,1], -Exponent[Expand[#[[1]]] /. sm->1/x, x]]& /@ rm;

  If[vb,Print["points: ",rs]];
  
  If[vb,Print[Length[rs]," points"]];

  (* aux function to determine rank of a set of points *)
  PRank[ps_List] := If[Length[ps] <= 1, 0, 
    MatrixRank @ Map[# - ps[[1]]&, Drop[ps,1]]];

  pr = PRank[rs];
  If[vb,Print["dimensionality of the convex hull: ",pr]];

  If[id != pr - 1,
    Return[{}];
(* A.P.: workaround, useful to find homogeneity degree
    rp = Map[# - rs[[1]]&, Drop[rs,1]] // NullSpace;
    rp = Map[Drop[#,1]&,rp];
    If[vb,Print["scaleless polynomial: null space basis: ",rp]];
    Return[{0, rp}]; 
*)
    ];

  (* attempt to build projection (works for alpha-representation) *)
  rp = Drop[#,id - nx]& /@ rs;

  (* if this does not work, must introduce local variables, etc *)
  If[pr != PRank[rp],
    Print["projection does not work, alternative methods not implemented."];
    Return[{}]];

  qh = QHull[rp,pr,ops]; (* convex hull of projection *)
  (* If[vb,Print["convex hull: ",qh]]; *)
  If[vb,Print["convex hull has ",Length[qh]," facets"]];

  (* normal vector composition: independent variables *)
  vs = v[#]& /@ Range[id];
  (* (nx+1)-dimensional   *)
  vv = Join[{1},ConstantArray[0,nx-id],vs];
  (* If[vb,Print["generic normal vector: ",vv]]; *)

  For[i = 1; re = {}; qb = {}, i <= Length[qh], i++,
    sy = Join[
      # == v& /@ (rs[[qh[[i]]]].vv), (* facet normal, pointing up *)
      # >= v& /@ (rs.vv)]; (* other points lie above facet *)
    sl = FindInstance[sy,Append[vs,v]];
    If[Length[sl] != 0,
      re = Append[re, vv /. sl[[1]]];
      qb = Append[qb,i]]];

  If[vb,Print["bottom facets numbers: ", qb]];

  (* remove artificial zeroth component '1' *)
  Return[Drop[#,1]& /@ re]]

(*******************************************************************************)
(* AlphaRepExpand[ks,ds,cs,hi]: expansion in alpha-representation
   Parameters:
     ks -- list of loop momenta (e.g. "{v1}")
     ds -- list of denominators (e.g. "{v1^2 + m^2, (v1+p0)^2 + m^2}")
     cs -- kinematic constraints (e.g. "{p0^2 -> -M^2}")
     hi -- scalings wrt small parameter x (e.g. "{m -> m*x^1, M -> M*x^0 }")
   Returns:
   Note:
     uses global symbols x (as "x,x[1],x[2],...")
*)

ClearAll[AlphaRepExpand];

AlphaRepExpand[ks_List,ds_List,cs_List,hi_List,ops___Rule] :=
  Module[{up,fp,nl,xs,re},

  (* produce alpha-parametrization *)
  {up,fp,nl} = UF[ks,ds,cs];

  (* alpha-parameters *)
  xs = Table[x[i],{i,1,Length[ds]}];

  (* attempt to expand using only F *)
  re = PExpand[fp /. hi, xs, x, ops, IntegralDim -> Length[ds]-1];

  (* if scaleless, attempt to expand the product *)
  If[Length[re] == 0,
    Print["F itself is scaleless, trying U*F"];
    re = PExpand[up * fp /. hi, xs, x, ops, IntegralDim -> Length[ds]-1],
    Print["F scaleful, append trivial region"];
    re = Append[re, ConstantArray[0,Length[ds]]]];

  (* only unique regions *)
  re = Union[re];

  Return[re]]

(*******************************************************************************)

