1# (C)2007-2008 by Marc Roeder,
2#   distribute under the terms of the GPL version 2.0 or later
3
4
5#############################################################################
6## This file provides functions to generate fundamental domains for
7## 3 dimensional Bieberbach groups and view them as colored polytopes in
8## JavaView.
9##
10## The functions you might want to use are in 3dimBieberbachFD.gap
11## Look there first.
12#############################################################################
13
14
15#############################################################################
16##
17## A square root approximation using continued fractions.
18##  Used for the Cholesky decomposition.
19
20
21    sqrt:=function(rat)
22        local   x,denom,  num;
23        denom:=DenominatorRat(rat);
24        num:=NumeratorRat(rat);
25        x:=Indeterminate(Rationals);
26        return ContinuedFractionApproximationOfRoot(denom*x^2-num,100);
27    end;
28
29
30#############################################################################
31##
32## A quick and dirty implementation of a Cholesky decomposition (approximation)
33##
34
35CholeskyDecomposition:=function(mat)
36    local   dim, L,  line,x,  col;
37
38    if not mat=TransposedMat(mat)
39       then
40        Error("matrix is not symmetric");
41    fi;
42    dim:=DimensionSquareMat(mat);
43    if not ForAll([1..dim],d->Determinant(mat{[1..d]}{[1..d]})>0)
44       then
45        Error("matrix is not positive definite");
46    fi;
47    L:=NullMat(dim,dim);
48    L[1][1]:=sqrt(mat[1][1]);
49    for line in [1..dim]
50      do
51        for col in [1..line-1]
52          do
53            if col>1
54               then
55                L[line][col]:=(mat[line][col]-(L[line]{[1..col-1]}*L[col]{[1..col-1]}))/L[col][col];
56            else
57                L[line][col]:=mat[line][col]/L[col][col];
58            fi;
59            if line>1
60               then
61                L[line][line]:=sqrt(mat[line][line]-L[line]{[1..line-1]}^2);
62            fi;
63        od;
64    od;
65    return L;
66end;
67
68
69
70#############################################################################
71## convert hsv color space to rgb color space.
72## hsv is a triple [h,s,v] where
73## h\in [0..259], s,v\in[0..255].
74##
75hsv2rgb:=function(hsv)
76    local   H,  S,  V,  Hi,  f,  p,  q,  t,  rgb;
77
78    H:=hsv[1];
79    S:=1/256*hsv[2];
80    V:=1/256*hsv[3];
81    Hi:=Int(H/60) mod 6;
82    f:=H/60-Hi;
83    p:=V*(1-S);
84    q:=V*(1-f*S);
85    t:=V*(1-(1-f)*S);
86    if Hi=0
87       then
88        rgb:=[V,t,p];
89    elif Hi=1
90      then
91        rgb:= [q,V,p];
92    elif Hi=2
93      then
94        rgb:= [p,V,t];
95    elif Hi=3
96      then
97        rgb:= [p,q,V];
98    elif Hi=4
99      then
100        rgb:= [t,p,V];
101    elif Hi=5
102      then
103        rgb:= [V,p,q];
104    else
105        Error("bad programming, Marc");
106    fi;
107    return List(rgb*256,Int);
108end;
109
110#############################################################################
111## create a list of length n which contains strings of the form
112##  "xxx yyy zzz" representing a point in rgb-colorspace.
113##
114ncolorStrings:=function(n)
115    local   colorlist,  Hvalues,  color,  rem,  H,  S,  V;
116
117    if n<=10
118       then
119        colorlist:=List([0..n-1]*360/n,i->[Int(i),255,255]);
120    else
121        colorlist:=[];
122        Hvalues:=Reversed([0..Int(n/3)+1]*360/(Int(n/3)+1));
123        for color in [1..n]
124          do
125#                        H:=(360/(Int(n/4)+1))*Int(color/4);
126            rem:=(color-1) mod 3;
127            if rem=0
128               then
129                H:=Remove(Hvalues);
130                S:=255;
131                V:=255;
132            elif rem=1
133              then
134                S:=127;
135                V:=255;
136            else
137                S:=255;
138                V:=127;
139            fi;
140           Add(colorlist,[H,S,V]);
141        od;
142    fi;
143    Apply(colorlist,hsv2rgb);
144    Apply(colorlist,c->List(c,String));
145    Apply(colorlist,c->JoinStringsWithSeparator(c," "));
146    return colorlist;
147end;
148
149
150##########################
151## the same, returning a slightly paler set of colors.
152##
153ncolorStringsPale:=function(n)
154    local   colorlist,  Hvalues,  color,  rem,  H,  S,  V;
155
156    if n<=10
157       then
158        colorlist:=List([0..n-1]*360/n,i->[Int(i),180,200]);
159    else
160        colorlist:=[];
161        Hvalues:=Reversed([0..Int(n/3)+1]*360/(Int(n/3)+1));
162        for color in [1..n]
163          do
164            #            H:=(360/(Int(n/4)+1))*Int(color/4);
165            rem:=(color-1) mod 3;
166            if rem=0
167               then
168                H:=Remove(Hvalues);
169                S:=127;
170                V:=255;
171            elif rem=1
172              then
173                S:=64;
174                V:=255;
175            else
176                S:=127;
177                V:=127;
178            fi;
179            Add(colorlist,[H,S,V]);
180        od;
181    fi;
182    Apply(colorlist,hsv2rgb);
183    Apply(colorlist,c->List(c,String));
184    Apply(colorlist,c->JoinStringsWithSeparator(c," "));
185    return colorlist;
186end;
187
188
189
190#############################################################################
191## Convert a rational number into a "floating point" string.
192## <precision> gives the number of digits after the first non-zero
193## digit after the point.
194## Example with precision=3: 1/7 is converted to "0.142"
195##                            1/700 is converted to "0.00142"
196##
197
198fraction2floatString:=function(q,precision)
199    local   sign,  signString,  magnitude,  beforepoint,  qright,
200            prezeros,  afterpoint,  returnstring;
201
202    if q=0
203       then
204        return "0";
205    fi;
206    sign:=SignRat(q);
207    if sign=-1
208       then
209        signString:="-";
210    else
211        signString:="";
212    fi;
213
214    if AbsoluteValue(q)>0
215       then
216        beforepoint:=String(sign*Int(q));
217        qright:=sign*(q-Int(q));
218    else
219        qright:=sign*q;
220        beforepoint:="0";
221    fi;
222
223    if qright<>0
224       then
225        magnitude:=LogInt(NumeratorRat(qright),10)-LogInt(DenominatorRat(qright),10);
226        if AbsoluteValue(qright)*10^(-magnitude)<1
227           then
228            magnitude:=magnitude-1;
229        fi;
230        prezeros:=Concatenation(List([1..-magnitude-1],i->"0"));
231        afterpoint:=String(Int(10^(precision-magnitude)*qright));
232        returnstring:=Concatenation([signString,beforepoint,".",prezeros,afterpoint]);
233    else
234        returnstring:=Concatenation(signString,beforepoint,".0");
235    fi;
236    while returnstring[Size(returnstring)]='0'
237      do
238        Unbind(returnstring[Size(returnstring)]);
239    od;
240    if returnstring[Size(returnstring)]='.'
241       then
242        Unbind(returnstring[Size(returnstring)]);
243    fi;
244    return returnstring;
245end;
246
247
248#############################################################################
249## This converts a number to a string and fills it up with leading zeros
250## if necessary.
251## Example with digits=3: 5->"005"
252##                        1234->"1234"
253##
254numberWithLeadingZeros:=function(n,digits)
255    local   string;
256
257    string:=String(n);
258    while Size(string)<digits
259       do
260        string:=Concatenation(["0",string]);
261    od;
262    return string;
263end;
264
265
266
267###################################
268
269vertexOrbitDecomposition:=function(vertexlist,group)
270    local   vertices,  vertexorbits,  vertex,  orbit;
271
272    vertices:=Set(vertexlist);
273    vertexorbits:=[];
274    while vertices<>[]
275      do
276        vertex:=vertices[1];
277        orbit:=Concatenation(
278                       OrbitPartInVertexSetsStandardSpaceGroup(group,[vertex],
279                               vertices));
280        Add(vertexorbits,orbit);
281        SubtractSet(vertices,orbit);
282    od;
283    Apply(vertexorbits,i->List(i,j->Position(vertexlist,j)));
284
285    return vertexorbits;
286end;
287
288###################################
289
290edgeOrbitDecomposition:=function(edges,vertexlist,group)
291    local   edgesasvectors,  edgeorbits,  edge,  orbit;
292
293    edgesasvectors:=Set(edges,i->Set(i,j->vertexlist[j]));
294    edgeorbits:=[];
295    while edgesasvectors<>[]
296      do
297        edge:=edgesasvectors[1];
298        orbit:=Set(OrbitPartInVertexSetsStandardSpaceGroup(group,edge,
299                       Set(vertexlist)));
300        SubtractSet(edgesasvectors,orbit);
301        Add(edgeorbits,orbit);
302    od;
303    edgeorbits:=Set(edgeorbits);
304    Apply(edgeorbits,i->Set(i,j->Set(j,k->Position(vertexlist,k))));
305    return edgeorbits;
306end;
307
308
309#############################################################################
310## Take a string and put <tag> and </tag> arund it
311##
312enclosedByTag:=function(string,tag)
313    return Concatenation(["<",tag,">",string,"</",tag,">\n"]);
314end;
315
316enclosedByTagNN:=function(string,tag)
317    return Concatenation(["<",tag,">",string,"</",tag,">"]);
318end;
319
320
321enclosedByTagWithOptions:=function(string,tag,options)
322    return Concatenation(["<",tag," ",options,">\n",string,"</",tag,">\n"]);
323end;
324
325####################
326
327javaviewWrappedDatastring:=function(title,abstract,detail,datastring,wrapper)
328    local   outstring,  stream,  line;
329
330    outstring:=[];
331    stream:=InputTextFile(wrapper);
332    while not IsEndOfStream(stream)
333      do
334        line:=ReadLine(stream);
335        if line="##Insert Title##\n"
336           then
337            line:=enclosedByTag(title,"title");
338            Append(line,"\n");
339        elif line ="##Insert Abstract##\n"
340          then
341            line:=enclosedByTag(abstract,"abstract");
342            Append(line,"\n");
343
344        elif line ="##Insert Detail##\n"
345          then
346            line:=enclosedByTag(detail,"detail");
347            Append(line,"\n");
348        elif line="##Insert Data String##\n"
349          then
350            line:=enclosedByTag(datastring,"geometries");
351            Append(line,"\n");
352        else
353        fi;
354        Append(outstring,line);
355    od;
356    CloseStream(stream);
357    return outstring;
358end;
359
360####################
361
362    isposdef:=function(mat)
363        return ForAll([1..Size(mat)],d->Determinant(mat{[1..d]}{[1..d]})>0);
364    end;
365
366
367#############################################################################
368## This generates the vertices for JavaView:
369##  The coloring is calculates separately.
370    ############
371
372    javaviewJustPointBlock:=function(vertices,group,precision)
373        local   vector2floatString,  cd,  coordstrings,  pstring;
374
375        vector2floatString:=function(v,precision)
376            local coordstring;
377            coordstring:= List(v,q->fraction2floatString(q,precision));
378            return JoinStringsWithSeparator(coordstring," ");
379        end;
380
381    cd:=CholeskyDecomposition(GramianOfAverageScalarProductFromFiniteMatrixGroup(PointGroup(group)));
382    coordstrings:=List(vertices,i->vector2floatString(i*cd,precision));
383    pstring:=Concatenation([List(coordstrings,s->enclosedByTag(s,"p")),
384                     [enclosedByTag("3","thickness")],
385                     ["\n"]]);
386    pstring:=Concatenation(pstring);
387    return enclosedByTag(pstring,"points");
388end;
389
390
391    javaviewPointColors:=function(vertices,group)
392        local   porbits,  pcolors,  pcolorlist;
393
394        porbits:=vertexOrbitDecomposition(vertices,group);
395
396        pcolors:=ncolorStrings(Size(porbits));
397        pcolorlist:=List([1..Size(vertices)],v->pcolors[PositionProperty(porbits,o->v in o)]);
398        return enclosedByTag(
399                          Concatenation(List(pcolorlist,s->enclosedByTag(s,"c"))),
400                      "colors");
401
402end;
403
404
405#############################################################################
406## And the edges
407################
408
409    javaviewEdgeBlock:=function(poly,group)
410        local   int2vec,  vec2int,  vertices,  facelattice,  edges,
411                edgeorbits,  orientededges,  orbit,  firstedge,  edge,
412                vecedge,  map,  edgeimage,  estring,  ecolors,
413                ecolorlist,  ecstring;
414
415        int2vec:=function(n)
416            return vertices[n];
417        end;
418        vec2int:=function(v)
419            return Position(vertices,v);
420        end;
421        vertices:=Polymake(poly,"VERTICES");
422        facelattice:=Polymake(poly,"FACE_LATTICE");
423        edges:=Set(facelattice[Size(facelattice)-1]);
424        edgeorbits:=edgeOrbitDecomposition(edges,vertices,group);
425
426        orientededges:=[];
427        for orbit in edgeorbits
428          do
429            firstedge:=Set(orbit[1],int2vec);
430            Add(orientededges,List(firstedge,vec2int));
431            for edge in orbit{[2..Size(orbit)]}
432              do
433                vecedge:=Set(edge,int2vec);
434                map:=RepresentativeActionOnRightOnSets(group,firstedge,vecedge);
435                edgeimage:=List(firstedge,v->(Concatenation(v,[1])*map));
436                Add(orientededges,List(edgeimage,i->vec2int(i{[1..3]})));
437            od;
438        od;
439
440        orientededges:=List(edges,e->First(orientededges,i->Set(i)=Set(e)));
441
442        estring:=List(orientededges,e->List(e-1,String));
443        Apply(estring,s->JoinStringsWithSeparator(s," "));
444        Apply(estring,s->enclosedByTag(s,"l"));
445        estring:=enclosedByTag(
446                         Concatenation(Concatenation(estring),
447                             enclosedByTag("4","thickness"))
448                             ,"lines");
449
450                         ecolors:=ncolorStrings(Size(edgeorbits));
451        ecolorlist:=List(edges,v->ecolors[PositionProperty(edgeorbits,o->v in o)]);
452        ecstring:=enclosedByTag(
453                          Concatenation(List(ecolorlist,s->enclosedByTag(s,"c"))),
454                          "colors");
455
456        return enclosedByTagWithOptions(Concatenation(["\n",estring,ecstring]),
457                       "lineSet", "line=\"show\" color=\"show\" arrow=\"show\"");
458     end;
459
460#############################################################################
461## FACETS
462
463
464javaviewFacetBlock:=function(poly,group)
465    local   crossProduct,  facets,  facelattice,  edges,  ofacets,
466            facet,  facetedges,  ofacet,  nextvertex,  nextedge,
467            vertices,  pos,  v1,  v2,  normal,  eq,  testpoint,
468            fstring,  facetorbits,  fcolors,  fcolorlist,  fcstring;
469
470    crossProduct:=function(x,y)
471        return [  x[2]*y[3]-x[3]*y[2],
472                  -(x[1]*y[3]-x[3]*y[1]),
473                  x[1]*y[2]-x[2]*y[1]];
474    end;
475
476    Polymake(poly,"VERTICES_IN_FACETS FACE_LATTICE");
477    facets:=Polymake(poly,"VERTICES_IN_FACETS");
478    ## first, generate "ordered" versions of the facets:
479    facelattice:=Polymake(poly,"FACE_LATTICE");
480    edges:=Set(facelattice[Size(facelattice)-1]);
481    ofacets:=[];
482    for facet in facets
483      do
484        facetedges:=Set(Filtered(edges,e->IsSubset(facet,e)));
485        ofacet:=List(facetedges[1]);
486        RemoveSet(facetedges,ofacet);
487        nextvertex:=ofacet[2];
488        while facetedges<>[]
489          do
490            nextedge:=First(facetedges,f->nextvertex in f);
491            RemoveSet(facetedges,nextedge);
492            nextvertex:=First(nextedge,i->i<>nextvertex);
493            Add(ofacet,nextvertex);
494        od;
495        Remove(ofacet);
496        Add(ofacets,ofacet);
497    od;
498
499    ## now check if the ordering must be reversed:
500    ## We're lucky this happens in 3-space.
501
502    vertices:=Polymake(poly,"VERTICES");
503   for pos in [1..Size(ofacets)]
504      do
505        v1:=vertices[ofacets[pos][2]]-vertices[ofacets[pos][1]];
506        v2:=vertices[ofacets[pos][3]]-vertices[ofacets[pos][1]];
507        normal:=crossProduct(v1,v2);
508        eq:=Concatenation([vertices[ofacets[pos][1]]*normal],normal);
509        testpoint:=vertices[Representative(Difference([1..Size(vertices)],ofacets[pos]))];
510        if WhichSideOfHyperplane(testpoint,eq)=-1
511           then
512            ofacets[pos]:=Reversed(ofacets[pos]);
513        fi;
514    od;
515
516    fstring:=List(ofacets,e->List(e-1,String));
517    Apply(fstring,s->JoinStringsWithSeparator(s," "));
518    Apply(fstring,s->enclosedByTag(s,"f"));
519    fstring:=enclosedByTag(Concatenation(fstring) ,"faces");
520    facetorbits:=edgeOrbitDecomposition(facets,vertices,group);
521
522    fcolors:=ncolorStringsPale(Size(facetorbits));
523    fcolorlist:=List(facets,v->fcolors[PositionProperty(facetorbits,o->v in o)]);
524    fcstring:=enclosedByTag(
525                      Concatenation(List(fcolorlist,s->enclosedByTag(s,"c"))),
526                      "colors");
527    return enclosedByTagWithOptions(Concatenation(["\n",fstring,fcstring]),
528            "faceSet", "face=\"show\" edge=\"show\" color=\"show\"");
529end;
530
531
532
533
534#############################################################################
535## Finally,
536## generate a nice javaview file.
537    ## It takes a number of affine matrices and generates the
538    ## image under each of them.
539    ###########################
540
541
542javaviewDatastring:=function(poly,maps,group,precision)
543    local   vertices,  pstring,  pcstring,  edgeblock,  facetblock,
544            allgeometries,  i,  gshow,  numberstring,  m,  pointblock,
545            facetgeometry,  edgegeometry;
546
547    vertices:=Polymake(poly,"VERTICES");
548    pstring:=javaviewJustPointBlock(vertices,group,precision);
549    pcstring:=javaviewPointColors(vertices,group);
550    edgeblock:=javaviewEdgeBlock(poly,group);
551    facetblock:=javaviewFacetBlock(poly,group);
552
553    allgeometries:=[];
554
555    for i in [1..Size(maps)]
556      do
557        if i=1
558           then
559            gshow:="\"show\"";
560            if Size(maps)=1
561               then
562                numberstring:="";
563            else
564                numberstring:=" FD";
565            fi;
566        else
567            gshow:="\"hide\"";
568            numberstring:=String(i-2);
569        fi;
570        m:=maps[i];
571        pstring:=javaviewJustPointBlock(List(vertices,v->
572                         List(Concatenation(v,[1])*m){[1..3]}),
573                         group,
574                         precision);
575        pointblock:=enclosedByTagWithOptions(
576                            Concatenation(["\n",pstring,pcstring]),
577                            "pointSet",
578                            "dim=\"3\" point=\"show\" color=\"show\" "
579                            );
580        facetgeometry:=enclosedByTagWithOptions(
581                               Concatenation(["\n",pointblock,"\n",facetblock,"\n"]),
582                           "geometry",
583                               Concatenation("name=\"points and facets ",numberstring,"\""," visible=",gshow)
584                               );
585        edgegeometry:=enclosedByTagWithOptions(
586                              Concatenation(["\n",pointblock,"\n",edgeblock,"\n"]),
587                              "geometry",
588                              Concatenation("name=\"points and edges ",numberstring,"\""," visible=",gshow)
589
590                                      );
591
592        Append(allgeometries,Concatenation([facetgeometry,"\n",edgegeometry,"\n\n"]));
593    od;
594    return allgeometries;
595end;
596
597
598
599#    ########### Here we are. The orbits are colored.
600#    ## Now we mark a vertex for each facet to give the whole thing some
601#    ## sort of orientation.
602#    ## Facet normals are calculated as well.
603#
604#    normals:=[];
605#    markingpoints:=[];
606#    marklines:=[];
607#    currentpointnr:=0;
608#
609#    for facetorbit in facetorbits
610#      do
611#        facet:=facetorbit[1];
612#        facetpos:=Position(facets,facet);
613#        v1:=vertices[ofacets[facetpos][1]];
614#        v2:=vertices[ofacets[facetpos][2]];
615#        v3:=vertices[ofacets[facetpos][3]];
616#        cp:=crossProduct((v2-v1)*cd,(v3-v1)*cd);
617#        normal:=cp;
618#
619#        maps:=List([1..Size(facetorbit)],i->
620#                   RepresentativeActionOnRightOnSets(group,
621#                           Set(facet,int2vec),
622#                           Set(facetorbit[i],int2vec)
623#                           )
624#                   );
625#        markedvertex:=facet[1];
626#        otherpointsformark:=Flat(Filtered(edges,e->markedvertex in e
627#                                    and IsSubset(facet,e)));
628#        otherpointsformark:=Filtered(Set(otherpointsformark),p->p<>markedvertex);
629#        Apply(otherpointsformark,p->1/3*(int2vec(p)-int2vec(markedvertex))+int2vec(markedvertex));
630#
631#
632#        for mapnr in [1..Size(maps)]
633#          do
634#            map:=maps[mapnr];
635#            #normals first:
636#            normalimage:=normal*LinearPartOfAffineMatOnRight(map)^cd;
637##            normalimage:=normalimage{[1..3]}*cd;
638#            normalimage:=normalimage*1/sqrt(normalimage*normalimage);
639#            facetpos:=Position(facets,facetorbit[mapnr]);
640#            normals[facetpos]:=enclosedByTag(vector2floatString(normalimage{[1..3]},precision),"n");
641#            #now the marked corners:
642#            imagepoints:=List(otherpointsformark,p->Concatenation(p,[1])*map);
643#            Apply(imagepoints,p->p{[1..3]}*cd);
644#
645#            for point in imagepoints
646#              do
647#                Add(markingpoints,enclosedByTag(vector2floatString(point,precision),"p"));
648#            od;
649#            Add(marklines,enclosedByTag(Concatenation([String(currentpointnr)," ",String(currentpointnr+1)]),"l"));
650#            currentpointnr:=currentpointnr+2;
651#        od;
652#    od;
653#
654#    normals:=enclosedByTag(Concatenation(normals),"normals");
655#
656#    Add(marklines,Concatenation([enclosedByTag("3","thickness"),"\n"]));
657#    marklines:=enclosedByTag(Concatenation(marklines),"lines");
658#    markingpoints:=enclosedByTag(Concatenation(markingpoints),"points");
659#    Add(markingpoints,'\n');
660#
661#    markblock:=Concatenation(
662#                       ["\n",
663#                        enclosedByTagWithOptions(markingpoints,
664#                                "pointSet",
665#                                "dim=\"3\" point=\"hide\""),
666#                        "\n",
667#                        enclosedByTagWithOptions(marklines,
668#                                "lineSet","line=\"show\"")
669#                                ]);
670#
671
672    #################
673#
674#    facetgeometry:=enclosedByTagWithOptions(Concatenation(["\n",pointblock,"\n",facetblock,"\n"]),"geometry","name=\"points and facets\"");
675#    edgegeometry:=enclosedByTagWithOptions(Concatenation(["\n",pointblock,"\n",edgeblock,"\n"]),"geometry","name=\"points and edges\"");
676#
677##    markgeometry:=enclosedByTagWithOptions(markblock,"geometry","geometry=\"show\" name=\"marks\"");
678#
679#    return Concatenation([facetgeometry,"\n",edgegeometry]);#,"\n",markgeometry]);
680
681#end;
682
683
684
685
686
687
688
689
690#############################################################################
691#############################################################################
692
693
694