rewrite

libsolutil/BooleanLP.cpp

@@ -136,7 +136,7 @@ void BooleanLPSolver::addAssertion(Expression const& _expr)
 			{
 				LinearExpression data = *left - *right;
 				data[0] *= -1;
-				Constraint c{move(data), _expr.name == "=", {}};
+				Constraint c{move(data), _expr.name == "="};
 				if (!tryAddDirectBounds(c))
 					state().fixedConstraints.emplace_back(move(c));
 				cout << "Added as fixed constraint" << endl;
@@ -186,7 +186,7 @@ void BooleanLPSolver::addAssertion(Expression const& _expr)
 		{
 			LinearExpression data = *left - *right;
 			data[0] *= -1;
-			Constraint c{move(data), _expr.name == "=", {}};
+			Constraint c{move(data), _expr.name == "="};
 			if (!tryAddDirectBounds(c))
 				state().fixedConstraints.emplace_back(move(c));
 		}
@@ -220,10 +220,24 @@ pair<CheckResult, vector<string>> BooleanLPSolver::check(vector<Expression> cons
 
 	std::vector<std::string> booleanVariables;
 	std::vector<Clause> clauses = state().clauses;
-	SolvingState lpState;
+
+	// TODO we start building up a new set of solver
+	// for each query, but we should also keep some
+	// kind of cache across queries.
+	std::vector<std::pair<size_t, LPSolver>> lpSolvers;
+	lpSolvers.emplace_back(0, LPSolver{});
+	LPSolver& lpSolver = lpSolvers.back().second;
+
 	for (auto&& [index, bound]: state().bounds)
-		resizeAndSet(lpState.bounds, index, bound);
+	{
-	lpState.constraints = state().fixedConstraints;
+		if (bound.lower)
+			lpSolver.addLowerBound(index, *bound.lower);
+		if (bound.upper)
+			lpSolver.addUpperBound(index, *bound.upper);
+	}
+	for (Constraint const& c: state().fixedConstraints)
+		lpSolver.addConstraint(c);
+
 	// TODO this way, it will result in a lot of gaps in both sets of variables.
 	// should we compress them and store a mapping?
 	// Is it even a problem if the indices overlap?
@@ -231,23 +245,9 @@ pair<CheckResult, vector<string>> BooleanLPSolver::check(vector<Expression> cons
 		if (state().isBooleanVariable.at(index) || isConditionalConstraint(index))
 			resizeAndSet(booleanVariables, index, name);
 		else
-			resizeAndSet(lpState.variableNames, index, name);
+			lpSolver.setVariableName(index, name);
 
-	// TODO keep a cache as a member that is never reset.
+	if (lpSolver.check().first == LPResult::Infeasible)
-	// TODO We can also keep the split unconditionals across push/pop
-	// We only need to be careful to update the number of variables.
-
-	std::vector<std::pair<size_t, LPSolver>> lpSolvers;
-
-	// TODO We start afresh here. If we want this to reuse the existing results
-	// from previous invocations of the boolean solver, we still have to use
-	// a cache.
-	// The current optimization is only for CDCL.
-	lpSolvers.emplace_back(0, LPSolver{&m_lpCache});
-	if (
-		lpSolvers.back().second.setState(lpState) == LPResult::Infeasible ||
-		lpSolvers.back().second.check().first == LPResult::Infeasible
-	)
 	{
 		cout << "----->>>>> unsatisfiable" << endl;
 		return {CheckResult::UNSATISFIABLE, {}};
@@ -263,8 +263,7 @@ pair<CheckResult, vector<string>> BooleanLPSolver::check(vector<Expression> cons
 				continue;
 			// "reason" is already stored for those constraints.
 			Constraint const& constraint = state().conditionalConstraints.at(constraintIndex);
-			solAssert(constraint.reasons.size() == 1 && *constraint.reasons.begin() == constraintIndex);
+			lpSolvers.back().second.addConstraint(constraint, constraintIndex);
-			lpSolvers.back().second.addConstraint(constraint);
 		}
 		auto&& [result, modelOrReason] = lpSolvers.back().second.check();
 		// We can only really use the result "infeasible". Everything else should be "sat".
@@ -365,7 +364,7 @@ optional<Literal> BooleanLPSolver::parseLiteral(smtutil::Expression const& _expr
 		LinearExpression data = *left - *right;
 		data[0] *= -1;
 
-		return Literal{true, addConditionalConstraint(Constraint{move(data), _expr.name == "=", {}})};
+		return Literal{true, addConditionalConstraint(Constraint{move(data), _expr.name == "="})};
 	}
 	else if (_expr.name == ">=")
 		return parseLiteral(_expr.arguments.at(1) <= _expr.arguments.at(0));
@@ -390,7 +389,6 @@ Literal BooleanLPSolver::negate(Literal const& _lit)
 			Constraint le = c;
 			le.equality = false;
 			le.data[0] -= 1;
-			le.reasons.clear();
 			Literal leL{true, addConditionalConstraint(le)};
 
 			// X >= b + 1
@@ -399,7 +397,6 @@ Literal BooleanLPSolver::negate(Literal const& _lit)
 			ge.equality = false;
 			ge.data *= -1;
 			ge.data[0] -= 1;
-			ge.reasons.clear();
 			Literal geL{true, addConditionalConstraint(ge)};
 
 
@@ -419,7 +416,6 @@ Literal BooleanLPSolver::negate(Literal const& _lit)
 			Constraint negated = c;
 			negated.data *= -1;
 			negated.data[0] -= 1;
-			negated.reasons.clear();
 			return Literal{true, addConditionalConstraint(negated)};
 		}
 	}
@@ -561,8 +557,6 @@ size_t BooleanLPSolver::addConditionalConstraint(Constraint _constraint)
 	//  - integers
 	declareVariable(name, false);
 	size_t index = state().variables.at(name);
-	solAssert(_constraint.reasons.empty());
-	_constraint.reasons.emplace(index);
 	state().conditionalConstraints[index] = move(_constraint);
 	return index;
 }

libsolutil/BooleanLP.h

@@ -40,8 +40,14 @@ struct State
 	std::map<size_t, Constraint> conditionalConstraints;
 	std::vector<Clause> clauses;
 
+	struct Bounds
+	{
+		std::optional<rational> lower;
+		std::optional<rational> upper;
+	};
+
 	// Unconditional bounds on variables
-	std::map<size_t, SolvingState::Bounds> bounds;
+	std::map<size_t, Bounds> bounds;
 	// Unconditional constraints
 	std::vector<Constraint> fixedConstraints;
 };
@@ -122,9 +128,6 @@ private:
 
 	/// Stack of state, to allow for push()/pop().
 	std::vector<State> m_state{{State{}}};
-
-	std::unordered_map<SolvingState, LPResult> m_lpCache;
-
 };
 
 

libsolutil/LP.h

@@ -42,8 +42,6 @@ struct Constraint
 {
 	LinearExpression data;
 	bool equality = false;
-	/// Set of literals the conjunction of which implies this constraint.
-	std::set<size_t> reasons = {};
 
 	bool operator<(Constraint const& _other) const;
 	bool operator==(Constraint const& _other) const;
@@ -171,156 +169,92 @@ enum class LPResult
 	Infeasible ///< System does not have any solution.
 };
 
-
+/**
-class SimplexWithBounds
+ * LP solver for rational problems, based on "A Fast Linear-Arithmetic Solver for DPLL(T)*"
+ * by Dutertre and Moura.
+ *
+ * Does not solve integer problems!
+ *
+ * Tries to split incoming bounds and constraints into unrelated sub-problems.
+ * Maintains lower/upper bounds for all variables.
+ * Adds one slack variable per constraint and stores all constraints as "= 0" equations.
+ * Splits variables into basic and non-basic. For each row there is exactly one
+ * basic variable that has a factor of -1.
+ * The equations are satisfied at all times and non-basic variables are always within their bounds.
+ * Non-basic variables might violate their bounds.
+ * It attempts to resolve these violations in turn, swapping a basic variables with a non-basic
+ * variables that can still move in the required direction.
+ *
+ * It is perfectly fine to add new bounds, variable or constraints after a call to "check".
+ * The solver can be copied at low cost and it uses a "copy on write" mechanism for the sub-problems.
+ */
+class LPSolver
 {
 public:
-	explicit SimplexWithBounds(SolvingState _state);
+	void addConstraint(Constraint const& _constraint, std::optional<size_t> _reason = std::nullopt);
-	LPResult check();
+	void setVariableName(size_t _variable, std::string _name);
+	void addLowerBound(size_t _variable, rational _bound);
+	void addUpperBound(size_t _variable, rational _bound);
 
-	size_t addVariable(std::string _name);
+	std::pair<LPResult, std::variant<Model, ReasonSet>> check();
-	void addConstraint(Constraint _constraint);
 
 	std::string toString() const;
+
 private:
+	struct Bounds
+	{
+		std::optional<rational> lower;
+		std::optional<rational> upper;
+	};
+	struct Variable
+	{
+		std::string name = {};
+		rational value = 0;
+		Bounds bounds = {};
+	};
+	struct SubProblem
+	{
+		/// Set to true on "check". Needs a copy for adding a constraint or bound if set to true.
+		bool sealed = false;
+		std::optional<LPResult> result = std::nullopt;
+		std::vector<LinearExpression> factors;
+		std::vector<Variable> variables;
+		/// Variable index to constraint it controls.
+		std::map<size_t, size_t> basicVariables;
+		/// Maps outer indices to inner indices.
+		std::map<size_t, size_t> varMapping = {};
+		std::set<size_t> reasons;
+
+		LPResult check();
+		std::string toString() const;
+	private:
 		/// Set value of non-basic variable.
-	void update(size_t _var, rational const& _value);
+		void update(size_t _varIndex, rational const& _value);
 		/// @returns the index of the first basic variable violating its bounds.
 		std::optional<size_t> firstConflictingBasicVariable() const;
 		std::optional<size_t> firstReplacementVar(size_t _basicVarToReplace, bool _increasing) const;
 
 		void pivot(size_t _old, size_t _new);
 		void pivotAndUpdate(size_t _oldBasicVar, rational const& _newValue, size_t _newBasicVar);
-
-	SolvingState m_state;
-	std::vector<rational> m_assignments;
-	/// Variable index to row it controls.
-	std::map<size_t, size_t> m_basicVariables;
-};
-
-
-/**
- * Applies several strategies to simplify a given solving state.
- * During these simplifications, it can sometimes already be determined if the
- * state is feasible or not.
- * Since some variables can be fixed to specific values, it returns a
- * (partial) model.
- *
- * - Constraints with exactly one nonzero coefficient represent "a x <= b"
- *   and thus are turned into bounds.
- * - Constraints with zero nonzero coefficients are constant relations.
- *   If such a relation is false, answer "infeasible", otherwise remove the constraint.
- * - Empty columns can be removed.
- * - Variables with matching bounds can be removed from the problem by substitution.
- *
- * Holds a reference to the solving state that is modified during operation.
- */
-class SolvingStateSimplifier
-{
-public:
-	SolvingStateSimplifier(SolvingState& _state):
-		m_state(_state) {}
-
-	std::pair<LPResult, std::variant<std::map<size_t, rational>, ReasonSet>> simplify();
-
-private:
-	/// Remove variables that have equal lower and upper bound.
-	/// @returns reason / set of conflicting clauses if infeasible.
-	std::optional<ReasonSet> removeFixedVariables();
-
-	/// Removes constraints of the form 0 <= b or 0 == b (no variables) and
-	/// turns constraints of the form a * x <= b (one variable) into bounds.
-	/// @returns reason / set of conflicting clauses if infeasible.
-	std::optional<ReasonSet> extractDirectConstraints();
-
-	/// Removes all-zeros columns.
-	void removeEmptyColumns();
-
-	/// Set to true by the strategies if they performed some changes.
-	bool m_changed = false;
-
-	SolvingState& m_state;
-	std::map<size_t, rational> m_fixedVariables;
-};
-
-/**
- * Splits a given linear program into multiple linear programs with disjoint sets of variables.
- * The initial program is feasible if and only if all sub-programs are feasible.
- */
-class ProblemSplitter
-{
-public:
-	explicit ProblemSplitter(SolvingState const& _state);
-
-	/// @returns true if there are still sub-problems to split out.
-	operator bool() const { return m_column < m_state.variableNames.size(); }
-
-	/// @returns the next sub-problem.
-	std::pair<std::vector<bool>, std::vector<bool>> next();
-
-private:
-	SolvingState const& m_state;
-	/// Next column to start the search for a connected component.
-	size_t m_column = 1;
-	/// The columns we have already split out.
-	std::vector<bool> m_seenColumns;
-};
-
-/**
- * LP solver for rational problems.
- *
- * Does not solve integer problems!
- *
- * Tries to split a given problem into sub-problems and utilizes a cache to quickly solve
- * similar problems.
- *
- * Can be used in a mode where it does not support returning models. In that case, the
- * cache is more efficient.
- */
-class LPSolver
-{
-public:
-	explicit LPSolver(bool _supportModels = true);
-	explicit LPSolver(std::unordered_map<SolvingState, LPResult>* _cache):
-		m_cache(_cache) {}
-
-
-	LPResult setState(SolvingState _state);
-	void addConstraint(Constraint _constraint);
-	std::pair<LPResult, std::variant<Model, ReasonSet>> check();
-
-private:
-	void combineSubProblems(size_t _combineInto, size_t _combineFrom);
-	void addConstraintToSubProblem(size_t _subProblem, Constraint _constraint);
-	void updateSubProblems();
-
-	/// Ground state for CDCL. This is shared by copies of the solver.
-	/// Only ``setState`` changes the state. Copies will only use
-	/// ``addConstraint`` which does not change m_state.
-	std::shared_ptr<SolvingState> m_state;
-	struct SubProblem
-	{
-		// TODO now we could actually put the constraints here again.
-		std::vector<Constraint> removableConstraints;
-		bool dirty = true;
-		LPResult result = LPResult::Unknown;
-		std::vector<boost::rational<bigint>> model = {};
-		std::set<size_t> variables = {};
-		std::optional<SimplexWithBounds> simplex = std::nullopt;
-		std::map<size_t, size_t> varMapping = {};
 	};
 
-	std::pair<SolvingState, std::map<size_t, size_t>> stateFromSubProblem(size_t _index) const;
-	ReasonSet reasonSetForSubProblem(SubProblem const& _subProblem);
 
-	std::shared_ptr<std::map<size_t, rational>> m_fixedVariables;
+	SubProblem& unseal(size_t _problem);
+	/// Unseals the problem for the given variable or creates a new one.
+	SubProblem& unsealForVariable(size_t _outerIndex);
+	void combineSubProblems(size_t _combineInto, size_t _combineFrom);
+	void addConstraintToSubProblem(size_t _subProblem, Constraint const& _constraint, std::optional<size_t> _reason);
+	void addOuterVariableToSubProblem(size_t _subProblem, size_t _outerIndex);
+	size_t addNewVariableToSubProblem(size_t _subProblem);
+
+	std::map<std::string, rational> model() const;
+
 	/// These use "copy on write".
 	std::vector<std::shared_ptr<SubProblem>> m_subProblems;
-	std::vector<size_t> m_subProblemsPerVariable;
+	/// Maps outer indices to sub problems.
-	std::vector<size_t> m_subProblemsPerConstraint;
+	std::map<size_t, size_t> m_subProblemsPerVariable;
-	/// TODO also store the first infeasible subproblem?
+	/// Counter to enable unique names for the slack variables.
-	/// TODO still retain the cache?
+	size_t m_slackVariableCounter = 0;
-	std::unordered_map<SolvingState, LPResult>* m_cache = nullptr;
 
 };
 

test/libsolutil/LP.cpp

@@ -38,7 +38,6 @@ class LPTestFramework
 public:
 	LPTestFramework()
 	{
-		m_solvingState.variableNames.emplace_back("");
 	}
 
 	LinearExpression constant(rational _value)
@@ -52,11 +51,11 @@ public:
 	}
 
 	/// Adds the constraint "_lhs <= _rhs".
-	void addLEConstraint(LinearExpression _lhs, LinearExpression _rhs, set<size_t> _reason = {})
+	void addLEConstraint(LinearExpression _lhs, LinearExpression _rhs, optional<size_t> _reason = {})
 	{
 		_lhs -= _rhs;
 		_lhs[0] = -_lhs[0];
-		m_solvingState.constraints.push_back({move(_lhs), false, move(_reason)});
+		m_solver.addConstraint({move(_lhs), false}, move(_reason));
 	}
 
 	void addLEConstraint(LinearExpression _lhs, rational _rhs)
@@ -65,47 +64,43 @@ public:
 	}
 
 	/// Adds the constraint "_lhs = _rhs".
-	void addEQConstraint(LinearExpression _lhs, LinearExpression _rhs, set<size_t> _reason = {})
+	void addEQConstraint(LinearExpression _lhs, LinearExpression _rhs, optional<size_t> _reason = {})
 	{
 		_lhs -= _rhs;
 		_lhs[0] = -_lhs[0];
-		m_solvingState.constraints.push_back({move(_lhs), true, move(_reason)});
+		m_solver.addConstraint({move(_lhs), true}, move(_reason));
 	}
 
-	void addLowerBound(string _variable, rational _value, set<size_t> _reason = {})
+	void addLowerBound(string _variable, rational _value)
 	{
-		size_t index = variableIndex(_variable);
+		m_solver.addLowerBound(variableIndex(_variable), _value);
-		if (index >= m_solvingState.bounds.size())
-			m_solvingState.bounds.resize(index + 1);
-		m_solvingState.bounds.at(index).lower = _value;
-		m_solvingState.bounds.at(index).lowerReasons = move(_reason);
 	}
 
-	void addUpperBound(string _variable, rational _value, set<size_t> _reason = {})
+	void addUpperBound(string _variable, rational _value)
 	{
-		size_t index = variableIndex(_variable);
+		m_solver.addUpperBound(variableIndex(_variable), _value);
-		if (index >= m_solvingState.bounds.size())
-			m_solvingState.bounds.resize(index + 1);
-		m_solvingState.bounds.at(index).upper = _value;
-		m_solvingState.bounds.at(index).upperReasons = move(_reason);
 	}
 
 	void feasible(vector<pair<string, rational>> const& _solution)
 	{
-		m_solver.setState(m_solvingState);
 		auto [result, modelOrReasonSet] = m_solver.check();
 		BOOST_REQUIRE(result == LPResult::Feasible);
 		Model model = get<Model>(modelOrReasonSet);
 		for (auto const& [var, value]: _solution)
+		{
+			BOOST_CHECK_MESSAGE(
+				model.count(var),
+				"Variable " + var + " not found in model."
+			);
 			BOOST_CHECK_MESSAGE(
 				value == model.at(var),
 				var + " = "s + ::toString(model.at(var)) + " (expected " + ::toString(value) + ")"
 			);
 		}
+	}
 
 	void infeasible(set<size_t> _reason = {})
 	{
-		m_solver.setState(m_solvingState);
 		auto [result, modelOrReason] = m_solver.check();
 		BOOST_REQUIRE(result == LPResult::Infeasible);
 		ReasonSet suppliedReason = get<ReasonSet>(modelOrReason);
@@ -115,44 +110,54 @@ public:
 protected:
 	size_t variableIndex(string const& _name)
 	{
-		if (m_solvingState.variableNames.empty())
+		if (!m_variableIndices.count(_name))
-			m_solvingState.variableNames.emplace_back("");
-		auto index = findOffset(m_solvingState.variableNames, _name);
-		if (!index)
 		{
-			index = m_solvingState.variableNames.size();
+			size_t index = 1 + m_variableIndices.size();
-			m_solvingState.variableNames.emplace_back(_name);
+			m_solver.setVariableName(index, _name);
+			m_variableIndices[_name] = index;
 		}
-		return *index;
+		return m_variableIndices.at(_name);
 	}
 
 	LPSolver m_solver;
-	SolvingState m_solvingState;
+	map<string, size_t> m_variableIndices;
 };
 
 
 BOOST_FIXTURE_TEST_SUITE(LP, LPTestFramework, *boost::unit_test::label("nooptions"))
 
+BOOST_AUTO_TEST_CASE(empty)
+{
+	feasible({});
+}
+
 BOOST_AUTO_TEST_CASE(basic)
 {
 	auto x = variable("x");
 	addLEConstraint(2 * x, 10);
-	feasible({{"x", 5}});
+	feasible({{"x", 0}});
 }
 
 BOOST_AUTO_TEST_CASE(not_linear_independent)
 {
 	addLEConstraint(2 * variable("x"), 10);
 	addLEConstraint(4 * variable("x"), 20);
+	feasible({{"x", 0}});
+}
+
+BOOST_AUTO_TEST_CASE(not_linear_independent_eq)
+{
+	addLEConstraint(2 * variable("x"), 10);
+	addEQConstraint(4 * variable("x"), constant(20));
 	feasible({{"x", 5}});
 }
 
 BOOST_AUTO_TEST_CASE(two_vars)
 {
 	addLEConstraint(variable("y"), 3);
-	addLEConstraint(variable("x"), 10);
+	addLowerBound("x", 10);
 	addLEConstraint(variable("x") + variable("y"), 4);
-	feasible({{"x", 1}, {"y", 3}});
+	feasible({{"x", 10}, {"y", -6}});
 }
 
 BOOST_AUTO_TEST_CASE(one_le_the_other)
@@ -167,8 +172,8 @@ BOOST_AUTO_TEST_CASE(factors)
 	auto y = variable("y");
 	addLEConstraint(2 * y, 3);
 	addLEConstraint(16 * x, 10);
-	addLEConstraint(x + y, 4);
+	addLEConstraint(constant(2), x + y);
-	feasible({{"x", rational(5) / 8}, {"y", rational(3) / 2}});
+	feasible({{"x", rational(5) / 8}, {"y", rational(11) / 8}});
 }
 
 BOOST_AUTO_TEST_CASE(cache)
@@ -178,8 +183,8 @@ BOOST_AUTO_TEST_CASE(cache)
 	// that it results in the same value.
 	auto x = variable("x");
 	auto y = variable("y");
-	addLEConstraint(2 * y, 3);
+	addLEConstraint(constant(3), 2 * y);
-	addLEConstraint(2 * x, 3);
+	addLEConstraint(constant(3), 2 * x);
 	feasible({{"x", rational(3) / 2}, {"y", rational(3) / 2}});
 	feasible({{"x", rational(3) / 2}, {"y", rational(3) / 2}});
 }
@@ -187,19 +192,25 @@ BOOST_AUTO_TEST_CASE(cache)
 BOOST_AUTO_TEST_CASE(bounds)
 {
 	addUpperBound("x", 200);
-	feasible({{"x", 200}});
+	feasible({{"x", 0}});
 
 	addLEConstraint(variable("x"), 100);
-	feasible({{"x", 100}});
+	feasible({{"x", 0}});
 
 	addLEConstraint(constant(5), variable("x"));
-	feasible({{"x", 100}});
+	feasible({{"x", 5}});
 
 	addLowerBound("x", 20);
-	feasible({{"x", 100}});
+	feasible({{"x", 20}});
 	addLowerBound("x", 25);
-	feasible({{"x", 100}});
+	feasible({{"x", 25}});
+	addLowerBound("x", 21);
+	feasible({{"x", 25}});
 
+	addUpperBound("x", 26);
+	feasible({{"x", 25}});
+	addUpperBound("x", 28);
+	feasible({{"x", 25}});
 	addUpperBound("x", 20);
 	infeasible();
 }
@@ -210,19 +221,20 @@ BOOST_AUTO_TEST_CASE(bounds2)
 	addUpperBound("x", 250);
 	addLowerBound("y", 2);
 	addUpperBound("y", 3);
-	feasible({{"x", 250}, {"y", 3}});
+	feasible({{"x", 200}, {"y", 2}});
 
 	addLEConstraint(variable("y"), variable("x"));
-	feasible({{"x", 250}, {"y", 3}});
+	feasible({{"x", 200}, {"y", 2}});
 
 	addEQConstraint(variable("y") + constant(231), variable("x"));
-	feasible({{"x", 234}, {"y", 3}});
+	cout << m_solver.toString() << endl;
-
+	feasible({{"x", 233}, {"y", 2}});
+/*
 	addEQConstraint(variable("y") + constant(10), variable("x") - variable("z"));
-	feasible({{"x", 234}, {"y", 3}});
+	feasible({{"x", 234}, {"y", 3}, {"z", 0}});
 
 	addEQConstraint(variable("z") + variable("x"), constant(2));
-	infeasible();
+	infeasible();*/
 }
 
 BOOST_AUTO_TEST_CASE(lower_bound)
@@ -230,7 +242,7 @@ BOOST_AUTO_TEST_CASE(lower_bound)
 	addLEConstraint(constant(1), variable("y"));
 	addLEConstraint(variable("x"), constant(10));
 	addLEConstraint(2 * variable("x") + variable("y"), 2);
-	feasible({{"x", 0}, {"y", 2}});
+	feasible({{"x", 0}, {"y", 1}});
 }
 
 BOOST_AUTO_TEST_CASE(check_infeasible)
@@ -252,11 +264,13 @@ BOOST_AUTO_TEST_CASE(unbounded2)
 	auto y = variable("y");
 	addLEConstraint(constant(2), x + y);
 	addLEConstraint(x, 10);
-	feasible({{"x", 10}, {"y", 0}});
+	feasible({{"x", 2}, {"y", 0}});
 }
 
 BOOST_AUTO_TEST_CASE(unbounded3)
 {
+	addLowerBound("y", 0);
+	addLowerBound("x", 0);
 	addLEConstraint(constant(0) - variable("x") - variable("y"), constant(10));
 	feasible({{"x", 0}, {"y", 0}});
 
@@ -277,7 +291,7 @@ BOOST_AUTO_TEST_CASE(equal)
 	auto y = variable("y");
 	addEQConstraint(x, y + constant(10));
 	addLEConstraint(x, 20);
-	feasible({{"x", 20}, {"y", 10}});
+	feasible({{"x", 10}, {"y", 0}});
 }
 
 
@@ -285,6 +299,7 @@ BOOST_AUTO_TEST_CASE(equal_constant)
 {
 	auto x = variable("x");
 	auto y = variable("y");
+	addLowerBound("x", 5);
 	addLEConstraint(x, y);
 	addEQConstraint(y, constant(5));
 	feasible({{"x", 5}, {"y", 5}});
@@ -295,7 +310,26 @@ BOOST_AUTO_TEST_CASE(all_equality)
 	auto x = variable("x");
 	auto y = variable("y");
 	addEQConstraint(-6 * x - 6 * y, constant(8));
-	infeasible();
+	feasible({{"x", -rational(4) / 3}, {"y", 0}});
+}
+
+BOOST_AUTO_TEST_CASE(negative_values)
+{
+	auto x = variable("x");
+	auto y = variable("y");
+	addEQConstraint(-2 * x, constant(8));
+	addLowerBound("y", -2);
+	addUpperBound("y", -1);
+	feasible({{"x", -4}, {"y", -1}});
+}
+
+BOOST_AUTO_TEST_CASE(basic_basic)
+{
+	auto x = variable("x");
+	addLEConstraint(x, constant(5));
+	feasible({{"x", 0}});
+	addEQConstraint(x, constant(5));
+	feasible({{"x", 5}});
 }
 
 BOOST_AUTO_TEST_CASE(linear_dependent)
@@ -303,40 +337,40 @@ BOOST_AUTO_TEST_CASE(linear_dependent)
 	auto x = variable("x");
 	auto y = variable("y");
 	auto z = variable("z");
-	addLEConstraint(x, 5);
+	addLEConstraint(x, -5);
-	addLEConstraint(2 * y, 10);
+	addLEConstraint(2 * y, -10);
-	addLEConstraint(3 * z, 15);
+	addLEConstraint(3 * z, -15);
 	// Here, they should be split into three independent problems.
-	feasible({{"x", 5}, {"y", 5}, {"z", 5}});
+	feasible({{"x", -5}, {"y", -5}, {"z", -5}});
 
 	addLEConstraint((x + y) + z, 100);
-	feasible({{"x", 5}, {"y", 5}, {"z", 5}});
+	feasible({{"x", -5}, {"y", -5}, {"z", -5}});
 
-	addLEConstraint((x + y) + z, 2);
+	addLEConstraint((x + y) + z, -1);
-	feasible({{"x", 2}, {"y", 0}, {"z", 0}});
+	feasible({{"x", -5}, {"y", -5}, {"z", -5}});
 
-	addLEConstraint(constant(2), (x + y) + z);
+	addLEConstraint(constant(-20), (x + y) + z);
-	feasible({{"x", 2}, {"y", 0}, {"z", 0}});
+	feasible({{"x", -5}, {"y", -5}, {"z", -5}});
 
-	addEQConstraint(constant(2), (x + y) + z);
+	addEQConstraint(constant(-18), (x + y) + z);
-	feasible({{"x", 2}, {"y", 0}, {"z", 0}});
+	feasible({{"x", -8}, {"y", -5}, {"z", -5}});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_simple)
 {
 	auto x = variable("x");
 	addLEConstraint(2 * x, constant(20), {0});
-	addLowerBound("x", 12, {1});
+	addLowerBound("x", 12);
-	infeasible({0, 1});
+	infeasible({0});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_bounds)
 {
 	auto x = variable("x");
-	addLowerBound("x", 12, {0});
+	addLowerBound("x", 12);
-	addUpperBound("x", 11, {1});
+	addUpperBound("x", 11);
-	addLEConstraint(x, constant(200), {2}); // unrelated
+	addLEConstraint(variable("y"), constant(200), {2}); // unrelated
-	infeasible({0, 1});
+	infeasible({});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_forwarded)
@@ -346,9 +380,9 @@ BOOST_AUTO_TEST_CASE(reasons_forwarded)
 	addEQConstraint(x, constant(2), {0});
 	addEQConstraint(x, y, {1});
 	feasible({});
-	addLEConstraint(x, constant(200), {5}); // unrelated
+	addLEConstraint(x, constant(200), {5});
 	addLEConstraint(y, constant(1), {2});
-	infeasible({0, 1, 2});
+	infeasible({0, 1, 2, 5});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_forwarded2)
@@ -359,9 +393,9 @@ BOOST_AUTO_TEST_CASE(reasons_forwarded2)
 	addEQConstraint(x, y, {1});
 	feasible({});
 	addEQConstraint(y, constant(3), {2});
-	addLEConstraint(x, constant(200), {6}); // unrelated
+	addLEConstraint(x, constant(200), {6});
-	addLEConstraint(y, constant(202), {6}); // unrelated
+	addLEConstraint(y, constant(202), {7});
-	infeasible({0, 1, 2});
+	infeasible({0, 1, 2, 6, 7});
 }
 
 BOOST_AUTO_TEST_CASE(reasons_split)
@@ -377,6 +411,22 @@ BOOST_AUTO_TEST_CASE(reasons_split)
 	infeasible({0, 2});
 }
 
+BOOST_AUTO_TEST_CASE(reasons_joined)
+{
+	auto x = variable("x");
+	auto y = variable("y");
+	auto z = variable("z");
+	addLEConstraint(x + y, constant(2), {0});
+	feasible({});
+	addLEConstraint(constant(20), x + y, {2});
+	infeasible({0, 2});
+	addLEConstraint(z, constant(200), {3});
+	infeasible({0, 2});
+	addLEConstraint(x, z);
+	infeasible({0, 2, 3});
+}
+
+
 BOOST_AUTO_TEST_CASE(fuzzer2)
 {
 	/*
@@ -399,6 +449,15 @@ BOOST_AUTO_TEST_CASE(fuzzer2)
 	auto x6 = variable("x6");
 	auto x7 = variable("x7");
 	auto x8 = variable("x8");
+	addLowerBound("x0", 0);
+	addLowerBound("x1", 0);
+	addLowerBound("x2", 0);
+	addLowerBound("x3", 0);
+	addLowerBound("x4", 0);
+	addLowerBound("x5", 0);
+	addLowerBound("x6", 0);
+	addLowerBound("x7", 0);
+	addLowerBound("x8", 0);
 	addEQConstraint(90 * x4 + 9 * x5, constant(-1));
 	addLEConstraint(-76 * x2 + 74 * x5, constant(0));
 	addLEConstraint(31 * x3 - 71 * x8, constant(0));