#include "LraViHelper.h"

#include "storm/modelchecker/helper/infinitehorizon/internal/ComponentUtility.h"


#include "storm/storage/MaximalEndComponent.h"
#include "storm/storage/StronglyConnectedComponent.h"

#include "storm/utility/graph.h"
#include "storm/utility/vector.h"
#include "storm/utility/macros.h"
#include "storm/utility/SignalHandler.h"

#include "storm/environment/solver/SolverEnvironment.h"
#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
#include "storm/environment/solver/MinMaxSolverEnvironment.h"


#include "storm/exceptions/UnmetRequirementException.h"


namespace storm {
    namespace modelchecker {
        namespace helper {
            namespace internal {
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                LraViHelper<ValueType, ComponentType, TransitionsType>::LraViHelper(ComponentType const& component, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, ValueType const& aperiodicFactor, storm::storage::BitVector const* timedStates, std::vector<ValueType> const* exitRates) : _component(component), _transitionMatrix(transitionMatrix), _timedStates(timedStates), _hasInstantStates(TransitionsType == LraViTransitionsType::DetTsNondetIs || TransitionsType == LraViTransitionsType::DetTsDetIs) {
                    // Run through the component and collect some data:
                    // We create two submodels, one consisting of the timed states of the component and one consisting of the instant states of the component.
                    // For this, we create a state index map that point from state indices of the input model to indices of the corresponding submodel of that state.
                    boost::container::flat_map<uint64_t, uint64_t> toSubModelStateMapping;
                    // We also obtain state and choices counts of the two submodels
                    uint64_t numTsSubModelStates(0), numTsSubModelChoices(0);
                    uint64_t numIsSubModelStates(0), numIsSubModelChoices(0);
                    // We will need to uniformize the timed MEC states by introducing a selfloop.
                    // For this, we need to find a uniformization rate which will be a little higher (given by aperiodicFactor) than the maximum rate occurring in the component.
                    _uniformizationRate = exitRates == nullptr ? storm::utility::one<ValueType>() : storm::utility::zero<ValueType>();
                    // Now run over the MEC and collect the required data.
                    for (auto const& element : _component) {
                        uint64_t const& componentState = getComponentElementState(element);
                        if (isTimedState(componentState)) {
                            toSubModelStateMapping.emplace(componentState, numTsSubModelStates);
                            ++numTsSubModelStates;
                            numTsSubModelChoices += getComponentElementChoiceCount(element);
                            STORM_LOG_ASSERT(nondetTs() || getComponentElementChoiceCount(element) == 1, "Timed state has multiple choices but only a single choice was expected.");
                            if (exitRates) {
                                _uniformizationRate = std::max(_uniformizationRate, (*exitRates)[componentState]);
                            }
                        } else {
                            toSubModelStateMapping.emplace(componentState, numIsSubModelStates);
                            ++numIsSubModelStates;
                            numIsSubModelChoices += getComponentElementChoiceCount(element);
                            STORM_LOG_ASSERT(nondetIs() || getComponentElementChoiceCount(element) == 1, "Instant state has multiple choices but only a single choice was expected.");
                        }
                    }
                    assert(numIsSubModelStates + numTsSubModelStates == _component.size());
                    assert(_hasInstantStates || numIsSubModelStates == 0);
                    STORM_LOG_ASSERT(nondetTs() || numTsSubModelStates == numTsSubModelChoices, "Unexpected choice count of deterministic timed submodel.");
                    STORM_LOG_ASSERT(nondetIs() || numIsSubModelStates == numIsSubModelChoices, "Unexpected choice count of deterministic instant submodel.");
                    _hasInstantStates = _hasInstantStates && numIsSubModelStates > 0;
                    STORM_LOG_THROW(numTsSubModelStates > 0, storm::exceptions::InvalidOperationException, "Bottom Component has no timed states. Computation of Long Run Average values not supported. Is this a Markov Automaton with Zeno behavior?");

                    // We make sure that every timed state gets a selfloop to make the model aperiodic
                    _uniformizationRate *= storm::utility::one<ValueType>() + aperiodicFactor;

                    // Now build the timed and the instant submodels.
                    // In addition, we also need the transitions between the two.
                    storm::storage::SparseMatrixBuilder<ValueType> tsTransitionsBuilder(numTsSubModelChoices, numTsSubModelStates, 0, true, nondetTs(), nondetTs() ? numTsSubModelStates : 0);
                    storm::storage::SparseMatrixBuilder<ValueType> tsToIsTransitionsBuilder, isTransitionsBuilder, isToTsTransitionsBuilder;
                    if (_hasInstantStates) {
                        tsToIsTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numTsSubModelChoices, numIsSubModelStates, 0, true, nondetTs(), nondetTs() ? numTsSubModelStates : 0);
                        isTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numIsSubModelChoices, numIsSubModelStates, 0, true, nondetIs(), nondetIs() ? numIsSubModelStates : 0);
                        isToTsTransitionsBuilder = storm::storage::SparseMatrixBuilder<ValueType>(numIsSubModelChoices, numTsSubModelStates, 0, true, nondetIs(), nondetIs() ? numIsSubModelStates : 0);
                        _IsChoiceValues.reserve(numIsSubModelChoices);
                    }
                    ValueType uniformizationFactor = storm::utility::one<ValueType>() / _uniformizationRate;
                    uint64_t currTsRow = 0;
                    uint64_t currIsRow = 0;
                    for (auto const& element : _component) {
                        uint64_t const& componentState = getComponentElementState(element);
                        if (isTimedState(componentState)) {
                            // The currently processed state is timed.
                            if (nondetTs()) {
                                tsTransitionsBuilder.newRowGroup(currTsRow);
                                if (_hasInstantStates) {
                                    tsToIsTransitionsBuilder.newRowGroup(currTsRow);
                                }
                            }
                            // We need to uniformize which means that a diagonal entry for the selfloop will be inserted.
                            // If there are exit rates, the uniformization factor needs to be updated.
                            if (exitRates) {
                                uniformizationFactor = (*exitRates)[componentState] / _uniformizationRate;
                            }
                            ValueType selfLoopProb = storm::utility::one<ValueType>() - uniformizationFactor;
                            uint64_t selfLoopColumn = toSubModelStateMapping[componentState];
                            for (auto componentChoiceIt = getComponentElementChoicesBegin(element); componentChoiceIt != getComponentElementChoicesEnd(element); ++componentChoiceIt) {
                                bool insertedDiagElement = false;
                                for (auto const& entry : this->_transitionMatrix.getRow(*componentChoiceIt)) {
                                    uint64_t subModelColumn = toSubModelStateMapping[entry.getColumn()];
                                    if (isTimedState(entry.getColumn())) {
                                        // We have a transition from a timed state to a timed state
                                        STORM_LOG_ASSERT(subModelColumn < numTsSubModelStates, "Invalid state for timed submodel");
                                        if (!insertedDiagElement && subModelColumn > selfLoopColumn) {
                                            // We passed the diagonal entry, so add it now before moving on to the next entry
                                            tsTransitionsBuilder.addNextValue(currTsRow, selfLoopColumn, selfLoopProb);
                                            insertedDiagElement = true;
                                        }
                                        if (!insertedDiagElement && subModelColumn == selfLoopColumn) {
                                            // The current entry is the diagonal (selfloop) entry
                                            tsTransitionsBuilder.addNextValue(currTsRow, selfLoopColumn, selfLoopProb + uniformizationFactor * entry.getValue());
                                            insertedDiagElement = true;
                                        } else {
                                            // The diagonal element either has been inserted already or still lies in front
                                            tsTransitionsBuilder.addNextValue(currTsRow, subModelColumn,  uniformizationFactor * entry.getValue());
                                        }
                                    } else {
                                        // We have a transition from a timed to a instant state
                                        STORM_LOG_ASSERT(subModelColumn < numIsSubModelStates, "Invalid state for instant submodel");
                                        tsToIsTransitionsBuilder.addNextValue(currTsRow, subModelColumn, uniformizationFactor * entry.getValue());
                                    }
                                }
                                // If the diagonal entry for the MS matrix still has not been set, we do that now
                                if (!insertedDiagElement) {
                                    tsTransitionsBuilder.addNextValue(currTsRow, selfLoopColumn, selfLoopProb);
                                }
                                ++currTsRow;
                            }
                        } else {
                            // The currently processed state is instant
                            if (nondetIs()) {
                                isTransitionsBuilder.newRowGroup(currIsRow);
                                isToTsTransitionsBuilder.newRowGroup(currIsRow);
                            }
                            for (auto componentChoiceIt = getComponentElementChoicesBegin(element); componentChoiceIt != getComponentElementChoicesEnd(element); ++componentChoiceIt) {
                                for (auto const& entry : this->_transitionMatrix.getRow(*componentChoiceIt)) {
                                    uint64_t subModelColumn = toSubModelStateMapping[entry.getColumn()];
                                    if (isTimedState(entry.getColumn())) {
                                        // We have a transition from an instant state to a timed state
                                        STORM_LOG_ASSERT(subModelColumn < numTsSubModelStates, "Invalid state for timed submodel");
                                        isToTsTransitionsBuilder.addNextValue(currIsRow, subModelColumn, entry.getValue());
                                    } else {
                                        // We have a transition from an instant to an instant state
                                        STORM_LOG_ASSERT(subModelColumn < numIsSubModelStates, "Invalid state for instant submodel");
                                        isTransitionsBuilder.addNextValue(currIsRow, subModelColumn, entry.getValue());
                                    }
                                }
                                ++currIsRow;
                            }
                        }
                    }
                    _TsTransitions = tsTransitionsBuilder.build();
                    if (_hasInstantStates) {
                        _TsToIsTransitions = tsToIsTransitionsBuilder.build();
                        _IsTransitions = isTransitionsBuilder.build();
                        _IsToTsTransitions = isToTsTransitionsBuilder.build();
                    }
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                ValueType LraViHelper<ValueType, ComponentType, TransitionsType>::performValueIteration(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector<ValueType> const* exitRates, storm::solver::OptimizationDirection const* dir, std::vector<uint64_t>* choices) {
                    
                    initializeNewValues(stateValueGetter, actionValueGetter, exitRates);
                    
                    ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().lra().getPrecision());
                    bool relative = env.solver().lra().getRelativeTerminationCriterion();
                    boost::optional<uint64_t> maxIter;
                    if (env.solver().lra().isMaximalIterationCountSet()) {
                        maxIter = env.solver().lra().getMaximalIterationCount();
                    }
                
                    // start the iterations
                    ValueType result = storm::utility::zero<ValueType>();
                    uint64_t iter = 0;
                    while (!maxIter.is_initialized() || iter < maxIter.get()) {
                        ++iter;
                        performIterationStep(env, dir);
                        
                        // Check if we are done
                        auto convergenceCheckResult = checkConvergence(relative, precision);
                        result = convergenceCheckResult.currentValue;
                        if (convergenceCheckResult.isPrecisionAchieved) {
                            break;
                        }
                        if (storm::utility::resources::isTerminate()) {
                            break;
                        }
                        // If there will be a next iteration, we have to prepare it.
                        prepareNextIteration(env);
                        
                    }
                    if (maxIter.is_initialized() && iter == maxIter.get()) {
                        STORM_LOG_WARN("LRA computation did not converge within " << iter << " iterations.");
                    } else if (storm::utility::resources::isTerminate()) {
                        STORM_LOG_WARN("LRA computation aborted after " << iter << " iterations.");
                    } else {
                        STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations.");
                    }
                    
                    if (choices) {
                        // We will be doing one more iteration step and track scheduler choices this time.
                        prepareNextIteration(env);
                        performIterationStep(env, dir, choices);
                    }
                    return result;
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                void LraViHelper<ValueType, ComponentType, TransitionsType>::initializeNewValues(ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, std::vector<ValueType> const* exitRates) {
                    // clear potential old values and reserve enough space for new values
                    _TsChoiceValues.clear();
                    _TsChoiceValues.reserve(_TsTransitions.getRowCount());
                    if (_hasInstantStates) {
                        _IsChoiceValues.clear();
                        _IsChoiceValues.reserve(_IsTransitions.getRowCount());
                    }
                    
                    // Set the new choice-based values
                    ValueType actionRewardScalingFactor = storm::utility::one<ValueType>() / _uniformizationRate;
                    for (auto const& element : _component) {
                        uint64_t const& componentState = getComponentElementState(element);
                        if (isTimedState(componentState)) {
                            if (exitRates) {
                                actionRewardScalingFactor = (*exitRates)[componentState] / _uniformizationRate;
                            }
                            for (auto componentChoiceIt = getComponentElementChoicesBegin(element); componentChoiceIt != getComponentElementChoicesEnd(element); ++componentChoiceIt) {
                                // Compute the values obtained for this choice.
                                _TsChoiceValues.push_back(stateValueGetter(componentState) / _uniformizationRate + actionValueGetter(*componentChoiceIt) * actionRewardScalingFactor);
                            }
                        } else {
                            for (auto componentChoiceIt = getComponentElementChoicesBegin(element); componentChoiceIt != getComponentElementChoicesEnd(element); ++componentChoiceIt) {
                                // Compute the values obtained for this choice.
                                // State values do not count here since no time passes in instant states.
                                _IsChoiceValues.push_back(actionValueGetter(*componentChoiceIt));
                            }
                        }
                    }

                    // Set-up new iteration vectors for timed states
                    _Tsx1.assign(_TsTransitions.getRowGroupCount(), storm::utility::zero<ValueType>());
                    _Tsx2 = _Tsx1;
                    
                    if (_hasInstantStates) {
                        // Set-up vectors for storing intermediate results for instant states.
                        _Isx.resize(_IsTransitions.getRowGroupCount(), storm::utility::zero<ValueType>());
                        _Isb = _IsChoiceValues;
                    }
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                void LraViHelper<ValueType, ComponentType, TransitionsType>::prepareSolversAndMultipliers(const Environment& env, storm::solver::OptimizationDirection const* dir) {
                    _TsMultiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _TsTransitions);
                    if (_hasInstantStates) {
                        if (_IsTransitions.getNonzeroEntryCount() > 0) {
                            // Set-up a solver for transitions within instant states
                            _IsSolverEnv = std::make_unique<storm::Environment>(env);
                            if (env.solver().isForceSoundness()) {
                                // To get correct results, the inner equation systems are solved exactly.
                                // TODO investigate how an error would propagate
                                _IsSolverEnv->solver().setForceExact(true);
                            }
                            bool isAcyclic = !storm::utility::graph::hasCycle(_IsTransitions);
                            if (isAcyclic) {
                                STORM_LOG_INFO("Instant transitions are acyclic.");
                                if (_IsSolverEnv->solver().minMax().isMethodSetFromDefault()) {
                                    _IsSolverEnv->solver().minMax().setMethod(storm::solver::MinMaxMethod::Acyclic);
                                }
                                if (_IsSolverEnv->solver().isLinearEquationSolverTypeSetFromDefaultValue()) {
                                   _IsSolverEnv->solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Acyclic);
                                }
                            }
                            if (nondetIs()) {
                                storm::solver::GeneralMinMaxLinearEquationSolverFactory<ValueType> factory;
                                _NondetIsSolver = factory.create(*_IsSolverEnv, _IsTransitions);
                                _NondetIsSolver->setHasUniqueSolution(true); // Assume non-zeno MA
                                _NondetIsSolver->setHasNoEndComponents(true); // assume non-zeno MA
                                _NondetIsSolver->setCachingEnabled(true);
                                auto req = _NondetIsSolver->getRequirements(*_IsSolverEnv, *dir);
                                req.clearUniqueSolution();
                                if (isAcyclic) {
                                    req.clearAcyclic();
                                }
                                // Computing a priori lower/upper bounds is not particularly easy, as there might be selfloops with high probabilities
                                // Which accumulate a lot of reward. Moreover, the right-hand-side of the equation system changes dynamically.
                                STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "The solver requirement " << req.getEnabledRequirementsAsString() << " has not been cleared.");
                                _NondetIsSolver->setRequirementsChecked(true);
                            } else {
                                storm::solver::GeneralLinearEquationSolverFactory<ValueType> factory;
                                if (factory.getEquationProblemFormat(*_IsSolverEnv) != storm::solver::LinearEquationSolverProblemFormat::FixedPointSystem) {
                                    // We need to convert the transition matrix connecting instant states
                                    // TODO: This could have been done already during construction of the matrix.
                                    // Insert diagonal entries.
                                    storm::storage::SparseMatrix<ValueType> converted(_IsTransitions, true);
                                    // Compute A' = 1-A
                                    converted.convertToEquationSystem();
                                    STORM_LOG_WARN("The selected equation solver requires to create a temporary " << converted.getDimensionsAsString());
                                    // Note that the solver has ownership of the converted matrix.
                                    _DetIsSolver = factory.create(*_IsSolverEnv, std::move(converted));
                                } else {
                                    _DetIsSolver = factory.create(*_IsSolverEnv, _IsTransitions);
                                }
                                _DetIsSolver->setCachingEnabled(true);
                                auto req = _DetIsSolver->getRequirements(*_IsSolverEnv);
                                if (isAcyclic) {
                                    req.clearAcyclic();
                                }
                                // A priori lower/upper bounds are hard (see MinMax version of this above)
                                STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "The solver requirement " << req.getEnabledRequirementsAsString() << " has not been cleared.");
                            }
                        }
                        
                        // Set up multipliers for transitions connecting timed and instant states
                        _TsToIsMultiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _TsToIsTransitions);
                        _IsToTsMultiplier = storm::solver::MultiplierFactory<ValueType>().create(env, _IsToTsTransitions);
                    }
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                void LraViHelper<ValueType, ComponentType, TransitionsType>::setInputModelChoices(std::vector<uint64_t>& choices, std::vector<uint64_t> const& localMecChoices, bool setChoiceZeroToTimedStates, bool setChoiceZeroToInstantStates) const {
                    // Transform the local choices (within this mec) to choice indices for the input model
                    uint64_t localState = 0;
                    for (auto const& element : _component) {
                        uint64_t elementState = getComponentElementState(element);
                        if ((setChoiceZeroToTimedStates && isTimedState(elementState)) || (setChoiceZeroToInstantStates && !isTimedState(elementState))) {
                            choices[elementState] = 0;
                        } else {
                            uint64_t choice = localMecChoices[localState];
                            STORM_LOG_ASSERT(choice < getComponentElementChoiceCount(element), "The selected choice does not seem to exist.");
                            uint64_t globalChoiceIndex = *(getComponentElementChoicesBegin(element) + choice);
                            choices[elementState] = globalChoiceIndex - _transitionMatrix.getRowGroupIndices()[elementState];
                            ++localState;
                        }
                    }
                    STORM_LOG_ASSERT(localState == localMecChoices.size(), "Did not traverse all component states.");
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                void LraViHelper<ValueType, ComponentType, TransitionsType>::performIterationStep(Environment const& env, storm::solver::OptimizationDirection const* dir, std::vector<uint64_t>* choices) {
                    STORM_LOG_ASSERT(!((nondetTs() || nondetIs()) && dir == nullptr), "No optimization direction provided for model with nondeterminism");
                    // Initialize value vectors, multiplers, and solver if this has not been done, yet
                    if (!_TsMultiplier) {
                        prepareSolversAndMultipliers(env, dir);
                    }
                    
                    // Compute new x values for the timed states
                    // Flip what is new and what is old
                    _Tsx1IsCurrent = !_Tsx1IsCurrent;
                    // At this point, xOld() points to what has been computed in the most recent call of performIterationStep (initially, this is the 0-vector).
                    // The result of this ongoing computation will be stored in xNew()
                    
                    // Compute the values obtained by a single uniformization step between timed states only
                    if (nondetTs()) {
                        if (choices == nullptr) {
                            _TsMultiplier->multiplyAndReduce(env, *dir, xOld(), &_TsChoiceValues, xNew());
                        } else {
                            // Also keep track of the choices made.
                            std::vector<uint64_t> tsChoices(_TsTransitions.getRowGroupCount());
                            _TsMultiplier->multiplyAndReduce(env, *dir, xOld(), &_TsChoiceValues, xNew(), &tsChoices);
                            // Note that nondeterminism within the timed states means that there can not be instant states (We either have MDPs or MAs)
                            // Hence, in this branch we don't have to care for choices at instant states.
                            STORM_LOG_ASSERT(!_hasInstantStates, "Nondeterministic timed states are only supported if there are no instant states.");
                            setInputModelChoices(*choices, tsChoices);
                        }
                    } else {
                        _TsMultiplier->multiply(env, xOld(), &_TsChoiceValues, xNew());
                    }
                    if (_hasInstantStates) {
                        // Add the values obtained by taking a single uniformization step that leads to an instant state followed by arbitrarily many instant steps.
                        // First compute the total values when taking arbitrarily many instant transitions (in no time)
                        if (_NondetIsSolver) {
                            // We might need to track the optimal choices.
                            if (choices == nullptr) {
                                _NondetIsSolver->solveEquations(*_IsSolverEnv, *dir, _Isx, _Isb);
                            } else {
                                _NondetIsSolver->setTrackScheduler();
                                _NondetIsSolver->solveEquations(*_IsSolverEnv, *dir, _Isx, _Isb);
                                setInputModelChoices(*choices, _NondetIsSolver->getSchedulerChoices(), true);
                            }
                        } else if (_DetIsSolver) {
                            _DetIsSolver->solveEquations(*_IsSolverEnv, _Isx, _Isb);
                        } else {
                            STORM_LOG_ASSERT(_IsTransitions.getNonzeroEntryCount() == 0, "If no solver was initialized, an empty matrix would have been expected.");
                            if (nondetIs()) {
                                if (choices == nullptr) {
                                    storm::utility::vector::reduceVectorMinOrMax(*dir, _Isb, _Isx, _IsTransitions.getRowGroupIndices());
                                } else {
                                    std::vector<uint64_t> psChoices(_IsTransitions.getRowGroupCount());
                                    storm::utility::vector::reduceVectorMinOrMax(*dir, _Isb, _Isx, _IsTransitions.getRowGroupIndices(), &psChoices);
                                    setInputModelChoices(*choices, psChoices, true);
                                }
                            } else {
                                // For deterministic instant states, there is nothing to reduce, i.e., we could just set _Isx = _Isb.
                                // For efficiency reasons, we do a swap instead:
                                _Isx.swap(_Isb);
                                // Note that at this point we have changed the contents of _Isb, but they will be overwritten anyway.
                                if (choices) {
                                    // Set choice 0 to all states.
                                    setInputModelChoices(*choices, {}, true, true);
                                }
                            }
                        }
                        // Now add the (weighted) values of the instant states to the values of the timed states.
                        _TsToIsMultiplier->multiply(env, _Isx, &xNew(), xNew());
                    }
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                typename LraViHelper<ValueType, ComponentType, TransitionsType>::ConvergenceCheckResult LraViHelper<ValueType, ComponentType, TransitionsType>::checkConvergence(bool relative, ValueType precision) const {
                    STORM_LOG_ASSERT(_TsMultiplier, "tried to check for convergence without doing an iteration first.");
                    // All values are scaled according to the uniformizationRate.
                    // We need to 'revert' this scaling when computing the absolute precision.
                    // However, for relative precision, the scaling cancels out.
                    ValueType threshold = relative ? precision : ValueType(precision / _uniformizationRate);
                    
                    ConvergenceCheckResult res = { true, storm::utility::one<ValueType>() };
                    // Now check whether the currently produced results are precise enough
                    STORM_LOG_ASSERT(threshold > storm::utility::zero<ValueType>(), "Did not expect a non-positive threshold.");
                    auto x1It = xOld().begin();
                    auto x1Ite = xOld().end();
                    auto x2It = xNew().begin();
                    ValueType maxDiff = (*x2It - *x1It);
                    ValueType minDiff = maxDiff;
                    // The difference between maxDiff and minDiff is zero at this point. Thus, it doesn't make sense to check the threshold now.
                    for (++x1It, ++x2It; x1It != x1Ite; ++x1It, ++x2It) {
                        ValueType diff = (*x2It - *x1It);
                        // Potentially update maxDiff or minDiff
                        bool skipCheck = false;
                        if (maxDiff < diff) {
                            maxDiff = diff;
                        } else if (minDiff > diff) {
                            minDiff = diff;
                        } else {
                            skipCheck = true;
                        }
                        // Check convergence
                        if (!skipCheck && (maxDiff - minDiff) > (relative ? (threshold * minDiff) : threshold)) {
                            res.isPrecisionAchieved = false;
                            break;
                        }
                    }
                    
                    // Compute the average of the maximal and the minimal difference.
                    ValueType avgDiff = (maxDiff + minDiff) / (storm::utility::convertNumber<ValueType>(2.0));
                    
                    // "Undo" the scaling of the values
                    res.currentValue = avgDiff * _uniformizationRate;
                    return res;
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                void LraViHelper<ValueType, ComponentType, TransitionsType>::prepareNextIteration(Environment const& env) {
                    // To avoid large (and numerically unstable) x-values, we substract a reference value.
                    ValueType referenceValue = xNew().front();
                    storm::utility::vector::applyPointwise<ValueType, ValueType>(xNew(), xNew(), [&referenceValue] (ValueType const& x_i) -> ValueType { return x_i - referenceValue; });
                    if (_hasInstantStates) {
                        // Update the RHS of the equation system for the instant states by taking the new values of timed states into account.
                        STORM_LOG_ASSERT(!nondetTs(), "Nondeterministic timed states not expected when there are also instant states.");
                        _IsToTsMultiplier->multiply(env, xNew(), &_IsChoiceValues, _Isb);
                    }
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                bool LraViHelper<ValueType, ComponentType, TransitionsType>::isTimedState(uint64_t const& inputModelStateIndex) const {
                    STORM_LOG_ASSERT(!_hasInstantStates || inputModelStateIndex < _timedStates->size(), "Unable to determine whether state " << inputModelStateIndex << " is timed.");
                    return !_hasInstantStates || _timedStates->get(inputModelStateIndex);
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                std::vector<ValueType>& LraViHelper<ValueType, ComponentType, TransitionsType>::xNew() {
                    return _Tsx1IsCurrent ? _Tsx1 : _Tsx2;
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                std::vector<ValueType> const& LraViHelper<ValueType, ComponentType, TransitionsType>::xNew() const {
                    return _Tsx1IsCurrent ? _Tsx1 : _Tsx2;
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                std::vector<ValueType>& LraViHelper<ValueType, ComponentType, TransitionsType>::xOld() {
                    return _Tsx1IsCurrent ? _Tsx2 : _Tsx1;
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                std::vector<ValueType> const& LraViHelper<ValueType, ComponentType, TransitionsType>::xOld() const {
                    return _Tsx1IsCurrent ? _Tsx2 : _Tsx1;
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                bool LraViHelper<ValueType, ComponentType, TransitionsType>::nondetTs() const {
                    return TransitionsType == LraViTransitionsType::NondetTsNoIs;
                }
                
                template <typename ValueType, typename ComponentType, LraViTransitionsType TransitionsType>
                bool LraViHelper<ValueType, ComponentType, TransitionsType>::nondetIs() const {
                    return TransitionsType == LraViTransitionsType::DetTsNondetIs;
                }
                
                template class LraViHelper<double, storm::storage::MaximalEndComponent, LraViTransitionsType::NondetTsNoIs>;
                template class LraViHelper<storm::RationalNumber, storm::storage::MaximalEndComponent, LraViTransitionsType::NondetTsNoIs>;
                template class LraViHelper<double, storm::storage::MaximalEndComponent, LraViTransitionsType::DetTsNondetIs>;
                template class LraViHelper<storm::RationalNumber, storm::storage::MaximalEndComponent, LraViTransitionsType::DetTsNondetIs>;
                
            }
        }
    }
}