@ -134,14 +134,23 @@ namespace storm {
void SparseMarkovAutomatonCslHelper : : computeBoundedReachabilityProbabilities ( OptimizationDirection dir , storm : : storage : : SparseMatrix < ValueType > const & transitionMatrix , std : : vector < ValueType > const & exitRates , storm : : storage : : BitVector const & goalStates , storm : : storage : : BitVector const & markovianNonGoalStates , storm : : storage : : BitVector const & probabilisticNonGoalStates , std : : vector < ValueType > & markovianNonGoalValues , std : : vector < ValueType > & probabilisticNonGoalValues , ValueType delta , uint64_t numberOfSteps , storm : : solver : : MinMaxLinearEquationSolverFactory < ValueType > const & minMaxLinearEquationSolverFactory ) {
STORM_LOG_THROW ( false , storm : : exceptions : : InvalidOperationException , " Computing bounded reachability probabilities is unsupported for this value type. " ) ;
}
template < typename ValueType , typename std : : enable_if < storm : : NumberTraits < ValueType > : : SupportsExponential , int > : : type >
void SparseMarkovAutomatonCslHelper : : printTransitions ( std : : vector < std : : vector < ValueType > > relReachability , std : : vector < ValueType > const & exitRateVector , storm : : storage : : SparseMatrix < ValueType > const & fullTransitionMatrix , storm : : storage : : BitVector const & markovianStates , storm : : storage : : BitVector const & psiStates , storm : : storage : : BitVector const & cycleStates , storm : : storage : : BitVector const & cycleGoalStates , std : : vector < std : : vector < ValueType > > & vd , std : : vector < std : : vector < ValueType > > & vu , std : : vector < std : : vector < ValueType > > & wu ) {
template < typename ValueType , typename std : : enable_if < storm : : NumberTraits < ValueType > : : SupportsExponential , int > : : type = 0 >
void SparseMarkovAutomatonCslHelper : : printTransitions (
storm : : storage : : SparseMatrix < ValueType > const & fullTransitionMatrix ,
std : : vector < ValueType > const & exitRateVector , storm : : storage : : BitVector const & markovianStates ,
storm : : storage : : BitVector const & psiStates , std : : vector < std : : vector < ValueType > > relReachability ,
const storage : : BitVector & cycleStates , const storage : : BitVector & cycleGoalStates ,
std : : vector < std : : vector < std : : vector < ValueType > > > & unifVectors ) {
std : : ofstream logfile ( " U+logfile.txt " , std : : ios : : app ) ;
auto const & rowGroupIndices = fullTransitionMatrix . getRowGroupIndices ( ) ;
auto numberOfStates = fullTransitionMatrix . getRowGroupCount ( ) ;
//Transition Matrix
logfile < < " number of states = num of row group count " < < numberOfStates < < " \n " ;
for ( uint_fast64_t i = 0 ; i < fullTransitionMatrix . getRowGroupCount ( ) ; i + + ) {
logfile < < " from node " < < i < < " " ;
@ -161,8 +170,8 @@ namespace storm {
logfile < < " \n " ;
logfile < < " probStates \t markovianStates \t goalStates \t cycleStates \t cycleGoalStates \n " ;
for ( int i = 0 ; i < markovianStates . size ( ) ; i + + ) {
logfile < < ( ~ markovianStates ) [ i ] < < " \t \t " < < markovianStates [ i ] < < " \t \t " < < psiStates [ i ] < < " \t \t " < < cycleStates [ i ] < < " \t \t " < < cycleGoalStates [ i ] < < " \n " ;
for ( int i = 0 ; i < markovianStates . size ( ) ; i + + ) {
logfile < < ( ~ markovianStates ) [ i ] < < " \t \t " < < markovianStates [ i ] < < " \t \t " < < psiStates [ i ] < < " \t \t " < < cycleStates [ i ] < < " \t \t " < < cycleGoalStates [ i ] < < " \n " ;
}
for ( int i = 0 ; i < relReachability . size ( ) ; i + + ) {
@ -174,25 +183,25 @@ namespace storm {
logfile < < " vd: \n " ;
for ( uint_fast 64_t i = 0 ; i < vd . size ( ) ; i + + ) {
for ( uint_fast 64_t j = 0 ; j < fullTransitionMatrix . getRowGroupCount ( ) ; j + + ) {
logfile < < vd [ i ] [ j ] < < " \t " ;
for ( uint64_t i = 0 ; i < unifVectors [ 0 ] . size ( ) ; i + + ) {
for ( uint64_t j = 0 ; j < unifVectors [ 0 ] [ i ] . size ( ) ; j + + ) {
logfile < < unifVectors [ 0 ] [ i ] [ j ] < < " \t " ;
}
logfile < < " \n " ;
}
logfile < < " \n vu: \n " ;
for ( uint_fast 64_t i = 0 ; i < v u. size ( ) ; i + + ) {
for ( uint_fast 64_t j = 0 ; j < fullTransitionMatrix . getRowGroupCount ( ) ; j + + ) {
logfile < < v u[ i ] [ j ] < < " \t " ;
for ( uint64_t i = 0 ; i < unifVectors [ 1 ] . size ( ) ; i + + ) {
for ( uint64_t j = 0 ; j < unifVectors [ 1 ] [ i ] . size ( ) ; j + + ) {
logfile < < unifVectors [ 1 ] [ i ] [ j ] < < " \t " ;
}
logfile < < " \n " ;
}
logfile < < " \n wu \n " ;
for ( uint_fast 64_t i = 0 ; i < w u. size ( ) ; i + + ) {
for ( uint_fast 64_t j = 0 ; j < fullTransitionMatrix . getRowGroupCount ( ) ; j + + ) {
logfile < < w u[ i ] [ j ] < < " \t " ;
for ( uint64_t i = 0 ; i < unifVectors [ 2 ] . size ( ) ; i + + ) {
for ( uint64_t j = 0 ; j < unifVectors [ 2 ] [ i ] . size ( ) ; j + + ) {
logfile < < unifVectors [ 2 ] [ i ] [ j ] < < " \t " ;
}
logfile < < " \n " ;
}
@ -212,256 +221,137 @@ namespace storm {
return res ;
}
template < typename ValueType , typename std : : enable_if < storm : : NumberTraits < ValueType > : : SupportsExponential , int > : : type >
void SparseMarkovAutomatonCslHelper : : calculateVu ( std : : vector < std : : vector < ValueType > > const & relativeReachability , OptimizationDirection dir , uint64_t k , uint64_t node , ValueType lambda , std : : vector < std : : vector < ValueType > > & vu , std : : vector < std : : vector < ValueType > > & wu , storm : : storage : : SparseMatrix < ValueType > const & fullTransitionMatrix , storm : : storage : : BitVector const & markovianStates , storm : : storage : : BitVector const & psiStates , std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > const & solver ) {
if ( v u[ k ] [ node ] ! = - 1 ) { return ; } //dynamic programming. avoiding multiple calculation.
uint64_t N = v u. size ( ) - 1 ;
void SparseMarkovAutomatonCslHelper : : calculateVu ( std : : vector < std : : vector < ValueType > > const & relativeReachability , OptimizationDirection dir , uint64_t k , uint64_t node , uint64_t const kind , ValueType lambda , uint64_t probSize , std : : vector < std : : vector < std : : vector < ValueType > > > & unifVectors , storm : : storage : : SparseMatrix < ValueType > const & fullTransitionMatrix , storm : : storage : : BitVector const & markovianStates , storm : : storage : : BitVector const & psiStates , std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > const & solver ) {
if ( unifVectors [ 1 ] [ k ] [ node ] ! = - 1 ) { return ; } //dynamic programming. avoiding multiple calculation.
uint64_t N = unifVectors [ 1 ] . size ( ) - 1 ;
auto rowGroupIndices = fullTransitionMatrix . getRowGroupIndices ( ) ;
ValueType res = 0 ;
for ( uint64_t i = k ; i < N ; i + + ) {
if ( wu [ N - 1 - ( i - k ) ] [ node ] = = - 1 ) {
calculateWu ( relativeReachability , dir , ( N - 1 - ( i - k ) ) , node , lambda , wu , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
if ( unifVectors [ 2 ] [ N - 1 - ( i - k ) ] [ node ] = = - 1 ) {
calculateUnifPlusVector ( N - 1 - ( i - k ) , node , 2 , lambda , probSize , relativeReachability , dir , unifVectors , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
//old: relativeReachability, dir, (N-1-(i-k)),node,lambda,wu,fullTransitionMatrix,markovianStates,psiStates, solver);
}
res + = poisson ( lambda , i ) * w u[ N - 1 - ( i - k ) ] [ node ] ;
res + = poisson ( lambda , i ) * unifVectors [ 2 ] [ N - 1 - ( i - k ) ] [ node ] ;
}
v u[ k ] [ node ] = res ;
unifVectors [ 1 ] [ k ] [ node ] = res ;
}
template < typename ValueType , typename std : : enable_if < storm : : NumberTraits < ValueType > : : SupportsExponential , int > : : type >
void SparseMarkovAutomatonCslHelper : : calculateWu ( std : : vector < std : : vector < ValueType > > const & relativeReachability , OptimizationDirection dir , uint64_t k , uint64_t node , ValueType lambda , std : : vector < std : : vector < ValueType > > & wu , storm : : storage : : SparseMatrix < ValueType > const & fullTransitionMatrix , storm : : storage : : BitVector const & markovianStates , storm : : storage : : BitVector const & psiStates , std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > const & solver ) {
if ( wu [ k ] [ node ] ! = - 1 ) { return ; } //dynamic programming. avoiding multiple calculation.
std : : ofstream logfile ( " U+logfile.txt " , std : : ios : : app ) ;
auto probabilisticStates = ~ markovianStates ;
uint64_t N = wu . size ( ) - 1 ;
auto const & rowGroupIndices = fullTransitionMatrix . getRowGroupIndices ( ) ;
ValueType res ;
if ( k = = N ) {
wu [ k ] [ node ] = 0 ;
return ;
}
if ( psiStates [ node ] ) {
wu [ k ] [ node ] = 1 ;
return ;
}
template < typename ValueType , typename std : : enable_if < storm : : NumberTraits < ValueType > : : SupportsExponential , int > : : type = 0 >
void SparseMarkovAutomatonCslHelper : : calculateUnifPlusVector ( uint64_t k , uint64_t node , uint64_t const kind , ValueType lambda , uint64_t probSize ,
std : : vector < std : : vector < ValueType > > const & relativeReachability ,
OptimizationDirection dir ,
std : : vector < std : : vector < std : : vector < ValueType > > > & unifVectors ,
storm : : storage : : SparseMatrix < ValueType > const & fullTransitionMatrix ,
storm : : storage : : BitVector const & markovianStates ,
storm : : storage : : BitVector const & psiStates ,
std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > const & solver ) {
if ( markovianStates [ node ] ) {
res = 0 ;
auto line = fullTransitionMatrix . getRow ( rowGroupIndices [ node ] ) ;
for ( auto & element : line ) {
uint64_t to = element . getColumn ( ) ;
if ( wu [ k + 1 ] [ to ] = = - 1 ) {
calculateWu ( relativeReachability , dir , k + 1 , to , lambda , wu , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
}
res + = element . getValue ( ) * wu [ k + 1 ] [ to ] ;
}
wu [ k ] [ node ] = res ;
} else {
// applying simpleMDP stuff
auto numberOfStates = relativeReachability [ 0 ] . size ( ) ;
typename storm : : storage : : SparseMatrix < ValueType > probMatrix = fullTransitionMatrix . getSubmatrix ( true , ~ markovianStates , ~ markovianStates , true ) ;
auto probSize = probMatrix . getRowGroupCount ( ) ;
std : : vector < ValueType > b ( probMatrix . getRowCount ( ) , 0 ) , x ( probSize , 0 ) ;
//calculate b
uint64_t lineCounter = 0 ;
for ( int i = 0 ; i < numberOfStates ; i + + ) {
if ( markovianStates [ i ] ) {
continue ;
}
auto rowStart = rowGroupIndices [ i ] ;
auto rowEnd = rowGroupIndices [ i + 1 ] ;
for ( auto j = rowStart ; j < rowEnd ; j + + ) {
ValueType res = 0 ;
for ( auto & element : fullTransitionMatrix . getRow ( j ) ) {
auto to = element . getColumn ( ) ;
if ( probabilisticStates [ to ] ) {
continue ;
}
if ( wu [ k ] [ to ] = = - 1 ) {
calculateWu ( relativeReachability , dir , k , to , lambda , wu , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
}
res = res + relativeReachability [ j ] [ to ] * wu [ k ] [ to ] ;
}
b [ lineCounter ] = res ;
lineCounter + + ;
}
}
logfile < < " \n wu reach vector b is \n " ;
for ( uint64_t i = 0 ; i < probSize ; i + + ) {
logfile < < b [ i ] < < " \n " ;
}
solver - > solveEquations ( dir , x , b ) ;
for ( uint64_t i = 0 ; i < probSize ; i + + ) {
auto trueI = transformIndice ( ~ markovianStates , i ) ;
wu [ k ] [ trueI ] = x [ i ] ;
}
/* unifPlus
res = - 1 ;
uint64_t rowStart = rowGroupIndices [ node ] ;
uint64_t rowEnd = rowGroupIndices [ node + 1 ] ;
for ( uint64_t i = rowStart ; i < rowEnd ; i + + ) {
auto line = fullTransitionMatrix . getRow ( i ) ;
ValueType between = 0 ;
for ( auto & element : line ) {
uint64_t to = element . getColumn ( ) ;
if ( to = = node ) {
continue ;
}
if ( wu [ k ] [ to ] = = - 1 ) {
calculateWu ( dir , k , to , lambda , wu , fullTransitionMatrix , markovianStates , psiStates ) ;
}
between + = element . getValue ( ) * wu [ k ] [ to ] ;
}
if ( maximize ( dir ) ) {
res = std : : max ( res , between ) ;
} else {
if ( res ! = - 1 ) {
res = std : : min ( res , between ) ;
} else {
res = between ;
}
}
} */
} // end no goal-prob state
}
template < typename ValueType , typename std : : enable_if < storm : : NumberTraits < ValueType > : : SupportsExponential , int > : : type >
void SparseMarkovAutomatonCslHelper : : calculateVd ( std : : vector < std : : vector < ValueType > > const & relativeReachability , OptimizationDirection dir , uint64_t k , uint64_t node , ValueType lambda , std : : vector < std : : vector < ValueType > > & vd , storm : : storage : : SparseMatrix < ValueType > const & fullTransitionMatrix , storm : : storage : : BitVector const & markovianStates , storm : : storage : : BitVector const & psiStates , std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > const & solver ) {
if ( unifVectors [ kind ] [ k ] [ node ] ! = - 1 ) { return ; }
std : : ofstream logfile ( " U+logfile.txt " , std : : ios : : app ) ;
logfile < < " calculating vector " < < kind < < " for k = " < < k < < " node " < < node < < " \t " ;
auto probabilisticStates = ~ markovianStates ;
if ( vd [ k ] [ node ] ! = - 1 ) { return ; } //dynamic programming. avoiding multiple calculation.
logfile < < " calculating vd for k = " < < k < < " node " < < node < < " \t " ;
uint64_t N = vd . size ( ) - 1 ;
auto numberOfStates = fullTransitionMatrix . getRowGroupCount ( ) ;
uint64_t N = unifVectors [ kind ] . size ( ) - 1 ;
auto const & rowGroupIndices = fullTransitionMatrix . getRowGroupIndices ( ) ;
ValueType res ;
// First Case, k==N, independent from kind of state
if ( k = = N ) {
logfile < < " k == N! res = 0 \n " ;
vd [ k ] [ node ] = 0 ;
unifVectors [ kind ] [ k ] [ node ] = 0 ;
return ;
}
//goal state
if ( psiStates [ node ] ) {
res = storm : : utility : : zero < ValueType > ( ) ;
for ( uint64_t i = k ; i < N ; i + + ) {
ValueType between = poisson ( lambda , i ) ;
res + = between ;
if ( kind = = 0 ) {
// Vd
res = storm : : utility : : zero < ValueType > ( ) ;
for ( uint64_t i = k ; i < N ; i + + ) {
ValueType between = poisson ( lambda , i ) ;
res + = between ;
}
unifVectors [ kind ] [ k ] [ node ] = res ;
} else {
// WU
unifVectors [ kind ] [ k ] [ node ] = 1 ;
}
vd [ k ] [ node ] = res ;
logfile < < " goal state node " < < node < < " res = " < < res < < " \n " ;
return ;
}
// no-goal markovian state
//markovian non-goal State
if ( markovianStates [ node ] ) {
logfile < < " markovian state: " ;
res = 0 ;
auto line = fullTransitionMatrix . getRow ( rowGroupIndices [ node ] ) ;
for ( auto & element : line ) {
uint64_t to = element . getColumn ( ) ;
if ( vd [ k + 1 ] [ to ] = = - 1 ) {
calculateVd ( relativeReachability , dir , k + 1 , to , lambda , vd , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
if ( unifVectors [ kind ] [ k + 1 ] [ to ] = = - 1 ) {
calculateUnifPlusVector ( k + 1 , to , kind , lambda , probSize , relativeReachability , dir , unifVectors , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
}
res + = element . getValue ( ) * vd [ k + 1 ] [ to ] ;
res + = element . getValue ( ) * unifVectors [ kind ] [ k + 1 ] [ to ] ;
}
vd [ k ] [ node ] = res ;
unifVectors [ kind ] [ k ] [ node ] = res ;
logfile < < " res = " < < res < < " \n " ;
} else { //no-goal prob state
// applying simpleMDP stuff
logfile < < " probState, applying MDP stuff " ;
auto numberOfStates = relativeReachability [ 0 ] . size ( ) ;
typename storm : : storage : : SparseMatrix < ValueType > probMatrix = fullTransitionMatrix . getSubmatrix ( true , ~ markovianStates , ~ markovianStates , true ) ;
auto probSize = probMatrix . getRowGroupCount ( ) ;
std : : vector < ValueType > b ( probMatrix . getRowCount ( ) , 0 ) , x ( probSize , 0 ) ;
return ;
}
//probabilistic non-goal State
if ( probabilisticStates [ node ] ) {
logfile < < " probabilistic state: " ;
std : : vector < ValueType > b ( probSize , 0 ) , x ( probabilisticStates . getNumberOfSetBits ( ) , 0 ) ;
//calculate b
uint64_t lineCounter = 0 ;
for ( int i = 0 ; i < numberOfStates ; i + + ) {
if ( markovianStates [ i ] ) {
for ( int i = 0 ; i < numberOfStates ; i + + ) {
if ( markovianStates [ i ] ) {
continue ;
}
auto rowStart = rowGroupIndices [ i ] ;
auto rowEnd = rowGroupIndices [ i + 1 ] ;
for ( auto j = rowStart ; j < rowEnd ; j + + ) {
ValueType res = 0 ;
for ( auto & element : fullTransitionMatrix . getRow ( j ) ) {
auto rowEnd = rowGroupIndices [ i + 1 ] ;
for ( auto j = rowStart ; j < rowEnd ; j + + ) {
res = 0 ;
for ( auto & element : fullTransitionMatrix . getRow ( j ) ) {
auto to = element . getColumn ( ) ;
if ( probabilisticStates [ to ] ) {
if ( probabilisticStates [ to ] ) {
continue ;
}
if ( vd [ k ] [ to ] = = - 1 ) {
calculateVd ( relativeReachability , dir , k , to , lambda , vd , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
if ( unifVectors [ kind ] [ k ] [ to ] = = - 1 ) {
calculateUnifPlusVector ( k , to , kind , lambda , probSize , relativeReachability , dir ,
unifVectors , fullTransitionMatrix , markovianStates ,
psiStates , solver ) ;
}
res = res + relativeReachability [ j ] [ to ] * vd [ k ] [ to ] ;
res = res + relativeReachability [ j ] [ to ] * unifVectors [ kind ] [ k ] [ to ] ;
}
b [ lineCounter ] = res ;
b [ lineCounter ] = res ;
lineCounter + + ;
}
}
logfile < < " \n vd reach vector b is \n " ;
for ( uint64_t i = 0 ; i < probMatrix . getRowCount ( ) ; i + + ) {
logfile < < b [ i ] < < " \n " ;
}
solver - > solveEquations ( dir , x , b ) ;
logfile < < " vd goal vector x is \n " ;
for ( int i = 0 ; i < probMatrix . getRowCount ( ) ; i + + ) {
logfile < < x [ i ] < < " \n " ;
}
for ( uint64_t i = 0 ; i < probSize ; i + + ) {
auto trueI = transformIndice ( ~ markovianStates , i ) ;
vd [ k ] [ trueI ] = x [ i ] ;
}
/* unif Plus
logfile < < " prob state: " ;
res = - 1 ;
uint64_t rowStart = rowGroupIndices [ node ] ;
uint64_t rowEnd = rowGroupIndices [ node + 1 ] ;
for ( uint64_t i = rowStart ; i < rowEnd ; i + + ) {
auto line = fullTransitionMatrix . getRow ( i ) ;
ValueType between = 0 ;
for ( auto & element : line ) {
uint64_t to = element . getColumn ( ) ;
if ( to = = node ) {
logfile < < " ignoring self loops for now " ;
continue ;
}
if ( vd [ k ] [ to ] = = - 1 ) {
calculateVd ( dir , k , to , lambda , vd , fullTransitionMatrix , markovianStates , psiStates ) ;
}
between + = element . getValue ( ) * vd [ k ] [ to ] ;
}
if ( maximize ( dir ) ) {
res = std : : max ( res , between ) ;
} else {
if ( res ! = - 1 ) {
res = std : : min ( res , between ) ;
} else {
res = between ;
}
}
auto trueI = transformIndice ( probabilisticStates , i ) ;
unifVectors [ kind ] [ k ] [ trueI ] = x [ i ] ;
}
*/
}
logfile < < " res = " < < unifVectors [ kind ] [ k ] [ node ] < < " but calculated more \n " ;
} //end probabilistic states
}
template < typename ValueType , typename std : : enable_if < storm : : NumberTraits < ValueType > : : SupportsExponential , int > : : type = 0 >
uint64_t SparseMarkovAutomatonCslHelper : : trajans ( storm : : storage : : SparseMatrix < ValueType > const & transitionMatrix , uint64_t node , std : : vector < uint64_t > & disc , std : : vector < uint64_t > & finish , uint64_t * counter ) {
auto const & rowGroupIndice = transitionMatrix . getRowGroupIndices ( ) ;
@ -585,24 +475,29 @@ namespace storm {
std : : vector < ValueType > SparseMarkovAutomatonCslHelper : : unifPlus ( OptimizationDirection dir , std : : pair < double , double > const & boundsPair , std : : vector < ValueType > const & exitRateVector , storm : : storage : : SparseMatrix < ValueType > const & transitionMatrix , storm : : storage : : BitVector const & markovianStates , storm : : storage : : BitVector const & psiStates , storm : : solver : : MinMaxLinearEquationSolverFactory < ValueType > const & minMaxLinearEquationSolverFactory ) {
STORM_LOG_TRACE ( " Using UnifPlus to compute bounded until probabilities. " ) ;
std : : ofstream logfile ( " U+logfile.txt " , std : : ios : : app ) ;
ValueType maxNorm = 0 ;
ValueType maxNorm = storm : : utility : : zero < ValueType > ( ) ;
//bitvectors to identify different kind of states
storm : : storage : : BitVector allStates ( markovianStates . size ( ) , true ) ;
auto probabilisticStates = ~ markovianStates ;
storm : : storage : : BitVector probabilisticStates = ~ markovianStates ;
//vectors to save calculation
std : : vector < std : : vector < ValueType > > vd , vu , wu ;
std : : vector < std : : vector < std : : vector < ValueType > > > unifVectors { } ;
//transition matrix with diagonal entries. The values can be changed during uniformisation
std : : vector < ValueType > exitRate { exitRateVector } ;
typename storm : : storage : : SparseMatrix < ValueType > fullTransitionMatrix = transitionMatrix . getSubmatrix ( true , allStates , allStates , true ) ;
typename storm : : storage : : SparseMatrix < ValueType > probMatrix = fullTransitionMatrix . getSubmatrix ( true , ~ markovianStates , ~ markovianStates , true ) ;
typename storm : : storage : : SparseMatrix < ValueType > probMatrix { } ;
uint64_t probSize = 0 ;
if ( probabilisticStates . getNumberOfSetBits ( ) ! = 0 ) { //work around in case there are no prob states
probMatrix = fullTransitionMatrix . getSubmatrix ( true , probabilisticStates , probabilisticStates , true ) ;
probSize = probMatrix . getRowCount ( ) ;
}
auto rowGroupIndices = fullTransitionMatrix . getRowGroupIndices ( ) ;
auto & rowGroupIndices = fullTransitionMatrix . getRowGroupIndices ( ) ;
//(1) define horizon, epsilon, kappa , N, lambda,
@ -616,13 +511,14 @@ namespace storm {
}
uint64_t N ;
//calculate relative ReachabilityVectors
std : : vector < ValueType > in ( numberOfStates , 0 ) ;
std : : vector < std : : vector < ValueType > > relReachability ( fullTransitionMatrix . getRowCount ( ) , in ) ;
printTransitions ( relReachability , exitRateVector , fullTransitionMatrix , markovianStates , psiStates , psiStates , psiStates , vd , vu , wu ) ; // TODO: delete when develepmont is finished
//calculate relative reachability
//calculate relative reachability
for ( uint64_t i = 0 ; i < numberOfStates ; i + + ) {
if ( markovianStates [ i ] ) {
@ -645,15 +541,14 @@ namespace storm {
requirements . clearBounds ( ) ;
STORM_LOG_THROW ( requirements . empty ( ) , storm : : exceptions : : UncheckedRequirementException , " Cannot establish requirements for solver. " ) ;
std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > solver = minMaxLinearEquationSolverFactory . create ( probMatrix ) ;
solver - > setHasUniqueSolution ( ) ;
solver - > setBounds ( storm : : utility : : zero < ValueType > ( ) , storm : : utility : : one < ValueType > ( ) ) ;
solver - > setRequirementsChecked ( ) ;
solver - > setCachingEnabled ( true ) ;
printTransitions ( relReachability , exitRateVector , fullTransitionMatrix , markovianStates , psiStates , psiStates , psiStates , vd , vu , wu ) ; // TODO: delete when develepmont is finished
std : : unique_ptr < storm : : solver : : MinMaxLinearEquationSolver < ValueType > > solver ;
if ( probSize ! = 0 ) {
solver = minMaxLinearEquationSolverFactory . create ( probMatrix ) ;
solver - > setHasUniqueSolution ( ) ;
solver - > setBounds ( storm : : utility : : zero < ValueType > ( ) , storm : : utility : : one < ValueType > ( ) ) ;
solver - > setRequirementsChecked ( ) ;
solver - > setCachingEnabled ( true ) ;
}
// while not close enough to precision:
do {
// (2) update parameter
@ -692,29 +587,37 @@ namespace storm {
vd = std : : vector < std : : vector < ValueType > > ( N + 1 , init ) ;
vu = std : : vector < std : : vector < ValueType > > ( N + 1 , init ) ;
wu = std : : vector < std : : vector < ValueType > > ( N + 1 , init ) ;
unifVectors . clear ( ) ;
unifVectors . push_back ( vd ) ;
unifVectors . push_back ( vd ) ;
unifVectors . push_back ( vd ) ;
printTransitions ( fullTransitionMatrix , exitRate , markovianStates , psiStates , relReachability , psiStates , psiStates , unifVectors ) ; //TODO remove
//define 0=vd 1=vu 2=wu
// (5) calculate vectors and maxNorm
for ( uint64_t i = 0 ; i < numberOfStates ; i + + ) {
for ( uint64_t k = N ; k < = N ; k - - ) {
calculateVd ( relReachability , dir , k , i , T * lambda , vd , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
calculateWu ( relReachability , dir , k , i , T * lambda , wu , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
calculateVu ( relReachability , dir , k , i , T * lambda , vu , wu , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
//also use iteration to keep maxNorm of vd and vu u p to date, so the loop-condition is easy to prove
ValueType diff = std : : abs ( vd [ k ] [ i ] - v u[ k ] [ i ] ) ;
calculateUnifPlusVector ( k , i , 0 , lambda , probSize , relReachability , dir , unifVectors , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
calculateUnifPlusVector ( k , i , 2 , lambda , probSize , relReachability , dir , unifVectors , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
calculateVu ( relReachability , dir , k , i , 1 , lambda , probSize , unifVectors , fullTransitionMatrix , markovianStates , psiStates , solver ) ;
//also use iteration to keep maxNorm of vd and vup to date, so the loop-condition is easy to prove
ValueType diff = std : : abs ( unifVectors [ 0 ] [ k ] [ i ] - unifVectors [ 1 ] [ k ] [ i ] ) ;
maxNorm = std : : max ( maxNorm , diff ) ;
}
}
printTransitions ( relReachability , exitRateVector , fullTransitionMatrix , markovianStates , psiStates , psiStates , psiStates , vd , vu , wu ) ; // TODO: delete when develepmont is finished
std : : cout < < " \n TBU was " < < unifVectors [ 0 ] [ 0 ] [ 0 ] < < " \n " ;
printTransitions ( fullTransitionMatrix , exitRate , markovianStates , psiStates , relReachability , psiStates , psiStates , unifVectors ) ; //TODO remove
// (6) double lambda
// (6) double lambda
lambda = 2 * lambda ;
} while ( maxNorm > epsilon * ( 1 - kappa ) ) ;
return vd [ 0 ] ;
return unifVectors [ 0 ] [ 0 ] ;
}
template < typename ValueType , typename std : : enable_if < storm : : NumberTraits < ValueType > : : SupportsExponential , int > : : type >
@ -1234,7 +1137,7 @@ namespace storm {
return v . front ( ) * uniformizationRate ;
}
template std : : vector < double > SparseMarkovAutomatonCslHelper : : computeBoundedUntilProbabilities ( OptimizationDirection dir , storm : : storage : : SparseMatrix < double > const & transitionMatrix , std : : vector < double > const & exitRateVector , storm : : storage : : BitVector const & markovianStates , storm : : storage : : BitVector const & psiStates , std : : pair < double , double > const & boundsPair , storm : : solver : : MinMaxLinearEquationSolverFactory < double > const & minMaxLinearEquationSolverFactory ) ;
template std : : vector < double > SparseMarkovAutomatonCslHelper : : computeUntilProbabilities ( OptimizationDirection dir , storm : : storage : : SparseMatrix < double > const & transitionMatrix , storm : : storage : : SparseMatrix < double > const & backwardTransitions , storm : : storage : : BitVector const & phiStates , storm : : storage : : BitVector const & psiStates , bool qualitative , storm : : solver : : MinMaxLinearEquationSolverFactory < double > const & minMaxLinearEquationSolverFactory ) ;