####  Polygons with interior hubs

###################### functions #####################################

# packages and numerical functions
# packages

 with(plots):
 with(numtheory):
 with(combinat):
 
# basic polygon functions
# polygon functions
 NumTriangles := proc(HP)
   RETURN((HP[1]+2*(nops(HP[2])+nops(HP[4])+nops(HP[5])))/3);
   end:
 
 NumBoundaryVertices := proc(HP)
   RETURN(HP[1]);
   end:
 
 NumInteriorVertices := proc(HP)
   RETURN(HP[3]);
   end:
 
 NumVertices := proc(HP)
   RETURN(NumBoundaryVertices(HP)+NumInteriorVertices(HP));
   end:
 
 NumPPdiags := proc(HP)
   RETURN(nops(HP[2]));
   end:
 
 BoundaryEdges := proc(HP)
   local n,k;
   n := NumBoundaryVertices(HP); 
   RETURN({seq({k,1+(k mod n)}, k=1..n)});
   end:
 
 PPdiags := proc(HP)
   RETURN(HP[2]);
   end:
 
 PQspokes := proc(HP)
   RETURN(HP[4]);
   end:
 
 QQspokes := proc(HP)
   RETURN(HP[5]);
   end:
 
 Pend := proc(e)
   RETURN(min(op(e)));
   end:
 
 Qend := proc(e)
   RETURN(max(op(e)));
   end:
 
 Alldiags := proc(HP)
   RETURN({op(PPdiags(HP)),op(PQspokes(HP)),op(QQspokes(HP))});
   end:
 
 
# geometric type 
 mu := proc(L)
   local x;
   RETURN(nops(L)-2-add(1/x,x=L));
   end:  
 
 hyperbolicpoly := proc(L)
   local m,vars,x,subsseq,mi;
   m := mu(L);
   vars := indets(m);
   subsseq := seq(1/x=0,x=vars); 
   mi := subs(x=x,subsseq,m);
   RETURN(evalb(mi >0));
   end:  
 
# polygon drawing 
 DrawHubPoly := proc()
   local HP,a,nP,nQ,Ppnts,Qpnts,Pedges,Pdiags,k,e,PQdiags,QQdiags,
         QX,QY,Qn,p,q,pic,labels,labelopt;
   HP := args(1);
   nP := NumBoundaryVertices(HP); 
   Ppnts := array(1..nP,
              evalf([seq([cos(2*k*Pi/nP),sin(2*k*Pi/nP)],k=1..nP)]));
   labelopt := false;
   for a in args[2..nargs] do
      if type(a,list) then Ppnts := a; fi;
      if type(a,relation) then labelopt := rhs(a); fi;
      od;      
   nQ := NumInteriorVertices(HP);
   if nQ >0 
     then 
      QX := array((1+nP)..(nP+nQ),[seq(0,k=1..nQ)]);
      QY := array((1+nP)..+(nP+nQ),[seq(0,k=1..nQ)]);
      Qn := array((1+nP)..+(nP+nQ),[seq(0,k=1..nQ)]);
      for e in PQspokes(HP) do
        p := Pend(e);
        q := Qend(e);
        Qn[q] := Qn[q]+1;
        QX[q] := QX[q]+Ppnts[p][1];
        QY[q] := QY[q]+Ppnts[p][2];
        od;
      Qpnts := array((1+nP)..(nP+nQ),
                [seq([QX[k]/Qn[k],QY[k]/Qn[k]],k=1+nP..nP+nQ)]);
     else Qpnts := {};
     fi;
   Pedges := {seq([Ppnts[e[1]],Ppnts[e[2]]], e=BoundaryEdges(HP))};
   Pdiags := {seq([Ppnts[e[1]],Ppnts[e[2]]], e=PPdiags(HP))};
   PQdiags := [seq([Ppnts[Pend(e)],Qpnts[Qend(e)]],e=PQspokes(HP))];
   QQdiags := [seq([Qpnts[e[1]],Qpnts[e[2]]],e=QQspokes(HP))];
   pic := plot({op(Pedges),op(Pdiags),op(PQdiags),op(QQdiags)},
                        view =[-1.2..1.2,-1.2..1.2],
                        scaling=constrained,colour = black,axes=none);
   labels := textplot({seq([1.1*Ppnts[i][1],1.1*Ppnts[i][2],``.i],
                       i=1..nP)},
                       view =[-1.2..1.2,-1.2..1.2],
                       scaling=constrained,colour = black,axes=none);
   if labelopt = true 
      then RETURN(display({pic,labels}));
      else RETURN(pic);
      fi;
   end: 
 
 DrawQuad := proc()
   local HP,QA,n,pnts,k,j,QC,QS,e1,e2,s;
   HP := args[1];
   QA := args[2];
   n := NumBoundaryVertices(HP); 
   QC := [];
   QS := array(1..4,[0,0,0,0]);
   k := 0;
   for j from 1 to nops(QA) do 
     if QA[j] <>1 
        then QC := [op(QC),j];
             k := k+1;
        else QS[(k+3 mod 4)+1] := QS[(k+3 mod 4)+1]+1;
        fi; 
     od; 
   QS[4] := n-4-add(QS[k],k=1..3);
   pnts:= array(1..n);
   pnts[QC[1]] := [1,1];
   pnts[QC[2]] := [-1,1]; 
   pnts[QC[3]] := [-1,-1];
   pnts[QC[4]] := [1,-1];
   for j from 1 to 4 do
     s := QS[j];
     e1 := pnts[QC[j]];
     e2 := pnts[QC[(j mod 4)+1]];
     for k from 1 to s do
         pnts[((QC[j]+k+n-1) mod n)+1] := 
                [(e1[1]*(s+1-k)+e2[1]*k)/(s+1),
                 (e1[2]*(s+1-k)+e2[2]*k)/(s+1)]; 
        od;
     od;   
   DrawHubPoly(HP,convert(pnts,list),args[3..nargs]);
   end: 
# basic hub creation
 MakeTriangle := proc();
   RETURN([3,{},0,{},{}]);
   end:
# 
 MakeHub := proc(n)
   local i;
   RETURN([n,{},1,{seq({i,n+1},i=1..n)},{}]);
   end:
 
 MakeHubData := proc()
   global _Hubs,_HubAttachData, _HubDataList,_H64,_H644,
          _H6444,_H84,_H844A,_H844B,_H86,_H104;
   local s,x;
   _H64 := [6,{},2,{{1,7},{2,7},{3,7},{4,7},{5,7},{1,8},{5,8},{6,8}},
            {{7, 8}}];
   _H644 := [6,{},3,{{1,7},{2,7},{4,7},{6,7},{2,8},
                     {3,8},{4,8},{4,9},{5,9},{6, 9}},
             {{7,8}, {7,9}}];          
   _H6444 := [6,{},4,{{1,7},{3,7},{5,7},{1,8},{2,8},{3,8},
                    {3,9},{4,9},{5,9},{1,10},{5,10},{6,10}}, 
              {{7,8},{7,9},{7,10}}];
   _H84 := [8,{},2,{{1,9},{2,9},{8,9},{2,10},{3,10},{4,10},
                    {5,10},{6,10},{7,10},{8,10}}, 
            {{9, 10}}];
   _H844A := [8,{},3,{{1,9},{2,9},{8,9},{2,10},{3,10},{4,10},
                       {5,10},{6,10},{8,10},{6,11},{7,11},{8,11}},
              {{9,10},{10,11}}];
   _H844B := [8,{},3,{{1,9},{2,9},{8,9},{2,10},{3,10},{4,10},
                      {6,10},{7,10},{8,10},{4,11},{5,11},{6,11}},
              {{9,10},{10,11}}];
 
   _H86 := [10,{},2,{{1,11},{2,11},{3,11},{9,11},{10,11},{3,12},
                    {4,12},{5,12},{6,12},{7,12},{8,12},{9,12}},
            {{11,12}}]; 
   _H104 := [10,{},2,{{2,11},{3,11},{4,11},{5,11},
                      {6,11},{7,11},{8,11},{9,11},{10,11},
                      {1,12},{2,12},{10,12}},
             {{11,12}}];
   _Hubs := [MakeTriangle(),MakeHub(4),MakeHub(6),MakeHub(8),
             MakeHub(10),MakeHub(12),_H64,_H644,_H6444,_H84,
             _H844A,_H844B,_H86,_H104];
   _HubAttachData :=[[MakeTriangle(),[1]],[MakeHub(4),[1]],
                 [MakeHub(6),[1]],[MakeHub(8),[1]],
                 [MakeHub(10),[1]],[MakeHub(12),[1]],
                 [_H64,[1,2,3,4,5,6]],[_H644,[1,2,3,4,5,6]],
                 [_H6444,[1,2]],[_H84,[1,2,3,4,5,6,7,8]],
                 [_H844A,[1,2,3,4,5,6,7,8]],[_H844B,[1,2,3,4]],
                 [_H86,[1]],[_H104,[1]]];
   _HubDataList := array(1..12,[seq([],i=1..12)]);
   for x in _HubAttachData do      
       s := NumTriangles(x[1]);
       _HubDataList[s] := [op(_HubDataList[s]),x];
       od; 
   RETURN(NULL):
   end:
# 
# polygon transformation and construction
# relabelling
 TranslateDiagsandSpokes := proc(HP,Ptab,Qtab)
   local newPdiags,newPQspokes,newQQspokes,e,nP;
   nP := NumBoundaryVertices(HP);
   newPdiags := {seq({Ptab[e[1]],Ptab[e[2]]},e=HP[2])};
   newPQspokes := {seq({Ptab[Pend(e)],Qtab[Qend(e)-nP]},e=HP[4])};
   newQQspokes := {seq({Qtab[e[1]-nP],Qtab[e[2]-nP]},e=HP[5])};
   RETURN([HP[1],newPdiags,HP[3],newPQspokes,newQQspokes]);
   end:
 
# construction
 JoinPolys := proc(HP1,k1,HP2,k2)
   local nP1,nQ1,Ptab1,Qtab1,THP1,nP2,nQ2,Ptab2,Qtab2,
         THP2,newnP,j,HP3,x;
   nP1 := NumBoundaryVertices(HP1);
   nP2 := NumBoundaryVertices(HP2);
   nQ1 := NumInteriorVertices(HP1);
   nQ2 := NumInteriorVertices(HP2);
   newnP := nP1+nP2-2; 
   Ptab1 := [seq(j,j=1..k1),seq(j+nP2-2,j=(k1+1)..nP1)];
   Qtab1 := [seq(j+newnP,j=1..nQ1)];
   Ptab2 := [seq(j-1+k1+nP2-k2, j=1..k2),seq(j-1+k1, j=1..(nP2-k2))];
   Ptab2 := [seq(x-1+newnP mod newnP +1,x=Ptab2)];
   Qtab2 := [seq(j+newnP+nQ1,j=1..nQ2)]; 
   THP1:= TranslateDiagsandSpokes(HP1,Ptab1,Qtab1);
   THP2:= TranslateDiagsandSpokes(HP2,Ptab2,Qtab2);
   HP3 := [nP1+nP2-2,
           {op(PPdiags(THP1)),op(PPdiags(THP2)),
            {k1,(k1+nP2-2) mod (newnP)+1}},
           nQ1+nQ2,
           {op(PQspokes(THP1)),op(PQspokes(THP2))},
           {op(QQspokes(THP1)),op(QQspokes(THP2))}
          ];
   RETURN(HP3);  
 end:
 
# dihedral action
 DihedralTable:= proc(n)
   local k,i,rots,refls;
   rots := [seq([seq(((i+k-1) mod n)+ 1, i=1..n)],k=0..n-1)];
   refls := [seq([seq(((-i+k-1) mod n)+ 1, i=1..n)],k=0..n-1)];
   RETURN([op(rots),op(refls)]);
   end:
 
 PermutePoly := proc(HP,Ptab)
   local Qtab,nP,nQ,j;
   nP := NumBoundaryVertices(HP);
   nQ := NumInteriorVertices(HP);
   Qtab := [seq(j+nP,j=1..nQ)];
   RETURN(TranslateDiagsandSpokes(HP,Ptab,Qtab)); 
   end:
   
 DihedralOrbit:= proc(HP,DHtab)
   local Ptab;
   RETURN({seq(PermutePoly(HP,Ptab),Ptab=DHtab)});
   end:
 
 DHOrbitSpace:= proc(HPL)
   local OS,LHPL,DHtab,n,HP,DHO;
   LHPL := {op(HPL)};
   OS := [];
   DHtab := DihedralTable(NumBoundaryVertices(LHPL[1]));
   while nops(LHPL) >0 do
     HP := LHPL[1];
     OS := [op(OS),HP];
     DHO := DihedralOrbit(HP,DHtab);
     LHPL := LHPL minus DHO;
     od;     
   RETURN(OS);
   end:
 
# polygon analysis 
# vertex type analysis
 ConnectedVertices := proc(edgeset,v)
   local VL,e;
   VL := {};
   for e in edgeset do
       if member(v,e) then VL := VL union e; fi;
       od;
   RETURN(VL minus {v});
   end:
 
 ThirdVertexSet := proc(edgeset,v1,v2)
   if not member({v1,v2},edgeset) then RETURN({}); fi;
   RETURN(ConnectedVertices(edgeset,v1) intersect
          ConnectedVertices(edgeset,v2));
   end:
 
 VertexTypes := proc(HP)
   local i,j,k,verts,vtypes,t,u,v,w,TVS,es,nV;
   es := {op(BoundaryEdges(HP)),op(Alldiags(HP))};
   nV := NumBoundaryVertices(HP)+NumInteriorVertices(HP);
   vtypes := array(1..nV);
   vtypes[1]:= 1; vtypes[2] := 2;
   verts := {1,2};
   while nops(verts) < nV do;       
         for u in verts do
         for v in verts minus {u} do
             TVS := ThirdVertexSet(es,u,v);
             t := 6-vtypes[u]-vtypes[v];
             for w in TVS do 
                 vtypes[w] := t;
                 verts := {op(verts),w}; 
                od;
             od;od;
         od;               
    RETURN([seq(vtypes[k],k=1..nV)]);
   end:
 
# angle functions
 NumAnglesatBoundaryVertex := proc(HP,n)
   local diag,s;
   s := 1;
   for diag in Alldiags(HP) do
       if member(n,diag) then s := s+1; fi;
       od;
   RETURN(s);
   end:
 
 NumAnglesatInteriorVertex := proc(HP,n)
   local diag,s;
   s := 0;
   for diag in Alldiags(HP) do
       if member(n,diag) then s := s+1; fi;
       od;
   RETURN(s);
   end:
 
 NumAnglesatVertexList := proc(HP)
   local n,nP,nQ;
   nP := NumBoundaryVertices(HP);
   nQ := NumInteriorVertices(HP);
   RETURN([seq(NumAnglesatBoundaryVertex(HP,n),n=1..nP),
           seq(NumAnglesatInteriorVertex(HP,n+nP),n=1..nQ)]);
   end:
 
# angle multiplicities and corner constraints
 VertexTypeMults := proc(HP)
   local nP,nQ,AML,VTL,s,t,i,j,k;
   nP := NumBoundaryVertices(HP);
   nQ := NumInteriorVertices(HP);
   AML := NumAnglesatVertexList(HP);
   VTL := VertexTypes(HP);
   s := array(1..6,[[],[],[]]); 
   t := array(1..6,[[],[],[]]); 
   for i from 1 to nP do
      j := VTL[i];
      s[j] := [op(s[j]),AML[i]]; 
      od;
    for i from 1+nP to nP+nQ do
      j := VTL[i];
      t[j] := [op(t[j]),AML[i]]; 
      od;
   RETURN([seq(s[k],k=1..3),seq(t[k],k=1..3)]);
   end:
 
 PossibleNumCorners := proc(BML,IML)
   local L,m,t,s;
   if nops(BML)= 0 then RETURN([0]); fi;
   if nops(IML) >0 
     then 
      m := 0;
      for s in BML do
        if IML[1]/s >2 then m := m+1; fi;
        od;
      RETURN([m]);
      fi;
   L := lcm(op(BML));
   if L=1 then RETURN([nops(BML)]); fi;
   m := 0;
   for t in BML do if t<L then m := m+1; fi; od;
   if m = nops(BML) 
      then RETURN([nops(BML)]);
      else RETURN([m,nops(BML)]); 
      fi;
    end:
 
 MinPolyCorners := proc(HP)
   local VTL,ML,L,i;
   VTL := VertexTypeMults(HP);
   L := [seq(PossibleNumCorners(VTL[i],VTL[i+3]),i=1..3)];
   RETURN(add(L[i][1],i=1..nops(L)));
   end:
 
# compatibilty and validity functions
 ValidMults := proc(VML)
   local nP,nQ,AML,VTL,s,t,i,j,k;
   for j from 4 to 6 do 
      if nops({op(VML[j])}) >1 then RETURN(false); fi; 
      od;
   for j from 1 to 3 do 
      for s in VML[j] do 
      for t in VML[j+3] do 
         if  (((t mod s) <>0) or (t=s)) then RETURN(false); fi; 
         od;
         od;
      od;
   RETURN(true);
   end:
 
 ValidPairVertexAnglesList := proc(HP,L)
   local nP,nQ,VTL,nAVL,AVL,i,eqns,sol,lL,vL;
   lL := convert(L,list);
   nP := NumBoundaryVertices(HP);
   nQ := NumInteriorVertices(HP);
   VTL := VertexTypes(HP);
   nAVL := NumAnglesatVertexList(HP); 
   AVL := [seq(nAVL[i]/lL[VTL[i]],i=1..nP+nQ)];
   eqns := {seq(AVL[i]=2,i=(1+nP)..(nP+nQ))};
   sol := solve(eqns);
   AVL := simplify(subs(sol,AVL));
   vL := simplify(subs(sol,lL));
   RETURN([vL,[seq(AVL[i],i=1..nP)]]);
   end:
 
 
 ValidPolys := proc(HP,parms)
   local lcmL,VTL,MML,VAL,VA,L,a,i,j,k,stuv,polypairs,lmn,vlmn,QA;
   VTL := VertexTypeMults(HP);
   if not ValidMults(VTL) then RETURN([]); fi; 
   MML := [seq(PossibleNumCorners(VTL[i],VTL[i+3]),i=1..3)];
   lcmL := [seq(lcm(op(VTL[i])),i=1..3)];
   lmn := array(1..3);
   polypairs := [];  
   for i from 1 to nops(MML[1]) do
   for j from 1 to nops(MML[2]) do
   for k from 1 to nops(MML[3]) do
     if (MML[1][i]+MML[2][j]+MML[3][k] = 4)
        then
             if i = nops(MML[1])
                then lmn[1] := lcmL[1]*parms[1]; 
                else lmn[1] := lcmL[1] 
                fi; 
             if j = nops(MML[2]) 
                then lmn[2] := lcmL[2]*parms[2]; 
                else lmn[2] := lcmL[2] 
                fi;
             if k = nops(MML[3]) 
                then lmn[3] := lcmL[3]*parms[3]; 
                else lmn[3] := lcmL[3] 
                fi; 
      VAL := ValidPairVertexAnglesList(HP,lmn);
      vlmn := VAL[1];
      VA := VAL[2];
             stuv := [];
             for a in VA do if a<>1 then stuv :=[op(stuv),1/a]; fi od;
             if (nops(stuv) = 4) and hyperbolicpoly(stuv)
                then polypairs := [op(polypairs),[vlmn,stuv,VA]];
                     fi;
             fi;
      od;od;od;
   RETURN(polypairs);
   end:


###################### driver ################
restart;   

# main loop - polygon creation
 
 MakeHubData():
 Hpolys.4 := {MakeHub(4)}:
 Horbspa.4 := Hpolys.4:
 save Horbspa.4,hpolyorbs.4.`.mpl`: 
 
 for K from 4 to top do
   print(`------------------------------------------------`);
   print(`numtriangs ` = K);
   print(`finding polys: `);
   st := time():
   Hpolys.K := [seq(X[1],X=_HubDataList[K])];
     for s from 1 to K-4 do
     for X in _HubDataList[s] do
       aHP := X[1]; SL := X[2];  
       for HP in Horbspa.(K-s) do
         nV := NumBoundaryVertices(HP);
         for k from 1 to nV do
         for l in SL do
            Hpolys.K := {op(Hpolys.K),JoinPolys(HP,k,aHP,l)};
            od; 
            od;
         od;
       od;
     od;
   print([`numpolys found ` = nops(Hpolys.K),
          `time taken `= time()-st]);
  
   print(`finding orbits `);
   st := time():
   Horbspa.K := [];
   vn := [op({seq(NumBoundaryVertices(HP2), HP2=Hpolys.K)})];
   print(`possible # of vertices `= vn);
   for u in vn do
     PS := [];
     for HP3 in Hpolys.K do
        if NumBoundaryVertices(HP3) = u then PS :=[op(PS),HP3];fi;
        od; 
     Horbspa.K := [op({op(Horbspa.K),op(DHOrbitSpace(PS))})];
     print([`orbits so far `= nops(Horbspa.K),`time taken `= time()-st]);
     od;
   save Horbspa.K,hpolyorbs.K.`.mpl`; 
   od: 
 
# main loop - quadrilateral search
  otop := top;
  print("Qrad Search");
  unassign('l','m','n');
 for K from 4 to otop do
 print(`------------------------------------------------`);
 print(`# triangs ` = K); 
   read hpolyorbs.K.`.mpl`;
   OS := Horbspa.K:
   VPL :=[];
   for HP in OS do 
        VP := ValidPolys(HP,[l,m,n]);
        if nops(VP)>0 
           then 
            print([`vertextypes` = VertexTypes(HP)]);     
            for P in VP do 
              print(P); 
              VPL := [op(VPL),seq([HP,op(x)],x=VP)];            
              od;          
        fi;
      od;
   save VPL, hquads.K.`.mpl`;
  od: