DGtal 1.4.2
Loading...
Searching...
No Matches
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
Aim: A helper class to build polytopes from digital sets and to check digital k-convexity and full co...
Definition DigitalConvexity.h:78
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}
DGtal::mmonomial
MPolynomial< 1, Ring, Alloc > mmonomial(unsigned int e)
Definition MPolynomial.h:1696
GIVEN
GIVEN("A cubical complex with random 3-cells")
Definition testCubicalComplex.cpp:70

References GIVEN(), DGtal::DigitalConvexity< TKSpace >::makeSimplex(), DGtal::mmonomial(), 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().