sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

some quick and dirty optimizations (that need to be reviewed) 

Former-commit-id: 86fbc66c90
main
TimQu 9 years ago
parent
b362047e4f
commit
71e181bd93
27 changed files with 314 additions and 431 deletions
Split View
  1. 4
      examples/multiobjective/ma/polling/polling.ma
  2. 16
      examples/multiobjective/ma/polling/untitled.txt
  3. 21
      examples/multiobjective/ma/stream/stream.ma
  4. 50
      examples/multiobjective/ma/stream/stream10.ma
  5. 50
      examples/multiobjective/ma/stream/stream100.ma
  6. 50
      examples/multiobjective/ma/stream/stream1000.ma
  7. 50
      examples/multiobjective/ma/stream/stream2.ma
  8. 50
      examples/multiobjective/ma/stream/stream50.ma
  9. 50
      examples/multiobjective/ma/stream/stream500.ma
  10. 50
      examples/multiobjective/ma/stream/stream5000.ma
  11. 1
      resources/3rdparty/carl
  12. 2
      src/logic/Bound.h
  13. 5
      src/logic/OperatorFormula.h
  14. 54
      src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp
  15. 148
      src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp
  16. 2
      src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h
  17. 40
      src/modelchecker/multiobjective/helper/SparseMultiObjectivePostprocessor.cpp
  18. 5
      src/modelchecker/multiobjective/helper/SparseMultiObjectiveResultData.h
  19. 30
      src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp
  20. 1
      src/models/sparse/MarkovAutomaton.cpp
  21. 19
      src/solver/GmmxxLinearEquationSolver.cpp
  22. 1
      src/solver/GmmxxLinearEquationSolver.h
  23. 12
      src/solver/StandardMinMaxLinearEquationSolver.cpp
  24. 3
      src/solver/StandardMinMaxLinearEquationSolver.h
  25. 6
      src/storage/SparseMatrix.cpp
  26. 21
      src/utility/vector.h
  27. 4
      test/functional/modelchecker/SparseMdpMultiObjectiveModelCheckerTest.cpp

4
examples/multiobjective/ma/polling/polling.ma
View File

@ -76,11 +76,11 @@ label "allqueuesfull" = q1Size=Q & q2Size=Q;
// Rewards adapted from Guck et al.: Modelling and Analysis of Markov Reward Automata
rewards "processedjobs1"
[copy1] true : 1;
[copy1] true : 0.1;
endrewards
rewards "processedjobs2"
[copy1] true : 1;
[copy1] true : 0.1;
endrewards
rewards "processedjobs"

16
examples/multiobjective/ma/polling/untitled.txt
View File

@ -0,0 +1,16 @@
mdp
const int J = 0;
module test
x : [0..1];
j : [0..1] init 1;
z : [0..1];
[] j=1 -> 1 : (x'=j);
[] j<1 -> 1 : (x'=0);
[] x=1 -> (z'=1);
endmodule

21
examples/multiobjective/ma/stream/stream.ma
View File

@ -4,7 +4,7 @@ ma
const int N; // num packages
const double inRate = 4;
const double processingRate = 5;
const double processingRate = 4;
module streamingclient
@ -12,12 +12,13 @@ module streamingclient
// 0: decide whether to start
// 1: buffering
// 2: running
// 3: done
// 3: success
n : [0..N]; // number of received packages
k : [0..N]; // number of processed packages
[buffer] s=0 & n<N -> 1 : (s'=1);
[buffer] s=0 & n<N & k=n -> 1 : (s'=1);
[buffer] s=0 & n<N & k<n -> 0.99: (s'=1) + 0.01 : (s'=2) & (k'=k+1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
@ -30,21 +31,15 @@ module streamingclient
<> s=3 -> 1 : true;
endmodule
// All packages received and buffer empty
label "underrun" = (s=0 & k>0);
label "running" = (s=2);
label "done" = (s=3);
rewards "buffering"
s=1 : 1;
endrewards
rewards "initialbuffering"
s=1 & k = 0: 1;
endrewards
rewards "intermediatebuffering"
s=1 & k > 0: 1;
endrewards
rewards "numRestarts"
rewards "numrestarts"
[start] k > 0 : 1;
endrewards

50
examples/multiobjective/ma/stream/stream10.ma
View File

@ -1,50 +0,0 @@
ma
const int N=10; // num packages
const double inRate = 4;
const double processingRate = 5;
module streamingclient
s : [0..3]; // current state:
// 0: decide whether to start
// 1: buffering
// 2: running
// 3: done
n : [0..N]; // number of received packages
k : [0..N]; // number of processed packages
[buffer] s=0 & n<N -> 1 : (s'=1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
<> s=2 & n<N & k<n -> inRate : (n'=n+1) + processingRate : (k'=k+1);
<> s=2 & n<N & k=n -> inRate : (n'=n+1) + processingRate : (s'=0);
<> s=2 & n=N & k<n -> processingRate : (k'=k+1);
<> s=2 & n=N & k=N -> processingRate : (s'=3);
<> s=3 -> 1 : true;
endmodule
// All packages received and buffer empty
label "done" = (s=3);
rewards "buffering"
s=1 : 1;
endrewards
rewards "initialbuffering"
s=1 & k = 0: 1;
endrewards
rewards "intermediatebuffering"
s=1 & k > 0: 1;
endrewards
rewards "numRestarts"
[start] k > 0 : 1;
endrewards

50
examples/multiobjective/ma/stream/stream100.ma
View File

@ -1,50 +0,0 @@
ma
const int N=100; // num packages
const double inRate = 4;
const double processingRate = 5;
module streamingclient
s : [0..3]; // current state:
// 0: decide whether to start
// 1: buffering
// 2: running
// 3: done
n : [0..N]; // number of received packages
k : [0..N]; // number of processed packages
[buffer] s=0 & n<N -> 1 : (s'=1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
<> s=2 & n<N & k<n -> inRate : (n'=n+1) + processingRate : (k'=k+1);
<> s=2 & n<N & k=n -> inRate : (n'=n+1) + processingRate : (s'=0);
<> s=2 & n=N & k<n -> processingRate : (k'=k+1);
<> s=2 & n=N & k=N -> processingRate : (s'=3);
<> s=3 -> 1 : true;
endmodule
// All packages received and buffer empty
label "done" = (s=3);
rewards "buffering"
s=1 : 1;
endrewards
rewards "initialbuffering"
s=1 & k = 0: 1;
endrewards
rewards "intermediatebuffering"
s=1 & k > 0: 1;
endrewards
rewards "numRestarts"
[start] k > 0 : 1;
endrewards

50
examples/multiobjective/ma/stream/stream1000.ma
View File

@ -1,50 +0,0 @@
ma
const int N=1000; // num packages
const double inRate = 4;
const double processingRate = 5;
module streamingclient
s : [0..3]; // current state:
// 0: decide whether to start
// 1: buffering
// 2: running
// 3: done
n : [0..N]; // number of received packages
k : [0..N]; // number of processed packages
[buffer] s=0 & n<N -> 1 : (s'=1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
<> s=2 & n<N & k<n -> inRate : (n'=n+1) + processingRate : (k'=k+1);
<> s=2 & n<N & k=n -> inRate : (n'=n+1) + processingRate : (s'=0);
<> s=2 & n=N & k<n -> processingRate : (k'=k+1);
<> s=2 & n=N & k=N -> processingRate : (s'=3);
<> s=3 -> 1 : true;
endmodule
// All packages received and buffer empty
label "done" = (s=3);
rewards "buffering"
s=1 : 1;
endrewards
rewards "initialbuffering"
s=1 & k = 0: 1;
endrewards
rewards "intermediatebuffering"
s=1 & k > 0: 1;
endrewards
rewards "numRestarts"
[start] k > 0 : 1;
endrewards

50
examples/multiobjective/ma/stream/stream2.ma
View File

@ -1,50 +0,0 @@
ma
const int N = 2; // num packages
const double inRate = 4;
const double processingRate = 5;
module streamingclient
s : [0..3]; // current state:
// 0: decide whether to start
// 1: buffering
// 2: running
// 3: done
n : [0..N]; // number of received packages
k : [0..N]; // number of processed packages
[buffer] s=0 & n<N -> 1 : (s'=1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
<> s=2 & n<N & k<n -> inRate : (n'=n+1) + processingRate : (k'=k+1);
<> s=2 & n<N & k=n -> inRate : (n'=n+1) + processingRate : (s'=0);
<> s=2 & n=N & k<n -> processingRate : (k'=k+1);
<> s=2 & n=N & k=N -> processingRate : (s'=3);
<> s=3 -> 1 : true;
endmodule
// All packages received and buffer empty
label "done" = (s=3);
rewards "buffering"
s=1 : 1;
endrewards
rewards "initialbuffering"
s=1 & k = 0: 1;
endrewards
rewards "intermediatebuffering"
s=1 & k > 0: 1;
endrewards
rewards "numRestarts"
[start] k > 0 : 1;
endrewards

50
examples/multiobjective/ma/stream/stream50.ma
View File

@ -1,50 +0,0 @@
ma
const int N=50; // num packages
const double inRate = 4;
const double processingRate = 5;
module streamingclient
s : [0..3]; // current state:
// 0: decide whether to start
// 1: buffering
// 2: running
// 3: done
n : [0..N]; // number of received packages
k : [0..N]; // number of processed packages
[buffer] s=0 & n<N -> 1 : (s'=1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
<> s=2 & n<N & k<n -> inRate : (n'=n+1) + processingRate : (k'=k+1);
<> s=2 & n<N & k=n -> inRate : (n'=n+1) + processingRate : (s'=0);
<> s=2 & n=N & k<n -> processingRate : (k'=k+1);
<> s=2 & n=N & k=N -> processingRate : (s'=3);
<> s=3 -> 1 : true;
endmodule
// All packages received and buffer empty
label "done" = (s=3);
rewards "buffering"
s=1 : 1;
endrewards
rewards "initialbuffering"
s=1 & k = 0: 1;
endrewards
rewards "intermediatebuffering"
s=1 & k > 0: 1;
endrewards
rewards "numRestarts"
[start] k > 0 : 1;
endrewards

50
examples/multiobjective/ma/stream/stream500.ma
View File

@ -1,50 +0,0 @@
ma
const int N=500; // num packages
const double inRate = 4;
const double processingRate = 5;
module streamingclient
s : [0..3]; // current state:
// 0: decide whether to start
// 1: buffering
// 2: running
// 3: done
n : [0..N]; // number of received packages
k : [0..N]; // number of processed packages
[buffer] s=0 & n<N -> 1 : (s'=1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
<> s=2 & n<N & k<n -> inRate : (n'=n+1) + processingRate : (k'=k+1);
<> s=2 & n<N & k=n -> inRate : (n'=n+1) + processingRate : (s'=0);
<> s=2 & n=N & k<n -> processingRate : (k'=k+1);
<> s=2 & n=N & k=N -> processingRate : (s'=3);
<> s=3 -> 1 : true;
endmodule
// All packages received and buffer empty
label "done" = (s=3);
rewards "buffering"
s=1 : 1;
endrewards
rewards "initialbuffering"
s=1 & k = 0: 1;
endrewards
rewards "intermediatebuffering"
s=1 & k > 0: 1;
endrewards
rewards "numRestarts"
[start] k > 0 : 1;
endrewards

50
examples/multiobjective/ma/stream/stream5000.ma
View File

@ -1,50 +0,0 @@
ma
const int N=5000; // num packages
const double inRate = 4;
const double processingRate = 5;
module streamingclient
s : [0..3]; // current state:
// 0: decide whether to start
// 1: buffering
// 2: running
// 3: done
n : [0..N]; // number of received packages
k : [0..N]; // number of processed packages
[buffer] s=0 & n<N -> 1 : (s'=1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
<> s=2 & n<N & k<n -> inRate : (n'=n+1) + processingRate : (k'=k+1);
<> s=2 & n<N & k=n -> inRate : (n'=n+1) + processingRate : (s'=0);
<> s=2 & n=N & k<n -> processingRate : (k'=k+1);
<> s=2 & n=N & k=N -> processingRate : (s'=3);
<> s=3 -> 1 : true;
endmodule
// All packages received and buffer empty
label "done" = (s=3);
rewards "buffering"
s=1 : 1;
endrewards
rewards "initialbuffering"
s=1 & k = 0: 1;
endrewards
rewards "intermediatebuffering"
s=1 & k > 0: 1;
endrewards
rewards "numRestarts"
[start] k > 0 : 1;
endrewards

1
resources/3rdparty/carl

@ -1 +0,0 @@
Subproject commit d67f986226cf846ba6366cdbe2abc21dff375542

2
src/logic/Bound.h
View File

@ -26,7 +26,7 @@ namespace storm {
template<typename ValueType>
std::ostream& operator<<(std::ostream& out, Bound<ValueType> const& bound) {
out << bound.comparisonType << bound.threshold;
out << bound.comparisonType << storm::utility::convertNumber<double>(bound.threshold);
return out;
}
}

5
src/logic/OperatorFormula.h
View File

@ -47,6 +47,9 @@ namespace storm {
virtual bool isOperatorFormula() const override;
OperatorInformation const& getOperatorInformation() const;
void setOperatorInformation(OperatorInformation newOperatorInformation) {
this->operatorInformation = newOperatorInformation;
}
virtual bool hasQualitativeResult() const override;
virtual bool hasQuantitativeResult() const override;
@ -59,4 +62,4 @@ namespace storm {
}
}
#endif /* STORM_LOGIC_OPERATORFORMULA_H_ */
#endif /* STORM_LOGIC_OPERATORFORMULA_H_ */

54
src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp
View File

@ -8,9 +8,11 @@
#include "src/transformer/EndComponentEliminator.h"
#include "src/utility/macros.h"
#include "src/utility/vector.h"
#include "src/solver/GmmxxLinearEquationSolver.h"
#include "src/exceptions/InvalidOperationException.h"
namespace storm {
namespace modelchecker {
namespace helper {
@ -67,6 +69,8 @@ namespace storm {
auto lowerTimeBoundIt = lowerTimeBounds.begin();
auto upperTimeBoundIt = upperTimeBounds.begin();
uint_fast64_t currentEpoch = std::max(lowerTimeBounds.empty() ? 0 : lowerTimeBoundIt->first, upperTimeBounds.empty() ? 0 : upperTimeBoundIt->first);
std::cout << "Checking " << currentEpoch << " epochs for current weight vector" << std::endl;
while(true) {
// Update the objectives that are considered at the current time epoch as well as the (weighted) reward vectors.
updateDataToCurrentEpoch(MS, PS, *minMax, consideredObjectives, currentEpoch, weightVector, lowerTimeBoundIt, lowerTimeBounds, upperTimeBoundIt, upperTimeBounds);
@ -78,10 +82,6 @@ namespace storm {
// Only perform such a step if there is time left.
if(currentEpoch>0) {
performMSStep(MS, PS, consideredObjectives);
// if(currentEpoch % 1000 == 0) {
// STORM_LOG_DEBUG(currentEpoch << " digitized time steps left. Current weighted value is " << PS.weightedSolutionVector[0]);
// std::cout << std::endl << currentEpoch << " digitized time steps left. Current weighted value is " << PS.weightedSolutionVector[0] << std::endl;
// }
--currentEpoch;
} else {
break;
@ -365,27 +365,47 @@ namespace storm {
template <class SparseMaModelType>
void SparseMaMultiObjectiveWeightVectorChecker<SparseMaModelType>::performPSStep(SubModel& PS, SubModel const& MS, MinMaxSolverData& minMax, LinEqSolverData& linEq, std::vector<uint_fast64_t>& optimalChoicesAtCurrentEpoch, storm::storage::BitVector const& consideredObjectives) const {
// compute a choice vector for the probabilistic states that is optimal w.r.t. the weighted reward vector
std::vector<uint_fast64_t> newOptimalChoices(PS.getNumberOfStates());
minMax.solver->solveEquations(minMax.x, minMax.b);
this->transformReducedSolutionToOriginalModel(minMax.matrix, minMax.x, minMax.solver->getScheduler()->getChoices(), minMax.toPSChoiceMapping, minMax.fromPSStateMapping, PS.toPS, PS.weightedSolutionVector, newOptimalChoices);
// check whether the linEqSolver needs to be updated, i.e., whether the scheduler has changed
if(newOptimalChoices != optimalChoicesAtCurrentEpoch) {
// std::cout << "X";
if(linEq.solver == nullptr || newOptimalChoices != optimalChoicesAtCurrentEpoch) {
optimalChoicesAtCurrentEpoch.swap(newOptimalChoices);
linEq.solver = nullptr;
storm::storage::SparseMatrix<ValueType> linEqMatrix = PS.toPS.selectRowsFromRowGroups(optimalChoicesAtCurrentEpoch, true);
linEqMatrix.convertToEquationSystem();
linEq.solver = linEq.factory.create(std::move(linEqMatrix));
} else {
// std::cout << " ";
// For a gmm solver, switch the method to jacobi since this is 'usually' faster.
storm::solver::GmmxxLinearEquationSolver<ValueType>* gmmSolver = dynamic_cast<storm::solver::GmmxxLinearEquationSolver<ValueType>*>(linEq.solver.get());
if(gmmSolver!=nullptr) {
STORM_LOG_INFO("Switching to Jacobi method");
gmmSolver->getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi);
gmmSolver->getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::None);
}
linEq.solver->allocateAuxMemory(storm::solver::LinearEquationSolverOperation::SolveEquations);
}
// Get the results for the individual objectives.
// Note that we do not consider an estimate for each objective (as done in the unbounded phase) since the results from the previous epoch are already pretty close
for(auto objIndex : consideredObjectives) {
PS.toMS.multiplyWithVector(MS.objectiveSolutionVectors[objIndex], PS.auxChoiceValues);
storm::utility::vector::addVectors(PS.auxChoiceValues, PS.objectiveRewardVectors[objIndex], PS.auxChoiceValues);
storm::utility::vector::selectVectorValues(linEq.b, optimalChoicesAtCurrentEpoch, PS.toPS.getRowGroupIndices(), PS.auxChoiceValues);
auto const& objectiveRewardVectorPS = PS.objectiveRewardVectors[objIndex];
auto const& objectiveSolutionVectorMS = MS.objectiveSolutionVectors[objIndex];
// compute rhs of equation system, i.e., PS.toMS * x + Rewards
// To safe some time, only do this for the obtained optimal choices
auto itGroupIndex = PS.toPS.getRowGroupIndices().begin();
auto itChoiceOffset = optimalChoicesAtCurrentEpoch.begin();
for(auto& bValue : linEq.b) {
uint_fast64_t row = (*itGroupIndex) + (*itChoiceOffset);
bValue = objectiveRewardVectorPS[row];
for(auto const& entry : PS.toMS.getRow(row)){
bValue += entry.getValue() * objectiveSolutionVectorMS[entry.getColumn()];
}
++itGroupIndex;
++itChoiceOffset;
}
linEq.solver->solveEquations(PS.objectiveSolutionVectors[objIndex], linEq.b);
}
}
@ -412,10 +432,10 @@ namespace storm {
template void SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::digitize<double>(SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::SubModel& subModel, double const& digitizationConstant) const;
template void SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::digitizeTimeBounds<double>(SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& lowerTimeBounds, SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& upperTimeBounds, double const& digitizationConstant);
#ifdef STORM_HAVE_CARL
template class SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>;
template storm::RationalNumber SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::getDigitizationConstant<storm::RationalNumber>() const;
template void SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::digitize<storm::RationalNumber>(SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const;
template void SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::digitizeTimeBounds<storm::RationalNumber>(SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& lowerTimeBounds, SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant);
// template class SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>;
// template storm::RationalNumber SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::getDigitizationConstant<storm::RationalNumber>() const;
// template void SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::digitize<storm::RationalNumber>(SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const;
// template void SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>::digitizeTimeBounds<storm::RationalNumber>(SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& lowerTimeBounds, SparseMaMultiObjectiveWeightVectorChecker<storm::models::sparse::MarkovAutomaton<double>>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant);
#endif
}

148
src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp
View File

@ -13,6 +13,8 @@
#include "src/exceptions/UnexpectedException.h"
#include "src/logic/CloneVisitor.h"
namespace storm {
namespace modelchecker {
namespace helper {
@ -23,6 +25,17 @@ namespace storm {
ResultData resultData;
resultData.overApproximation() = storm::storage::geometry::Polytope<RationalNumberType>::createUniversalPolytope();
resultData.underApproximation() = storm::storage::geometry::Polytope<RationalNumberType>::createEmptyPolytope();
resultData.stopWatchWeightVectorChecker.pause();
resultData.stopWatchWeightVectorChecker.reset();
// Set the maximum gap between lower and upper bound of the weightVectorChecker result.
// This is the maximal edge length of the box we have to consider around each computed point
// We pick the gap such that the maximal distance between two points within this box is less than the given precision divided by two.
SparseModelValueType gap = storm::utility::convertNumber<SparseModelValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
gap /= (storm::utility::one<SparseModelValueType>() + storm::utility::one<SparseModelValueType>());
gap /= storm::utility::sqrt(static_cast<SparseModelValueType>(preprocessorData.objectives.size()));
weightVectorChecker->setMaximumLowerUpperBoundGap(gap);
if(!checkIfPreprocessingWasConclusive(preprocessorData)) {
switch(preprocessorData.queryType) {
case PreprocessorData::QueryType::Achievability:
@ -57,8 +70,6 @@ namespace storm {
template <class SparseModelType, typename RationalNumberType>
void SparseMultiObjectiveHelper<SparseModelType, RationalNumberType>::achievabilityQuery(PreprocessorData const& preprocessorData, WeightVectorCheckerType weightVectorChecker, ResultData& resultData) {
//Set the maximum gap between lower and upper bound of the weightVectorChecker result we initially allow for this query type. Note that we could also try other values
weightVectorChecker->setMaximumLowerUpperBoundGap(storm::utility::convertNumber<SparseModelValueType>(0.1));
// Get a point that represents the thresholds
Point thresholds;
thresholds.reserve(preprocessorData.objectives.size());
@ -84,8 +95,6 @@ namespace storm {
template <class SparseModelType, typename RationalNumberType>
void SparseMultiObjectiveHelper<SparseModelType, RationalNumberType>::numericalQuery(PreprocessorData const& preprocessorData, WeightVectorCheckerType weightVectorChecker, ResultData& resultData) {
STORM_LOG_ASSERT(preprocessorData.indexOfOptimizingObjective, "Detected numerical query but index of optimizing objective is not set.");
// Set the maximum gap between lower and upper bound of the weightVectorChecker result we initially allow for this query type. Note that we could also try other values
weightVectorChecker->setMaximumLowerUpperBoundGap(storm::utility::convertNumber<SparseModelValueType>(0.1));
// initialize some data
uint_fast64_t optimizingObjIndex = *preprocessorData.indexOfOptimizingObjective;
Point thresholds;
@ -131,12 +140,6 @@ namespace storm {
// Note that we do not have to care whether a threshold is strict anymore, because the resulting optimum should be
// the supremum over all strategies. Hence, one could combine a scheduler inducing the optimum value (but possibly violating strict
// thresholds) and (with very low probability) a scheduler that satisfies all (possibly strict) thresholds.
// Pick the maximum gap between lower and upper bound of the weightVectorChecker result.
SparseModelValueType gap = storm::utility::convertNumber<SparseModelValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
gap /= (storm::utility::one<SparseModelValueType>() + storm::utility::one<SparseModelValueType>());
weightVectorChecker->setMaximumLowerUpperBoundGap(gap);
while(true) {
std::pair<Point, bool> optimizationRes = resultData.underApproximation()->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective);
STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "The underapproximation is either unbounded or empty.");
@ -148,6 +151,7 @@ namespace storm {
if(optimizationRes.second) {
resultData.setPrecisionOfResult(optimizationRes.first[optimizingObjIndex] - thresholds[optimizingObjIndex]);
STORM_LOG_DEBUG("Solution can be improved by at most " << resultData.template getPrecisionOfResult<double>());
std::cout << "Current precision of the numerical result is " << resultData.template getPrecisionOfResult<double>() << std::endl;
}
if(targetPrecisionReached(resultData) || maxStepsPerformed(resultData)) {
resultData.setOptimumIsAchievable(checkIfThresholdsAreSatisfied(resultData.underApproximation(), thresholds, strictThresholds));
@ -161,19 +165,63 @@ namespace storm {
template <class SparseModelType, typename RationalNumberType>
void SparseMultiObjectiveHelper<SparseModelType, RationalNumberType>::paretoQuery(PreprocessorData const& preprocessorData, WeightVectorCheckerType weightVectorChecker, ResultData& resultData) {
// Set the maximum gap between lower and upper bound of the weightVectorChecker result.
// This is the maximal edge length of the box we have to consider around each computed point
// We pick the gap such that the maximal distance between two points within this box is less than the given precision divided by two.
SparseModelValueType gap = storm::utility::convertNumber<SparseModelValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
gap /= (storm::utility::one<SparseModelValueType>() + storm::utility::one<SparseModelValueType>());
gap /= storm::utility::sqrt(static_cast<SparseModelValueType>(preprocessorData.objectives.size()));
weightVectorChecker->setMaximumLowerUpperBoundGap(gap);
//First consider the objectives individually
for(uint_fast64_t objIndex = 0; objIndex<preprocessorData.objectives.size() && !maxStepsPerformed(resultData); ++objIndex) {
WeightVector direction(preprocessorData.objectives.size(), storm::utility::zero<RationalNumberType>());
direction[objIndex] = storm::utility::one<RationalNumberType>();
performRefinementStep(std::move(direction), preprocessorData.produceSchedulers, weightVectorChecker, resultData);
//TODO remove this code..
std::cout << "multi(" << preprocessorData.originalFormula.getSubformula(0);
for(uint_fast64_t subFIndex = 1; subFIndex < preprocessorData.originalFormula.getNumberOfSubformulas(); ++subFIndex) {
std::cout << ", ";
RationalNumberType thresholdValue =
storm::utility::convertNumber<RationalNumberType>(resultData.refinementSteps().back().getLowerBoundPoint()[subFIndex]) * storm::utility::convertNumber<RationalNumberType>(preprocessorData.objectives[subFIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber<RationalNumberType>(preprocessorData.objectives[subFIndex].toOriginalValueTransformationOffset);
auto const& subF = preprocessorData.originalFormula.getSubformula(subFIndex);
if(subF.isOperatorFormula() && subF.asOperatorFormula().hasOptimalityType()) {
storm::logic::CloneVisitor cv;
auto copySubF = cv.clone(subF);
if(subF.asOperatorFormula().getOptimalityType() == storm::solver::OptimizationDirection::Maximize) {
storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound<storm::RationalNumber>(storm::logic::ComparisonType::GreaterEqual, thresholdValue));
copySubF->asOperatorFormula().setOperatorInformation(opInfo);
} else {
storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound<storm::RationalNumber>(storm::logic::ComparisonType::LessEqual, thresholdValue));
copySubF->asOperatorFormula().setOperatorInformation(opInfo);
}
std::cout << *copySubF;
}
}
std::cout << ")" << std::endl;
//TODO remove this code..
std::cout << "multi(";
for(uint_fast64_t subFIndex = 0; subFIndex < preprocessorData.originalFormula.getNumberOfSubformulas(); ++subFIndex) {
if(subFIndex != 0) {
std::cout << ", ";
}
RationalNumberType thresholdValue =
storm::utility::convertNumber<RationalNumberType>(resultData.refinementSteps().back().getLowerBoundPoint()[subFIndex]) * storm::utility::convertNumber<RationalNumberType>(preprocessorData.objectives[subFIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber<RationalNumberType>(preprocessorData.objectives[subFIndex].toOriginalValueTransformationOffset);
auto const& subF = preprocessorData.originalFormula.getSubformula(subFIndex);
if(subF.isOperatorFormula() && subF.asOperatorFormula().hasOptimalityType()) {
storm::logic::CloneVisitor cv;
auto copySubF = cv.clone(subF);
if(subF.asOperatorFormula().getOptimalityType() == storm::solver::OptimizationDirection::Maximize) {
storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound<storm::RationalNumber>(storm::logic::ComparisonType::GreaterEqual, thresholdValue));
copySubF->asOperatorFormula().setOperatorInformation(opInfo);
} else {
storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound<storm::RationalNumber>(storm::logic::ComparisonType::LessEqual, thresholdValue));
copySubF->asOperatorFormula().setOperatorInformation(opInfo);
}
std::cout << *copySubF;
}
}
std::cout << ")" << std::endl;
}
while(true) {
@ -192,6 +240,65 @@ namespace storm {
}
}
resultData.setPrecisionOfResult(farestDistance);
std::cout << "Current precision of the approximation of the pareto curve is " << resultData.template getPrecisionOfResult<double>() << std::endl;
//TODO remove this code..
std::cout << "multi(" << preprocessorData.originalFormula.getSubformula(0);
for(uint_fast64_t subFIndex = 1; subFIndex < preprocessorData.originalFormula.getNumberOfSubformulas(); ++subFIndex) {
std::cout << ", ";
RationalNumberType thresholdValue =
storm::utility::convertNumber<RationalNumberType>(resultData.refinementSteps().back().getLowerBoundPoint()[subFIndex]) * storm::utility::convertNumber<RationalNumberType>(preprocessorData.objectives[subFIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber<RationalNumberType>(preprocessorData.objectives[subFIndex].toOriginalValueTransformationOffset);
auto const& subF = preprocessorData.originalFormula.getSubformula(subFIndex);
if(subF.isOperatorFormula() && subF.asOperatorFormula().hasOptimalityType()) {
storm::logic::CloneVisitor cv;
auto copySubF = cv.clone(subF);
if(subF.asOperatorFormula().getOptimalityType() == storm::solver::OptimizationDirection::Maximize) {
storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound<storm::RationalNumber>(storm::logic::ComparisonType::GreaterEqual, thresholdValue));
copySubF->asOperatorFormula().setOperatorInformation(opInfo);
} else {
storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound<storm::RationalNumber>(storm::logic::ComparisonType::LessEqual, thresholdValue));
copySubF->asOperatorFormula().setOperatorInformation(opInfo);
}
std::cout << *copySubF;
}
}
std::cout << ")" << std::endl;
//TODO remove this code..
std::cout << "multi(";
for(uint_fast64_t subFIndex = 0; subFIndex < preprocessorData.originalFormula.getNumberOfSubformulas(); ++subFIndex) {
if(subFIndex != 0) {
std::cout << ", ";
}
RationalNumberType thresholdValue =
storm::utility::convertNumber<RationalNumberType>(resultData.refinementSteps().back().getLowerBoundPoint()[subFIndex]) * storm::utility::convertNumber<RationalNumberType>(preprocessorData.objectives[subFIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber<RationalNumberType>(preprocessorData.objectives[subFIndex].toOriginalValueTransformationOffset);
auto const& subF = preprocessorData.originalFormula.getSubformula(subFIndex);
if(subF.isOperatorFormula() && subF.asOperatorFormula().hasOptimalityType()) {
storm::logic::CloneVisitor cv;
auto copySubF = cv.clone(subF);
if(subF.asOperatorFormula().getOptimalityType() == storm::solver::OptimizationDirection::Maximize) {
storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound<storm::RationalNumber>(storm::logic::ComparisonType::GreaterEqual, thresholdValue));
copySubF->asOperatorFormula().setOperatorInformation(opInfo);
} else {
storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound<storm::RationalNumber>(storm::logic::ComparisonType::LessEqual, thresholdValue));
copySubF->asOperatorFormula().setOperatorInformation(opInfo);
}
std::cout << *copySubF;
}
}
std::cout << ")" << std::endl;
STORM_LOG_DEBUG("Current precision of the approximation of the pareto curve is " << resultData.template getPrecisionOfResult<double>());
if(targetPrecisionReached(resultData) || maxStepsPerformed(resultData)) {
break;
@ -264,7 +371,11 @@ namespace storm {
break;
}
}
result.stopWatchWeightVectorChecker.unpause();
weightVectorChecker->check(storm::utility::vector::convertNumericVector<SparseModelValueType>(direction));
result.stopWatchWeightVectorChecker.pause();
if(oldMaximumLowerUpperBoundGap) {
// Reset the precision back to the previous values
weightVectorChecker->setMaximumLowerUpperBoundGap(*oldMaximumLowerUpperBoundGap);
@ -283,6 +394,7 @@ namespace storm {
}
updateOverApproximation(result.refinementSteps(), result.overApproximation());
updateUnderApproximation(result.refinementSteps(), result.underApproximation());
std::cout << "number of performed refinement Steps is: " << result.refinementSteps().size() << std::endl;
}
template <class SparseModelType, typename RationalNumberType>

2
src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h
View File

@ -52,7 +52,7 @@ namespace storm {
static WeightVector findSeparatingVector(Point const& pointToBeSeparated, std::shared_ptr<storm::storage::geometry::Polytope<RationalNumberType>> const& underApproximation, storm::storage::BitVector& individualObjectivesToBeChecked);
/*
* Refines the current result w.r.t. the given direction vector
* Refines the current result w.r.t. the given direction vector.
*/
static void performRefinementStep(WeightVector&& direction, bool saveScheduler, WeightVectorCheckerType weightVectorChecker, ResultData& resultData);

40
src/modelchecker/multiobjective/helper/SparseMultiObjectivePostprocessor.cpp
View File

@ -14,6 +14,7 @@
#include "src/settings//SettingsManager.h"
#include "src/settings/modules/MultiObjectiveSettings.h"
#include "src/settings/modules/GeneralSettings.h"
#include "src/settings/modules/IOSettings.h"
#include "src/utility/export.h"
#include "src/utility/macros.h"
#include "src/utility/vector.h"
@ -40,7 +41,8 @@ namespace storm {
STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unknown Query Type");
}
STORM_LOG_INFO(getInfo(result, preprocessorData, resultData, preprocessorStopwatch, helperStopwatch, postprocessorStopwatch));
//STORM_LOG_INFO(getInfo(result, preprocessorData, resultData, preprocessorStopwatch, helperStopwatch, postprocessorStopwatch));
std:: cout << getInfo(result, preprocessorData, resultData, preprocessorStopwatch, helperStopwatch, postprocessorStopwatch);
exportPlot(result, preprocessorData, resultData, preprocessorStopwatch, helperStopwatch, postprocessorStopwatch);
@ -191,7 +193,12 @@ namespace storm {
template<typename SparseModelType, typename RationalNumberType>
std::shared_ptr<storm::storage::geometry::Polytope<RationalNumberType>> SparseMultiObjectivePostprocessor<SparseModelType, RationalNumberType>::transformToOriginalValues(std::shared_ptr<storm::storage::geometry::Polytope<RationalNumberType>> const& polytope, PreprocessorData const& preprocessorData) {
if(polytope->isEmpty()) {
return storm::storage::geometry::Polytope<RationalNumberType>::createEmptyPolytope();
}
if(polytope->isUniversal()) {
return storm::storage::geometry::Polytope<RationalNumberType>::createUniversalPolytope();
}
uint_fast64_t numObjectives = preprocessorData.objectives.size();
std::vector<std::vector<RationalNumberType>> transformationMatrix(numObjectives, std::vector<RationalNumberType>(numObjectives, storm::utility::zero<RationalNumberType>()));
std::vector<RationalNumberType> transformationVector;
@ -299,7 +306,7 @@ namespace storm {
combinedTime.addToAccumulatedSeconds(preprocessorStopwatch->getAccumulatedSeconds());
}
if(helperStopwatch) {
statistics << "\t Value Iterations: " << std::setw(8) << *helperStopwatch << std::endl;
statistics << "\t Iterations: " << std::setw(8) << *helperStopwatch << std::endl;
combinedTime.addToAccumulatedSeconds(helperStopwatch->getAccumulatedSeconds());
}
if(postprocessorStopwatch) {
@ -317,7 +324,32 @@ namespace storm {
statistics << "Convergence precision for iterative solvers: " << settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() << std::endl;
statistics << std::endl;
statistics << "---------------------------------------------------------------------------------------------------------------------------------------" << std::endl;
statistics << " Data in CSV format " << std::endl;
statistics << "---------------------------------------------------------------------------------------------------------------------------------------" << std::endl;
statistics << "model_Header;Constants;States;Choices;Transitions;Properties;" << std::endl;
statistics << "model_data;"
<< storm::settings::modules::IOSettings().getConstantDefinitionString() << ";"
<< preprocessorData.originalModel.getNumberOfStates() << ";"
<< preprocessorData.originalModel.getNumberOfChoices() << ";"
<< preprocessorData.originalModel.getNumberOfTransitions() << ";"
<< preprocessorData.originalFormula << ";"
<< std::endl;
statistics << "result_Header;Iterations;Time Combined;Accuracy;Time Iterations;Time Computation Optimal Points" << std::endl;
statistics << "result_data;"
<< resultData.refinementSteps().size() << ";"
<< combinedTime << ";";
if(resultData.isPrecisionOfResultSet()) {
statistics << storm::utility::convertNumber<double>(resultData.getPrecisionOfResult()) << ";";
} else {
statistics << ";";
}
if(helperStopwatch) {
statistics << *helperStopwatch << ";";
} else {
statistics << ";";
}
statistics << resultData.stopWatchWeightVectorChecker << ";";
statistics << std::endl;
return statistics.str();
}

5
src/modelchecker/multiobjective/helper/SparseMultiObjectiveResultData.h
View File

@ -8,6 +8,7 @@
#include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveRefinementStep.h"
#include "src/storage/geometry/Polytope.h"
#include "src/utility/macros.h"
#include "src/utility/Stopwatch.h"
#include "src/exceptions/InvalidStateException.h"
@ -98,6 +99,9 @@ namespace storm {
return maxStepsPerformed;
}
//Keeps track of the time we spent with weight vector checking (i.e., computation of optimal points)
storm::utility::Stopwatch stopWatchWeightVectorChecker;
private:
enum class Tribool { False, True, Indeterminate };
@ -126,6 +130,7 @@ namespace storm {
//Stores whether the computation was aborted due to performing too many refinement steps
bool maxStepsPerformed = false;
};
}
}

30
src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp
View File

@ -132,10 +132,13 @@ namespace storm {
storm::storage::SparseMatrix<ValueType> deterministicBackwardTransitions = deterministicMatrix.transpose();
std::vector<ValueType> deterministicStateRewards(deterministicMatrix.getRowCount());
storm::solver::GeneralLinearEquationSolverFactory<ValueType> linearEquationSolverFactory;
//TODO check if all but one entry of weightVector is zero
//Also only compute values for objectives with weightVector != zero,
//one check can be omitted as the result can be computed back from the weighed result and the results from the remaining objectives
for(uint_fast64_t objIndex = 0; objIndex < data.objectives.size(); ++objIndex) {
// We compute an estimate for the results of the individual objectives which is obtained from the weighted result and the result of the objectives computed so far.
// Note that weightedResult = Sum_{i=1}^{n} w_i * objectiveResult_i.
std::vector<ValueType> weightedSumOfUncheckedObjectives = weightedResult;
ValueType sumOfWeightsOfUncheckedObjectives = storm::utility::vector::sum_if(weightVector, objectivesWithNoUpperTimeBound);
for(uint_fast64_t const& objIndex : storm::utility::vector::getSortedIndices(weightVector)) {
if(objectivesWithNoUpperTimeBound.get(objIndex)){
offsetsToLowerBound[objIndex] = storm::utility::zero<ValueType>();
offsetsToUpperBound[objIndex] = storm::utility::zero<ValueType>();
@ -144,13 +147,28 @@ namespace storm {
// As target states, we pick the states from which no reward is reachable.
storm::storage::BitVector targetStates = ~storm::utility::graph::performProbGreater0(deterministicBackwardTransitions, storm::storage::BitVector(deterministicMatrix.getRowCount(), true), statesWithRewards);
//TODO: we could give the solver some hint for the result (e.g., weightVector[ObjIndex] * weightedResult or (weightedResult - sum_i=0^objIndex-1 objectiveResult) * weightVector[objIndex]/ sum_i=objIndex^n weightVector[i] )
// Compute the estimate for this objective
if(!storm::utility::isZero(weightVector[objIndex])) {
objectiveResults[objIndex] = weightedSumOfUncheckedObjectives;
storm::utility::vector::scaleVectorInPlace(objectiveResults[objIndex], storm::utility::one<ValueType>() / sumOfWeightsOfUncheckedObjectives);
}
// Make sure that the objectiveResult is initialized in some way
objectiveResults[objIndex].resize(data.preprocessedModel.getNumberOfStates(), storm::utility::zero<ValueType>());
// Invoke the linear equation solver
objectiveResults[objIndex] = storm::modelchecker::helper::SparseDtmcPrctlHelper<ValueType>::computeReachabilityRewards(deterministicMatrix,
deterministicBackwardTransitions,
deterministicStateRewards,
targetStates,
false, //no qualitative checking,
linearEquationSolverFactory);
linearEquationSolverFactory,
objectiveResults[objIndex]);
// Update the estimate for the next objectives.
if(!storm::utility::isZero(weightVector[objIndex])) {
storm::utility::vector::addScaledVector(weightedSumOfUncheckedObjectives, objectiveResults[objIndex], -weightVector[objIndex]);
sumOfWeightsOfUncheckedObjectives -= weightVector[objIndex];
}
} else {
objectiveResults[objIndex] = std::vector<ValueType>(data.preprocessedModel.getNumberOfStates(), storm::utility::zero<ValueType>());
}

1
src/models/sparse/MarkovAutomaton.cpp
View File

@ -360,6 +360,7 @@ namespace storm {
this->printModelInformationHeaderToStream(out);
out << "Choices: \t" << this->getNumberOfChoices() << std::endl;
out << "Markovian St.: \t" << this->getMarkovianStates().getNumberOfSetBits() << std::endl;
out << "Max. Rate.: \t" << this->getMaximalExitRate() << std::endl;
this->printModelInformationFooterToStream(out);
}

19
src/solver/GmmxxLinearEquationSolver.cpp
View File

@ -125,6 +125,7 @@ namespace storm {
localA.reset();
this->A = &A;
gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
jacobiDecomposition.reset();
}
template<typename ValueType>
@ -213,23 +214,25 @@ namespace storm {
this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
}
// Get a Jacobi decomposition of the matrix A.
std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
// Convert the LU matrix to gmm++'s format.
std::unique_ptr<gmm::csr_matrix<ValueType>> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.first));
// Get a Jacobi decomposition of the matrix A (if not already available).
if(jacobiDecomposition == nullptr) {
std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> nativeJacobiDecomposition = A.getJacobiDecomposition();
// Convert the LU matrix to gmm++'s format.
jacobiDecomposition = std::make_unique<std::pair<gmm::csr_matrix<ValueType>, std::vector<ValueType>>>(*storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(nativeJacobiDecomposition.first)),
std::move(nativeJacobiDecomposition.second));
}
std::vector<ValueType>* currentX = &x;
std::vector<ValueType>* nextX = auxiliaryJacobiMemory.get();
std::vector<ValueType>* nextX = auxiliaryJacobiMemory.get();
// Set up additional environment variables.
uint_fast64_t iterationCount = 0;
bool converged = false;
while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
// Compute D^-1 * (b - LU * x) and store result in nextX.
gmm::mult_add(*gmmLU, gmm::scaled(*currentX, -storm::utility::one<ValueType>()), b, *nextX);
storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX);
gmm::mult_add(jacobiDecomposition->first, gmm::scaled(*currentX, -storm::utility::one<ValueType>()), b, *nextX);
storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition->second, *nextX, *nextX);
// Now check if the process already converged within our precision.
converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion());

1
src/solver/GmmxxLinearEquationSolver.h
View File

@ -138,6 +138,7 @@ namespace storm {
// Auxiliary storage for the Jacobi method.
mutable std::unique_ptr<std::vector<ValueType>> auxiliaryJacobiMemory;
mutable std::unique_ptr<std::pair<gmm::csr_matrix<ValueType>, std::vector<ValueType>>> jacobiDecomposition;
};
template<typename ValueType>

12
src/solver/StandardMinMaxLinearEquationSolver.cpp
View File

@ -107,7 +107,7 @@ namespace storm {
storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A.getRowGroupIndices(), b);
// Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration.
auto solver = linearEquationSolverFactory->create(std::move(submatrix));
solver = linearEquationSolverFactory->create(std::move(submatrix));
solver->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
Status status = Status::InProgress;
@ -165,6 +165,8 @@ namespace storm {
if (this->isTrackSchedulerSet()) {
this->scheduler = std::make_unique<storm::storage::TotalScheduler>(std::move(scheduler));
}
solver.reset();
if(status == Status::Converged || status == Status::TerminatedEarly) {
return true;
@ -202,7 +204,9 @@ namespace storm {
template<typename ValueType>
bool StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
if(solver == nullptr) {
solver = linearEquationSolverFactory->create(A);
}
bool allocatedAuxMemory = !this->hasAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations);
if (allocatedAuxMemory) {
this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations);
@ -274,7 +278,9 @@ namespace storm {
this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly);
}
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
if(solver == nullptr) {
solver = linearEquationSolverFactory->create(A);
}
for (uint64_t i = 0; i < n; ++i) {
solver->multiply(x, b, *auxiliaryRepeatedMultiplyMemory);

3
src/solver/StandardMinMaxLinearEquationSolver.h
View File

@ -68,6 +68,7 @@ namespace storm {
/// The factory used to obtain linear equation solvers.
std::unique_ptr<LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
mutable std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver;
// If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
// when the solver is destructed.
@ -129,4 +130,4 @@ namespace storm {
};
}
}
}

6
src/storage/SparseMatrix.cpp
View File

@ -695,6 +695,9 @@ namespace storm {
*writeIt = std::move(intermediateRowEntry);
++writeIt;
}
} else {
// skip the intermediate rows
writeIt = getRow(smallerRow).begin();
}
// write the larger row
for(auto& largerRowEntry : largerRowContents) {
@ -719,6 +722,9 @@ namespace storm {
*writeIt = std::move(*intermediateRowEntryIt);
--writeIt;
}
} else {
// skip the intermediate rows
writeIt = getRow(smallerRow).end() - 1;
}
// write the larger row
for(auto largerRowEntryIt = largerRowContents.rbegin(); largerRowEntryIt != largerRowContents.rend(); ++largerRowEntryIt) {

21
src/utility/vector.h
View File

@ -101,7 +101,22 @@ namespace storm {
}
/*!
* Retrieves a list of indices such that the first index points to the entry of the given vector
* with the highest value, the second index points to the entry with the second highest value, ...
* Example: v={3,8,4,5,1} yields res={1,3,2,0,4}
*/
template<typename T>
std::vector<uint_fast64_t> getSortedIndices(std::vector<T> const& v){
std::vector<uint_fast64_t> res = buildVectorForRange(0, v.size());
std::sort(res.begin(), res.end(), [&v](uint_fast64_t index1, uint_fast64_t index2) { return v[index1] > v[index2];});
return res;
}
/*!
* Selects the elements from a vector at the specified positions and writes them consecutively into another vector.
* @param vector The vector into which the selected elements are to be written.
@ -383,7 +398,7 @@ namespace storm {
}
/*!
* Adds each element of the first vector and (the corresponding element of the second vector times the given factor) and writes the result into the first vector.
* Computes x:= x + a*y, i.e., adds each element of the first vector and (the corresponding element of the second vector times the given factor) and writes the result into the first vector.
*
* @param firstOperand The first operand.
* @param secondOperand The second operand
@ -802,7 +817,7 @@ namespace storm {
bool hasNonZeroEntry(std::vector<T> const& v){
return std::any_of(v.begin(), v.end(), [](T value){return !storm::utility::isZero(value);});
}
/*!
* Output vector as string.
*

4
test/functional/modelchecker/SparseMdpMultiObjectiveModelCheckerTest.cpp
View File

@ -30,11 +30,11 @@ TEST(SparseMdpMultiObjectiveModelCheckerTest, probEqual1Objective) {
std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(storm::modelchecker::CheckTask<storm::logic::Formula, double>(*formulas[0], true));
ASSERT_TRUE(result->isExplicitQuantitativeCheckResult());
EXPECT_NEAR(7.647058824, result->asExplicitQuantitativeCheckResult<double>()[initState], storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision());
EXPECT_NEAR(7.647057, result->asExplicitQuantitativeCheckResult<double>()[initState], storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision());
result = checker.check(storm::modelchecker::CheckTask<storm::logic::Formula, double>(*formulas[1], true));
ASSERT_TRUE(result->isExplicitQuantitativeCheckResult());
EXPECT_NEAR(7.647058824, result->asExplicitQuantitativeCheckResult<double>()[initState], storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision());
EXPECT_NEAR(7.647057, result->asExplicitQuantitativeCheckResult<double>()[initState], storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision());
}

Write Preview
|||||||
100:0
Loading…
Cancel
Save