Merge remote-tracking branch 'origin/sylvanRationalFunctions' into menu_games

Former-commit-id: 258428ac4b

resources/3rdparty/sylvan/src/storm_function_wrapper.cpp

@@ -3,9 +3,16 @@
 #include <cstring>
 #include <iostream>
 #include <sstream>
+#include <set>
+#include <map>
 #include "src/adapters/CarlAdapter.h"
 #include "sylvan_storm_rational_function.h"
 
+#include <sylvan_config.h>
+#include <sylvan.h>
+#include <sylvan_common.h>
+#include <sylvan_mtbdd.h>
+
 #undef DEBUG_STORM_FUNCTION_WRAPPER
 
 #ifdef DEBUG_STORM_FUNCTION_WRAPPER
@@ -206,16 +213,37 @@ void print_storm_rational_function_to_file(storm_rational_function_ptr a, FILE*
 	fprintf(out, "%s", s.c_str());
 }
 
-MTBDD storm_rational_function_leaf_parameter_replacement(uint64_t node_value, uint32_t node_type, void* context) {
-	if (node_type != sylvan_storm_rational_function_get_type()) {
-		//
-	} else {
-		//
+MTBDD testiTest(storm::RationalFunction const& currentFunction, std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>> const& replacements) {
+	if (currentFunction.isConstant()) {
+		return mtbdd_storm_rational_function((storm_rational_function_ptr)&currentFunction);
+	}
+
+	std::set<storm::RationalFunctionVariable> variablesInFunction = currentFunction.gatherVariables();
+	std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>>::const_iterator it = replacements.cbegin();
+	std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>>::const_iterator end = replacements.cend();
+
+	// Walking the (ordered) map enforces an ordering on the MTBDD
+	for (; it != end; ++it) {
+		if (variablesInFunction.find(it->first) != variablesInFunction.cend()) {
+			std::map<storm::RationalFunctionVariable, storm::RationalNumber> highReplacement = {{it->first, it->second.second.first}};
+			std::map<storm::RationalFunctionVariable, storm::RationalNumber> lowReplacement = {{it->first, it->second.second.second}};
+			MTBDD high = testiTest(currentFunction.substitute(highReplacement), replacements);
+			MTBDD low = testiTest(currentFunction.substitute(lowReplacement), replacements);
+			LACE_ME
+			return mtbdd_ite(mtbdd_ithvar(it->second.first), high, low);
+		}
+	}
+
+	return mtbdd_storm_rational_function((storm_rational_function_ptr)&currentFunction);
+}
+
+
+MTBDD storm_rational_function_leaf_parameter_replacement(MTBDD dd, storm_rational_function_ptr a, void* context) {
+	storm::RationalFunction& srf_a = *(storm::RationalFunction*)a;
+	if (srf_a.isConstant()) {
+		return dd;
 	}
 	
-	(void)node_value;
-	(void)node_type;
-	(void)context;
-	
-	return mtbdd_invalid;
+	std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>>* replacements = (std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>>*)context;
+	return testiTest(srf_a, *replacements);
 }

resources/3rdparty/sylvan/src/storm_function_wrapper.h

@@ -31,7 +31,7 @@ void print_storm_rational_function_to_file(storm_rational_function_ptr a, FILE*
 
 int storm_rational_function_is_zero(storm_rational_function_ptr a);
 
-MTBDD storm_rational_function_leaf_parameter_replacement(uint64_t node_value, uint32_t node_type, void* context);
+MTBDD storm_rational_function_leaf_parameter_replacement(MTBDD dd, storm_rational_function_ptr a, void* context);
 
 #ifdef __cplusplus
 }

resources/3rdparty/sylvan/src/sylvan_obj_mtbdd_storm.hpp

@@ -30,6 +30,8 @@
     Mtbdd DivideRF(const Mtbdd &other) const;
 	
 	Mtbdd AbstractPlusRF(const BddSet &variables) const;
+	
+	Mtbdd ReplaceLeavesRF(void* context) const;
 #endif
 	
 	/**

resources/3rdparty/sylvan/src/sylvan_obj_storm.cpp

@@ -57,6 +57,11 @@ Mtbdd Mtbdd::AbstractPlusRF(const BddSet &variables) const {
     return sylvan_storm_rational_function_abstract_plus(mtbdd, variables.set.bdd);
 }
 
+Mtbdd Mtbdd::ReplaceLeavesRF(void* context) const {
+	LACE_ME;
+    return sylvan_storm_rational_function_replace_leaves(mtbdd, (size_t)context);
+}
+
 #endif
 
 Mtbdd

resources/3rdparty/sylvan/src/sylvan_storm_rational_function.c

@@ -380,7 +380,7 @@ TASK_IMPL_2(MTBDD, sylvan_storm_rational_function_op_neg, MTBDD, dd, size_t, p)
 /**
  * Operation "replace leaves" for one storm::RationalFunction MTBDD
  */
-TASK_IMPL_2(MTBDD, sylvan_storm_rational_function_op_replace_leaves, MTBDD, dd, void*, context)
+TASK_IMPL_2(MTBDD, sylvan_storm_rational_function_op_replace_leaves, MTBDD, dd, size_t, context)
 {
 	LOG_I("task_impl_2 op_replace")
     /* Handle partial functions */
@@ -388,7 +388,12 @@ TASK_IMPL_2(MTBDD, sylvan_storm_rational_function_op_replace_leaves, MTBDD, dd,
 
     /* Compute result for leaf */
     if (mtbdd_isleaf(dd)) {
-		return storm_rational_function_leaf_parameter_replacement(mtbdd_getvalue(dd), mtbdd_gettype(dd), context);
+		if (mtbdd_gettype(dd) != sylvan_storm_rational_function_type) {
+			assert(0);
+		}
+		
+		storm_rational_function_ptr mdd = (storm_rational_function_ptr)mtbdd_getvalue(dd);
+		return storm_rational_function_leaf_parameter_replacement(dd, mdd, (void*)context);
     }
 
     return mtbdd_invalid;

resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h

@@ -116,12 +116,12 @@ typedef MTBDD (*mtbddLeaveReplacementFunction)(uint64_t, uint32_t, void*);
 /**
  * Operation "replace" for one storm::RationalFunction MTBDD
  */
-TASK_DECL_2(MTBDD, sylvan_storm_rational_function_op_replace_leaves, MTBDD, void*)
+TASK_DECL_2(MTBDD, sylvan_storm_rational_function_op_replace_leaves, MTBDD, size_t)
 
 /**
  * Compute the MTBDD that arises from a after calling the mtbddLeaveReplacementFunction on each leaf.
  */
-#define sylvan_storm_rational_function_replace_leaves(a, func, ctx) mtbdd_uapply(a, TASK(sylvan_storm_rational_function_op_replace_leaves), ctx)
+#define sylvan_storm_rational_function_replace_leaves(a, ctx) mtbdd_uapply(a, TASK(sylvan_storm_rational_function_op_replace_leaves), ctx)
 
 #ifdef __cplusplus
 }

src/adapters/CarlAdapter.h

@@ -15,6 +15,7 @@
 #include <carl/core/FactorizedPolynomial.h>
 #include <carl/core/Relation.h>
 #include <carl/core/SimpleConstraint.h>
+#include <carl/util/stringparser.h>
 
 namespace carl {
     // Define hash values for all polynomials and rational function.
@@ -60,6 +61,7 @@ namespace storm {
     typedef carl::Variable RationalFunctionVariable;
     typedef carl::MultivariatePolynomial<RationalNumber> RawPolynomial;
     typedef carl::FactorizedPolynomial<RawPolynomial> Polynomial;
+	typedef carl::Cache<carl::PolynomialFactorizationPair<RawPolynomial>> RawPolynomialCache;
     typedef carl::Relation CompareRelation;
     
     typedef carl::RationalFunction<Polynomial, true> RationalFunction;

src/storage/dd/Add.cpp

@@ -784,6 +784,34 @@ namespace storm {
             return internalAdd;
         }
         
+#ifdef STORM_HAVE_CARL
+		template<DdType LibraryType, typename ValueType>
+        Add<LibraryType, ValueType> Add<LibraryType, ValueType>::replaceLeaves(std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>> const& replacementMap) const {
+			STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Not yet implemented: replaceLeaves");
+		}
+
+		template<>
+        Add<storm::dd::DdType::Sylvan, storm::RationalFunction> Add<storm::dd::DdType::Sylvan, storm::RationalFunction>::replaceLeaves(std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>> const& replacementMap) const {
+			std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>> internalReplacementMap;
+			std::set<storm::expressions::Variable> containedMetaVariables = this->getContainedMetaVariables();
+			
+			std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>>::const_iterator it = replacementMap.cbegin();
+			std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>>::const_iterator end = replacementMap.cend();
+			
+			for (; it != end; ++it) {
+				DdMetaVariable<storm::dd::DdType::Sylvan> const& metaVariable = this->getDdManager().getMetaVariable(it->second.first);
+				STORM_LOG_THROW(metaVariable.getNumberOfDdVariables() == 1, storm::exceptions::InvalidArgumentException, "Cannot use MetaVariable with more then one internal DD variable.");
+				
+				auto const& ddVariable = metaVariable.getDdVariables().at(0);
+				
+				internalReplacementMap.insert(std::make_pair(it->first, std::make_pair(ddVariable.getIndex(), it->second.second)));
+				containedMetaVariables.insert(it->second.first);
+			}
+			
+			return Add<storm::dd::DdType::Sylvan, storm::RationalFunction>(this->getDdManager(), internalAdd.replaceLeaves(internalReplacementMap), containedMetaVariables);
+		}
+#endif
+		
         template class Add<storm::dd::DdType::CUDD, double>;
         template class Add<storm::dd::DdType::CUDD, uint_fast64_t>;
 

src/storage/dd/Add.h

@@ -14,6 +14,9 @@
 #include "src/storage/dd/cudd/CuddAddIterator.h"
 #include "src/storage/dd/sylvan/SylvanAddIterator.h"
 
+#include "storm-config.h"
+#include "src/adapters/CarlAdapter.h"
+
 namespace storm {
     namespace dd {
         template<DdType LibraryType>
@@ -243,6 +246,16 @@ namespace storm {
              */
             Add<LibraryType, ValueType> maximum(Add<LibraryType, ValueType> const& other) const;
             
+#ifdef STORM_HAVE_CARL
+			/*!
+             * Replaces the leaves in this MTBDD, using the supplied variable replacement map.
+             *
+             * @param replacementMap The variable replacement map.
+             * @return The resulting function represented as an ADD.
+             */
+            Add<LibraryType, ValueType> replaceLeaves(std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>> const& replacementMap) const;
+#endif
+			
             /*!
              * Sum-abstracts from the given meta variables.
              *

src/storage/dd/sylvan/InternalSylvanAdd.cpp

@@ -289,6 +289,19 @@ namespace storm {
         }
 #endif
 
+        template<typename ValueType>
+        InternalAdd<DdType::Sylvan, ValueType> InternalAdd<DdType::Sylvan, ValueType>::replaceLeaves(std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>> const& replacementMap) const {
+            STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Not yet implemented: replaceLeaves");
+        }
+
+#ifdef STORM_HAVE_CARL
+		template<>
+        InternalAdd<DdType::Sylvan, storm::RationalFunction> InternalAdd<DdType::Sylvan, storm::RationalFunction>::replaceLeaves(std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>> const& replacementMap) const {
+            return InternalAdd<DdType::Sylvan, storm::RationalFunction>(ddManager, this->sylvanMtbdd.ReplaceLeavesRF((void*)&replacementMap));
+        }
+#endif
+
+
         template<typename ValueType>
         InternalAdd<DdType::Sylvan, ValueType> InternalAdd<DdType::Sylvan, ValueType>::sumAbstract(InternalBdd<DdType::Sylvan> const& cube) const {
             return InternalAdd<DdType::Sylvan, ValueType>(ddManager, this->sylvanMtbdd.AbstractPlus(cube.sylvanBdd));

src/storage/dd/sylvan/InternalSylvanAdd.h

@@ -248,6 +248,16 @@ namespace storm {
              */
             InternalAdd<DdType::Sylvan, ValueType> maximum(InternalAdd<DdType::Sylvan, ValueType> const& other) const;
             
+#ifdef STORM_HAVE_CARL
+			/*!
+             * Replaces the leaves in this MTBDD, using the supplied variable replacement map.
+             *
+             * @param replacementMap The variable replacement map.
+             * @return The resulting function represented as an ADD.
+             */
+            InternalAdd<DdType::Sylvan, ValueType> replaceLeaves(std::map<storm::RationalFunctionVariable, std::pair<uint32_t, std::pair<storm::RationalNumber, storm::RationalNumber>>> const& replacementMap) const;
+#endif
+			
             /*!
              * Sum-abstracts from the given cube.
              *

test/functional/storage/SylvanDdTest.cpp

@@ -163,7 +163,7 @@ TEST(SylvanDd, BddExistAbstractRepresentative) {
 	EXPECT_TRUE(bddX0Y0Z0 == representative_xyz);
 }
 
-TEST(SylvanDd, AddMinExistAbstractRepresentative) {
+TEST(SylvanDd, AddMinAbstractRepresentative) {
 	std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
 
 	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddZero;
@@ -268,7 +268,7 @@ TEST(SylvanDd, AddMinExistAbstractRepresentative) {
 	EXPECT_TRUE(representative_complex_xyz == comparison_complex_xyz);
 }
 
-TEST(SylvanDd, AddMaxExistAbstractRepresentative) {
+TEST(SylvanDd, AddMaxAbstractRepresentative) {
 	std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
 
 	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddZero;
@@ -414,6 +414,214 @@ TEST(SylvanDd, EncodingTest) {
 }
 
 #ifdef STORM_HAVE_CARL
+TEST(SylvanDd, RationalFunctionLeaveReplacementNonVariable) {
+	std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> zero;
+    ASSERT_NO_THROW(zero = manager->template getAddZero<storm::RationalFunction>());
+	
+	std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>> replacementMap;
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> zeroReplacementResult = zero.replaceLeaves(replacementMap);
+	
+	EXPECT_EQ(0ul, zeroReplacementResult.getNonZeroCount());
+    EXPECT_EQ(1ul, zeroReplacementResult.getLeafCount());
+    EXPECT_EQ(1ul, zeroReplacementResult.getNodeCount());
+	EXPECT_TRUE(zeroReplacementResult == zero);
+}
+
+TEST(SylvanDd, RationalFunctionLeaveReplacementSimpleVariable) {
+	std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
+	
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> function;
+	
+	std::shared_ptr<storm::RawPolynomialCache> cache = std::make_shared<storm::RawPolynomialCache>();
+	
+	carl::StringParser parser;
+	parser.setVariables({"x"});
+	
+	storm::RawPolynomial polyX = parser.parseMultivariatePolynomial<storm::RationalNumber>("x");	
+	storm::RationalFunction variableX = storm::RationalFunction(storm::Polynomial(polyX, cache));
+	
+	ASSERT_NO_THROW(function = manager->template getConstant<storm::RationalFunction>(variableX));
+	
+	std::pair<storm::expressions::Variable, storm::expressions::Variable> xExpr;
+	ASSERT_NO_THROW(xExpr = manager->addMetaVariable("x", 0, 1));
+	
+	std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>> replacementMap;
+	storm::RationalNumber rnOneThird = storm::RationalNumber(1) / storm::RationalNumber(3);
+	storm::RationalNumber rnTwoThird = storm::RationalNumber(2) / storm::RationalNumber(3);
+	replacementMap.insert(std::make_pair(parser.variables().find("x")->second, std::make_pair(xExpr.first, std::make_pair(rnOneThird, rnTwoThird))));
+	
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> replacedAddSimpleX = function.replaceLeaves(replacementMap);
+	
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddX0 = manager->getEncoding(xExpr.first, 0);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddX1 = manager->getEncoding(xExpr.first, 1);
+
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> complexAdd = 
+          (bddX0.template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(rnTwoThird)))
+        + (bddX1.template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(rnOneThird)));
+	
+	EXPECT_EQ(2ul, replacedAddSimpleX.getNonZeroCount());
+    EXPECT_EQ(2ul, replacedAddSimpleX.getLeafCount());
+    EXPECT_EQ(3ul, replacedAddSimpleX.getNodeCount());
+	EXPECT_TRUE(replacedAddSimpleX == complexAdd);
+}
+
+TEST(SylvanDd, RationalFunctionLeaveReplacementTwoVariables) {
+	std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
+	
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> function;
+	std::shared_ptr<storm::RawPolynomialCache> cache = std::make_shared<storm::RawPolynomialCache>();
+	carl::StringParser parser;
+	parser.setVariables({"x", "y"});
+	
+	storm::RationalFunction variableX(storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("x"), cache));
+	storm::RationalFunction variableY(storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("y"), cache));
+	ASSERT_NO_THROW(function = manager->template getConstant<storm::RationalFunction>(variableX * variableY));
+	
+	std::pair<storm::expressions::Variable, storm::expressions::Variable> xExpr;
+	std::pair<storm::expressions::Variable, storm::expressions::Variable> yExpr;
+	ASSERT_NO_THROW(xExpr = manager->addMetaVariable("x", 0, 1));
+	ASSERT_NO_THROW(yExpr = manager->addMetaVariable("y", 0, 1));
+	
+	std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>> replacementMap;
+	storm::RationalNumber rnOneThird = storm::RationalNumber(1) / storm::RationalNumber(3);
+	storm::RationalNumber rnTwoThird = storm::RationalNumber(2) / storm::RationalNumber(3);
+	storm::RationalNumber rnOne = storm::RationalNumber(1);
+	storm::RationalNumber rnTen = storm::RationalNumber(10);
+	replacementMap.insert(std::make_pair(parser.variables().find("x")->second, std::make_pair(xExpr.first, std::make_pair(rnOneThird, rnTwoThird))));
+	replacementMap.insert(std::make_pair(parser.variables().find("y")->second, std::make_pair(yExpr.first, std::make_pair(rnOne, rnTen))));
+	
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> replacedAdd = function.replaceLeaves(replacementMap);
+	
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddX0 = manager->getEncoding(xExpr.first, 0);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddX1 = manager->getEncoding(xExpr.first, 1);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddY0 = manager->getEncoding(yExpr.first, 0);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddY1 = manager->getEncoding(yExpr.first, 1);
+
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> complexAdd = 
+          ((bddX0 && bddY0).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(rnTwoThird * rnTen)))
+        + ((bddX0 && bddY1).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(rnTwoThird)))
+		+ ((bddX1 && bddY0).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(rnOneThird * rnTen)))
+		+ ((bddX1 && bddY1).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(rnOneThird)));
+	
+	EXPECT_EQ(4ul, replacedAdd.getNonZeroCount());
+    EXPECT_EQ(4ul, replacedAdd.getLeafCount());
+    EXPECT_EQ(7ul, replacedAdd.getNodeCount());
+	EXPECT_TRUE(replacedAdd == complexAdd);
+}
+
+TEST(SylvanDd, RationalFunctionCarlSubstituteTest) {
+	std::shared_ptr<storm::RawPolynomialCache> cache = std::make_shared<storm::RawPolynomialCache>();
+	carl::StringParser parser;
+	parser.setVariables({"x", "y", "z"});
+	storm::RationalFunction zHalfed = storm::RationalFunction(storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("z"), cache), storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("2"), cache));
+	storm::RationalFunction rationalFunction = storm::RationalFunction(storm::Polynomial(parser.parseMultivariatePolynomial<storm::RationalNumber>("2*x+x*y"), cache), storm::Polynomial(parser.parseMultivariatePolynomial<storm::RationalNumber>("1"), cache)) + zHalfed;
+	
+	std::map<storm::RationalFunctionVariable, storm::RationalNumber> replacement = {{parser.variables().find("x")->second, storm::RationalNumber(2)}};
+	storm::RationalFunction subX = rationalFunction.substitute(replacement);
+	
+	storm::RationalFunction cmp = storm::RationalFunction(storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("4+2*y"), cache), storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("1"), cache)) + zHalfed;
+	EXPECT_EQ(subX, cmp);
+	
+	storm::RawPolynomial poly(parser.template parseMultivariatePolynomial<storm::RationalNumber>("2*x+x*y") + parser.template parseMultivariatePolynomial<storm::RationalNumber>("z") / storm::RationalNumber(2));
+	storm::RawPolynomial polySub = poly.substitute(replacement);
+	EXPECT_EQ(polySub, parser.template parseMultivariatePolynomial<storm::RationalNumber>("4+2*y") + parser.template parseMultivariatePolynomial<storm::RationalNumber>("z") / storm::RationalNumber(2));
+}
+
+TEST(SylvanDd, RationalFunctionLeaveReplacementComplexFunction) {
+	std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
+	
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> function;
+	std::shared_ptr<storm::RawPolynomialCache> cache = std::make_shared<storm::RawPolynomialCache>();
+	carl::StringParser parser;
+	parser.setVariables({"x", "y", "z"});
+	storm::RationalFunction zDivTwoY = storm::RationalFunction(storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("z"), cache), storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("2*y"), cache));
+
+	//storm::RationalFunction rationalFunction(two * x + x*y + constantOneDivTwo * z / y);
+	storm::RationalFunction rationalFunction = storm::RationalFunction(storm::Polynomial(parser.parseMultivariatePolynomial<storm::RationalNumber>("2*x+x*y"), cache), storm::Polynomial(parser.parseMultivariatePolynomial<storm::RationalNumber>("1"), cache)) + zDivTwoY;	
+	ASSERT_NO_THROW(function = manager->template getConstant<storm::RationalFunction>(rationalFunction));
+	
+	std::pair<storm::expressions::Variable, storm::expressions::Variable> xExpr;
+	std::pair<storm::expressions::Variable, storm::expressions::Variable> yExpr;
+	std::pair<storm::expressions::Variable, storm::expressions::Variable> zExpr;
+	ASSERT_NO_THROW(xExpr = manager->addMetaVariable("x", 0, 1));
+	ASSERT_NO_THROW(yExpr = manager->addMetaVariable("y", 0, 1));
+	ASSERT_NO_THROW(zExpr = manager->addMetaVariable("z", 0, 1));
+	
+	std::map<storm::RationalFunctionVariable, std::pair<storm::expressions::Variable, std::pair<storm::RationalNumber, storm::RationalNumber>>> replacementMap;
+	storm::RationalNumber rnTwo(2);
+	storm::RationalNumber rnThree(3);
+	storm::RationalNumber rnFive(5);
+	storm::RationalNumber rnSeven(7);
+	storm::RationalNumber rnEleven(11);
+	storm::RationalNumber rnThirteen(13);
+	replacementMap.insert(std::make_pair(parser.variables().find("x")->second, std::make_pair(xExpr.first, std::make_pair(rnTwo, rnSeven))));
+	replacementMap.insert(std::make_pair(parser.variables().find("y")->second, std::make_pair(yExpr.first, std::make_pair(rnThree, rnEleven))));
+	replacementMap.insert(std::make_pair(parser.variables().find("z")->second, std::make_pair(zExpr.first, std::make_pair(rnFive, rnThirteen))));
+	
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> replacedAdd = function.replaceLeaves(replacementMap);
+	
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddX0 = manager->getEncoding(xExpr.first, 0);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddX1 = manager->getEncoding(xExpr.first, 1);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddY0 = manager->getEncoding(yExpr.first, 0);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddY1 = manager->getEncoding(yExpr.first, 1);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddZ0 = manager->getEncoding(zExpr.first, 0);
+	storm::dd::Bdd<storm::dd::DdType::Sylvan> bddZ1 = manager->getEncoding(zExpr.first, 1);
+
+	auto f = [&](bool x, bool y, bool z) {
+		storm::RationalNumber result(2);
+		if (!x) {
+			result *= rnSeven;
+		} else {
+			result *= rnTwo;
+		}
+		
+		storm::RationalNumber partTwo(1);
+		if (!x) {
+			partTwo *= rnSeven;
+		} else {
+			partTwo *= rnTwo;
+		}
+		if (!y) {
+			partTwo *= rnEleven;
+		} else {
+			partTwo *= rnThree;
+		}
+		
+		storm::RationalNumber partThree(1);
+		if (!z) {
+			partThree *= rnThirteen;
+		} else {
+			partThree *= rnFive;
+		}
+		if (!y) {
+			partThree /= storm::RationalNumber(2) * rnEleven;
+		} else {
+			partThree /= storm::RationalNumber(2) * rnThree;
+		}
+		
+		return result + partTwo + partThree;
+	};
+	
+	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> complexAdd = 
+          ((bddX0 && (bddY0 && bddZ0)).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(f(false, false, false))))
+		+ ((bddX0 && (bddY0 && bddZ1)).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(f(false, false, true))))
+		+ ((bddX0 && (bddY1 && bddZ0)).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(f(false, true, false))))
+		+ ((bddX0 && (bddY1 && bddZ1)).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(f(false, true, true))))
+		+ ((bddX1 && (bddY0 && bddZ0)).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(f(true, false, false))))
+		+ ((bddX1 && (bddY0 && bddZ1)).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(f(true, false, true))))
+		+ ((bddX1 && (bddY1 && bddZ0)).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(f(true, true, false))))
+		+ ((bddX1 && (bddY1 && bddZ1)).template toAdd<storm::RationalFunction>() * manager->template getConstant<storm::RationalFunction>(storm::RationalFunction(f(true, true, true))));
+	
+	EXPECT_EQ(8ul, replacedAdd.getNonZeroCount());
+    EXPECT_EQ(8ul, replacedAdd.getLeafCount());
+    EXPECT_EQ(15ul, replacedAdd.getNodeCount());
+	EXPECT_TRUE(replacedAdd == complexAdd);
+	
+	replacedAdd.exportToDot("sylvan_replacedAddC.dot");
+	complexAdd.exportToDot("sylvan_complexAddC.dot");
+}
+
 TEST(SylvanDd, RationalFunctionConstants) {
     std::shared_ptr<storm::dd::DdManager<storm::dd::DdType::Sylvan>> manager(new storm::dd::DdManager<storm::dd::DdType::Sylvan>());
     storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> zero;
@@ -439,33 +647,15 @@ TEST(SylvanDd, RationalFunctionConstants) {
     EXPECT_EQ(1ul, two.getLeafCount());
     EXPECT_EQ(1ul, two.getNodeCount());
 	
-	// The cache that is used in case the underlying type needs a cache.
-	std::shared_ptr<carl::Cache<carl::PolynomialFactorizationPair<storm::RawPolynomial>>> cache = std::make_shared<carl::Cache<carl::PolynomialFactorizationPair<storm::RawPolynomial>>>();
-	
 	storm::dd::Add<storm::dd::DdType::Sylvan, storm::RationalFunction> function;
-	carl::Variable x = carl::freshRealVariable("x");
-	carl::Variable y = carl::freshRealVariable("y");
-	carl::Variable z = carl::freshRealVariable("z");
+	std::shared_ptr<storm::RawPolynomialCache> cache = std::make_shared<storm::RawPolynomialCache>();
+	carl::StringParser parser;
+	parser.setVariables({"x", "y", "z"});
 	
-	storm::RationalFunction constantOne(1);
+	storm::RationalFunction partA = storm::RationalFunction(storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("2*x+x*y"), cache));
+	storm::RationalFunction partB = storm::RationalFunction(storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("z"), cache), storm::Polynomial(parser.template parseMultivariatePolynomial<storm::RationalNumber>("2*y"), cache));
 
-	storm::RationalFunction variableX = storm::RationalFunction(typename storm::RationalFunction::PolyType(typename storm::RationalFunction::PolyType::PolyType(x), cache));
-	storm::RationalFunction variableY = storm::RationalFunction(typename storm::RationalFunction::PolyType(typename storm::RationalFunction::PolyType::PolyType(y), cache));
-	storm::RationalFunction variableZ = storm::RationalFunction(typename storm::RationalFunction::PolyType(typename storm::RationalFunction::PolyType::PolyType(z), cache));
-	
-	storm::RationalFunction constantOneDivTwo(constantOne / constantTwo);
-	storm::RationalFunction tmpFunctionA(constantOneDivTwo);
-	tmpFunctionA *= variableZ;
-	tmpFunctionA /= variableY;
-	storm::RationalFunction tmpFunctionB(variableX);
-	tmpFunctionB *= variableY;
-	
-	
-	//storm::RationalFunction rationalFunction(two * x + x*y + constantOneDivTwo * z / y);
-	storm::RationalFunction rationalFunction(constantTwo);
-	rationalFunction *= variableX;
-	rationalFunction += tmpFunctionB;
-	rationalFunction += tmpFunctionA;
+	storm::RationalFunction rationalFunction = storm::RationalFunction(partA + partB);
 	
     ASSERT_NO_THROW(function = manager->template getConstant<storm::RationalFunction>(rationalFunction));