1############################################################################# 2## 3## NCone.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## 15## Reps 16## 17#################################### 18 19DeclareRepresentation( "IsExternalConeRep", 20 IsCone and IsExternalFanRep, 21 [ ] ); 22 23DeclareRepresentation( "IsConvexConeRep", 24 IsExternalConeRep, 25 [ ] ); 26 27DeclareRepresentation( "IsInternalConeRep", 28 IsCone and IsInternalFanRep, 29 [ ] ); 30 31#################################### 32## 33## Types and Families 34## 35#################################### 36 37 38BindGlobal( "TheFamilyOfCones", 39 NewFamily( "TheFamilyOfCones" , IsCone ) ); 40 41BindGlobal( "TheTypeExternalCone", 42 NewType( TheFamilyOfCones, 43 IsCone and IsExternalConeRep ) ); 44 45BindGlobal( "TheTypeConvexCone", 46 NewType( TheFamilyOfCones, 47 IsConvexConeRep ) ); 48 49BindGlobal( "TheTypeInternalCone", 50 NewType( TheFamilyOfCones, 51 IsInternalConeRep ) ); 52 53 54##################################### 55## 56## Property Computation 57## 58##################################### 59 60## 61InstallMethod( IsPointed, 62 "for homalg cones.", 63 [ IsCone ], 64 65 function( cone ) 66 67 return Cdd_IsPointed( ExternalCddCone ( cone ) ); 68 69end ); 70 71## 72InstallMethod( InteriorPoint, 73 [ IsConvexObject and IsCone ], 74 function( cone ) 75 local point, denominators; 76 point := Cdd_InteriorPoint( ExternalCddCone( cone ) ); 77 denominators := List( point, DenominatorRat ); 78 if DuplicateFreeList( denominators ) = [ 1 ] then 79 return point; 80 else 81 return Lcm( denominators )*point; 82 fi; 83end ); 84 85## 86InstallMethod( IsComplete, 87 " for cones", 88 [ IsCone ], 89 90 function( cone ) 91 local rays; 92 93 if IsPointed( cone ) or not IsFullDimensional( cone ) then 94 95 return false; 96 97 fi; 98 99 rays:= RayGenerators( cone ); 100 101 return IsFullDimensional( cone ) and ForAll( rays, i-> -i in rays ); 102 103end ); 104 105## Let N= Z^{1 \times n}. Then N is free Z-modue. Let r_1,...,r_k \in N be the generating rays of the 106## cone. Let A= [ r_1,..., r_k ]^T \in Z^{ k \times n }. Let M= N/ Z^{1 \times k}A. Now let B the smith 107## normal form of A. Then B \in Z^{ k \times n } and there exists l<=k, l<=n with B_{i,i} \neq 0 and B_{i-1,i}| B_{i-1,i} for 108## all 1<i<= l and B_{i,i}=0 for i>l. We have now M= Z/Zb_{1,1} \oplus ... \oplus Z/Zb_{l,l} \oplus Z^{1 \times max{n,k}-l }. 109## If all B_{i,i}=1 for i<=l, then M= Z^{1 \times max{n,k}-l }. i.e. M is free. Thus there is H \subset N with N= H \oplus Z^{1 \times k}A. 110## ( Corollary 7.55, Advanced modern algebra, J.rotman ). 111## 112InstallMethod( IsSmooth, 113 "for cones", 114 [ IsCone ], 115 116 function( cone ) 117 local rays, smith; 118 119 rays := RayGenerators( cone ); 120 121 if RankMat( rays ) <> Length( rays ) then 122 123 return false; 124 125 fi; 126 127 smith := SmithNormalFormIntegerMat(rays); 128 129 smith := List( Flat( smith ), i -> AbsInt( i ) ); 130 131 if Maximum( smith ) <> 1 then 132 133 return false; 134 135 fi; 136 137 return true; 138 139end ); 140 141## 142InstallMethod( IsRegularCone, 143 "for cones.", 144 [ IsCone ], 145 function( cone ) 146 147 return IsSmooth( cone ); 148 149end ); 150 151 152## 153InstallMethod( IsSimplicial, 154 " for cones", 155 [ IsCone ], 156 function( cone ) 157 local rays; 158 159 rays:= RayGenerators( cone ); 160 161 return Length( rays )= RankMat( rays ); 162 163end ); 164 165## 166InstallMethod( IsFullDimensional, 167 "for cones", 168 [ IsCone ], 169 function( cone ) 170 171 return RankMat( RayGenerators( cone ) ) = AmbientSpaceDimension( cone ); 172 173end ); 174 175## 176InstallTrueMethod( HasConvexSupport, IsCone ); 177 178## 179InstallMethod( IsRay, 180 "for cones.", 181 [ IsCone ], 182 183 function( cone ) 184 185 return Dimension( cone ) = 1; 186 187end ); 188 189##################################### 190## 191## Attribute Computation 192## 193##################################### 194 195InstallMethod( RayGenerators, 196 [ IsCone ], 197 198 function( cone ) 199# local nmz_cone, l, r; 200 201 return Cdd_GeneratingRays( ExternalCddCone( cone ) ); 202 203# nmz_cone := ExternalNmzCone( cone ); 204 205# r := NmzGenerators( nmz_cone ); 206 207# l := NmzMaximalSubspace( nmz_cone ); 208 209# return Concatenation( r, l, -l ); 210 211 end ); 212 213## 214InstallMethod( DualCone, 215 "for external cones", 216 [ IsCone ], 217 218 function( cone ) 219 local dual; 220 221 if RayGenerators( cone ) = [ ] then 222 dual := ConeByInequalities( [ List( [ 1 .. AmbientSpaceDimension( cone ) ], i -> 0 ) ] ); 223 else 224 dual := ConeByInequalities( RayGenerators( cone ) ); 225 fi; 226 227 SetDualCone( dual, cone ); 228 229 SetContainingGrid( dual, ContainingGrid( cone ) ); 230 231 return dual; 232 233end ); 234 235## 236InstallMethod( DefiningInequalities, 237 [ IsCone ], 238 239 function( cone ) 240 local inequalities, new_inequalities, equalities, i, u; 241 242 inequalities:= ShallowCopy( Cdd_Inequalities( ExternalCddCone( cone ) ) ); 243 244 equalities:= ShallowCopy( Cdd_Equalities( ExternalCddCone( cone ) ) ); 245 246 for i in equalities do 247 248 Append( inequalities, [ i,-i ] ); 249 250 od; 251 252new_inequalities:= [ ]; 253 254 for i in inequalities do 255 256 u:= ShallowCopy( i ); 257 258 Remove( u , 1 ); 259 260 Add(new_inequalities, u ); 261 262 od; 263 264 return new_inequalities; 265 266end ); 267 268 269## 270InstallMethod( FactorConeEmbedding, 271 "for cones", 272 [ IsCone ], 273 274 function( cone ) 275 276 return TheIdentityMorphism( ContainingGrid( cone ) ); 277 278end ); 279 280## 281InstallMethod( EqualitiesOfCone, 282 "for external Cone", 283 [ IsCone ], 284 285 function( cone ) 286 local equalities, new_equalities, u, i; 287 288 equalities:= Cdd_Equalities( ExternalCddCone( cone ) ); 289 290 new_equalities:= [ ]; 291 292 for i in equalities do 293 294 u:= ShallowCopy( i ); 295 296 Remove( u , 1 ); 297 298 Add(new_equalities, u ); 299 300 od; 301 302 return new_equalities; 303 304end ); 305 306## 307InstallMethod( RelativeInteriorRay, 308 "for a cone", 309 [ IsCone ], 310 311 function( cone ) 312 local rays, rand_mat; 313 314 rays := RayGenerators( cone ); 315 316 rand_mat := RandomMat( Length( rays ), 1 ); 317 318 rand_mat := Concatenation( rand_mat ); 319 320 rand_mat := List( rand_mat, i -> AbsInt( i ) + 1 ); 321 322 return Sum( [ 1 .. Length( rays ) ], i -> rays[ i ] * rand_mat[ i ] ); 323 324end ); 325 326## 327InstallMethod( HilbertBasisOfDualCone, 328 "for cone", 329 [ IsCone ], 330 331 function( cone ) 332 333 return HilbertBasis( DualCone( cone ) ); 334 335end ); 336 337## This method is commented since it returns a wrong answer for not-pinted cones 338## C := Cone( [ e1 ] ); 339## 340# InstallMethod( HilbertBasisOfDualCone, 341# "for cone", 342# [ IsCone ], 343 344# function( cone ) 345# local ray_generators, d, i, dim, V, D, max, v, I, b, DpD, d1, d2, Dgens, 346# zero_element, entry; 347 348# ray_generators := RayGenerators( cone ); 349 350# dim := AmbientSpaceDimension( cone ); 351 352# max := Maximum( List( Concatenation( ray_generators ), AbsInt ) ); 353 354# D := [ ]; 355 356# ## This needs to be done smarter 357# I:= Cartesian( List( [ 1 .. dim ], i -> [ -max .. max ] ) ); 358 359# for v in I do 360 361# if ForAll( ray_generators, i -> i * v >= 0 ) then 362 363# Add( D, v ); 364 365# fi; 366 367# od; 368 369# DpD := [ ]; 370 371# for d1 in D do 372 373# if d1 * d1 <> 0 then 374 375# for d2 in D do 376 377# if d2 * d2 <> 0 then 378 379# Add( DpD, d1 + d2 ); 380 381# fi; 382 383# od; 384 385# fi; 386 387# od; 388 389# Dgens := [ ]; 390 391# for d in D do 392 393# if not d in DpD then 394 395# Add( Dgens, d ); 396 397# fi; 398 399# od; 400 401# if not Dgens = [ ] then 402 403# zero_element := ListWithIdenticalEntries( Length( Dgens[ 1 ] ), 0 ); 404 405# i := Position( Dgens, zero_element ); 406 407# if i <> fail then 408 409# Remove( Dgens, i ); 410 411# fi; 412 413# fi; 414 415# entry := ToDoListEntry( [ [ cone, "DualCone" ] ], [ DualCone, cone ], "HilbertBasis", Dgens ); 416 417# AddToToDoList( entry ); 418 419# return Dgens; 420 421# end); 422 423## 424InstallMethod( AmbientSpaceDimension, 425 "for cones", 426 [ IsCone ], 427 428 function( cone ) 429 430 if Length( RayGenerators( cone ) ) > 0 then 431 432 return Length( RayGenerators( cone )[ 1 ] ); 433 434 else 435 436 return 1; 437 438 fi; 439 440end ); 441 442## 443InstallMethod( Dimension, 444 "for cones", 445 [ IsCone and HasRayGenerators ], 446 447 function( cone ) 448 449 if Length( RayGenerators( cone ) ) > 0 then 450 451 return RankMat( RayGenerators( cone ) ); 452 453 else 454 455 return 0; 456 457 fi; 458 459 TryNextMethod(); 460 461end ); 462 463## 464InstallMethod( Dimension, 465 "for cones", 466 [ IsCone ], 467 function( cone ) 468 469 return Cdd_Dimension( ExternalCddCone( cone ) ); 470 471end ); 472 473## 474InstallMethod( HilbertBasis, 475 "for cones", 476 [ IsCone ], 477 478 function( cone ) 479 local ineq, const; 480 481 if not IsPointed( cone ) then 482 483 Error( "Hilbert basis for not-pointed cones is not yet implemented, you can use the command 'LatticePointsGenerators' " ); 484 485 fi; 486 487 488 if IsPackageMarkedForLoading( "NormalizInterface", ">=1.1.0" ) then 489 490 return Set( ValueGlobal( "NmzHilbertBasis" )( ExternalNmzCone( cone ) ) ); 491 492 elif IsPackageMarkedForLoading( "4ti2Interface", ">=2018.07.06" ) then 493 494 ineq := DefiningInequalities( cone ); 495 496 const := ListWithIdenticalEntries( Length( ineq ), 0 ); 497 498 return Set( ValueGlobal( "4ti2Interface_zsolve_equalities_and_inequalities" )( [ ], [ ], ineq, const )[ 2 ]: precision := "gmp" ); 499 500 else 501 502 Error( "4ti2Interface or NormalizInterface should be loaded!" ); 503 504 fi; 505 506end ); 507 508## 509InstallMethod( RaysInFacets, 510 " for cones", 511 [ IsCone ], 512 513 function( cone ) 514 local external_cone, list_of_facets, generating_rays, list, current_cone, current_list, current_ray_generators, i; 515 516 external_cone := Cdd_H_Rep ( ExternalCddCone ( cone ) ); 517 518 list_of_facets:= Cdd_Facets( external_cone ); 519 520 generating_rays:= RayGenerators( cone ); 521 522 list:= [ ]; 523 524 for i in list_of_facets do 525 526 current_cone := Cdd_ExtendLinearity( external_cone, i ); 527 528 current_ray_generators := Cdd_GeneratingRays( current_cone ) ; 529 530 current_list:= List( [1..Length( generating_rays )], 531 532 function(j) 533 534 if generating_rays[j] in Cone( current_cone ) then 535 return 1; 536 else 537 return 0; 538 fi; 539 540 end ); 541 542 Add( list, current_list ); 543 544 od; 545 546return list; 547 548end ); 549 550## 551InstallMethod( RaysInFaces, 552 " for cones", 553 [ IsCone ], 554 555 function( cone ) 556 local external_cone, list_of_faces, generating_rays, list, current_cone, current_list, current_ray_generators, i,j; 557 558 external_cone := Cdd_H_Rep( ExternalCddCone ( cone ) ); 559 560 list_of_faces:= Cdd_Faces( external_cone ); 561 562 generating_rays:= RayGenerators( cone ); 563 564 list:= [ ]; 565 566 for i in list_of_faces do 567 568 if i[ 1 ] <> 0 then 569 570 current_cone := Cdd_ExtendLinearity( external_cone, i[ 2 ] ); 571 572 current_ray_generators := Cdd_GeneratingRays( current_cone ) ; 573 574 current_list:= List( [ 1 .. Length( generating_rays ) ], 575 576 function(j) 577 578 if generating_rays[j] in Cone( current_cone ) then 579 return 1; 580 else 581 return 0; 582 fi; 583 584 end ); 585 586 Add( list, current_list ); 587 588 fi; 589 590 od; 591 592return list; 593 594end ); 595 596## 597InstallMethod( Facets, 598 "for external cones", 599 [ IsCone ], 600 601 function( cone ) 602 local raylist, rays, conelist, i, lis, j; 603 604 raylist := RaysInFacets( cone ); 605 606 rays := RayGenerators( cone ); 607 608 conelist := [ ]; 609 610 for i in [ 1..Length( raylist ) ] do 611 612 lis := [ ]; 613 614 for j in [ 1 .. Length( raylist[ i ] ) ] do 615 616 if raylist[ i ][ j ] = 1 then 617 618 lis := Concatenation( lis, [ rays[ j ] ] ); 619 620 fi; 621 622 od; 623 624 conelist := Concatenation( conelist, [ lis ] ); 625 626 od; 627 628 if conelist = [ [ ] ] then 629 return [ Cone( [ List( [ 1 .. AmbientSpaceDimension( cone ) ], i->0 ) ] ) ]; 630 fi; 631 632 return List( conelist, Cone ); 633 634end ); 635 636## 637InstallMethod( FacesOfCone, 638 " for external cones", 639 [ IsCone ], 640 641 function( cone ) 642 local raylist, rays, conelist, i, lis, j; 643 644 raylist := RaysInFaces( cone ); 645 646 rays := RayGenerators( cone ); 647 648 conelist := [ ]; 649 650 for i in [ 1..Length( raylist ) ] do 651 652 lis := [ ]; 653 654 for j in [ 1 .. Length( raylist[ i ] ) ] do 655 656 if raylist[ i ][ j ] = 1 then 657 658 lis := Concatenation( lis, [ rays[ j ] ] ); 659 660 fi; 661 662 od; 663 664 conelist := Concatenation( conelist, [ lis ] ); 665 666 od; 667 668 lis := [ Cone( [ List( [ 1 .. AmbientSpaceDimension( cone ) ], i-> 0 ) ] ) ]; 669 670 return Concatenation( lis, List( conelist, Cone ) ); 671 672end ); 673 674## 675InstallMethod( FVector, 676 "for cones", 677 [ IsCone ], 678 function( cone ) 679 local external_cone, faces; 680 681 external_cone := Cdd_H_Rep( ExternalCddCone( cone ) ); 682 683 faces := Cdd_Faces( external_cone ); 684 685 return List( [ 1 .. Dimension( cone ) ], 686 i -> Length( PositionsProperty( faces, face -> face[ 1 ] = i ) ) ); 687 end ); 688 689## 690InstallMethod( LinearSubspaceGenerators, 691 [ IsCone ], 692 693 function( cone ) 694 local inequalities, basis, new_basis; 695 696 inequalities:= DefiningInequalities( cone ); 697 698 basis := NullspaceMat( TransposedMat( inequalities ) ); 699 700 new_basis := List( basis, i-> LcmOfDenominatorRatInList( i )*i ); 701 702 return new_basis; 703 704end ); 705 706## 707InstallMethod( LinealitySpaceGenerators, 708 [ IsCone ], 709 710 LinearSubspaceGenerators 711 712 ); 713 714## 715InstallMethod( GridGeneratedByCone, 716 " for homalg cones.", 717 [ IsCone ], 718 719 function( cone ) 720 local rays, M, grid; 721 722 rays := RayGenerators( cone ); 723 724 M := HomalgMatrix( rays, HOMALG_MATRICES.ZZ ); 725 726 M := HomalgMap( M, "free", ContainingGrid( cone ) ); 727 728 grid := ImageSubobject( M ); 729 730 SetFactorGrid( cone, FactorObject( grid ) ); 731 732 SetFactorGridMorphism( cone, CokernelEpi( M ) ); 733 734 return grid; 735 736end ); 737 738## 739InstallMethod( GridGeneratedByOrthogonalCone, 740 " for homalg cones.", 741 [ IsCone ], 742 743 function( cone ) 744 local rays, M; 745 746 rays := RayGenerators( cone ); 747 748 M := HomalgMatrix( rays, HOMALG_MATRICES.ZZ ); 749 750 M := Involution( M ); 751 752 M := HomalgMap( M, ContainingGrid( cone ), "free" ); 753 754 return KernelSubobject( M ); 755 756end ); 757 758## 759InstallMethod( FactorGrid, 760 " for homalg cones.", 761 [ IsCone ], 762 763 function( cone ) 764 local rays, M, grid; 765 766 rays := RayGenerators( cone ); 767 768 M := HomalgMatrix( rays, HOMALG_MATRICES.ZZ ); 769 770 M := HomalgMap( M, "free", ContainingGrid( cone ) ); 771 772 grid := ImageSubobject( M ); 773 774 SetGridGeneratedByCone( cone, grid ); 775 776 SetFactorGridMorphism( cone, CokernelEpi( M ) ); 777 778 return FactorObject( grid ); 779 780end ); 781 782## 783InstallMethod( FactorGridMorphism, 784 " for homalg cones.", 785 [ IsCone ], 786 787 function( cone ) 788 local rays, grid, M; 789 790 rays := RayGenerators( cone ); 791 792 M := HomalgMatrix( rays, HOMALG_MATRICES.ZZ ); 793 794 M := HomalgMap( M, "free", ContainingGrid( cone ) ); 795 796 grid := ImageSubobject( M ); 797 798 SetGridGeneratedByCone( cone, grid ); 799 800 SetFactorGrid( cone, FactorObject( grid ) ); 801 802 return CokernelEpi( M ); 803 804end ); 805 806InstallMethod( LatticePointsGenerators, 807 [ IsCone ], 808 809 function( cone ) 810 local n; 811 812 n:= AmbientSpaceDimension( cone ); 813 814 return LatticePointsGenerators( Polyhedron( [ List( [ 1 .. n ], i -> 0 ) ], cone ) ); 815 816end ); 817 818## 819InstallMethod( StarSubdivisionOfIthMaximalCone, 820 " for homalg cones and fans", 821 [ IsFan, IsInt ], 822 823 function( fan, noofcone ) 824 local maxcones, cone, ray, cone2; 825 826 maxcones := MaximalCones( fan ); 827 828 if Length( maxcones ) < noofcone then 829 830 Error( " not enough maximal cones" ); 831 832 fi; 833 834 if not IsSmooth( maxcones[ noofcone ] ) then 835 836 Error( " the specified cone is not smooth!" ); 837 838 fi; 839 840 maxcones := List( maxcones, RayGenerators ); 841 842 cone := maxcones[ noofcone ]; 843 844 ray := Sum( cone ); 845 846 cone2 := Concatenation( cone, [ ray ] ); 847 848 cone2 := UnorderedTuples( cone2, Length( cone2 ) - 1 ); 849 850 Apply( cone2, Set ); 851 852 Apply( maxcones, Set ); 853 854 maxcones := Concatenation( maxcones, cone2 ); 855 856 maxcones := Difference( maxcones, [ Set( cone ) ] ); 857 858 maxcones := Fan( maxcones ); 859 860 SetContainingGrid( maxcones, ContainingGrid( fan ) ); 861 862 return maxcones; 863 864end ); 865 866 867## 868InstallMethod( StarFan, 869 " for homalg cones in fans", 870 [ IsCone and HasSuperFan ], 871 872 function( cone ) 873 874 return StarFan( cone, SuperFan( cone ) ); 875 876end ); 877 878## 879InstallMethod( StarFan, 880 " for homalg cones", 881 [ IsCone, IsFan ], 882 883 function( cone, fan ) 884 local maxcones; 885 886 maxcones := MaximalCones( fan ); 887 888 maxcones := Filtered( maxcones, i -> Contains( i, cone ) ); 889 890 maxcones := List( maxcones, HilbertBasis ); 891 892 ## FIXME: THIS IS BAD CODE! REPAIR IT! 893 maxcones := List( maxcones, i -> List( i, j -> HomalgMap( HomalgMatrix( [ j ], HOMALG_MATRICES.ZZ ), 1 * HOMALG_MATRICES.ZZ, ContainingGrid( cone ) ) ) ); 894 895 maxcones := List( maxcones, i -> List( i, j -> UnderlyingListOfRingElementsInCurrentPresentation( ApplyMorphismToElement( ByASmallerPresentation( FactorGridMorphism( cone ) ), HomalgElement( j ) ) ) ) ); 896 897 maxcones := Fan( maxcones ); 898 899 return maxcones; 900 901end ); 902 903############################## 904## 905## Methods 906## 907############################## 908 909## 910InstallMethod( FourierProjection, 911 [ IsCone, IsInt ], 912 function( cone, n ) 913 local ray_generators, new_rays, i, j; 914 915 ray_generators:= RayGenerators( cone ); 916 917 new_rays:= [ ]; 918 919 for i in ray_generators do 920 921 j:= ShallowCopy( i ); 922 Remove(j, n ); 923 Add( new_rays, j ); 924 925 od; 926 927 return Cone( new_rays ); 928 929end ); 930 931InstallMethod( \*, 932 [ IsInt, IsCone ], 933 function( n, cone ) 934 935 if n>0 then 936 return cone; 937 elif n=0 then 938 return Cone( [ List([ 1 .. AmbientSpaceDimension( cone ) ], i-> 0 ) ] ); 939 else 940 return Cone( -RayGenerators( cone ) ); 941 fi; 942 943end ); 944 945## 946InstallMethod( \*, 947 " cartesian product for cones.", 948 [ IsCone, IsCone ], 949 950 function( cone1, cone2 ) 951 local rays1, rays2, i, j, raysnew; 952 953 rays1 := RayGenerators( cone1 ); 954 955 rays2 := RayGenerators( cone2 ); 956 957 rays1 := Concatenation( rays1, [ List( [ 1 .. Length( rays1[ 1 ] ) ], i -> 0 ) ] ); 958 959 rays2 := Concatenation( rays2, [ List( [ 1 .. Length( rays2[ 1 ] ) ], i -> 0 ) ] ); 960 961 raysnew := [ 1 .. Length( rays1 ) * Length( rays2 ) ]; 962 963 for i in [ 1 .. Length( rays1 ) ] do 964 965 for j in [ 1 .. Length( rays2 ) ] do 966 967 raysnew[ (i-1)*Length( rays2 ) + j ] := Concatenation( rays1[ i ], rays2[ j ] ); 968 969 od; 970 971 od; 972 973 raysnew := Cone( raysnew ); 974 975 SetContainingGrid( raysnew, ContainingGrid( cone1 ) + ContainingGrid( cone2 ) ); 976 977 return raysnew; 978 979end ); 980 981## 982InstallMethod( NonReducedInequalities, 983 "for a cone", 984 [ IsCone ], 985 986 function( cone ) 987 local inequalities; 988 989 if HasDefiningInequalities( cone ) then 990 991 return DefiningInequalities( cone ); 992 993 fi; 994 995 inequalities := [ ]; 996 997 if IsBound( cone!.input_equalities ) then 998 999 inequalities := Concatenation( inequalities, cone!.input_equalities, - cone!.input_equalities ); 1000 1001 fi; 1002 1003 if IsBound( cone!.input_inequalities ) then 1004 1005 inequalities := Concatenation( inequalities, cone!.input_inequalities ); 1006 1007 fi; 1008 1009 if inequalities <> [ ] then 1010 1011 return inequalities; 1012 1013 fi; 1014 1015 return DefiningInequalities( cone ); 1016 1017end ); 1018 1019InstallMethod( IntersectionOfCones, 1020 "for homalg cones", 1021 [ IsCone and HasDefiningInequalities, IsCone and HasDefiningInequalities ], 1022 1023 function( cone1, cone2 ) 1024 local inequalities, cone; 1025 1026 if not Rank( ContainingGrid( cone1 ) )= Rank( ContainingGrid( cone2 ) ) then 1027 1028 Error( "cones are not from the same grid" ); 1029 1030 fi; 1031 1032 inequalities := Unique( Concatenation( [ NonReducedInequalities( cone1 ), NonReducedInequalities( cone2 ) ] ) ); 1033 1034 cone := ConeByInequalities( inequalities ); 1035 1036 SetContainingGrid( cone, ContainingGrid( cone1 ) ); 1037 1038 return cone; 1039 1040end ); 1041 1042## 1043InstallMethod( IntersectionOfCones, 1044 "for homalg cones", 1045 [ IsCone, IsCone ], 1046 1047 function( cone1, cone2 ) 1048 local cone, ext_cone; 1049 1050 if not Rank( ContainingGrid( cone1 ) )= Rank( ContainingGrid( cone2 ) ) then 1051 1052 Error( "cones are not from the same grid" ); 1053 1054 fi; 1055 1056 ext_cone := Cdd_Intersection( ExternalCddCone( cone1), ExternalCddCone( cone2 ) ); 1057 1058 cone := Cone( ext_cone ); 1059 1060 SetContainingGrid( cone, ContainingGrid( cone1 ) ); 1061 1062 return cone; 1063 1064end ); 1065 1066## 1067InstallMethod( IntersectionOfCones, 1068 "for a list of convex cones", 1069 [ IsList ], 1070 1071 function( list_of_cones ) 1072 local grid, inequalities, equalities, i, cone; 1073 1074 if list_of_cones = [ ] then 1075 1076 Error( "Cannot create an empty cone\n" ); 1077 1078 fi; 1079 1080 if not ForAll( list_of_cones, IsCone ) then 1081 1082 Error( "not all list elements are cones\n" ); 1083 1084 fi; 1085 1086 grid := ContainingGrid( list_of_cones[ 1 ] ); 1087 1088 inequalities := [ ]; 1089 1090 equalities := [ ]; 1091 1092 for i in list_of_cones do 1093 1094 if HasDefiningInequalities( i ) then 1095 1096 Add( inequalities, DefiningInequalities( i ) ); 1097 1098 elif IsBound( i!.input_inequalities ) then 1099 1100 Add( inequalities, i!.input_inequalities ); 1101 1102 if IsBound( i!.input_equalities ) then 1103 1104 Add( equalities, i!.input_equalities ); 1105 1106 fi; 1107 1108 else 1109 1110 Add( inequalities, DefiningInequalities( i ) ); 1111 1112 fi; 1113 1114 od; 1115 1116 if equalities <> [ ] then 1117 1118 cone := ConeByEqualitiesAndInequalities( Concatenation( equalities ), Concatenation( inequalities ) ); 1119 1120 else 1121 1122 cone := ConeByInequalities( Concatenation( inequalities ) ); 1123 1124 fi; 1125 1126 SetContainingGrid( cone, grid ); 1127 1128 return cone; 1129 1130end ); 1131 1132## 1133InstallMethod( Contains, 1134 " for homalg cones", 1135 [ IsCone, IsCone ], 1136 1137 function( ambcone, cone ) 1138 local ineq; 1139 1140 if HasRayGenerators( ambcone ) and HasRayGenerators( cone ) then 1141 1142 if IsSubset( Set( RayGenerators( ambcone ) ), Set( RayGenerators( cone ) ) ) then 1143 1144 return true; 1145 1146 fi; 1147 1148 fi; 1149 1150 if HasDefiningInequalities( ambcone ) and HasDefiningInequalities( cone ) then 1151 1152 if IsSubset( Set( DefiningInequalities( cone ) ), Set( DefiningInequalities( ambcone ) ) ) then 1153 1154 return true; 1155 1156 fi; 1157 1158 fi; 1159 1160 ineq := NonReducedInequalities( ambcone ); 1161 1162 cone := RayGenerators( cone ); 1163 1164 ineq := List( cone, i -> ineq * i ); 1165 1166 ineq := Flat( ineq ); 1167 1168 return ForAll( ineq, i -> i >= 0 ); 1169 1170end ); 1171 1172## 1173InstallMethod( \*, 1174 "for a matrix and a cone", 1175 [ IsHomalgMatrix, IsCone ], 1176 1177 function( matrix, cone ) 1178 local ray_list, multiplied_rays, ring, new_cone; 1179 1180 ring := HomalgRing( matrix ); 1181 1182 ray_list := RayGenerators( cone ); 1183 1184 ray_list := List( ray_list, i -> HomalgMatrix( i, ring ) ); 1185 1186 multiplied_rays := List( ray_list, i -> matrix * i ); 1187 1188 multiplied_rays := List( multiplied_rays, i -> EntriesOfHomalgMatrix( i ) ); 1189 1190 new_cone := Cone( multiplied_rays ); 1191 1192 SetContainingGrid( new_cone, ContainingGrid( cone ) ); 1193 1194 return new_cone; 1195 1196end ); 1197 1198 1199 1200## 1201InstallMethod( \=, 1202 "for two cones", 1203 [ IsCone, IsCone ], 1204 1205 function( cone1, cone2 ) 1206 1207 return Contains( cone1, cone2 ) and Contains( cone2, cone1 ); 1208 1209end ); 1210 1211BindGlobal( "RayGeneratorContainedInCone", \in ); 1212 1213## 1214InstallMethod( \in, 1215 "for cones", 1216 [ IsList, IsCone ], 1217 1218 function( raygen, cone ) 1219 local ineq; 1220 1221 ineq := NonReducedInequalities( cone ); 1222 1223 ineq := List( ineq, i -> Sum( [ 1 .. Length( i ) ], j -> i[ j ]*raygen[ j ] ) ); 1224 1225 return ForAll( ineq, i -> i >= 0 ); 1226 1227end ); 1228 1229## 1230InstallMethod( \in, 1231 "for cones", 1232 [ IsList, IsCone and IsSimplicial ], 1233 1234 function( raygen, cone ) 1235 local ray_generators, matrix; 1236 1237 ray_generators := RayGenerators( cone ); 1238 1239 ##FIXME: One can use homalg for this, but at the moment 1240 ## we do not want the overhead. 1241 matrix := SolutionMat( ray_generators, raygen ); 1242 1243 return ForAll( matrix, i -> i >= 0 ); 1244 1245end ); 1246 1247## 1248InstallMethod( IsRelativeInteriorRay, 1249 "for cones", 1250 [ IsList, IsCone ], 1251 1252 function( ray, cone ) 1253 local ineq, mineq, equations; 1254 1255 ineq := NonReducedInequalities( cone ); 1256 1257 mineq := -ineq; 1258 1259 equations := Intersection( ineq, mineq ); 1260 1261 ineq := Difference( ineq, equations ); 1262 1263 ineq := List( ineq, i -> i * ray ); 1264 1265 equations := List( equations, i -> i * ray ); 1266 1267 return ForAll( ineq, i -> i > 0 ) and ForAll( equations, i -> i = 0 ); 1268 1269end ); 1270 1271 1272####################################### 1273## 1274## Methods to construct external cones 1275## 1276####################################### 1277 1278InstallMethod( ExternalCddCone, 1279 [ IsCone ], 1280 1281 function( cone ) 1282 1283 local list, new_list, number_of_equalities, linearity, i, u ; 1284 1285 new_list:= [ ]; 1286 if IsBound( cone!.input_rays ) and Length( cone!.input_rays )= 1 and IsZero( cone!.input_rays ) then 1287 1288 new_list:= [ Concatenation( [ 1 ], cone!.input_rays[ 1 ] ) ]; 1289 1290 return Cdd_PolyhedronByGenerators( new_list ); 1291 1292 fi; 1293 1294 if IsBound( cone!.input_rays ) then 1295 1296 list := cone!.input_rays; 1297 1298 for i in [1..Length( list ) ] do 1299 1300 u:= ShallowCopy( list[ i ] ); 1301 1302 Add( u, 0, 1 ); 1303 1304 Add( new_list, u ); 1305 1306 od; 1307 1308 return Cdd_PolyhedronByGenerators( new_list ); 1309 1310 fi; 1311 1312 1313 if IsBound( cone!.input_equalities ) then 1314 1315 list := StructuralCopy( cone!.input_equalities ); 1316 1317 number_of_equalities:= Length( list ); 1318 1319 linearity := [1..number_of_equalities]; 1320 1321 Append( list, StructuralCopy( cone!.input_inequalities ) ); 1322 1323 for i in [1..Length( list ) ] do 1324 1325 u:= ShallowCopy( list[ i ] ); 1326 1327 Add( u, 0, 1 ); 1328 1329 Add( new_list, u ); 1330 1331 od; 1332 1333 return Cdd_PolyhedronByInequalities( new_list, linearity ); 1334 1335 else 1336 1337 list:= StructuralCopy( cone!.input_inequalities ); 1338 1339 for i in [1..Length( list ) ] do 1340 1341 u:= ShallowCopy( list[ i ] ); 1342 1343 Add( u, 0, 1 ); 1344 1345 Add( new_list, u ); 1346 1347 od; 1348 1349 return Cdd_PolyhedronByInequalities( new_list ); 1350 1351 fi; 1352 1353end ); 1354 1355 1356InstallMethod( ExternalNmzCone, 1357 [ IsCone ], 1358 function( cone ) 1359 local a, list, equalities, i; 1360 1361 list:= []; 1362 1363 if IsBound( cone!.input_rays ) then 1364 1365 list := StructuralCopy( cone!.input_rays ); 1366 1367 if ForAll( Concatenation( list ), IsInt ) then 1368 1369 return ValueGlobal( "NmzCone" )( [ "integral_closure", list ] ); 1370 1371 else 1372 1373 a := DuplicateFreeList( List( Concatenation( list ), l -> DenominatorRat( l ) ) ); 1374 1375 list := Lcm( a ) * list; 1376 1377 return ValueGlobal( "NmzCone" )( [ "integral_closure", list ] ); 1378 1379 fi; 1380 1381 fi; 1382 1383 list:= StructuralCopy( cone!.input_inequalities ); 1384 1385 if IsBound( cone!.input_equalities ) then 1386 1387 equalities:= StructuralCopy( cone!.input_equalities ); 1388 1389 for i in equalities do 1390 1391 Append( list, [ i, -i ] ); 1392 1393 od; 1394 1395 fi; 1396 1397 if ForAll( Concatenation( list ), IsInt ) then 1398 1399 return ValueGlobal( "NmzCone" )( ["inequalities", list ] ); 1400 1401 else 1402 1403 a := DuplicateFreeList( List( Concatenation( list ), l -> DenominatorRat( l ) ) ); 1404 1405 list := Lcm( a ) * list; 1406 1407 return ValueGlobal( "NmzCone" )( ["inequalities", list ] ); 1408 1409 fi; 1410 1411 Error( "The cone should be defined by vertices or inequalities!" ); 1412 1413end ); 1414 1415##################################### 1416## 1417## Constructors 1418## 1419##################################### 1420 1421## 1422InstallMethod( ConeByInequalities, 1423 "constructor for Cones by inequalities", 1424 [ IsList ], 1425 1426 function( inequalities ) 1427 local cone, newgens, i, vals; 1428 1429 if Length( inequalities ) = 0 then 1430 1431 Error( "someone must set a dimension\n" ); 1432 1433 fi; 1434 1435 cone := rec( input_inequalities := inequalities ); 1436 1437 ObjectifyWithAttributes( 1438 cone, TheTypeConvexCone 1439 ); 1440 1441 SetAmbientSpaceDimension( cone, Length( inequalities[ 1 ] ) ); 1442 1443 return cone; 1444 1445end ); 1446 1447## 1448InstallMethod( ConeByEqualitiesAndInequalities, 1449 "constructor for cones by equalities and inequalities", 1450 [ IsList, IsList ], 1451 1452 function( equalities, inequalities ) 1453 local cone, newgens, i, vals; 1454 1455 if Length( equalities ) = 0 and Length( inequalities ) = 0 then 1456 1457 Error( "someone must set a dimension\n" ); 1458 1459 fi; 1460 1461 cone := rec( input_inequalities := inequalities, input_equalities := equalities ); 1462 1463 ObjectifyWithAttributes( 1464 cone, TheTypeConvexCone 1465 ); 1466 1467 if Length( equalities ) > 0 then 1468 1469 SetAmbientSpaceDimension( cone, Length( equalities[ 1 ] ) ); 1470 1471 else 1472 1473 SetAmbientSpaceDimension( cone, Length( inequalities[ 1 ] ) ); 1474 1475 fi; 1476 1477 return cone; 1478 1479end ); 1480 1481## 1482InstallMethod( ConeByGenerators, 1483 "constructor for cones by generators", 1484 [ IsList ], 1485 1486 function( raylist ) 1487 local cone, newgens, i, vals; 1488 1489 if Length( raylist ) = 0 then 1490 1491 # Think again about this 1492 return ConeByGenerators( [ [ ] ] ); 1493 1494 fi; 1495 1496 newgens := [ ]; 1497 1498 for i in raylist do 1499 1500 if IsList( i ) then 1501 1502 Add( newgens, i ); 1503 1504 elif IsCone( i ) then 1505 1506 Append( newgens, RayGenerators( i ) ); 1507 1508 else 1509 1510 Error( " wrong rays" ); 1511 1512 fi; 1513 1514 od; 1515 1516 cone := rec( input_rays := newgens ); 1517 1518 ObjectifyWithAttributes( 1519 cone, TheTypeConvexCone 1520 ); 1521 1522 if Length( raylist ) =1 and IsZero( raylist[ 1 ] ) then 1523 1524 SetIsZero( cone, true ); 1525 1526 fi; 1527 1528 newgens := Set( newgens ); 1529 1530 SetAmbientSpaceDimension( cone, Length( newgens[ 1 ] ) ); 1531 1532 if Length( newgens ) = 1 and not Set( newgens[ 1 ] ) = [ 0 ] then 1533 1534 SetIsRay( cone, true ); 1535 1536 else 1537 1538 SetIsRay( cone, false ); 1539 1540 fi; 1541 1542 return cone; 1543 1544end ); 1545 1546## 1547InstallMethod( Cone, 1548 "construct cone from list of generating rays", 1549 [ IsList ], 1550 1551 ConeByGenerators 1552 1553 ); 1554 1555InstallMethod( Cone, 1556 "Construct cone from Cdd cone", 1557 [ IsCddPolyhedron ], 1558 1559 function( cdd_cone ) 1560 local inequalities, equalities, 1561 new_inequalities, new_equalities, u, i; 1562 1563 if cdd_cone!.rep_type = "H-rep" then 1564 1565 inequalities:= Cdd_Inequalities( cdd_cone ); 1566 1567 new_inequalities:= [ ]; 1568 1569 for i in inequalities do 1570 1571 u:= ShallowCopy( i ); 1572 1573 Remove( u , 1 ); 1574 1575 Add(new_inequalities, u ); 1576 1577 od; 1578 1579 if cdd_cone!.linearity <> [] then 1580 1581 equalities:= Cdd_Equalities( cdd_cone ); 1582 1583 new_equalities:= [ ]; 1584 1585 for i in equalities do 1586 1587 u:= ShallowCopy( i ); 1588 1589 Remove( u , 1 ); 1590 1591 Add(new_equalities, u ); 1592 1593 od; 1594 1595 return ConeByEqualitiesAndInequalities( new_equalities, new_inequalities); 1596 1597 fi; 1598 1599 return ConeByInequalities( new_inequalities ); 1600 1601 else 1602 1603 return ConeByGenerators( Cdd_GeneratingRays( cdd_cone ) ); 1604 1605 fi; 1606 1607end ); 1608 1609 1610 1611################################ 1612## 1613## Displays and Views 1614## 1615################################ 1616 1617## 1618InstallMethod( ViewObj, 1619 "for homalg cones", 1620 [ IsCone ], 1621 1622 function( cone ) 1623 local str; 1624 1625 Print( "<A" ); 1626 1627 if HasIsSmooth( cone ) then 1628 1629 if IsSmooth( cone ) then 1630 1631 Print( " smooth" ); 1632 1633 fi; 1634 1635 fi; 1636 1637 if HasIsPointed( cone ) then 1638 1639 if IsPointed( cone ) then 1640 1641 Print( " pointed" ); 1642 1643 fi; 1644 1645 fi; 1646 1647 if HasIsSimplicial( cone ) then 1648 1649 if IsSimplicial( cone ) then 1650 1651 Print( " simplicial" ); 1652 1653 fi; 1654 1655 fi; 1656 1657 if HasIsRay( cone ) and IsRay( cone ) then 1658 1659 Print( " ray" ); 1660 1661 else 1662 1663 Print( " cone" ); 1664 1665 fi; 1666 1667 Print( " in |R^" ); 1668 1669 Print( String( AmbientSpaceDimension( cone ) ) ); 1670 1671 if HasDimension( cone ) then 1672 1673 Print( " of dimension ", String( Dimension( cone ) ) ); 1674 1675 fi; 1676 1677 if HasRayGenerators( cone ) then 1678 1679 Print( " with ", String( Length( RayGenerators( cone ) ) )," ray generators" ); 1680 1681 fi; 1682 1683 Print( ">" ); 1684 1685end ); 1686 1687## 1688InstallMethod( Display, 1689 "for homalg cones", 1690 [ IsCone ], 1691 1692 function( cone ) 1693 local str; 1694 1695 Print( "A" ); 1696 1697 if HasIsSmooth( cone ) then 1698 1699 if IsSmooth( cone ) then 1700 1701 Print( " smooth" ); 1702 1703 fi; 1704 1705 fi; 1706 1707 if HasIsPointed( cone ) then 1708 1709 if IsPointed( cone ) then 1710 1711 Print( " pointed" ); 1712 1713 fi; 1714 1715 fi; 1716 1717 if IsRay( cone ) then 1718 1719 Print( " ray" ); 1720 1721 else 1722 1723 Print( " cone" ); 1724 1725 fi; 1726 1727 Print( " in |R^" ); 1728 1729 Print( String( AmbientSpaceDimension( cone ) ) ); 1730 1731 if HasDimension( cone ) then 1732 1733 Print( " of dimension ", String( Dimension( cone ) ) ); 1734 1735 fi; 1736 1737 if HasRayGenerators( cone ) then 1738 1739 Print( " with ray generators ", String( RayGenerators( cone ) ) ); 1740 1741 fi; 1742 1743 Print( ".\n" ); 1744 1745end ); 1746