From b7492d543acaefd3f8cfe98bb863965f5a1667a4 Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Mon, 15 Dec 2014 17:32:09 +0100
Subject: [PATCH] Further work regarding rewards in parameterized models. Note:
 this includes some debug output.

Former-commit-id: ac65f020a578b21191bef3b8a21df713e4759c15
---
 src/adapters/ExplicitModelAdapter.h           |  1 +
 .../reachability/SparseSccModelChecker.cpp    | 58 +++++++++++++++----
 ...ministicModelBisimulationDecomposition.cpp | 14 +++++
 src/stormParametric.cpp                       |  2 +-
 src/utility/ConstantsComparator.cpp           | 17 +++++-
 src/utility/ConstantsComparator.h             |  8 ++-
 6 files changed, 85 insertions(+), 15 deletions(-)

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<ValueType> buildStateRewards(std::vector<storm::prism::StateReward> const& rewards, StateInformation const& stateInformation, expressions::ExpressionEvaluation<ValueType>& eval) {
                 std::vector<ValueType> 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<ValueType>() / (storm::utility::constantOne<ValueType>() - 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<typename ValueType>
@@ -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<ValueType> oneStepProbabilities = dtmc.getTransitionMatrix().getConstrainedRowSumVector(maybeStates, psiStates);
                 
-                // Create a vector that holds the one-step rewards.
-                std::vector<ValueType> oneStepRewards = std::vector<ValueType>(maybeStates.getNumberOfSetBits(), storm::utility::zero<ValueType>());
-                
                 // Project the state reward vector to all maybe-states.
                 boost::optional<std::vector<ValueType>> 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<std::size_t> 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<typename ValueType>
@@ -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<ValueType> 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<storm::models::AbstractDeterministicModel<ValueType>>(new ModelType(builder.build(), std::move(newLabeling), std::move(stateRewards)));
+            
+            std::cout << quotient->getTransitionMatrix() << std::endl;
         }
         
         template<typename ValueType>
@@ -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<ValueType>::infinity();
         }
         
+        template<typename ValueType>
+        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<storm::Polynomial>::isConstant(storm::Polynomial const& value) const {
             return value.isConstant();
         }
-        
 #endif
+
         template<typename IndexType, typename ValueType>
         storm::storage::MatrixEntry<IndexType, ValueType> simplify(storm::storage::MatrixEntry<IndexType, ValueType> 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<storm::storage::sparse::state_type, RationalFunction>& simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, RationalFunction>& matrixEntry);
         template storm::storage::MatrixEntry<storm::storage::sparse::state_type, RationalFunction>&& simplify(storm::storage::MatrixEntry<storm::storage::sparse::state_type, RationalFunction>&& 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<typename ValueType>
         ValueType infinity();
         
+        template<typename ValueType>
+        ValueType pow(ValueType const& value, uint_fast64_t exponent);
+        
+        template<>
+        RationalFunction pow(RationalFunction const& value, uint_fast64_t exponent);
+
         template<typename ValueType>
         ValueType simplify(ValueType value);
 
@@ -92,8 +98,6 @@ namespace storm {
             
             bool isConstant(storm::Polynomial const& value) const;
         };
-        
-        
 #endif
         
         template<typename IndexType, typename ValueType>