DGtal 1.3.0
Searching...
No Matches
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
40using namespace std;
41using namespace DGtal;
42
43
45// Functions for testing class EhrhartPolynomial.
47
48SCENARIO( "EhrhartPolynomial< Z2 > unit tests", "[ehrhart_polynomial][2d]" )
49{
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
97SCENARIO( "EhrhartPolynomial< Z3 > unit tests", "[ehrhart_polynomial][3d]" )
98{
101 typedef Space::Point Point;
102 typedef DGtal::int64_t Integer;
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
145SCENARIO( "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;
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
192SCENARIO( "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;
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 unsigned int nb_not_cvx_ok = 0;
212 for ( auto a : domain )
213 for ( auto b : domain )
214 for ( auto c : domain )
215 {
216 if ( ! ( ( a < b ) && ( b < c ) ) ) continue;
217 if ( ! dconv.isSimplexFullDimensional( { a, b, c } ) ) continue;
218 const auto P = dconv.makeSimplex( { a, b, c } );
219 const bool fcvx = dconv.isFullyConvex( P );
220 const auto P0 = P + UnitSegment( 0 );
221 const auto P1 = P + UnitSegment( 1 );
222 const auto D = P.getDomain();
223 const auto W = D.upperBound() - D.lowerBound();
224 const auto E_P = Ehrhart( P );
225 const auto E_P0 = Ehrhart( P0 );
226 const auto E_P1 = Ehrhart( P1 );
227 const auto LP = E_P.numerator();
228 const auto LP0 = E_P0.numerator();
229 const auto LP1 = E_P1.numerator();
230 const auto LD = E_P.denominator();
231 const auto LD0 = E_P0.denominator();
232 const auto LD1 = E_P1.denominator();
233 const bool c2_0 = ( LP[ 2 ] + Integer( W[ 1 ] ) * LD ) == LP0[ 2 ];
234 const bool c1_0 = ( LP[ 1 ] + Integer( 1 ) * LD ) == LP0[ 1 ];
235 const bool c0_0 = LP[ 0 ] == LP0[ 0 ];
236 const bool d_0 = LD == LD0;
237 const bool c2_1 = ( LP[ 2 ] + Integer( W[ 0 ] ) * LD ) == LP1[ 2 ];
238 const bool c1_1 = ( LP[ 1 ] + Integer( 1 ) * LD ) == LP1[ 1 ];
239 const bool c0_1 = LP[ 0 ] == LP1[ 0 ];
240 const bool d_1 = LD == LD1;
241 nb0_ok += ( c2_0 && c1_0 && c0_0 && d_0 ) ? 1 : 0;
242 nb1_ok += ( c2_1 && c1_1 && c0_1 && d_1 ) ? 1 : 0;
243 nbsimplex += 1;
244 nb_cvx_ok += fcvx ? 1 : 0;
245 if ( ! ( c2_0 && c1_0 && c0_0 && d_0 ) ){
246 trace.info() << "------------------------------------" << std::endl;
247 trace.info() << ( fcvx ? "FULLCVX" : "NOTFULLCVX" ) << a << b << c << std::endl;
248 trace.info() << "W[0] = " << W[ 0 ] << "W[1] = " << W[ 1 ] << std::endl;
249 trace.info() << "LP = (" << LP << ")/" << E_P.denominator() << std::endl;
250 trace.info() << "LP0= (" << LP0 << ")/" << E_P0.denominator() << std::endl;
251 break;
252 }
253 if ( ! ( c2_1 && c1_1 && c0_1 && d_1 ) ){
254 trace.info() << "------------------------------------" << std::endl;
255 trace.info() << ( fcvx ? "FULLCVX" : "NOTFULLCVX" ) << a << b << c << std::endl;
256 trace.info() << "W[0] = " << W[ 0 ] << "W[1] = " << W[ 1 ] << std::endl;
257 trace.info() << "LP = (" << LP << ")/" << E_P.denominator() << std::endl;
258 trace.info() << "LP1= (" << LP1 << ")/" << E_P1.denominator() << std::endl;
259 break;
260 }
261 }
262 THEN( "Not all triangles are fully convex." ) {
263 REQUIRE( nb_cvx_ok < nbsimplex );
264 }
265 THEN( "Ehrhart polynomials are predictable." ) {
266 CAPTURE( nb_cvx_ok );
267 REQUIRE( nb0_ok == nbsimplex );
268 REQUIRE( nb1_ok == nbsimplex );
269 }
270 }
271}
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()
Point::Coordinate Integer
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:154
STL namespace.
MyPointD Point
Definition: testClone2.cpp:383
CAPTURE(thicknessHV)
GIVEN("A cubical complex with random 3-cells")
Domain domain
HyperRectDomain< Space > Domain
REQUIRE(domain.isInside(aPoint))
SCENARIO("UnorderedSetByBlock< PointVector< 2, int > unit tests with 32 bits blocks", "[unorderedsetbyblock][2d]")