SparseParameterLiftingModelChecker.cpp

Go to the documentation of this file.

#include "storm-pars/modelchecker/region/SparseParameterLiftingModelChecker.h"
 
#include <storm-pars/analysis/MonotonicityChecker.h>
#include <boost/container/flat_set.hpp>
#include <queue>
 
#include "storm/adapters/RationalFunctionAdapter.h"
#include "storm/logic/FragmentSpecification.h"
#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h"
#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
#include "storm/models/sparse/Dtmc.h"
#include "storm/models/sparse/Mdp.h"
#include "storm/models/sparse/StandardRewardModel.h"
#include "storm/utility/Stopwatch.h"
#include "storm/utility/vector.h"
 
#include "storm/exceptions/InvalidArgumentException.h"
#include "storm/exceptions/NotSupportedException.h"
 
namespace storm {
namespace modelchecker {
 
template<typename SparseModelType, typename ConstantType>
SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::SparseParameterLiftingModelChecker() {
    // Intentionally left empty
}
 
template<typename SparseModelType, typename ConstantType>
void SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyFormula(
    Environment const& env, storm::modelchecker::CheckTask<storm::logic::Formula, typename SparseModelType::ValueType> const& checkTask) {
    currentFormula = checkTask.getFormula().asSharedPointer();
    currentCheckTask = std::make_unique<storm::modelchecker::CheckTask<storm::logic::Formula, ConstantType>>(
        checkTask.substituteFormula(*currentFormula).template convertValueType<ConstantType>());
 
    if (currentCheckTask->getFormula().isProbabilityOperatorFormula()) {
        auto const& probOpFormula = currentCheckTask->getFormula().asProbabilityOperatorFormula();
        if (probOpFormula.getSubformula().isBoundedUntilFormula()) {
            specifyBoundedUntilFormula(currentCheckTask->substituteFormula(probOpFormula.getSubformula().asBoundedUntilFormula()));
        } else if (probOpFormula.getSubformula().isUntilFormula()) {
            specifyUntilFormula(env, currentCheckTask->substituteFormula(probOpFormula.getSubformula().asUntilFormula()));
        } else if (probOpFormula.getSubformula().isEventuallyFormula()) {
            specifyReachabilityProbabilityFormula(env, currentCheckTask->substituteFormula(probOpFormula.getSubformula().asEventuallyFormula()));
        } else {
            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property.");
        }
    } else if (currentCheckTask->getFormula().isRewardOperatorFormula()) {
        auto const& rewOpFormula = currentCheckTask->getFormula().asRewardOperatorFormula();
        if (rewOpFormula.getSubformula().isEventuallyFormula()) {
            specifyReachabilityRewardFormula(env, currentCheckTask->substituteFormula(rewOpFormula.getSubformula().asEventuallyFormula()));
        } else if (rewOpFormula.getSubformula().isCumulativeRewardFormula()) {
            specifyCumulativeRewardFormula(currentCheckTask->substituteFormula(rewOpFormula.getSubformula().asCumulativeRewardFormula()));
        }
    }
}
 
template<typename SparseModelType, typename ConstantType>
RegionResult SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::analyzeRegion(
    Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region, RegionResultHypothesis const& hypothesis,
    RegionResult const& initialResult, bool sampleVerticesOfRegion,
    std::shared_ptr<storm::analysis::LocalMonotonicityResult<typename RegionModelChecker<typename SparseModelType::ValueType>::VariableType>>
        localMonotonicityResult) {
    typedef typename RegionModelChecker<typename SparseModelType::ValueType>::VariableType VariableType;
    typedef typename storm::analysis::MonotonicityResult<VariableType>::Monotonicity Monotonicity;
    typedef typename storm::utility::parametric::Valuation<typename SparseModelType::ValueType> Valuation;
    typedef typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::CoefficientType CoefficientType;
    STORM_LOG_THROW(this->currentCheckTask->isOnlyInitialStatesRelevantSet(), storm::exceptions::NotSupportedException,
                    "Analyzing regions with parameter lifting requires a property where only the value in the initial states is relevant.");
    STORM_LOG_THROW(this->currentCheckTask->isBoundSet(), storm::exceptions::NotSupportedException,
                    "Analyzing regions with parameter lifting requires a bounded property.");
    STORM_LOG_THROW(this->parametricModel->getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException,
                    "Analyzing regions with parameter lifting requires a model with a single initial state.");
 
    RegionResult result = initialResult;
 
    // Check if we need to check the formula on one point to decide whether to show AllSat or AllViolated
    if (hypothesis == RegionResultHypothesis::Unknown && result == RegionResult::Unknown) {
        result = getInstantiationChecker()
                         .check(env, region.getCenterPoint())
                         ->asExplicitQualitativeCheckResult()[*this->parametricModel->getInitialStates().begin()]
                     ? RegionResult::CenterSat
                     : RegionResult::CenterViolated;
    }
 
    bool existsSat = (hypothesis == RegionResultHypothesis::AllSat || result == RegionResult::ExistsSat || result == RegionResult::CenterSat);
    bool existsViolated =
        (hypothesis == RegionResultHypothesis::AllViolated || result == RegionResult::ExistsViolated || result == RegionResult::CenterViolated);
 
    // Here we check on global monotonicity
    if (localMonotonicityResult != nullptr && localMonotonicityResult->isDone()) {
        // Try to check it with a global monotonicity result
        auto monRes = localMonotonicityResult->getGlobalMonotonicityResult();
        bool lowerBound = isLowerBound(this->currentCheckTask->getBound().comparisonType);
 
        if (monRes->isDone() && monRes->isAllMonotonicity()) {
            // Build valuations
            auto monMap = monRes->getMonotonicityResult();
            Valuation valuationToCheckSat;
            Valuation valuationToCheckViolated;
            for (auto var : region.getVariables()) {
                auto monVar = monMap[var];
                if (monVar == Monotonicity::Constant) {
                    valuationToCheckSat.insert(std::pair<VariableType, CoefficientType>(var, region.getLowerBoundary(var)));
                    valuationToCheckViolated.insert(std::pair<VariableType, CoefficientType>(var, region.getLowerBoundary(var)));
                } else if (monVar == Monotonicity::Decr) {
                    if (lowerBound) {
                        valuationToCheckSat.insert(std::pair<VariableType, CoefficientType>(var, region.getUpperBoundary(var)));
                        valuationToCheckViolated.insert(std::pair<VariableType, CoefficientType>(var, region.getLowerBoundary(var)));
                    } else {
                        valuationToCheckSat.insert(std::pair<VariableType, CoefficientType>(var, region.getLowerBoundary(var)));
                        valuationToCheckViolated.insert(std::pair<VariableType, CoefficientType>(var, region.getUpperBoundary(var)));
                    }
                } else if (monVar == Monotonicity::Incr) {
                    if (lowerBound) {
                        valuationToCheckSat.insert(std::pair<VariableType, CoefficientType>(var, region.getLowerBoundary(var)));
                        valuationToCheckViolated.insert(std::pair<VariableType, CoefficientType>(var, region.getUpperBoundary(var)));
                    } else {
                        valuationToCheckSat.insert(std::pair<VariableType, CoefficientType>(var, region.getUpperBoundary(var)));
                        valuationToCheckViolated.insert(std::pair<VariableType, CoefficientType>(var, region.getLowerBoundary(var)));
                    }
                }
            }
 
            // Check for result
            if (existsSat && getInstantiationCheckerSAT()
                                 .check(env, valuationToCheckSat)
                                 ->asExplicitQualitativeCheckResult()[*this->parametricModel->getInitialStates().begin()]) {
                STORM_LOG_INFO("Region " << region << " is AllSat, discovered with instantiation checker on " << valuationToCheckSat
                                         << " and help of monotonicity\n");
                RegionModelChecker<typename SparseModelType::ValueType>::numberOfRegionsKnownThroughMonotonicity++;
                return RegionResult::AllSat;
            }
 
            if (existsViolated && !getInstantiationCheckerVIO()
                                       .check(env, valuationToCheckViolated)
                                       ->asExplicitQualitativeCheckResult()[*this->parametricModel->getInitialStates().begin()]) {
                STORM_LOG_INFO("Region " << region << " is AllViolated, discovered with instantiation checker on " << valuationToCheckViolated
                                         << " and help of monotonicity\n");
                RegionModelChecker<typename SparseModelType::ValueType>::numberOfRegionsKnownThroughMonotonicity++;
                return RegionResult::AllViolated;
            }
 
            return RegionResult::ExistsBoth;
        }
    }
 
    // Try to prove AllSat or AllViolated, depending on the hypothesis or the current result
    if (existsSat) {
        // show AllSat:
        storm::solver::OptimizationDirection parameterOptimizationDirection = isLowerBound(this->currentCheckTask->getBound().comparisonType)
                                                                                  ? storm::solver::OptimizationDirection::Minimize
                                                                                  : storm::solver::OptimizationDirection::Maximize;
        auto checkResult = this->check(env, region, parameterOptimizationDirection, localMonotonicityResult);
        if (checkResult->asExplicitQualitativeCheckResult()[*this->parametricModel->getInitialStates().begin()]) {
            result = RegionResult::AllSat;
        } else if (sampleVerticesOfRegion) {
            result = sampleVertices(env, region, result);
        }
    } else if (existsViolated) {
        // show AllViolated:
        storm::solver::OptimizationDirection parameterOptimizationDirection = isLowerBound(this->currentCheckTask->getBound().comparisonType)
                                                                                  ? storm::solver::OptimizationDirection::Maximize
                                                                                  : storm::solver::OptimizationDirection::Minimize;
        auto checkResult = this->check(env, region, parameterOptimizationDirection, localMonotonicityResult);
        if (!checkResult->asExplicitQualitativeCheckResult()[*this->parametricModel->getInitialStates().begin()]) {
            result = RegionResult::AllViolated;
        } else if (sampleVerticesOfRegion) {
            result = sampleVertices(env, region, result);
        }
    } else {
        STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "When analyzing a region, an invalid initial result was given: " << initialResult);
    }
    return result;
}
 
template<typename SparseModelType, typename ConstantType>
RegionResult SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::sampleVertices(
    Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region, RegionResult const& initialResult) {
    RegionResult result = initialResult;
 
    if (result == RegionResult::AllSat || result == RegionResult::AllViolated) {
        return result;
    }
 
    bool hasSatPoint = result == RegionResult::ExistsSat || result == RegionResult::CenterSat;
    bool hasViolatedPoint = result == RegionResult::ExistsViolated || result == RegionResult::CenterViolated;
 
    // Check if there is a point in the region for which the property is satisfied
    auto vertices = region.getVerticesOfRegion(region.getVariables());
    auto vertexIt = vertices.begin();
    while (vertexIt != vertices.end() && !(hasSatPoint && hasViolatedPoint)) {
        if (getInstantiationChecker().check(env, *vertexIt)->asExplicitQualitativeCheckResult()[*this->parametricModel->getInitialStates().begin()]) {
            hasSatPoint = true;
        } else {
            hasViolatedPoint = true;
        }
        ++vertexIt;
    }
 
    if (hasSatPoint) {
        if (hasViolatedPoint) {
            result = RegionResult::ExistsBoth;
        } else if (result != RegionResult::CenterSat) {
            result = RegionResult::ExistsSat;
        }
    } else if (hasViolatedPoint && result != RegionResult::CenterViolated) {
        result = RegionResult::ExistsViolated;
    }
 
    return result;
}
 
template<typename SparseModelType, typename ConstantType>
std::unique_ptr<CheckResult> SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::check(
    Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region,
    storm::solver::OptimizationDirection const& dirForParameters,
    std::shared_ptr<storm::analysis::LocalMonotonicityResult<typename RegionModelChecker<typename SparseModelType::ValueType>::VariableType>>
        localMonotonicityResult) {
    auto quantitativeResult = computeQuantitativeValues(env, region, dirForParameters, localMonotonicityResult);
    lastValue = quantitativeResult->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()];
    if (currentCheckTask->getFormula().hasQuantitativeResult()) {
        return quantitativeResult;
    } else {
        return quantitativeResult->template asExplicitQuantitativeCheckResult<ConstantType>().compareAgainstBound(
            this->currentCheckTask->getFormula().asOperatorFormula().getComparisonType(),
            this->currentCheckTask->getFormula().asOperatorFormula().template getThresholdAs<ConstantType>());
    }
}
 
template<typename SparseModelType, typename ConstantType>
std::unique_ptr<QuantitativeCheckResult<ConstantType>> SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::getBound(
    Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region,
    storm::solver::OptimizationDirection const& dirForParameters,
    std::shared_ptr<storm::analysis::LocalMonotonicityResult<typename RegionModelChecker<typename SparseModelType::ValueType>::VariableType>>
        localMonotonicityResult) {
    STORM_LOG_WARN_COND(this->currentCheckTask->getFormula().hasQuantitativeResult(), "Computing quantitative bounds for a qualitative formula...");
    return std::make_unique<ExplicitQuantitativeCheckResult<ConstantType>>(std::move(
        computeQuantitativeValues(env, region, dirForParameters, localMonotonicityResult)->template asExplicitQuantitativeCheckResult<ConstantType>()));
}
 
template<typename SparseModelType, typename ConstantType>
typename SparseModelType::ValueType SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::getBoundAtInitState(
    Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region,
    storm::solver::OptimizationDirection const& dirForParameters) {
    STORM_LOG_THROW(this->parametricModel->getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException,
                    "Getting a bound at the initial state requires a model with a single initial state.");
    return storm::utility::convertNumber<typename SparseModelType::ValueType>(
        getBound(env, region, dirForParameters)
            ->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()]);
}
 
template<typename SparseModelType, typename ConstantType>
storm::modelchecker::SparseInstantiationModelChecker<SparseModelType, ConstantType>&
SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::getInstantiationCheckerSAT() {
    return getInstantiationChecker();
}
 
template<typename SparseModelType, typename ConstantType>
storm::modelchecker::SparseInstantiationModelChecker<SparseModelType, ConstantType>&
SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::getInstantiationCheckerVIO() {
    return getInstantiationChecker();
}
 
template<typename SparseModelType, typename ConstantType>
struct RegionBound {
    typedef typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::VariableType VariableType;
 
    RegionBound(storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& r, std::shared_ptr<storm::analysis::Order> o,
                std::shared_ptr<storm::analysis::LocalMonotonicityResult<VariableType>> l, boost::optional<ConstantType> const& b)
        : region(r), order(o), localMonRes(l), bound(b) {}
    storm::storage::ParameterRegion<typename SparseModelType::ValueType> region;
    std::shared_ptr<storm::analysis::Order> order;
    std::shared_ptr<storm::analysis::LocalMonotonicityResult<VariableType>> localMonRes;
    boost::optional<ConstantType> bound;
};
 
template<typename SparseModelType, typename ConstantType>
std::pair<ConstantType, typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation>
SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::computeExtremalValue(
    Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dir,
    typename SparseModelType::ValueType const& precision, bool absolutePrecision, boost::optional<ConstantType> const& initialValue,
    std::optional<storm::logic::Bound> const& terminationBound) {
    typedef typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::CoefficientType CoefficientType;
    typedef typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation Valuation;
    STORM_LOG_THROW(this->parametricModel->getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException,
                    "Getting extremal values at the initial state requires a model with a single initial state.");
    STORM_LOG_THROW(!this->currentCheckTask->isBoundSet(), storm::exceptions::NotSupportedException,
                    "Computing extremal values with parameter lifting requires no bound on the operator.");
    bool const useMonotonicity = this->isUseMonotonicitySet();
    bool const minimize = storm::solver::minimize(dir);
 
    // Comparator for the region queue
    auto cmp = storm::solver::minimize(dir) ? [](RegionBound<SparseModelType, ConstantType> const& lhs,
                                                 RegionBound<SparseModelType, ConstantType> const& rhs) { return lhs.bound > rhs.bound; }
                                            : [](RegionBound<SparseModelType, ConstantType> const& lhs, RegionBound<SparseModelType, ConstantType> const& rhs) {
                                                  return lhs.bound < rhs.bound;
                                              };
    std::priority_queue<RegionBound<SparseModelType, ConstantType>, std::vector<RegionBound<SparseModelType, ConstantType>>, decltype(cmp)> regionQueue(cmp);
    storm::utility::Stopwatch initialWatch(true);
 
    storm::utility::Stopwatch boundsWatch(false);
    auto numberOfPLACallsBounds = 0;
    boost::optional<ConstantType> initBound;
    if (useMonotonicity) {
        if (this->isUseBoundsSet()) {
            numberOfPLACallsBounds++;
            numberOfPLACallsBounds++;
            auto minBound = getBound(env, region, storm::solver::OptimizationDirection::Minimize, nullptr)
                                ->template asExplicitQuantitativeCheckResult<ConstantType>()
                                .getValueVector();
            auto maxBound = getBound(env, region, storm::solver::OptimizationDirection::Maximize, nullptr)
                                ->template asExplicitQuantitativeCheckResult<ConstantType>()
                                .getValueVector();
            if (minimize) {
                initBound = minBound[*this->parametricModel->getInitialStates().begin()];
            } else {
                initBound = maxBound[*this->parametricModel->getInitialStates().begin()];
            }
            orderExtender->setMinValuesInit(minBound);
            orderExtender->setMaxValuesInit(maxBound);
        }
        auto order = this->extendOrder(nullptr, region);
        auto monRes = std::shared_ptr<storm::analysis::LocalMonotonicityResult<VariableType>>(
            new storm::analysis::LocalMonotonicityResult<VariableType>(order->getNumberOfStates()));
        storm::utility::Stopwatch monotonicityWatch(true);
        this->extendLocalMonotonicityResult(region, order, monRes);
        monotonicityWatch.stop();
        STORM_LOG_INFO("\nTotal time for monotonicity checking: " << monotonicityWatch << ".\n\n");
 
        regionQueue.emplace(region, order, monRes, initBound);
    } else {
        regionQueue.emplace(region, nullptr, nullptr, initBound);
    }
 
    // The results
    boost::optional<ConstantType> value;
    Valuation valuation;
    if (!initialValue) {
        auto init = getGoodInitialPoint(env, region, dir, regionQueue.top().localMonRes);
        value = storm::utility::convertNumber<ConstantType>(init.first);
        valuation = std::move(init.second);
    } else {
        value = initialValue;
    }
 
    initialWatch.stop();
    STORM_LOG_INFO("\nTotal time for initial points: " << initialWatch << ".\n\n");
    if (!initialValue) {
        STORM_LOG_INFO("Initial value: " << value.get() << " at " << valuation);
    } else {
        STORM_LOG_INFO("Initial value: " << value.get() << " as provided by the user");
    }
 
    if (terminationBound != std::nullopt && !terminationBound.value().isSatisfied(value.get())) {
        return std::make_pair(storm::utility::convertNumber<ConstantType>(value.get()), valuation);
    }
 
    auto numberOfSplits = 0;
    auto numberOfPLACalls = 0;
    auto numberOfOrderCopies = 0;
    auto numberOfMonResCopies = 0;
    storm::utility::Stopwatch loopWatch(true);
    if (!(useMonotonicity && regionQueue.top().localMonRes->getGlobalMonotonicityResult()->isDone() &&
          regionQueue.top().localMonRes->getGlobalMonotonicityResult()->isAllMonotonicity())) {
        // Doing the extremal computation, only when we don't use monotonicity or there are possibly not monotone variables.
        auto totalArea = storm::utility::convertNumber<ConstantType>(region.area());
        auto coveredArea = storm::utility::zero<ConstantType>();
        while (!regionQueue.empty()) {
            assert(value);
            auto currRegion = regionQueue.top().region;
            auto order = regionQueue.top().order;
            auto localMonotonicityResult = regionQueue.top().localMonRes;
            auto currBound = regionQueue.top().bound;
            STORM_LOG_INFO("Currently looking at region: " << currRegion);
 
            std::vector<storm::storage::ParameterRegion<typename SparseModelType::ValueType>> newRegions;
 
            // Check whether this region needs further investigation based on the bound of the parent region
            bool investigateBounds = !currBound;
 
            if (!investigateBounds && absolutePrecision) {
                investigateBounds = (minimize && currBound.get() < value.get() - storm::utility::convertNumber<ConstantType>(precision)) ||
                                    (!minimize && currBound.get() > value.get() + storm::utility::convertNumber<ConstantType>(precision));
            } else if (!investigateBounds && !absolutePrecision) {
                investigateBounds = (minimize && (currBound.get() * (1 + storm::utility::convertNumber<ConstantType>(precision)) < value.get())) ||
                                    (!minimize && (currBound.get() * (1 - storm::utility::convertNumber<ConstantType>(precision)) > value.get()));
            }
 
            if (investigateBounds) {
                numberOfPLACalls++;
                auto bounds =
                    getBound(env, currRegion, dir, localMonotonicityResult)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
                // TODO this looks dangerous
                currBound = bounds[*this->parametricModel->getInitialStates().begin()];
                // Check whether this region needs further investigation based on the bound of this region
                bool lookAtRegion;
                if (absolutePrecision) {
                    lookAtRegion = (minimize && currBound.get() < value.get() - storm::utility::convertNumber<ConstantType>(precision)) ||
                                   (!minimize && currBound.get() > value.get() + storm::utility::convertNumber<ConstantType>(precision));
                } else {
                    lookAtRegion = (minimize && (currBound.get() * (1 + storm::utility::convertNumber<ConstantType>(precision)) < value.get())) ||
                                   (!minimize && (currBound.get() * (1 - storm::utility::convertNumber<ConstantType>(precision)) > value.get()));
                }
                if (lookAtRegion) {
                    if (useMonotonicity) {
                        // Continue extending order/monotonicity result
                        bool changedOrder = false;
                        if (!order->getDoneBuilding() && orderExtender->isHope(order)) {
                            if (numberOfCopiesOrder[order] != 1) {
                                numberOfCopiesOrder[order]--;
                                order = copyOrder(order);
                                numberOfOrderCopies++;
                            } else {
                                assert(numberOfCopiesOrder[order] == 1);
                            }
                            this->extendOrder(order, currRegion);
                            changedOrder = true;
                        }
                        if (changedOrder) {
                            assert(!localMonotonicityResult->isDone());
 
                            if (numberOfCopiesMonRes[localMonotonicityResult] != 1) {
                                numberOfCopiesMonRes[localMonotonicityResult]--;
                                localMonotonicityResult = localMonotonicityResult->copy();
                                numberOfMonResCopies++;
                            } else {
                                assert(numberOfCopiesMonRes[localMonotonicityResult] == 1);
                            }
                            this->extendLocalMonotonicityResult(currRegion, order, localMonotonicityResult);
                            STORM_LOG_INFO("Order and monotonicity result got extended");
                        }
                    }
 
                    // Check whether this region contains a new 'good' value and set this value
                    auto point =
                        useMonotonicity ? currRegion.getPoint(dir, *(localMonotonicityResult->getGlobalMonotonicityResult())) : currRegion.getCenterPoint();
                    auto currValue = getInstantiationChecker()
                                         .check(env, point)
                                         ->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()];
                    if (!value || (minimize ? currValue <= value.get() : currValue >= value.get())) {
                        value = currValue;
                        valuation = point;
                        if (terminationBound != std::nullopt && !terminationBound.value().isSatisfied(value.get())) {
                            return std::make_pair(storm::utility::convertNumber<ConstantType>(value.get()), valuation);
                        }
                    }
 
                    bool splitRegion;
                    if (absolutePrecision) {
                        splitRegion = (minimize && currBound.get() < value.get() - storm::utility::convertNumber<ConstantType>(precision)) ||
                                      (!minimize && currBound.get() > value.get() + storm::utility::convertNumber<ConstantType>(precision));
                    } else {
                        splitRegion = (minimize && (currBound.get() * (1 + storm::utility::convertNumber<ConstantType>(precision)) < value.get())) ||
                                      (!minimize && (currBound.get() * (1 - storm::utility::convertNumber<ConstantType>(precision)) > value.get()));
                    }
                    if (splitRegion) {
                        // We will split the region in this case, but first we set the bounds to extend the order for the new regions.
                        if (useMonotonicity && this->isUseBoundsSet() && !order->getDoneBuilding()) {
                            boundsWatch.start();
                            numberOfPLACallsBounds++;
                            if (minimize) {
                                orderExtender->setMinMaxValues(
                                    order, bounds,
                                    getBound(env, currRegion, storm::solver::OptimizationDirection::Maximize, localMonotonicityResult)
                                        ->template asExplicitQuantitativeCheckResult<ConstantType>()
                                        .getValueVector());
                            } else {
                                orderExtender->setMinMaxValues(
                                    order,
                                    getBound(env, currRegion, storm::solver::OptimizationDirection::Maximize, localMonotonicityResult)
                                        ->template asExplicitQuantitativeCheckResult<ConstantType>()
                                        .getValueVector(),
                                    bounds);
                            }
                            boundsWatch.stop();
                        }
                        // Now split the region
                        if (useMonotonicity) {
                            this->splitSmart(currRegion, newRegions, *(localMonotonicityResult->getGlobalMonotonicityResult()), true);
                        } else if (this->isRegionSplitEstimateSupported()) {
                            auto empty = storm::analysis::MonotonicityResult<VariableType>();
                            this->splitSmart(currRegion, newRegions, empty, true);
                        } else {
                            currRegion.split(currRegion.getCenterPoint(), newRegions);
                        }
                    }
                }
            }
 
            if (newRegions.empty()) {
                // When the newRegions is empty we are done with the current region
                coveredArea += storm::utility::convertNumber<ConstantType>(currRegion.area());
                if (order != nullptr) {
                    numberOfCopiesOrder[order]--;
                    numberOfCopiesMonRes[localMonotonicityResult]--;
                }
                regionQueue.pop();
            } else {
                regionQueue.pop();
                STORM_LOG_INFO("Splitting region " << currRegion << " into " << newRegions.size());
                numberOfSplits++;
                // Add the new regions to the queue
                if (useMonotonicity) {
                    for (auto& r : newRegions) {
                        r.setBoundParent(storm::utility::convertNumber<CoefficientType>(currBound.get()));
                        regionQueue.emplace(r, order, localMonotonicityResult, currBound.get());
                    }
                    if (numberOfCopiesOrder.find(order) != numberOfCopiesOrder.end()) {
                        numberOfCopiesOrder[order] += newRegions.size();
                        numberOfCopiesMonRes[localMonotonicityResult] += newRegions.size();
                    } else {
                        numberOfCopiesOrder[order] = newRegions.size();
                        numberOfCopiesMonRes[localMonotonicityResult] = newRegions.size();
                    }
                } else {
                    for (auto& r : newRegions) {
                        r.setBoundParent(storm::utility::convertNumber<CoefficientType>(currBound.get()));
                        regionQueue.emplace(r, nullptr, nullptr, currBound.get());
                    }
                }
            }
 
            STORM_LOG_INFO("Covered " << (coveredArea * storm::utility::convertNumber<ConstantType>(100.0) / totalArea) << "% of the region.\n");
            STORM_LOG_INFO("Best value: " << value.get() << ". Regions queued: " << regionQueue.size() << "\n");
        }
        loopWatch.stop();
    }
 
    STORM_LOG_INFO("Total number of splits: " << numberOfSplits << '\n');
    STORM_LOG_INFO("Total number of pla calls: " << numberOfPLACalls << '\n');
    if (useMonotonicity) {
        STORM_LOG_INFO("Total number of pla calls for bounds for monotonicity checking: " << numberOfPLACallsBounds << '\n');
        STORM_LOG_INFO("Total number of copies of the order: " << numberOfOrderCopies << '\n');
        STORM_LOG_INFO("Total number of copies of the local monotonicity result: " << numberOfMonResCopies << '\n');
    }
    STORM_LOG_INFO("\nTotal time for region refinement: " << loopWatch << ".\n\n");
    STORM_LOG_INFO("\nTotal time for additional bounds: " << boundsWatch << ".\n\n");
 
    return std::make_pair(storm::utility::convertNumber<ConstantType>(value.get()), valuation);
}
 
template<typename SparseModelType, typename ConstantType>
std::pair<typename SparseModelType::ValueType, typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation>
SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::computeExtremalValue(
    Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region, storm::solver::OptimizationDirection const& dir,
    typename SparseModelType::ValueType const& precision, bool absolutePrecision, std::optional<storm::logic::Bound> const& terminationBound) {
    auto res = computeExtremalValue(env, region, dir, precision, absolutePrecision, boost::none, terminationBound);
    return {storm::utility::convertNumber<typename SparseModelType::ValueType>(res.first), res.second};
}
 
template<typename SparseModelType, typename ConstantType>
bool SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::verifyRegion(
    Environment const& env, storm::storage::ParameterRegion<typename SparseModelType::ValueType> const& region, storm::logic::Bound const& bound) {
    // Use the bound from the formula.
    ConstantType valueToCheck = storm::utility::convertNumber<ConstantType>(bound.threshold.evaluateAsDouble());
    // We will try to violate the bound.
    storm::solver::OptimizationDirection dir =
        isLowerBound(bound.comparisonType) ? storm::solver::OptimizationDirection::Minimize : storm::solver::OptimizationDirection::Maximize;
    // We pass the bound as an invariant; as soon as it is obtained, we can stop the search.
    auto res = computeExtremalValue(env, region, dir, storm::utility::zero<typename SparseModelType::ValueType>(), false, boost::none, bound).first;
    STORM_LOG_DEBUG("Reported extremal value " << res);
    // TODO use termination bound instead of initial value?
    return storm::solver::minimize(dir) ? storm::utility::convertNumber<ConstantType>(res) >= valueToCheck
                                        : storm::utility::convertNumber<ConstantType>(res) <= valueToCheck;
}
 
template<typename SparseModelType, typename ConstantType>
SparseModelType const& SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::getConsideredParametricModel() const {
    return *parametricModel;
}
 
template<typename SparseModelType, typename ConstantType>
CheckTask<storm::logic::Formula, ConstantType> const& SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::getCurrentCheckTask() const {
    return *currentCheckTask;
}
 
template<typename SparseModelType, typename ConstantType>
void SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyBoundedUntilFormula(
    const CheckTask<logic::BoundedUntilFormula, ConstantType>& checkTask) {
    STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property.");
}
 
template<typename SparseModelType, typename ConstantType>
void SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyUntilFormula(Environment const& env,
                                                                                            CheckTask<logic::UntilFormula, ConstantType> const& checkTask) {
    STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property.");
}
 
template<typename SparseModelType, typename ConstantType>
void SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyReachabilityProbabilityFormula(
    Environment const& env, CheckTask<logic::EventuallyFormula, ConstantType> const& checkTask) {
    // transform to until formula
    auto untilFormula =
        std::make_shared<storm::logic::UntilFormula const>(storm::logic::Formula::getTrueFormula(), checkTask.getFormula().getSubformula().asSharedPointer());
    specifyUntilFormula(env, currentCheckTask->substituteFormula(*untilFormula));
}
 
template<typename SparseModelType, typename ConstantType>
void SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyReachabilityRewardFormula(
    Environment const& env, CheckTask<logic::EventuallyFormula, ConstantType> const& checkTask) {
    STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property.");
}
 
template<typename SparseModelType, typename ConstantType>
void SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::specifyCumulativeRewardFormula(
    const CheckTask<logic::CumulativeRewardFormula, ConstantType>& checkTask) {
    STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Parameter lifting is not supported for the given property.");
}
 
template<typename SparseModelType, typename ConstantType>
std::shared_ptr<storm::analysis::Order> SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::copyOrder(
    std::shared_ptr<storm::analysis::Order> order) {
    auto res = order->copy();
    if (orderExtender) {
        orderExtender->setUnknownStates(order, res);
        orderExtender->copyMinMax(order, res);
    }
    return res;
}
 
template<typename SparseModelType, typename ConstantType>
std::pair<typename SparseModelType::ValueType, typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation>
SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::checkForPossibleMonotonicity(
    Environment const& env, const storage::ParameterRegion<typename SparseModelType::ValueType>& region, std::set<VariableType>& possibleMonotoneIncrParameters,
    std::set<VariableType>& possibleMonotoneDecrParameters, std::set<VariableType>& possibleNotMonotoneParameters,
    std::set<VariableType> const& consideredVariables, storm::solver::OptimizationDirection const& dir) {
    typedef typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation Valuation;
    typedef typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::CoefficientType CoefficientType;
    ConstantType value = storm::solver::minimize(dir) ? 1 : 0;
    Valuation valuation;
    for (auto& var : consideredVariables) {
        ConstantType previousCenter = -1;
        bool monDecr = true;
        bool monIncr = true;
        auto valuationCenter = region.getCenterPoint();
 
        valuationCenter[var] = region.getLowerBoundary(var);
        // TODO: make cmdline argument or 1/precision
        int numberOfSamples = 50;
        auto stepSize = (region.getUpperBoundary(var) - region.getLowerBoundary(var)) / (numberOfSamples - 1);
 
        while (valuationCenter[var] <= region.getUpperBoundary(var)) {
            // Create valuation
            ConstantType valueCenter = getInstantiationChecker()
                                           .check(env, valuationCenter)
                                           ->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()];
            if (storm::solver::minimize(dir) ? valueCenter <= value : valueCenter >= value) {
                value = valueCenter;
                valuation = valuationCenter;
            }
            // Calculate difference with result for previous valuation
            ConstantType diffCenter = previousCenter - valueCenter;
            assert(previousCenter == -1 || (diffCenter >= -1 && diffCenter <= 1));
            if (previousCenter != -1) {
                assert(previousCenter != -1 && previousCenter != -1);
                monDecr &= diffCenter > 0 && diffCenter > 0 && diffCenter > 0;  // then previous value is larger than the current value from the initial states
                monIncr &= diffCenter < 0 && diffCenter < 0 && diffCenter < 0;
            }
            previousCenter = valueCenter;
            if (!monDecr && !monIncr) {
                break;
            }
            valuationCenter[var] += stepSize;
        }
        if (monIncr) {
            possibleMonotoneParameters.insert(var);
            possibleMonotoneIncrParameters.insert(var);
        } else if (monDecr) {
            possibleMonotoneParameters.insert(var);
            possibleMonotoneDecrParameters.insert(var);
        } else {
            possibleNotMonotoneParameters.insert(var);
        }
    }
    return std::make_pair(storm::utility::convertNumber<typename SparseModelType::ValueType>(value), std::move(valuation));
}
 
template<typename SparseModelType, typename ConstantType>
std::pair<typename SparseModelType::ValueType, typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation>
SparseParameterLiftingModelChecker<SparseModelType, ConstantType>::getGoodInitialPoint(
    const Environment& env, const storage::ParameterRegion<typename SparseModelType::ValueType>& region, const OptimizationDirection& dir,
    std::shared_ptr<storm::analysis::LocalMonotonicityResult<VariableType>> localMonRes) {
    typedef typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::Valuation Valuation;
    typedef typename storm::storage::ParameterRegion<typename SparseModelType::ValueType>::CoefficientType CoefficientType;
    ConstantType value = storm::solver::minimize(dir) ? 1 : 0;
    Valuation valuation;
    std::set<VariableType> monIncr, monDecr, notMon, notMonFirst;
    STORM_LOG_INFO("Number of parameters: " << region.getVariables().size() << '\n';);
 
    if (localMonRes != nullptr) {
        localMonRes->getGlobalMonotonicityResult()->splitBasedOnMonotonicity(region.getVariables(), monIncr, monDecr, notMonFirst);
 
        auto numMon = monIncr.size() + monDecr.size();
        STORM_LOG_INFO("Number of monotone parameters: " << numMon << '\n';);
 
        if (numMon < region.getVariables().size()) {
            checkForPossibleMonotonicity(env, region, monIncr, monDecr, notMon, notMonFirst, dir);
            STORM_LOG_INFO("Number of possible monotone parameters: " << (monIncr.size() + monDecr.size() - numMon) << '\n';);
            STORM_LOG_INFO("Number of definitely not monotone parameters: " << notMon.size() << '\n';);
        }
 
        valuation = region.getPoint(dir, monIncr, monDecr);
    } else {
        valuation = region.getCenterPoint();
    }
    value = getInstantiationChecker()
                .check(env, valuation)
                ->template asExplicitQuantitativeCheckResult<ConstantType>()[*this->parametricModel->getInitialStates().begin()];
 
    return std::make_pair(storm::utility::convertNumber<typename SparseModelType::ValueType>(value), std::move(valuation));
}
 
template class SparseParameterLiftingModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>, double>;
template class SparseParameterLiftingModelChecker<storm::models::sparse::Mdp<storm::RationalFunction>, double>;
template class SparseParameterLiftingModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>, storm::RationalNumber>;
template class SparseParameterLiftingModelChecker<storm::models::sparse::Mdp<storm::RationalFunction>, storm::RationalNumber>;
}  // namespace modelchecker
}  // namespace storm