diff --git a/src/adapters/CarlAdapter.h b/src/adapters/CarlAdapter.h index d31802717..3100b13f8 100644 --- a/src/adapters/CarlAdapter.h +++ b/src/adapters/CarlAdapter.h @@ -36,7 +36,8 @@ namespace carl { namespace storm { typedef carl::Variable Variable; - typedef carl::MultivariatePolynomial RawPolynomial; + typedef cln::cl_RA Coefficient; + typedef carl::MultivariatePolynomial RawPolynomial; typedef carl::FactorizedPolynomial Polynomial; typedef carl::CompareRelation CompareRelation; typedef carl::RationalFunction RationalFunction; diff --git a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp index 023607f44..3a4e51999 100644 --- a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp +++ b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.cpp @@ -174,7 +174,7 @@ namespace storm { template - SparseDtmcRegionModelChecker::SparseDtmcRegionModelChecker(storm::models::sparse::Dtmc const& model) : model(model), eliminationModelChecker(model), smtSolver(nullptr), probabilityOperatorFormula(nullptr) { + SparseDtmcRegionModelChecker::SparseDtmcRegionModelChecker(storm::models::sparse::Dtmc const& model) : model(model), eliminationModelChecker(model), smtSolver(nullptr), probabilityOperatorFormula(nullptr), sampleDtmc(nullptr), approxMdp(nullptr) { // Intentionally left empty. } @@ -255,6 +255,8 @@ namespace storm { eliminateStatesConstSucc(this->subsystem, this->flexibleTransitions, this->flexibleBackwardTransitions, this->oneStepProbabilities, this->hasOnlyLinearFunctions, this->initialState); STORM_LOG_DEBUG("Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl); std::cout << "Eliminated " << subsystem.size() - subsystem.getNumberOfSetBits() << " of " << subsystem.size() << " states that had constant outgoing transitions." << std::endl; + + //eliminate the remaining states to get the reachability probability function this->sparseTransitions = flexibleTransitions.getSparseMatrix(); this->sparseBackwardTransitions = this->sparseTransitions.transpose(); @@ -270,6 +272,7 @@ namespace storm { std::chrono::high_resolution_clock::time_point timeInitialStateEliminationEnd = std::chrono::high_resolution_clock::now(); + initializeSampleDtmcAndApproxMdp(this->sampleDtmc,this->sampleDtmcMapping,this->approxMdp,this->approxMdpMapping,this->approxMdpSubstitutions,this->subsystem, this->sparseTransitions,this->oneStepProbabilities, this->initialState); initializeSMTSolver(this->smtSolver, this->reachProbFunction,*this->probabilityOperatorFormula); //some information for statistics... @@ -282,6 +285,9 @@ namespace storm { this->numOfRegionsSolvedThroughApproximation=0; this->numOfRegionsSolvedThroughSubsystemSmt=0; this->numOfRegionsSolvedThroughFullSmt=0; + this->numOfRegionsExistsBoth=0; + this->numOfRegionsAllSat=0; + this->numOfRegionsAllViolated=0; this->timeCheckRegion=std::chrono::high_resolution_clock::duration::zero(); this->timeSampling=std::chrono::high_resolution_clock::duration::zero(); this->timeApproximation=std::chrono::high_resolution_clock::duration::zero(); @@ -325,6 +331,203 @@ namespace storm { subsys.set(initialState, true); } + template + void SparseDtmcRegionModelChecker::initializeSampleDtmcAndApproxMdp( + std::shared_ptr>& sampleDtmc, + std::vector&>> & sampleDtmcMapping, + std::shared_ptr> & approxMdp, + std::vector&, std::size_t>> & approxMdpMapping, + std::vector> & approxMdpSubstitutions, + storm::storage::BitVector const& subsys, + storm::storage::SparseMatrix const& transitions, + std::vector const& oneStepProbs, + storm::storage::sparse::state_type const& initState) { + + //the matrix "transitions" might have empty rows where states have been eliminated. + //The DTMC/ MDP matrices should not have such rows. We therefore leave them out, but we have to change the indices of the states accordingly. + std::vector newStateIndexMap(transitions.getRowCount(), transitions.getRowCount()); //initialize with some illegal index to easily check if a transition leads to an unselected state + storm::storage::sparse::state_type newStateIndex=0; + for(auto const& state : subsys){ + newStateIndexMap[state]=newStateIndex; + ++newStateIndex; + } + //We need to add a target state to which the oneStepProbabilities will lead as well as a sink state to which the "missing" probability will lead + storm::storage::sparse::state_type numStates=newStateIndex+2; + storm::storage::sparse::state_type targetState=numStates-2; + storm::storage::sparse::state_type sinkState= numStates-1; + + //We can now fill in the dummy data of the transition matrices for dtmc and mdp + std::vector oneStepToSink(oneStepProbs.size(), storm::utility::zero()); // -2 because not needed for target and sink state + storm::storage::SparseMatrixBuilder dtmcMatrixBuilder(numStates, numStates); + storm::storage::SparseMatrixBuilder mdpMatrixBuilder(0, numStates, 0, true, true, numStates); + std::vector > mdpRowSubstitutions; + storm::storage::sparse::state_type currentMdpRow=0; + //go through rows: + for(storm::storage::sparse::state_type oldStateIndex : subsys){ + ParametricType valueToSinkState=storm::utility::regions::getNewFunction(storm::utility::one()); + // the dtmc and the mdp rows need to sum up to one, therefore the first entry that we add has value one. + ConstantType dummyEntry=storm::utility::one(); + // store the columns and values that we have added because we need the same information for the next rows in the mdp + std::vector> rowInformation; + mdpMatrixBuilder.newRowGroup(currentMdpRow); + std::set occurringParameters; + //go through columns: + for(auto const& entry: transitions.getRow(oldStateIndex)){ + valueToSinkState-=entry.getValue(); + storm::utility::regions::gatherOccurringVariables(entry.getValue(), occurringParameters); + dtmcMatrixBuilder.addNextValue(newStateIndexMap[oldStateIndex],newStateIndexMap[entry.getColumn()],dummyEntry); + mdpMatrixBuilder.addNextValue(currentMdpRow, newStateIndexMap[entry.getColumn()], dummyEntry); + rowInformation.emplace_back(std::make_pair(newStateIndexMap[entry.getColumn()], dummyEntry)); + dummyEntry=storm::utility::zero(); + STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]!=transitions.getRowCount(), storm::exceptions::UnexpectedException, "There is a transition to a state that should have been eliminated."); + } + //transition to target state + if(!this->parametricTypeComparator.isZero(oneStepProbs[oldStateIndex])){ + valueToSinkState-=oneStepProbs[oldStateIndex]; + storm::utility::regions::gatherOccurringVariables(oneStepProbs[oldStateIndex], occurringParameters); + + dtmcMatrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], targetState, dummyEntry); + mdpMatrixBuilder.addNextValue(currentMdpRow, targetState, dummyEntry); + rowInformation.emplace_back(std::make_pair(targetState, dummyEntry)); + dummyEntry=storm::utility::zero(); + } + //transition to sink state + if(!this->parametricTypeComparator.isZero(valueToSinkState)){ + dtmcMatrixBuilder.addNextValue(newStateIndexMap[oldStateIndex], sinkState, dummyEntry); + mdpMatrixBuilder.addNextValue(currentMdpRow, sinkState, dummyEntry); + rowInformation.emplace_back(std::make_pair(sinkState, dummyEntry)); + oneStepToSink[oldStateIndex]=valueToSinkState; + } + // Find the different combinations of occuring variables and lower/upper bounds (mathematically speaking we want to have the set 'occuringVariables times {u,l}' ) + std::size_t numOfSubstitutions=std::pow(2,occurringParameters.size()); //1 substitution = 1 combination of lower and upper bounds + for(uint_fast64_t substitutionId=0; substitutionId>parameterIndex)%2==0 ){ + mdpRowSubstitutions.back().insert(std::make_pair(parameter, TypeOfBound::LOWER)); + } + else{ + mdpRowSubstitutions.back().insert(std::make_pair(parameter, TypeOfBound::UPPER)); + } + ++parameterIndex; + } + } + //We already added one row in the MDP. Now add the remaining 2^(occuringParameters.size())-1 rows + for(++currentMdpRow; currentMdpRow()); + dtmcMatrixBuilder.addNextValue(sinkState, sinkState, storm::utility::one()); + mdpMatrixBuilder.newRowGroup(currentMdpRow); + mdpMatrixBuilder.addNextValue(currentMdpRow, targetState, storm::utility::one()); + ++currentMdpRow; + mdpMatrixBuilder.newRowGroup(currentMdpRow); + mdpMatrixBuilder.addNextValue(currentMdpRow, sinkState, storm::utility::one()); + ++currentMdpRow; + //The DTMC labeling + storm::models::sparse::StateLabeling dtmcStateLabeling(numStates); + storm::storage::BitVector initLabel(numStates, false); + initLabel.set(newStateIndexMap[this->initialState], true); + dtmcStateLabeling.addLabel("init", std::move(initLabel)); + storm::storage::BitVector targetLabel(numStates, false); + targetLabel.set(numStates-2, true); + dtmcStateLabeling.addLabel("target", std::move(targetLabel)); + storm::storage::BitVector sinkLabel(numStates, false); + sinkLabel.set(numStates-1, true); + dtmcStateLabeling.addLabel("sink", std::move(sinkLabel)); + // The MDP labeling (is the same) + storm::models::sparse::StateLabeling mdpStateLabeling(dtmcStateLabeling); + // other ingredients + boost::optional> dtmcNoStateRewards, mdpNoStateRewards; + boost::optional> dtmcNoTransitionRewards, mdpNoTransitionRewards; + boost::optional>> dtmcNoChoiceLabeling, mdpNoChoiceLabeling; + // the final dtmc and mdp + sampleDtmc=std::make_shared>(dtmcMatrixBuilder.build(), std::move(dtmcStateLabeling), std::move(dtmcNoStateRewards), std::move(dtmcNoTransitionRewards), std::move(dtmcNoChoiceLabeling)); + approxMdp=std::make_shared>(mdpMatrixBuilder.build(), std::move(mdpStateLabeling), std::move(mdpNoStateRewards), std::move(mdpNoTransitionRewards), std::move(mdpNoChoiceLabeling)); + + //build mapping vectors and the substitution vector. + sampleDtmcMapping=std::vector&>>(); + sampleDtmcMapping.reserve(sampleDtmc->getTransitionMatrix().getEntryCount()); + approxMdpMapping=std::vector&, std::size_t>>(); + approxMdpMapping.reserve(approxMdp->getTransitionMatrix().getEntryCount()); + approxMdpSubstitutions = std::vector >(); + + currentMdpRow=0; + //go through rows: + for(storm::storage::sparse::state_type oldStateIndex : subsys){ + //Get the index of the substitution for the current mdp row. (add it to approxMdpSubstitutions if we see this substitution the first time) + std::size_t substitutionIndex = std::find(approxMdpSubstitutions.begin(),approxMdpSubstitutions.end(), mdpRowSubstitutions[currentMdpRow]) - approxMdpSubstitutions.begin(); + if(substitutionIndex==approxMdpSubstitutions.size()){ + approxMdpSubstitutions.push_back(mdpRowSubstitutions[currentMdpRow]); + } + typename storm::storage::SparseMatrix::iterator dtmcEntryIt=sampleDtmc->getTransitionMatrix().getRow(newStateIndexMap[oldStateIndex]).begin(); + typename storm::storage::SparseMatrix::iterator mdpEntryIt=approxMdp->getTransitionMatrix().getRow(currentMdpRow).begin(); + //go through columns to fill the dtmc and the first mdp row (of this group): + for(auto const& entry: transitions.getRow(oldStateIndex)){ + STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]==dtmcEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entries of the DTMC matrix and the pDtmc Transitionmatrix do not match"); + STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entries of the MDP matrix and the pDtmc Transitionmatrix do not match (" << newStateIndexMap[entry.getColumn()] << " vs " << mdpEntryIt->getColumn() << ")"); + sampleDtmcMapping.emplace_back(entry.getValue(),*dtmcEntryIt); + approxMdpMapping.emplace_back(entry.getValue(),*mdpEntryIt, substitutionIndex); + ++dtmcEntryIt; + ++mdpEntryIt; + } + //transition to target state + if(!this->parametricTypeComparator.isZero(oneStepProbs[oldStateIndex])){ + STORM_LOG_THROW(targetState==dtmcEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the DTMC matrix should match the target state but it does not."); + STORM_LOG_THROW(targetState==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the MDP matrix should match the target state but it does not."); + sampleDtmcMapping.emplace_back(oneStepProbs[oldStateIndex],*dtmcEntryIt); + approxMdpMapping.emplace_back(oneStepProbs[oldStateIndex], *mdpEntryIt, substitutionIndex); + ++dtmcEntryIt; + ++mdpEntryIt; + } + //transition to sink state + if(!this->parametricTypeComparator.isZero(oneStepToSink[oldStateIndex])){ + STORM_LOG_THROW(sinkState==dtmcEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the DTMC matrix should match the sink state but it does not."); + STORM_LOG_THROW(sinkState==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the MDP matrix should match the sink state but it does not."); + sampleDtmcMapping.emplace_back(oneStepToSink[oldStateIndex],*dtmcEntryIt); + approxMdpMapping.emplace_back(oneStepToSink[oldStateIndex], *mdpEntryIt, substitutionIndex); + ++dtmcEntryIt; + ++mdpEntryIt; + } + //fill in the remaining rows of the mdp + storm::storage::sparse::state_type startOfNextRowGroup = currentMdpRow + approxMdp->getTransitionMatrix().getRowGroupSize(newStateIndexMap[oldStateIndex]); + for(++currentMdpRow;currentMdpRowgetTransitionMatrix().getRow(currentMdpRow).begin(); + for(auto const& entry: transitions.getRow(oldStateIndex)){ + STORM_LOG_THROW(newStateIndexMap[entry.getColumn()]==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entries of the MDP matrix and the pDtmc Transitionmatrix do not match"); + approxMdpMapping.emplace_back(entry.getValue(),*mdpEntryIt, substitutionIndex); + ++mdpEntryIt; + } + //transition to target state + if(!this->parametricTypeComparator.isZero(oneStepProbs[oldStateIndex])){ + STORM_LOG_THROW(targetState==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the MDP matrix should match the target state but it does not."); + approxMdpMapping.emplace_back(oneStepProbs[oldStateIndex], *mdpEntryIt, substitutionIndex); + ++mdpEntryIt; + } + //transition to sink state + if(!this->parametricTypeComparator.isZero(oneStepToSink[oldStateIndex])){ + STORM_LOG_THROW(sinkState==mdpEntryIt->getColumn(), storm::exceptions::UnexpectedException, "The entry of the MDP matrix should match the sink state but it does not."); + approxMdpMapping.emplace_back(oneStepToSink[oldStateIndex], *mdpEntryIt, substitutionIndex); + ++mdpEntryIt; + } + } + } + } + template ParametricType SparseDtmcRegionModelChecker::computeReachProbFunction( @@ -444,7 +647,7 @@ namespace storm { STORM_LOG_THROW(this->probabilityOperatorFormula!=nullptr, storm::exceptions::InvalidStateException, "Tried to analyze a region although no property has been specified" ); STORM_LOG_DEBUG("Analyzing the region " << region.toString()); - std::cout << "Analyzing the region " << region.toString() << std::endl; + // std::cout << "Analyzing the region " << region.toString() << std::endl; //switches for the different steps. bool doApproximation=this->hasOnlyLinearFunctions; // this approach is only correct if the model has only linear functions @@ -457,13 +660,13 @@ namespace storm { std::vector upperBounds; if(doApproximation){ STORM_LOG_DEBUG("Checking approximative probabilities..."); - std::cout << " Checking approximative probabilities..." << std::endl; + // std::cout << " Checking approximative probabilities..." << std::endl; if(checkApproximativeProbabilities(region, lowerBounds, upperBounds)){ ++this->numOfRegionsSolvedThroughApproximation; STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through approximation."); doSampling=false; doSubsystemSmt=false; - doFullSmt=true; //TODO set this back to false.. this is just for testing + doFullSmt=false; } } std::chrono::high_resolution_clock::time_point timeApproximationEnd = std::chrono::high_resolution_clock::now(); @@ -471,13 +674,13 @@ namespace storm { std::chrono::high_resolution_clock::time_point timeSamplingStart = std::chrono::high_resolution_clock::now(); if(doSampling){ STORM_LOG_DEBUG("Checking sample points..."); - std::cout << " Checking sample points..." << std::endl; + // std::cout << " Checking sample points..." << std::endl; if(checkSamplePoints(region)){ ++this->numOfRegionsSolvedThroughSampling; STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through sampling."); doApproximation=false; doSubsystemSmt=false; - doFullSmt=true;//TODO set this back to false.. this is just for testing + doFullSmt=false; } } std::chrono::high_resolution_clock::time_point timeSamplingEnd = std::chrono::high_resolution_clock::now(); @@ -491,7 +694,7 @@ namespace storm { std::chrono::high_resolution_clock::time_point timeFullSmtStart = std::chrono::high_resolution_clock::now(); if(doFullSmt){ STORM_LOG_DEBUG("Checking with Smt Solving..."); - std::cout << " Checking with Smt Solving..." << std::endl; + // std::cout << " Checking with Smt Solving..." << std::endl; if(checkFullSmt(region)){ ++this->numOfRegionsSolvedThroughFullSmt; STORM_LOG_DEBUG("Result '" << region.checkResultToString() <<"' obtained through Smt Solving."); @@ -510,14 +713,38 @@ namespace storm { this->timeApproximation += timeApproximationEnd - timeApproximationStart; this->timeSubsystemSmt += timeSubsystemSmtEnd - timeSubsystemSmtStart; this->timeFullSmt += timeFullSmtEnd - timeFullSmtStart; + switch(region.getCheckResult()){ + case RegionCheckResult::EXISTSBOTH: + ++this->numOfRegionsExistsBoth; + break; + case RegionCheckResult::ALLSAT: + ++this->numOfRegionsAllSat; + break; + case RegionCheckResult::ALLVIOLATED: + ++this->numOfRegionsAllViolated; + break; + default: + break; + } } template void SparseDtmcRegionModelChecker::checkRegions(std::vector& regions) { + STORM_LOG_DEBUG("Checking " << regions.size() << "regions."); + std::cout << "Checking " << regions.size() << "regions. Progress: "; + std::cout.flush(); + + uint_fast64_t progress=0; + uint_fast64_t checkedRegions=0; for(auto& region : regions){ this->checkRegion(region); + if((checkedRegions++)*10/regions.size()==progress){ + std::cout << progress++; + std::cout.flush(); + } } + std::cout << std::endl; } template @@ -525,7 +752,7 @@ namespace storm { auto samplingPoints = region.getVerticesOfRegion(region.getVariables()); //test the 4 corner points for (auto const& point : samplingPoints){ - if(checkPoint(region, point, true)){ + if(checkPoint(region, point, false)){ return true; } } @@ -540,7 +767,26 @@ namespace storm { valueInBoundOfFormula = this->valueIsInBoundOfFormula(storm::utility::regions::evaluateFunction(this->reachProbFunction, point)); } else{ - STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "sample point check with model instantiation not yet implemented"); + //put the values into the dtmc matrix + for( std::pair&>& mappingPair : this->sampleDtmcMapping){ + mappingPair.second.setValue(storm::utility::regions::convertNumber( + storm::utility::regions::evaluateFunction(mappingPair.first, point) + ) + ); + } + std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); + storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); + storm::modelchecker::SparseDtmcPrctlModelChecker modelChecker(*this->sampleDtmc); + std::unique_ptr resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false); + valueInBoundOfFormula=this->valueIsInBoundOfFormula(resultPtr->asExplicitQuantitativeCheckResult().getValueVector()[*this->sampleDtmc->getInitialStates().begin()]); + + //Delete from here + // ConstantType result=resultPtr->asExplicitQuantitativeCheckResult().getValueVector()[*this->sampleDtmc->getInitialStates().begin()]; + // ConstantType otherresult=storm::utility::regions::convertNumber(storm::utility::regions::evaluateFunction(this->reachProbFunction, point)); + + // STORM_LOG_THROW((std::abs(result - otherresult) <= 0.01),storm::exceptions::UnexpectedException, "The results of new DTMC algorithm does not match: " << result << " vs. " << otherresult); + //To here + } if(valueInBoundOfFormula){ @@ -571,12 +817,14 @@ namespace storm { STORM_LOG_THROW(this->hasOnlyLinearFunctions, storm::exceptions::UnexpectedException, "Tried to generate bounds on the probability (only applicable if all functions are linear), but there are nonlinear functions."); //build the mdp and a reachability formula and create a modelchecker std::chrono::high_resolution_clock::time_point timeMDPBuildStart = std::chrono::high_resolution_clock::now(); - storm::models::sparse::Mdp mdp = buildMdpForApproximation(region); + buildMdpForApproximation(region); std::chrono::high_resolution_clock::time_point timeMDPBuildEnd = std::chrono::high_resolution_clock::now(); this->timeMDPBuild += timeMDPBuildEnd-timeMDPBuildStart; std::shared_ptr targetFormulaPtr(new storm::logic::AtomicLabelFormula("target")); storm::logic::EventuallyFormula eventuallyFormula(targetFormulaPtr); - storm::modelchecker::SparseMdpPrctlModelChecker modelChecker(mdp); + storm::modelchecker::SparseMdpPrctlModelChecker modelChecker(*this->approxMdp); + + //Decide whether we should compute lower or upper bounds first. //This does not matter if the current result is unknown. However, let us assume that it is more likely that the final result will be ALLSAT. So we test this first. @@ -619,11 +867,20 @@ namespace storm { for(int i=1; i<=2; ++i){ //perform model checking on the mdp std::unique_ptr resultPtr = modelChecker.computeEventuallyProbabilities(eventuallyFormula,false,opType); - + + //DELETE FROM HERE + // storm::models::sparse::Mdp othermdp = buildMdpForApproximation2(region); + // storm::modelchecker::SparseMdpPrctlModelChecker othermodelChecker(othermdp); + // std::unique_ptr otherresultPtr = othermodelChecker.computeEventuallyProbabilities(eventuallyFormula,false,opType); + // ConstantType origResult=resultPtr->asExplicitQuantitativeCheckResult().getValueVector()[*this->approxMdp->getInitialStates().begin()]; + // ConstantType otherResult=otherresultPtr->asExplicitQuantitativeCheckResult().getValueVector()[*othermdp.getInitialStates().begin()]; + // STORM_LOG_THROW(this->constantTypeComparator.isEqual(origResult,otherResult),storm::exceptions::UnexpectedException, "The results of new mdp algorithm does not match: " << origResult << " vs. " << otherResult); + //TO HERE + //check if the approximation yielded a conclusive result if(opType==storm::logic::OptimalityType::Maximize){ upperBounds = resultPtr->asExplicitQuantitativeCheckResult().getValueVector(); - if(valueIsInBoundOfFormula(upperBounds[*mdp.getInitialStates().begin()])){ + if(valueIsInBoundOfFormula(upperBounds[*this->approxMdp->getInitialStates().begin()])){ if((this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::Less) || (this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::LessEqual)){ region.setCheckResult(RegionCheckResult::ALLSAT); @@ -639,11 +896,11 @@ namespace storm { } //flip the optype for the next iteration opType=storm::logic::OptimalityType::Minimize; - if(i==1) std::cout << " Requiring a second model checker run (with Minimize)" << std::endl; + // if(i==1) std::cout << " Requiring a second model checker run (with Minimize)" << std::endl; } else if(opType==storm::logic::OptimalityType::Minimize){ lowerBounds = resultPtr->asExplicitQuantitativeCheckResult().getValueVector(); - if(valueIsInBoundOfFormula(lowerBounds[*mdp.getInitialStates().begin()])){ + if(valueIsInBoundOfFormula(lowerBounds[*this->approxMdp->getInitialStates().begin()])){ if((this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::Greater) || (this->probabilityOperatorFormula->getComparisonType()== storm::logic::ComparisonType::GreaterEqual)){ region.setCheckResult(RegionCheckResult::ALLSAT); @@ -659,7 +916,7 @@ namespace storm { } //flip the optype for the next iteration opType=storm::logic::OptimalityType::Maximize; - if(i==1) std::cout << " Requiring a second model checker run (with Maximize)" << std::endl; + // if(i==1) std::cout << " Requiring a second model checker run (with Maximize)" << std::endl; } } @@ -668,7 +925,37 @@ namespace storm { } template - storm::models::sparse::Mdp SparseDtmcRegionModelChecker::buildMdpForApproximation(const ParameterRegion& region) { + void SparseDtmcRegionModelChecker::buildMdpForApproximation(const ParameterRegion& region) { + //instantiate the substitutions for the given region + std::vector> substitutions(this->approxMdpSubstitutions.size()); + for(uint_fast64_t substitutionIndex=0; substitutionIndexapproxMdpSubstitutions.size(); ++substitutionIndex){ + for(std::pair const& sub : this->approxMdpSubstitutions[substitutionIndex]){ + switch(sub.second){ + case TypeOfBound::LOWER: + substitutions[substitutionIndex].insert(std::make_pair(sub.first, region.getLowerBound(sub.first))); + break; + case TypeOfBound::UPPER: + substitutions[substitutionIndex].insert(std::make_pair(sub.first, region.getUpperBound(sub.first))); + break; + default: + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected Type of Bound"); + } + } + } + + //now put the values into the mdp matrix + for( std::tuple&, size_t>& mappingTriple : this->approxMdpMapping){ + std::get<1>(mappingTriple).setValue(storm::utility::regions::convertNumber( + storm::utility::regions::evaluateFunction(std::get<0>(mappingTriple),substitutions[std::get<2>(mappingTriple)]) + ) + ); + } + + } + + //DELETEME + template + storm::models::sparse::Mdp SparseDtmcRegionModelChecker::buildMdpForApproximation2(const ParameterRegion& region) { //We are going to build a (non parametric) MDP which has an action for the lower bound and an action for the upper bound of every parameter //The matrix (and the Choice labeling) @@ -877,16 +1164,94 @@ namespace storm { } } - - - + template + template + bool SparseDtmcRegionModelChecker::valueIsInBoundOfFormula(ValueType value){ + STORM_LOG_THROW(this->probabilityOperatorFormula!=nullptr, storm::exceptions::InvalidStateException, "Tried to compare a value to the bound of a formula, but no formula specified."); + double valueAsDouble = storm::utility::regions::convertNumber(value); + // std::cout << "checked whether value: " << value << " (= " << valueAsDouble << " double ) is in bound of formula: " << *this->probabilityOperatorFormula << std::endl; + // std::cout << " (valueAsDouble >= this->probabilityOperatorFormula->getBound()) is " << (valueAsDouble >= this->probabilityOperatorFormula->getBound()) << std::endl; + switch (this->probabilityOperatorFormula->getComparisonType()) { + case storm::logic::ComparisonType::Greater: + return (valueAsDouble > this->probabilityOperatorFormula->getBound()); + case storm::logic::ComparisonType::GreaterEqual: + return (valueAsDouble >= this->probabilityOperatorFormula->getBound()); + case storm::logic::ComparisonType::Less: + return (valueAsDouble < this->probabilityOperatorFormula->getBound()); + case storm::logic::ComparisonType::LessEqual: + return (valueAsDouble <= this->probabilityOperatorFormula->getBound()); + default: + STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported"); + } + } + template + void SparseDtmcRegionModelChecker::printStatisticsToStream(std::ostream& outstream) { + + if(this->probabilityOperatorFormula==nullptr){ + outstream << "Statistic Region Model Checker Statistics Error: No formula specified." << std::endl; + return; + } + + std::chrono::milliseconds timePreprocessingInMilliseconds = std::chrono::duration_cast(this->timePreprocessing); + std::chrono::milliseconds timeInitialStateEliminationInMilliseconds = std::chrono::duration_cast(this->timeInitialStateElimination); + std::chrono::milliseconds timeComputeReachProbFunctionInMilliseconds = std::chrono::duration_cast(this->timeComputeReachProbFunction); + std::chrono::milliseconds timeCheckRegionInMilliseconds = std::chrono::duration_cast(this->timeCheckRegion); + std::chrono::milliseconds timeSammplingInMilliseconds = std::chrono::duration_cast(this->timeSampling); + std::chrono::milliseconds timeApproximationInMilliseconds = std::chrono::duration_cast(this->timeApproximation); + std::chrono::milliseconds timeMDPBuildInMilliseconds = std::chrono::duration_cast(this->timeMDPBuild); + std::chrono::milliseconds timeFullSmtInMilliseconds = std::chrono::duration_cast(this->timeFullSmt); + + std::chrono::high_resolution_clock::duration timeOverall = timePreprocessing + timeCheckRegion; // + ... + std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast(timeOverall); + + std::size_t subsystemTransitions = this->sparseTransitions.getNonzeroEntryCount(); + for(auto const& transition : this->oneStepProbabilities){ + if(!this->parametricTypeComparator.isZero(transition)){ + ++subsystemTransitions; + } + } + + uint_fast64_t numOfSolvedRegions= this->numOfRegionsExistsBoth + this->numOfRegionsAllSat + this->numOfRegionsAllViolated; + + outstream << std::endl << "Statistics Region Model Checker Statistics:" << std::endl; + outstream << "-----------------------------------------------" << std::endl; + outstream << "Model: " << this->model.getNumberOfStates() << " states, " << this->model.getNumberOfTransitions() << " transitions." << std::endl; + outstream << "Reduced model: " << this->subsystem.getNumberOfSetBits() << " states, " << subsystemTransitions << "transitions" << std::endl; + outstream << "Formula: " << *this->probabilityOperatorFormula << std::endl; + outstream << (this->hasOnlyLinearFunctions ? "A" : "Not a") << "ll occuring functions in the model are linear" << std::endl; + outstream << "Number of checked regions: " << this->numOfCheckedRegions << std::endl; + outstream << " Number of solved regions: " << numOfSolvedRegions << "(" << numOfSolvedRegions*100/this->numOfCheckedRegions << "%)" << std::endl; + outstream << " AllSat: " << this->numOfRegionsAllSat << "(" << this->numOfRegionsAllSat*100/this->numOfCheckedRegions << "%)" << std::endl; + outstream << " AllViolated: " << this->numOfRegionsAllViolated << "(" << this->numOfRegionsAllViolated*100/this->numOfCheckedRegions << "%)" << std::endl; + outstream << " ExistsBoth: " << this->numOfRegionsExistsBoth << "(" << this->numOfRegionsExistsBoth*100/this->numOfCheckedRegions << "%)" << std::endl; + outstream << " Unsolved: " << this->numOfCheckedRegions - numOfSolvedRegions << "(" << (this->numOfCheckedRegions - numOfSolvedRegions)*100/this->numOfCheckedRegions << "%)" << std::endl; + outstream << " -- " << std::endl; + outstream << " " << this->numOfRegionsSolvedThroughSampling << " regions solved through Sampling" << std::endl; + outstream << " " << this->numOfRegionsSolvedThroughApproximation << " regions solved through Approximation" << std::endl; + outstream << " " << this->numOfRegionsSolvedThroughSubsystemSmt << " regions solved through SubsystemSmt" << std::endl; + outstream << " " << this->numOfRegionsSolvedThroughFullSmt << " regions solved through FullSmt" << std::endl; + outstream << std::endl; + outstream << "Running times:" << std::endl; + outstream << " " << timeOverallInMilliseconds.count() << "ms overall" << std::endl; + outstream << " " << timePreprocessingInMilliseconds.count() << "ms Preprocessing including... " << std::endl; + outstream << " " << timeInitialStateEliminationInMilliseconds.count() << "ms Initial state elimination including..." << std::endl; + outstream << " " << timeComputeReachProbFunctionInMilliseconds.count() << "ms to compute the reachability probability function" << std::endl; + outstream << " " << timeCheckRegionInMilliseconds.count() << "ms Region Check including... " << std::endl; + outstream << " " << timeSammplingInMilliseconds.count() << "ms Sampling " << std::endl; + outstream << " " << timeApproximationInMilliseconds.count() << "ms Approximation including... " << std::endl; + outstream << " " << timeMDPBuildInMilliseconds.count() << "ms to build the MDP" << std::endl; + outstream << " " << timeFullSmtInMilliseconds.count() << "ms Full Smt solving" << std::endl; + outstream << "-----------------------------------------------" << std::endl; + + } + +#ifdef STORM_HAVE_CARL - - + //DELETEME template<> std::pair,std::vector>> SparseDtmcRegionModelChecker::instantiateFlexibleMatrix(FlexibleMatrix const& matrix, std::vector> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector const& oneStepProbabilities, bool addSelfLoops) const { @@ -984,13 +1349,14 @@ namespace storm { return std::pair, std::vector>>(matrixBuilder.build(), std::move(choiceLabeling)); } + //DELETEME template std::pair,std::vector>> SparseDtmcRegionModelChecker::instantiateFlexibleMatrix(FlexibleMatrix const& matrix, std::vector> const& substitutions, storm::storage::BitVector const& filter, bool addSinkState, std::vector const& oneStepProbabilities, bool addSelfLoops) const{ STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Instantiation of flexible matrix is not supported for this type"); } -#ifdef STORM_HAVE_CARL - + + //DELETEME template<> void SparseDtmcRegionModelChecker::eliminateStates(storm::storage::BitVector& subsystem, FlexibleMatrix& flexibleMatrix, std::vector& oneStepProbabilities, FlexibleMatrix& flexibleBackwardTransitions, storm::storage::BitVector const& initialstates, storm::storage::SparseMatrix const& forwardTransitions, boost::optional> const& statePriorities){ @@ -1057,7 +1423,7 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "elimination of states not suported for this type"); } - + //OBSOLETE ... template<> void SparseDtmcRegionModelChecker::formulateModelWithSMT(storm::solver::Smt2SmtSolver& solver, std::vector& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities){ carl::VariablePool& varPool = carl::VariablePool::getInstance(); @@ -1092,6 +1458,7 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "SMT formulation is not supported for this type"); } + //DELETEME template<> void SparseDtmcRegionModelChecker::restrictProbabilityVariables(storm::solver::Smt2SmtSolver& solver, std::vector const& stateProbVars, storm::storage::BitVector const& subsystem, FlexibleMatrix const& flexibleMatrix, std::vector const& oneStepProbabilities, ParameterRegion const& region, storm::logic::ComparisonType const& compType){ //We are going to build a new (non parametric) MDP which has an action for the lower bound and an action for the upper bound of every parameter @@ -1165,6 +1532,7 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "restricting Probability Variables is not supported for this type"); } + //DELETEME template<> bool SparseDtmcRegionModelChecker::checkRegionOld(storm::logic::Formula const& formula, std::vector parameterRegions){ //Note: this is an 'experimental' implementation @@ -1340,79 +1708,10 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Region check is not supported for this type"); } - template - template - bool SparseDtmcRegionModelChecker::valueIsInBoundOfFormula(ValueType value){ - STORM_LOG_THROW(this->probabilityOperatorFormula!=nullptr, storm::exceptions::InvalidStateException, "Tried to compare a value to the bound of a formula, but no formula specified."); - double valueAsDouble = storm::utility::regions::convertNumber(value); - switch (this->probabilityOperatorFormula->getComparisonType()) { - case storm::logic::ComparisonType::Greater: - return (valueAsDouble > this->probabilityOperatorFormula->getBound()); - case storm::logic::ComparisonType::GreaterEqual: - return (valueAsDouble >= this->probabilityOperatorFormula->getBound()); - case storm::logic::ComparisonType::Less: - return (valueAsDouble < this->probabilityOperatorFormula->getBound()); - case storm::logic::ComparisonType::LessEqual: - return (valueAsDouble <= this->probabilityOperatorFormula->getBound()); - default: - STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "the comparison relation of the formula is not supported"); - } - } + #endif - template - void SparseDtmcRegionModelChecker::printStatisticsToStream(std::ostream& outstream) { - - if(this->probabilityOperatorFormula==nullptr){ - outstream << "Statistic Region Model Checker Statistics Error: No formula specified." << std::endl; - return; - } - - std::chrono::milliseconds timePreprocessingInMilliseconds = std::chrono::duration_cast(this->timePreprocessing); - std::chrono::milliseconds timeInitialStateEliminationInMilliseconds = std::chrono::duration_cast(this->timeInitialStateElimination); - std::chrono::milliseconds timeComputeReachProbFunctionInMilliseconds = std::chrono::duration_cast(this->timeComputeReachProbFunction); - std::chrono::milliseconds timeCheckRegionInMilliseconds = std::chrono::duration_cast(this->timeCheckRegion); - std::chrono::milliseconds timeSammplingInMilliseconds = std::chrono::duration_cast(this->timeSampling); - std::chrono::milliseconds timeApproximationInMilliseconds = std::chrono::duration_cast(this->timeApproximation); - std::chrono::milliseconds timeMDPBuildInMilliseconds = std::chrono::duration_cast(this->timeMDPBuild); - std::chrono::milliseconds timeFullSmtInMilliseconds = std::chrono::duration_cast(this->timeFullSmt); - - std::chrono::high_resolution_clock::duration timeOverall = timePreprocessing + timeCheckRegion; // + ... - std::chrono::milliseconds timeOverallInMilliseconds = std::chrono::duration_cast(timeOverall); - - std::size_t subsystemTransitions = this->sparseTransitions.getNonzeroEntryCount(); - for(auto const& transition : this->oneStepProbabilities){ - if(!this->parametricTypeComparator.isZero(transition)){ - ++subsystemTransitions; - } - } - - outstream << std::endl << "Statistics Region Model Checker Statistics:" << std::endl; - outstream << "-----------------------------------------------" << std::endl; - outstream << "Model: " << this->model.getNumberOfStates() << " states, " << this->model.getNumberOfTransitions() << " transitions." << std::endl; - outstream << "Reduced model: " << this->subsystem.getNumberOfSetBits() << " states, " << subsystemTransitions << "transitions" << std::endl; - outstream << "Formula: " << *this->probabilityOperatorFormula << std::endl; - outstream << (this->hasOnlyLinearFunctions ? "A" : "Not a") << "ll occuring functions in the model are linear" << std::endl; - outstream << "Number of checked regions: " << this->numOfCheckedRegions << " with..." << std::endl; - outstream << " " << this->numOfRegionsSolvedThroughSampling << " regions solved through Sampling" << std::endl; - outstream << " " << this->numOfRegionsSolvedThroughApproximation << " regions solved through Approximation" << std::endl; - outstream << " " << this->numOfRegionsSolvedThroughSubsystemSmt << " regions solved through SubsystemSmt" << std::endl; - outstream << " " << this->numOfRegionsSolvedThroughFullSmt << " regions solved through FullSmt" << std::endl; - outstream << "Running times:" << std::endl; - outstream << " " << timeOverallInMilliseconds.count() << "ms overall" << std::endl; - outstream << " " << timePreprocessingInMilliseconds.count() << "ms Preprocessing including... " << std::endl; - outstream << " " << timeInitialStateEliminationInMilliseconds.count() << "ms Initial state elimination including..." << std::endl; - outstream << " " << timeComputeReachProbFunctionInMilliseconds.count() << "ms to compute the reachability probability function" << std::endl; - outstream << " " << timeCheckRegionInMilliseconds.count() << "ms Region Check including... " << std::endl; - outstream << " " << timeSammplingInMilliseconds.count() << "ms Sampling " << std::endl; - outstream << " " << timeApproximationInMilliseconds.count() << "ms Approximation including... " << std::endl; - outstream << " " << timeMDPBuildInMilliseconds.count() << "ms to build the MDP" << std::endl; - outstream << " " << timeFullSmtInMilliseconds.count() << "ms Full Smt solving" << std::endl; - outstream << "-----------------------------------------------" << std::endl; - - } -#endif #ifdef STORM_HAVE_CARL template class SparseDtmcRegionModelChecker; diff --git a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h index 09b2d6928..5f2de241a 100644 --- a/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h +++ b/src/modelchecker/reachability/SparseDtmcRegionModelChecker.h @@ -20,7 +20,7 @@ namespace storm { //The type of variables and bounds depends on the template arguments typedef typename std::conditional<(std::is_same::value), storm::Variable,std::nullptr_t>::type VariableType; - typedef typename std::conditional<(std::is_same::value), storm::RationalFunction::CoeffType,std::nullptr_t>::type BoundType; + typedef typename std::conditional<(std::is_same::value), storm::Coefficient,std::nullptr_t>::type BoundType; /*! * The possible results for a single region @@ -149,6 +149,9 @@ namespace storm { * Prints statistical information (mostly running times) to the given stream. */ void printStatisticsToStream(std::ostream& outstream); + + + private: typedef typename storm::modelchecker::SparseDtmcEliminationModelChecker::FlexibleSparseMatrix FlexibleMatrix; @@ -207,6 +210,44 @@ namespace storm { storm::storage::sparse::state_type const& initState ); + + + enum class TypeOfBound { + LOWER, + UPPER + }; + + /*! + * Initializes a DTMC that can be used to get the probability result for a certain parameter evaluation. + * To quickly insert different evaluations, we provide a dtmc with some dummy values in the transition entries. + * Furthermore, another matrix is given that is used to connect matrix entries and transition functions. + * In a similar way an MDP is initialized. The Mdp can be used to approximate the probabilities. + * In addition to matrix entries and transition functions, we specify which variables and which bounds needs to be substituted + * + * + * @param sampleDtmc This dtmc can represent different instantiations of the considered pDtmc + * @param sampleDtmcMapping Connects the entries of the sampleDtmc Matrix with the corresponding transitionFunctions + * @param approxMdp This mdp will later be used to approximate the reachability probabilities + * @param approxMdpMapping Connects the entries of the approxMdp Matrix with the corresponding transitionFunctions and a substitution (given as an index in the approxMdpSubstitutionsVector) + * @param approxMdpSubstitutions contains the substitutions of the parameters + * + * @param subsys the states of the flexTransitions that are still part of the pDTMC (i.e. that have not been eliminated) + * @param flexTransitions the transitions of the pDTMC + * @param oneStepProbs the probabilities to move to a target state + * @param initState the initial state of the pDtmc + */ + void initializeSampleDtmcAndApproxMdp( + std::shared_ptr> & sampleDtmc, + std::vector&>> & sampleDtmcMapping, + std::shared_ptr> & approxMdp, + std::vector&, std::size_t>> & approxMdpMapping, + std::vector> & approxMdpSubstitutions, + storm::storage::BitVector const& subsys, + storm::storage::SparseMatrix const& transitions, + std::vector const& oneStepProbs, + storm::storage::sparse::state_type const& initState + ); + //Computes the reachability probability function by performing state elimination ParametricType computeReachProbFunction( storm::storage::BitVector const& subsys, @@ -254,9 +295,16 @@ namespace storm { /*! - * Actually builds the mdp that is used to obtain bounds on the maximal/minimal reachability probability + * Builds the mdp that is used to obtain bounds on the maximal/minimal reachability probability + * The result is stored in this->approxMdp + */ + void buildMdpForApproximation(ParameterRegion const& region); + + /*! + * Builds the mdp that is used to obtain bounds on the maximal/minimal reachability probability + */ - storm::models::sparse::Mdp buildMdpForApproximation(ParameterRegion const& region); + storm::models::sparse::Mdp buildMdpForApproximation2(ParameterRegion const& region); /*! * Starts the SMTSolver to get the result. @@ -303,6 +351,20 @@ namespace storm { storm::storage::BitVector subsystem; // a flag that is true if there are only linear functions at transitions of the model bool hasOnlyLinearFunctions; + // the dtmc that can be instantiated to check the value at a certain point + std::shared_ptr> sampleDtmc; + // a vector that links entries of the dtmc matrix with the corresponding transition functions (for fast instantiation) + std::vector&>> sampleDtmcMapping; + // the Mdp that is used to approximate the probability values + std::shared_ptr> approxMdp; + // a vector that links entries of the mdp matrix with the corresponding transition functions and substitutions (given as an index of this->approxMdpSubstitutions vector) + std::vector&, std::size_t>> approxMdpMapping; + // stores the different substitutions of the variables + std::vector> approxMdpSubstitutions; + + + + // The function for the reachability probability in the initial state ParametricType reachProbFunction; @@ -310,11 +372,16 @@ namespace storm { // runtimes and other information for statistics. uint_fast64_t numOfCheckedRegions; + uint_fast64_t numOfRegionsSolvedThroughSampling; uint_fast64_t numOfRegionsSolvedThroughApproximation; uint_fast64_t numOfRegionsSolvedThroughSubsystemSmt; uint_fast64_t numOfRegionsSolvedThroughFullSmt; + uint_fast64_t numOfRegionsExistsBoth; + uint_fast64_t numOfRegionsAllSat; + uint_fast64_t numOfRegionsAllViolated; + std::chrono::high_resolution_clock::duration timePreprocessing; std::chrono::high_resolution_clock::duration timeInitialStateElimination; std::chrono::high_resolution_clock::duration timeComputeReachProbFunction; diff --git a/src/utility/constants.cpp b/src/utility/constants.cpp index d4625b578..fde1f5bfd 100644 --- a/src/utility/constants.cpp +++ b/src/utility/constants.cpp @@ -253,6 +253,9 @@ namespace storm { template RationalFunction pow(RationalFunction const& value, uint_fast64_t exponent); + template Coefficient one(); + template Coefficient zero(); + template Polynomial one(); template Polynomial zero(); template RationalFunction simplify(RationalFunction value); diff --git a/src/utility/regions.cpp b/src/utility/regions.cpp index fd1717014..14215b0d0 100644 --- a/src/utility/regions.cpp +++ b/src/utility/regions.cpp @@ -102,7 +102,7 @@ namespace storm { template<> - storm::RationalFunction::CoeffType convertNumber(double const& number, bool const& roundDown, double const& precision){ + storm::Coefficient convertNumber(double const& number, bool const& roundDown, double const& precision){ double actualPrecision = (precision==0.0 ? storm::settings::generalSettings().getPrecision() : precision); uint_fast64_t denominator = 1.0/actualPrecision; uint_fast64_t numerator; @@ -111,14 +111,14 @@ namespace storm { } else{ numerator = std::ceil(number*denominator); } - storm::RationalFunction::CoeffType result(numerator); + storm::Coefficient result(numerator); result = result/denominator; return result; } template<> storm::RationalFunction convertNumber(double const& number, bool const& roundDown, double const& precision){ - return storm::RationalFunction(convertNumber(number, roundDown, precision)); + return storm::RationalFunction(convertNumber(number, roundDown, precision)); } template<> @@ -255,6 +255,19 @@ namespace storm { solver->add(variable, value); } + template<> + storm::RationalFunction getNewFunction(storm::Coefficient initialValue) { + std::shared_ptr>> cache(new carl::Cache>()); + return storm::RationalFunction(storm::RationalFunction::PolyType(storm::RationalFunction::PolyType::PolyType(initialValue), cache)); + } + + template<> + storm::RationalFunction getNewFunction(storm::Variable initialValue) { + std::shared_ptr>> cache(new carl::Cache>()); + return storm::RationalFunction(storm::RationalFunction::PolyType(storm::RationalFunction::PolyType::PolyType(initialValue), cache)); + } + + //explicit instantiations template double convertNumber(double const& number, bool const& roundDown, double const& precision); @@ -262,8 +275,8 @@ namespace storm { template class RegionParser; template storm::RationalFunction convertNumber(double const& number, bool const& roundDown, double const& precision); - template storm::RationalFunction::CoeffType convertNumber(double const& number, bool const& roundDown, double const& precision); - template double convertNumber(storm::RationalFunction::CoeffType const& number, bool const& roundDown, double const& precision); + template storm::Coefficient convertNumber(double const& number, bool const& roundDown, double const& precision); + template double convertNumber(storm::Coefficient const& number, bool const& roundDown, double const& precision); template cln::cl_RA convertNumber(cln::cl_RA const& number, bool const& roundDown, double const& precision); template storm::Variable getVariableFromString(std::string variableString); @@ -277,6 +290,8 @@ namespace storm { template void addParameterBoundsToSmtSolver(std::shared_ptr solver, storm::Variable const& variable, storm::logic::ComparisonType relation, cln::cl_RA const& bound); template void addBoolVariableToSmtSolver(std::shared_ptr solver, storm::Variable const& variable, bool value); + template storm::RationalFunction getNewFunction(storm::Coefficient initialValue); + template storm::RationalFunction getNewFunction(storm::Variable initialValue); #endif } } diff --git a/src/utility/regions.h b/src/utility/regions.h index bc1a76307..a77bc1bdc 100644 --- a/src/utility/regions.h +++ b/src/utility/regions.h @@ -155,6 +155,14 @@ namespace storm { */ template void addBoolVariableToSmtSolver(std::shared_ptr solver, VariableType const& variable, bool value); + + + /*! + * Returns a new function that initially has the given value (which might be a constant or a variable) + * + */ + template + ParametricType getNewFunction(ValueType initialValue); } } }