@ -23,15 +23,15 @@ namespace storm {
STORM_LOG_DEBUG ( " Invoked QuickHull on " < < points . size ( ) < < " points " ) ;
const uint_fast64_t dimension = points . front ( ) . rows ( ) ;
if ( dimension = = 1 ) {
if ( dimension = = 1 ) {
handle1DPoints ( points , generateRelevantVerticesAndVertexSets ) ;
} else {
// Generate initial set of d+1 affine independent points (if such a set exists)
std : : vector < uint_fast64_t > vertexIndices ;
if ( this - > findInitialVertices ( points , vertexIndices ) ) {
if ( this - > findInitialVertices ( points , vertexIndices ) ) {
// compute point inside initial facet
EigenVector insidePoint ( EigenVector : : Zero ( dimension ) ) ;
for ( uint_fast64_t vertexIndex : vertexIndices ) {
for ( uint_fast64_t vertexIndex : vertexIndices ) {
insidePoint + = points [ vertexIndex ] ;
}
insidePoint / = storm : : utility : : convertNumber < ValueType > ( ( uint_fast64_t ) vertexIndices . size ( ) ) ;
@ -79,13 +79,13 @@ namespace storm {
ValueType maxValue = points . front ( ) ( 0 ) ;
uint_fast64_t minIndex = 0 ;
uint_fast64_t maxIndex = 0 ;
for ( uint_fast64_t pointIndex = 1 ; pointIndex < points . size ( ) ; + + pointIndex ) {
for ( uint_fast64_t pointIndex = 1 ; pointIndex < points . size ( ) ; + + pointIndex ) {
ValueType const & pointValue = points [ pointIndex ] ( 0 ) ;
if ( pointValue < minValue ) {
if ( pointValue < minValue ) {
minValue = pointValue ;
minIndex = pointIndex ;
}
if ( pointValue > maxValue ) {
if ( pointValue > maxValue ) {
maxValue = pointValue ;
maxIndex = pointIndex ;
}
@ -96,11 +96,11 @@ namespace storm {
resultVector = EigenVector ( 2 ) ;
resultVector ( 0 ) = - minValue ;
resultVector ( 1 ) = maxValue ;
if ( generateRelevantVerticesAndVertexSets ) {
if ( generateRelevantVerticesAndVertexSets ) {
relevantVertices . push_back ( points [ minIndex ] ) ;
std : : vector < uint_fast64_t > minVertexSet ( 1 , 0 ) ;
vertexSets . push_back ( std : : move ( minVertexSet ) ) ;
if ( minIndex ! = maxIndex & & points [ minIndex ] ! = points [ maxIndex ] ) {
if ( minIndex ! = maxIndex & & points [ minIndex ] ! = points [ maxIndex ] ) {
relevantVertices . push_back ( points [ maxIndex ] ) ;
std : : vector < uint_fast64_t > maxVertexSet ( 1 , 1 ) ;
vertexSets . push_back ( std : : move ( maxVertexSet ) ) ;
@ -117,7 +117,7 @@ namespace storm {
vectorMatrix . col ( i ) < < points [ subset [ i ] ] , storm : : utility : : one < ValueType > ( ) ;
}
vectorMatrix . col ( subset . size ( ) ) < < points [ item ] , storm : : utility : : one < ValueType > ( ) ;
return ( vectorMatrix . fullPivLu ( ) . rank ( ) > subset . size ( ) ) ;
return ( vectorMatrix . fullPivLu ( ) . rank ( ) > ( int_fast64_t ) subset . size ( ) ) ;
}
template < typename ValueType >
@ -129,13 +129,13 @@ namespace storm {
const uint_fast64_t dimension = points . front ( ) . rows ( ) ;
EigenVector refPoint = points . front ( ) ;
EigenMatrix constraints ( points . size ( ) - 1 , dimension ) ;
for ( unsigned row = 1 ; row < points . size ( ) ; + + row ) {
for ( unsigned row = 1 ; row < points . size ( ) ; + + row ) {
constraints . row ( row - 1 ) = points [ row ] - refPoint ;
}
EigenVector normal = constraints . fullPivLu ( ) . kernel ( ) . col ( 0 ) ;
// Eigen returns the column vector 0...0 if the kernel is empty (i.e., there is no such hyperplane)
if ( normal . isZero ( ) ) {
if ( normal . isZero ( ) ) {
pointsAffineDependent = false ;
} else {
points . push_back ( refPoint + normal ) ;
@ -156,7 +156,7 @@ namespace storm {
newA . row ( row ) = resultMatrix . row ( row ) ;
newb ( row ) = resultVector ( row ) ;
}
for ( ; row < numOfRegularConstraints + numOfAdditionalConstraints ; + + row ) {
for ( ; row < numOfRegularConstraints + numOfAdditionalConstraints ; + + row ) {
newA . row ( row ) = std : : move ( additionalConstraints [ row - numOfRegularConstraints ] . first ) ;
newb ( row ) = additionalConstraints [ row - numOfRegularConstraints ] . second ;
}
@ -165,9 +165,9 @@ namespace storm {
// clear the additionally added points. Note that the order of the points might have changed
storm : : storage : : BitVector keptPoints ( points . size ( ) , true ) ;
for ( uint_fast64_t pointIndex = 0 ; pointIndex < points . size ( ) ; + + pointIndex ) {
for ( int_fast64_t row = 0 ; row < resultMatrix . rows ( ) ; + + row ) {
if ( ( resultMatrix . row ( row ) * points [ pointIndex ] ) ( 0 ) > resultVector ( row ) ) {
for ( uint_fast64_t pointIndex = 0 ; pointIndex < points . size ( ) ; + + pointIndex ) {
for ( int_fast64_t row = 0 ; row < resultMatrix . rows ( ) ; + + row ) {
if ( ( resultMatrix . row ( row ) * points [ pointIndex ] ) ( 0 ) > resultVector ( row ) ) {
keptPoints . set ( pointIndex , false ) ;
break ;
}
@ -175,11 +175,11 @@ namespace storm {
}
points = storm : : utility : : vector : : filterVector ( points , keptPoints ) ;
if ( generateRelevantVerticesAndVertexSets ) {
if ( generateRelevantVerticesAndVertexSets ) {
storm : : storage : : BitVector keptVertices ( relevantVertices . size ( ) , true ) ;
for ( uint_fast64_t vertexIndex = 0 ; vertexIndex < relevantVertices . size ( ) ; + + vertexIndex ) {
for ( int_fast64_t row = 0 ; row < resultMatrix . rows ( ) ; + + row ) {
if ( ( resultMatrix . row ( row ) * relevantVertices [ vertexIndex ] ) ( 0 ) > resultVector ( row ) ) {
for ( uint_fast64_t vertexIndex = 0 ; vertexIndex < relevantVertices . size ( ) ; + + vertexIndex ) {
for ( int_fast64_t row = 0 ; row < resultMatrix . rows ( ) ; + + row ) {
if ( ( resultMatrix . row ( row ) * relevantVertices [ vertexIndex ] ) ( 0 ) > resultVector ( row ) ) {
keptVertices . set ( vertexIndex , false ) ;
break ;
}
@ -195,22 +195,22 @@ namespace storm {
template < typename ValueType >
bool QuickHull < ValueType > : : findInitialVertices ( std : : vector < EigenVector > & points , std : : vector < uint_fast64_t > & verticesOfInitialPolytope ) const {
const uint_fast64_t dimension = points . front ( ) . rows ( ) ;
if ( points . size ( ) < dimension + 1 ) {
if ( points . size ( ) < dimension + 1 ) {
//not enough points to obtain a (non-degenerated) polytope
return false ;
}
// We first find some good candidates to get a (hopefully) large initial mesh.
storm : : storage : : BitVector notGoodCandidates ( points . size ( ) , true ) ;
for ( uint_fast64_t dim = 0 ; dim < dimension ; + + dim ) {
if ( ! notGoodCandidates . empty ( ) ) {
for ( uint_fast64_t dim = 0 ; dim < dimension ; + + dim ) {
if ( ! notGoodCandidates . empty ( ) ) {
uint_fast64_t minIndex = * notGoodCandidates . begin ( ) ;
uint_fast64_t maxIndex = * notGoodCandidates . begin ( ) ;
for ( uint_fast64_t pointIndex : notGoodCandidates ) {
if ( points [ minIndex ] ( dim ) > points [ pointIndex ] ( dim ) ) {
for ( uint_fast64_t pointIndex : notGoodCandidates ) {
if ( points [ minIndex ] ( dim ) > points [ pointIndex ] ( dim ) ) {
minIndex = pointIndex ;
}
if ( points [ maxIndex ] ( dim ) < points [ pointIndex ] ( dim ) ) {
if ( points [ maxIndex ] ( dim ) < points [ pointIndex ] ( dim ) ) {
maxIndex = pointIndex ;
}
}
@ -222,8 +222,8 @@ namespace storm {
//Found candidates. Now swap them to the front.
const uint_fast64_t numOfGoodCandidates = goodCandidates . getNumberOfSetBits ( ) ;
for ( auto const & goodCandidate : goodCandidates ) {
if ( goodCandidate > = numOfGoodCandidates ) {
for ( auto const & goodCandidate : goodCandidates ) {
if ( goodCandidate > = numOfGoodCandidates ) {
uint_fast64_t notGoodCandidate = * notGoodCandidates . begin ( ) ;
assert ( notGoodCandidate < numOfGoodCandidates ) ;
std : : swap ( points [ notGoodCandidate ] , points [ goodCandidate ] ) ;
@ -232,7 +232,7 @@ namespace storm {
}
storm : : storage : : geometry : : SubsetEnumerator < std : : vector < EigenVector > > subsetEnum ( points . size ( ) , dimension + 1 , points , affineFilter ) ;
if ( subsetEnum . setToFirstSubset ( ) ) {
if ( subsetEnum . setToFirstSubset ( ) ) {
verticesOfInitialPolytope = subsetEnum . getCurrentSubset ( ) ;
return true ;
} else {
@ -247,7 +247,7 @@ namespace storm {
std : : vector < Facet > result ;
result . reserve ( dimension + 1 ) ;
storm : : storage : : geometry : : SubsetEnumerator < > subsetEnum ( verticesOfInitialPolytope . size ( ) , dimension ) ;
if ( ! subsetEnum . setToFirstSubset ( ) ) {
if ( ! subsetEnum . setToFirstSubset ( ) ) {
STORM_LOG_THROW ( false , storm : : exceptions : : UnexpectedException , " Could not find an initial subset. " ) ;
}
do {
@ -255,13 +255,13 @@ namespace storm {
// set the points that lie on the new facet
std : : vector < uint_fast64_t > const & subset ( subsetEnum . getCurrentSubset ( ) ) ;
newFacet . points . reserve ( subset . size ( ) ) ;
for ( uint_fast64_t i : subset ) {
for ( uint_fast64_t i : subset ) {
newFacet . points . push_back ( verticesOfInitialPolytope [ i ] ) ;
}
//neighbors: these are always the remaining facets
newFacet . neighbors . reserve ( dimension ) ;
for ( uint_fast64_t i = 0 ; i < dimension + 1 ; + + i ) {
if ( i ! = result . size ( ) ) { //initFacets.size() will be the index of this new facet!
for ( uint_fast64_t i = 0 ; i < dimension + 1 ; + + i ) {
if ( i ! = result . size ( ) ) { //initFacets.size() will be the index of this new facet!
newFacet . neighbors . push_back ( i ) ;
}
}
@ -280,14 +280,14 @@ namespace storm {
assert ( facet . points . size ( ) = = dimension ) ;
EigenVector const & refPoint = points [ facet . points . back ( ) ] ;
EigenMatrix constraints ( dimension - 1 , dimension ) ;
for ( unsigned row = 0 ; row < dimension - 1 ; + + row ) {
for ( unsigned row = 0 ; row < dimension - 1 ; + + row ) {
constraints . row ( row ) = ( points [ facet . points [ row ] ] - refPoint ) ;
}
facet . normal = constraints . fullPivLu ( ) . kernel ( ) . col ( 0 ) ;
facet . offset = facet . normal . dot ( refPoint ) ;
// invert the plane if the insidePoint is not contained in it
if ( facet . normal . dot ( insidePoint ) > facet . offset ) {
if ( facet . normal . dot ( insidePoint ) > facet . offset ) {
facet . normal = - facet . normal ;
facet . offset = - facet . offset ;
}
@ -306,19 +306,19 @@ namespace storm {
computeOutsideSetOfFacet ( facets [ facetIndex ] , currentOutsidePoints , points ) ;
}
for ( uint_fast64_t facetCount = currentFacets . getNextSetIndex ( 0 ) ; facetCount ! = currentFacets . size ( ) ; facetCount = currentFacets . getNextSetIndex ( facetCount + 1 ) ) {
for ( uint_fast64_t facetCount = currentFacets . getNextSetIndex ( 0 ) ; facetCount ! = currentFacets . size ( ) ; facetCount = currentFacets . getNextSetIndex ( facetCount + 1 ) ) {
// set all points to false to get rid of points that lie within the polytope after each iteration
currentOutsidePoints . clear ( ) ;
// Find a facet with a non-empty outside set
if ( ! facets [ facetCount ] . outsideSet . empty ( ) ) {
if ( ! facets [ facetCount ] . outsideSet . empty ( ) ) {
uint_fast64_t numberOfNewFacets = 0 ;
// Now we compute the enlarged mesh
uint_fast64_t farAwayPointIndex = facets [ facetCount ] . maxOutsidePointIndex ;
// get Visible set from maxOutsidePoint of the current facet
std : : set < uint_fast64_t > visibleSet = getVisibleSet ( facets , facetCount , points [ farAwayPointIndex ] ) ;
std : : set < uint_fast64_t > invisibleSet = getInvisibleNeighbors ( facets , visibleSet ) ;
for ( auto invisFacetIt = invisibleSet . begin ( ) ; invisFacetIt ! = invisibleSet . end ( ) ; + + invisFacetIt ) {
for ( auto visFacetIt = visibleSet . begin ( ) ; visFacetIt ! = visibleSet . end ( ) ; + + visFacetIt ) {
for ( auto invisFacetIt = invisibleSet . begin ( ) ; invisFacetIt ! = invisibleSet . end ( ) ; + + invisFacetIt ) {
for ( auto visFacetIt = visibleSet . begin ( ) ; visFacetIt ! = visibleSet . end ( ) ; + + visFacetIt ) {
if ( std : : find ( facets [ * invisFacetIt ] . neighbors . begin ( ) , facets [ * invisFacetIt ] . neighbors . end ( ) , * visFacetIt ) ! = facets [ * invisFacetIt ] . neighbors . end ( ) ) {
Facet newFacet ;
// Set points of Facet
@ -339,14 +339,14 @@ namespace storm {
}
}
for ( auto visibleFacet : visibleSet ) {
for ( uint_fast64_t outsidePoint : facets [ visibleFacet ] . outsideSet ) {
for ( auto visibleFacet : visibleSet ) {
for ( uint_fast64_t outsidePoint : facets [ visibleFacet ] . outsideSet ) {
currentOutsidePoints . set ( outsidePoint , true ) ;
}
currentFacets . set ( visibleFacet , false ) ;
}
// compute new outside sets
for ( uint_fast64_t facetIndex : currentFacets ) {
for ( uint_fast64_t facetIndex : currentFacets ) {
computeOutsideSetOfFacet ( facets [ facetIndex ] , currentOutsidePoints , points ) ;
}
@ -360,18 +360,18 @@ namespace storm {
void QuickHull < ValueType > : : getPolytopeFromMesh ( std : : vector < EigenVector > const & points , std : : vector < Facet > const & facets , storm : : storage : : BitVector const & currentFacets , bool generateRelevantVerticesAndVertexSets ) {
storm : : storage : : geometry : : HyperplaneCollector < ValueType > hyperplaneCollector ;
for ( auto facet : currentFacets ) {
for ( auto facet : currentFacets ) {
hyperplaneCollector . insert ( std : : move ( facets [ facet ] . normal ) , std : : move ( facets [ facet ] . offset ) , generateRelevantVerticesAndVertexSets ? & facets [ facet ] . points : nullptr ) ;
}
if ( generateRelevantVerticesAndVertexSets ) {
if ( generateRelevantVerticesAndVertexSets ) {
//Get the mapping from a hyperplane to the set of vertices that lie on that plane, erase the duplicates, and count for each vertex the number of hyperplanes on which that vertex lies
vertexSets = hyperplaneCollector . getIndexLists ( ) ;
std : : vector < uint_fast64_t > hyperplanesOnVertexCounter ( points . size ( ) , 0 ) ;
for ( auto & vertexVector : vertexSets ) {
for ( auto & vertexVector : vertexSets ) {
std : : set < uint_fast64_t > vertexSet ;
for ( auto const & i : vertexVector ) {
if ( vertexSet . insert ( i ) . second ) {
for ( auto const & i : vertexVector ) {
if ( vertexSet . insert ( i ) . second ) {
+ + hyperplanesOnVertexCounter [ i ] ;
}
}
@ -382,8 +382,8 @@ namespace storm {
//Therefore, we additionally store the old indices for every vertex to be able to translate from old to new indices
std : : unordered_map < EigenVector , std : : vector < uint_fast64_t > > relevantVerticesMap ;
relevantVerticesMap . reserve ( points . size ( ) ) ;
for ( uint_fast64_t vertexIndex = 0 ; vertexIndex < hyperplanesOnVertexCounter . size ( ) ; + + vertexIndex ) {
if ( hyperplanesOnVertexCounter [ vertexIndex ] > = points . front ( ) . rows ( ) ) {
for ( uint_fast64_t vertexIndex = 0 ; vertexIndex < hyperplanesOnVertexCounter . size ( ) ; + + vertexIndex ) {
if ( ( int_fast64_t ) hyperplanesOnVertexCounter [ vertexIndex ] > = points . front ( ) . rows ( ) ) {
auto mapEntry = relevantVerticesMap . insert ( typename std : : unordered_map < EigenVector , std : : vector < uint_fast64_t > > : : value_type ( points [ vertexIndex ] , std : : vector < uint_fast64_t > ( ) ) ) . first ;
mapEntry - > second . push_back ( vertexIndex ) ;
}
@ -392,17 +392,17 @@ namespace storm {
std : : vector < uint_fast64_t > oldToNewIndexMapping ( points . size ( ) , points . size ( ) ) ; //Initialize with some illegal value
relevantVertices . clear ( ) ;
relevantVertices . reserve ( relevantVerticesMap . size ( ) ) ;
for ( auto const & mapEntry : relevantVerticesMap ) {
for ( auto const & oldIndex : mapEntry . second ) {
for ( auto const & mapEntry : relevantVerticesMap ) {
for ( auto const & oldIndex : mapEntry . second ) {
oldToNewIndexMapping [ oldIndex ] = relevantVertices . size ( ) ;
}
relevantVertices . push_back ( mapEntry . first ) ;
}
//Actually translate and erase duplicates
for ( auto & vertexVector : vertexSets ) {
for ( auto & vertexVector : vertexSets ) {
std : : set < uint_fast64_t > vertexSet ;
for ( auto const & oldIndex : vertexVector ) {
if ( hyperplanesOnVertexCounter [ oldIndex ] > = points . front ( ) . rows ( ) ) {
for ( auto const & oldIndex : vertexVector ) {
if ( ( int_fast64_t ) hyperplanesOnVertexCounter [ oldIndex ] > = points . front ( ) . rows ( ) ) {
vertexSet . insert ( oldToNewIndexMapping [ oldIndex ] ) ;
}
}
@ -424,14 +424,14 @@ namespace storm {
std : : set < uint_fast64_t > visibleSet ;
facetsChecked . insert ( startIndex ) ;
visibleSet . insert ( startIndex ) ;
for ( uint_fast64_t i = 0 ; i < facets [ startIndex ] . neighbors . size ( ) ; + + i ) {
for ( uint_fast64_t i = 0 ; i < facets [ startIndex ] . neighbors . size ( ) ; + + i ) {
facetsToCheck . insert ( facets [ startIndex ] . neighbors [ i ] ) ;
}
while ( ! facetsToCheck . empty ( ) ) {
auto elementIt = facetsToCheck . begin ( ) ;
if ( point . dot ( facets [ * elementIt ] . normal ) > facets [ * elementIt ] . offset ) {
visibleSet . insert ( * elementIt ) ;
for ( uint_fast64_t i = 0 ; i < facets [ * elementIt ] . neighbors . size ( ) ; + + i ) {
for ( uint_fast64_t i = 0 ; i < facets [ * elementIt ] . neighbors . size ( ) ; + + i ) {
if ( facetsChecked . find ( facets [ * elementIt ] . neighbors [ i ] ) = = facetsChecked . end ( ) ) {
facetsToCheck . insert ( facets [ * elementIt ] . neighbors [ i ] ) ;
}
@ -445,8 +445,8 @@ namespace storm {
template < typename ValueType >
void QuickHull < ValueType > : : setNeighborhoodOfNewFacets ( std : : vector < Facet > & facets , uint_fast64_t firstNewFacet , uint_fast64_t dimension ) const {
for ( uint_fast64_t currentFacet = firstNewFacet ; currentFacet < facets . size ( ) ; + + currentFacet ) {
for ( uint_fast64_t otherFacet = currentFacet + 1 ; otherFacet < facets . size ( ) ; + + otherFacet ) {
for ( uint_fast64_t currentFacet = firstNewFacet ; currentFacet < facets . size ( ) ; + + currentFacet ) {
for ( uint_fast64_t otherFacet = currentFacet + 1 ; otherFacet < facets . size ( ) ; + + otherFacet ) {
if ( getCommonPoints ( facets [ currentFacet ] , facets [ otherFacet ] ) . size ( ) > = dimension - 1 ) {
facets [ currentFacet ] . neighbors . push_back ( otherFacet ) ;
facets [ otherFacet ] . neighbors . push_back ( currentFacet ) ;
@ -466,9 +466,9 @@ namespace storm {
template < typename ValueType >
void QuickHull < ValueType > : : computeOutsideSetOfFacet ( Facet & facet , storm : : storage : : BitVector & currentOutsidePoints , std : : vector < EigenVector > const & points ) const {
ValueType maxMultiplicationResult = facet . offset ;
for ( uint_fast64_t pointIndex : currentOutsidePoints ) {
for ( uint_fast64_t pointIndex : currentOutsidePoints ) {
ValueType multiplicationResult = points [ pointIndex ] . dot ( facet . normal ) ;
if ( multiplicationResult > facet . offset ) {
if ( multiplicationResult > facet . offset ) {
currentOutsidePoints . set ( pointIndex , false ) ; // we already know that the point lies outside so it can be ignored for future facets
facet . outsideSet . push_back ( pointIndex ) ;
if ( multiplicationResult > maxMultiplicationResult ) {
@ -482,8 +482,8 @@ namespace storm {
template < typename ValueType >
std : : vector < uint_fast64_t > QuickHull < ValueType > : : getCommonPoints ( Facet const & lhs , Facet const & rhs ) const {
std : : vector < uint_fast64_t > commonPoints ;
for ( uint_fast64_t lhsPoint : lhs . points ) {
for ( uint_fast64_t rhsPoint : rhs . points ) {
for ( uint_fast64_t lhsPoint : lhs . points ) {
for ( uint_fast64_t rhsPoint : rhs . points ) {
if ( lhsPoint = = rhsPoint ) {
commonPoints . push_back ( lhsPoint ) ;
}
@ -495,8 +495,8 @@ namespace storm {
template < typename ValueType >
std : : set < uint_fast64_t > QuickHull < ValueType > : : getInvisibleNeighbors ( std : : vector < Facet > & facets , std : : set < uint_fast64_t > const & visibleSet ) const {
std : : set < uint_fast64_t > invisibleNeighbors ;
for ( auto currentFacetIt = visibleSet . begin ( ) ; currentFacetIt ! = visibleSet . end ( ) ; + + currentFacetIt ) {
for ( uint_fast64_t currentNeighbor = 0 ; currentNeighbor < facets [ * currentFacetIt ] . neighbors . size ( ) ; + + currentNeighbor ) {
for ( auto currentFacetIt = visibleSet . begin ( ) ; currentFacetIt ! = visibleSet . end ( ) ; + + currentFacetIt ) {
for ( uint_fast64_t currentNeighbor = 0 ; currentNeighbor < facets [ * currentFacetIt ] . neighbors . size ( ) ; + + currentNeighbor ) {
if ( visibleSet . find ( facets [ * currentFacetIt ] . neighbors [ currentNeighbor ] ) = = visibleSet . end ( ) ) {
invisibleNeighbors . insert ( facets [ * currentFacetIt ] . neighbors [ currentNeighbor ] ) ;
}
TimQu
8 years ago
xxxxxxxxxx