StrongBisimulationTest.cpp

Go to the documentation of this file.

#include "BisimulationTestHelper.h"
 
#include <algorithm>
#include <cstdint>
#include <numeric>
#include <random>
#include <set>
 
#include "storm/adapters/RationalNumberAdapter.h"
#include "storm/environment/Environment.h"
#include "storm/exceptions/UnexpectedException.h"
#include "storm/transformer/StatePermuter.h"
 
namespace {
 
using storm::test::bisimulation::buildFromPrism;
using storm::test::bisimulation::buildMarkovAutomaton;
using storm::test::bisimulation::buildModel;
using storm::test::bisimulation::BuildOptions;
using storm::test::bisimulation::checkFormula;
using storm::test::bisimulation::Options;
using storm::test::bisimulation::StateLabelPreservation;
using storm::test::bisimulation::strongOptions;
 
using ValueType = double;

void checkChoiceMapping(storm::models::sparse::Model<ValueType> const& model, storm::bisimulation::ReturnType<ValueType> const& result) {
    ASSERT_TRUE(result.toQuotientChoiceMapping.has_value());
    ASSERT_EQ(model.getNumberOfChoices(), result.toQuotientChoiceMapping->size());
    for (uint64_t state = 0; state < model.getNumberOfStates(); ++state) {
        auto const quotientChoices = result.quotient->getTransitionMatrix().getRowGroupIndices(result.toQuotientStateMapping[state]);
        std::set<uint64_t> const expectedChoices(quotientChoices.begin(), quotientChoices.end());
        std::set<uint64_t> mappedChoices;
        for (uint64_t const choice : model.getTransitionMatrix().getRowGroupIndices(state)) {
            mappedChoices.insert((*result.toQuotientChoiceMapping)[choice]);
        }
        EXPECT_EQ(expectedChoices, mappedChoices) << "unexpected quotient choices for the choices of state " << state;
    }
}

void testQuotient(std::string const& prismFile, std::string const& formulaString, uint64_t expectedModelStates, uint64_t expectedStates,
                  uint64_t expectedTransitions, uint64_t expectedChoices, Options const options = strongOptions(), BuildOptions buildOptions = {}) {
#ifndef STORM_HAVE_Z3
    GTEST_SKIP() << "Z3 not available.";
#endif
    if (options.stateLabelPreservation == StateLabelPreservation::All) {
        buildOptions.allLabels = true;
    }
    auto const input = buildFromPrism<ValueType>(prismFile, formulaString, buildOptions);
    ASSERT_EQ(expectedModelStates, input.model->getNumberOfStates());
 
    Options quotientOptions = options;
    quotientOptions.createQuotientChoiceMapping = true;  // Does not affect the quotient, but lets us check the mapping.
    auto const result = storm::bisimulation::performBisimulationMinimization<ValueType>(*input.model, input.formulas, quotientOptions);
    auto const& quotient = result.quotient;
    checkChoiceMapping(*input.model, result);
    EXPECT_EQ(input.model->getType(), quotient->getType());
    EXPECT_EQ(expectedStates, quotient->getNumberOfStates());
    EXPECT_EQ(expectedTransitions, quotient->getNumberOfTransitions());
    EXPECT_EQ(expectedChoices, quotient->getNumberOfChoices());
    // The quotient keeps the choice labels and the choice origins of the model (but not necessarily an empty choice labeling).
    auto const choiceLabels = [](auto const& model) {
        std::set<std::string> result;
        if (model->hasChoiceLabeling()) {
            auto const labels = model->getChoiceLabeling().getLabels();
            result.insert(labels.begin(), labels.end());
        }
        return result;
    };
    EXPECT_EQ(choiceLabels(input.model), choiceLabels(quotient));
    EXPECT_EQ(input.model->hasChoiceOrigins(), quotient->hasChoiceOrigins());
 
    // The formulas were parsed for the PRISM program, so unlike the formula string they can be checked even if they contain atomic expressions.
    for (auto const& formula : input.formulas) {
        EXPECT_NEAR(checkFormula<ValueType>(quotient, formula), checkFormula<ValueType>(input.model, formula), 1e-9) << *formula;
    }
}

Options allLabelOptions() {
    Options options = strongOptions();
    options.stateLabelPreservation = StateLabelPreservation::All;
    return options;
}

Options approximateOptions() {
    Options options = strongOptions();
    options.tolerance = storm::utility::convertNumber<storm::RationalNumber>(1e-9);
    return options;
}
 
// ------------------------------------------------------------
// Deterministic models
// ------------------------------------------------------------

TEST(StrongBisimulationTest, SilentChainIsNotCollapsed) {
    // 0 -> 1 -> 2 -> {3, 4}; 3 and 4 absorbing, 3 labeled "goal".
    storm::storage::SparseMatrixBuilder<ValueType> builder(5, 5);
    builder.addNextValue(0, 1, 1.0);
    builder.addNextValue(1, 2, 1.0);
    builder.addNextValue(2, 3, 0.5);
    builder.addNextValue(2, 4, 0.5);
    builder.addNextValue(3, 3, 1.0);
    builder.addNextValue(4, 4, 1.0);
    auto const model = buildModel<storm::models::sparse::Dtmc<ValueType>>(builder.build(), {{"goal", {3}}});
 
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient;
    EXPECT_EQ(5ull, quotient->getNumberOfStates());
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "P=? [F \"goal\"]"), 0.5, 1e-12);
}

TEST(StrongBisimulationTest, IdenticalDistributions) {
    // 1 and 2 have the same distribution over {3} and {4}, so they are merged; 0 keeps them apart only through its own distribution.
    storm::storage::SparseMatrixBuilder<ValueType> builder(5, 5);
    builder.addNextValue(0, 1, 0.5);
    builder.addNextValue(0, 2, 0.5);
    builder.addNextValue(1, 3, 0.25);
    builder.addNextValue(1, 4, 0.75);
    builder.addNextValue(2, 3, 0.25);
    builder.addNextValue(2, 4, 0.75);
    builder.addNextValue(3, 3, 1.0);
    builder.addNextValue(4, 4, 1.0);
    auto const model = buildModel<storm::models::sparse::Dtmc<ValueType>>(builder.build(), {{"goal", {3}}});
 
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient;
    EXPECT_EQ(4ull, quotient->getNumberOfStates());  // {0}, {1,2}, {3}, {4}
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "P=? [F \"goal\"]"), 0.25, 1e-12);
}

TEST(StrongBisimulationTest, CtmcInternalRateIsObservable) {
    // 0 -> 1 with rate 100 and 0 -> 2 with rate 1; 1 -> 2 with rate 1; 2 absorbing and labeled "goal".
    storm::storage::SparseMatrixBuilder<ValueType> builder(3, 3);
    builder.addNextValue(0, 1, 100.0);
    builder.addNextValue(0, 2, 1.0);
    builder.addNextValue(1, 2, 1.0);
    builder.addNextValue(2, 2, 1.0);
    auto const model = buildModel<storm::models::sparse::Ctmc<ValueType>>(builder.build(), {{"goal", {2}}});
 
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient;
    EXPECT_EQ(storm::models::ModelType::Ctmc, quotient->getType());
    EXPECT_EQ(3ull, quotient->getNumberOfStates());  // weak bisimulation merges 0 and 1 here, strong bisimulation does not
}

TEST(StrongBisimulationTest, CtmcExitRateIsObservable) {
    // 0 and 1 both move to 2 only, but at different speeds.
    storm::storage::SparseMatrixBuilder<ValueType> builder(3, 3);
    builder.addNextValue(0, 2, 1.0);
    builder.addNextValue(1, 2, 2.0);
    builder.addNextValue(2, 2, 1.0);
    auto const model = buildModel<storm::models::sparse::Ctmc<ValueType>>(builder.build(), {{"goal", {2}}});
 
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient;
    EXPECT_EQ(3ull, quotient->getNumberOfStates());
}

TEST(StrongBisimulationTest, PropositionalSubformulas) {
    // 0 moves uniformly to 1, 2 and 3, which are absorbing. 1 is labeled with "a" and "b", 2 is unlabeled and 3 is labeled with "a".
    storm::storage::SparseMatrixBuilder<ValueType> builder(4, 4);
    builder.addNextValue(0, 1, 1.0 / 3);
    builder.addNextValue(0, 2, 1.0 / 3);
    builder.addNextValue(0, 3, 1.0 / 3);
    builder.addNextValue(1, 1, 1.0);
    builder.addNextValue(2, 2, 1.0);
    builder.addNextValue(3, 3, 1.0);
    auto const model = buildModel<storm::models::sparse::Dtmc<ValueType>>(builder.build(), {{"a", {1, 3}}, {"b", {1}}});
    Options individualOptions = strongOptions();
    individualOptions.stateLabelPreservation = StateLabelPreservation::FormulaIndividual;
 
    storm::parser::FormulaParser formulaParser;
    for (std::string const formulaString : {"P=? [F \"a\" & !\"b\"]", "P=? [F P>=0.5 [F \"a\" & !\"b\"]]"}) {
        std::vector<std::shared_ptr<storm::logic::Formula const>> const formulas{formulaParser.parseSingleFormulaFromString(formulaString)};
        auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, formulas, strongOptions()).quotient;
        EXPECT_EQ(3ull, quotient->getNumberOfStates()) << formulaString;  // {0}, {1, 2} and {3}
        EXPECT_NEAR(checkFormula<ValueType>(quotient, formulaString), 1.0 / 3, 1e-12) << formulaString;
 
        auto const individualQuotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, formulas, individualOptions).quotient;
        EXPECT_EQ(4ull, individualQuotient->getNumberOfStates()) << formulaString;
        EXPECT_NEAR(checkFormula<ValueType>(individualQuotient, formulaString), 1.0 / 3, 1e-12) << formulaString;
    }
}

TEST(StrongBisimulationTest, InitLabelInPropositionalSubformula) {
    // 0 and 1 both move to the absorbing state 2, which is labeled "a". Only 1 is initial, so merging 0 and 1 would make 0 the representative of the initial
    // block and thereby lose the initial state.
    storm::storage::SparseMatrixBuilder<ValueType> builder(3, 3);
    builder.addNextValue(0, 2, 1.0);
    builder.addNextValue(1, 2, 1.0);
    builder.addNextValue(2, 2, 1.0);
    auto const model = buildModel<storm::models::sparse::Dtmc<ValueType>>(builder.build(), {{"a", {2}}});
    model->getStateLabeling().removeLabelFromState("init", 0);
    model->getStateLabeling().addLabelToState("init", 1);
 
    std::string const formulaString = "P=? [F \"a\" & !\"init\"]";
    storm::parser::FormulaParser formulaParser;
    std::vector<std::shared_ptr<storm::logic::Formula const>> const formulas{formulaParser.parseSingleFormulaFromString(formulaString)};
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, formulas, strongOptions()).quotient;
    EXPECT_EQ(3ull, quotient->getNumberOfStates());
    ASSERT_EQ(1ull, quotient->getInitialStates().getNumberOfSetBits());
    EXPECT_NEAR(checkFormula<ValueType>(quotient, formulaString), 1.0, 1e-12);
}

TEST(StrongBisimulationTest, InitLabelCombinedWithOtherLabel) {
    // 0 (labeled "a") and 1 (the only initial state) both move to the absorbing state 2. No state satisfies "init" & "a", but a block {0, 1, 2} would be
    // represented by state 0 and thus yield a quotient state that is both initial and labeled "a".
    storm::storage::SparseMatrixBuilder<ValueType> builder(3, 3);
    builder.addNextValue(0, 2, 1.0);
    builder.addNextValue(1, 2, 1.0);
    builder.addNextValue(2, 2, 1.0);
    auto const model = buildModel<storm::models::sparse::Dtmc<ValueType>>(builder.build(), {{"a", {0}}});
    model->getStateLabeling().removeLabelFromState("init", 0);
    model->getStateLabeling().addLabelToState("init", 1);
 
    std::string const formulaString = "P=? [F \"init\" & \"a\"]";
    storm::parser::FormulaParser formulaParser;
    std::vector<std::shared_ptr<storm::logic::Formula const>> const formulas{formulaParser.parseSingleFormulaFromString(formulaString)};
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, formulas, strongOptions()).quotient;
    EXPECT_EQ(2ull, quotient->getNumberOfStates());  // {1} and {0, 2}
    EXPECT_NEAR(checkFormula<ValueType>(quotient, formulaString), 0.0, 1e-12);
}

TEST(StrongBisimulationTest, FloatingPointRoundingNeedsTolerance) {
    // 1 reaches {3, 4} via two transitions (0.1 and 0.2), 2 via a single one (0.3); both sums equal 0.3 mathematically but not as doubles.
    storm::storage::SparseMatrixBuilder<ValueType> builder(6, 6);
    builder.addNextValue(0, 1, 0.5);
    builder.addNextValue(0, 2, 0.5);
    builder.addNextValue(1, 3, 0.1);
    builder.addNextValue(1, 4, 0.2);
    builder.addNextValue(1, 5, 0.7);
    builder.addNextValue(2, 3, 0.3);
    builder.addNextValue(2, 5, 0.7);
    builder.addNextValue(3, 3, 1.0);
    builder.addNextValue(4, 4, 1.0);
    builder.addNextValue(5, 5, 1.0);
    auto const model = buildModel<storm::models::sparse::Dtmc<ValueType>>(builder.build(), {{"target", {3, 4}}});
 
    EXPECT_EQ(5ull, storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient->getNumberOfStates());
    EXPECT_EQ(4ull, storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, approximateOptions()).quotient->getNumberOfStates());
    // When no tolerance is set, the model's stochastic tolerance is used -- the default model tolerance (see
    // storm::Environment::DefaultModelTolerance), which also identifies the two states.
    EXPECT_EQ(storm::utility::convertNumber<ValueType>(storm::Environment::DefaultModelTolerance), model->getStochasticTolerance());
    EXPECT_EQ(4ull, storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, Options{}).quotient->getNumberOfStates());
}
 
// ------------------------------------------------------------
// Nondeterministic models
// ------------------------------------------------------------

TEST(StrongBisimulationTest, ChoiceSetsAreCompared) {
    // States 0 and 1 offer the same two distributions over {2} and {3}, but in a different order, and 1 offers one of them twice.
    storm::storage::SparseMatrixBuilder<ValueType> builder(7, 4, 0, false, true, 4);
    builder.newRowGroup(0);
    builder.addNextValue(0, 2, 1.0);  // state 0, choice 0
    builder.addNextValue(1, 2, 0.5);  // state 0, choice 1
    builder.addNextValue(1, 3, 0.5);
    builder.newRowGroup(2);
    builder.addNextValue(2, 2, 0.5);  // state 1, choice 0
    builder.addNextValue(2, 3, 0.5);
    builder.addNextValue(3, 2, 1.0);  // state 1, choice 1
    builder.addNextValue(4, 2, 1.0);  // state 1, choice 2 (a duplicate of choice 1)
    builder.newRowGroup(5);
    builder.addNextValue(5, 2, 1.0);  // state 2, absorbing and labeled "goal"
    builder.newRowGroup(6);
    builder.addNextValue(6, 3, 1.0);  // state 3, absorbing
    auto const model = buildModel<storm::models::sparse::Mdp<ValueType>>(builder.build(), {{"goal", {2}}});
 
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient;
    EXPECT_EQ(3ull, quotient->getNumberOfStates());  // {0,1}, {2}, {3}
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "Pmin=? [F \"goal\"]"), 0.5, 1e-12);
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "Pmax=? [F \"goal\"]"), 1.0, 1e-12);
}

TEST(StrongBisimulationTest, ChoiceOrderMattersIfActionSensitive) {
    // States 0 and 1 offer the same two distributions over {2} and {3}, but in a different order.
    storm::storage::SparseMatrixBuilder<ValueType> builder(6, 4, 0, false, true, 4);
    builder.newRowGroup(0);
    builder.addNextValue(0, 2, 1.0);  // state 0, choice 0
    builder.addNextValue(1, 2, 0.5);  // state 0, choice 1
    builder.addNextValue(1, 3, 0.5);
    builder.newRowGroup(2);
    builder.addNextValue(2, 2, 0.5);  // state 1, choice 0
    builder.addNextValue(2, 3, 0.5);
    builder.addNextValue(3, 2, 1.0);  // state 1, choice 1
    builder.newRowGroup(4);
    builder.addNextValue(4, 2, 1.0);  // state 2, absorbing and labeled "goal"
    builder.newRowGroup(5);
    builder.addNextValue(5, 3, 1.0);  // state 3, absorbing
    auto const model = buildModel<storm::models::sparse::Mdp<ValueType>>(builder.build(), {{"goal", {2}}});
 
    EXPECT_EQ(3ull, storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient->getNumberOfStates());
 
    Options options = strongOptions();
    options.actionSensitive = true;
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, options).quotient;
    EXPECT_EQ(4ull, quotient->getNumberOfStates());
    EXPECT_EQ(6ull, quotient->getNumberOfChoices());
}

TEST(StrongBisimulationTest, ApproximateChoicePartnersAtDifferentPositions) {
    // States 0 and 1 have three choices each over "goal" (2) and "sink" (3). Every choice has a partner within 0.01 in the other state, but the middle
    // choices are 0.011 apart.
    storm::storage::SparseMatrixBuilder<ValueType> builder(8, 4, 0, false, true, 4);
    builder.newRowGroup(0);
    builder.addNextValue(0, 2, 0.001);
    builder.addNextValue(0, 3, 0.999);
    builder.addNextValue(1, 2, 0.019);
    builder.addNextValue(1, 3, 0.981);
    builder.addNextValue(2, 2, 0.039);
    builder.addNextValue(2, 3, 0.961);
    builder.newRowGroup(3);
    builder.addNextValue(3, 2, 0.010);
    builder.addNextValue(3, 3, 0.990);
    builder.addNextValue(4, 2, 0.030);
    builder.addNextValue(4, 3, 0.970);
    builder.addNextValue(5, 2, 0.048);
    builder.addNextValue(5, 3, 0.952);
    builder.newRowGroup(6);
    builder.addNextValue(6, 2, 1.0);
    builder.newRowGroup(7);
    builder.addNextValue(7, 3, 1.0);
    auto const model = buildModel<storm::models::sparse::Mdp<ValueType>>(builder.build(), {{"goal", {2}}, {"sink", {3}}});
 
    EXPECT_EQ(4ull, storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient->getNumberOfStates());
    Options options = strongOptions();
    options.tolerance = storm::utility::convertNumber<storm::RationalNumber>(0.02);  // i.e., halfTolerance 0.01
    EXPECT_EQ(3ull, storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, options).quotient->getNumberOfStates());
 
    // The choices 0.001 and 0.019 of state 0 both only have the partner 0.010 in state 1, so the quotient choices cannot be matched with the choices of
    // state 1 one-to-one.
    options.createQuotientChoiceMapping = true;
    STORM_SILENT_EXPECT_THROW(storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, options), storm::exceptions::UnexpectedException);
}

TEST(StrongBisimulationTest, ApproximateChoiceMappingIsOneToOne) {
    // The states 0 and 1 (and, analogously, 4 and 5) have three choices each over "goal" (2) and "sink" (3). With tolerance 0.2, the two states of such a
    // pair are bisimilar (their choices are pairwise 0.05 apart), while the choices of a single state are more than 0.1 apart. Some choices are within the
    // tolerance of two choices of the other state, e.g. 0.15 of state 1 is within 0.2 of both, 0.1 and 0.32 of state 0. The second pair mirrors the first one
    // (p vs. 1-p), so that one of them is affected no matter how the signatures are sorted and which state of a pair is the representative.
    storm::storage::SparseMatrixBuilder<ValueType> builder(14, 6, 0, true, true, 6);  // States 4 and 5 have no predecessors, so we force the dimensions.
    auto const addChoice = [&builder](uint64_t const row, ValueType const goalProbability) {
        builder.addNextValue(row, 2, goalProbability);
        builder.addNextValue(row, 3, storm::utility::one<ValueType>() - goalProbability);
    };
    builder.newRowGroup(0);  // state 0
    addChoice(0, 0.1);
    addChoice(1, 0.32);
    addChoice(2, 0.5);
    builder.newRowGroup(3);  // state 1
    addChoice(3, 0.15);
    addChoice(4, 0.27);
    addChoice(5, 0.45);
    builder.newRowGroup(6);
    builder.addNextValue(6, 2, 1.0);  // state 2, absorbing and labeled "goal"
    builder.newRowGroup(7);
    builder.addNextValue(7, 3, 1.0);  // state 3, absorbing and labeled "sink"
    builder.newRowGroup(8);           // state 4
    addChoice(8, 0.5);
    addChoice(9, 0.68);
    addChoice(10, 0.9);
    builder.newRowGroup(11);  // state 5
    addChoice(11, 0.55);
    addChoice(12, 0.73);
    addChoice(13, 0.85);
    auto const model = buildModel<storm::models::sparse::Mdp<ValueType>>(builder.build(), {{"goal", {2}}, {"sink", {3}}});
 
    Options options = strongOptions();
    options.tolerance = storm::utility::convertNumber<storm::RationalNumber>(0.2);  // i.e., halfTolerance 0.1
    options.createQuotientChoiceMapping = true;
    auto const result = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, options);
    ASSERT_EQ(4ull, result.quotient->getNumberOfStates());  // {0,1}, {2}, {3}, {4,5}
    ASSERT_EQ(8ull, result.quotient->getNumberOfChoices());
    checkChoiceMapping(*model, result);
 
    // Within a pair, the i'th choices of the two states are only close to each other, so they share their quotient choice.
    auto const& mapping = *result.toQuotientChoiceMapping;
    for (uint64_t i = 0; i < 3; ++i) {
        EXPECT_EQ(mapping[i], mapping[3 + i]) << "unexpected quotient choice for the choices " << i << " and " << (3 + i);
        EXPECT_EQ(mapping[8 + i], mapping[11 + i]) << "unexpected quotient choice for the choices " << (8 + i) << " and " << (11 + i);
    }
}

TEST(StrongBisimulationTest, MdpQuotientIsNotMeasureDriven) {
    // A chain 0 -> 1 -> 2 -> 2 in which the goal is reached with probability one, plus a second choice in state 0 that takes a detour through state 3.
    storm::storage::SparseMatrixBuilder<ValueType> builder(5, 4, 0, false, true, 4);
    builder.newRowGroup(0);
    builder.addNextValue(0, 1, 1.0);  // state 0, choice 0
    builder.addNextValue(1, 3, 1.0);  // state 0, choice 1
    builder.newRowGroup(2);
    builder.addNextValue(2, 2, 1.0);  // state 1
    builder.newRowGroup(3);
    builder.addNextValue(3, 2, 1.0);  // state 2, absorbing and labeled "goal"
    builder.newRowGroup(4);
    builder.addNextValue(4, 1, 1.0);  // state 3
    auto const model = buildModel<storm::models::sparse::Mdp<ValueType>>(builder.build(), {{"goal", {2}}});
 
    storm::parser::FormulaParser formulaParser;
    std::vector<std::shared_ptr<storm::logic::Formula const>> const formulas{formulaParser.parseSingleFormulaFromString("Pmax=? [F \"goal\"]")};
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, formulas, strongOptions()).quotient;
 
    // Every state reaches the goal with probability one, so a measure-driven initial partition would have merged all four states into a single one.
    EXPECT_EQ(4ull, quotient->getNumberOfStates());
    // No quotient state may carry both "init" and "goal", because no state of the original model does.
    auto const& labeling = quotient->getStateLabeling();
    ASSERT_TRUE(labeling.containsLabel("init"));
    ASSERT_TRUE(labeling.containsLabel("goal"));
    EXPECT_TRUE(labeling.getStates("init").isDisjointFrom(labeling.getStates("goal")));
    // A step-bounded property distinguishes the states that a measure-driven partition would have merged.
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "Pmax=? [F<=1 \"goal\"]"), checkFormula<ValueType>(model, "Pmax=? [F<=1 \"goal\"]"), 1e-12);
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "Pmax=? [F<=2 \"goal\"]"), checkFormula<ValueType>(model, "Pmax=? [F<=2 \"goal\"]"), 1e-12);
}

TEST(StrongBisimulationTest, MarkovAutomaton) {
    auto const model = storm::api::buildExplicitDRNModel<ValueType>(STORM_TEST_RESOURCES_DIR "/ma/jobscheduler.drn");
    ASSERT_EQ(storm::models::ModelType::MarkovAutomaton, model->getType());
 
    storm::parser::FormulaParser formulaParser;
    std::vector<std::shared_ptr<storm::logic::Formula const>> const formulas{formulaParser.parseSingleFormulaFromString("Tmin=? [F \"all_jobs_finished\"]")};
 
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, formulas, strongOptions()).quotient;
    EXPECT_EQ(storm::models::ModelType::MarkovAutomaton, quotient->getType());
    EXPECT_LE(quotient->getNumberOfStates(), model->getNumberOfStates());
    auto const markovAutomaton = quotient->template as<storm::models::sparse::MarkovAutomaton<ValueType>>();
    // Markovian and probabilistic states must not be mixed, and a Markovian state has exactly one choice.
    for (uint64_t state = 0; state < quotient->getNumberOfStates(); ++state) {
        if (markovAutomaton->isMarkovianState(state)) {
            EXPECT_EQ(1ull, quotient->getTransitionMatrix().getRowGroupSize(state));
        }
    }
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "Tmin=? [F \"all_jobs_finished\"]"), checkFormula<ValueType>(model, "Tmin=? [F \"all_jobs_finished\"]"),
                1e-9);
    // Expected times are preserved in both optimization directions, cf. https://github.com/stormchecker/storm/issues/498.
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "Tmax=? [F \"all_jobs_finished\"]"), checkFormula<ValueType>(model, "Tmax=? [F \"all_jobs_finished\"]"),
                1e-9);
}

TEST(StrongBisimulationTest, MarkovAutomatonExpectedTime) {
    // The probabilistic state 0 chooses between the Markovian states 1 and 2, which both move on to the "done" state 3 but at different rates.
    storm::storage::SparseMatrixBuilder<ValueType> builder(5, 4, 0, false, true, 4);
    builder.newRowGroup(0);
    builder.addNextValue(0, 1, 1.0);  // state 0 is probabilistic, so its rows hold probabilities: choice 0
    builder.addNextValue(1, 2, 1.0);  // state 0, choice 1
    builder.newRowGroup(2);
    builder.addNextValue(2, 3, 1.0);  // state 1, Markovian with exit rate 1
    builder.newRowGroup(3);
    builder.addNextValue(3, 3, 2.0);  // state 2, Markovian with exit rate 2
    builder.newRowGroup(4);
    builder.addNextValue(4, 3, 1.0);  // state 3, Markovian, absorbing and labeled "done"
    auto const model = buildMarkovAutomaton<ValueType>(builder.build(), storm::storage::BitVector(4, {1, 2, 3}), {{"done", {3}}});
    ASSERT_NEAR(0.5, checkFormula<ValueType>(model, "Tmin=? [F \"done\"]"), 1e-12);
    ASSERT_NEAR(1.0, checkFormula<ValueType>(model, "Tmax=? [F \"done\"]"), 1e-12);
 
    storm::parser::FormulaParser formulaParser;
    std::vector<std::shared_ptr<storm::logic::Formula const>> const formulas{formulaParser.parseSingleFormulaFromString("Tmin=? [F \"done\"]")};
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, formulas, strongOptions()).quotient;
 
    // States 1 and 2 have the same distribution over the blocks, so only their exit rates keep them apart.
    EXPECT_EQ(4ull, quotient->getNumberOfStates());
    EXPECT_NEAR(0.5, checkFormula<ValueType>(quotient, "Tmin=? [F \"done\"]"), 1e-12);
    EXPECT_NEAR(1.0, checkFormula<ValueType>(quotient, "Tmax=? [F \"done\"]"), 1e-12);
}

TEST(StrongBisimulationTest, MarkovAutomatonHybridStates) {
    // The hybrid states 0 and 1 behave the same: both move to the Markovian state 2 (which moves on to the "goal" state 3), either through their Markovian
    // choice or, without letting time pass, through their probabilistic one.
    storm::storage::SparseMatrixBuilder<ValueType> builder(6, 4, 0, false, true, 4);
    builder.newRowGroup(0);
    builder.addNextValue(0, 2, 2.0);  // The Markovian choice of state 0, holding a rate.
    builder.addNextValue(1, 2, 1.0);  // The probabilistic choice of state 0, with the same distribution.
    builder.newRowGroup(2);
    builder.addNextValue(2, 2, 2.0);
    builder.addNextValue(3, 2, 1.0);
    builder.newRowGroup(4);
    builder.addNextValue(4, 3, 4.0);
    builder.newRowGroup(5);
    builder.addNextValue(5, 3, 1.0);
    storm::storage::BitVector markovianStates(4);
    markovianStates.set(0);
    markovianStates.set(1);
    markovianStates.set(2);
    auto const model = buildMarkovAutomaton<ValueType>(builder.build(), std::move(markovianStates), {{"goal", {3}}});
    ASSERT_FALSE(model->isClosed());
 
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, strongOptions()).quotient;
    EXPECT_EQ(3ull, quotient->getNumberOfStates());   // {0, 1}, {2} and {3}
    EXPECT_EQ(4ull, quotient->getNumberOfChoices());  // The merged hybrid state keeps both of its choices.
    auto const markovAutomaton = quotient->template as<storm::models::sparse::MarkovAutomaton<ValueType>>();
    ASSERT_FALSE(markovAutomaton->isClosed());
    auto const hybridState = *quotient->getInitialStates().begin();
    EXPECT_TRUE(markovAutomaton->isMarkovianState(hybridState));
    EXPECT_EQ(2ull, quotient->getTransitionMatrix().getRowGroupSize(hybridState));
    EXPECT_EQ(2.0, markovAutomaton->getExitRate(hybridState));
 
    // Closing both models applies the maximal progress assumption, i.e. it keeps the probabilistic choice of a hybrid state. This also checks that the
    // Markovian choice is the first choice of the quotient state, since that is the one close() removes.
    model->close();
    markovAutomaton->close();
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "Tmin=? [F \"goal\"]"), 0.25, 1e-9);
    EXPECT_NEAR(checkFormula<ValueType>(model, "Tmin=? [F \"goal\"]"), 0.25, 1e-9);
}
 
// ------------------------------------------------------------
// Benchmark models
// ------------------------------------------------------------
 
TEST(StrongBisimulationTest, Die) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm", "P=? [F \"one\"]", 13ull, 5ull, 8ull, 5ull);
}
 
TEST(StrongBisimulationTest, DieAllLabels) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm", "P=? [F \"one\"]", 13ull, 11ull, 17ull, 11ull, allLabelOptions());
}

TEST(StrongBisimulationTest, DieUnnamedRewardModel) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm", "R=? [F \"done\"]", 13ull, 5ull, 7ull, 5ull);
}

TEST(StrongBisimulationTest, DiePropositionalSubformulas) {
    std::string const formulaString = "P=? [F (s=7 & d=1) | \"two\"]";
    testQuotient(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm", formulaString, 13ull, 6ull, 10ull, 6ull);
    Options options = strongOptions();
    options.stateLabelPreservation = StateLabelPreservation::FormulaIndividual;
    testQuotient(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm", formulaString, 13ull, 7ull, 11ull, 7ull, options);
}

TEST(StrongBisimulationTest, DieStateValuations) {
#ifndef STORM_HAVE_Z3
    GTEST_SKIP() << "Z3 not available.";
#endif
    auto const input = buildFromPrism<ValueType>(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm", "P=? [F \"one\"]", {.stateValuations = true});
    ASSERT_TRUE(input.model->hasStateValuations());
 
    auto const result = storm::bisimulation::performBisimulationMinimization<ValueType>(*input.model, input.formulas, strongOptions());
    ASSERT_TRUE(result.quotient->hasStateValuations());
    ASSERT_EQ(result.quotient->getNumberOfStates(), result.quotient->getStateValuations().getNumberOfEntities());
    // Every quotient state carries the valuation of one of the states that it represents.
    std::set<uint64_t> quotientStatesWithRepresentedValuation;
    for (uint64_t state = 0; state < input.model->getNumberOfStates(); ++state) {
        uint64_t const quotientState = result.toQuotientStateMapping[state];
        if (input.model->getStateValuations().toString(state) == result.quotient->getStateValuations().toString(quotientState)) {
            quotientStatesWithRepresentedValuation.insert(quotientState);
        }
    }
    EXPECT_EQ(result.quotient->getNumberOfStates(), quotientStatesWithRepresentedValuation.size());
}
 
TEST(StrongBisimulationTest, Crowds) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/dtmc/crowds5_5.pm", "P=? [F \"observe0Greater1\"]", 7403ull, 65ull, 105ull, 65ull);
}
 
TEST(StrongBisimulationTest, CrowdsAllLabels) {
    // The model is larger than above because the target states of the formula are no longer made absorbing.
    testQuotient(STORM_TEST_RESOURCES_DIR "/dtmc/crowds5_5.pm", "P=? [F \"observe0Greater1\"]", 8607ull, 2149ull, 3912ull, 2149ull, allLabelOptions());
}
 
TEST(StrongBisimulationTest, CtmcEmbedded) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm", "P=? [F<=10000 \"down\"]", 2076ull, 310ull, 1752ull, 310ull, approximateOptions());
}
 
TEST(StrongBisimulationTest, CtmcCluster) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/ctmc/cluster2.sm", "P=? [F<=100 !\"minimum\"]", 276ull, 147ull, 569ull, 147ull);
}
 
TEST(StrongBisimulationTest, CtmcClusterChoiceLabelsAndOrigins) {
    // The actions of the left and the right cluster have different names and stem from different modules, which rules out any reduction.
    testQuotient(STORM_TEST_RESOURCES_DIR "/ctmc/cluster2.sm", "P=? [F<=100 !\"minimum\"]", 276ull, 276ull, 1120ull, 276ull, strongOptions(),
                 {.choiceLabels = true, .choiceOrigins = true});
}
 
TEST(StrongBisimulationTest, DtmcBrpChoiceLabelsAndOrigins) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/dtmc/brp-16-2.pm", "P=? [F \"target\"]", 613ull, 412ull, 572ull, 412ull, strongOptions(),
                 {.choiceLabels = true, .choiceOrigins = true});
}

TEST(StrongBisimulationTest, CtmcEmbeddedExact) {
#ifndef STORM_HAVE_Z3
    GTEST_SKIP() << "Z3 not available.";
#endif
    std::string const formulaString = "P=? [F<=10000 \"down\"]";
    auto const exactInput = buildFromPrism<storm::RationalNumber>(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm", formulaString, {.allLabels = true});
    ASSERT_EQ(3478ull, exactInput.model->getNumberOfStates());
    ASSERT_EQ(14639ull, exactInput.model->getNumberOfTransitions());
 
    Options options = strongOptions();
    options.stateLabelPreservation = StateLabelPreservation::All;
    auto const labeledQuotient =
        storm::bisimulation::performBisimulationMinimization<storm::RationalNumber>(*exactInput.model, exactInput.formulas, options).quotient;
    EXPECT_EQ(1127ull, labeledQuotient->getNumberOfStates());
    EXPECT_EQ(5730ull, labeledQuotient->getNumberOfTransitions());
 
    // Without the labels the quotient is much coarser. This is the configuration in which the reported sizes differed the most (136 on macOS vs 180 on
    // Linux), since the initial partition consists of a single block there.
    options.stateLabelPreservation = StateLabelPreservation::None;
    auto const unlabeledQuotient = storm::bisimulation::performBisimulationMinimization<storm::RationalNumber>(*exactInput.model, {}, options).quotient;
    EXPECT_EQ(98ull, unlabeledQuotient->getNumberOfStates());
    EXPECT_EQ(539ull, unlabeledQuotient->getNumberOfTransitions());
}

TEST(StrongBisimulationTest, CtmcEmbeddedWithTolerance) {
#ifndef STORM_HAVE_Z3
    GTEST_SKIP() << "Z3 not available.";
#endif
    std::string const formulaString = "P=? [F<=10000 \"down\"]";
    auto const input = buildFromPrism<double>(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm", formulaString, {.allLabels = true});
 
    Options options = approximateOptions();
    options.stateLabelPreservation = StateLabelPreservation::All;
    auto const labeledQuotient = storm::bisimulation::performBisimulationMinimization<double>(*input.model, input.formulas, options).quotient;
    EXPECT_EQ(1127ull, labeledQuotient->getNumberOfStates());
    EXPECT_EQ(5730ull, labeledQuotient->getNumberOfTransitions());
 
    options.stateLabelPreservation = StateLabelPreservation::None;
    auto const unlabeledQuotient = storm::bisimulation::performBisimulationMinimization<double>(*input.model, {}, options).quotient;
    EXPECT_EQ(98ull, unlabeledQuotient->getNumberOfStates());
    EXPECT_EQ(539ull, unlabeledQuotient->getNumberOfTransitions());
}
 
TEST(StrongBisimulationTest, MdpTwoDice) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm", "Pmin=? [F \"two\"]", 169ull, 11ull, 26ull, 14ull);
}
 
TEST(StrongBisimulationTest, MdpTwoDiceAllLabels) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm", "Pmin=? [F \"two\"]", 169ull, 77ull, 183ull, 97ull, allLabelOptions());
}
 
TEST(StrongBisimulationTest, MdpCoin) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/coin2-2.nm", "Pmin=? [F \"finished\"]", 272ull, 55ull, 96ull, 78ull);
}

TEST(StrongBisimulationTest, MdpCoinQuotientIsNotMeasureDriven) {
#ifndef STORM_HAVE_Z3
    GTEST_SKIP() << "Z3 not available.";
#endif
    std::string const formulaString = "Pmax=? [F \"finished\"]";
    auto const input = buildFromPrism<ValueType>(STORM_TEST_RESOURCES_DIR "/mdp/coin2-2.nm", formulaString);
    ASSERT_NEAR(1.0, checkFormula<ValueType>(input.model, formulaString), 1e-9);
 
    auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*input.model, input.formulas, strongOptions()).quotient;
    EXPECT_EQ(55ull, quotient->getNumberOfStates());
    EXPECT_NEAR(checkFormula<ValueType>(quotient, formulaString), 1.0, 1e-9);
    // The quotient preserves all of PCTL, in particular the step-bounded variant of the formula it was built for.
    EXPECT_NEAR(checkFormula<ValueType>(quotient, "Pmax=? [F<=10 \"finished\"]"), checkFormula<ValueType>(input.model, "Pmax=? [F<=10 \"finished\"]"), 1e-9);
}
 
TEST(StrongBisimulationTest, MdpTwoDiceActionSensitive) {
    Options options = strongOptions();
    options.actionSensitive = true;
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm", "Pmin=? [F \"two\"]", 169ull, 33ull, 97ull, 58ull, options);
}
 
TEST(StrongBisimulationTest, MdpCoinChoiceOrigins) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/coin2-2.nm", "Pmin=? [F \"finished\"]", 272ull, 195ull, 399ull, 321ull, strongOptions(),
                 {.choiceOrigins = true});
}
 
TEST(StrongBisimulationTest, MdpLeaderChoiceLabels) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/leader3.nm", "Pmin=? [F \"elected\"]", 364ull, 169ull, 325ull, 262ull, strongOptions(), {.choiceLabels = true});
}
 
TEST(StrongBisimulationTest, MdpLeaderChoiceLabelsAndOriginsActionSensitive) {
    Options options = strongOptions();
    options.actionSensitive = true;
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/leader3.nm", "Pmin=? [F \"elected\"]", 364ull, 169ull, 325ull, 262ull, options,
                 {.choiceLabels = true, .choiceOrigins = true});
}

TEST(StrongBisimulationTest, MdpTwoDiceApproximative) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm", "Pmin=? [F \"two\"]", 169ull, 11ull, 26ull, 14ull, approximateOptions());
}
 
TEST(StrongBisimulationTest, MdpCoinApproximative) {
    testQuotient(STORM_TEST_RESOURCES_DIR "/mdp/coin2-2.nm", "Pmin=? [F \"finished\"]", 272ull, 55ull, 96ull, 78ull, approximateOptions());
}
 
// ------------------------------------------------------------
// Consistency between the refinement strategies and the value types
// ------------------------------------------------------------

template<typename VT>
void testRefinementStrategiesAgree(std::string const& prismFile, std::string const& formulaString) {
#ifndef STORM_HAVE_Z3
    GTEST_SKIP() << "Z3 not available.";
#endif
    auto const input = buildFromPrism<VT>(prismFile, formulaString);
 
    auto splitterOptions = strongOptions();
    splitterOptions.preferSignatureRefinement = false;
    auto const splitterQuotient = storm::bisimulation::performBisimulationMinimization<VT>(*input.model, input.formulas, splitterOptions).quotient;
 
    auto signatureOptions = splitterOptions;
    signatureOptions.preferSignatureRefinement = true;
    auto const signatureQuotient = storm::bisimulation::performBisimulationMinimization<VT>(*input.model, input.formulas, signatureOptions).quotient;
 
    EXPECT_EQ(splitterQuotient->getNumberOfStates(), signatureQuotient->getNumberOfStates());
    EXPECT_EQ(splitterQuotient->getNumberOfTransitions(), signatureQuotient->getNumberOfTransitions());
}
 
TEST(StrongBisimulationTest, RefinementStrategiesAgreeOnDie) {
    testRefinementStrategiesAgree<ValueType>(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm", "P=? [F \"one\"]");
}
 
TEST(StrongBisimulationTest, RefinementStrategiesAgreeOnCrowds) {
    testRefinementStrategiesAgree<ValueType>(STORM_TEST_RESOURCES_DIR "/dtmc/crowds5_5.pm", "P=? [F \"observe0Greater1\"]");
}

TEST(StrongBisimulationTest, RefinementStrategiesAgreeOnCtmc) {
    testRefinementStrategiesAgree<storm::RationalNumber>(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm", "P=? [F<=10000 \"down\"]");
}

void testExactAgreesWithApproximate(std::string const& prismFile, std::string const& formulaString) {
#ifndef STORM_HAVE_Z3
    GTEST_SKIP() << "Z3 not available.";
#endif
    auto const doubleInput = buildFromPrism<double>(prismFile, formulaString);
    auto const exactInput = buildFromPrism<storm::RationalNumber>(prismFile, formulaString);
 
    auto const doubleQuotient =
        storm::bisimulation::performBisimulationMinimization<double>(*doubleInput.model, doubleInput.formulas, strongOptions()).quotient;
    auto const exactQuotient =
        storm::bisimulation::performBisimulationMinimization<storm::RationalNumber>(*exactInput.model, exactInput.formulas, strongOptions()).quotient;
 
    EXPECT_EQ(doubleQuotient->getNumberOfStates(), exactQuotient->getNumberOfStates());
    EXPECT_EQ(doubleQuotient->getNumberOfTransitions(), exactQuotient->getNumberOfTransitions());
}
 
TEST(StrongBisimulationTest, ExactAgreesWithApproximateOnDie) {
    testExactAgreesWithApproximate(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm", "P=? [F \"one\"]");
}
 
TEST(StrongBisimulationTest, ExactAgreesWithApproximateOnTwoDice) {
    testExactAgreesWithApproximate(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm", "Pmin=? [F \"two\"]");
}

template<typename VT>
void testPermutationInvariance(std::string const& prismFile, std::string const& formulaString, uint64_t const numSeeds, Options options = strongOptions()) {
#ifndef STORM_HAVE_Z3
    GTEST_SKIP() << "Z3 not available.";
#endif
    options.stateLabelPreservation = StateLabelPreservation::All;
    auto const input = buildFromPrism<VT>(prismFile, formulaString, {.allLabels = true});
    auto const quotient = storm::bisimulation::performBisimulationMinimization<VT>(*input.model, input.formulas, options).quotient;
 
    std::vector<uint64_t> permutation(input.model->getNumberOfStates());
    std::iota(permutation.begin(), permutation.end(), 0ull);
    for (uint64_t seed = 0; seed < numSeeds; ++seed) {
        std::mt19937_64 rng(seed);
        std::shuffle(permutation.begin(), permutation.end(), rng);
        auto const permutedModel = storm::transformer::permuteStates(*input.model, permutation);
        auto const permutedQuotient = storm::bisimulation::performBisimulationMinimization<VT>(*permutedModel, input.formulas, options).quotient;
        EXPECT_EQ(quotient->getNumberOfStates(), permutedQuotient->getNumberOfStates()) << "Seed " << seed << ".";
        EXPECT_EQ(quotient->getNumberOfTransitions(), permutedQuotient->getNumberOfTransitions()) << "Seed " << seed << ".";
    }
}
 
TEST(StrongBisimulationTest, PermutationInvarianceExactCtmc) {
    testPermutationInvariance<storm::RationalNumber>(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm", "P=? [F<=10000 \"down\"]", 3ull);
}
 
TEST(StrongBisimulationTest, PermutationInvarianceCtmc) {
    testPermutationInvariance<double>(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm", "P=? [F<=10000 \"down\"]", 3ull, approximateOptions());
}
 
TEST(StrongBisimulationTest, PermutationInvarianceMdp) {
    testPermutationInvariance<double>(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm", "Pmin=? [F \"two\"]", 3ull, approximateOptions());
}
 
// ------------------------------------------------------------
// Randomized self-checks
// ------------------------------------------------------------

template<typename ModelType>
std::shared_ptr<ModelType> buildRandomModel(uint64_t const numStates, uint64_t const seed) {
    bool constexpr isNondeterministic = std::is_same_v<ModelType, storm::models::sparse::Mdp<ValueType>>;
    std::mt19937_64 rng(seed);
    storm::storage::SparseMatrixBuilder<ValueType> builder(0, numStates, 0, false, isNondeterministic, isNondeterministic ? numStates : 0);
    uint64_t row = 0;
    for (uint64_t state = 0; state < numStates; ++state) {
        if constexpr (isNondeterministic) {
            builder.newRowGroup(row);
        }
        uint64_t const numChoices = isNondeterministic ? 1 + rng() % 2 : 1;
        for (uint64_t choice = 0; choice < numChoices; ++choice, ++row) {
            std::map<uint64_t, ValueType> distribution;
            uint64_t const numSuccessors = 1 + rng() % 3;
            for (uint64_t i = 0; i < numSuccessors; ++i) {
                // Small values keep the rows short and make coinciding distributions (and thus actual merges) reasonably likely.
                distribution[rng() % numStates] += static_cast<ValueType>(1 + rng() % 4);
            }
            ValueType sum = storm::utility::zero<ValueType>();
            for (auto const& [_, value] : distribution) {
                sum += value;
            }
            for (auto const& [column, value] : distribution) {
                builder.addNextValue(row, column, value / sum);
            }
        }
    }
    std::map<std::string, std::vector<uint64_t>> labels{{"a", {}}, {"b", {}}};
    for (uint64_t state = 0; state < numStates; ++state) {
        if (rng() % 3 == 0) {
            labels["a"].push_back(state);
        }
        if (rng() % 4 == 0) {
            labels["b"].push_back(state);
        }
    }
    return buildModel<ModelType>(builder.build(row, numStates, numStates), labels);
}

template<typename ModelType>
void testIdempotence(uint64_t const numStates, uint64_t const numSeeds, std::string const& formulaString) {
    for (uint64_t seed = 0; seed < numSeeds; ++seed) {
        auto const model = buildRandomModel<ModelType>(numStates, seed);
        Options const options = strongOptions();
        auto const quotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, options).quotient;
        ASSERT_LE(quotient->getNumberOfStates(), model->getNumberOfStates()) << "Seed " << seed << ".";
        auto const quotientOfQuotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*quotient, {}, options).quotient;
        EXPECT_EQ(quotient->getNumberOfStates(), quotientOfQuotient->getNumberOfStates()) << "Not idempotent for seed " << seed << ".";
        EXPECT_EQ(quotient->getNumberOfTransitions(), quotientOfQuotient->getNumberOfTransitions()) << "Not idempotent for seed " << seed << ".";
        // The comparison is against the default precision of the underlying equation solver, not against the exactness of the quotient.
        EXPECT_NEAR(checkFormula<ValueType>(quotient, formulaString), checkFormula<ValueType>(model, formulaString), 1e-4) << "Seed " << seed << ".";
    }
}
 
TEST(StrongBisimulationTest, RandomDtmcs) {
    testIdempotence<storm::models::sparse::Dtmc<ValueType>>(12, 500, "P=? [F \"a\"]");
    testIdempotence<storm::models::sparse::Dtmc<ValueType>>(60, 200, "P=? [(!\"b\") U \"a\"]");
}
 
TEST(StrongBisimulationTest, RandomCtmcs) {
    testIdempotence<storm::models::sparse::Ctmc<ValueType>>(12, 500, "P=? [F \"a\"]");
    testIdempotence<storm::models::sparse::Ctmc<ValueType>>(60, 200, "P=? [(!\"b\") U \"a\"]");
}
 
TEST(StrongBisimulationTest, RandomMdps) {
    testIdempotence<storm::models::sparse::Mdp<ValueType>>(12, 500, "Pmin=? [F \"a\"]");
    testIdempotence<storm::models::sparse::Mdp<ValueType>>(60, 200, "Pmax=? [(!\"b\") U \"a\"]");
}

template<typename ModelType>
void testRefinementStrategiesAgreeOnRandomModels(uint64_t const numStates, uint64_t const numSeeds) {
    for (uint64_t seed = 0; seed < numSeeds; ++seed) {
        auto const model = buildRandomModel<ModelType>(numStates, seed);
        auto splitterOptions = strongOptions();
        splitterOptions.preferSignatureRefinement = false;
        auto const splitterQuotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, splitterOptions).quotient;
        auto signatureOptions = splitterOptions;
        signatureOptions.preferSignatureRefinement = true;
        auto const signatureQuotient = storm::bisimulation::performBisimulationMinimization<ValueType>(*model, {}, signatureOptions).quotient;
        EXPECT_EQ(splitterQuotient->getNumberOfStates(), signatureQuotient->getNumberOfStates()) << "Seed " << seed << ".";
        EXPECT_EQ(splitterQuotient->getNumberOfTransitions(), signatureQuotient->getNumberOfTransitions()) << "Seed " << seed << ".";
    }
}
 
TEST(StrongBisimulationTest, RefinementStrategiesAgreeOnRandomDtmcs) {
    testRefinementStrategiesAgreeOnRandomModels<storm::models::sparse::Dtmc<ValueType>>(12, 500);
    testRefinementStrategiesAgreeOnRandomModels<storm::models::sparse::Dtmc<ValueType>>(60, 200);
}
 
TEST(StrongBisimulationTest, RefinementStrategiesAgreeOnRandomCtmcs) {
    testRefinementStrategiesAgreeOnRandomModels<storm::models::sparse::Ctmc<ValueType>>(12, 500);
    testRefinementStrategiesAgreeOnRandomModels<storm::models::sparse::Ctmc<ValueType>>(60, 200);
}
 
}  // namespace