1#############################################################################
2##
3##  NFan.gi         NConvex package                 Sebastian Gutsche
4##                                                  Kamal Saleh
5##
6##  Copyright 2019 Mathematics Faculty, Siegen University, Germany
7##
8##  Fans for NConvex package.
9##
10#############################################################################
11
12####################################
13##
14## Reps
15##
16####################################
17
18DeclareRepresentation( "IsExternalFanRep",
19                       IsFan and IsExternalConvexObjectRep,
20                       [ ]
21                      );
22
23DeclareRepresentation( "IsConvexFanRep",
24                       IsExternalFanRep,
25                       [ ]
26                      );
27
28DeclareRepresentation( "IsInternalFanRep",
29                       IsFan and IsInternalConvexObjectRep,
30                       [ ]
31                      );
32
33####################################
34##
35## Types and Families
36##
37####################################
38
39
40BindGlobal( "TheFamilyOfFans",
41        NewFamily( "TheFamilyOfFans" , IsFan ) );
42
43BindGlobal( "TheTypeExternalFan",
44        NewType( TheFamilyOfFans,
45                 IsFan and IsExternalFanRep ) );
46
47BindGlobal( "TheTypeConvexFan",
48        NewType( TheFamilyOfFans,
49                 IsConvexFanRep ) );
50
51BindGlobal( "TheTypeInternalFan",
52        NewType( TheFamilyOfFans,
53                 IsInternalFanRep ) );
54
55####################################
56##
57## Constructors
58##
59####################################
60
61##
62InstallMethod( Fan,
63               " for cones",
64               [ IsList ],
65
66  function( cones )
67    local newgens, i, point, extobj, type;
68
69    if Length( cones ) = 0 then
70
71        Error( " no empty fans allowed." );
72
73    fi;
74
75    newgens := [ ];
76
77    for i in cones do
78
79        if IsCone( i ) then
80
81            if IsBound( i!.input_rays ) then
82
83                Add( newgens, i!.input_rays );
84
85            else
86
87                Add( newgens, RayGenerators( i ) );
88
89            fi;
90
91        elif IsList( i ) then
92
93            Add( newgens, i );
94
95        else
96
97            Error( " wrong cones inserted" );
98
99        fi;
100
101    od;
102
103    point := rec( input_cone_list := newgens );
104
105    ObjectifyWithAttributes(
106        point, TheTypeInternalFan
107        );
108
109    SetAmbientSpaceDimension( point, Length( newgens[ 1 ][ 1 ] ) );
110
111    return point;
112
113end );
114
115##
116InstallMethod( Fan,
117               " for homalg fans",
118               [ IsFan ],
119
120  IdFunc
121
122);
123
124
125##
126InstallMethod( Fan,
127               "for rays and cones.",
128               [ IsList, IsList ],
129
130  function( rays, cones )
131    local point, indices;
132
133    if Length( cones ) = 0 or Length( rays ) = 0 then
134
135        Error( "fan has to have the trivial cone.\n" );
136
137    fi;
138
139    indices := Set( Flat( cones ) );
140
141    if not ForAll( indices, i -> i > 0 and i<= Length( rays ) ) then
142
143      Error( "wrong cones inserted \n" );
144
145    fi;
146
147    point := rec( input_rays := rays, input_cones := cones );
148
149    ObjectifyWithAttributes(
150        point, TheTypeConvexFan
151        );
152
153    SetAmbientSpaceDimension( point, Length( rays[ 1 ] ) );
154
155    return point;
156
157end );
158
159##
160InstallMethod( FanWithFixedRays,
161               "for rays and cones.",
162               [ IsList, IsList ],
163
164  Fan );
165
166##
167InstallMethod( DeriveFansFromTriangulation,
168               "for a list of rays.",
169               [ IsList, IsBool ],
170  function( rays, single_fan_desired )
171    local triangulations, i, j, fans;
172
173    # (0) Check if TopcomInterface is available
174    if TestPackageAvailability( "TopcomInterface", ">=2019.06.15" ) = fail then
175      Error( "The package TopcomInterface is not available " );
176    fi;
177
178    # (1) and marked to be loaded as suggested package
179    if not IsPackageMarkedForLoading( "TopcomInterface", ">=2019.06.15" ) then
180      Error( "The package TopcomInterface has not been marked to be loaded by NConvex if available " );
181    fi;
182
183    # (2) Check that the given rays are valid input
184
185    # (2a) Are the rays all of the same length?
186    if Length( DuplicateFreeList( List( [ 1 .. Length( rays ) ], i -> Length( rays[ i ] ) ) ) ) > 1 then
187      Error( "The rays must be lists of the equal lengths " );
188      return;
189    fi;
190
191    # (2b) Are the rays lists of integers?
192    for i in [ 1 .. Length( rays ) ] do
193      for j in [ 1 .. Length( rays[ i ] ) ] do
194        if not IsInt( rays[ i ][ j ] ) then
195          Error( "The rays must be lists of integers " );
196          return;
197        fi;
198      od;
199    od;
200
201    # (3) compute all fine and regular triangulations
202    triangulations := ValueGlobal( "points2allfinetriangs" )( rays, [], ["regular"] );
203
204    # to match the conventions of NConvex, the counting of rays must start at 1
205    # whilst topcom (in C++-standard) starts the counting at 0
206    Apply( triangulations, i -> i + 1 );
207
208    # (4) iterate over the obtained triangulations and turn them into fans
209    if single_fan_desired then
210      fans := Fan( rays, triangulations[ 1 ] );
211    else
212      fans := List( [ 1 .. Length( triangulations ) ], i -> Fan( rays, triangulations[ i ] ) );
213    fi;
214
215    # return the result
216    return fans;
217
218end );
219
220##
221InstallMethod( FansFromTriangulation,
222               "for a list of rays.",
223               [ IsList ],
224  function( rays )
225
226    return DeriveFansFromTriangulation( rays, false );
227
228end );
229
230##
231InstallMethod( FanFromTriangulation,
232               "for a list of rays.",
233               [ IsList ],
234  function( rays )
235
236    return DeriveFansFromTriangulation( rays, true );
237
238end );
239
240##############################
241##
242##  Attributes
243##
244##############################
245
246
247InstallMethod( RayGenerators,
248               "for fans",
249               [ IsFan ],
250  function( fan )
251
252    return RayGenerators( CanonicalizeFan( fan ) );
253
254end );
255
256InstallMethod( RaysInMaximalCones,
257               [ IsFan ],
258  function( fan )
259
260    return RaysInMaximalCones( CanonicalizeFan( fan ) );
261
262end );
263
264##
265InstallMethod( GivenRayGenerators,
266               "for external fans.",
267               [ IsFan ],
268
269  function( fan )
270
271    if IsBound( fan!.input_rays ) then
272
273        return fan!.input_rays;
274
275    elif IsBound( fan!.input_cone_list ) then
276
277        #return DuplicateFreeList( Concatenation( fan!.input_cone_list ) );
278        return List( Set( Union( fan!.input_cone_list ) ) );
279
280    else
281
282        Error( "Something went wrong." );
283
284    fi;
285
286end );
287
288##
289InstallMethod( Rays,
290               "for fans.",
291               [ IsFan ],
292
293  function( fan )
294    local rays;
295
296    rays := RayGenerators( fan );
297
298    rays := List( rays, i -> Cone( [ i ] ) );
299
300    List( rays, function( i ) SetContainingGrid( i, ContainingGrid( fan ) ); return 0; end );
301
302    return rays;
303
304end );
305
306## This function ignore the small cones which live in some other cone.
307
308InstallMethod( RaysInTheGivenMaximalCones,
309               "for fans",
310               [ IsFan ],
311
312  function( fan )
313    local rays, cones, i, j;
314
315    if IsBound( fan!.input_cones ) and IsBound( fan!.input_rays ) then
316
317        rays := GivenRayGenerators( fan );
318
319        cones := List( [ 1 .. Length( fan!.input_cones ) ], i -> List( [ 1 .. Length( rays ) ], j -> 0 ) );
320
321        for i in [ 1 .. Length( fan!.input_cones ) ] do
322
323          for j in fan!.input_cones[ i ] do
324
325            cones[ i ][ j ] := 1;
326
327          od;
328
329        od;
330
331        return ListOfMaximalConesInList( cones );
332
333    fi;
334
335    if IsBound( fan!.input_cone_list ) then
336
337      rays := GivenRayGenerators( fan );
338
339      ## Dont use ListWithIdenticalEntries here since it has new sideeffects.
340      cones := List( [ 1 .. Length( fan!.input_cone_list ) ], i -> List( [ 1 .. Length( rays ) ], j -> 0 ) );
341
342      for i in [ 1 .. Length( fan!.input_cone_list ) ] do
343
344        for j in [ 1 .. Length( rays ) ] do
345
346          if rays[ j ] in fan!.input_cone_list[ i ] then
347
348            cones[ i ][ j ] := 1;
349
350          fi;
351
352        od;
353
354      od;
355
356      return ListOfMaximalConesInList( cones );
357
358    fi;
359
360    if IsCone( fan ) then
361
362      return RaysInMaximalCones( Fan( [ fan ] ) );
363
364    fi;
365
366    TryNextMethod(  );
367
368end );
369
370InstallMethod( MaximalCones,
371               "for external fans.",
372               [ IsFan ],
373
374  function( fan )
375    local raylist, rays, conelist, i, lis, j;
376
377    raylist := RaysInMaximalCones( fan );
378
379    rays := RayGenerators( fan );
380
381    conelist := [ ];
382
383    for i in [ 1..Length( raylist ) ] do
384
385      lis := [ ];
386
387      for j in [ 1 .. Length( raylist[ i ] ) ] do
388
389        if raylist[ i ][ j ] = 1 then
390
391          lis := Concatenation( lis, [ rays[ j ] ] );
392
393        fi;
394
395      od;
396
397      conelist := Concatenation( conelist, [ lis ] );
398
399    od;
400
401    conelist := List( conelist, Cone );
402
403    Perform( conelist, function( i ) SetContainingGrid( i, ContainingGrid( fan ) ); return 0; end );
404
405    Perform( conelist, function( i ) SetSuperFan( i, fan ); return 0; end );
406
407    return conelist;
408
409end );
410
411##
412InstallMethod( MaximalCones,
413               [ IsFan, IsInt],
414
415  function( fan, n )
416    local all_max_cones, new_list, i;
417
418    all_max_cones:= MaximalCones( fan );
419
420    new_list:= [ ];
421
422    for i in all_max_cones do
423
424      if Dimension( i ) = n then Add( new_list, i ); fi;
425
426    od;
427
428    return new_list;
429
430end );
431
432##
433InstallMethod( GivenMaximalCones,
434               "for external fans.",
435               [ IsFan ],
436
437  function( fan )
438    local raylist, rays, conelist, i, lis, j;
439
440    raylist := RaysInTheGivenMaximalCones( fan );
441
442    rays := GivenRayGenerators( fan );
443
444    conelist := [ ];
445
446    for i in [ 1..Length( raylist ) ] do
447
448      lis := [ ];
449
450      for j in [ 1 .. Length( raylist[ i ] ) ] do
451
452        if raylist[ i ][ j ] = 1 then
453
454          lis := Concatenation( lis, [ rays[ j ] ] );
455
456        fi;
457
458      od;
459
460      conelist := Concatenation( conelist, [ lis ] );
461
462    od;
463
464    conelist := List( conelist, Cone );
465
466    Perform( conelist, function( i ) SetContainingGrid( i, ContainingGrid( fan ) ); return 0; end );
467
468    Perform( conelist, function( i ) SetSuperFan( i, fan ); return 0; end );
469
470    return conelist;
471
472end );
473
474##
475InstallMethod( CanonicalizeFan,
476               [ IsFan ],
477
478  function( fan )
479    local list_of_max, new_gen, cones,i,j, F, max_cones, rays_in_max_cones;
480
481    list_of_max := GivenMaximalCones( fan );
482
483    new_gen := [ ];
484
485    for i in list_of_max do
486
487      Append( new_gen, RayGenerators( i ) );
488
489    od;
490
491    #new_gen := DuplicateFreeList( new_gen );
492
493    new_gen:= List( Set( new_gen ) );
494
495    cones := List( [ 1 .. Length( list_of_max ) ], i -> List( [ 1 .. Length( new_gen ) ], j -> 0 ) );
496
497    for i in [ 1 .. Length( list_of_max ) ] do
498
499      for j in [ 1 .. Length( new_gen ) ] do
500
501        if new_gen[ j ] in RayGenerators( list_of_max[ i ] ) then
502
503          cones[ i ][ j ] := 1;
504
505        fi;
506
507      od;
508
509    od;
510
511    max_cones:= ListOfMaximalConesInList( cones );
512
513    rays_in_max_cones:= [ ];
514
515    for i in [ 1 .. Length( max_cones ) ] do
516
517      Add( rays_in_max_cones, [ ] );
518
519      for j in [ 1..Length( new_gen ) ] do
520
521           if max_cones[ i ][ j ] =1 then Add( rays_in_max_cones[ i ], j ); fi;
522
523      od;
524
525    od;
526
527    F := Fan( new_gen, rays_in_max_cones );
528
529    SetRayGenerators( F, new_gen );
530
531    SetRaysInMaximalCones( F, max_cones );
532
533    return F;
534
535end );
536
537##
538InstallMethod( RaysInAllCones,
539            [ IsFan ],
540
541  function( fan )
542    local max_cones, cones, current_list_of_faces, rays, L, i;
543
544    if IsCone( fan ) then
545
546      max_cones := [ fan ];
547
548    else
549
550      max_cones:= MaximalCones( fan );
551
552    fi;
553
554    cones := [ ];
555
556    for i in max_cones do
557
558      current_list_of_faces:= FacesOfCone( i );
559
560      cones := Concatenation( cones, List( current_list_of_faces, RayGenerators ) );
561
562    od;
563
564    cones := DuplicateFreeList( cones );
565
566    rays := RayGenerators( fan );
567
568    if not ForAll( DuplicateFreeList( Concatenation( cones ) ), r -> r = [ ] or r in rays ) then
569
570      # If this error happens, it means that r is a positive multiple of some element in rays,
571      # which is not problem and can be fixed very easily.
572      Error( "This should not happen, please report this error!" );
573
574    fi;
575
576    L := List( cones, cone ->
577        List( rays, function( r )
578                      if r in cone then
579                        return 1;
580                      else
581                        return 0;
582                      fi;
583                    end ) );
584
585    return DuplicateFreeList( L );
586
587end );
588
589##
590InstallMethod( AllCones,
591        [ IsFan ],
592  function( fan )
593    local n, rays, rays_in_all_cones;
594
595    n := AmbientSpaceDimension( fan );
596
597    rays := RayGenerators( fan );
598
599    rays_in_all_cones := RaysInAllCones( fan );
600
601    rays_in_all_cones := List( rays_in_all_cones, r -> rays{ Positions( r, 1 ) } );
602
603    rays_in_all_cones[ Position( rays_in_all_cones, [  ] ) ] := [ ListWithIdenticalEntries( n, 0 ) ];
604
605    return List( rays_in_all_cones, Cone );
606
607end );
608
609##
610InstallMethod( FVector,
611               [ IsFan ],
612  function( fan )
613    local dim_of_cones;
614
615    dim_of_cones := List( AllCones( fan ), Dimension );
616
617    return List( [ 1.. Dimension( fan ) ], i-> Length( Positions( dim_of_cones, i ) ) );
618
619end );
620
621##
622InstallMethod( Dimension,
623               "for fans",
624               [ IsFan ],
625
626  function( fan )
627
628    return Maximum( List(MaximalCones( fan ), i-> Dimension( i ) ) );
629
630end );
631
632##
633InstallMethod( AmbientSpaceDimension,
634               "for fans",
635               [ IsFan ],
636
637  function( fan )
638
639    return Length( RayGenerators( fan )[ 1 ] );
640
641end );
642
643##
644InstallMethod( PrimitiveCollections,
645               "for fans",
646               [ IsFan ],
647
648  function( fan )
649    local rays, max_cones, facets, all_points, d_max, primitive_collections, d, checked, facet,
650         I_minus_j, scanner, j, I;
651
652    # collect data of the fan
653    rays := RayGenerators( fan );
654
655    max_cones := MaximalCones( fan );
656
657    facets := List( [ 1 .. Length( max_cones ) ],
658                        i -> List( RayGenerators( max_cones[ i ] ), k -> Position( rays, k ) ) );
659
660    all_points := [ 1 .. Length( rays ) ];
661
662    d_max := Maximum( List( facets, i -> Length( i ) ) ) + 1;
663
664    # identify and return the primitive collections
665    primitive_collections := [];
666
667    for d in [ 1 .. d_max ] do
668      checked := [];
669
670      for facet in facets do
671
672        for I_minus_j in Combinations( facet, d ) do
673
674          scanner := Difference( all_points, Flat( I_minus_j ) );
675
676          for j in scanner do
677
678            I := Concatenation( I_minus_j, [ j ] );
679
680            if not I in checked then
681
682              Add( checked, I );
683
684              # (1) I is contained in the primitive collections iff it is not contained in any facet
685              if ForAll( [ 1 .. Length( facets ) ],
686                            i -> not IsSubsetSet( facets[ i ], I ) ) then
687
688                # (2) I is generator of the primitive collections iff
689                # primitive_collections does not contain a "smaller" generator
690                if ForAll( [ 1 .. Length( primitive_collections ) ],
691                              i -> not IsSubsetSet( I, primitive_collections[i] ) ) then
692
693                  # then add this new generator
694                  Append( primitive_collections, [ I ] );
695
696                fi;
697
698              fi;
699
700            fi;
701          od;
702
703        od;
704
705      od;
706
707    od;
708
709    return primitive_collections;
710
711end );
712
713#########################
714##
715##  Properties
716##
717#########################
718
719InstallMethod( IsWellDefinedFan,
720               [ IsFan ],
721
722  function( fan )
723    local max_cones, combi, n;
724
725    max_cones := MaximalCones( fan );
726
727    n := Length( max_cones );
728
729    combi := Combinations( [ 1 .. n ], 2 );
730
731    return ForAll( combi,
732      function( two_indices )
733        local C1, C2, U;
734
735        C1 := max_cones[ two_indices[ 1 ] ];
736
737        C2 := max_cones[ two_indices[ 2 ] ];
738
739        U := IntersectionOfCones( C1, C2 );
740
741        if U in FacesOfCone( C1 ) and U in FacesOfCone( C2 ) then
742
743          return true ;
744
745        else
746
747          return false;
748
749        fi;
750
751      end );
752
753end );
754
755##
756InstallMethod( IsComplete,
757              [ IsFan ],
758  function( fan )
759    local list_of_cones, facets, facets_without_duplicates, positions;
760
761    if Dimension( fan ) < AmbientSpaceDimension( fan ) then
762
763      return false;
764
765    fi;
766
767    list_of_cones := MaximalCones( fan );
768
769    if not ForAll( list_of_cones, IsFullDimensional ) then
770
771      return false;
772
773    fi;
774
775    facets := Concatenation( List( list_of_cones, Facets ) );
776
777    facets_without_duplicates := DuplicateFreeList( facets );
778
779    positions := List( facets_without_duplicates,
780                  facet -> Length( Positions( facets, facet ) ) );
781
782    if Set( positions ) = Set( [ 2 ] ) then
783
784      return true;
785
786    elif Set( positions ) = Set( [ 1, 2 ] ) then
787
788      return false;
789
790    else
791
792      Print( "This should not happen, This may be caused by a not well-defined fan!\n" );
793      return false;
794
795    fi;
796
797end );
798
799##
800InstallMethod( IsPointed,
801               "for fans",
802               [ IsFan ],
803
804  function( fan )
805
806    return ForAll( MaximalCones( fan ), IsPointed );
807
808end );
809
810##
811InstallMethod( IsSmooth,
812               "for fans",
813               [ IsFan ],
814
815  function( fan )
816
817    return ForAll( MaximalCones( fan ), IsSmooth );
818
819end );
820
821##
822InstallMethod( IsRegularFan,
823               "whether a fan is a normalfan or not",
824               [ IsFan ],
825
826  function( fan )
827    local max_cones, ambient_dim, rays, max_cones_ineqs, embed, nr_rays, nd, equations,
828      inequations, r, L1, L0, i, hyper_surface, cone, index_rays;
829
830    if not IsComplete( fan ) then
831
832      return false;
833
834    fi;
835
836    if AmbientSpaceDimension( fan ) <= 2 then
837
838      return true;
839
840    fi;
841
842    ## Algorithm is taken from the Maple Convex package.
843    rays := RayGenerators( fan );
844
845    ambient_dim := AmbientSpaceDimension( fan );
846
847    max_cones := MaximalCones( fan );
848
849    max_cones_ineqs := List( max_cones, DefiningInequalities );
850
851    nr_rays := Length( rays );
852
853    nd := ambient_dim * Length( max_cones );
854
855    embed := function( a, b, c, d, e )
856              local return_list, e1, d1;
857              if e < c then
858                 e1 := e;
859                 e := c;
860                 c := e1;
861                 d1 := d;
862                 d := b;
863                 b := d1;
864              fi;
865              return_list := ListWithIdenticalEntries( c, 0 );
866              return_list := Concatenation( return_list, b );
867              return_list := Concatenation( return_list,
868              ListWithIdenticalEntries( e - Length( b ) - c, 0 ) );
869              return_list := Concatenation( return_list, d );
870              return Concatenation( return_list,
871                ListWithIdenticalEntries( a - Length( return_list ), 0 ) );
872             end;
873
874    ## FIXME: Our convention is to handle only pointed fans.
875    ## convex handles fans with lineality spaces, so the lines differ.
876    equations := List( [ 1 .. Length( max_cones ) ],
877                       i -> List( EqualitiesOfCone( max_cones[ i ] ),
878                                  r -> embed( nd, r, ambient_dim * ( i - 1 ), [ ], 0 ) ) );
879
880    equations := Concatenation( equations );
881
882    inequations := [ ];
883
884    index_rays := [ 1 .. nr_rays ];
885
886    for r in [ 1 .. nr_rays ] do
887
888        L0 := [];
889
890        L1 := [];
891
892        for i in [ 1 .. Length( max_cones ) ] do
893
894          if rays[ r ] in max_cones[ i ] then
895
896            Add( L1, i );
897
898          else
899
900            Add( L0, i );
901
902          fi;
903
904        od;
905
906        i := ambient_dim * ( L1[ 1 ] - 1 );
907
908        index_rays[ r ] := i;
909
910        Remove( L1, L1[ 1 ] );
911
912        equations := Concatenation( equations,
913                          List( L1,
914                            j -> embed( nd, rays[ r ], i, - rays[ r ], ambient_dim * ( j - 1 ) ) ) );
915
916        inequations := Concatenation( inequations,
917                          List( L0,
918                            j -> embed( nd, rays[ r ], i, - rays[ r ], ambient_dim * ( j - 1 ) ) ) );
919
920    od;
921
922    hyper_surface := ConeByEqualitiesAndInequalities( equations, [ ] );
923
924    i := AmbientSpaceDimension( hyper_surface ) - Dimension( hyper_surface );
925
926    cone := ConeByEqualitiesAndInequalities( equations, inequations );
927
928    r := AmbientSpaceDimension( cone ) - Dimension( cone );
929
930    return i = r;
931
932end );
933
934##
935## This is implementation of the Shephard's criterion
936## (Theorem 4.7, Combinatorial convexity and algebraic geometry, Ewald, Guenter)
937##
938InstallMethod( IsNormalFan,
939          [ IsFan ],
940
941  function( fan )
942      local rays, cones, mat, G, polytopes, P, M;
943
944      if not IsComplete( fan ) then
945
946        return false;
947
948      fi;
949
950      if AmbientSpaceDimension( fan ) <= 2 then
951
952        return true;
953
954      fi;
955
956      if HasIsRegularFan( fan ) then
957
958        return IsRegularFan( fan );
959
960      fi;
961
962      rays := RayGenerators( fan );
963
964      cones := ShallowCopy( RaysInAllCones( fan ) );
965
966      Remove( cones, PositionProperty( cones, IsZero ) );
967
968      mat := HomalgMatrix( rays, HOMALG_RATIONALS );
969
970      G := GaleTransform( mat );
971
972      if IsZero( G ) then
973
974        return true;
975
976      fi;
977
978      G := EntriesOfHomalgMatrixAsListList( G );
979
980      polytopes := List( cones, cone -> Set( DuplicateFreeList( G{ Positions( cone, 0 ) } ) ) );
981
982      polytopes := DuplicateFreeList( polytopes );
983
984      polytopes := List( polytopes, l -> Polytope( l ) );
985
986      P := Iterated( polytopes, IntersectionOfPolytopes );
987
988      if Dimension( P ) = -1 then
989
990        return false;
991
992      elif Dimension( P ) = 0 then
993
994        M := Vertices( P )[ 1 ];
995
996        return ForAll( polytopes, P -> IsInteriorPoint( M, P ) );
997
998      else
999
1000        M := RandomInteriorPoint( P );
1001
1002        if ForAll( polytopes, P -> IsInteriorPoint( M, P ) ) then
1003
1004          return true;
1005
1006        else
1007
1008          return false;
1009
1010        fi;
1011
1012      fi;
1013
1014end );
1015
1016##
1017InstallMethod( IsRegularFan, [ IsFan ], IsNormalFan );
1018
1019##
1020InstallMethod( IsFullDimensional,
1021               "for fans",
1022               [ IsFan ],
1023
1024  function( fan )
1025
1026    return ForAny( MaximalCones( fan ), i -> Dimension( i ) = AmbientSpaceDimension( i ) );
1027
1028end );
1029
1030
1031##
1032InstallMethod( IsSimplicial,
1033               " for homalg fans",
1034               [ IsFan ],
1035
1036  function( fan )
1037
1038    fan := MaximalCones( fan );
1039
1040    return ForAll( fan, IsSimplicial );
1041
1042end );
1043
1044##
1045InstallMethod( IsFanoFan,
1046            [ IsFan ],
1047  function( fan )
1048    local polyt;
1049
1050    if HasIsComplete( fan ) and not IsComplete( fan ) then
1051
1052      return false;
1053
1054    fi;
1055
1056    polyt := Polytope( RayGenerators( fan ) );
1057
1058      if not IsFullDimensional( polyt ) or
1059            not IsInteriorPoint(
1060              ListWithIdenticalEntries( AmbientSpaceDimension( polyt ), 0 ), polyt ) then
1061
1062        return false;
1063
1064      fi;
1065
1066    return fan = NormalFan( PolarPolytope( polyt ) );
1067
1068end );
1069
1070#########################
1071##
1072##  Methods
1073##
1074#########################
1075
1076##
1077InstallMethod( \=,
1078            [ IsFan, IsFan ],
1079  function( fan1, fan2 )
1080
1081    if RayGenerators( fan1 ) = RayGenerators( fan2 ) and
1082        Set( RaysInMaximalCones( fan1 ) ) = Set( RaysInMaximalCones( fan2 ) ) then
1083
1084          return true;
1085
1086    fi;
1087
1088    if RayGenerators( fan1 ) <> RayGenerators( fan2 ) and
1089        Set( RayGenerators( fan1 ) ) = Set( RayGenerators( fan2 ) ) then
1090
1091          Error( "This should not happen! Please report this error." );
1092
1093    fi;
1094
1095    return false;
1096
1097end );
1098
1099##
1100InstallMethod( \*,
1101               "for fans.",
1102               [ IsFan, IsFan ],
1103
1104  function( fan1, fan2 )
1105    local rays1, rays2, m1, m2, new_m, new_rays, cones1, cones2, i, j, k, new_cones, akt_cone, new_fan;
1106
1107    rays1 := RayGenerators( fan1 );
1108
1109    rays2 := RayGenerators( fan2 );
1110
1111    m1 := Rank( ContainingGrid( fan1 ) );
1112
1113    m2 := Rank( ContainingGrid( fan2 ) );
1114
1115    m1 := List( [ 1 .. m1 ], i -> 0 );
1116
1117    m2 := List( [ 1 .. m2 ], i -> 0 );
1118
1119    rays1 := List( rays1, i -> Concatenation( i, m2 ) );
1120
1121    rays2 := List( rays2, i -> Concatenation( m1, i ) );
1122
1123    new_rays := Concatenation( rays1, rays2 );
1124
1125    cones1 := RaysInMaximalCones( fan1 );
1126
1127    cones2 := RaysInMaximalCones( fan2 );
1128
1129    new_cones := [ ];
1130
1131    m1 := Length( rays1 );
1132
1133    m2 := Length( rays2 );
1134
1135    for i in cones1 do
1136
1137      for j in cones2 do
1138
1139        akt_cone := [ ];
1140
1141        for k in [ 1 .. m1 ] do
1142
1143          if i[ k ] = 1 then
1144
1145            Add( akt_cone, k );
1146
1147          fi;
1148
1149        od;
1150
1151        for k in [ 1 .. m2 ] do
1152
1153          if j[ k ] = 1 then
1154
1155            Add( akt_cone, k + m1 );
1156
1157          fi;
1158
1159        od;
1160
1161        Add( new_cones, akt_cone );
1162
1163      od;
1164
1165    od;
1166
1167    new_fan := FanWithFixedRays( new_rays, new_cones );
1168
1169    SetContainingGrid( new_fan, ContainingGrid( fan1 ) + ContainingGrid( fan2 ) );
1170
1171    return new_fan;
1172
1173end );
1174
1175##
1176InstallMethod( ToricStarFan,
1177               "for fans",
1178               [ IsFan, IsCone ],
1179
1180  function( fan, cone )
1181    local maximal_cones, rays_of_cone, defining_inequalities, value_list, cone_list, i, j, breaker;
1182
1183    maximal_cones := MaximalCones( fan );
1184
1185    rays_of_cone := RayGenerators( cone );
1186
1187    cone_list := [ ];
1188
1189    breaker := false;
1190
1191    for i in maximal_cones do
1192
1193      defining_inequalities := DefiningInequalities( i );
1194
1195      for j in rays_of_cone do
1196
1197        value_list := List( defining_inequalities, k -> k * j );
1198
1199        if not ForAll( value_list, k -> k >= 0 ) or not 0 in value_list then
1200
1201          breaker := true;
1202
1203          continue;
1204
1205        fi;
1206
1207      od;
1208
1209      if breaker then
1210
1211        breaker := false;
1212
1213        continue;
1214
1215      fi;
1216
1217      Add( cone_list, cone );
1218
1219    od;
1220
1221    cone_list := Fan( cone_list );
1222
1223    SetContainingGrid( cone_list, ContainingGrid( fan ) );
1224
1225end );
1226##
1227InstallMethod( \*,
1228               "for homalg fans.",
1229               [ IsCone, IsFan ],
1230
1231  function( cone, fan )
1232
1233    return Fan( [ cone ] ) * fan;
1234
1235end );
1236
1237##
1238InstallMethod( \*,
1239               "for homalg fans.",
1240               [ IsFan, IsCone ],
1241
1242  function( fan, cone )
1243
1244    return fan * Fan( [ cone ] );
1245
1246end );
1247
1248##
1249InstallMethod( ToricStarFan,
1250               "for fans",
1251               [ IsFan, IsCone ],
1252
1253  function( fan, cone )
1254    local maximal_cones, rays_of_cone, defining_inequalities, value_list, cone_list, i, j, breaker;
1255
1256    maximal_cones := MaximalCones( fan );
1257
1258    rays_of_cone := RayGenerators( cone );
1259
1260    cone_list := [ ];
1261
1262    breaker := false;
1263
1264    for i in maximal_cones do
1265
1266      defining_inequalities := DefiningInequalities( i );
1267
1268      for j in rays_of_cone do
1269
1270        value_list := List( defining_inequalities, k -> k * j );
1271
1272        if not ForAll( value_list, k -> k >= 0 ) or not 0 in value_list then
1273
1274          breaker := true;
1275
1276          continue;
1277
1278        fi;
1279
1280      od;
1281
1282      if breaker then
1283
1284        breaker := false;
1285
1286        continue;
1287
1288      fi;
1289
1290      Add( cone_list, cone );
1291
1292    od;
1293
1294    cone_list := Fan( cone_list );
1295
1296    SetContainingGrid( cone_list, ContainingGrid( fan ) );
1297
1298    return cone_list;
1299
1300end );
1301
1302#########################
1303##
1304##  Simple functions
1305##
1306#########################
1307
1308InstallMethod( FirstLessTheSecond,
1309               [ IsList, IsList],
1310
1311  function( u, v )
1312
1313    if Length( u ) <> Length( v ) then
1314
1315      Error( "The two lists should have the same length!" );
1316
1317    fi;
1318
1319    return ForAll( [ 1 .. Length( u ) ], i-> u[ i ] <= v[ i ] );
1320
1321end );
1322
1323InstallMethod( OneMaximalConeInList,
1324              [ IsList ],
1325  function( u )
1326    local list, max, new_u, i;
1327
1328    #new_u:= DuplicateFreeList( ShallowCopy( u ) );
1329    new_u:= List( Set( u ) );
1330
1331    max := new_u[ 1 ];
1332
1333    for i in new_u do
1334
1335      if FirstLessTheSecond( max, i ) then
1336
1337        max := i;
1338
1339      fi;
1340
1341    od;
1342
1343    list := [ ];
1344
1345    for i in new_u do
1346
1347      if FirstLessTheSecond( i, max ) then
1348
1349        Add( list, i );
1350
1351      fi;
1352
1353    od;
1354
1355    return [ max, list ];
1356
1357 end );
1358
1359##
1360InstallMethod( ListOfMaximalConesInList,
1361               [ IsList ],
1362
1363  function( L )
1364    local l, list_of_max, current, new_L;
1365
1366    list_of_max := [ ];
1367
1368    #new_L:= DuplicateFreeList( L );
1369    new_L := Set( L );
1370
1371    while Length( new_L )<> 0 do
1372
1373      current := OneMaximalConeInList( new_L );
1374
1375      Add( list_of_max, current[ 1 ] );
1376
1377      SubtractSet( new_L, current[ 2] );
1378
1379    od;
1380
1381    return list_of_max;
1382
1383end );
1384
1385####################################
1386##
1387## Display Methods
1388##
1389####################################
1390
1391##
1392InstallMethod( ViewObj,
1393               "for homalg fans",
1394               [ IsFan ],
1395
1396  function( fan )
1397    local str;
1398
1399    Print( "<A" );
1400
1401    if HasIsComplete( fan ) then
1402
1403      if IsComplete( fan ) then
1404
1405          Print( " complete" );
1406
1407      fi;
1408
1409    fi;
1410
1411    if HasIsPointed( fan ) then
1412
1413      if IsPointed( fan ) then
1414
1415          Print( " pointed" );
1416
1417      fi;
1418
1419    fi;
1420
1421    if HasIsSmooth( fan ) then
1422
1423      if IsSmooth( fan ) then
1424
1425          Print( " smooth" );
1426
1427      fi;
1428
1429    fi;
1430
1431    Print( " fan in |R^" );
1432
1433    Print( String( AmbientSpaceDimension( fan ) ) );
1434
1435    if HasRays( fan ) then
1436
1437      Print( " with ", String( Length( Rays( fan ) ) )," rays" );
1438
1439    fi;
1440
1441    Print( ">" );
1442
1443end );
1444
1445##
1446InstallMethod( Display,
1447               "for homalg polytopes",
1448               [ IsFan ],
1449
1450  function( fan )
1451    local str;
1452
1453    Print( "A" );
1454
1455    if HasIsComplete( fan ) then
1456
1457      if IsComplete( fan ) then
1458
1459        Print( " complete" );
1460
1461      fi;
1462
1463    fi;
1464
1465    Print( " fan in |R^" );
1466
1467    Print( String( AmbientSpaceDimension( fan ) ) );
1468
1469    if HasRays( fan ) then
1470
1471      Print( " with ", String( Length( Rays( fan ) ) )," rays" );
1472
1473    fi;
1474
1475    Print( ".\n" );
1476
1477end );
1478
1479