diff --git a/src/adapters/ExplicitModelAdapter.h b/src/adapters/ExplicitModelAdapter.h index ea5a1491c..67c1bb805 100644 --- a/src/adapters/ExplicitModelAdapter.h +++ b/src/adapters/ExplicitModelAdapter.h @@ -739,6 +739,7 @@ namespace storm { static std::vector buildStateRewards(std::vector const& rewards, StateInformation const& stateInformation, expressions::ExpressionEvaluation& eval) { std::vector result(stateInformation.reachableStates.size()); for (uint_fast64_t index = 0; index < stateInformation.reachableStates.size(); index++) { + std::cout << index << ": " << *stateInformation.reachableStates[index] << std::endl; result[index] = ValueType(0); for (auto const& reward : rewards) { // Add this reward to the state if the state is included in the state reward. diff --git a/src/modelchecker/reachability/SparseSccModelChecker.cpp b/src/modelchecker/reachability/SparseSccModelChecker.cpp index 99622b0b9..d239b8b3e 100644 --- a/src/modelchecker/reachability/SparseSccModelChecker.cpp +++ b/src/modelchecker/reachability/SparseSccModelChecker.cpp @@ -66,8 +66,28 @@ namespace storm { } // Finally eliminate initial state. - STORM_PRINT_AND_LOG("Eliminating initial state " << *initialStates.begin() << "." << std::endl); - eliminateState(flexibleMatrix, oneStepProbabilities, *initialStates.begin(), flexibleBackwardTransitions, stateRewards); + if (!stateRewards) { + // If we are computing probabilities, then we can simply call the state elimination procedure. It + // will scale the transition row of the initial state with 1/(1-loopProbability). + STORM_PRINT_AND_LOG("Eliminating initial state " << *initialStates.begin() << "." << std::endl); + eliminateState(flexibleMatrix, oneStepProbabilities, *initialStates.begin(), flexibleBackwardTransitions, stateRewards); + } else { + // If we are computing rewards, we cannot call the state elimination procedure for technical reasons. + // Instead, we need to get rid of a potential loop in this state explicitly. + + // Start by finding the self-loop element. Since it can only be the only remaining outgoing transition + // of the initial state, this amounts to checking whether the outgoing transitions of the initial + // state are non-empty. + if (!flexibleMatrix.getRow(*initialStates.begin()).empty()) { + STORM_LOG_ASSERT(flexibleMatrix.getRow(*initialStates.begin()).size() == 1, "At most one outgoing transition expected at this point, but found more."); + STORM_LOG_ASSERT(flexibleMatrix.getRow(*initialStates.begin()).front().getColumn() == *initialStates.begin(), "Remaining entry should be a self-loop, but it is not."); + ValueType loopProbability = flexibleMatrix.getRow(*initialStates.begin()).front().getValue(); + loopProbability = storm::utility::constantOne() / (storm::utility::constantOne() - loopProbability); + loopProbability = storm::utility::pow(loopProbability, 2); + STORM_PRINT_AND_LOG("Scaling the transition reward of the initial state."); + stateRewards.get()[(*initialStates.begin())] *= loopProbability; + } + } // Make sure that we have eliminated all transitions from the initial state. STORM_LOG_ASSERT(flexibleMatrix.getRow(*initialStates.begin()).empty(), "The transitions of the initial states are non-empty."); @@ -98,7 +118,11 @@ namespace storm { } // Now, we return the value for the only initial state. - return storm::utility::simplify(oneStepProbabilities[*initialStates.begin()]); + if (stateRewards) { + return storm::utility::simplify(stateRewards.get()[*initialStates.begin()]); + } else { + return storm::utility::simplify(oneStepProbabilities[*initialStates.begin()]); + } } template @@ -134,12 +158,16 @@ namespace storm { // Create a vector for the probabilities to go to a state with probability 1 in one step. std::vector oneStepProbabilities = dtmc.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, psiStates); - // Create a vector that holds the one-step rewards. - std::vector oneStepRewards = std::vector(maybeStates.getNumberOfSetBits(), storm::utility::zero()); - // Project the state reward vector to all maybe-states. boost::optional> stateRewards(maybeStates.getNumberOfSetBits()); storm::utility::vector::selectVectorValues(stateRewards.get(), maybeStates, dtmc.getStateRewardVector()); + for (uint_fast64_t i = 0; i < dtmc.getStateRewardVector().size(); ++i) { + std::cout << i << ": " << dtmc.getStateRewardVector()[i] << ", "; + } + std::cout << std::endl; + std::cout << maybeStates << std::endl; + std::cout << psiStates << std::endl; + std::cout << infinityStates << std::endl; // Determine the set of initial states of the sub-DTMC. storm::storage::BitVector newInitialStates = dtmc.getInitialStates() % maybeStates; @@ -152,7 +180,7 @@ namespace storm { // impose ordering constraints later. std::vector statePriorities = getStatePriorities(submatrix, submatrixTransposed, newInitialStates, oneStepProbabilities); - return computeReachabilityValue(submatrix, oneStepRewards, submatrixTransposed, newInitialStates, phiStates, psiStates, stateRewards, statePriorities); + return computeReachabilityValue(submatrix, oneStepProbabilities, submatrixTransposed, newInitialStates, phiStates, psiStates, stateRewards, statePriorities); } template @@ -402,7 +430,9 @@ namespace storm { entry.setValue(storm::utility::simplify(entry.getValue() * loopProbability)); } } - oneStepProbabilities[state] = oneStepProbabilities[state] * loopProbability; + if (!stateRewards) { + oneStepProbabilities[state] = oneStepProbabilities[state] * loopProbability; + } } STORM_LOG_DEBUG((hasSelfLoop ? "State has self-loop." : "State does not have a self-loop.")); @@ -479,9 +509,15 @@ namespace storm { // Now move the new transitions in place. predecessorForwardTransitions = std::move(newSuccessors); - // Add the probabilities to go to a target state in just one step. - oneStepProbabilities[predecessor] += storm::utility::simplify(multiplyFactor * oneStepProbabilities[state]); - STORM_LOG_DEBUG("Fixed new next-state probabilities of predecessor states."); + if (!stateRewards) { + // Add the probabilities to go to a target state in just one step if we have to compute probabilities. + oneStepProbabilities[predecessor] += storm::utility::simplify(multiplyFactor * oneStepProbabilities[state]); + STORM_LOG_DEBUG("Fixed new next-state probabilities of predecessor states."); + } else { + // If we are computing rewards, we basically scale the state reward of the state to eliminate and + // add the result to the state reward of the predecessor. + stateRewards.get()[predecessor] += storm::utility::simplify(multiplyElement->getValue() * storm::utility::pow(loopProbability, 2) * stateRewards.get()[state]); + } } // Finally, we need to add the predecessor to the set of predecessors of every successor. diff --git a/src/storage/DeterministicModelBisimulationDecomposition.cpp b/src/storage/DeterministicModelBisimulationDecomposition.cpp index 930f6cda2..4c542a16e 100644 --- a/src/storage/DeterministicModelBisimulationDecomposition.cpp +++ b/src/storage/DeterministicModelBisimulationDecomposition.cpp @@ -640,6 +640,9 @@ namespace storm { // (b) the new labeling, // (c) the new reward structures. + std::cout << model.getTransitionMatrix() << std::endl; + std::cout << model.getLabeledStates("end") << std::endl;; + // Prepare a matrix builder for (a). storm::storage::SparseMatrixBuilder builder(this->size(), this->size()); @@ -739,6 +742,13 @@ namespace storm { // If the model has state rewards, we simply copy the state reward of the representative state, because // all states in a block are guaranteed to have the same state reward. if (keepRewards && model.hasStateRewards()) { + std::cout << "is block absorbing? " << oldBlock.isAbsorbing() << std::endl; + std::cout << "setting reward for block " << blockIndex << " to " << model.getStateRewardVector()[representativeState] << std::endl; + std::cout << "block: " << std::endl; + for (auto element : this->blocks[blockIndex]) { + std::cout << element << ", "; + } + std::cout << std::endl; stateRewards.get()[blockIndex] = model.getStateRewardVector()[representativeState]; } } @@ -751,6 +761,8 @@ namespace storm { // Finally construct the quotient model. this->quotient = std::shared_ptr>(new ModelType(builder.build(), std::move(newLabeling), std::move(stateRewards))); + + std::cout << quotient->getTransitionMatrix() << std::endl; } template @@ -1225,6 +1237,8 @@ namespace storm { this->initializeSilentProbabilities(model, partition); } + partition.print(); + return partition; } diff --git a/src/stormParametric.cpp b/src/stormParametric.cpp index d8d9d4f78..20a32cda1 100644 --- a/src/stormParametric.cpp +++ b/src/stormParametric.cpp @@ -235,7 +235,7 @@ int main(const int argc, const char** argv) { } } - STORM_LOG_ASSERT(parameters == valueFunction.gatherVariables(), "Parameters in result and program definition do not coincide."); +// STORM_LOG_ASSERT(parameters == valueFunction.gatherVariables(), "Parameters in result and program definition do not coincide."); if(storm::settings::parametricSettings().exportResultToFile()) { storm::utility::exportParametricMcResult(valueFunction, constraintCollector); diff --git a/src/utility/ConstantsComparator.cpp b/src/utility/ConstantsComparator.cpp index 9a0e4b31b..48c9a74ed 100644 --- a/src/utility/ConstantsComparator.cpp +++ b/src/utility/ConstantsComparator.cpp @@ -20,6 +20,11 @@ namespace storm { return std::numeric_limits::infinity(); } + template + ValueType pow(ValueType const& value, uint_fast64_t exponent) { + return std::pow(value, exponent); + } + template<> double simplify(double value) { // In the general case, we don't to anything here, but merely return the value. If something else is @@ -67,6 +72,11 @@ namespace storm { } #ifdef PARAMETRIC_SYSTEMS + template<> + RationalFunction pow(RationalFunction const& value, uint_fast64_t exponent) { + return carl::pow(exponent, value); + } + template<> RationalFunction simplify(RationalFunction value) { value.simplify(); @@ -116,8 +126,8 @@ namespace storm { bool ConstantsComparator::isConstant(storm::Polynomial const& value) const { return value.isConstant(); } - #endif + template storm::storage::MatrixEntry simplify(storm::storage::MatrixEntry matrixEntry) { simplify(matrixEntry.getValue()); @@ -141,6 +151,8 @@ namespace storm { template double one(); template double zero(); template double infinity(); + + template double pow(double const& value, uint_fast64_t exponent); template double simplify(double value); @@ -155,6 +167,8 @@ namespace storm { template RationalFunction one(); template RationalFunction zero(); + template RationalFunction pow(RationalFunction const& value, uint_fast64_t exponent); + template Polynomial one(); template Polynomial zero(); template RationalFunction simplify(RationalFunction value); @@ -165,5 +179,6 @@ namespace storm { template storm::storage::MatrixEntry& simplify(storm::storage::MatrixEntry& matrixEntry); template storm::storage::MatrixEntry&& simplify(storm::storage::MatrixEntry&& matrixEntry); #endif + } } \ No newline at end of file diff --git a/src/utility/ConstantsComparator.h b/src/utility/ConstantsComparator.h index 8c9ee9de7..00f2c14f9 100644 --- a/src/utility/ConstantsComparator.h +++ b/src/utility/ConstantsComparator.h @@ -27,6 +27,12 @@ namespace storm { template ValueType infinity(); + template + ValueType pow(ValueType const& value, uint_fast64_t exponent); + + template<> + RationalFunction pow(RationalFunction const& value, uint_fast64_t exponent); + template ValueType simplify(ValueType value); @@ -92,8 +98,6 @@ namespace storm { bool isConstant(storm::Polynomial const& value) const; }; - - #endif template