From 71e181bd938ec0a9a315ee0b9b5da35a510dfa91 Mon Sep 17 00:00:00 2001 From: TimQu Date: Tue, 1 Nov 2016 19:29:03 +0100 Subject: [PATCH] some quick and dirty optimizations (that need to be reviewed) Former-commit-id: 86fbc66c903d7305cb5fc723c6b049de8eeeec84 --- examples/multiobjective/ma/polling/polling.ma | 4 +- .../multiobjective/ma/polling/untitled.txt | 16 ++ examples/multiobjective/ma/stream/stream.ma | 21 +-- examples/multiobjective/ma/stream/stream10.ma | 50 ------ .../multiobjective/ma/stream/stream100.ma | 50 ------ .../multiobjective/ma/stream/stream1000.ma | 50 ------ examples/multiobjective/ma/stream/stream2.ma | 50 ------ examples/multiobjective/ma/stream/stream50.ma | 50 ------ .../multiobjective/ma/stream/stream500.ma | 50 ------ .../multiobjective/ma/stream/stream5000.ma | 50 ------ resources/3rdparty/carl | 1 - src/logic/Bound.h | 2 +- src/logic/OperatorFormula.h | 5 +- ...rseMaMultiObjectiveWeightVectorChecker.cpp | 54 +++++-- .../helper/SparseMultiObjectiveHelper.cpp | 148 +++++++++++++++--- .../helper/SparseMultiObjectiveHelper.h | 2 +- .../SparseMultiObjectivePostprocessor.cpp | 40 ++++- .../helper/SparseMultiObjectiveResultData.h | 5 + ...parseMultiObjectiveWeightVectorChecker.cpp | 30 +++- src/models/sparse/MarkovAutomaton.cpp | 1 + src/solver/GmmxxLinearEquationSolver.cpp | 19 ++- src/solver/GmmxxLinearEquationSolver.h | 1 + .../StandardMinMaxLinearEquationSolver.cpp | 12 +- .../StandardMinMaxLinearEquationSolver.h | 3 +- src/storage/SparseMatrix.cpp | 6 + src/utility/vector.h | 21 ++- ...parseMdpMultiObjectiveModelCheckerTest.cpp | 4 +- 27 files changed, 314 insertions(+), 431 deletions(-) create mode 100644 examples/multiobjective/ma/polling/untitled.txt delete mode 100644 examples/multiobjective/ma/stream/stream10.ma delete mode 100644 examples/multiobjective/ma/stream/stream100.ma delete mode 100644 examples/multiobjective/ma/stream/stream1000.ma delete mode 100644 examples/multiobjective/ma/stream/stream2.ma delete mode 100644 examples/multiobjective/ma/stream/stream50.ma delete mode 100644 examples/multiobjective/ma/stream/stream500.ma delete mode 100644 examples/multiobjective/ma/stream/stream5000.ma delete mode 160000 resources/3rdparty/carl diff --git a/examples/multiobjective/ma/polling/polling.ma b/examples/multiobjective/ma/polling/polling.ma index a7523d2c6..5c8a45ad8 100644 --- a/examples/multiobjective/ma/polling/polling.ma +++ b/examples/multiobjective/ma/polling/polling.ma @@ -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" diff --git a/examples/multiobjective/ma/polling/untitled.txt b/examples/multiobjective/ma/polling/untitled.txt new file mode 100644 index 000000000..2be3cb7df --- /dev/null +++ b/examples/multiobjective/ma/polling/untitled.txt @@ -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 + + diff --git a/examples/multiobjective/ma/stream/stream.ma b/examples/multiobjective/ma/stream/stream.ma index 820934a8c..9a657dd04 100644 --- a/examples/multiobjective/ma/stream/stream.ma +++ b/examples/multiobjective/ma/stream/stream.ma @@ -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 1 : (s'=1); + [buffer] s=0 & n 1 : (s'=1); + [buffer] s=0 & n 0.99: (s'=1) + 0.01 : (s'=2) & (k'=k+1); [start] s=0 & k 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 diff --git a/examples/multiobjective/ma/stream/stream10.ma b/examples/multiobjective/ma/stream/stream10.ma deleted file mode 100644 index c7fc431f9..000000000 --- a/examples/multiobjective/ma/stream/stream10.ma +++ /dev/null @@ -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 1 : (s'=1); - [start] s=0 & k 1 : (s'=2) & (k'=k+1); - - <> s=1 -> inRate : (n'=n+1) & (s'=0); - - <> s=2 & n inRate : (n'=n+1) + processingRate : (k'=k+1); - <> s=2 & n inRate : (n'=n+1) + processingRate : (s'=0); - <> s=2 & n=N & k 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 diff --git a/examples/multiobjective/ma/stream/stream100.ma b/examples/multiobjective/ma/stream/stream100.ma deleted file mode 100644 index 8eee818ba..000000000 --- a/examples/multiobjective/ma/stream/stream100.ma +++ /dev/null @@ -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 1 : (s'=1); - [start] s=0 & k 1 : (s'=2) & (k'=k+1); - - <> s=1 -> inRate : (n'=n+1) & (s'=0); - - <> s=2 & n inRate : (n'=n+1) + processingRate : (k'=k+1); - <> s=2 & n inRate : (n'=n+1) + processingRate : (s'=0); - <> s=2 & n=N & k 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 diff --git a/examples/multiobjective/ma/stream/stream1000.ma b/examples/multiobjective/ma/stream/stream1000.ma deleted file mode 100644 index 1860928f7..000000000 --- a/examples/multiobjective/ma/stream/stream1000.ma +++ /dev/null @@ -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 1 : (s'=1); - [start] s=0 & k 1 : (s'=2) & (k'=k+1); - - <> s=1 -> inRate : (n'=n+1) & (s'=0); - - <> s=2 & n inRate : (n'=n+1) + processingRate : (k'=k+1); - <> s=2 & n inRate : (n'=n+1) + processingRate : (s'=0); - <> s=2 & n=N & k 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 diff --git a/examples/multiobjective/ma/stream/stream2.ma b/examples/multiobjective/ma/stream/stream2.ma deleted file mode 100644 index 56378a490..000000000 --- a/examples/multiobjective/ma/stream/stream2.ma +++ /dev/null @@ -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 1 : (s'=1); - [start] s=0 & k 1 : (s'=2) & (k'=k+1); - - <> s=1 -> inRate : (n'=n+1) & (s'=0); - - <> s=2 & n inRate : (n'=n+1) + processingRate : (k'=k+1); - <> s=2 & n inRate : (n'=n+1) + processingRate : (s'=0); - <> s=2 & n=N & k 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 diff --git a/examples/multiobjective/ma/stream/stream50.ma b/examples/multiobjective/ma/stream/stream50.ma deleted file mode 100644 index 190122cde..000000000 --- a/examples/multiobjective/ma/stream/stream50.ma +++ /dev/null @@ -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 1 : (s'=1); - [start] s=0 & k 1 : (s'=2) & (k'=k+1); - - <> s=1 -> inRate : (n'=n+1) & (s'=0); - - <> s=2 & n inRate : (n'=n+1) + processingRate : (k'=k+1); - <> s=2 & n inRate : (n'=n+1) + processingRate : (s'=0); - <> s=2 & n=N & k 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 diff --git a/examples/multiobjective/ma/stream/stream500.ma b/examples/multiobjective/ma/stream/stream500.ma deleted file mode 100644 index 90933bd63..000000000 --- a/examples/multiobjective/ma/stream/stream500.ma +++ /dev/null @@ -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 1 : (s'=1); - [start] s=0 & k 1 : (s'=2) & (k'=k+1); - - <> s=1 -> inRate : (n'=n+1) & (s'=0); - - <> s=2 & n inRate : (n'=n+1) + processingRate : (k'=k+1); - <> s=2 & n inRate : (n'=n+1) + processingRate : (s'=0); - <> s=2 & n=N & k 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 diff --git a/examples/multiobjective/ma/stream/stream5000.ma b/examples/multiobjective/ma/stream/stream5000.ma deleted file mode 100644 index 32dea35a6..000000000 --- a/examples/multiobjective/ma/stream/stream5000.ma +++ /dev/null @@ -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 1 : (s'=1); - [start] s=0 & k 1 : (s'=2) & (k'=k+1); - - <> s=1 -> inRate : (n'=n+1) & (s'=0); - - <> s=2 & n inRate : (n'=n+1) + processingRate : (k'=k+1); - <> s=2 & n inRate : (n'=n+1) + processingRate : (s'=0); - <> s=2 & n=N & k 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 diff --git a/resources/3rdparty/carl b/resources/3rdparty/carl deleted file mode 160000 index d67f98622..000000000 --- a/resources/3rdparty/carl +++ /dev/null @@ -1 +0,0 @@ -Subproject commit d67f986226cf846ba6366cdbe2abc21dff375542 diff --git a/src/logic/Bound.h b/src/logic/Bound.h index e56a2c04f..fbfc2fcf8 100644 --- a/src/logic/Bound.h +++ b/src/logic/Bound.h @@ -26,7 +26,7 @@ namespace storm { template std::ostream& operator<<(std::ostream& out, Bound const& bound) { - out << bound.comparisonType << bound.threshold; + out << bound.comparisonType << storm::utility::convertNumber(bound.threshold); return out; } } diff --git a/src/logic/OperatorFormula.h b/src/logic/OperatorFormula.h index a1247d431..fd5495da6 100644 --- a/src/logic/OperatorFormula.h +++ b/src/logic/OperatorFormula.h @@ -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_ */ \ No newline at end of file +#endif /* STORM_LOGIC_OPERATORFORMULA_H_ */ diff --git a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp index a891f1b1c..688ed3709 100644 --- a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp @@ -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 void SparseMaMultiObjectiveWeightVectorChecker::performPSStep(SubModel& PS, SubModel const& MS, MinMaxSolverData& minMax, LinEqSolverData& linEq, std::vector& 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 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 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* gmmSolver = dynamic_cast*>(linEq.solver.get()); + if(gmmSolver!=nullptr) { + STORM_LOG_INFO("Switching to Jacobi method"); + gmmSolver->getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi); + gmmSolver->getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings::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>::digitize(SparseMaMultiObjectiveWeightVectorChecker>::SubModel& subModel, double const& digitizationConstant) const; template void SparseMaMultiObjectiveWeightVectorChecker>::digitizeTimeBounds(SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& upperTimeBounds, double const& digitizationConstant); #ifdef STORM_HAVE_CARL - template class SparseMaMultiObjectiveWeightVectorChecker>; - template storm::RationalNumber SparseMaMultiObjectiveWeightVectorChecker>::getDigitizationConstant() const; - template void SparseMaMultiObjectiveWeightVectorChecker>::digitize(SparseMaMultiObjectiveWeightVectorChecker>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const; - template void SparseMaMultiObjectiveWeightVectorChecker>::digitizeTimeBounds(SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant); +// template class SparseMaMultiObjectiveWeightVectorChecker>; + // template storm::RationalNumber SparseMaMultiObjectiveWeightVectorChecker>::getDigitizationConstant() const; + // template void SparseMaMultiObjectiveWeightVectorChecker>::digitize(SparseMaMultiObjectiveWeightVectorChecker>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const; +// template void SparseMaMultiObjectiveWeightVectorChecker>::digitizeTimeBounds(SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant); #endif } diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp index 627db4eba..fb63c037b 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp @@ -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::createUniversalPolytope(); resultData.underApproximation() = storm::storage::geometry::Polytope::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(storm::settings::getModule().getPrecision()); + gap /= (storm::utility::one() + storm::utility::one()); + gap /= storm::utility::sqrt(static_cast(preprocessorData.objectives.size())); + weightVectorChecker->setMaximumLowerUpperBoundGap(gap); + if(!checkIfPreprocessingWasConclusive(preprocessorData)) { switch(preprocessorData.queryType) { case PreprocessorData::QueryType::Achievability: @@ -57,8 +70,6 @@ namespace storm { template void SparseMultiObjectiveHelper::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(0.1)); // Get a point that represents the thresholds Point thresholds; thresholds.reserve(preprocessorData.objectives.size()); @@ -84,8 +95,6 @@ namespace storm { template void SparseMultiObjectiveHelper::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(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(storm::settings::getModule().getPrecision()); - gap /= (storm::utility::one() + storm::utility::one()); - weightVectorChecker->setMaximumLowerUpperBoundGap(gap); - while(true) { std::pair 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()); + std::cout << "Current precision of the numerical result is " << resultData.template getPrecisionOfResult() << std::endl; } if(targetPrecisionReached(resultData) || maxStepsPerformed(resultData)) { resultData.setOptimumIsAchievable(checkIfThresholdsAreSatisfied(resultData.underApproximation(), thresholds, strictThresholds)); @@ -161,19 +165,63 @@ namespace storm { template void SparseMultiObjectiveHelper::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(storm::settings::getModule().getPrecision()); - gap /= (storm::utility::one() + storm::utility::one()); - gap /= storm::utility::sqrt(static_cast(preprocessorData.objectives.size())); - weightVectorChecker->setMaximumLowerUpperBoundGap(gap); - //First consider the objectives individually for(uint_fast64_t objIndex = 0; objIndex()); direction[objIndex] = storm::utility::one(); 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(resultData.refinementSteps().back().getLowerBoundPoint()[subFIndex]) * storm::utility::convertNumber(preprocessorData.objectives[subFIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber(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::logic::ComparisonType::GreaterEqual, thresholdValue)); + copySubF->asOperatorFormula().setOperatorInformation(opInfo); + } else { + storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound(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(resultData.refinementSteps().back().getLowerBoundPoint()[subFIndex]) * storm::utility::convertNumber(preprocessorData.objectives[subFIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber(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::logic::ComparisonType::GreaterEqual, thresholdValue)); + copySubF->asOperatorFormula().setOperatorInformation(opInfo); + } else { + storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound(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() << 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(resultData.refinementSteps().back().getLowerBoundPoint()[subFIndex]) * storm::utility::convertNumber(preprocessorData.objectives[subFIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber(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::logic::ComparisonType::GreaterEqual, thresholdValue)); + copySubF->asOperatorFormula().setOperatorInformation(opInfo); + } else { + storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound(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(resultData.refinementSteps().back().getLowerBoundPoint()[subFIndex]) * storm::utility::convertNumber(preprocessorData.objectives[subFIndex].toOriginalValueTransformationFactor) + storm::utility::convertNumber(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::logic::ComparisonType::GreaterEqual, thresholdValue)); + copySubF->asOperatorFormula().setOperatorInformation(opInfo); + } else { + storm::logic::OperatorInformation opInfo(boost::none, storm::logic::Bound(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()); if(targetPrecisionReached(resultData) || maxStepsPerformed(resultData)) { break; @@ -264,7 +371,11 @@ namespace storm { break; } } + + result.stopWatchWeightVectorChecker.unpause(); weightVectorChecker->check(storm::utility::vector::convertNumericVector(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 diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h index c47aefaa8..0afb224b7 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.h @@ -52,7 +52,7 @@ namespace storm { static WeightVector findSeparatingVector(Point const& pointToBeSeparated, std::shared_ptr> 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); diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectivePostprocessor.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectivePostprocessor.cpp index 0bfea1603..b7307ccb8 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectivePostprocessor.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectivePostprocessor.cpp @@ -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 std::shared_ptr> SparseMultiObjectivePostprocessor::transformToOriginalValues(std::shared_ptr> const& polytope, PreprocessorData const& preprocessorData) { - + if(polytope->isEmpty()) { + return storm::storage::geometry::Polytope::createEmptyPolytope(); + } + if(polytope->isUniversal()) { + return storm::storage::geometry::Polytope::createUniversalPolytope(); + } uint_fast64_t numObjectives = preprocessorData.objectives.size(); std::vector> transformationMatrix(numObjectives, std::vector(numObjectives, storm::utility::zero())); std::vector 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().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(resultData.getPrecisionOfResult()) << ";"; + } else { + statistics << ";"; + } + if(helperStopwatch) { + statistics << *helperStopwatch << ";"; + } else { + statistics << ";"; + } + statistics << resultData.stopWatchWeightVectorChecker << ";"; + statistics << std::endl; return statistics.str(); } diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveResultData.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveResultData.h index b9ebba998..75f025fcd 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveResultData.h +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveResultData.h @@ -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; + }; } } diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp index b6b98ab2c..3c748b16e 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp @@ -132,10 +132,13 @@ namespace storm { storm::storage::SparseMatrix deterministicBackwardTransitions = deterministicMatrix.transpose(); std::vector deterministicStateRewards(deterministicMatrix.getRowCount()); storm::solver::GeneralLinearEquationSolverFactory 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 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(); offsetsToUpperBound[objIndex] = storm::utility::zero(); @@ -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() / sumOfWeightsOfUncheckedObjectives); + } + + // Make sure that the objectiveResult is initialized in some way + objectiveResults[objIndex].resize(data.preprocessedModel.getNumberOfStates(), storm::utility::zero()); + + // Invoke the linear equation solver objectiveResults[objIndex] = storm::modelchecker::helper::SparseDtmcPrctlHelper::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(data.preprocessedModel.getNumberOfStates(), storm::utility::zero()); } diff --git a/src/models/sparse/MarkovAutomaton.cpp b/src/models/sparse/MarkovAutomaton.cpp index 887d3a13b..770d8f028 100644 --- a/src/models/sparse/MarkovAutomaton.cpp +++ b/src/models/sparse/MarkovAutomaton.cpp @@ -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); } diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index a1f57beec..fd15e7e27 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -125,6 +125,7 @@ namespace storm { localA.reset(); this->A = &A; gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + jacobiDecomposition.reset(); } template @@ -213,23 +214,25 @@ namespace storm { this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations); } - // Get a Jacobi decomposition of the matrix A. - std::pair, std::vector> jacobiDecomposition = A.getJacobiDecomposition(); - - // Convert the LU matrix to gmm++'s format. - std::unique_ptr> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.first)); + // Get a Jacobi decomposition of the matrix A (if not already available). + if(jacobiDecomposition == nullptr) { + std::pair, std::vector> nativeJacobiDecomposition = A.getJacobiDecomposition(); + // Convert the LU matrix to gmm++'s format. + jacobiDecomposition = std::make_unique, std::vector>>(*storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(nativeJacobiDecomposition.first)), + std::move(nativeJacobiDecomposition.second)); + } std::vector* currentX = &x; - std::vector* nextX = auxiliaryJacobiMemory.get(); + std::vector* 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()), b, *nextX); - storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX); + gmm::mult_add(jacobiDecomposition->first, gmm::scaled(*currentX, -storm::utility::one()), 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()); diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 7d4c09ad1..825c35a3d 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -138,6 +138,7 @@ namespace storm { // Auxiliary storage for the Jacobi method. mutable std::unique_ptr> auxiliaryJacobiMemory; + mutable std::unique_ptr, std::vector>> jacobiDecomposition; }; template diff --git a/src/solver/StandardMinMaxLinearEquationSolver.cpp b/src/solver/StandardMinMaxLinearEquationSolver.cpp index 6e0313cf9..714356b91 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/solver/StandardMinMaxLinearEquationSolver.cpp @@ -107,7 +107,7 @@ namespace storm { storm::utility::vector::selectVectorValues(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(std::move(scheduler)); } + + solver.reset(); if(status == Status::Converged || status == Status::TerminatedEarly) { return true; @@ -202,7 +204,9 @@ namespace storm { template bool StandardMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - std::unique_ptr> 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> solver = linearEquationSolverFactory->create(A); + if(solver == nullptr) { + solver = linearEquationSolverFactory->create(A); + } for (uint64_t i = 0; i < n; ++i) { solver->multiply(x, b, *auxiliaryRepeatedMultiplyMemory); diff --git a/src/solver/StandardMinMaxLinearEquationSolver.h b/src/solver/StandardMinMaxLinearEquationSolver.h index d6d36b272..86a960694 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/solver/StandardMinMaxLinearEquationSolver.h @@ -68,6 +68,7 @@ namespace storm { /// The factory used to obtain linear equation solvers. std::unique_ptr> linearEquationSolverFactory; + mutable std::unique_ptr> 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 { }; } -} \ No newline at end of file +} diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index eb0738580..dd996bed4 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -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) { diff --git a/src/utility/vector.h b/src/utility/vector.h index 05d61f781..f6680cf01 100644 --- a/src/utility/vector.h +++ b/src/utility/vector.h @@ -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 + std::vector getSortedIndices(std::vector const& v){ + std::vector 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 const& v){ return std::any_of(v.begin(), v.end(), [](T value){return !storm::utility::isZero(value);}); } - + /*! * Output vector as string. * diff --git a/test/functional/modelchecker/SparseMdpMultiObjectiveModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpMultiObjectiveModelCheckerTest.cpp index ab17dc3a5..f99939612 100644 --- a/test/functional/modelchecker/SparseMdpMultiObjectiveModelCheckerTest.cpp +++ b/test/functional/modelchecker/SparseMdpMultiObjectiveModelCheckerTest.cpp @@ -30,11 +30,11 @@ TEST(SparseMdpMultiObjectiveModelCheckerTest, probEqual1Objective) { std::unique_ptr result = checker.check(storm::modelchecker::CheckTask(*formulas[0], true)); ASSERT_TRUE(result->isExplicitQuantitativeCheckResult()); - EXPECT_NEAR(7.647058824, result->asExplicitQuantitativeCheckResult()[initState], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(7.647057, result->asExplicitQuantitativeCheckResult()[initState], storm::settings::getModule().getPrecision()); result = checker.check(storm::modelchecker::CheckTask(*formulas[1], true)); ASSERT_TRUE(result->isExplicitQuantitativeCheckResult()); - EXPECT_NEAR(7.647058824, result->asExplicitQuantitativeCheckResult()[initState], storm::settings::getModule().getPrecision()); + EXPECT_NEAR(7.647057, result->asExplicitQuantitativeCheckResult()[initState], storm::settings::getModule().getPrecision()); }