DGtal  1.5.beta
Functions
testEhrhartPolynomial.cpp File Reference
DGtal Test Files
#include <iostream>
#include <vector>
#include <algorithm>
#include "DGtal/base/Common.h"
#include "DGtal/geometry/volumes/EhrhartPolynomial.h"
#include "DGtal/geometry/volumes/DigitalConvexity.h"
#include "DGtalCatch.h"
Include dependency graph for testEhrhartPolynomial.cpp:

Go to the source code of this file.

Functions

 SCENARIO ("EhrhartPolynomial< Z2 > unit tests", "[ehrhart_polynomial][2d]")
 
 SCENARIO ("EhrhartPolynomial< Z3 > unit tests", "[ehrhart_polynomial][3d]")
 
 SCENARIO ("EhrhartPolynomial< Z4 > unit tests", "[ehrhart_polynomial][4d]")
 
 SCENARIO ("EhrhartPolynomial< Z2 > triangle tests", "[ehrhart_polynomial][2d]")
 

Detailed Description

This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Author
Jacques-Olivier Lachaud (jacqu.nosp@m.es-o.nosp@m.livie.nosp@m.r.la.nosp@m.chaud.nosp@m.@uni.nosp@m.v-sav.nosp@m.oie..nosp@m.fr ) Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
Date
2021/08/24

Functions for testing class EhrhartPolynomial.

This file is part of the DGtal library.

Definition in file testEhrhartPolynomial.cpp.

Function Documentation

◆ SCENARIO() [1/4]

SCENARIO ( "EhrhartPolynomial< Z2 > triangle tests"  ,
""  [ehrhart_polynomial][2d] 
)

Definition at line 192 of file testEhrhartPolynomial.cpp.

193 {
194  typedef SpaceND< 2, int > Space;
195  typedef HyperRectDomain< Space > Domain;
196  typedef KhalimskySpaceND< 2, int > KSpace;
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 }
DGtal::DigitalConvexity::isSimplexFullDimensional
static bool isSimplexFullDimensional(PointIterator itB, PointIterator itE)
DGtal::DigitalConvexity::makeSimplex
static LatticePolytope makeSimplex(PointIterator itB, PointIterator itE)
DGtal::DigitalConvexity::isFullyConvex
bool isFullyConvex(const PointRange &X, bool convex0=false) const
DGtal::EhrhartPolynomial
Aim: This class implements the class Ehrhart Polynomial which is related to lattice point enumeration...
Definition: EhrhartPolynomial.h:67
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:394
DGtal::Trace::info
std::ostream & info()
DGtal::int64_t
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
DGtal::trace
Trace trace
Definition: Common.h:153
Point
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE
CAPTURE(thicknessHV)
domain
Domain domain
Definition: testProjection.cpp:88
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
REQUIRE
REQUIRE(domain.isInside(aPoint))

References CAPTURE(), domain, DGtal::Trace::info(), DGtal::DigitalConvexity< TKSpace >::isFullyConvex(), DGtal::DigitalConvexity< TKSpace >::isSimplexFullDimensional(), DGtal::DigitalConvexity< TKSpace >::makeSimplex(), REQUIRE(), and DGtal::trace.

◆ SCENARIO() [2/4]

SCENARIO ( "EhrhartPolynomial< Z2 > unit tests"  ,
""  [ehrhart_polynomial][2d] 
)

Definition at line 48 of file testEhrhartPolynomial.cpp.

49 {
50  typedef SpaceND< 2, int > Space;
51  typedef KhalimskySpaceND< 2, int > KSpace;
52  typedef Space::Point Point;
53  typedef DGtal::int64_t Integer;
54  typedef EhrhartPolynomial< Space, Integer > Ehrhart;
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 }
GIVEN
GIVEN("A cubical complex with random 3-cells")
Definition: testCubicalComplex.cpp:70

References GIVEN(), DGtal::DigitalConvexity< TKSpace >::makeSimplex(), and REQUIRE().

◆ SCENARIO() [3/4]

SCENARIO ( "EhrhartPolynomial< Z3 > unit tests"  ,
""  [ehrhart_polynomial][3d] 
)

Definition at line 97 of file testEhrhartPolynomial.cpp.

98 {
99  typedef SpaceND< 3, int > Space;
100  typedef KhalimskySpaceND< 3, int > KSpace;
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 }
DGtal::DigitalConvexity::makePolytope
LatticePolytope makePolytope(const PointRange &X, bool make_minkowski_summable=false) const

References CAPTURE(), GIVEN(), DGtal::DigitalConvexity< TKSpace >::makePolytope(), and REQUIRE().

◆ SCENARIO() [4/4]

SCENARIO ( "EhrhartPolynomial< Z4 > unit tests"  ,
""  [ehrhart_polynomial][4d] 
)

Definition at line 145 of file testEhrhartPolynomial.cpp.

146 {
147  typedef SpaceND< 4, int > Space;
148  typedef KhalimskySpaceND< 4, int > KSpace;
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 }

References CAPTURE(), GIVEN(), DGtal::DigitalConvexity< TKSpace >::makePolytope(), and REQUIRE().