// Regular models by Tim Dokchitser. Beta version, March 2018, please do not distribute. 
// Based on ``Models of curves over DVRs'', to be available soon.

freeze;   // tested by rmtest.m, rmtest3.m, test.m

///////////////////////
// 1. Common functions ///
///////////////////////

Z:=Integers();
Q:=Rationals();
PR:=PolynomialRing;
RFF:=RationalFunctionField;
Left:=func<s,n|S[[1..Min(#S,n)]] where S is Sprint(s)>;
Right:=func<s,n|S[[Max(1,#S-n+1)..#S]] where S is Sprint(s)>;
Count:=func<S,x|#[z: z in S | z eq x]>;   // count occurences of x in a list
DelCRs:=func<s|&cat([x: x in Eltseq(Sprint(s)) | x ne "\n"])>; // delete \n's
DelSpaces:=func<s|&cat([x: x in Eltseq(Sprint(s)) | (x ne " ") and (x ne "\n")])>;
RightCosetRepresentatives:=func<G,H|DoubleCosetRepresentatives(G,H,sub<G|>)>;
vec:=func<v|Type(v) eq ModTupFldElt select v else 
  (Type(v) eq SeqEnum select Vector(Q,v) else Vector(Q,Eltseq(v)))>;


function DelSymbols(s,delsymb) 
  if Type(s) ne MonStgElt then s:=Sprint(s); end if;
  if Type(delsymb) eq MonStgElt then delsymb:=Eltseq(delsymb); end if;  
  s:=[c: c in Eltseq(s) | not (c in delsymb)];
  return IsEmpty(s) select "" else &cat s;
end function;

function MLoCase(s)
  if (#s eq 0) or (s[1] eq ToLower(s[1])) then return s; end if;
  return ToLower(s[1]) cat s[[2..#s]];
end function;

function ReplaceStringFunc(s,fs,ts)
  if Type(s) ne MonStgElt then s:=Sprint(s); end if;
  if Type(fs) eq SeqEnum then
    for i:=1 to #fs do
      s:=ReplaceStringFunc(s,fs[i],ts[i]);
    end for;
    return s;
  end if;
  return SubstituteString(s,fs,ts);
end function;

procedure ReplaceString(~s,fs,ts)
  s:=ReplaceStringFunc(s,fs,ts);
end procedure;

function Last(v)
  if Type(v) in [SeqEnum,List,SetIndx,Tup] then return v[#v]; end if;
  if HasSignature("Eltseq",[Type(v)]) then 
    w:=Eltseq(v); return w[#w]; 
  end if;
  error "Last: unrecognised type";
end function;

function trim(s) 
  n:=#s;
  while (n ge 2) and (s[n] in [" ","\n"]) do n-:=1; end while;
  f:=1;
  while (f le #s) and (s[f] in [" ","\n"]) do f+:=1; end while;
  if (f eq 1) and (n eq #s) then return s; end if;
  return s[[f..n]];
end function;

procedure SortBy(~L,f: sign:=1)
  if Type(L) eq SetEnum then
    L:=SetToSequence(L);
  end if;
  Sort(~L,func<a,b|fa lt fb select -sign else (fa eq fb select 0 else sign)
                   where fa is f(a) where fb is f(b)>);
end procedure;                

function SortSet(X)           return Sort(SetToSequence(Set(X))); end function;

function SetAndMultiplicities(X)
  S:=SortSet(Set(X));
  M:=[Multiplicity(X,x): x in S];
  return S,M;
end function;

function VectorContent(v)
  den:=LCM([Denominator(x): x in v]);
  v:=[Z!(den*x): x in v];
  gcd:=GCD(v);
  if gcd eq 0 then return 0; end if;
  v:=[x div gcd: x in v];
  return gcd/den,v;
end function;

function PrintSequence(L: sep:=", ", mult:="^^%o", prfun:=Sprint)    
  if IsEmpty(L) then return ""; end if;
  if Type(L) eq SetMulti then 
    S,M:=SetAndMultiplicities(L);
    L:=[prfun(S[i])*(M[i] eq 1 select "" else Sprintf(mult,M[i])): i in [1..#S]];
    prfun:=Sprint;
  end if;
  return &cat[prfun(L[i]) * (i eq #L select "" else sep): i in [1..#L]];
end function;

function PrintPolynomial(cfs: Reverse:=false, Var:="x")
  s:="";
  for step:=0 to #cfs-1 do
    i:=Reverse select step else #cfs-1-step;
    c:=Sprint(cfs[i+1]);
    
    if c eq "0" then continue; end if;

    if i eq 0 then ;
      elif i eq 1 then c cat:="*"*Var;
      else c cat:= Sprintf("*%o^%o",Var,i);
    end if;
    if Left(c,2) eq "1*" then c:=c[[3..#c]]; end if;
    if Left(c,3) eq "-1*" then c:="-"*c[[4..#c]]; end if;
    if (c[[1]] in ["-","+"]) and (#s ne 0) and (s[[#s]] eq "+") then
      s:=Left(s,#s-1);
    end if;
    s cat:= c cat "+";
  end for;
  s:=Left(s,#s-1);
  if #s eq 0 then s:="0"; end if;
  return s;
end function;  

function Express(L,x: R:=Z)
  return Eltseq(Solution(Matrix(R,L),Vector(R,x)));
end function;

function TeXPolynomial(f: vnames:="default", vorder:="default", cfprint:=Sprint)
  
  if Type(f) in [FldFunRatUElt,FldFunRatMElt] then
    num:=TeXPolynomial(Numerator(f): vnames:=vnames, vorder:=vorder, cfprint:=cfprint);
    den:=TeXPolynomial(Denominator(f): vnames:=vnames, vorder:=vorder, cfprint:=cfprint);
    if den eq "1" then return num; end if;
    if exists{c: c in Eltseq(num) | c in ["*","+","/","-"]} then num:="("*num*")"; end if;
    if exists{c: c in Eltseq(den) | c in ["*","+","/","-"]} then num:="("*den*")"; end if;
    return Sprintf("%o/%o",num,den);
  end if;
          
  R<[v]>:=Parent(f);
  if vorder cmpeq "default" then 
    vorder:=[1..#v];
  end if;
  if vnames cmpeq "default" then 
    vnames:=[Sprint(R.i): i in [1..#v]];
  end if;
  s:="";
  x:=vorder[1];
  for i:=Degree(f,x) to 0 by -1 do
    c:=Coefficient(f,x,i);
    if c eq 0 then continue; end if;
    if #vorder eq 1 
      then cs:=cfprint(c);
      else cs:=TeXPolynomial(c: vnames:=vnames, vorder:=vorder[[2..#vorder]], cfprint:=cfprint);
    end if;
    if (Position(cs,"+") ne 0 or Position(cs[[2..#cs]],"-") ne 0) and (i gt 0) then
      cs:="("*cs*")";
    end if;
    if cs eq "1" then cs:=""; end if;
    if cs eq "-1" then cs:="-"; end if;
    if i ne 0 then cs*:=vnames[x]; end if;
    lbr:=i gt 9 or i lt 0 select "{" else "";
    rbr:=i gt 9 or i lt 0 select "}" else "";
    if i gt 1 then cs*:="^"*lbr*Sprint(i)*rbr; end if;
    if cs eq "" then cs:="1"; end if;
    if cs eq "-" then cs:="-1"; end if;
    if s eq "" then s:=cs; else s*:=(cs[1] eq "-" select "" else "+")*cs; end if;
  end for;
  return s;
end function;

function Match(s,m)
  R:=RealField();
  s:=Eltseq(s cat CodeToString(127));
  m:=Eltseq(m cat CodeToString(127));

  // Convert "HEL%%LO %S(%I)" into [* "HEL%LO ","%S","(","%I",")" *]

  M:=[""];
  skip:=false;
  for i:=1 to #m do
    c:=m[i];
    if c ne "%" then
      l:=#(M[#M]);
      cmd:=(l ne 0) and (M[#M][1] eq "%");
      if cmd and (l eq 1) and (StringToCode(c) in [StringToCode("a")..StringToCode("z")]) then
        c:=CodeToString(StringToCode(c)-StringToCode("a")+StringToCode("A"));
      end if;
      M[#M] cat:= c;
      if not cmd then continue; end if;
      if (c in ["O","M"]) and (l eq 1) then continue; end if;
      Append(~M,"");
      continue;
    end if;
    if skip then skip:=false; continue; end if;
    if m[i+1] eq "%" then
      M[#M] cat:= c; skip:=true; continue;
    end if;
    Append(~M,"%");
  end for;
  M:=[s: s in M | #s ne 0];

  buf:=[* *];
  for i:=1 to #M do
    c:=M[i];
    if (#c eq 1) or (c[1] ne "%") then
      if (#s lt #c) or (&cat(s[1..#c]) ne c) then return false,buf; end if;
      for j:=1 to #c do Remove(~s,1); end for;
    elif c in ["%S","%E"] then
      if i eq #M then
        Append(~buf,&cat s); 
        if c eq "%E" then try buf[#buf]:=eval buf[#buf]; catch e return false, buf; end try; end if;
        return true,buf;
      end if;
      stop:=M[i+1];
      l:=0;
      if (#stop eq 3) and (stop[1] eq "%") and (stop[2] in ["M","O"]) then
        stop:=stop[3];
      end if;
      if (#stop eq 1) or (stop[1] ne "%") then
        for j:=1 to #s do
          if &cat(s[j..Min(#s,j+#stop-1)]) eq stop then l:=j; break; end if;
        end for;
      elif stop in ["%W","%I","%R"] then
        for j:=1 to #s do
          if (s[j] in ["0","1","2","3","4","5","6","7","8","9"])
             or ((stop ne "%W") and (s[j] eq "-"))
           then l:=j; break; end if;
        end for;
      else
        l:=1;
      end if;
      if l eq 0 then return false,buf; end if;
      Append(~buf,"");
      for j:=1 to l-1 do
        buf[#buf] cat:= s[1];
        Remove(~s,1);
      end for;
      if c eq "%E" then try buf[#buf]:=eval buf[#buf]; catch e return false, buf; end try; end if;
    elif (c[2] in ["M","O"]) and (#c eq 3) then
      ok:=false;
      while s[1] eq c[3] do Remove(~s,1); ok:=true; end while;
      if (c[2] eq "M") and (not ok) then return false,buf; end if;
    elif c[2] in ["W","I","R"] then
      l:=0;
      W:=s[1];
      if (W eq "-") and (c[2] eq "W") then return false,buf; end if;
      ok:=true;
      repeat
        try
          assert Last(W) ne " ";
          w:=eval W; 
          W cat:= s[l+2];
          error if not IsCoercible(c[2] eq "R" select R else Q,w),"";
          res:=w;
          l+:=1; 
        catch e
          ok:=false;
        end try;
      until not ok;
      if l eq 0 then return false,buf; end if;
      Append(~buf,res);
      for j:=1 to l do Remove(~s,1); end for;
    else
      error "Match: "*c*" is not implemented";
    end if;
  end for;
  return #s eq 0,buf;
end function;


// output

intrinsic Rewrite(filename::MonStgElt)
{Open file for writing, clearing its contents}
  Flush(Open(filename,"w"));
end intrinsic;


intrinsic write(filename::MonStgElt, str: con:=true, rewrite:=false, newline:="\n")
{Write a string to a file, and flush the results. By default (con:=true) prints on screen as well}
  if con then str; end if;
  F:=Open(filename,rewrite select "w" else "a");
  WriteBytes(F,[StringToCode(c): c in Eltseq(Sprint(str)*newline)]);
  Flush(F);
end intrinsic;

intrinsic writeq(filename::MonStgElt, str: rewrite:=false, newline:="\n")
{Quiet write a string to a file, and flush the results. (Quiet means does not print on the screen.)}
  write(filename,str: con:=false, rewrite:=rewrite, newline:=newline);
end intrinsic;



///////////////////////
// 2. Drawing /////////
///////////////////////


declare type TBox;                   // TeX Box

declare attributes TBox: 
  c,parent,     // children, parent
  banchors,     // [x1,x2,...] anchors on the bottom
  tanchors,     // anchors on the top, possibly empty
  width,height, // width and height
  left,right,   // left, right edges = [[x,y1,y2],[x,y1,y2],...] set of segments
  // all of the above are measured from 0,0
  x1,y1,        // where currently placed in the parent
  tex,          // function(B) -> draw everything necessary except children, returns tikz string
  // additional data for components
  style;


Line1Margins:=<2/5,5/5,1/2,1/2>;     // left, right, top, bottom
HLine2Margins:=<2/5,2/5,1/3,1/3>;
VLine2Margins:=<2/5,2/5,1/3,1/3>;
P1LineLength:=2/3;
RedBulletWidth:=1/4;
RedBulletHeight:=1/4;
RedBulletMargins:=<1/8,1/8,1/8,1/8>;
AvoidMinXSpace:=3/5;


intrinsic IsCoercible(M::TBox, y::.) -> BoolElt, .     {Coerce a TeX Box.}
  return false, _; end intrinsic;
intrinsic 'in'(x::TBox, y::.) -> BoolElt               {"in" function for a TeX Box.}
  return false; end intrinsic;


intrinsic TBoxCreate() -> TBox                         {Create a TeX Box (hackobj function).}
  B:=New(TBox); 
  B`c:=[]; B`left:=[]; B`right:=[]; 
  B`banchors:=[0]; B`tanchors:=[]; 
  B`parent:=false;
  B`x1:=0; B`y1:=0; B`width:=0; B`height:=0;
  B`tex:=func<B|"">;
  B`style:=<>;
  return B;
end intrinsic;


intrinsic Print(B::TBox, level::MonStgElt)             {Print a TeX Box.}
  printf "TBox %o %ox%o %o->%o",#B`c,B`width,B`height,DelSpaces(B`banchors),DelSpaces(B`tanchors);
end intrinsic;


intrinsic Box(width, height, banchors::SeqEnum, tanchors::SeqEnum, margins::Tup, tex::UserProgram) -> TBox
{Create a TBox encompassing a TeX primitive}
  B:=TBoxCreate();
  leftmargin,rightmargin,topmargin,bottommargin:=Explode(margins);
  B`left:=[[-leftmargin,-bottommargin,height+topmargin]];
  B`right:=[[width+rightmargin,-bottommargin,height+topmargin]];
  B`banchors:=banchors;
  B`tanchors:=tanchors;
  B`tex:=tex;
  B`width:=width;
  B`height:=height;
  return B;
end intrinsic;


intrinsic Duplicate(b::TBox) -> TBox
{Duplicate a box with all its children}
  B:=Box(b`width,b`height,b`banchors,b`tanchors,<0,0,0,0>,b`tex);
  B`x1:=b`x1;
  B`y1:=b`y1;
  B`left:=b`left;
  B`right:=b`right;
  B`parent:=b`parent;
  B`style:=b`style;
  B`c:=[Duplicate(ch): ch in b`c];
  for i:=1 to #B`c do 
    B`c[i]`parent:=B;
  end for;
  return B;
end intrinsic;


function RefineMargin(S,mm)   // mm = "Min" or "Max",
  assert mm in ["Min","Max"];
  f:=mm eq "Min" select Min else Max;
  Ys:=SortSet(&join[{d[2],d[3]}: d in S]);
  Xs:=[];
  
  // compute breakpoints
  for d in S do
    x,y1,y2:=Explode(d);
    i1:=Position(Ys,y1);
    i2:=Position(Ys,y2);
    for i:=i1 to i2 do
      if IsDefined(Xs,i) then Xs[i]:=f(Xs[i],x); else Xs[i]:=x; end if;
    end for;
  end for;

  // convert to intervals - currently not correct, too greedy for chains
  breaks:=[<Xs[1],Ys[1]>];
  for i:=2 to #Xs-1 do
    if mm eq "Max" then 
      if Xs[i] gt Xs[i-1] then Append(~breaks,<Xs[i],Ys[i]>); end if;
      if Xs[i] gt Xs[i+1] then Append(~breaks,<Xs[i+1],Ys[i]>); end if;
    else
      if Xs[i] lt Xs[i-1] then Append(~breaks,<Xs[i],Ys[i]>); end if;
      if Xs[i] lt Xs[i+1] then Append(~breaks,<Xs[i+1],Ys[i]>); end if;
    end if;
  end for;
  Append(~breaks,<0,Last(Ys)>);    
  return [[breaks[i][1],breaks[i][2],breaks[i+1][2]]: i in [1..#breaks-1]];
end function;


function LeftMarginWRXY(B)
  x:=B`x1;
  y:=B`y1;
  return [[d[1]+x,d[2]+y,d[3]+y]: d in B`left];
end function;


function RightMarginWRXY(B)
  x:=B`x1;
  y:=B`y1;
  return [[d[1]+x,d[2]+y,d[3]+y]: d in B`right];
end function;


intrinsic DimensionsFromChildren(~B::TBox)
{Compute dimensions and margins of a box from its children}
  c:=B`c;
  // box size
  B`width:=Max([b`x1+b`width: b in c]);
  B`height:=Max([b`y1+b`height: b in c]);
  // margins
  B`left:=RefineMargin(&cat [LeftMarginWRXY(b): b in c],"Min");
  B`right:=RefineMargin(&cat [RightMarginWRXY(b): b in c],"Max");
end intrinsic;


intrinsic Box(c::SeqEnum[TBox]) -> TBox
{Create TBox from a sequence of children}

  B:=TBoxCreate();

  // translate the children so that the box starts at 0,0; and set parent
  minx1:=Min([b`x1: b in c]);
  miny1:=Min([b`y1: b in c]);
  for i:=1 to #c do
    c[i]`x1-:=minx1; c[i]`y1-:=miny1; c[i]`parent:=B;
  end for;
  B`c:=c;
  DimensionsFromChildren(~B);
  return B;
end intrinsic;


function PrintReal(x)
  s:=Sprint(RealField()!x);
  p:=Position(s,".");
  assert p ne 0;
  return s[[1..Min(#s,p+2)]];
end function;


intrinsic AbsX(B::TBox) -> FldRatElt
{Absolute x1 position of a box}
  x:=0;
  while Type(B) eq TBox do
    x:=x+B`x1;
    B:=B`parent;
  end while;
  return x;
end intrinsic;


intrinsic AbsY(B::TBox) -> FldRatElt
{Absolute y1 position of a box}
  y:=0;
  while Type(B) eq TBox do
    y:=y+B`y1;
    B:=B`parent;
  end while;
  return y;
end intrinsic;


function TeXPosition(B,x,y)
  return Sprintf("(%o,%o)",PrintReal(x+AbsX(B)),PrintReal(y+AbsY(B)));
end function;


/*
    ok,buf:=Match(d,"%s at %s");
    at:="";
    if ok then d:=buf[1]; at:=buf[2]; end if;
*/

function HLineTeX(B)
  color,label,style,labelpos,v,source:=Explode(B`style);
  color:=#color eq 0 select "" else ","*color;
  main:=style eq "l1";
  if label ne "" then label:=Sprintf("node[%ol] {%o} ",labelpos,label); end if;
  if source ne "" then 
    ok,buf:=Match(source,"%s@%s");
    assert ok;
    srcstyle:=trim(buf[2]); 
    source:=Sprintf("node[facel%o%o] {%o}",#srcstyle eq 0 select "" else ",",
      srcstyle,trim(buf[1])); 
  end if;
  if v lt 0 then
    x1:=TeXPosition(B,B`width,0);
    x2:=TeXPosition(B,0,B`height);
  else
    x1:=TeXPosition(B,0,0);
    x2:=TeXPosition(B,B`width,B`height);
  end if;
  if labelpos eq "main"
    // below right of line
    then return Sprintf("\\draw[%o%o] %o--%o %o%o;\n",style,color,x1,x2,label,source);
    else return Sprintf("\\draw[%o%o] %o--%o%o%o;\n",style,color,x1,label,x2,source);
    // left of line
    //then return Sprintf("\\draw[%o%o] %o%o--%o %o;\n",style,color,x1,source,x2,label);
    //else return Sprintf("\\draw[%o%o] %o%o--%o%o;\n",style,color,x1,source,label,x2);
  end if;
end function;


function VLineTeX(B)
  color,label,style,labelpos:=Explode(B`style);
  color:=#color eq 0 select "" else ","*color;
  if label ne "" then label:=Sprintf("node[%ol] {%o} ",labelpos,label); end if;
  return Sprintf("\\draw[%o%o] %o--%o%o;\n",style,color, 
    TeXPosition(B,0,0),label,TeXPosition(B,0,B`height));
end function;


intrinsic HLine(width, label::MonStgElt, labelpos::MonStgElt, color::MonStgElt, style::MonStgElt, source::MonStgElt: vshift:=0, left:=false) -> TBox
{TBox with a horizontal (unless vshift>0 or <0 is given) line and given initial width, color and style}
  if Left(style,5) eq "l1" then // main component 
    margins:=Line1Margins;
  elif Left(style,5) eq "l2" then // chain of P1s
    margins:=HLine2Margins;
  else 
    error "HLine: Unknown line style (for margins)";
  end if;
  banchors:=left select [width] else [0];
  tanchors:=left select [0] else [width];
  B:=Box(width,Abs(vshift),banchors,tanchors,margins,HLineTeX);
  B`style:=<color,label,style,labelpos,Sign(vshift),source>;
  return B;
end intrinsic;


intrinsic VLine(height, label::MonStgElt, labelpos::MonStgElt, color::MonStgElt, style::MonStgElt) -> TBox
{TBox with a vertical line and given initial height}
  B:=Box(0,height,[0],[0],VLine2Margins,VLineTeX);
  B`style:=<color,label,style,labelpos>;
  return B;
end intrinsic;


function RedBulletTeX(B: style:="default", tex:="default")
  if style eq "default" then style:="[red]"; end if;
  if tex eq "default" then tex:="$\\bullet$"; end if;
  return Sprintf("\\node%o at %o {%o};\n",style,TeXPosition(B,0,0),tex);
end function;


intrinsic RedBullet() -> TBox
{TBox with a red bullet for a singular point}
   return Box(RedBulletWidth, RedBulletHeight, [0], [], RedBulletMargins, RedBulletTeX);
end intrinsic;


intrinsic SpecialPoint(style::MonStgElt, tex::MonStgElt: width:=RedBulletWidth, 
  height:=RedBulletHeight, margins:=RedBulletMargins) -> TBox
{TBox with a special point (e.g. for a singular point); defaults 
style="default", tex="default" give a standard red bullet}
  return Box(width,height,[0],[],margins,func<B|RedBulletTeX(B: style:=style, tex:=tex)>);
end intrinsic;


function SelfP1Chain(m,color,style)

  if #m eq 0 then                                      // 0 : loop
    if style eq "qmark" then
      s:="\\path[draw,thick%o] %o++(0,-0.1)--(-0.2,0.35)++(0.0,-0.2)--(0.0,0.6)++(-0.15,-0.1)
   --(0.025,0.52)++(0.075,-0.44) node[qmark]{} ++(0.075,0.44)--(0.35,0.51)++(-0.15,0.09)
   --(0.4,0.15)++(0.0,0.2)--(0.2,-0.1)++(-0.15,0.63)--(0.07,0.53)++(0.02,0.0)
   --(0.11,0.53)++(0.02,0.0)--(0.15,0.53);";
    else  
      s:="\\path[draw,thick%o] %o edge[white,line width=2] ++(0.6,0) 
          edge[out=0,in=-90] ++(0.6,0.4) ++(0.6,0.4) edge[out=90,in=0] ++(-0.3,0.3) ++(-0.3,0.3)
          edge[out=180,in=90] ++(-0.3,-0.3) ++(-0.3,-0.3) edge[out=-90,in=180] ++(0.6,-0.4);";
    end if;
    color:=color eq "" select "" else ","*color;
    tex:=func<B|Sprintf(s,color,TeXPosition(B,0,0))>;
    B:=Box(0.6, 1, [0], [], <1/3,1/3,1/3,1/3>, tex);
    B`banchors:=[0,0.6];
    return B;
  end if;   

  if #m eq 1 then                                      // 1: arc component
    s:="\\draw[-,shorten <=-0.5em,shorten >=-0.5em%o] %o ++(1,0) arc (0:180:0.5) node[midway,abovel] {%o};";
    color:=color eq "" select "" else ","*color;
    tex:=func<B|Sprintf(s,color,TeXPosition(B,0,0),m[1])>;
    B:=Box(1, 1, [0], [], <1/4,1/4,1/4,1/4>, tex);
    B`banchors:=[0,1];
    return B;
  end if;    
   
  r8:=#m mod 8;
  d8:=#m div 8 - (#m mod 8 in [4,5,6,7] select 0 else 1);
  URUL:=&cat[["U","R","U","L"]: i in [1..d8]];
  ULUR:=&cat[["U","L","U","R"]: i in [1..d8]];
  LDRD:=&cat[["L","D","R","D"]: i in [1..d8]];
  RDLD:=&cat[["R","D","L","D"]: i in [1..d8]];

  if #m in [2..9] then dir:=
    [["LSU","LSD"],                                     // 2 : chains
     ["U","R","D"],                                     // 3
     ["LSU","RU","RD","LSD"],                           // 4
     ["LSU","RSU","R","RSDA","LSD"],                    // 5
     ["LSU","U","RU","RD","D","LSD"],                   // 6
     ["LSU","U","RU","R","RD","D","LSD"],               // 7
     ["LSU","RU","RSD","RSU","RSD","RSUB","RD","LSD"],  // 8
     ["U","LSU","U","RSU","R","RSDA","D","LSD","D"]]    // 9
   [#m-1];
  elif r8 eq 0 then
    dir:=URUL cat ["U","RSU","RSD","RSU","RSD","RSU","RSDA","D"] cat LDRD;
  elif r8 eq 1 then
    dir:=URUL cat ["LSU","RSU","RSD","RSU","R","RSD","RSU","RSD","LSD"] cat LDRD;
  elif r8 eq 2 then
    dir:=ULUR cat ["U","L","U","R","RSU","RSDA","R","D","L","D"] cat RDLD;
  elif r8 eq 3 then
    dir:=ULUR cat ["U","L","U","R","U","R","D","R","D","L","D"] cat RDLD;
  elif r8 eq 4 then
    dir:=["RSU"] cat ULUR cat ["RSU","RSDA"] cat RDLD cat ["RSDA"];
  elif r8 eq 5 then
    dir:=["LSU"] cat ULUR cat ["RSU","R","RSDA"] cat RDLD cat ["LSD"];
  elif r8 eq 6 then
    dir:=URUL cat ["U","R","RSU","RSDA","R","D"] cat LDRD;
  else
    dir:=URUL cat ["U","R","U","R","D","R","D"] cat LDRD;
  end if;

  x:=0;
  y:=0;
  c:=[];

  for i:=1 to #m do
    d:=dir[i];
    l:=P1LineLength;
    l2:=l/2;

    if   d eq "L"   then bx1:=x-l;  by1:=y;    x-:=l;          vshift:=0;   label:="above";
    elif d eq "U"   then bx1:=x;    by1:=y;    y+:=l;          vshift:=0;   label:="left";
    elif d eq "D"   then bx1:=x;    by1:=y-l;  y-:=l;          vshift:=0;   label:="right";
    elif d eq "R"   then bx1:=x;    by1:=y;    x+:=l;          vshift:=0;   label:="above";
    elif d eq "RU"  then bx1:=x;    by1:=y;    x+:=l;  y+:=l;  vshift:=l;   label:="aboveleft";
    elif d eq "RD"  then bx1:=x;    by1:=y-l;  x+:=l;  y-:=l;  vshift:=-l;  label:="aboveright";
    elif d eq "RSU" then bx1:=x;    by1:=y;    x+:=l2; y+:=l;  vshift:=l;   label:="aboveleft"; l:=l2;
    elif d eq "RSUB" then bx1:=x;   by1:=y;    x+:=l2; y+:=l;  vshift:=l;   label:="belowright"; l:=l2;
    elif d eq "RSD" then bx1:=x;    by1:=y-l;  x+:=l2; y-:=l;  vshift:=-l;  label:="belowleft"; l:=l2;
    elif d eq "RSDA" then bx1:=x;   by1:=y-l;  x+:=l2; y-:=l;  vshift:=-l;  label:="aboveright"; l:=l2;
    elif d eq "LU"  then bx1:=x-l;  by1:=y;    x-:=l;  y+:=l;  vshift:=-l;  label:="belowleft"; 
    elif d eq "LD"  then bx1:=x-l;  by1:=y-l;  x-:=l;  y-:=l;  vshift:=l;   label:="belowright";  
    elif d eq "LSU" then bx1:=x-l2; by1:=y;    x-:=l2; y+:=l;  vshift:=-l;  label:="aboveright"; l:=l2;
    elif d eq "LSD" then bx1:=x-l2; by1:=y-l;  x-:=l2; y-:=l;  vshift:=l;   label:="aboveleft";  l:=l2;
    else error "Unrecognized direction "*d; 
    end if;

    b:=d[1] in ["U","D"] 
       select VLine(l,Sprint(m[i]),label,color,style/*: down:=d[1] eq "D"*/)
       else   HLine(l,Sprint(m[i]),label,color,style,"": vshift:=vshift, left:=d in ["L","RD","LU"]);
    b`x1:=bx1;
    b`y1:=by1;

    Append(~c,b);

  end for;

  B:=Box(c);
  if #m eq 2 
    then B`banchors:=[c[#c]`x1+c[#c]`banchors[1],c[1]`x1+c[1]`banchors[1]];
    else B`banchors:=[c[1]`x1+c[1]`banchors[1],c[#c]`x1+c[#c]`banchors[1]];
  end if;
  return B;
end function;


intrinsic Box(m::SeqEnum, linksup::BoolElt: color:="", style:="", self:=false) -> TBox
{TBox from chain of multiplicities}
  if style eq "" then style:="l2"; end if;
  if self then return
    SelfP1Chain(m,color,style);
  end if;
  if IsEmpty(m) then
    vprint Newton,2: "Box: empty P1 chain (contractible link)";
    return VLine(P1LineLength,"","right","bluedash",style);
  end if;

  dir:=["U","L","U","R"];
  x:=0;
  y:=0;
  c:=[];
  l:=P1LineLength;
  v:=0; 
  hlabel:="above";
  for i:=1 to #m do
    r:=(i-1) mod 4 + 1;

    d:=dir[r];
    if d eq "L" and linksup and i eq #m then v:=-l; hlabel:="belowleft"; end if;
    if d eq "R" and linksup and i eq #m then v:=l; hlabel:="belowright"; end if;
    vlabel:=r eq 3 select "left" else "right";
    b:=d eq "U" select VLine(l,Sprint(m[i]),vlabel,color,style)
                else   HLine(l,Sprint(m[i]),hlabel,color,style,"": vshift:=v, left:=d eq "L");
    if   d eq "U" then b`x1:=x;   b`y1:=y;  y+:=l; 
    elif d eq "R" then b`x1:=x;   b`y1:=y;  x+:=l; 
    else               b`x1:=x-l; b`y1:=y;  x-:=l; 
    end if;
    Append(~c,b);

  end for;
  B:=Box(c);
  B`banchors:=[c[1]`x1+c[1]`banchors[1]];
  if linksup then 
    B`tanchors:=[c[#c]`x1+c[#c]`tanchors[1]];
  end if;
  return B;
end intrinsic;


intrinsic TeX(B::TBox) -> MonStgElt
{TeX a TBox with all its children}
  s:=B`tex(B);        // Box itself
  for b in B`c do
    s*:=TeX(b);     // followed by all the children
  end for;
  return s;
end intrinsic;


function VerticalOverlap(ya1,ya2,yb1,yb2)
  if ya1 ge yb2 then return false; end if;
  if yb1 ge ya2 then return false; end if;
  return true; 
end function;


intrinsic Place(~B::TBox, U::TBox, cmds: first:=false, xpos:=0, ypos:=0) 
{Place a new child box U into a box B; commands [<"R",B1>,<"A2",C>,...]; 
first prepends otherwise appends}

  B`c:=first select [U] cat B`c else B`c cat [U];
  U`parent:=B;
  U`x1:=xpos;
  U`y1:=ypos;
  
  for c in cmds do 

    if c[1] eq "R" then       // "R": right of a node    
      F:=c[2];
      shift:=0;
      for dF in F`right, dU in F`left do
        if VerticalOverlap(F`y1+dF[2],F`y1+dF[3],U`y1+dU[2],U`y1+dU[3]) then 
          shift:=Max(shift,(F`x1+dF[1])-(U`x1+dU[1]));
        end if;
      end for;
      U`x1+:=shift;
    elif c[1] eq "A2" then    // "A": avoid links up from a given chain
      avoid:=[AbsX(b)+t: t in b`tanchors, b in c[2]];
      x:=AbsX(U)+U`banchors[1];

      vprintf Newton,2: "A: %o from %o",x,DelSpaces(avoid);
      while exists(xa){xa: xa in avoid | Abs(xa-x) lt AvoidMinXSpace} do
        U`x1+:=Abs(xa-x)+AvoidMinXSpace;
        x:=AbsX(U);
      end while;        
    else
      error Sprintf("Place: %o is not implemented",c[1]);    
    end if;
  end for;

  DimensionsFromChildren(~B);
end intrinsic;


procedure PlaceBoxesAtHeight(~B, H, C, avoid)
  for c in C do
    if Type(c) ne TBox then continue; end if;
    Place(~B,c,[*<"R",B>,<"A2",[b: b in avoid | Type(b) eq TBox]>*]: ypos:=H);
  end for;
end procedure;


intrinsic VerticalScale(~B::TBox, r)
{Rescale a box vertically by a factor of r}
  if r eq 1 then return; end if;
  B`y1*:=r;                                // position in parent and total height
  B`height*:=r;                            
  for i:=1 to #B`c do                      // children
    VerticalScale(~B`c[i],r);
  end for;
  for i:=1 to #B`left do                   // left margin
    x,y1,y2:=Explode(B`left[i]);
    B`left[i]:=[x,r*y1,r*y2];
  end for;
  for i:=1 to #B`right do                  // right margin
    x,y1,y2:=Explode(B`right[i]);
    B`right[i]:=[x,r*y1,r*y2];
  end for;
end intrinsic;


intrinsic VerticalScaleTo(~B::TBox, H)
{Rescale a box vertically to height H}
  if B`height in [0,H] then return; end if;
  VerticalScale(~B,H/B`height);
end intrinsic;


intrinsic VerticalScaleTo(~C::SeqEnum[TBox], H)
{Rescale sequence of boxes vertically to height H, not touching ones with no top anchors}
  for i:=1 to #C do
    if IsEmpty(C[i]`tanchors) then continue; end if;
    VerticalScaleTo(~C[i],H);
  end for;
end intrinsic;


intrinsic VerticalScaleTo(~C::List, H)
{Rescale sequence of boxes vertically to height H, not touching ones with no top anchors}
  S:=[c: c in C | Type(c) eq TBox];
  VerticalScaleTo(~S,H);
end intrinsic;


intrinsic HorizontalScale(~B::TBox, r)
{Rescale a box horizontally by a factor of r}
  if r eq 1 then return; end if;
  B`x1*:=r;                                // position in parent and total height
  B`width*:=r;                            
  B`banchors:=[r*x: x in B`banchors];
  B`tanchors:=[r*x: x in B`tanchors];
  for i:=1 to #B`c do                      // children
    HorizontalScale(~B`c[i],r);
  end for;
  for i:=1 to #B`left do                   // left margin
    x,y1,y2:=Explode(B`left[i]);
    B`left[i]:=[r*x,y1,y2];
  end for;
  for i:=1 to #B`right do                  // right margin
    x,y1,y2:=Explode(B`right[i]);
    B`right[i]:=[r*x,y1,r*y2];
  end for;
end intrinsic;


function VertexData(ch1)
  singular:=false;
  mult:=1;
  genus:=0;
  source:="";
  for c in [* c: c in ch1 | Type(c) eq Tup *] do     // vertex commands
    if   c[1] eq "singular" then singular:=c[2]; 
    elif c[1] eq "multiplicity" then mult:=c[2]; 
    elif c[1] eq "genus" then genus:=c[2]; 
    elif c[1] eq "source" then source:=c[2]*"@"*c[3]; 
    else error Sprintf("Tup command %o not recognised (not singular, multiplicity, genus, source)",c[1]);
    end if;
  end for;
  ch1:=[* x: x in ch1 | Type(x) eq TBox *];          // vertex boxes
  return singular,mult,genus,source,ch1;
end function;    


function HLineMainComponent(B,genus,mult,singular,source,childrenabove,childrenbelow)
  if IsEmpty(childrenabove) then 
    xpos1:=0;
    xpos2:=0;
  else
    xpos1:=childrenabove[1]`x1 + childrenabove[1]`banchors[1];
    xpos2:=childrenabove[#childrenabove]`x1 + Last(childrenabove[#childrenabove]`banchors);
  end if;
  lastchainxpos:=[t+AbsX(b)-AbsX(B): t in b`tanchors, b in childrenbelow];
  xpos1:=Min(lastchainxpos cat [xpos1]);
  xpos2:=Max(lastchainxpos cat [xpos2]);
  color:=singular select "l1sing" else "";
  smult:=Sprint(mult);
  if genus ne 0 then smult*:=Sprintf(" g%o",genus); end if;
  multlen:=1/4*(#smult+1);
  ReplaceString(~smult,"g","\\smash{g}");

  H:=HLine(xpos2-xpos1+2/3,smult,"main",color,"l1",source);
  return H,xpos1;
end function;


procedure BottomAnchorsFromChildren(~B,ch); 
  B`banchors:=[b`x1+a: a in b`banchors, b in ch];
end procedure;


procedure TopAnchorsFromChildren(~B,ch); 
  B`tanchors:=[b`x1+a: a in b`tanchors, b in ch];
end procedure;


procedure BoxChain(~G,~D,P: fromleft:=false, placetop:=false, nolinkup:=false, removechain:=true)
  // P = sequence of vertices v1->...->vn
  // places all v2..v(n-1) with their children in a box
  // fromleft: lonely children followed by chains left to right

  LP:=[Label(v): v in P];
  if placetop then Append(~LP,0); end if;
  
  B:=TBoxCreate();

  for i:=1 to #LP-1 do                                    // intermediate steps
    Y:=B`height;

    ch1:=i eq 1 select [**]  else D[[LP[i]]];               // lonely children of P[i]
    ch2:=LP[i+1] eq 0 select [**] else D[[LP[i],LP[i+1]]];  // chains P[i]->P[i+1]
    if i ne 1 then
      singular,mult,genus,source,ch1:=VertexData(D[[LP[i]]]);
    end if;
    if not removechain then
      ch1:=[* Duplicate(c): c in ch1 *];
      ch2:=[* Duplicate(c): c in ch2 *];
    end if;

    ch:=fromleft select (ch1 cat ch2) else (ch2 cat ch1);
   
    if #ch1 eq 0 
      then MaxCh1Height:=0;
      else MaxCh1Height:=Max([b`height: b in ch1]) + 1/2;
    end if;

    H:=Max([b`height: b in ch2] cat [MaxCh1Height]);
    for j:=1 to #ch2 do
      VerticalScaleTo(~ch2[j],H);
    end for;
    t:=B`tanchors;

    cmds:=IsEmpty(t) select [**] else [*<"A2",t>*]; // avoid links up
    for j:=1 to #ch do
      //cmdsj:=j eq 1 select [**] else [*<"R",B>*];        // right of the whole box
      cmdsj:=/*j eq 1 select [**] else*/ [*<"R",ch[k]>: k in [1..j-1]*];        // right of previous children
      Place(~B,ch[j],cmdsj cat cmds: ypos:=Y);
    end for;

    if i eq 1 then
      BottomAnchorsFromChildren(~B,ch); 
    else
      H,xpos1:=HLineMainComponent(B,genus,mult,singular,source,ch,t);
      Place(~B,H,[]: xpos:=xpos1, ypos:=Y, first);
    end if;

    B`tanchors:=ch2;
  end for;

  if nolinkup 
    then B`tanchors:=[];
    else TopAnchorsFromChildren(~B,B`tanchors); 
  end if;

  chlab:=placetop select [LP[1]] else [LP[1],LP[#P]];
  if not IsDefined(D,chlab) then D[chlab]:=[**]; end if;

  Append(~D[chlab],B);

  if not placetop then 
    G+:={P[1],P[#P]};
  end if;

  if removechain then
    G:=G-{P[i]: i in [2..#LP-1]};
  end if;
  
end procedure;


function LiftPathToG(G,P)       // lift a sequence of vertices from a subgraph to a graph
  L:=[Label(v): v in P];        // via their labels
  GP:=[];
  count:=0;
  for v in Vertices(G) do
    p:=Position(L,Label(v));
    if p eq 0 then continue; end if;
    GP[p]:=v;
    count+:=1;
  end for;
  assert count eq #P;
  return GP;
end function;


function LongestTail(G)
  maxlen:=0; 
  maxP:=[];
  for v in VertexSet(G) do
  for cv in Components(G-{v}) do
    c:=sub<G|Set(LiftPathToG(G,SetToSequence(cv)))>;
    N:=NumberOfVertices(c);
    if Diameter(c) ne N-1 then continue; end if;  // not a chain
    P:=LiftPathToG(G,DiameterPath(c));
    if #P le maxlen then continue; end if;
    fromv:=[w: w in P | w in Neighbours(v)];
    if (#fromv ne 1) or (fromv[1] ne P[1]) and (fromv[1] ne P[#P]) then continue; end if;
    if fromv[1] eq P[#P] then P:=Reverse(P); end if;    // path touching v at P[1] only
    maxlen:=N;
    maxP:=[v] cat P;
  end for;
  end for;
  return maxlen,maxP;
end function;


function LongestChain(G)
  maxlen:=0; 
  maxP:=[];
  for v1,v2 in VertexSet(G) do
  if Label(v1) ge Label(v2) then continue; end if;
  for cv in Components(G-{v1,v2}) do
    c:=sub<G|Set(LiftPathToG(G,SetToSequence(cv)))>;
    N:=NumberOfVertices(c);
    if Diameter(c) ne N-1 then continue; end if;        // not a chain
    P:=LiftPathToG(G,DiameterPath(c));
    if #P le maxlen then continue; end if;
    fromv1:=[w: w in P | w in Neighbours(v1)];
    fromv2:=[w: w in P | w in Neighbours(v2)];
    if (#fromv1 ne 1) or (#fromv2 ne 1) or {fromv1[1],fromv2[1]} ne {P[1],P[#P]} then continue; end if;
    if fromv1[1] eq P[#P] then P:=Reverse(P); end if;    // path touching v1 at P[1] only, v2 at P[#P]
    maxlen:=N;
    maxP:=[v1] cat P cat [v2];
  end for;
  end for;
  return maxlen,maxP;
end function;


function MaxBoxHeight(ch)
  return Max([0] cat [b`height + (IsEmpty(b`tanchors) select 1/2 else 0): b in ch | Type(b) eq TBox]);
end function;


procedure PlaceHLineMainComponent(~B,y,d,chb,chu: first:=true)
  chb:=[c: c in chb | Type(c) eq TBox];
  chu:=[c: c in chu | Type(c) eq TBox];
  singular,mult,genus,source:=VertexData(d);
  H,xpos1:=HLineMainComponent(B,genus,mult,singular,source,chu,chb);
  Place(~B,H,[]: xpos:=xpos1, ypos:=y, first:=first);
end procedure;


function Styles(s,L)
  style:="";
  for d in Reverse(L) do
    k,v:=Explode(d);
    if Regexp("[\[\{,=]"*k*"[\]\},=]",s*style) then
      style:=Sprintf("  %o/.style={%o},\n",k,v)*style;
    end if;
  end for;
  return Prune(Prune(style));
end function;


function GetStyle(Style,key: bool:=false, forcelist:=false)
  ok:=exists(val){d[2]: d in Style | d[1] eq key};
  error if not ok, Sprintf("GetStyle: key %o not found",key);
  if bool then forcelist:=["true","false"]; end if;
  error if forcelist cmpne false and val notin forcelist,
  Sprintf("GetStyle: value %o=%o is not in %o",key,val,DelSpaces(forcelist));
  return bool select (eval val) else val;  
end function;


function TikzScale(s)
  isone:=func<s|s in ["","1","1.0","1.00","1.000"]>;
  if isone(s) then return ""; end if;                         // no scaling
  ok,buf:=Match(s,"%s,%s");
  if not ok then return Sprintf("scale=%o",s); end if;
  xscale:=isone(buf[1]) select "" else "xscale="*buf[1];
  yscale:=isone(buf[2]) select "" else "yscale="*buf[2];
  return PrintSequence([u: u in [xscale,yscale] | u ne ""]: sep:=",");
end function;


intrinsic TeX(G::GrphUnd, D::Assoc, Style::SeqEnum[SeqEnum]) -> MonStgElt
{Returns a tikz picture drawing the dual graph}

  ZSeq:=Parent([2]);

  // Check for multiply assigned chains
  for k in Keys(D) do
    error if (#Set(k) eq 2) and IsDefined(D,Reverse(k)), Sprintf("D[%o] and D[%o] both assigned",k,Reverse(k)); 
  end for;

  // Now duplicate and reverse them [strings are not reversed, only integer chains]
  for k in Keys(D) do
  if (#k eq 2) and not IsDefined(D,Reverse(k)) then
    D[Reverse(k)]:=[* IsCoercible(ZSeq,d) select 
      Reverse(ZSeq!d) else d: d in D[k] *];
  end if;
  end for;

  for k in Keys(D) do
  for i:=1 to #D[k] do
    d:=D[k][i];
    if (#k eq 1) and (d cmpeq [-1]) then      // Singular point on a component
      D[k][i]:=RedBullet();
    elif (#k eq 1) and (Type(d) eq List) and (#d eq 2) then    // Special point on a component
      D[k][i]:=SpecialPoint(d[1],d[2]);                        // [* style, tex *]
    elif (#k eq 1) and (Type(d) eq List) and (#d eq 1) and d[1] cmpeq "bluenode" then  
      D[k][i]:=SelfP1Chain([],"blue","qmark");
    elif (#k eq 2) and (#Set(k) eq 1) /*and IsCoercible(ZSeq,d)*/ then  // non-singular chain of P1s
      D[k][i]:=Box(/*ZSeq!*/d,false: self);                // from a component to itself
    elif IsCoercible(ZSeq,d) then                      // chain of P1s between two
      chain:=ZSeq!d;                        
      color:="";
      if chain eq [-1] then chain:=[""]; color:="l2sing"; end if;   // singular chain
      D[k][i]:=Box(chain,#k eq 2: color:=color);         
    elif ExtendedType(d) eq SeqEnum[MonStgElt] then
      color:="";
      style:="";
      if Left(d[1],6) eq "style=" then 
        style:=d[1][[7..#d[1]]]; d:=d[[2..#d]];
      end if;
      if d eq ["-1"] then d:=[""]; color:="l2sing"; end if; // singular chain
      D[k][i]:=Box(d,#k eq 2: color:=color, style:=style); 
    elif Type(d) eq Tup then 
      ;       // processed separately
    else
      error Sprintf("D[%o] = %o not understood",DelSpaces(k),DelSpaces(D[k][i]));
    end if;
  end for;
  end for;

  for k in Keys(D) do                   // Move self-chains into node data
  if (#k eq 2) and (#Set(k) eq 1) then
    D[[k[1]]] cat:= D[k];
    Remove(~D,k);
  end if;
  end for;
  
  repeat
    changed:=false;

    repeat            // find and box all the tails v1(->...->vn)
      vprint Newton,1: "Finding longest tail";
      maxlen,maxP:=LongestTail(G);
      if maxlen eq 0 then break; end if;
      changed:=true;
      vprint Newton,1: "Boxing a chain",maxP;
      BoxChain(~G,~D,maxP: fromleft, placetop);
      vprint Newton,1: "Done";
    until false;
    
    repeat            // find and box all the chains v1(->...->)vn
      maxlen,maxP:=LongestChain(G);
      if maxlen eq 0 then break; end if;
      changed:=true;
      BoxChain(~G,~D,maxP: fromleft, removechain:=false);
      BoxChain(~G,~D,Reverse(maxP): fromleft);
    until false;
  
  until not changed;
  
  V:=VertexSet(G);
  E:=EdgeSet(G);
  dd:=func<s|Type(s) eq RngIntElt select D[[Label(V.s)]] else D[[Label(V.i): i in s]]>;
      
  if (#V eq 4) and (#E eq 6) then    // K4 
    H1:=0;                                               // Height to place box #1 = 0 
    ch1:=dd([1,2]) cat dd(1); 
    H2:=MaxBoxHeight(ch1);                               // #2
    ch2A:=dd([1,3]);
    ch2B:=dd(2) cat dd([2,3]);
    H3:=Max(MaxBoxHeight(ch2A),H2+MaxBoxHeight(ch2B));   // #3
    ch3A:=dd([1,4]);
    ch3B:=dd([2,4]);
    ch3C:=dd(3) cat dd([3,4]);
    H4:=Max([MaxBoxHeight(ch3A),H2+MaxBoxHeight(ch3B),H3+MaxBoxHeight(ch3C)]); // #4

    VerticalScaleTo(~ch2A,H3);                            // Rescale all children accordingly
    VerticalScaleTo(~ch2B,H3-H2);
    VerticalScaleTo(~ch3A,H4);
    VerticalScaleTo(~ch3B,H4-H2);
    VerticalScaleTo(~ch3C,H4-H3);

    B:=TBoxCreate();
    PlaceBoxesAtHeight(~B,H1,ch3A cat ch2A cat ch1,[]);   //! could avoid stuff below #1 if necessary
    PlaceBoxesAtHeight(~B,H2,ch2B,ch1);
    PlaceBoxesAtHeight(~B,H3,ch3C,ch2A cat ch2B);
    PlaceHLineMainComponent(~B,H1,dd(1),[],ch3A cat ch2A cat ch1);
    PlaceHLineMainComponent(~B,H3,dd(3),ch2A cat ch2B,ch3C);  //!
    PlaceBoxesAtHeight(~B,H2,ch3B,ch1);
    PlaceHLineMainComponent(~B,H2,dd(2),ch1,ch2B cat ch3B);
    PlaceHLineMainComponent(~B,H4,dd(4),ch3A cat ch3B cat ch3C,[]);

  elif #V eq 1 then                // One remaining vertex

    B:=TBoxCreate();
    ch:=D[[Label(V.1)]];
    singular,mult,genus,source,ch:=VertexData(ch);
    for j:=1 to #ch do
      cmdsj:=j eq 1 select [**] else [*<"R",ch[j-1]>*];        // right of previous child
      Place(~B,ch[j],cmdsj);
    end for;
    H,xpos1:=HLineMainComponent(B,genus,mult,singular,source,ch,[]);
    Place(~B,H,[]: xpos:=xpos1, ypos:=0, first);
  else 
    return Sprintf("\\footnotesize Graph possibly unconnected or not planar --- cannot draw (V=%o E=%o)\\\\%o",#V,#E,DelCRs(Sprint(G,"Maximal")));
  end if;  
  
  out:=TeX(B);
  sty:=Styles(out,Style);
  if #sty ne 0 then sty:=",\n"*sty; end if;
  return Sprintf("\\begin{tikzpicture}[%o%o]\n%o\\end{tikzpicture}",
    TikzScale(GetStyle(Style,"ModelScale")),sty,out);

end intrinsic;


//////////////////////////////
// 3. Polytope and printing //
//////////////////////////////


function Select(L,x) 
  return [i: i in [1..#L] | L[i] eq x];
end function;


LastCoordinate:=func<vlist|[v.3: v in vlist]>;
FiniteFaces:=func<F,inf|[f: f in F | inf notin LastCoordinate(Vertices(f))]>;


function Down(P)  // Project a 3D polytope down to 2D
  assert Dimension(Ambient(P)) eq 3;
  return Polytope(SetToSequence({Prune(Eltseq(u)): u in Vertices(P)}));
end function;

function DeterminantOneFractions(f1,f2)       
  // Fill in [f1,f2] with internal fractions of small denominator a la Hirzebruch-Jung
  // f1<f2 and consecutive dets = -1 for historical reasons
  assert f1 lt f2;
  N:=Numerator;
  D:=Denominator;
  // integers between f1 and f2
  L:={f1,f2} join {n: n in [Ceiling(f1)..Floor(f2)]};
  // possible small fractions close to f1 and f2
  for d:=1 to D(f1)-1 do
    n:=Ceiling(d*f1);
    if n/d lt f2 then Include(~L,n/d); end if;
  end for;
  for d:=1 to D(f2)-1 do
    n:=Floor(d*f2);
    if n/d gt f1 then Include(~L,n/d); end if;
  end for;
  // Sort and prune unnecessary fractions
  L:=SortSet(L);
  repeat
    changed:=false;
    for i:=#L-1 to 2 by -1 do
    if D(L[i]) gt Max(D(L[i-1]),D(L[i+1])) then
      Remove(~L,i);
      changed:=true;    
    end if;
    end for;
  until not changed;
  // Check that everything is tipsy-dosy
  for i:=1 to #L-1 do
    assert N(L[i])*D(L[i+1])-N(L[i+1])*D(L[i]) eq -1;
  end for;
  return L;
end function;


// Shift a list of vectors (monomial exponents) into the positive cone, minimally

procedure ShiftMonomialsToZero(~Mvecs)
  shift:=[Min([v[j]: v in Mvecs]): j in [1..3]];
  for i in [1..#Mvecs], j in [1..3] do 
    Mvecs[i][j]-:=shift[j]; 
  end for;
end procedure;

// Transform an equation given by a list of monomial exponents and coefficients
// via a matrix Minv. Then reduce coefficients using red and Evaluate at the
// variables vars. Returns evaluated reduced equation, non-evaluated, non-reduced

function TransformPolynomial(f,Minv)
  Rk:=Parent(f);
  k:=BaseRing(Rk);
  cfs,mons:=CoefficientsAndMonomials(f);
  monvectors:=[Exponents(m): m in mons];
  Mvecs:=[Eltseq(Minv*Transpose(Matrix(Z,[m]))): m in monvectors];
  ShiftMonomialsToZero(~Mvecs);
  ftrans:=&+[cfs[i]*Monomial(Rk,Mvecs[i]): i in [1..#cfs]];
  return ftrans;
end function;


function TransformAndReduce(monvectors,cfs,red,Minv,vars)
  K:=Parent(cfs[1]);
  k:=Parent(red(1));
  RK:=PR(K,3);
  Rk:=PR(k,3);
  Mvecs:=[Eltseq(Minv*m): m in monvectors];
  ShiftMonomialsToZero(~Mvecs);
  ftrans:=&+[cfs[i]*Monomial(RK,Mvecs[i]): i in [1..#cfs]];
  ftransred:=&+[red(cfs[i])*Monomial(Rk,Mvecs[i]): i in [1..#cfs]];
  fred:=Evaluate(ftransred,vars);
  return fred,ftransred,ftrans;
end function;


function FactorisationMPolSplittingField(f)
   
  fct:=Factorisation(f);        // Reducible over the ground field - go in recursively

  if (#fct gt 1) or (#fct eq 1 and fct[1][2] gt 1) then
    out:=[];
    for d in fct do
      fmp:=FactorisationMPolSplittingField(d[1]);
      for j:=1 to #fmp do fmp[j][2]*:=d[2]; end for;
      out cat:= fmp;
    end for;
    return out;
  end if;

  // Use curve machinery to factor over the splitting field
  k:=BaseRing(Parent(f));
  //assert IsFinite(k);
  kk:=FieldOfGeometricIrreducibility(Curve(AffineSpace(k,2),f));
  return kk eq k select [<f,1>] else Factorisation(PR(kk,2)!f);
end function; 

procedure AddComponent(~comps,name,i,d1,d2,sing,M)
  Append(~comps,<name,i,d1,d2,sing,[**],[**],M>);     
end procedure; 

procedure AddChain(~comp,toC,mults,sing,src,~S)
  vprintf Newton,1: "%o->%o %o%o (src %o)\n",comp[1],toC,DelSpaces(mults),
    sing select " sing" else "",DelSpaces(src);

  if not IsEmpty(mults) then    // Style can overwrite (non-empty) multiplicities with "Ch%i"
    if exists(i){i: i in [1..#S] | S[i,1] eq "_chaincount"} then   // Keep chain count
      chaincount:=(eval S[i,2])+1; S[i,2]:=Sprint(chaincount);
    else
      chaincount:=1; Append(~S,["_chaincount","1"]);
    end if;
    mults:=[Sprint(x): x in mults];
    if GetStyle(S,"ChainNumbers": bool) then    // Print chain numbers if requested
      mults[#mults]:=Sprintf("\\llap{\\scalebox{0.6}{\\color{blue}%o}\\,}%o",chaincount,mults[#mults]);
    end if;
    if exists(newmults){d[2]: d in S | d[1] eq "Ch"*Sprint(chaincount)} then
      mults:=Split(newmults," ");
    end if;    
  end if;
  
  Append(~comp[6],<toC,mults,sing,src>);
end procedure;

procedure AddSingularPoint(~comp,src,name)
  Append(~comp[7],<src,name>);
end procedure;

// Orient a sequence of vertices for drawing

function Orient(S)  
  if #S le 2 then return S; end if;
  C:=ComplexField();
  avg:=&+S / #S;
  args:=[Arg(C!Eltseq(v-avg)): v in S];
  Sort(~args,~perm);
  return S[Eltseq(perm)];
end function;

function angle(v)
  arg:=Round(180*Arg(ComplexField()!Eltseq(v))/Pi(RealField()));
  return arg lt 0 select arg+360 else arg;
end function;

und:=func<x|#s le 1 select s else Sprintf("{%o}",s) where s is Sprint(x)>;

function FaceNameTeX(s,Style: dollar:="$")
  if exists(d){d[2]: d in Style | d[1] eq s} then
    L:=Split(d,"@": IncludeEmpty);
    return dollar*L[1]*dollar,L[[2..#L]];
  end if;
  ok,buf:=Match(s,"FL%i");
  if ok then return Sprintf(dollar*"F^L_%o"*dollar,und(buf[1])),[]; end if;
  ok,buf:=Match(s,"F%i");
  if ok then return Sprintf(dollar*"F_%o"*dollar,und(buf[1])),[]; end if;
  ok,buf:=Match(s,"F%i%s");
  if ok then return Sprintf(dollar*"F_%o^%o"*dollar,und(buf[1]),und(buf[2])),[]; end if;
  error "FaceNameTeX: unrecognized face name "*s;
end function;


function printpower(v,n)
  if n eq 0 then return ""; end if;
  if n eq 1 then return v; end if;
  n:=Sprint(n);
  return #n eq 1 select Sprintf("%o^%o",v,n) else Sprintf("%o^{%o}",v,n);
end function;


function PrintMatrixAsMonomialTransformation(m,v)
  L:=[];
  for r in Rows(Transpose(m)) do
    a,b,c:=Explode(Eltseq(r));
    s:=Sprintf("%o%o%o",printpower(v[1],a),printpower(v[2],b),printpower(v[3],c));
    if s eq "" then s:="1"; end if;
    Append(~L,s);
  end for;
  return L;
end function;


intrinsic DVRData(K::FldRat, p::RngIntElt) -> .
{DVRData K,k,v,red,lift,pi for K=Q, p=prime}
  require IsPrime(p): "p is not a prime number";
  k:=GF(p);
  v:=func<x|Valuation(x,p)>;
  red:=func<x|k!x>;
  lift:=func<x|K!Z!x>;
  pi:=p;
  return K,k,v,red,lift,pi;
end intrinsic;


intrinsic DVRData(Z::RngInt, p::RngIntElt) -> .
{DVRData K,k,v,red,lift,pi for O=Z, p=prime}
  return DVRData(Q,p);
end intrinsic;


intrinsic DVRData(K::FldNum, p::RngOrdIdl) -> .
{DVRData K,k,v,red,lift,pi for K=number field, p=prime ideal}
  require IsPrime(p): "p is not a prime ideal";
  v:=func<x|Valuation(x,p)>;
  OK:=Integers(K);
  k,m:=ResidueClassField(p);
  red:=func<x|m(OK!x)>;
  lift:=func<x|K!(x@@m)>;
  pi:=UniformizingElement(p);
  return K,k,v,red,lift,pi;
end intrinsic;


intrinsic DVRData(K::FldNum, p::PlcNumElt) -> .
{DVRData K,k,v,red,lift,pi for K=number field, p=place}
  return DVRData(K,Ideal(p));
end intrinsic;


intrinsic DVRData(O::RngOrd, p::RngOrdIdl) -> .
{DVRData K,k,v,red,lift,pi for O=integers in a number field, p=prime ideal}
  return DVRData(FieldOfFractions(O),p);
end intrinsic;


intrinsic DVRData(K::FldPad) -> .
{DVRData K,k,v,red,lift,pi for K=p-adic field}
  v:=func<x|Valuation(x)>;
  OK:=Integers(K);   
  k,m:=ResidueClassField(OK);
  pi:=UniformizingElement(OK);  
  red:=func<x|m(OK!x)>;
  lift:=func<x|K!(x@@m)>;
  return K,k,v,red,lift,pi;
end intrinsic;


intrinsic DVRData(O::RngPad) -> .
{DVRData K,k,v,red,lift,pi for O=integers in a p-adic field}
  return DVRData(FieldOfFractions(O));
end intrinsic;


intrinsic DVRData(K::FldFunRat, p::FldFunRatUElt) -> .
{DVRData K,k,v,red,lift,pi for rational function field in one variable}
  assert Rank(K) eq 1;
  OK:=Integers(K);
  pi:=OK!p;  
  v:=func<x|Valuation(x,pi)>;
  k,m:=ResidueClassField(OK,ideal<OK|pi>);
  red:=func<x|m(OK!x)>;
  lift:=func<x|K!(x@@m)>;
  return K,k,v,red,lift,pi;
end intrinsic;


intrinsic DVRData(K::FldFunRat) -> .
{DVRData K,k,v,red,lift,pi for rational function field in one variable}
  return DVRData(K,Integers(K)!(K.1));
end intrinsic;


intrinsic DVRData(K::RngSerLaur) -> .
{DVRData K,k,v,red,lift,pi for Laurent series ring in one variable}
  OK:=Integers(K);
  pi:=OK!(K.1);  
  v:=func<x|Valuation(x)>;
  k,m:=ResidueClassField(OK);
  red:=func<x|m(OK!x)>;
  lift:=func<x|K!(x@@m)>;
  return K,k,v,red,lift,pi;
end intrinsic;


intrinsic DVRData(O::RngSerPow) -> .
{DVRData K,k,v,red,lift,pi for power series ring in one variable}
  return DVRData(FieldOfFractions(O));
end intrinsic;


intrinsic DefaultSettings(DeltaScale,ModelScale) -> Assoc
{Default settings for Model(...)}
  S:=AssociativeArray();
  S["ShowEquation"]:=false;
  S["ShowCharts"]:=false;
  S["DeltaScale"]:=DeltaScale;
  S["ModelScale"]:=ModelScale;
  return S;
end intrinsic;


intrinsic Model(f::RngMPolElt,P: Style:=AssociativeArray(), DVR:="default") -> MonStgElt, Assoc
{TeX model for a curve C given by f=0 (main function)}

  DefaultStyle:=[
    // Formats
    ["OutputFormat","%o\\hfill%o\\hfill\\hfill\n%o\\hfill%o"],
    ["DeltaFormat","\\pbox[c]{20cm}{\n%o\n}"],
    ["ModelFormat","\\pbox[c]{20cm}{\n%o\n}"],
    ["ChartsFormat","\n\\par\\bigskip{\\small $\\begin{array}{lllll}\n%o\\end{array}\n$}"],
    ["EquationFormat","{\\small $f=%o$}\n\n\\vskip 7pt\n"],
    // Main parameters
    ["SmallValuations","true"],
    ["LargeValuations","true"],
    ["ContractFaces","true"],
    ["RemoveFaces","true"],
    ["BakerRegulars","true"],
    ["ShowCharts","false"],
    ["FaceNames","true"],
    ["ChainNumbers","false"],
    ["Warnings","true"],
    ["ShowEquation","false"],    // true, false or an integer (max number of terms)
    ["DeltaScale","0.8,0.6"],    // e.g. "1" (no scaling), "0.7" or "0.8,0.6" (vert,horiz)
    ["ModelScale","0.8,0.7"],    //      same
    // Delta_v tikz style 
    ["sml","scale=0.55"], 
    ["lrg","scale=0.9,inner sep=0.2em"], 
    ["fname","blue,scale=0.55"], 
    ["lname","scale=0.55,sloped,red,above=-0.07em,near end"],
    ["lin","-,shorten <=-0.07em,shorten >=-0.07em"],
    ["rem","black!20,thin"],
    // Regular model tikz style 
    ["l1","shorten >=-1.3em,shorten <=-0.5em,thick"],
    ["l2","shorten >=-0.3em,shorten <=-0.3em"],
    ["l1sing","very thick,red!70"],
    ["l2sing","l2,very thick,red!70"],
    ["lfnt","font=\\tiny"],
    ["rlab","xshift=15,yshift=6"],     // label to the right of the component (rather than below)
    ["qmark","coordinate,label={[scale=0.5]?}"],   // question mark for self-node of unknown depth
    ["zigz","l2,snake=zigzag,segment length=2,segment amplitude=1,blue!70!black,shorten >=0,shorten <=0"],
    ["qmleft","lfnt,above left=0.15em and -0.06em,scale=0.7"],
    ["leftl","left=-3pt,lfnt"],
    ["rightl","right=-3pt,lfnt"],
    ["mainl","scale=0.8,above left=-0.17em and -1.5em"],
    ["abovel","above=-2.5pt,lfnt"],
    ["belowl","below=-2.5pt,lfnt"],
    ["aboveleftl","above left=-4.5pt,lfnt"],
    ["aboverightl","above right=-4.5pt,lfnt"],
    ["belowleftl","below left=-4.5pt,lfnt"],
    ["facel","scale=0.5,blue,below right=-0.5pt and 6pt"],
    ["redbull","red,label={[red,scale=0.6,above=-0.17]#1}"],
    ["bluedash","blue!80!white,dashed"],
    ["belowrightl","below right=-4.5pt,lfnt"]
  ];  

  OldStyle:=Style;
  NewStyle:=DefaultStyle;      // Default + Style associative array -> Style list
  for k in Keys(Style) do
    ok:=exists(i){i: i in [1..#NewStyle] | NewStyle[i][1] eq k};  
    if not ok and k[1] in ["F","V","L","C"] then     // Face/vertex/(line)/chain styles
      Append(~NewStyle,[k,Style[k]]); continue;
    end if;
    error if not ok, Sprintf("Style: key %o is not recognized",k);
    NewStyle[i][2]:=Sprint(Style[k]);
  end for;
  Style:=NewStyle;

  // Process options

  ShowEquation:=GetStyle(Style,"ShowEquation");
  ok,val:=Match(ShowEquation,"%i");
  if ok then                        // ShowEquation=n: show equation if f has <=n terms
    ShowEquation:=#Monomials(Numerator(f)) le val[1]; 
  else
    ShowEquation:=GetStyle(Style,"ShowEquation": bool);
  end if;

  SmallValuations := GetStyle(Style,"SmallValuations": bool);
  LargeValuations := GetStyle(Style,"LargeValuations": bool);
  ContractFaces   := GetStyle(Style,"ContractFaces": bool);
  RemoveFaces     := GetStyle(Style,"RemoveFaces": bool);
  BakerRegulars   := GetStyle(Style,"BakerRegulars": bool);
  ShowCharts      := GetStyle(Style,"ShowCharts": bool);
  FaceNames       := GetStyle(Style,"FaceNames": bool);
  Warnings        := GetStyle(Style,"Warnings": bool);
  
  K:=BaseRing(Parent(f));
  if DVR cmpne "default" then
    assert #DVR eq 6;
    K,k,vK,red,lift,pi:=Explode(DVR);
  elif Type(K) in [FldPad,RngPad,RngSerLaur,RngSerPow] then 
    error if P cmpne "default", "P should not be specified for K:"*Sprint(Type(K));
    K,k,vK,red,lift,pi:=DVRData(K); 
  else 
    error if P cmpeq "default", "P must be specified for K:"*Sprint(Type(K));
    K,k,vK,red,lift,pi:=DVRData(K,P);
  end if;
  
  R<x,y>:=PR(K,2);
  forg:=f;

  procedure Warning(s)
    if Warnings then "Warning: "*s; end if;
  end procedure;

  if IsFinite(k) then
    RootsFunction:=RootsInSplittingField;
    PointsFunction:=PointsOverSplittingField;
    FactorizationFunction:=FactorisationMPolSplittingField;
  else
    Warning("Residue field is not finite. Do not have an algorithm to factor/find roots\n"*
            "for schemes in extensions. May miss components in the regular model.\n"*
            "Use at your own risk");
    RootsFunction:=Roots;
    FactorizationFunction:=Factorisation;
    PointsFunction:=Points;
  end if;
  
  if Rank(Parent(f)) eq 3 
    then f:=Numerator(Evaluate(f,[x,y,pi]));
    else f:=Numerator(Evaluate(f,[x,y]));
  end if;

  
  cfs,monsP:=CoefficientsAndMonomials(f);
  error if #monsP eq 1, "One monomial equation: empty Newton polygon"; 
  monsP:=[Exponents(m): m in monsP];
  
  mons:=[];
  for i:=1 to #cfs do
    val:=vK(cfs[i]);
    cfs[i]:=cfs[i]/pi^val;
    mons[i]:=monsP[i] cat [val];
  end for;
  monvectors:=[Transpose(Matrix(Z,[m])): m in mons];   // monomials as column vectors
  Delta:=Polytope([Prune(m): m in mons]);
  AllP:=SortSet(Points(Delta));
  DeltaCenter:=&+AllP/#AllP;
  InnP:=SortSet(InteriorPoints(Delta));
  
  inf:=1+Max([Last(m): m in mons]);               // 0,1,2-dim faces of Delta tilde
  DeltaTop:=[[m[1],m[2],inf]: m in mons];         
  P:=Polytope(mons cat DeltaTop);
  F:=Faces(P);
  F:=[FiniteFaces(FList,inf): FList in F];
  Vtilde,Ltilde,Ftilde,_:=Explode(F);
  if assigned FaceReorder then Ftilde:=Ftilde[FaceReorder]; end if;
  if assigned LineReorder then Ltilde:=Ltilde[LineReorder]; end if;
  F2H:=[Last(Inequalities(F)): F in Ftilde];
  
  Vs:=[Vertices(Down(F))[1]: F in Vtilde];
  Ls:=[Down(F): F in Ltilde];
  Fs:=[Down(F): F in Ftilde];
  
  function vF(j,x,y)            // z-coordinate function on jth 2-face evaluated at [x,y]
    Heqn,Hv:=Explode(F2H[j]);
    Hx,Hy,Hz:=Explode(Eltseq(Heqn));
    return (Hv - Hx*x - Hy*y)/Hz;
  end function;
  
  function FacePositions(S)
    if IsEmpty(S) then return []; end if;
    if Dimension(S[1]) eq 0 then return [Position(Vs,Vertices(x)[1]): x in S]; end if;
    if Dimension(S[1]) eq 1 then return [Position(Ls,x): x in S]; end if;
    if Dimension(S[1]) eq 2 then return [Position(Fs,x): x in S]; end if;
    error "FacePositions: unrecognized dimension";
  end function;
  
  vP:=AssociativeArray();                        // height vP: Delta -> Q
  for j in [1..#Ftilde], P in Points(Fs[j]) do   // defined on integer points
    xP,yP:=Explode(Eltseq(P));
    if not IsDefined(vP,P) then
      vP[P]:=vF(j,xP,yP);
    end if;
  end for;
  
  AllFP:=[SortSet(Points(F)): F in Fs];
  AllLP:=[SortSet(Points(L)): L in Ls];
  
  InnFP:=[SortSet(InteriorPoints(F)): F in Fs];   // AllLP = L meet Z^2, AllFP = F meet Z^2
  InnLP:=[SortSet(InteriorPoints(L)): L in Ls];   // and inner versions InnLP, InnFP
  
  FIsRemovable:=[false: F in Fs];    // Repeatedly remove removable faces
  newleft:=#FIsRemovable;
  if RemoveFaces then
  repeat
    oldleft:=newleft;
    FIsRemovable:=[IsEmpty(InnFP[j]) and
      #(Set(AllFP[j]) meet &join{Set(AllFP[i]): i in [1..#Fs] | 
        not FIsRemovable[i] and (i ne j)}) le 2: j in [1..#Fs]];
    newleft:=Count(FIsRemovable,false);
  until oldleft eq newleft;
  end if;
  
  deltaL:=[Denominator(VectorContent([vP[P]: P in AllLP[j]])): j in [1..#Ls]];  // denominators
  deltaF:=[Denominator(VectorContent([vP[P]: P in AllFP[j]])): j in [1..#Fs]];
  
  degreeL:=[Z| (#Points(Ls[i])-1)/deltaL[i] : i in [1..#Ls]];
  
  FL:=[FacePositions(Faces(F,1)): F in Fs];             // list of 1-faces on a given 2-face
  FV:=[FacePositions(Faces(F,0)): F in Fs];             // list of 0-faces on a given 2-face
  FC:=[(&+Vs[FV[j]])/#(FV[j]): j in [1..#Fs]];          // centers of 2-faces
  
  LF:=[[j: j in [1..#Fs] | (i in FL[j])                 // list of 2-faces bounding a given 1-face
    and not FIsRemovable[j]]: i in [1..#Ls]];           //   [only the non-removable ones]
  
  LSquarefree:=[];                    // 1-faces with squarefree reduction of f_L
  for i:=1 to #Ls do
    lcf:=[];
    for P in AllLP[i] do
      if vP[P] notin Z then continue; end if;
      j:=Position(monsP,Eltseq(P));
      Append(~lcf,(j eq 0) or (mons[j][3] ne vP[P]) select 0 else red(cfs[j]));
    end for;
    fL:=Polynomial(lcf);
    Append(~LSquarefree,Degree(GCD(fL,Derivative(fL))) eq 0);
  end for;  
  
  FIsContractible:=[false: F in Fs];    
    // Contractible faces: triangles, two 1-faces with delta=1 and length 1, 
    // meeting at an interior point, and one outer face with delta=delta_F
    // or one triangular In-like face [Riemann-Hurwitz]
    // + special case for genus 1?
  //! Should check face is non-singular in the In case
  cc:=[];
  if #InnP eq 0 then 
    Warning("Contractible faces in genus 0 are not implemented");
  else
  for i:=1 to #Fs do 
    if FIsRemovable[i] then continue; end if;                         // want not removable 
    if exists{P: P in InnFP[i] | vP[P] in Z} then continue; end if;   //  and not of positive genus
    l:=FL[i];
    if #l ne 3 then continue; end if;
    if not forall{j: j in l | LSquarefree[j]} then continue; end if;
    chains:=&cat[ [l[j]: count in [1..degreeL[l[j]]]] : j in [1..3] | (deltaL[l[j]] lt deltaF[i]) or (#LF[l[j]] gt 1)];
    error if #chains lt 2 and RemoveFaces, "Did not expect a removable face here";
    if #chains gt 2 then continue; end if;
    if not forall{j: j in [1..3] | (#LF[l[j]] eq 1) or LSquarefree[l[j]]} 
      then continue; end if;                                          // and non-singular inner face
    vprintf Newton,1: "Contractible face: #%o %o (%o)\n",i,PrintSequence(Vertices(Fs[i]): d),DelSpaces(chains);
    if {i: i in [1..#Fs] | not FIsRemovable[i] and not FIsContractible[i]} eq {i} 
      then continue; end if; // last remaining not contractible face  
    FIsContractible[i]:=true;
    Append(~cc,chains);
  end for;
  end if;
  
//  if forall{i: i in [1..#Fs] | FIsRemovable[i] or FIsContractible[i]} then
//    [&+[&+Eltseq(v): v in Vs[FV[i]]]: i in [1..#Fs]];
//    assert exists(i){i: i in [1..#Fs] | not FIsRemovable[i]};
//    FIsContractible[i]:=false;
//    error "All faces are contractible!";  
//  end if;
    
  if RemoveFaces and ContractFaces then
    assert forall{i: i in Select(FIsContractible,true) | #(Set(AllFP[i]) meet Set(InnP)) eq 1};
    assert forall{i: i in Select(FIsRemovable,true) | #(Set(AllFP[i]) meet Set(InnP)) eq 0};
  end if;
  
  // merge contractible chains
  while exists(d){[i,j]: i,j in [1..#cc] | (i lt j) and #(Set(cc[i]) meet Set(cc[j])) eq 1} do
    c1:=cc[d[1]];
    c2:=cc[d[2]];
    if c1[1] eq c2[1]         then c:=Prune(Reverse(c1)) cat c2; 
    elif c1[1] eq Last(c2)    then c:=Prune(c2) cat c1; 
    elif Last(c1) eq c2[1]    then c:=Prune(c1) cat c2; 
    elif Last(c1) eq Last(c2) then c:=Prune(c1) cat Reverse(c2); 
      else error "Could not merge chains"; 
    end if;
    vprintf Newton,0: "Merging chains %o %o -> %o\n",c1,c2,c;
    cc[d[1]]:=c;
    Remove(~cc,d[2]);
  end while;
  ContractibleChains:=cc;
    
  if not ContractFaces then
    FIsContractible:=[false: F in Fs];    
    ContractibleChains:=[];
  end if;
    
  M:=[[]: l in Ls];                             // chart transformation matrices
  MInv:=[[]: l in Ls];                          // and their inverses
  
  LSlopes:=[];          // intermediate det 1 fractions from sF1 to sF2, inclusive
  ki:=[];               // list of [k_0,k_1,...,k_{r+1}] for every 1-face L
  
  LStarPoints:=[];      // list of [P0,P1,P2] for every 1-face L
  for i:=1 to #Ls do    // P1-P0 on L, counterclockwise on F1, 
    if IsEmpty(LF[i]) then continue; end if;  // removable
    F1j:=LF[i][1];      // P1-P0, P2-P1 form a determinant one basis 
    cF1:=FC[F1j];       // (second one points into F1)
    F2j:=#LF[i] eq 2 select LF[i][2] else 0;    // 0 if outer
    lp:=AllLP[i];
    P0:=lp[1]; 
    P1:=lp[2]; 
    v:=P1-P0;
    vx:=v.1;
    vy:=v.2;
    if -vy*(cF1.1-P0.1)+vx*(cF1.2-P0.2) lt 0 then  // want P1-P0 counterclockwise on F1
      P0:=lp[#lp];
      P1:=lp[#lp-1];
      v:=P1-P0;
    end if;
    vx:=v.1;
    vy:=v.2;
    vz:=vF(F1j,P1.1,P1.2)-vF(F1j,P0.1,P0.2);    // nu vector = (v.1,v.2,vz)
    g,a,b:=XGCD(vx,vy);
    assert -vy*(cF1.1-P0.1)+vx*(cF1.2-P0.2) gt 0;
    assert g eq 1;
    P2:=P0 + Parent(v)![-b,a];
    wx:=P2.1-P0.1;
    wy:=P2.2-P0.2;
    assert Determinant(Matrix(Q,[Eltseq(P1-P0),Eltseq(P2-P0)])) eq 1;
    LStarPoints[i]:=[P0,P1,P2];
  
    delta:=deltaL[i];
    
    sF1:=delta*(vF(F1j,P2.1,P2.2)-vF(F1j,P0.1,P0.2));    // s_F1
    outer:=F2j eq 0;
    inner:=not outer;
    sF2:=outer select Floor(sF1-1) else                  // s_F2 in the outer case
         delta*(vF(F2j,P2.1,P2.2)-vF(F2j,P0.1,P0.2));    // s_F2 in the inner case
    LSlopes[i]:=[-x: x in DeterminantOneFractions(-sF1,-sF2)];   
  
    deltavzinv:= delta eq 1 select 0 else (Integers(delta)!Z!(delta*vz))^(-1);
    
    ki[i]:=[Z|-Numerator(x)*deltavzinv: x in LSlopes[i]];
  
    for j:=1 to #LSlopes[i]-1 do
      di:=Denominator(LSlopes[i][j]);
      dip:=Denominator(LSlopes[i][j+1]);
      ni:=Numerator(LSlopes[i][j]);
      nip:=Numerator(LSlopes[i][j+1]);
      kii:=ki[i][j];
      kip:=ki[i][j+1];
      eps:=outer and (j+1 eq #LSlopes[i]) select 1 else 0;      // Modify in the outer case for the last j
      if eps eq 1 then assert di eq dip; end if;
  
      M[i][j]:=Matrix(Z, [ 
         [delta*vx,  di*wx+kii*vx,    -dip*wx-kip*vx  +   eps*(di*wx+kii*vx)  ],
         [delta*vy,  di*wy+kii*vy,    -dip*wy-kip*vy  +   eps*(di*wy+kii*vy) ],
         [delta*vz,  ni/delta+kii*vz, -nip/delta-kip*vz + eps*(ni/delta+kii*vz) ] ]);
      assert Determinant(M[i][j]) eq 1;
      MInv[i][j]:=M[i][j]^(-1);
    end for;
  end for;
  
  LPtsRec:=recformat<L, root, mult, sing, name>; 
  FIsSingular:=[];
  
  // now construct the dual graph,
  // extending k repeatedly if necessary
  
  newk:=k;
  for trial:=1 to 40 do   
  
    if newk ne k then                                   // extend k
      vprint Newton,1: "Extending",k,"to",newk;
      k:=newk; oldred:=red;   red:=func<x|k!oldred(x)>; 
    end if;  
  
    comps:=[];                   // [<name,Findex,redeqn,mult,sing,chains,singpts,M=transmatrix>,...]
    Lpts:=[];                    // [<Lindex,root,mult,sing,name>,...]
    singptcount:=0;
  
    Rk1<t>:=PR(k);
    Rk2<X2,Y2>:=PR(k,2);
    A:=AffineSpace(Rk2);
  
    for i:=1 to #Ls do
      if IsEmpty(LF[i]) then continue; end if;  // removable
      tr:=MInv[i][1];
      ffred:=TransformAndReduce(monvectors,cfs,red,tr,[X2,Y2,0]);
      m:=Min([Degree(m,Y2): m in Monomials(ffred)]);
      flred:=Evaluate(ffred div Y2^m,[t,0]);
      rts:=[r: r in RootsFunction(flred) | r[1] ne 0];
      if not IsEmpty(rts) then 
        newk:=Parent(rts[1][1]); 
        if newk ne k then continue trial; end if;
      end if;
      for r in rts do 
        sing:=r[2] gt 1;
        if sing and (Evaluate(Derivative(ffred,2),[r[1],0]) ne 0) and     // Baker regular
           BakerRegulars and (#LF[i] eq 1) and (deltaF[LF[i][1]] eq 1)    // on outer L with delta_F=1
        then     
          vprint Newton,1: "Baker regular",i;
          sing:=false;
        end if;
        if sing then 
          name:=CodeToString(97+singptcount);
          singptcount+:=1;
        else 
          name:="";
        end if;
        Append(~Lpts,rec<LPtsRec|L:=i,root:=r[1],mult:=r[2],sing:=sing,name:=name>);   
      end for;    
      vprint Newton,1: Sprintf("L%o %o",i,PrintSequence([DelSpaces(r): r in rts]));
    end for;
  
    Lvanishing:=[[]: P in Lpts]; // [<compindex,orderofvanishing>]
    
    for i:=1 to #Fs do
      FIsSingular[i]:=false;
      if FIsRemovable[i] /*or FIsContractible[i]*/ then continue; end if;
      F:=Fs[i];
      Li:=FL[i][1];
      FisF1:=LF[Li][1] eq i;
      j:=FisF1 select 1 else #LSlopes[Li]-1;
      minv:=MInv[Li][j];
  
      // if the chain is not from F but to F, swap Y<->Pi and X<->X^-1
      if not FisF1 then minv:=Matrix(Z,[[-1,0,0],[0,0,1],[0,1,0]])*minv; end if;
      m:=minv^(-1);
      
      fred,ftransred,ftrans:=TransformAndReduce(monvectors,cfs,red,minv,[X2,Y2,0]);
      fct:=[d: d in FactorisationMPolSplittingField(fred) | d[1] notin [X2,Y2]];
  
      fcompseqns:=[d[1]: d in fct];
      for d in fct do                                  // comps (dim>0)
        //if d[1] in [X2,Y2] then continue; end if;      // from factors of fred on F
        if d[2] gt 1 then FIsSingular[i]:=true; end if;
        newk:=BaseRing(Parent(d[1]));
        if newk ne k then continue trial; end if;
        name:=Sprintf("F%o%o",i,#fct eq 1 select "" else CodeToString(96+Position(fct,d)));
        vprintf Newton,1: "%o %o [%o] delta=%o\n",name,DelSpaces(d[1]),d[2],deltaF[i];
        AddComponent(~comps,name,i,d[1],d[2],d[2] gt 1,m);
      end for;
  
      U:=Curve(A,fred);
      Sing:=SingularSubscheme(U);
      for X in IrreducibleComponents(Sing) do   // singular points (dim=0)
      if Dimension(X) eq 0 then
        S0:=PointsOverSplittingField(X);
        for P in [P: P in S0 | P[1] ne 0 and P[2] ne 0] do 
          FIsSingular[i]:=true;
          EP:=Eltseq(P);
          newk:=Universe(EP);
          if newk ne k then continue trial; end if;
          node:=IsNode(U,P);
          Co:=[j: j in [1..#comps] | (comps[j][2] eq i) and Evaluate(comps[j][3],Eltseq(P)) eq 0];
          vprint Newton,1: Sprintf("%o on %o %o",DelSpaces(P),PrintSequence([comps[j][1]: j in Co]: sep:=","),node select "node" else "not a node");        
          if node then           // ordinary double point
            assert #Co le 2;     // on one or two comps -> chain 
                                 //! unknown length: xy=pi^n with n not implemented
            chain:=["style=zigz",
              Sprintf("$\\!\\!\\!\\!\\!?\\>\\>\\,$%o",deltaF[i],deltaF[i])];
            AddChain(~comps[Co[1]],Co[#Co],chain,true,P,~Style);
          elif #Co eq 1 then
            AddSingularPoint(~comps[Co[1]],P,DelSpaces(P));
          elif #Co eq 2 then
            AddChain(~comps[Co[1]],Co[2],[-1],true,P,~Style);
          else
            name:=Sprintf("F%oS",i);          
            AddComponent(~comps,name,i,1,1,true,DiagonalMatrix(Z,[0,0,0]));     // show as extra component
            for j in Co do 
              AddChain(~comps[j],#comps,[-1],true,P,~Style);
            end for;          
          end if;
        end for;
      end if;
      end for;   
  
      // which 2-face comps degenerate where on 1-faces
   
      Rk3<X3,Y3,P3>:=PR(k,3);
      Rk<t>:=PR(k);
      for L in FL[i] do
      for fci:=1 to #comps do
        fc:=comps[fci];
        if fc[2] ne i then continue; end if;
        fe:=fc[3];
        fv:=Evaluate(fe,[X3,Y3]);
        ft:=Evaluate(TransformPolynomial(fv,MInv[L][1]*m),[t,0,0]);
        for Pj:=1 to #Lpts do
          P:=Lpts[Pj];
          if P`L ne L then continue; end if;
          vprint Newton,1: Sprintf("(%o) %o %o %o -> %o*%o=%o",fci,fc[1],DelSpaces(ft),
            DelSpaces(<P`name,P`L,P`root,P`mult,P`sing>),
            Valuation(ft,t-P`root),fc[4],Valuation(ft,t-P`root)*fc[4]
          );
          Append(~Lvanishing[Pj],<fci,Valuation(ft,t-P`root)*fc[4]>);
        end for;
      end for;
      end for;
  
    end for;
  
    break;
  end for;
  
  // check that orders of vanishing add up
  
  mults1:=[&+[d[2]: d in L]: L in Lvanishing];  
  mults2:=[#LF[d`L]*d`mult: d in Lpts];
  assert mults1 eq mults2;
  
  // chains between comps
  
  for Pi:=1 to #Lpts do
    L,r,mult,sing,name:=Explode([*D`L,D`root,D`mult,D`sing,D`name*]) where D is Lpts[Pi];
    co:=[d: d in Lvanishing[Pi] | d[2] gt 0];
    outer:=#LF[L] eq 1;
    F1:=co[1][1];
    F2:=#co eq 1 select 0 else co[2][1];
    if comps[F1][2] ne LF[L][1] then
      assert F2 ne 0;
      t:=F1; F1:=F2; F2:=t;
    end if;  
    if not sing then 
      assert #LF[L] eq #co; 
      slopes:=LSlopes[L];
      slopes:=slopes[2..#slopes-1];
      slopes:=[Z|delta*Denominator(x): x in slopes] where delta is deltaL[L];
      AddChain(~comps[F1],F2,slopes,false,r,~Style);
      continue;
    end if;
    vprint Newton,1: "singular point:",L,outer,r,mult,co;
    if   #co eq 1 then AddSingularPoint(~comps[F1],r,name);
    elif #co eq 2 then AddChain(~comps[F1],F2,["style=l2sing",name],true,r,~Style);
    else 
      name:=Sprintf("FL%o",L);          
      AddComponent(~comps,name,-Pi,r,mult,true,DiagonalMatrix(Z,[0,0,0]));     // show as extra component
      for c in co do      
        AddChain(~comps[c[1]],#comps,[-1],true,r,~Style);
      end for;
    end if;
  end for;
  
  ////////////////////////////// Dual graph
  
  G:=Graph<#comps|>;
  
  V:=VertexSet(G);
  for i:=1 to #V do AssignLabel(V!i,i); end for;     // vertex labels
  
  D:=AssociativeArray();   // chains and other data
  
  procedure AddData(~D,i,j,c)
    if     IsDefined(D,[i,j]) then Append(~D[[i,j]],c); 
      elif IsDefined(D,[j,i]) then Append(~D[[j,i]],Reverse(c)); 
      else D[[i,j]]:=[* c *];
    end if;
  end procedure;
  
  for i:=1 to #comps do
    name,F,eqn,mult,sing,chains,singpts:=Explode(comps[i]);
    if F gt 0 then mult*:=deltaF[F]; end if;
    if (F le 0) or sing or IsEmpty([P: P in InteriorPoints(Fs[F]) | vP[Eltseq(P)] in Z])
      then g:=0;
      else vprint Newton,1: "Computing geometric genus";
           g:=GeometricGenus(Curve(A,eqn));    
           //! make faster by checking full non-singularity?
    end if;
    data:=[* <"singular",sing>, <"genus",g> *];
    if FaceNames then
      src,srcsty:=FaceNameTeX(name,Style);
      Append(~data,<"source",src,#srcsty lt 2 select "" else srcsty[2]>); 
    end if;
    if (not sing) and (mult ne 1) 
      then Append(~data, <"multiplicity",mult>); end if;
    for P in singpts do
      pn:=Sprint(P[2]);
      repeat
        ok,buf:=Match(pn,"%S$.1^%I%S");
        if not ok then break; end if;
        pn:=Sprintf("%o$r^{%o}$%o",buf[1],buf[2],buf[3]);
      until false;    
      Append(~data,[* Sprintf("[redbull=%o]",#pn eq 1 select pn else "{"*pn*"}"),"$\\bullet$" *]);    
    end for;
    for d in chains do
      toi,chmults,chsing,src:=Explode(d);
      chmultsstr:=[Sprint(x): x in chmults];
      if (toi eq 0) and not IsEmpty(chmults) then 
        if sing then chmults:=["style=l2sing"] cat chmultsstr; end if;
        data cat:= [* chmults *];   // src not used
      elif toi eq i then 
        Append(~data,[* "bluenode" *]);    //! node: chain component to self
      elif toi ne 0 then 
        G+:={i,toi}; 
        if ExtendedType(chmults) ne SeqEnum[MonStgElt] then
          if chsing then chmults:=["style=l2sing"] cat chmultsstr; end if;
        end if;
        AddData(~D,i,toi,chmults);
      end if;
    end for;
    D[[i]]:=data;
  end for;
  
  // contracting chains in the dual graph, corresponding to contractible faces
  
  vprint Newton,1: "Contracting chains";
  
  cmult:=func<Fi|deltaF[comps[Fi][2]]*comps[Fi][4]>;    // main component multiplicity
  for i:=1 to #comps do
    name,F,eqn,mult,sing,chains,singpts:=Explode(comps[i]);
    if F le 0 or not FIsContractible[F] or
       (#[F: F in comps | F[2] eq comps[i][2]] ne 1) then continue; end if; 
    assert not sing;
    ch:=[**];
    for j in [1..#comps] do
//      if IsDefined(D,[i,j]) then ch cat:=[<i,j,eval Sprint(d)>: d in D[[i,j]]]; end if; //! 
//      if IsDefined(D,[j,i]) then ch cat:=[<j,i,eval Sprint(d)>: d in D[[j,i]]]; end if;
/*
"ch",ch,[Sprint(ExtendedType(c[3])): c in ch];
if IsDefined(D,[i,j]) then
  "ij",D[[i,j]],[Sprint(ExtendedType(d)): d in D[[i,j]]];
end if;
if IsDefined(D,[j,i]) then
  "ji",D[[j,i]],[Sprint(ExtendedType(d)): d in D[[j,i]]];
end if;
*/
      if IsDefined(D,[i,j]) then ch cat:=[*<i,j,d>: d in D[[i,j]]*]; end if; //! 
      if IsDefined(D,[j,i]) then ch cat:=[*<j,i,d>: d in D[[j,i]]*]; end if;
    end for;
    assert #ch eq 2; 
    d1,d2:=Explode(ch);
    fromc:=d1[1]+d1[2]-i;
    toc:=d2[1]+d2[2]-i;
    d1chain:=d1[3]; d1s:=not IsEmpty(d1chain) and Type(Universe(d1chain)) eq MonStg;
    d2chain:=d2[3]; d2s:=not IsEmpty(d2chain) and Type(Universe(d2chain)) eq MonStg;
    if d1s and exists{d: d in d1chain | not Match(d,"%i")} then
      ch:=d1chain;    // not an integer in d1 or d2 -> brutally take one of the chains
    elif d2s and exists{d: d in d2chain | not Match(d,"%i")} then
      ch:=d2chain;
    else 
      if d1s then d1chain:=[eval s: s in d1chain]; end if;
      if d2s then d2chain:=[eval s: s in d2chain]; end if;
      ch0:=(d1[2] eq i select d1chain else Reverse(d1chain)) cat [cmult(i)] cat
           (d2[1] eq i select d2chain else Reverse(d2chain));
      ch:=[cmult(fromc)] cat ch0 cat [cmult(toc)];
      assert forall{j: j in [2..#ch-1] | (ch[j-1]+ch[j+1]) mod ch[j] eq 0};  // check self-intersections
      while exists(j){j: j in [2..#ch-1] | ch[j] eq ch[j-1]+ch[j+1]} do    
        Remove(~ch,j);    // blow down
      end while;
      ch:=ch[[2..#ch-1]];
      vprintf Newton,1: "Contracting in D: %o->%o %o %o\n",fromc,toc,DelSpaces(ch0),
        ch eq ch0 select "" else DelSpaces(ch);
    end if;    
    Remove(~D,[i]); 
    Remove(~D,[d1[1],d1[2]]); 
    Remove(~D,[d2[1],d2[2]]);
    G-:=Position(VertexLabels(G),i);
    if fromc ne toc then
      G+:={Position(VertexLabels(G),fromc),Position(VertexLabels(G),toc)};
    end if;
    AddData(~D,fromc,toc,ch);
  end for;
  
  ////////////////////// TeX for Delta_v
  
  vprint Newton,1: "TeX for Delta_v";

  s:="";
 
  // shading: fill singular faces
  
  for Fi in Select(FIsSingular,true) do          
    s*:=Sprintf("\\fill[red!5] %o--cycle;\n",
      PrintSequence([DelSpaces(v): v in Orient(Vertices(Fs[Fi]))]: sep:="--"));
  end for;
  
  // name face components
  
  if FaceNames then
  for i:=1 to #Fs do
    if FIsRemovable[i] or FIsContractible[i] then continue; end if;
    pts:=Points(Fs[i]);
    lat:=Parent(Representative(pts));
    ptsv:={P: P in pts | P+lat![0,1] in pts};
    ptsh:={P: P in pts | P+lat![1,0] in pts};
    ptsq:={P: P in ptsv | P+lat![1,0] in ptsv};
    ptsv:=[P+lat![0,1/2]: P in ptsv | not exists{l : l in Ls[FL[i]] | P+lat![0,1/2] in l}];
    ptsh:=[P+lat![1/2,0]: P in ptsh | not exists{l : l in Ls[FL[i]] | P+lat![1/2,0] in l}];
    names:=[];
    namestyle:="";
    for d in comps do
      if d[2] ne i then continue; end if;
      dname,style:=FaceNameTeX(d[1],Style);
      Append(~names,dname);
      if (#style ne 0) and (style[1] ne "") then namestyle:=","*style[1]; end if;
    end for;
    names:=PrintSequence(names);
    v:=Vs[FV[i]]; 
    c:=&+v/#v;
    d:=&+[(x-c)/  (1+Norm(x-c)) : x in v]/#v;
    if not IsEmpty(ptsq) then 
      ptsq:=[P+lat![1/2,1/2]: P in ptsq];
      _,j:=Min([Norm(c-P): P in ptsq]); c:=ptsq[j];
    elif not IsEmpty(ptsv) then 
      _,j:=Min([Norm(c-P): P in ptsv]); c:=ptsv[j];  
    elif not IsEmpty(ptsh) then 
      _,j:=Min([Norm(c-P): P in ptsh]); c:=ptsh[j];  
    else
      for i:=1 to 5 do
        cn:=c+i*d; 
        if Max([Abs(x-Round(x)): x in Eltseq(cn)]) ge 1/4 then c:=cn; break; end if;  
      end for;
    end if;
    s*:=Sprintf("\\node[fname%o] at (%o,%o) {%o};\n",namestyle,RealField(2)!c.1,RealField(2)!c.2,names);
  end for;
  end if;
  
  // valuations (large or small) at integer points
  
  AllPNonRem:=&join{Set(AllFP[i]): i in [1..#Fs] | not FIsRemovable[i] and not FIsContractible[i]};
  S:=&join[Set(InnFP[i]): i in [1..#Fs] | FIsSingular[i]];
  
  for i:=1 to #AllP do
    P:=AllP[i];
    x,y:=Explode(Eltseq(P));
    z:=vP[P];

    if exists(d){d[2]: d in Style | d[1] eq Sprintf("V%o%o",x,y)} 
      then zstr:=d;
      else zstr:=z;
    end if;

    if     [x,y] in S then             nodestyle:=",red!70!black";
      elif [x,y] notin AllPNonRem then nodestyle:=",rem";
      else nodestyle:="";
    end if;

    s*:=Sprintf("\\node[%o%o] at (%o,%o) (%o) {%o};\n",
      [x,y,z] in mons select (LargeValuations select "lrg" else "coordinate") else 
        SmallValuations select "sml" else "coordinate",
      nodestyle,x,y,i,zstr);   // lrg / sml / coordinate
  end for; 
  
  // paths in Delta_v
  
  vprint Newton,1: "Paths";
  
  s*:="\\draw[lin]\n";
  for i:=1 to #Ls do
    L:=Ls[i];
    delta:=deltaL[i];
    pt:=AllLP[i];
    pti:=[Position(AllP,p): p in pt];
    label:=["": i in [1..#pt-1]];
    
    if forall{f: f in LF[i] | FIsRemovable[f] or FIsContractible[f]} 
        then color:=["[rem]": i in [1..#pt-1]]; else                 // removable/contractible
      BadSegments:=[Z|P`mult*delta: P in Lpts | (P`L eq i) and P`sing];
      BadSegmentNames:=[P`name: P in Lpts | (P`L eq i) and P`sing];
      GoodSegments:=[Z|delta: i in [1..(#pt-1-&+BadSegments) div delta]];  // alternate singular/non-singular
      color:=[];
      goodturn:=#GoodSegments gt #BadSegments;
      repeat
        if IsEmpty(BadSegments) and IsEmpty(GoodSegments) then break; end if;
        if goodturn then
          if not IsEmpty(GoodSegments) then
            for j:=1 to Last(GoodSegments) do Append(~color,""); end for;
            GoodSegments:=Prune(GoodSegments);
          end if;
        else
          if not IsEmpty(BadSegments) then
            for j:=1 to Last(BadSegments) do 
              Append(~color,Sprintf("[red%o]",
                j eq 1 select ",(-" else 
                (j eq Last(BadSegments) select ",-)" 
                else "")));  
              if j eq 1 then label[#color]:=Sprintf("node[lname] {%o} ",Last(BadSegmentNames)); end if;
            end for;
            BadSegments:=Prune(BadSegments);
            BadSegmentNames:=Prune(BadSegmentNames);
          end if;
        end if;
        goodturn:=not goodturn;
      until false;
    end if;
    assert #color eq #pt-1;  
    for i:=1 to #pt-1 do 
      s*:=Sprintf(" (%o) edge%o %o(%o)",pti[i],color[i],label[i],pti[i+1]);
    end for;
    s*:="\n";
  end for;  
  s*:=";\n";
  
  vprint Newton,1: "Contractible chains";
  
  // draw arcs for contractible chains in Delta_v

  for c in ContractibleChains do
    if c[1] eq Last(c) 
      then v:={P: P in InnLP[c[1]] | vP[P] in Z};    // I_n-like from a line to itself
      else v:=&meet[Set(Vertices(L)): L in Ls[c]];   // Rays from one central point
    end if;
    assert #v eq 1;  
    v:=Representative(v);
    v1:=Representative(Set(Vertices(Ls[c[1]])) diff {v});
    v2:=Representative(Set(Vertices(Ls[Last(c)])) diff {v,v1});
    arg1,arg2:=Explode(Sort([angle(w-v): w in [v1,v2]]));
    assert exists(i){i: i in [1..#Fs] | not FIsContractible[i] and not FIsRemovable[i]};
    argc:=angle(FC[i]-v);  
    if (argc ge arg1) and (argc le arg2) then arg2-:=360; end if;
    s*:=Sprintf("\\draw [-] (%o) ++(%o:1.2em) arc (%o:%o:1.2em);\n",Position(AllP,Eltseq(v)),arg1,arg1,arg2); 
  end for;

  sty:=Styles(s,Style);
  if #sty ne 0 then sty:=",\n"*sty; end if;
  charttex:=Sprintf("\\begin{tikzpicture}[%o%o]\n%o\\end{tikzpicture}",
    TikzScale(GetStyle(Style,"DeltaScale")),sty,s);
  
  // Defining equation
  
  vprint Newton,1: "TexPolynomial for defining equation";
  
  if ShowEquation then
    if Rank(Parent(forg)) eq 2 then
      forgstr:=TeXPolynomial(forg: vnames:=["x","y"], vorder:=[2,1]);
    else
      forgstr:=TeXPolynomial(forg: vnames:=["x","y","p"], vorder:=[2,1,3]);
    end if;
    forgstr:=Sprintf(GetStyle(Style,"EquationFormat"),forgstr);
  else 
    forgstr:="";
  end if;
  
  // Charts and transformation matrices
  
  vprint Newton,1: "Charts and transformation matrices";

  if ShowCharts then  
    s:="";
    for i:=1 to #comps do
      name,F,eqn,mult,sing,chains,singpts,m:=Explode(comps[i]);
      if IsZero(m) then continue; end if;
      if /*sing or*/ (F le 0) or FIsContractible[F] or FIsRemovable[F] then continue; end if;
      tr1:=PrintMatrixAsMonomialTransformation(m^(-1),["X","Y","Z"]);
      tr2:=PrintMatrixAsMonomialTransformation(m,["x","y","p"]);
      eqns:=TeXPolynomial(eqn: vorder:=[2,1], vnames:=["X","Y"]);
      s*:=Sprintf("  %o & x=%o & X=%o && %o=0 \\cr\n"* 
                  "     & y=%o & Y=%o && Z^%o=0\\cr\n"* 
                  "     & p=%o & Z=%o \\cr\n",
         FaceNameTeX(name,Style: dollar:=""),tr1[1],tr2[1],eqns,tr1[2],tr2[2],mult*deltaF[F],tr1[3],tr2[3]);
    end for;
    for P in Lpts do
      if not P`sing then continue; end if;
      s*:=Sprintf(" %o & L=%o & r=[%o]^{%o} \\cr\n",P`name,P`L,P`root,P`mult);  
    end for; 
    transformations:=Sprintf(GetStyle(Style,"ChartsFormat"),s);
  else 
    transformations:="";
  end if;
  
  // TeX dual graph 
  vprint Newton,1: "TeX dual graph";
  texregmodel:=TeX(G,D,Style);    

  // Final output  
  charttex:=Sprintf(GetStyle(Style,"DeltaFormat"),charttex);
  texregmodel:=Sprintf(GetStyle(Style,"ModelFormat"),texregmodel);
  texmodel:=Sprintf(GetStyle(Style,"OutputFormat"),forgstr,charttex,texregmodel,transformations);
  
  // Data
  S:=OldStyle;
  S["Fs"]:=Fs;
  S["Ls"]:=Ls;
  S["Vs"]:=Vs;
  S["vF"]:=vF;
  S["FIsRemovable"]:=FIsRemovable;
  S["FIsContractible"]:=FIsContractible;
  S["FIsSingular"]:=FIsSingular;
  S["comps"]:=comps;
  S["Lpts"]:=Lpts;
  S["Lvanishing"]:=Lvanishing;
  S["AllFP"]:=AllFP;
  S["AllLP"]:=AllLP;
  S["InnFP"]:=InnFP;
  S["InnLP"]:=InnLP;
  S["deltaF"]:=deltaF;
  S["deltaL"]:=deltaL;
  S["DVRData"]:=<K,k,vK,red,lift,pi>;
  
  return texmodel,S;
end intrinsic;


intrinsic Model(f::RngMPolElt: Style:=AssociativeArray(), DVR:="default") -> MonStgElt, Assoc
{TeX model for a curve C given by f=0 (main function)}
  return Model(f,"default": Style:=Style, DVR:=DVR);
end intrinsic;

  
intrinsic Model(C::CrvPln, P: Style:=AssociativeArray(), DVR:="default") -> MonStgElt, Assoc
{TeX model for a plane curve}
  F:=DefiningPolynomials(C);
  require #F eq 1: "Curve must be defined by one equation";
  return Model(F[1],P: Style:=Style, DVR:=DVR);
end intrinsic;

  
intrinsic Model(C::CrvHyp, P: Style:=AssociativeArray(), DVR:="default") -> MonStgElt, Assoc
{TeX model for a hyperelliptic curve}
  g,h:=HyperellipticPolynomials(C);
  R<x,y>:=PR(BaseField(C),2);
  f:=y^2+Evaluate(h,x)*y-Evaluate(g,x);
  return Model(f,P: Style:=Style, DVR:=DVR);
end intrinsic;


intrinsic Model(C::CrvEll, P: Style:=AssociativeArray(), DVR:="default") -> MonStgElt, Assoc
{TeX model for an elliptic curve}
  R<x,y>:=PR(BaseField(C),2);
  f:=Evaluate(DefiningPolynomial(C),[x,y,1]);
  return Model(f,P: Style:=Style, DVR:=DVR);
end intrinsic;
