DGtal  1.5.beta
testEhrhartPolynomial.cpp
Go to the documentation of this file.
1 
31 #include <iostream>
32 #include <vector>
33 #include <algorithm>
34 #include "DGtal/base/Common.h"
35 #include "DGtal/geometry/volumes/EhrhartPolynomial.h"
36 #include "DGtal/geometry/volumes/DigitalConvexity.h"
37 #include "DGtalCatch.h"
39 
40 using namespace std;
41 using namespace DGtal;
42 
43 
45 // Functions for testing class EhrhartPolynomial.
47 
48 SCENARIO( "EhrhartPolynomial< Z2 > unit tests", "[ehrhart_polynomial][2d]" )
49 {
50  typedef SpaceND< 2, int > Space;
52  typedef Space::Point Point;
53  typedef DGtal::int64_t Integer;
55 
56  GIVEN( "Simplex (0,0), (1,0), (2,1)" ) {
57  std::vector< Point > T = { Point(0,0), Point(1,0), Point(2,1) };
58  DigitalConvexity< KSpace > dconv( Point::diagonal( -100 ), Point::diagonal( 100 ) );
59  auto P = dconv.makeSimplex( T.cbegin(), T.cend() );
60  Ehrhart E( P );
61  THEN( "Its Ehrhart polynomial is 1/2( 2+ 3t+t^2) ") {
62  auto expP = mmonomial<Integer>( 2 ) + 3 * mmonomial<Integer>( 1 ) + 2;
63  REQUIRE( E.numerator() == expP );
64  REQUIRE( E.denominator() == 2 );
65  }
66  THEN( "Its Ehrhart polynomial counts lattice points" ) {
67  auto P0 = 0 * P;
68  auto n0 = E.count( 0 );
69  REQUIRE( P0.count() == n0 );
70  auto P1 = 1 * P;
71  auto n1 = E.count( 1 );
72  REQUIRE( P1.count() == n1 );
73  auto P2 = 2 * P;
74  auto n2 = E.count( 2 );
75  REQUIRE( P2.count() == n2 );
76  auto P3 = 3 * P;
77  auto n3 = E.count( 3 );
78  REQUIRE( P3.count() == n3 );
79  }
80  THEN( "Its Ehrhart polynomial counts interior lattice points" ) {
81  auto P1 = 1 * P;
82  auto n1 = E.countInterior( 1 );
83  REQUIRE( P1.countInterior() == n1 );
84  auto P2 = 2 * P;
85  auto n2 = E.countInterior( 2 );
86  REQUIRE( P2.countInterior() == n2 );
87  auto P3 = 3 * P;
88  auto n3 = E.countInterior( 3 );
89  REQUIRE( P3.countInterior() == n3 );
90  auto P4 = 4 * P;
91  auto n4 = E.countInterior( 4 );
92  REQUIRE( P4.countInterior() == n4 );
93  }
94  }
95 }
96 
97 SCENARIO( "EhrhartPolynomial< Z3 > unit tests", "[ehrhart_polynomial][3d]" )
98 {
99  typedef SpaceND< 3, int > Space;
101  typedef Space::Point Point;
102  typedef DGtal::int64_t Integer;
103  typedef EhrhartPolynomial< Space, Integer > Ehrhart;
104 
105  GIVEN( "A convex polytope" ) {
106  std::vector< Point > T = { Point(0,0,0), Point(1,0,1), Point(2,1,0), Point(0,0,2),
107  Point(0,3,1) };
108  DigitalConvexity< KSpace > dconv( Point::diagonal( -100 ), Point::diagonal( 100 ) );
109  auto P = dconv.makePolytope( T );
110  Ehrhart E( P );
111  CAPTURE( E.numerator() );
112  CAPTURE( E.denominator() );
113  THEN( "Its Ehrhart polynomial counts lattice points" ) {
114  auto P0 = 0 * P;
115  auto n0 = E.count( 0 );
116  REQUIRE( P0.count() == n0 );
117  auto P1 = 1 * P;
118  auto n1 = E.count( 1 );
119  REQUIRE( P1.count() == n1 );
120  auto P2 = 2 * P;
121  auto n2 = E.count( 2 );
122  REQUIRE( P2.count() == n2 );
123  auto P3 = 3 * P;
124  auto n3 = E.count( 3 );
125  REQUIRE( P3.count() == n3 );
126  }
127  THEN( "Its Ehrhart polynomial counts interior lattice points" ) {
128  auto P1 = 1 * P;
129  auto n1 = E.countInterior( 1 );
130  REQUIRE( P1.countInterior() == n1 );
131  auto P2 = 2 * P;
132  auto n2 = E.countInterior( 2 );
133  REQUIRE( P2.countInterior() == n2 );
134  auto P3 = 3 * P;
135  auto n3 = E.countInterior( 3 );
136  REQUIRE( P3.countInterior() == n3 );
137  auto P4 = 4 * P;
138  auto n4 = E.countInterior( 4 );
139  REQUIRE( P4.countInterior() == n4 );
140  }
141  }
142 }
143 
144 
145 SCENARIO( "EhrhartPolynomial< Z4 > unit tests", "[ehrhart_polynomial][4d]" )
146 {
147  typedef SpaceND< 4, int > Space;
149  typedef Space::Point Point;
150  typedef DGtal::int64_t Integer;
151  typedef EhrhartPolynomial< Space, Integer > Ehrhart;
152 
153  GIVEN( "A convex polytope" ) {
154  std::vector< Point > T = { Point(0,0,0,0), Point(1,0,1,-1), Point(2,1,0,2),
155  Point(0,0,2,0), Point(0,4,1,3), Point(3,0,-1,1) };
156  DigitalConvexity< KSpace > dconv( Point::diagonal( -100 ), Point::diagonal( 100 ) );
157  auto P = dconv.makePolytope( T );
158  Ehrhart E( P );
159  CAPTURE( E.numerator() );
160  CAPTURE( E.denominator() );
161  THEN( "Its Ehrhart polynomial counts lattice points" ) {
162  auto P0 = 0 * P;
163  auto n0 = E.count( 0 );
164  REQUIRE( P0.count() == n0 );
165  auto P1 = 1 * P;
166  auto n1 = E.count( 1 );
167  REQUIRE( P1.count() == n1 );
168  auto P2 = 2 * P;
169  auto n2 = E.count( 2 );
170  REQUIRE( P2.count() == n2 );
171  auto P3 = 3 * P;
172  auto n3 = E.count( 3 );
173  REQUIRE( P3.count() == n3 );
174  }
175  THEN( "Its Ehrhart polynomial counts interior lattice points" ) {
176  auto P1 = 1 * P;
177  auto n1 = E.countInterior( 1 );
178  REQUIRE( P1.countInterior() == n1 );
179  auto P2 = 2 * P;
180  auto n2 = E.countInterior( 2 );
181  REQUIRE( P2.countInterior() == n2 );
182  auto P3 = 3 * P;
183  auto n3 = E.countInterior( 3 );
184  REQUIRE( P3.countInterior() == n3 );
185  auto P4 = 4 * P;
186  auto n4 = E.countInterior( 4 );
187  REQUIRE( P4.countInterior() == n4 );
188  }
189  }
190 }
191 
192 SCENARIO( "EhrhartPolynomial< Z2 > triangle tests", "[ehrhart_polynomial][2d]" )
193 {
194  typedef SpaceND< 2, int > Space;
197  typedef Space::Point Point;
198  typedef DGtal::int64_t Integer;
199  typedef EhrhartPolynomial< Space, Integer > Ehrhart;
200  typedef Ehrhart::LatticePolytope Polytope;
201  typedef Polytope::UnitSegment UnitSegment;
202 
203  Domain domain( Point( 0, 0 ), Point( 4, 4 ) );
204  DigitalConvexity<KSpace> dconv( Point( -1, -1 ), Point( 5, 5 ) );
205 
206  WHEN( "Computing all triangles in domain (0,0)-(4,4)." ) {
207  unsigned int nbsimplex = 0;
208  unsigned int nb0_ok = 0;
209  unsigned int nb1_ok = 0;
210  unsigned int nb_cvx_ok = 0;
211  for ( auto a : domain )
212  for ( auto b : domain )
213  for ( auto c : domain )
214  {
215  if ( ! ( ( a < b ) && ( b < c ) ) ) continue;
216  if ( ! dconv.isSimplexFullDimensional( { a, b, c } ) ) continue;
217  const auto P = dconv.makeSimplex( { a, b, c } );
218  const bool fcvx = dconv.isFullyConvex( P );
219  const auto P0 = P + UnitSegment( 0 );
220  const auto P1 = P + UnitSegment( 1 );
221  const auto D = P.getDomain();
222  const auto W = D.upperBound() - D.lowerBound();
223  const auto E_P = Ehrhart( P );
224  const auto E_P0 = Ehrhart( P0 );
225  const auto E_P1 = Ehrhart( P1 );
226  const auto LP = E_P.numerator();
227  const auto LP0 = E_P0.numerator();
228  const auto LP1 = E_P1.numerator();
229  const auto LD = E_P.denominator();
230  const auto LD0 = E_P0.denominator();
231  const auto LD1 = E_P1.denominator();
232  const bool c2_0 = ( LP[ 2 ] + Integer( W[ 1 ] ) * LD ) == LP0[ 2 ];
233  const bool c1_0 = ( LP[ 1 ] + Integer( 1 ) * LD ) == LP0[ 1 ];
234  const bool c0_0 = LP[ 0 ] == LP0[ 0 ];
235  const bool d_0 = LD == LD0;
236  const bool c2_1 = ( LP[ 2 ] + Integer( W[ 0 ] ) * LD ) == LP1[ 2 ];
237  const bool c1_1 = ( LP[ 1 ] + Integer( 1 ) * LD ) == LP1[ 1 ];
238  const bool c0_1 = LP[ 0 ] == LP1[ 0 ];
239  const bool d_1 = LD == LD1;
240  nb0_ok += ( c2_0 && c1_0 && c0_0 && d_0 ) ? 1 : 0;
241  nb1_ok += ( c2_1 && c1_1 && c0_1 && d_1 ) ? 1 : 0;
242  nbsimplex += 1;
243  nb_cvx_ok += fcvx ? 1 : 0;
244  if ( ! ( c2_0 && c1_0 && c0_0 && d_0 ) ){
245  trace.info() << "------------------------------------" << std::endl;
246  trace.info() << ( fcvx ? "FULLCVX" : "NOTFULLCVX" ) << a << b << c << std::endl;
247  trace.info() << "W[0] = " << W[ 0 ] << "W[1] = " << W[ 1 ] << std::endl;
248  trace.info() << "LP = (" << LP << ")/" << E_P.denominator() << std::endl;
249  trace.info() << "LP0= (" << LP0 << ")/" << E_P0.denominator() << std::endl;
250  break;
251  }
252  if ( ! ( c2_1 && c1_1 && c0_1 && d_1 ) ){
253  trace.info() << "------------------------------------" << std::endl;
254  trace.info() << ( fcvx ? "FULLCVX" : "NOTFULLCVX" ) << a << b << c << std::endl;
255  trace.info() << "W[0] = " << W[ 0 ] << "W[1] = " << W[ 1 ] << std::endl;
256  trace.info() << "LP = (" << LP << ")/" << E_P.denominator() << std::endl;
257  trace.info() << "LP1= (" << LP1 << ")/" << E_P1.denominator() << std::endl;
258  break;
259  }
260  }
261  THEN( "Not all triangles are fully convex." ) {
262  REQUIRE( nb_cvx_ok < nbsimplex );
263  }
264  THEN( "Ehrhart polynomials are predictable." ) {
265  CAPTURE( nb_cvx_ok );
266  REQUIRE( nb0_ok == nbsimplex );
267  REQUIRE( nb1_ok == nbsimplex );
268  }
269  }
270 }
static bool isSimplexFullDimensional(PointIterator itB, PointIterator itE)
static LatticePolytope makeSimplex(PointIterator itB, PointIterator itE)
bool isFullyConvex(const PointRange &X, bool convex0=false) const
LatticePolytope makePolytope(const PointRange &X, bool make_minkowski_summable=false) const
Aim: This class implements the class Ehrhart Polynomial which is related to lattice point enumeration...
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
std::ostream & info()
DGtal is the top-level namespace which contains all DGtal functions and types.
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
Trace trace
Definition: Common.h:153
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE(thicknessHV)
GIVEN("A cubical complex with random 3-cells")
SCENARIO("EhrhartPolynomial< Z2 > unit tests", "[ehrhart_polynomial][2d]")
Domain domain
HyperRectDomain< Space > Domain
REQUIRE(domain.isInside(aPoint))