@ -29,27 +29,44 @@ namespace storm {
template < typename ModelType >
QuantileHelper < ModelType > : : QuantileHelper ( ModelType const & model , storm : : logic : : QuantileFormula const & quantileFormula ) : model ( model ) , quantileFormula ( quantileFormula ) {
// Intentionally left empty
// Do all kinds of sanity check.
std : : set < storm : : expressions : : Variable > quantileVariables ;
for ( auto const & quantileVariable : quantileFormula . getBoundVariables ( ) ) {
STORM_LOG_THROW ( quantileVariables . count ( quantileVariable . second ) = = 0 , storm : : exceptions : : NotSupportedException , " Quantile formula considers the same bound variable twice. " ) ;
quantileVariables . insert ( quantileVariable . second ) ;
}
STORM_LOG_THROW ( quantileFormula . getSubformula ( ) . isProbabilityOperatorFormula ( ) , storm : : exceptions : : NotSupportedException , " Quantile formula needs probability operator inside. The formula " < < quantileFormula < < " is not supported. " ) ;
auto const & probOpFormula = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) ;
STORM_LOG_THROW ( probOpFormula . getBound ( ) . comparisonType = = storm : : logic : : ComparisonType : : Greater | | probOpFormula . getBound ( ) . comparisonType = = storm : : logic : : ComparisonType : : LessEqual , storm : : exceptions : : NotSupportedException , " Probability operator inside quantile formula needs to have bound > or <=. " ) ;
STORM_LOG_THROW ( probOpFormula . hasBound ( ) , storm : : exceptions : : InvalidOperationException , " Probability operator inside quantile formula needs to have a bound. " ) ;
STORM_LOG_THROW ( ! model . isNondeterministicModel ( ) | | probOpFormula . hasOptimalityType ( ) , storm : : exceptions : : InvalidOperationException , " Probability operator inside quantile formula needs to have an optimality type. " ) ;
STORM_LOG_WARN_COND ( probOpFormula . getBound ( ) . comparisonType = = storm : : logic : : ComparisonType : : Greater | | probOpFormula . getBound ( ) . comparisonType = = storm : : logic : : ComparisonType : : LessEqual , " Probability operator inside quantile formula needs to have bound > or <=. The specified comparison type might lead to non-termination. " ) ; // This has to do with letting bound variables approach infinity, e.g., Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
bool lowerBounded = storm : : logic : : isLowerBound ( probOpFormula . getBound ( ) . comparisonType ) ;
STORM_LOG_THROW ( probOpFormula . getSubformula ( ) . isBoundedUntilFormula ( ) , storm : : exceptions : : NotSupportedException , " Quantile formula needs bounded until probability operator formula as subformula. The formula " < < quantileFormula < < " is not supported. " ) ;
auto const & boundedUntilFormula = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getSubformula ( ) . asBoundedUntilFormula ( ) ;
auto const & boundedUntilFormula = probOpFormula . getSubformula ( ) . asBoundedUntilFormula ( ) ;
std : : set < storm : : expressions : : Variable > boundVariables ;
for ( uint64_t dim = 0 ; dim < boundedUntilFormula . getDimension ( ) ; + + dim ) {
storm : : expressions : : Expression boundExpression ;
if ( boundedUntilFormula . hasUpperBound ( dim ) ) {
STORM_LOG_THROW ( ! boundedUntilFormula . hasLowerBound ( dim ) , storm : : exceptions : : NotSupportedException , " Interval bounds are not supported within quantile formulas. " ) ;
STORM_LOG_THROW ( ! boundedUntilFormula . isUpperBoundStrict ( dim ) , storm : : exceptions : : NotSupportedException , " Only non-strict upper reward bounds are supported for quantiles. " ) ;
boundExpression = boundedUntilFormula . getUpperBound ( dim ) ;
} else if ( boundedUntilFormula . hasLowerBound ( dim ) ) {
STORM_LOG_THROW ( ! boundedUntilFormula . isLowerBoundStrict ( dim ) , storm : : exceptions : : NotSupportedException , " Only non-strict lower reward bounds are supported for quantiles. " ) ;
boundExpression = boundedUntilFormula . getLowerBound ( dim ) ;
}
if ( boundExpression . isInitialized ( ) & & boundExpression . containsVariables ( ) ) {
STORM_LOG_THROW ( boundExpression . isVariable ( ) , storm : : exceptions : : NotSupportedException , " Non-trivial bound expressions such as ' " < < boundExpression < < " ' are not supported. Either specify a constant or a quantile variable. " ) ;
storm : : expressions : : Variable const & boundVariable = boundExpression . getBaseExpression ( ) . asVariableExpression ( ) . getVariable ( ) ;
STORM_LOG_THROW ( boundVariables . count ( boundVariable ) = = 0 , storm : : exceptions : : NotSupportedException , " Variable " < < boundExpression < < " occurs at multiple reward bounds. " ) ;
boundVariables . insert ( boundVariable ) ;
STORM_LOG_THROW ( quantileVariables . count ( boundVariable ) = = 1 , storm : : exceptions : : NotSupportedException , " The formula contains undefined constant ' " < < boundExpression < < " '. " ) ;
}
}
// Only > and <= are supported for upper bounds. This is to make sure that Pr>0.7 [F "goal"] holds iff Pr>0.7 [F<=B "goal"] holds for some B.
// Only >= and < are supported for lower bounds. (EC construction..)
// TODO
// Fix precision increasing
// Multiple quantile formulas in the same file yield constants def clash
// Bounds are either constants or variables that are declared in the quantile formula.
// Prop op has optimality type
// No Prmin with lower cost bounds: Ec construction fails. In single obj we would do 1-Prmax[F "nonRewOrNonGoalEC"] but this invalidates other lower/upper cost bounds.
/* Todo: Current assumptions:
* Subformula is always prob operator with bounded until
* Each bound variable occurs at most once ( change this ? )
* cost bounds are assumed to be always non - strict ( the epochs returned by the reward unfolding assume strict lower and non - strict upper cost bounds , though . . . )
* ' reasonable ' quantile ( e . g . not quantile ( max B , Pmax > 0.5 [ F < = B G ] ) ( we could also filter out these things early on )
*/
// ignore optimization direction for quantile variables
}
enum class BoundTransformation {
@ -58,7 +75,7 @@ namespace storm {
GreaterEqualZero ,
LessEqualZero
} ;
std : : shared_ptr < storm : : logic : : ProbabilityOperatorFormula > transformBoundedUntilOperator ( storm : : logic : : ProbabilityOperatorFormula const & boundedUntilOperator , std : : vector < BoundTransformation > const & transformations ) {
std : : shared_ptr < storm : : logic : : ProbabilityOperatorFormula > transformBoundedUntilOperator ( storm : : logic : : ProbabilityOperatorFormula const & boundedUntilOperator , std : : vector < BoundTransformation > const & transformations , bool complementQuery = false ) {
auto const & origBoundedUntil = boundedUntilOperator . getSubformula ( ) . asBoundedUntilFormula ( ) ;
STORM_LOG_ASSERT ( transformations . size ( ) = = origBoundedUntil . getDimension ( ) , " Tried to replace the bound of a dimension that is higher than the number of dimensions of the formula. " ) ;
std : : vector < std : : shared_ptr < storm : : logic : : Formula const > > leftSubformulas , rightSubformulas ;
@ -107,7 +124,11 @@ namespace storm {
} else {
newBoundedUntil = std : : make_shared < storm : : logic : : BoundedUntilFormula > ( origBoundedUntil . getLeftSubformula ( ) . asSharedPointer ( ) , origBoundedUntil . getRightSubformula ( ) . asSharedPointer ( ) , lowerBounds , upperBounds , timeBoundReferences ) ;
}
return std : : make_shared < storm : : logic : : ProbabilityOperatorFormula > ( newBoundedUntil , boundedUntilOperator . getOperatorInformation ( ) ) ;
storm : : logic : : OperatorInformation newOpInfo ( boundedUntilOperator . getOperatorInformation ( ) . optimalityType , boundedUntilOperator . getBound ( ) ) ;
if ( complementQuery ) {
newOpInfo . bound - > comparisonType = storm : : logic : : invert ( newOpInfo . bound - > comparisonType ) ;
}
return std : : make_shared < storm : : logic : : ProbabilityOperatorFormula > ( newBoundedUntil , newOpInfo ) ;
}
/// Increases the precision of solver results
@ -174,55 +195,26 @@ namespace storm {
}
template < typename ModelType >
storm : : storage : : BitVector QuantileHelper < ModelType > : : getDimensionsForVariable ( storm : : expressions : : Variable const & var ) const {
auto const & boundedUntil = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getSubformula ( ) . asBoundedUntilFormula ( ) ;
storm : : storage : : BitVector result ( boundedUntil . getDimension ( ) , false ) ;
for ( uint64_t dim = 0 ; dim < boundedUntil . getDimension ( ) ; + + dim ) {
if ( boundedUntil . hasLowerBound ( dim ) & & boundedUntil . getLowerBound ( dim ) . isVariable ( ) & & boundedUntil . getLowerBound ( dim ) . getBaseExpression ( ) . asVariableExpression ( ) . getVariable ( ) = = var ) {
result . set ( dim , true ) ;
}
if ( boundedUntil . hasUpperBound ( dim ) & & boundedUntil . getUpperBound ( dim ) . isVariable ( ) & & boundedUntil . getUpperBound ( dim ) . getBaseExpression ( ) . asVariableExpression ( ) . getVariable ( ) = = var ) {
result . set ( dim , true ) ;
}
}
return result ;
}
std : : vector < std : : vector < typename ModelType : : ValueType > > QuantileHelper < ModelType > : : computeQuantile ( Environment const & env ) {
numCheckedEpochs = 0 ;
numPrecisionRefinements = 0 ;
cachedSubQueryResults . clear ( ) ;
template < typename ModelType >
std : : vector < std : : vector < typename ModelType : : ValueType > > QuantileHelper < ModelType > : : computeQuantile ( Environment const & env ) const {
std : : vector < std : : vector < ValueType > > result ;
Environment envCpy = env ; // It might be necessary to increase the precision during the computation
numCheckedEpochs = 0 ; numPrecisionRefinements = 0 ;
if ( getOpenDimensions ( ) . getNumberOfSetBits ( ) = = 1 ) {
uint64_t dimension = * getOpenDimensions ( ) . begin ( ) ;
// Call the internal recursive function
auto internalResult = computeQuantile ( envCpy , getOpenDimensions ( ) , false ) ;
bool zeroSatisfiesFormula = checkLimitValue ( envCpy , storm : : storage : : BitVector ( getDimension ( ) , false ) ) ;
bool infSatisfiesFormula = checkLimitValue ( envCpy , getOpenDimensions ( ) ) ;
if ( zeroSatisfiesFormula ! = infSatisfiesFormula ) {
while ( true ) {
auto quantileRes = computeQuantileForDimension ( envCpy , dimension ) ;
if ( quantileRes ) {
result = { { storm : : utility : : convertNumber < ValueType > ( quantileRes - > first ) * quantileRes - > second } } ;
break ;
}
increasePrecision ( envCpy ) ;
+ + numPrecisionRefinements ;
}
} else if ( zeroSatisfiesFormula ) {
// thus also infSatisfiesFormula is true, i.e., all bound values satisfy the formula
if ( storm : : solver : : minimize ( getOptimizationDirForDimension ( dimension ) ) ) {
result = { { storm : : utility : : zero < ValueType > ( ) } } ;
} else {
result = { { storm : : utility : : infinity < ValueType > ( ) } } ;
}
} else {
// i.e., no bound value satisfies the formula
result = { { } } ;
}
} else if ( getOpenDimensions ( ) . getNumberOfSetBits ( ) = = 2 ) {
result = computeTwoDimensionalQuantile ( envCpy ) ;
// Translate the result by applying the scaling factors.
for ( auto const & costLimits : internalResult . first . getGenerator ( ) ) {
std : : vector < ValueType > resultPoint ( costLimits . size ( ) ) ;
storm : : utility : : vector : : applyPointwise < CostLimit , ValueType , ValueType > ( costLimits , internalResult . second , resultPoint , [ ] ( CostLimit const & costLimit , ValueType const & factor ) - > ValueType {
if ( costLimit . isInfinity ( ) ) {
return storm : : utility : : infinity < ValueType > ( ) ;
} else {
STORM_LOG_THROW ( false , storm : : exceptions : : NotSupportedException , " The quantile formula considers " < < getOpenDimensions ( ) . getNumberOfSetBits ( ) < < " open dimensions. Only one- or two-dimensional quantiles are supported. " ) ;
return storm : : utility : : convertNumber < ValueType > ( costLimit . get ( ) ) * factor ;
} } ) ;
result . push_back ( resultPoint ) ;
}
if ( storm : : settings : : getModule < storm : : settings : : modules : : CoreSettings > ( ) . isShowStatisticsSet ( ) ) {
std : : cout < < " Number of checked epochs: " < < numCheckedEpochs < < std : : endl ;
@ -231,20 +223,26 @@ namespace storm {
return result ;
}
template < typename ModelType >
std : : pair < CostLimitClosure , std : : vector < typename QuantileHelper < ModelType > : : ValueType > > QuantileHelper < ModelType > : : computeQuantile ( Environment & env , storm : : storage : : BitVector const & consideredDimensions ) const {
// Todo: ask the cache first
std : : pair < CostLimitClosure , std : : vector < typename QuantileHelper < ModelType > : : ValueType > > QuantileHelper < ModelType > : : computeQuantile ( Environment & env , storm : : storage : : BitVector const & consideredDimensions , bool complementaryQuery ) {
STORM_LOG_ASSERT ( consideredDimensions . isSubsetOf ( getOpenDimensions ( ) ) , " Considered dimensions for a quantile query should be a subset of the set of dimensions without a fixed bound. " ) ;
storm : : storage : : BitVector cacheKey = consideredDimensions ;
cacheKey . resize ( cacheKey . size ( ) + 1 , complementaryQuery ) ;
auto cacheIt = cachedSubQueryResults . find ( cacheKey ) ;
if ( cacheIt ! = cachedSubQueryResults . end ( ) ) {
return cacheIt - > second ;
}
auto boundedUntilOp = transformBoundedUntilOperator ( quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) , std : : vector < BoundTransformation > ( getDimension ( ) , BoundTransformation : : None ) ) ;
auto boundedUntilOp = transformBoundedUntilOperator ( quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) , std : : vector < BoundTransformation > ( getDimension ( ) , BoundTransformation : : None ) , complementaryQuery ) ;
std : : set < storm : : expressions : : Variable > infinityVariables ;
storm : : storage : : BitVector lowerBoundedDimensions ( getDimension ( ) ) ;
storm : : storage : : BitVector downwardClosedDimensions ( getDimension ( ) ) ;
bool hasLowerValueBound = storm : : logic : : isLowerBound ( boundedUntilOp - > getComparisonType ( ) ) ;
for ( auto const & d : getOpenDimensions ( ) ) {
if ( consideredDimensions . get ( d ) ) {
bool hasLowerCostBound = boundedUntilOp - > getSubformula ( ) . asBoundedUntilFormula ( ) . hasLowerBound ( d ) ;
lowerBoundedDimensions . set ( d , hasLowerCostBound ) ;
bool hasLowerValueBound = storm : : logic : : isLowerBound ( boundedUntilOp - > getComparisonType ( ) ) ;
downwardClosedDimensions . set ( d , hasLowerCostBound = = hasLowerValueBound ) ;
} else {
infinityVariables . insert ( getVariableForDimension ( d ) ) ;
@ -253,23 +251,57 @@ namespace storm {
downwardClosedDimensions = downwardClosedDimensions % consideredDimensions ;
CostLimitClosure satCostLimits ( downwardClosedDimensions ) , unsatCostLimits ( ~ downwardClosedDimensions ) ;
// Initialize the (un)sat cost limits to guarantee termination
bool onlyUpperCostBounds = lowerBoundedDimensions . empty ( ) ;
bool onlyLowerCostBounds = lowerBoundedDimensions = = consideredDimensions ;
if ( onlyUpperCostBounds | | onlyLowerCostBounds ) {
for ( auto const & k : consideredDimensions ) {
storm : : storage : : BitVector subQueryDimensions = consideredDimensions ;
subQueryDimensions . set ( k , false ) ;
bool subQueryComplement = complementaryQuery ! = ( ( onlyUpperCostBounds & & hasLowerValueBound ) | | ( onlyLowerCostBounds & & ! hasLowerValueBound ) ) ;
auto subQueryResult = computeQuantile ( env , subQueryDimensions , subQueryComplement ) ;
for ( auto const & subQueryCostLimit : subQueryResult . first . getGenerator ( ) ) {
CostLimits initPoint ;
uint64_t i = 0 ;
for ( auto const & dim : consideredDimensions ) {
if ( dim = = k ) {
initPoint . push_back ( CostLimit : : infinity ( ) ) ;
} else {
initPoint . push_back ( subQueryCostLimit [ i ] ) ;
+ + i ;
}
}
if ( subQueryComplement ) {
unsatCostLimits . insert ( initPoint ) ;
} else {
satCostLimits . insert ( initPoint ) ;
}
}
}
} else {
STORM_LOG_WARN ( " Quantile formula considers mixtures of upper and lower reward-bounds. Termination is not guaranteed. " ) ;
}
// Loop until the goal precision is reached.
STORM_LOG_DEBUG ( " Computing quantile for dimensions: " < < consideredDimensions ) ;
while ( true ) {
// initialize Reward unfolding and data that will be needed for each epoch
// initialize r eward unfolding and data that will be needed for each epoch
MultiDimensionalRewardUnfolding < ValueType , true > rewardUnfolding ( model , boundedUntilOp , infinityVariables ) ;
if ( computeQuantile ( env , consideredDimensions , lowerBoundedDimensions , satCostLimits , unsatCostLimits , rewardUnfolding ) ) {
if ( computeQuantile ( env , consideredDimensions , * boundedUntilOp , lowerBoundedDimensions , satCostLimits , unsatCostLimits , rewardUnfolding ) ) {
std : : vector < ValueType > scalingFactors ;
for ( auto const & dim : consideredDimensions ) {
scalingFactors . push_back ( rewardUnfolding . getDimension ( dim ) . scalingFactor ) ;
}
return { satCostLimits , scalingFactors } ;
std : : pair < CostLimitClosure , std : : vector < ValueType > > result ( satCostLimits , scalingFactors ) ;
cachedSubQueryResults . emplace ( cacheKey , result ) ;
return result ;
}
STORM_LOG_WARN ( " Restarting quantile computation due to insufficient precision. " ) ;
+ + numPrecisionRefinements ;
increasePrecision ( env ) ;
}
}
bool getNextCandidateCostLimit ( CostLimit const & maxCostLimit , CostLimits & current ) {
if ( maxCostLimit . get ( ) = = 0 ) {
return false ;
@ -296,12 +328,6 @@ namespace storm {
return true ;
}
bool unionContainsAllCostLimits ( CostLimitClosure const & lhs , CostLimitClosure const & rhs ) {
// TODO
assert ( false ) ;
}
bool translateEpochToCostLimits ( EpochManager : : Epoch const & epoch , EpochManager : : Epoch const & startEpoch , storm : : storage : : BitVector const & consideredDimensions , storm : : storage : : BitVector const & lowerBoundedDimensions , EpochManager const & epochManager , CostLimits & epochAsCostLimits ) {
for ( uint64_t dim = 0 ; dim < consideredDimensions . size ( ) ; + + dim ) {
if ( consideredDimensions . get ( dim ) ) {
@ -331,9 +357,8 @@ namespace storm {
return true ;
}
template < typename ModelType >
bool QuantileHelper < ModelType > : : computeQuantile ( Environment & env , storm : : storage : : BitVector const & consideredDimensions , storm : : storage : : BitVector const & lowerBoundedDimensions , CostLimitClosure & satCostLimits , CostLimitClosure & unsatCostLimits , MultiDimensionalRewardUnfolding < ValueType , true > & rewardUnfolding ) const {
bool QuantileHelper < ModelType > : : computeQuantile ( Environment & env , storm : : storage : : BitVector const & consideredDimensions , storm : : logic : : ProbabilityOperatorFormula const & boundedUntilOperator , storm : : storage : : BitVector const & lowerBoundedDimensions , CostLimitClosure & satCostLimits , CostLimitClosure & unsatCostLimits , MultiDimensionalRewardUnfolding < ValueType , true > & rewardUnfolding ) {
auto lowerBound = rewardUnfolding . getLowerObjectiveBound ( ) ;
auto upperBound = rewardUnfolding . getUpperObjectiveBound ( ) ;
@ -341,12 +366,14 @@ namespace storm {
std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > minMaxSolver ;
std : : set < EpochManager : : Epoch > checkedEpochs ;
CostLimit currentMaxCostLimit ( 0 ) ;
for ( CostLimit currentMaxCostLimit ( 0 ) ; ! unionContainsAllCostLimits ( satCostLimits , unsatCostLimits ) ; + + currentMaxCostLimit . get ( ) ) {
bool progress = true ;
for ( CostLimit currentMaxCostLimit ( 0 ) ; progress ; + + currentMaxCostLimit . get ( ) ) {
CostLimits currentCandidate ( satCostLimits . dimension ( ) , currentMaxCostLimit ) ;
// We can only stop the exploration, if the upward closure of the point in the 'top right corner' is contained in the (un)satCostlLimits.
progress = ! satCostLimits . containsUpwardClosure ( currentCandidate ) & & ! unsatCostLimits . containsUpwardClosure ( currentCandidate ) ;
do {
if ( ! satCostLimits . contains ( currentCandidate ) & & ! unsatCostLimits . contains ( currentCandidate ) ) {
progress = true ;
// Transform candidate cost limits to an appropriate start epoch
auto startEpoch = rewardUnfolding . getStartEpoch ( true ) ;
auto costLimitIt = currentCandidate . begin ( ) ;
@ -367,7 +394,7 @@ namespace storm {
if ( checkedEpochs . count ( epoch ) = = 0 ) {
checkedEpochs . insert ( epoch ) ;
auto & epochModel = rewardUnfolding . setCurrentEpoch ( epoch ) ;
rewardUnfolding . setSolutionForCurrentEpoch ( epochModel . analyzeSingleObjective ( env , quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getOptimalityType ( ) , x , b , minMaxSolver , lowerBound , upperBound ) ) ;
rewardUnfolding . setSolutionForCurrentEpoch ( epochModel . analyzeSingleObjective ( env , boundedUntilOperator . getOptimalityType ( ) , x , b , minMaxSolver , lowerBound , upperBound ) ) ;
+ + numCheckedEpochs ;
CostLimits epochAsCostLimits ;
@ -376,15 +403,14 @@ namespace storm {
bool propertySatisfied ;
if ( env . solver ( ) . isForceSoundness ( ) ) {
auto lowerUpperValue = getLowerUpperBound ( env , currValue ) ;
propertySatisfied = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getBound ( ) . isSatisfied ( lowerUpperValue . first ) ;
if ( propertySatisfied ! = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getBound ( ) . isSatisfied ( lowerUpperValue . second ) ) {
propertySatisfied = boundedUntilOperator . getBound ( ) . isSatisfied ( lowerUpperValue . first ) ;
if ( propertySatisfied ! = boundedUntilOperator . getBound ( ) . isSatisfied ( lowerUpperValue . second ) ) {
// unclear result due to insufficient precision.
return false ;
}
} else {
propertySatisfied = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getBound ( ) . isSatisfied ( currValue ) ;
propertySatisfied = boundedUntilOperator . getBound ( ) . isSatisfied ( currValue ) ;
}
if ( propertySatisfied ) {
satCostLimits . insert ( epochAsCostLimits ) ;
} else {
@ -400,382 +426,6 @@ namespace storm {
}
///////
template < typename ValueType >
void filterDominatedPoints ( std : : vector < std : : vector < ValueType > > & points , std : : vector < storm : : solver : : OptimizationDirection > const & dirs ) {
std : : vector < std : : vector < ValueType > > result ;
// Note: this is slow and not inplace but also most likely not performance critical
for ( auto const & p1 : points ) {
bool p1Dominated = false ;
for ( auto const & p2 : points ) {
assert ( p1 . size ( ) = = p2 . size ( ) ) ;
bool p2DominatesP1 = false ;
for ( uint64_t i = 0 ; i < dirs . size ( ) ; + + i ) {
if ( storm : : solver : : minimize ( dirs [ i ] ) ? p2 [ i ] < = p1 [ i ] : p2 [ i ] > = p1 [ i ] ) {
if ( p2 [ i ] ! = p1 [ i ] ) {
p2DominatesP1 = true ;
}
} else {
p2DominatesP1 = false ;
break ;
}
}
if ( p2DominatesP1 ) {
p1Dominated = true ;
break ;
}
}
if ( ! p1Dominated ) {
result . push_back ( p1 ) ;
}
}
points = std : : move ( result ) ;
}
template < typename ModelType >
std : : vector < std : : vector < typename ModelType : : ValueType > > QuantileHelper < ModelType > : : computeTwoDimensionalQuantile ( Environment & env ) const {
std : : vector < std : : vector < ValueType > > result ;
auto const & probOpFormula = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) ;
auto const & boundedUntilFormula = probOpFormula . getSubformula ( ) . asBoundedUntilFormula ( ) ;
std : : vector < storm : : storage : : BitVector > dimensionsAsBitVector ;
std : : vector < uint64_t > dimensions ;
std : : vector < storm : : solver : : OptimizationDirection > optimizationDirections ;
std : : vector < bool > lowerCostBounds ;
for ( auto const & dirVar : quantileFormula . getBoundVariables ( ) ) {
dimensionsAsBitVector . push_back ( getDimensionsForVariable ( dirVar . second ) ) ;
STORM_LOG_THROW ( dimensionsAsBitVector . back ( ) . getNumberOfSetBits ( ) = = 1 , storm : : exceptions : : NotSupportedException , " There is not exactly one reward bound referring to quantile variable ' " < < dirVar . second . getName ( ) < < " '. " ) ;
dimensions . push_back ( * dimensionsAsBitVector . back ( ) . begin ( ) ) ;
lowerCostBounds . push_back ( boundedUntilFormula . hasLowerBound ( dimensions . back ( ) ) ) ;
optimizationDirections . push_back ( dirVar . first ) ;
}
STORM_LOG_ASSERT ( dimensions . size ( ) = = 2 , " Expected to have exactly two open dimensions. " ) ;
if ( optimizationDirections [ 0 ] = = optimizationDirections [ 1 ] ) {
if ( lowerCostBounds [ 0 ] = = lowerCostBounds [ 1 ] ) {
// TODO: Assert that we have a reasonable probability bound. STORM_LOG_THROW(storm::solver::minimize(optimizationDirections[0])
bool infInfProbSatisfiesFormula = checkLimitValue ( env , storm : : storage : : BitVector ( getDimension ( ) , false ) ) ;
bool zeroZeroProbSatisfiesFormula = checkLimitValue ( env , dimensionsAsBitVector [ 0 ] | dimensionsAsBitVector [ 1 ] ) ;
if ( infInfProbSatisfiesFormula ! = zeroZeroProbSatisfiesFormula ) {
std : : vector < std : : pair < int64_t , typename ModelType : : ValueType > > singleQuantileValues ;
for ( uint64_t i = 0 ; i < 2 ; + + i ) {
// find out whether the bounds B_i = 0 and B_1-i = infinity satisfy the formula
bool zeroInfSatisfiesFormula = checkLimitValue ( env , dimensionsAsBitVector [ 1 - i ] ) ;
if ( zeroInfSatisfiesFormula = = storm : : solver : : minimize ( optimizationDirections [ 0 ] ) ) {
// There is bound value b such that the point B_i=0 and B_1-i = b is part of the result
singleQuantileValues . emplace_back ( 0 , storm : : utility : : zero < ValueType > ( ) ) ;
} else {
// Compute quantile where 1-i is set to infinity
singleQuantileValues . push_back ( computeQuantileForDimension ( env , dimensions [ i ] ) . get ( ) ) ;
if ( ! storm : : solver : : minimize ( optimizationDirections [ i ] ) ) {
// When maximizing bounds, the computed single dimensional quantile value is sat for all values of the other bound.
// Increase the computed value so that there is at least one unsat assignment of bounds with B[i] = singleQuantileValues[i]
+ + singleQuantileValues . back ( ) . first ;
}
}
// Decrease the value for lower cost bounds to convert >= to >
if ( lowerCostBounds [ i ] ) {
- - singleQuantileValues [ i ] . first ;
}
}
std : : vector < int64_t > currentEpochValues = { ( int64_t ) singleQuantileValues [ 0 ] . first , ( int64_t ) singleQuantileValues [ 1 ] . first } ; // signed int is needed to allow lower bounds >-1 (aka >=0).
while ( ! exploreTwoDimensionalQuantile ( env , singleQuantileValues , currentEpochValues , result ) ) {
increasePrecision ( env ) ;
+ + numPrecisionRefinements ;
}
} else if ( infInfProbSatisfiesFormula ) {
// then also zeroZeroProb satisfies the formula, i.e., all bound values satisfy the formula
if ( storm : : solver : : minimize ( optimizationDirections [ 0 ] ) ) {
result = { { storm : : utility : : zero < ValueType > ( ) , storm : : utility : : zero < ValueType > ( ) } } ;
} else {
result = { { storm : : utility : : infinity < ValueType > ( ) , storm : : utility : : infinity < ValueType > ( ) } } ;
}
} else {
// i.e., no bound value satisfies the formula
result = { { } } ;
}
} else {
// TODO: this is an "unreasonable" case
}
} else {
// TODO: find reasonable min/max cases
}
return result ;
}
template < typename ModelType >
bool QuantileHelper < ModelType > : : exploreTwoDimensionalQuantile ( Environment const & env , std : : vector < std : : pair < int64_t , typename ModelType : : ValueType > > const & startEpochValues , std : : vector < int64_t > & currentEpochValues , std : : vector < std : : vector < ValueType > > & resultPoints ) const {
// Initialize some data for easy access
auto const & probOpFormula = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) ;
auto const & boundedUntilFormula = probOpFormula . getSubformula ( ) . asBoundedUntilFormula ( ) ;
std : : vector < storm : : storage : : BitVector > dimensionsAsBitVector ;
std : : vector < uint64_t > dimensions ;
std : : vector < storm : : solver : : OptimizationDirection > optimizationDirections ;
std : : vector < bool > lowerCostBounds ;
for ( auto const & dirVar : quantileFormula . getBoundVariables ( ) ) {
dimensionsAsBitVector . push_back ( getDimensionsForVariable ( dirVar . second ) ) ;
STORM_LOG_THROW ( dimensionsAsBitVector . back ( ) . getNumberOfSetBits ( ) = = 1 , storm : : exceptions : : NotSupportedException , " There is not exactly one reward bound referring to quantile variable ' " < < dirVar . second . getName ( ) < < " '. " ) ;
dimensions . push_back ( * dimensionsAsBitVector . back ( ) . begin ( ) ) ;
lowerCostBounds . push_back ( boundedUntilFormula . hasLowerBound ( dimensions . back ( ) ) ) ;
optimizationDirections . push_back ( dirVar . first ) ;
}
STORM_LOG_ASSERT ( dimensions . size ( ) = = 2 , " Expected to have exactly two open dimensions. " ) ;
MultiDimensionalRewardUnfolding < ValueType , true > rewardUnfolding ( model , transformBoundedUntilOperator ( quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) , std : : vector < BoundTransformation > ( getDimension ( ) , BoundTransformation : : None ) ) ) ;
// initialize data that will be needed for each epoch
auto lowerBound = rewardUnfolding . getLowerObjectiveBound ( ) ;
auto upperBound = rewardUnfolding . getUpperObjectiveBound ( ) ;
std : : vector < ValueType > x , b ;
std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > minMaxSolver ;
if ( currentEpochValues [ 0 ] < 0 & & currentEpochValues [ 1 ] < 0 ) {
// This case can only happen in these cases:
assert ( lowerCostBounds [ 0 ] ) ;
assert ( currentEpochValues [ 0 ] = = - 1 ) ;
assert ( currentEpochValues [ 1 ] = = - 1 ) ;
// This case has been checked already, so we can skip it.
// Skipping this is actually necessary, since the rewardUnfolding will handle formulas like [F{"a"}>A,{"b"}>B "init"] incorrectly if A=B=-1.
+ + currentEpochValues [ 1 ] ;
}
std : : set < typename EpochManager : : Epoch > exploredEpochs ;
while ( true ) {
auto currEpoch = rewardUnfolding . getStartEpoch ( true ) ;
for ( uint64_t i = 0 ; i < 2 ; + + i ) {
if ( currentEpochValues [ i ] > = 0 ) {
rewardUnfolding . getEpochManager ( ) . setDimensionOfEpoch ( currEpoch , dimensions [ i ] , ( uint64_t ) currentEpochValues [ i ] ) ;
} else {
rewardUnfolding . getEpochManager ( ) . setBottomDimension ( currEpoch , dimensions [ i ] ) ;
}
}
auto epochOrder = rewardUnfolding . getEpochComputationOrder ( currEpoch ) ;
for ( auto const & epoch : epochOrder ) {
if ( exploredEpochs . count ( epoch ) = = 0 ) {
auto & epochModel = rewardUnfolding . setCurrentEpoch ( epoch ) ;
rewardUnfolding . setSolutionForCurrentEpoch ( epochModel . analyzeSingleObjective ( env , quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getOptimalityType ( ) , x , b , minMaxSolver , lowerBound , upperBound ) ) ;
+ + numCheckedEpochs ;
exploredEpochs . insert ( epoch ) ;
}
}
ValueType currValue = rewardUnfolding . getInitialStateResult ( currEpoch ) ;
bool propertySatisfied ;
if ( env . solver ( ) . isForceSoundness ( ) ) {
auto lowerUpperValue = getLowerUpperBound ( env , currValue ) ;
propertySatisfied = probOpFormula . getBound ( ) . isSatisfied ( lowerUpperValue . first ) ;
if ( propertySatisfied ! = probOpFormula . getBound ( ) . isSatisfied ( lowerUpperValue . second ) ) {
// unclear result due to insufficient precision.
return false ;
}
} else {
propertySatisfied = probOpFormula . getBound ( ) . isSatisfied ( currValue ) ;
}
// If the reward bounds should be as small as possible, we should stop as soon as the property is satisfied.
// If the reward bounds should be as large as possible, we should stop as soon as the property is violated and then go a step backwards
if ( storm : : solver : : minimize ( optimizationDirections [ 0 ] ) = = propertySatisfied ) {
// We found another point for the result! Translate it to the original domain
auto point = currentEpochValues ;
std : : vector < ValueType > convertedPoint ;
for ( uint64_t i = 0 ; i < 2 ; + + i ) {
if ( lowerCostBounds [ i ] ) {
// Translate > to >=
+ + point [ i ] ;
}
if ( i = = 1 & & storm : : solver : : maximize ( optimizationDirections [ i ] ) ) {
// When maximizing, we actually searched for each x-value the smallest y-value that leads to a property violation. Hence, decreasing y by one means property satisfaction
- - point [ i ] ;
}
if ( point [ i ] < 0 ) {
// Skip this point
convertedPoint . clear ( ) ;
continue ;
}
convertedPoint . push_back ( storm : : utility : : convertNumber < ValueType > ( point [ i ] ) * rewardUnfolding . getDimension ( dimensions [ i ] ) . scalingFactor ) ;
}
if ( ! convertedPoint . empty ( ) ) {
resultPoints . push_back ( std : : move ( convertedPoint ) ) ;
}
if ( currentEpochValues [ 1 ] = = startEpochValues [ 1 ] . first ) {
break ;
} else {
+ + currentEpochValues [ 0 ] ;
currentEpochValues [ 1 ] = startEpochValues [ 1 ] . first ;
}
} else {
+ + currentEpochValues [ 1 ] ;
}
}
// When maximizing, there are border cases where one dimension can be arbitrarily large
for ( uint64_t i = 0 ; i < 2 ; + + i ) {
if ( storm : : solver : : maximize ( optimizationDirections [ i ] ) & & startEpochValues [ i ] . first > 0 ) {
std : : vector < ValueType > newResultPoint ( 2 ) ;
newResultPoint [ i ] = storm : : utility : : convertNumber < ValueType > ( startEpochValues [ i ] . first - 1 ) * startEpochValues [ i ] . second ;
newResultPoint [ 1 - i ] = storm : : utility : : infinity < ValueType > ( ) ;
resultPoints . push_back ( newResultPoint ) ;
}
}
filterDominatedPoints ( resultPoints , optimizationDirections ) ;
return true ;
}
template < typename ModelType >
boost : : optional < std : : pair < uint64_t , typename ModelType : : ValueType > > QuantileHelper < ModelType > : : computeQuantileForDimension ( Environment const & env , uint64_t dimension ) const {
// We assume that there is one bound value that violates the quantile and one bound value that satisfies it.
// Let all other open bounds approach infinity
std : : set < storm : : expressions : : Variable > infinityVariables ;
for ( auto const & d : getOpenDimensions ( ) ) {
if ( d ! = dimension ) {
infinityVariables . insert ( getVariableForDimension ( d ) ) ;
}
}
auto transformedFormula = transformBoundedUntilOperator ( quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) , std : : vector < BoundTransformation > ( getDimension ( ) , BoundTransformation : : None ) ) ;
MultiDimensionalRewardUnfolding < ValueType , true > rewardUnfolding ( model , transformedFormula , infinityVariables ) ;
bool upperCostBound = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getSubformula ( ) . asBoundedUntilFormula ( ) . hasUpperBound ( dimension ) ;
bool minimizingRewardBound = storm : : solver : : minimize ( getOptimizationDirForDimension ( dimension ) ) ;
// initialize data that will be needed for each epoch
auto lowerBound = rewardUnfolding . getLowerObjectiveBound ( ) ;
auto upperBound = rewardUnfolding . getUpperObjectiveBound ( ) ;
std : : vector < ValueType > x , b ;
std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > minMaxSolver ;
uint64_t epochValue = 0 ;
std : : set < typename EpochManager : : Epoch > exploredEpochs ;
while ( true ) {
auto currEpoch = rewardUnfolding . getStartEpoch ( true ) ;
rewardUnfolding . getEpochManager ( ) . setDimensionOfEpoch ( currEpoch , dimension , ( uint64_t ) epochValue ) ;
auto epochOrder = rewardUnfolding . getEpochComputationOrder ( currEpoch ) ;
for ( auto const & epoch : epochOrder ) {
if ( exploredEpochs . count ( epoch ) = = 0 ) {
auto & epochModel = rewardUnfolding . setCurrentEpoch ( epoch ) ;
rewardUnfolding . setSolutionForCurrentEpoch ( epochModel . analyzeSingleObjective ( env , quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getOptimalityType ( ) , x , b , minMaxSolver , lowerBound , upperBound ) ) ;
+ + numCheckedEpochs ;
exploredEpochs . insert ( epoch ) ;
}
}
ValueType currValue = rewardUnfolding . getInitialStateResult ( currEpoch ) ;
bool propertySatisfied ;
if ( env . solver ( ) . isForceSoundness ( ) ) {
auto lowerUpperValue = getLowerUpperBound ( env , currValue ) ;
propertySatisfied = transformedFormula - > getBound ( ) . isSatisfied ( lowerUpperValue . first ) ;
if ( propertySatisfied ! = transformedFormula - > getBound ( ) . isSatisfied ( lowerUpperValue . second ) ) {
// unclear result due to insufficient precision.
return boost : : none ;
}
} else {
propertySatisfied = transformedFormula - > getBound ( ) . isSatisfied ( currValue ) ;
}
// If the reward bound should be as small as possible, we should stop as soon as the property is satisfied.
// If the reward bound should be as large as possible, we should stop as soon as sthe property is violated and then go a step backwards
if ( minimizingRewardBound & & propertySatisfied ) {
// We found a satisfying value!
if ( ! upperCostBound ) {
// The rewardunfolding assumes strict lower cost bounds while we assume non-strict ones. Hence, >B becomes >=(B+1).
+ + epochValue ;
}
break ;
} else if ( ! minimizingRewardBound & & ! propertySatisfied ) {
// We found a non-satisfying value. Go one step back to get the largest satisfying value.
// ... however, lower cost bounds need to be converted from strict bounds >B to non strict bounds >=(B+1).
// Hence, no operation is necessary in this case.
if ( upperCostBound ) {
STORM_LOG_ASSERT ( epochValue > 0 , " The property does not seem to be satisfiable. This case should have been excluded earlier. " ) ;
- - epochValue ;
}
break ;
}
+ + epochValue ;
}
return std : : make_pair ( epochValue , rewardUnfolding . getDimension ( dimension ) . scalingFactor ) ;
}
template < typename ModelType >
bool QuantileHelper < ModelType > : : checkLimitValue ( Environment & env , storm : : storage : : BitVector const & infDimensions ) const {
auto const & probabilityBound = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getBound ( ) ;
// Increase the precision until we get a conclusive result
while ( true ) {
ValueType numericResult = computeLimitValue ( env , infDimensions ) ;
if ( env . solver ( ) . isForceSoundness ( ) ) {
auto lowerUpper = getLowerUpperBound ( env , numericResult ) ;
bool lowerSat = probabilityBound . isSatisfied ( lowerUpper . first ) ;
bool upperSat = probabilityBound . isSatisfied ( lowerUpper . second ) ;
if ( lowerSat = = upperSat ) {
return lowerSat ;
} else {
increasePrecision ( env ) ;
+ + numPrecisionRefinements ;
}
} else {
return probabilityBound . isSatisfied ( numericResult ) ;
}
}
}
template < typename ModelType >
typename ModelType : : ValueType QuantileHelper < ModelType > : : computeLimitValue ( Environment const & env , storm : : storage : : BitVector const & infDimensions ) const {
// To compute the limit for an upper bounded dimension, we can simply drop the bound
// To compute the limit for a lower bounded dimension, we only consider epoch steps that lie in an end component and check for the bound >0 instead.
// Notice, however, that this approach becomes problematic if, at the same time, we consider an upper bounded dimension with bound value zero.
std : : vector < BoundTransformation > bts ( getDimension ( ) , BoundTransformation : : None ) ;
std : : set < storm : : expressions : : Variable > infinityVariables ;
for ( auto const & d : getOpenDimensions ( ) ) {
bool upperCostBound = quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getSubformula ( ) . asBoundedUntilFormula ( ) . hasUpperBound ( d ) ;
if ( infDimensions . get ( d ) ) {
infinityVariables . insert ( getVariableForDimension ( d ) ) ;
} else {
if ( upperCostBound ) {
bts [ d ] = BoundTransformation : : LessEqualZero ;
} else {
bts [ d ] = BoundTransformation : : GreaterEqualZero ;
}
}
}
auto transformedFormula = transformBoundedUntilOperator ( quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) , bts ) ;
MultiDimensionalRewardUnfolding < ValueType , true > rewardUnfolding ( model , transformedFormula , infinityVariables ) ;
if ( ! rewardUnfolding . getProb1Objectives ( ) . empty ( ) ) {
assert ( rewardUnfolding . getProb1Objectives ( ) . size ( ) = = 1 ) ;
// The probability is one.
return storm : : utility : : one < ValueType > ( ) ;
}
// Get lower and upper bounds for the solution.
auto lowerBound = rewardUnfolding . getLowerObjectiveBound ( ) ;
auto upperBound = rewardUnfolding . getUpperObjectiveBound ( ) ;
// Initialize epoch models
auto initEpoch = rewardUnfolding . getStartEpoch ( ) ;
auto epochOrder = rewardUnfolding . getEpochComputationOrder ( initEpoch ) ;
// initialize data that will be needed for each epoch
std : : vector < ValueType > x , b ;
std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > minMaxSolver ;
for ( auto const & epoch : epochOrder ) {
auto & epochModel = rewardUnfolding . setCurrentEpoch ( epoch ) ;
rewardUnfolding . setSolutionForCurrentEpoch ( epochModel . analyzeSingleObjective ( env , quantileFormula . getSubformula ( ) . asProbabilityOperatorFormula ( ) . getOptimalityType ( ) , x , b , minMaxSolver , lowerBound , upperBound ) ) ;
+ + numCheckedEpochs ;
}
ValueType numericResult = rewardUnfolding . getInitialStateResult ( initEpoch ) ;
STORM_LOG_TRACE ( " Limit probability for infinity dimensions " < < infDimensions < < " is " < < numericResult < < " . " ) ;
return numericResult ;
}
template class QuantileHelper < storm : : models : : sparse : : Mdp < double > > ;
template class QuantileHelper < storm : : models : : sparse : : Mdp < storm : : RationalNumber > > ;
Tim Quatmann
6 years ago