Repository: plaans/aries (language: python)

Issue:
Fix reification of constraints with optional variables
Expressions with optional variables should be reified to optional literals

For instance a constraint `a < b` should be reified to a literal `l` such that

- `prez(l) <=> prez(a) && prez(b)`
- `l = true` iff `prez(a) && prez(b) & int(a) < int(b)`
- `l = false` iff `prez(a) && prez(b) & int(a) >= int(b)`
- `l = absent` iff `!prez(a) || !prez(b)`

Currently the value of `l` is unspecified if one of the variables is absent which is unsatisfactory and requires the addition of specialized constraints such as `OptLeq`


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml
index 75325ffd..59ca5fd5 100644
--- a/.github/workflows/rust.yml
+++ b/.github/workflows/rust.yml
@@ -39,8 +39,13 @@ jobs:
   solving:
     name: Solving
     runs-on: ubuntu-20.04
-    env:
-      RUST_BACKTRACE: 1 # have backtraces in the CI environment
+    env: # Make sure we have fast binaries and debuggable binaries
+      CARGO_PROFILE_RELEASE_DEBUG: false
+      CARGO_PROFILE_DEV_DEBUG: false
+      CARGO_PROFILE_RELEASE_DEBUG_ASSERTIONS: true
+      CARGO_PROFILE_DEV_OPT_LEVEL: 3
+      CARGO_PROFILE_DEV_LTO: thin
+      RUST_BACKTRACE: false
     steps:
       - uses: actions/checkout@v2
       - uses: actions-rs/toolchain@v1
@@ -53,7 +58,7 @@ jobs:
       - name: Get problems from LFS
         run: git lfs pull
       - name: SAT solving
-        run: ./ci/sat.py
+        run: ./ci/sat.py debug
       - name: Jobshop solving
         run: ./ci/jobshop.py
       - name: GG solving
diff --git a/Cargo.lock b/Cargo.lock
index 4a63fc77..a171689a 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -47,6 +47,18 @@ dependencies = [
  "vec_map",
 ]
 
+[[package]]
+name = "aries_core"
+version = "0.1.0"
+dependencies = [
+ "anyhow",
+ "aries_backtrack",
+ "aries_collections",
+ "aries_utils",
+ "lazy_static",
+ "streaming-iterator",
+]
+
 [[package]]
 name = "aries_gg"
 version = "0.1.0"
@@ -63,6 +75,7 @@ version = "0.1.0"
 dependencies = [
  "anyhow",
  "aries_backtrack",
+ "aries_core",
  "aries_model",
  "aries_solver",
  "aries_tnet",
@@ -76,10 +89,9 @@ dependencies = [
  "anyhow",
  "aries_backtrack",
  "aries_collections",
+ "aries_core",
  "aries_utils",
- "criterion",
  "lazy_static",
- "rand",
  "streaming-iterator",
 ]
 
@@ -90,6 +102,7 @@ dependencies = [
  "anyhow",
  "aries_backtrack",
  "aries_collections",
+ "aries_core",
  "aries_model",
  "aries_planning",
  "aries_solver",
@@ -105,6 +118,7 @@ version = "0.1.0"
 dependencies = [
  "anyhow",
  "aries_collections",
+ "aries_core",
  "aries_model",
  "aries_utils",
  "env_param",
@@ -119,6 +133,7 @@ name = "aries_sat"
 version = "0.1.0"
 dependencies = [
  "anyhow",
+ "aries_core",
  "aries_model",
  "aries_solver",
  "structopt",
@@ -132,9 +147,11 @@ name = "aries_smt"
 version = "0.1.0"
 dependencies = [
  "aries_backtrack",
+ "aries_core",
  "aries_model",
  "aries_solver",
  "aries_tnet",
+ "itertools",
 ]
 
 [[package]]
@@ -143,6 +160,7 @@ version = "0.1.0"
 dependencies = [
  "aries_backtrack",
  "aries_collections",
+ "aries_core",
  "aries_model",
  "env_param",
  "itertools",
@@ -156,6 +174,7 @@ version = "0.1.0"
 dependencies = [
  "aries_backtrack",
  "aries_collections",
+ "aries_core",
  "aries_model",
  "aries_solver",
  "criterion",
diff --git a/Cargo.toml b/Cargo.toml
index 6fa00505..9050d8a1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,7 +1,7 @@
 [workspace]
 
 members = [
-    "planning", "planners", "tnet", "collections", "model", "solver", "env_param", "utils", "examples/sat", "examples/jobshop", "examples/smt", "examples/gg"
+    "core", "planning", "planners", "tnet", "collections", "model", "solver", "env_param", "utils", "examples/sat", "examples/jobshop", "examples/smt", "examples/gg"
 ]
 
 [profile.dev]
diff --git a/core/Cargo.toml b/core/Cargo.toml
new file mode 100644
index 00000000..284b402e
--- /dev/null
+++ b/core/Cargo.toml
@@ -0,0 +1,34 @@
+[package]
+name = "aries_core"
+version = "0.1.0"
+authors = ["Arthur Bit-Monnot <abitmonnot@laas.fr>"]
+edition = "2018"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
+anyhow = "1.0.35"
+streaming-iterator = "0.1.5"
+aries_backtrack = { path = "../backtrack" }
+aries_collections = { path = "../collections" }
+aries_utils = { path = "../utils" }
+lazy_static = "1.4.0"
+
+# TODO: reactivate benchmarks
+#[dev-dependencies]
+#criterion = "0.3"
+#rand = "0.8.3"
+#
+#
+#[[bench]]
+#name = "bounds"
+#harness = false
+#
+#[[bench]]
+#name = "watches"
+#harness = false
+#
+#
+#[[bench]]
+#name = "domains"
+#harness = false
\ No newline at end of file
diff --git a/model/benches/bounds.rs b/core/benches/bounds.rs
similarity index 78%
rename from model/benches/bounds.rs
rename to core/benches/bounds.rs
index 7919ed78..4901b123 100644
--- a/model/benches/bounds.rs
+++ b/core/benches/bounds.rs
@@ -1,5 +1,5 @@
-use aries_model::bounds::Lit;
-use aries_model::Model;
+use aries_core::state::Domains;
+use aries_core::*;
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
 use rand::prelude::SliceRandom;
 use rand::rngs::StdRng;
@@ -21,10 +21,10 @@ fn entailment(xs: &[Lit], ys: &[Lit]) -> u64 {
 pub fn criterion_benchmark(c: &mut Criterion) {
     let mut rng = StdRng::seed_from_u64(2398248538438434234);
 
-    let mut model = Model::new();
+    let mut model = Domains::new();
     let mut bounds = Vec::new();
     for _ in 0..50 {
-        let var = model.new_ivar(0, 100, "");
+        let var = model.new_var(0, 100);
         for v in -20..20 {
             bounds.push(Lit::leq(var, v));
             bounds.push(Lit::geq(var, v));
@@ -33,13 +33,13 @@ pub fn criterion_benchmark(c: &mut Criterion) {
 
     bounds.shuffle(&mut rng);
 
-    c.bench_function("bounds-entail-many-vars", |b| {
+    c.bench_function("literals-entail-many-vars", |b| {
         b.iter(|| entailment(black_box(&bounds), black_box(&bounds)))
     });
 
     let mut bounds = Vec::new();
     for _ in 0..5 {
-        let var = model.new_ivar(0, 100, "");
+        let var = model.new_var(0, 100);
         for v in -20..20 {
             bounds.push(Lit::leq(var, v));
             bounds.push(Lit::geq(var, v));
@@ -48,13 +48,13 @@ pub fn criterion_benchmark(c: &mut Criterion) {
 
     bounds.shuffle(&mut rng);
 
-    c.bench_function("bounds-entail-few-vars", |b| {
+    c.bench_function("literals-entail-few-vars", |b| {
         b.iter(|| entailment(black_box(&bounds), black_box(&bounds)))
     });
 
     let mut bounds = Vec::new();
 
-    let var = model.new_ivar(0, 100, "");
+    let var = model.new_var(0, 100);
     for v in -40..40 {
         bounds.push(Lit::leq(var, v));
         bounds.push(Lit::geq(var, v));
@@ -62,7 +62,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
 
     bounds.shuffle(&mut rng);
 
-    c.bench_function("bounds-entail-one-var", |b| {
+    c.bench_function("literals-entail-one-var", |b| {
         b.iter(|| entailment(black_box(&bounds), black_box(&bounds)))
     });
 }
diff --git a/model/benches/domains.rs b/core/benches/domains.rs
similarity index 92%
rename from model/benches/domains.rs
rename to core/benches/domains.rs
index a8d44ba5..182017c7 100644
--- a/model/benches/domains.rs
+++ b/core/benches/domains.rs
@@ -1,13 +1,15 @@
-use aries_backtrack::Backtrack;
-use aries_model::bounds::Lit;
-use aries_model::extensions::assignments::AssignmentExt;
-use aries_model::lang::IVar;
-use aries_model::state::domains::Domains;
-use aries_model::{Model, WriterId};
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
 use rand::prelude::{SliceRandom, StdRng};
 use rand::{Rng, SeedableRng};
 
+use aries_backtrack::Backtrack;
+use aries_core::literals::Lit;
+use aries_core::state::domains::Domains;
+use aries_core::WriterId;
+use aries_model::extensions::assignments::AssignmentExt;
+use aries_model::lang::IVar;
+use aries_model::Model;
+
 fn count_entailed(xs: &[Lit], model: &Model) -> usize {
     let mut count = 0;
     for &x in xs {
@@ -60,14 +62,14 @@ pub fn read_benchmark(c: &mut Criterion) {
     literals.shuffle(&mut rng);
     variables.shuffle(&mut rng);
 
-    c.bench_function("model-count-entailed-bounds", |b| {
+    c.bench_function("model-count-entailed-literals", |b| {
         b.iter(|| count_entailed(black_box(&literals), black_box(&model)))
     });
 
-    c.bench_function("model-count-lower-bounds", |b| {
+    c.bench_function("model-count-lower-literals", |b| {
         b.iter(|| count_lower_bounds(black_box(&variables), black_box(&model)))
     });
-    c.bench_function("model-count-upper-bounds", |b| {
+    c.bench_function("model-count-upper-literals", |b| {
         b.iter(|| count_upper_bounds(black_box(&variables), black_box(&model)))
     });
 }
diff --git a/model/benches/watches.rs b/core/benches/watches.rs
similarity index 96%
rename from model/benches/watches.rs
rename to core/benches/watches.rs
index fea7c5b5..8383a250 100644
--- a/model/benches/watches.rs
+++ b/core/benches/watches.rs
@@ -1,4 +1,4 @@
-use aries_model::bounds::{Lit, Watches};
+use aries_core::literals::{Lit, Watches};
 use aries_model::Model;
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
 use rand::prelude::{SliceRandom, StdRng};
@@ -43,7 +43,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
         b.iter(|| count_watches(black_box(&bounds), black_box(&watches)))
     });
 
-    // shuffle bounds
+    // shuffle literals
 
     bounds.shuffle(&mut rng);
 
diff --git a/model/src/bounds/bound_value.rs b/core/src/bound_value.rs
similarity index 93%
rename from model/src/bounds/bound_value.rs
rename to core/src/bound_value.rs
index d256706b..74c547e6 100644
--- a/model/src/bounds/bound_value.rs
+++ b/core/src/bound_value.rs
@@ -1,4 +1,4 @@
-use crate::lang::IntCst;
+use crate::*;
 
 #[derive(Copy, Clone, Debug, Eq, PartialEq, Ord, PartialOrd, Hash)]
 pub struct BoundValue(i32);
@@ -23,7 +23,7 @@ impl BoundValue {
     /// represent an upper bound, returns true if the two are incompatible.
     ///
     /// ```
-    /// use aries_model::bounds::BoundValue;
+    /// use aries_core::*;
     /// assert!(!BoundValue::lb(5).compatible_with_symmetric(BoundValue::ub(4)));
     /// assert!(BoundValue::lb(5).compatible_with_symmetric(BoundValue::ub(5)));
     /// assert!(BoundValue::lb(5).compatible_with_symmetric(BoundValue::ub(6)));
@@ -44,7 +44,7 @@ impl BoundValue {
     /// This should be called with a lower and an upper bound.
     ///
     /// ```
-    /// use aries_model::bounds::BoundValue;
+    /// use aries_core::BoundValue;
     /// assert!(!BoundValue::lb(5).equal_to_symmetric(BoundValue::ub(4)));
     /// assert!(BoundValue::lb(5).equal_to_symmetric(BoundValue::ub(5)));
     /// assert!(!BoundValue::lb(5).equal_to_symmetric(BoundValue::ub(6)));
@@ -126,7 +126,7 @@ impl std::ops::Sub<BoundValueAdd> for BoundValue {
 /// with the one of the bound value.
 ///
 /// ```
-/// use aries_model::bounds::{BoundValue, BoundValueAdd};
+/// use aries_core::*;
 /// let ub_add = BoundValueAdd::on_ub(5);
 /// let lb_add = BoundValueAdd::on_ub(-4);
 /// assert_eq!(BoundValue::ub(3) + BoundValueAdd::on_ub(5), BoundValue::ub(8));
@@ -140,12 +140,18 @@ impl std::ops::Sub<BoundValueAdd> for BoundValue {
 pub struct BoundValueAdd(IntCst);
 
 impl BoundValueAdd {
-    /// The zero value addition, independently of whether it represents applies on lower or upper bounds.
+    /// The zero value addition, independently of whether it represents applies on lower or upper literals.
     pub const ZERO: BoundValueAdd = BoundValueAdd(0);
 
+    /// Adding this to a [BoundValue] is equivalent to relaxing it to the next value.
+    pub const RELAXATION: BoundValueAdd = BoundValueAdd(1);
+
+    /// Adding this to a [BoundValue] is equivalent to tightening it to the next value.
+    pub const TIGHTENING: BoundValueAdd = BoundValueAdd(-1);
+
     /// The maximum representable value. Not that using it anything else than a default value for comparison
     /// is likely to result in an overflow.
-    pub const MAX: BoundValueAdd = BoundValueAdd(IntCst::MAX);
+    pub const MAX: BoundValueAdd = BoundValueAdd(INT_CST_MAX);
 
     /// Construct the BVA that represents an update on a lower bound
     pub fn on_lb(increment: IntCst) -> Self {
@@ -192,7 +198,7 @@ impl std::ops::Add<BoundValueAdd> for BoundValueAdd {
 
 #[cfg(test)]
 mod test {
-    use crate::bounds::{BoundValue, BoundValueAdd};
+    use crate::*;
 
     #[test]
     fn test_compatibility() {
diff --git a/core/src/lib.rs b/core/src/lib.rs
new file mode 100644
index 00000000..66fcc66f
--- /dev/null
+++ b/core/src/lib.rs
@@ -0,0 +1,26 @@
+use crate::state::Cause;
+pub use bound_value::*;
+pub use lit::*;
+pub use var_bound::*;
+pub use variable::*;
+
+mod bound_value;
+mod lit;
+pub mod literals;
+pub mod state;
+mod var_bound;
+mod variable;
+
+/// Identifies an external writer to the model.
+#[derive(Ord, PartialOrd, PartialEq, Eq, Copy, Clone, Hash, Debug)]
+pub struct WriterId(pub u8);
+
+impl WriterId {
+    pub fn new(num: impl Into<u8>) -> WriterId {
+        WriterId(num.into())
+    }
+
+    pub fn cause(&self, cause: impl Into<u32>) -> Cause {
+        Cause::inference(*self, cause)
+    }
+}
diff --git a/model/src/bounds/bound.rs b/core/src/lit.rs
similarity index 64%
rename from model/src/bounds/bound.rs
rename to core/src/lit.rs
index 4c64adf6..d0ce6414 100644
--- a/model/src/bounds/bound.rs
+++ b/core/src/lit.rs
@@ -1,11 +1,6 @@
-use crate::bounds::var_bound::VarBound;
-use crate::bounds::BoundValue;
-use crate::lang::{BVar, ConversionError};
-use crate::lang::{IntCst, VarRef};
+use crate::*;
 use core::convert::{From, Into};
 use std::cmp::Ordering;
-use std::convert::TryFrom;
-use std::mem::transmute;
 
 /// A literal `Lit` represents a lower or upper bound on a discrete variable
 /// (i.e. an integer, boolean or symbolic variable).
@@ -19,19 +14,16 @@ use std::mem::transmute;
 /// The `unpack()` method extract all fields into a tuple.
 ///
 /// ```
-/// use aries_model::Model;
-/// use aries_model::lang::VarRef;
-/// use aries_model::bounds::{Lit, Relation};
-/// let mut model = Model::<&'static str>::new();
-/// let x = model.new_bvar("X");
-/// let x_is_true: Lit = x.true_lit();
-/// let x_is_false: Lit = x.false_lit();
-/// let y = model.new_ivar(0, 10, "Y");
+/// use aries_core::*;
+/// use aries_core::state::IntDomains;
+/// let mut state = IntDomains::new();
+/// let x = state.new_var(0, 1);
+/// let x_is_true: Lit = x.geq(1);
+/// let x_is_false: Lit = !x_is_true;
+/// let y = state.new_var(0, 10);
 /// let y_geq_5 = Lit::geq(y, 5);
 ///
 /// // the `<=` is internally converted into a `<`
-/// // the variable is converted into a `VarRef`
-/// let y: VarRef = y.into();
 /// assert_eq!(y_geq_5.variable(), y);
 /// assert_eq!(y_geq_5.relation(), Relation::Gt);
 /// assert_eq!(y_geq_5.value(), 4);
@@ -47,47 +39,36 @@ use std::mem::transmute;
 ///  - value
 ///
 /// As a result, ordering a vector of `Lit`s will group them by variable, then among literals on the same variable by relation.
-/// An important invariant is that, in a sorted list, a bound can only entail the bounds immediately following it.
+/// An important invariant is that, in a sorted list, a bound can only entail the literals immediately following it.
 ///
 /// ```
-/// use aries_model::Model;
-/// use aries_model::bounds::Lit;
-/// let mut model = Model::<&'static str>::new();
-/// let x = model.new_ivar(0, 10, "X");
-/// let y = model.new_ivar(0, 10, "Y");
-/// let mut bounds = vec![Lit::geq(y, 4), Lit::geq(x,1), Lit::leq(x, 3), Lit::leq(x, 4), Lit::leq(x, 6), Lit::geq(x,2)];
-/// bounds.sort();
-/// assert_eq!(bounds, vec![Lit::geq(x,2), Lit::geq(x,1), Lit::leq(x, 3), Lit::leq(x, 4), Lit::leq(x, 6), Lit::geq(y, 4)]);
+/// use aries_core::*;
+/// let x = VarRef::from_u32(1);
+/// let y = VarRef::from_u32(2);
+/// let mut literals = vec![Lit::geq(y, 4), Lit::geq(x,1), Lit::leq(x, 3), Lit::leq(x, 4), Lit::leq(x, 6), Lit::geq(x,2)];
+/// literals.sort();
+/// assert_eq!(literals, vec![Lit::geq(x,2), Lit::geq(x,1), Lit::leq(x, 3), Lit::leq(x, 4), Lit::leq(x, 6), Lit::geq(y, 4)]);
 /// ```
 #[derive(Copy, Clone, Eq, PartialEq, Hash, Ord, PartialOrd)]
 pub struct Lit {
-    /// Union of the variable (highest 31 bits) and the relation (lowest bit)
+    /// Either the lower bound or the upper bound ot hte affected variable.
+    ///
+    /// Implemented as the union of the variable (highest 31 bits) and the relation (lowest bit)
     /// This encoding allows:
-    ///  - to very efficiently check whether two bounds have the same `(variable, relation)` part
+    ///  - to very efficiently check whether two literals have the same `(variable, relation)` part
     ///    which is one of the critical operation in `entails`.
     ///  - to use as an index in a table: each variable will have two slots: one of the LEQ relation
     ///    and one for the GT relation
-    //  TODO: use a VarBound
-    pub(in crate::bounds) var_rel: u32,
+    var_bound: VarBound,
     /// +/- the value of the relation. The value of a GT relation is negated before being stored.
     /// This design allows to test entailment without testing the relation of the Bound
-    pub(in crate::bounds) raw_value: BoundValue,
+    raw_value: BoundValue,
 }
 
 #[derive(Ord, PartialOrd, Eq, PartialEq, Debug, Copy, Clone)]
-#[repr(u8)]
 pub enum Relation {
-    Gt = 0,
-    Leq = 1,
-}
-
-impl std::ops::Not for Relation {
-    type Output = Self;
-
-    #[inline]
-    fn not(self) -> Self::Output {
-        unsafe { transmute((self as u8) ^ 0x1) }
-    }
+    Gt,
+    Leq,
 }
 
 impl std::fmt::Display for Relation {
@@ -99,9 +80,6 @@ impl std::fmt::Display for Relation {
     }
 }
 
-const REL_MASK: u32 = 0x1;
-const VAR_MASK: u32 = !REL_MASK;
-
 impl Lit {
     /// A literal that is always true. It is defined by stating that the special variable [VarRef::ZERO] is
     /// lesser than or equal to 0, which is always true.
@@ -112,41 +90,41 @@ impl Lit {
     #[inline]
     pub const fn from_parts(var_bound: VarBound, value: BoundValue) -> Self {
         Lit {
-            var_rel: var_bound.to_u32(),
+            var_bound,
             raw_value: value,
         }
     }
 
     #[inline]
     pub const fn new(variable: VarRef, relation: Relation, value: IntCst) -> Self {
-        let var_part = variable.to_u32() << 1;
-        let relation_part = relation as u32;
-        let raw_value = match relation {
-            Relation::Leq => BoundValue::ub(value),
-            Relation::Gt => BoundValue::lb(value + 1),
-        };
-        Lit {
-            var_rel: var_part | relation_part,
-            raw_value,
+        match relation {
+            Relation::Leq => Lit {
+                var_bound: VarBound::ub(variable),
+                raw_value: BoundValue::ub(value),
+            },
+            Relation::Gt => Lit {
+                var_bound: VarBound::lb(variable),
+                raw_value: BoundValue::lb(value + 1),
+            },
         }
     }
 
     #[inline]
     pub fn variable(self) -> VarRef {
-        let var_part = self.var_rel & VAR_MASK;
-        let var = var_part >> 1;
-        VarRef::from(var)
+        self.var_bound.variable()
     }
 
     #[inline]
-    pub fn relation(self) -> Relation {
-        let rel_part = self.var_rel & REL_MASK;
-        let rel = rel_part as u8;
-        unsafe { transmute(rel) }
+    pub const fn relation(self) -> Relation {
+        if self.var_bound.is_ub() {
+            Relation::Leq
+        } else {
+            Relation::Gt
+        }
     }
 
     #[inline]
-    pub fn value(self) -> IntCst {
+    pub const fn value(self) -> IntCst {
         match self.relation() {
             Relation::Leq => self.raw_value.as_ub(),
             Relation::Gt => self.raw_value.as_lb() - 1,
@@ -154,12 +132,12 @@ impl Lit {
     }
 
     #[inline]
-    pub fn affected_bound(self) -> VarBound {
-        VarBound::from_raw(self.var_rel)
+    pub const fn affected_bound(self) -> VarBound {
+        self.var_bound
     }
 
     #[inline]
-    pub fn bound_value(self) -> BoundValue {
+    pub const fn bound_value(self) -> BoundValue {
         self.raw_value
     }
 
@@ -182,31 +160,24 @@ impl Lit {
         Lit::new(var.into(), Relation::Gt, val)
     }
 
-    pub fn is_true(v: BVar) -> Lit {
-        Lit::geq(v, 1)
-    }
-    pub fn is_false(v: BVar) -> Lit {
-        Lit::leq(v, 0)
-    }
-
     #[inline]
     pub const fn not(self) -> Self {
         Lit {
-            var_rel: self.var_rel ^ 0x1,
+            var_bound: self.var_bound.symmetric_bound(),
             raw_value: self.raw_value.neg(),
         }
     }
 
     #[inline]
     pub fn entails(self, other: Lit) -> bool {
-        self.var_rel == other.var_rel && self.raw_value.stronger(other.raw_value)
+        self.var_bound == other.var_bound && self.raw_value.stronger(other.raw_value)
     }
 
     pub fn unpack(self) -> (VarRef, Relation, IntCst) {
         (self.variable(), self.relation(), self.value())
     }
 
-    /// An ordering that will group bounds by (given from highest to lowest priority):
+    /// An ordering that will group literals by (given from highest to lowest priority):
     ///  - variable
     ///  - affected bound (lower, upper)
     ///  - by value of the bound
@@ -224,12 +195,6 @@ impl std::ops::Not for Lit {
     }
 }
 
-impl From<BVar> for Lit {
-    fn from(v: BVar) -> Self {
-        v.true_lit()
-    }
-}
-
 impl From<bool> for Lit {
     fn from(b: bool) -> Self {
         if b {
@@ -240,18 +205,6 @@ impl From<bool> for Lit {
     }
 }
 
-impl TryFrom<Lit> for bool {
-    type Error = ConversionError;
-
-    fn try_from(value: Lit) -> Result<Self, Self::Error> {
-        match value {
-            Lit::TRUE => Ok(true),
-            Lit::FALSE => Ok(false),
-            _ => Err(ConversionError::NotConstant),
-        }
-    }
-}
-
 impl std::fmt::Debug for Lit {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         match *self {
diff --git a/core/src/literals.rs b/core/src/literals.rs
new file mode 100644
index 00000000..a90727c3
--- /dev/null
+++ b/core/src/literals.rs
@@ -0,0 +1,9 @@
+pub use disjunction::*;
+pub use implication_graph::*;
+pub use lit_set::*;
+pub use watches::*;
+
+mod disjunction;
+mod implication_graph;
+mod lit_set;
+mod watches;
diff --git a/model/src/bounds/disjunction.rs b/core/src/literals/disjunction.rs
similarity index 98%
rename from model/src/bounds/disjunction.rs
rename to core/src/literals/disjunction.rs
index 130a3545..0b72c994 100644
--- a/model/src/bounds/disjunction.rs
+++ b/core/src/literals/disjunction.rs
@@ -1,4 +1,4 @@
-use crate::bounds::{Lit, Relation};
+use crate::*;
 use std::borrow::Borrow;
 
 /// A set of literals representing a disjunction.
@@ -126,7 +126,6 @@ impl AsRef<[Lit]> for Disjunction {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::lang::{IntCst, VarRef};
 
     fn leq(var: VarRef, val: IntCst) -> Lit {
         Lit::leq(var, val)
diff --git a/core/src/literals/implication_graph.rs b/core/src/literals/implication_graph.rs
new file mode 100644
index 00000000..523ad204
--- /dev/null
+++ b/core/src/literals/implication_graph.rs
@@ -0,0 +1,146 @@
+use crate::literals::{LitSet, Watches};
+use crate::*;
+
+/// An implication in the form of a 2-SAT network.
+///
+/// It allows answering in polynomial time whether a given literal implies another one (directly or indirectly)
+///
+/// The network accounts for implication between literals on the same variable:
+/// `(X < 0) => (X < 1)` is implicit and always true.
+/// Thus, with an explicit implication `(Y < 0) => (X < 0)`, the network is able to infer the following facts:
+///  - `(Y < 0)  => (X < 0)`
+///  - `(Y < 0)  => (X < 1)`
+///  - `(Y < -1) => (X < 0)`
+///  - `(Y < 0)  => (Y < 1)`   (directly accessible with [`Lit::entails`])
+///
+/// # Limitations
+///
+/// The network does not check contradiction in the form of `x => !x`.
+/// It also does not check for edge duplication, which might a cause for inefficiencies.
+///
+/// # Example
+/// ```
+/// use aries_core::*;
+/// use aries_core::literals::ImplicationGraph;
+/// let mut set = ImplicationGraph::empty();
+/// let v1 = VarRef::from_u32(3); // arbitrary variable
+/// let v2 = VarRef::from_u32(4); // arbitrary variable
+/// assert!(!set.implies(v1.leq(0), v2.leq(0)));
+/// set.add_implication(v1.leq(0), v2.leq(0));
+/// assert!(set.implies(v1.leq(0), v2.leq(0)));
+/// assert!(set.implies(v1.leq(0), v2.leq(1)));
+/// assert!(set.implies(v1.leq(-1), v2.leq(0)));
+/// assert!(!set.implies(v1.leq(1), v2.leq(0)));
+/// ```
+#[derive(Clone, Default)]
+pub struct ImplicationGraph {
+    edges: Watches<Lit>,
+    num_edges: usize,
+}
+
+impl ImplicationGraph {
+    /// Creates an empty implication graph
+    pub fn empty() -> Self {
+        Self::default()
+    }
+
+    /// Record the fact that `from` implies `to`.
+    pub fn add_implication(&mut self, from: Lit, to: Lit) {
+        if to == Lit::TRUE || from == Lit::FALSE || from.entails(to) {
+            return;
+        }
+        self.num_edges += 1;
+        self.edges.add_watch(to, from);
+        self.edges.add_watch(!from, !to);
+        debug_assert!(self.implies(from, to));
+        debug_assert!(self.implies(!to, !from));
+    }
+
+    /// Return true if there is a direct or indirect implication `x => y`.
+    pub fn implies(&self, x: Lit, y: Lit) -> bool {
+        if y == Lit::TRUE || x == Lit::FALSE || x.entails(y) {
+            return true;
+        }
+        if self.edges.watches_on(!y).next().is_none() {
+            // fail fast: no incoming edges to y, which thus is not reachable
+            // this is possible to check because, for each  (x -> y) edge we have a (!y -> !x) edge.
+            return false;
+        }
+        // starting from `x`, do a depth first search in the implication graph,
+        // looking for a literal that entails `y`
+
+        // list of literals that were previously encountered in search
+        let mut visited = LitSet::with_capacity(64);
+        let mut queue = Vec::with_capacity(64);
+        queue.push(x);
+
+        // dfs through implications
+        while let Some(curr) = queue.pop() {
+            if visited.contains(curr) {
+                continue;
+            }
+            visited.insert(curr);
+            for next in self.edges.watches_on(curr) {
+                if next.entails(y) {
+                    return true;
+                } else {
+                    queue.push(next);
+                }
+            }
+        }
+        false
+    }
+
+    pub fn direct_implications_of(&self, lit: Lit) -> impl Iterator<Item = Lit> + '_ {
+        self.edges.watches_on(lit)
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use crate::literals::ImplicationGraph;
+    use crate::*;
+
+    const A: VarRef = VarRef::from_u32(0);
+    const B: VarRef = VarRef::from_u32(1);
+    const C: VarRef = VarRef::from_u32(2);
+    const D: VarRef = VarRef::from_u32(3);
+
+    #[test]
+    fn test_implications() {
+        let mut g = ImplicationGraph::empty();
+
+        assert!(g.implies(A.leq(0), A.leq(0)));
+        assert!(g.implies(A.leq(0), A.leq(1)));
+        assert!(!g.implies(A.leq(0), B.leq(0)));
+        assert!(!g.implies(A.leq(0), A.leq(-1)));
+
+        g.add_implication(A.leq(1), B.leq(1));
+        assert!(g.implies(A.leq(1), B.leq(1)));
+        assert!(g.implies(A.leq(0), B.leq(1)));
+        assert!(g.implies(A.leq(1), B.leq(2)));
+        assert!(g.implies(A.leq(0), B.leq(2)));
+        assert!(!g.implies(A.leq(1), B.leq(0)));
+        assert!(!g.implies(A.leq(1), B.leq(0)));
+
+        g.add_implication(B.leq(2), C.leq(2));
+        assert!(g.implies(A.leq(1), C.leq(2)));
+        assert!(g.implies(A.leq(1), C.leq(3)));
+        assert!(!g.implies(A.leq(1), C.leq(1)));
+        assert!(g.implies(A.leq(0), C.leq(2)));
+        assert!(!g.implies(A.leq(2), C.leq(2)));
+    }
+
+    #[test]
+    fn test_implication_cycle() {
+        let mut g = ImplicationGraph::empty();
+
+        g.add_implication(A.leq(0), B.leq(0));
+        g.add_implication(B.leq(0), A.leq(0));
+
+        g.add_implication(C.leq(0), D.leq(0));
+        g.add_implication(D.leq(0), C.leq(0));
+
+        assert!(!g.implies(A.leq(0), C.leq(0)))
+    }
+}
diff --git a/core/src/literals/lit_set.rs b/core/src/literals/lit_set.rs
new file mode 100644
index 00000000..364ee0a1
--- /dev/null
+++ b/core/src/literals/lit_set.rs
@@ -0,0 +1,240 @@
+use crate::*;
+use std::array::TryFromSliceError;
+use std::convert::{TryFrom, TryInto};
+
+/// Set of literals.
+///
+/// The set is said to contain a literal `l` iff
+///  - `l` was previously inserted in the set; or
+///  - another literal `l2` was previously inserted such that `l2.entails(l)`.
+///
+/// # Limitations
+///
+/// The implementation is optimized for small, dynamic sets as it stores literals in an unsorted vector.
+/// A specialized implementation for large sets could be implemented with a `Map<VarBound, BoundValue>`
+///
+/// # Example
+/// ```
+/// use aries_core::literals::LitSet;
+/// use aries_core::VarRef;
+/// let mut set = LitSet::empty();
+/// let var = VarRef::from_u32(3); // arbitrary variable
+/// assert!(!set.contains(var.leq(0)));
+/// set.insert(var.leq(0));
+/// assert!(set.contains(var.leq(0)));
+/// assert!(set.contains(var.leq(1))); // present because entailed by `var.leq(0)`
+/// assert!(!set.contains(var.leq(-1))); // not present as it is not entailed
+/// ```
+#[derive(Default, Clone, Debug)]
+pub struct LitSet {
+    elements: Vec<Lit>,
+}
+
+impl LitSet {
+    pub const EMPTY: LitSet = Self { elements: vec![] };
+
+    pub fn empty() -> Self {
+        Self::default()
+    }
+
+    pub fn with_capacity(capacity: usize) -> Self {
+        LitSet {
+            elements: Vec::with_capacity(capacity),
+        }
+    }
+
+    pub fn into_sorted(self) -> StableLitSet {
+        StableLitSet::new(self)
+    }
+
+    pub fn literals(&self) -> impl Iterator<Item = Lit> + '_ {
+        self.elements.iter().copied()
+    }
+
+    pub fn contains(&self, elem: Lit) -> bool {
+        self.elements.iter().any(|l| l.entails(elem))
+    }
+
+    /// Insert a literal `lit` into the set.
+    ///
+    /// Note that all literals directly implied by `lit` are also implicitly inserted.
+    pub fn insert(&mut self, lit: Lit) {
+        if lit == Lit::TRUE {
+            return;
+        }
+        for l in self.elements.iter_mut() {
+            if l.entails(lit) {
+                return;
+            } else if lit.entails(*l) {
+                *l = lit;
+                return;
+            }
+        }
+        self.elements.push(lit)
+    }
+
+    /// Remove a literal `rm` from the set.
+    ///
+    /// Removal needs to account for implicit literals.
+    ///
+    /// `{ (X <= 1) } is equivalent to { (X <= 1), (X <= 2), (X <= 3), ...}`
+    /// Hence, if we remove (X <= 2) from this set, we are left with
+    /// `{ (X <= 3), (X <= 4), ... }`, which is normalized to `{ (X <= 3) }`
+    ///
+    /// When doing this operation, it is likely that the resulting literal would be tautological.
+    /// For instance, consider is a boolean variable `B` whose value is either `0` or `1`.
+    /// Removing `(B <= 0)` from the set `{ (B <= 0) }` would result in the set `{ (B <= 1) }`.
+    /// The literal `(B <= 1)` is tautological and can be ignored in the set.
+    /// The method will use the method `tautological` to determine which literals are always true.
+    ///
+    pub fn remove(&mut self, rm: Lit, tautology: impl Fn(Lit) -> bool) {
+        for (i, l) in self.elements.iter_mut().enumerate() {
+            if l.entails(rm) {
+                let new_lit = Lit::from_parts(rm.affected_bound(), rm.bound_value() + BoundValueAdd::RELAXATION);
+                if tautology(new_lit) {
+                    // the literal that we would insert is a tautology, simply remove the literal
+                    self.elements.swap_remove(i);
+                } else {
+                    // non tautological, insert the literal.
+                    *l = new_lit;
+                }
+                // by construction, there is at most one literal on a VarBound
+                return;
+            }
+        }
+    }
+}
+
+impl<T: IntoIterator<Item = Lit>> From<T> for LitSet {
+    fn from(lits: T) -> Self {
+        let mut set = LitSet::empty();
+        for l in lits {
+            set.insert(l);
+        }
+        set
+    }
+}
+
+/// A set of literal in a canonical form and which can thus be used for comparison.
+#[derive(Clone, Ord, PartialOrd, Eq, PartialEq, Hash, Debug)]
+pub struct StableLitSet {
+    elements: Vec<Lit>,
+}
+
+impl StableLitSet {
+    pub const EMPTY: Self = Self { elements: vec![] };
+
+    pub fn new(mut set: LitSet) -> Self {
+        set.elements.sort();
+        Self { elements: set.elements }
+    }
+
+    pub fn literals(&self) -> impl Iterator<Item = Lit> + '_ {
+        self.elements.iter().copied()
+    }
+
+    pub fn len(&self) -> usize {
+        self.elements.len()
+    }
+
+    pub fn is_empty(&self) -> bool {
+        self.elements.is_empty()
+    }
+}
+
+impl<T: IntoIterator<Item = Lit>> From<T> for StableLitSet {
+    fn from(lits: T) -> Self {
+        Self::new(lits.into())
+    }
+}
+
+impl<const N: usize> TryFrom<&StableLitSet> for [Lit; N] {
+    type Error = TryFromSliceError;
+
+    fn try_from(value: &StableLitSet) -> Result<Self, Self::Error> {
+        value.elements.as_slice().try_into()
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    const A: VarRef = VarRef::from_u32(1);
+    const B: VarRef = VarRef::from_u32(2);
+    const C: VarRef = VarRef::from_u32(3);
+
+    #[test]
+    fn test_lit_set() {
+        let mut set = LitSet::empty();
+
+        assert!(!set.contains(A.leq(1)));
+        assert!(!set.contains(A.geq(1)));
+        assert_eq!(set.elements.len(), 0);
+
+        set.insert(A.leq(1));
+        assert_eq!(set.elements.len(), 1);
+        assert!(set.contains(A.leq(1)));
+        assert!(set.contains(A.leq(2)));
+        assert!(!set.contains(A.leq(0)));
+
+        // insert (A <= 2) that should (implicitly) already be in the set
+        set.insert(A.leq(2));
+        assert_eq!(set.elements.len(), 1);
+        assert!(set.contains(A.leq(1)));
+        assert!(set.contains(A.leq(2)));
+        assert!(!set.contains(A.leq(0)));
+        assert!(!set.contains(B.leq(10)));
+        assert!(!set.contains(C.leq(10)));
+
+        // insert (A <= 0) that should subsume the (A <= 1) literal present in the set
+        set.insert(A.leq(0));
+        assert_eq!(set.elements.len(), 1);
+        assert!(set.contains(A.leq(1)));
+        assert!(set.contains(A.leq(2)));
+        assert!(set.contains(A.leq(0)));
+        assert!(!set.contains(A.leq(-1)));
+        assert!(!set.contains(B.leq(10)));
+        assert!(!set.contains(C.leq(10)));
+
+        set.insert(B.geq(5));
+        assert_eq!(set.elements.len(), 2);
+        set.insert(C.geq(5));
+        assert_eq!(set.elements.len(), 3);
+        set.insert(C.geq(3));
+        assert!(set.contains(A.leq(0)));
+        assert!(!set.contains(A.leq(-1)));
+        assert!(set.contains(B.geq(5)));
+        assert!(!set.contains(B.geq(6)));
+        assert!(set.contains(C.geq(5)));
+        assert!(!set.contains(C.geq(6)));
+        assert!(!set.contains(B.leq(10)));
+        assert!(!set.contains(C.leq(10)));
+    }
+
+    #[test]
+    fn test_lit_set_removal() {
+        let mut set = LitSet::empty();
+
+        let tauto = |l| A.leq(4).entails(l);
+
+        assert!(!set.contains(A.leq(1)));
+        assert!(!set.contains(A.geq(1)));
+        assert_eq!(set.elements.len(), 0);
+
+        set.insert(A.leq(1));
+        assert_eq!(set.elements.len(), 1);
+        assert!(set.contains(A.leq(1)));
+        assert!(set.contains(A.leq(2)));
+        assert!(!set.contains(A.leq(0)));
+
+        set.remove(A.leq(1), tauto);
+        assert_eq!(set.elements.len(), 1);
+        assert!(!set.contains(A.leq(1)));
+        assert!(set.contains(A.leq(2)));
+        assert!(!set.contains(A.leq(0)));
+
+        set.remove(A.leq(3), tauto);
+        assert_eq!(set.elements.len(), 0);
+    }
+}
diff --git a/model/src/bounds/watches.rs b/core/src/literals/watches.rs
similarity index 92%
rename from model/src/bounds/watches.rs
rename to core/src/literals/watches.rs
index 2ff40465..dd971178 100644
--- a/model/src/bounds/watches.rs
+++ b/core/src/literals/watches.rs
@@ -1,9 +1,7 @@
+use crate::*;
 use aries_collections::ref_store::RefVec;
 
-use crate::bounds::var_bound::VarBound;
-use crate::bounds::{BoundValue, Lit};
-
-/// A set of bounds watches on bound changes.
+/// A set of literals watches on bound changes.
 /// The event watches are all on the same bound (i.e. the lower or the upper bound) of a single variable.
 #[derive(Clone)]
 pub struct WatchSet<Watcher> {
@@ -17,7 +15,7 @@ impl<Watcher> WatchSet<Watcher> {
     pub fn add_watch(&mut self, watcher: Watcher, literal: Lit) {
         self.watches.push(Watch {
             watcher,
-            guard: literal.raw_value,
+            guard: literal.bound_value(),
         })
     }
 
@@ -42,7 +40,7 @@ impl<Watcher> WatchSet<Watcher> {
     {
         self.watches
             .iter()
-            .any(|w| w.watcher == watcher && literal.raw_value.stronger(w.guard))
+            .any(|w| w.watcher == watcher && literal.bound_value().stronger(w.guard))
     }
 
     pub fn watches_on(&self, literal: Lit) -> impl Iterator<Item = Watcher> + '_
@@ -50,7 +48,7 @@ impl<Watcher> WatchSet<Watcher> {
         Watcher: Copy,
     {
         self.watches.iter().filter_map(move |w| {
-            if literal.raw_value.stronger(w.guard) {
+            if literal.bound_value().stronger(w.guard) {
                 Some(w.watcher)
             } else {
                 None
@@ -65,7 +63,7 @@ impl<Watcher> WatchSet<Watcher> {
     pub fn move_watches_to(&mut self, literal: Lit, out: &mut WatchSet<Watcher>) {
         let mut i = 0;
         while i < self.watches.len() {
-            if literal.raw_value.stronger(self.watches[i].guard) {
+            if literal.bound_value().stronger(self.watches[i].guard) {
                 let w = self.watches.swap_remove(i);
                 out.watches.push(w);
             } else {
@@ -88,10 +86,7 @@ pub struct Watch<Watcher> {
 }
 impl<Watcher> Watch<Watcher> {
     pub fn to_lit(&self, var_bound: VarBound) -> Lit {
-        Lit {
-            var_rel: u32::from(var_bound),
-            raw_value: self.guard,
-        }
+        Lit::from_parts(var_bound, self.guard)
     }
 }
 
@@ -169,14 +164,12 @@ impl<Watcher> Default for Watches<Watcher> {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::bounds::Lit;
-    use crate::Model;
+    use crate::Lit;
 
     #[test]
     fn test_watches() {
-        let mut model: Model<&'static str> = Model::new();
-        let a = model.new_ivar(0, 10, "a");
-        let b = model.new_ivar(0, 10, "b");
+        let a = VarRef::from_u32(1);
+        let b = VarRef::from_u32(2);
 
         let watches = &mut Watches::new();
 
diff --git a/model/src/state.rs b/core/src/state.rs
similarity index 89%
rename from model/src/state.rs
rename to core/src/state.rs
index f8ea3d7a..bff46524 100644
--- a/model/src/state.rs
+++ b/core/src/state.rs
@@ -4,7 +4,6 @@ mod domains;
 mod event;
 mod explanation;
 mod int_domains;
-mod presence_graph;
 
 pub use cause::*;
 pub use domain::*;
@@ -13,7 +12,7 @@ pub use event::*;
 pub use explanation::*;
 pub use int_domains::*;
 
-use crate::bounds::Lit;
+use crate::Lit;
 
 /// Represents a triggered event of setting a conflicting literal.
 #[derive(PartialEq, Eq, Debug, Copy, Clone)]
diff --git a/model/src/state/cause.rs b/core/src/state/cause.rs
similarity index 99%
rename from model/src/state/cause.rs
rename to core/src/state/cause.rs
index 271a81d8..ac517209 100644
--- a/model/src/state/cause.rs
+++ b/core/src/state/cause.rs
@@ -1,4 +1,4 @@
-use crate::bounds::Lit;
+use crate::Lit;
 use crate::WriterId;
 
 /// Cause of an event that originates from outside of the solver.
diff --git a/model/src/state/domain.rs b/core/src/state/domain.rs
similarity index 95%
rename from model/src/state/domain.rs
rename to core/src/state/domain.rs
index 55eb7eaf..a7118217 100644
--- a/model/src/state/domain.rs
+++ b/core/src/state/domain.rs
@@ -1,4 +1,4 @@
-use crate::lang::IntCst;
+use crate::IntCst;
 use std::fmt::Formatter;
 
 #[derive(Clone, Ord, PartialOrd, Eq, PartialEq, Debug)]
@@ -39,9 +39,9 @@ impl std::fmt::Display for IntDomain {
 /// Represents the domain of an optional variable
 #[derive(Eq, PartialEq, Copy, Clone)]
 pub enum OptDomain {
-    /// The variable is necessarily present and must take a value in the given bounds.
+    /// The variable is necessarily present and must take a value in the given literals.
     Present(IntCst, IntCst),
-    /// It is unknown whether the variable is present but if it is it must take a value in the given bounds.
+    /// It is unknown whether the variable is present but if it is it must take a value in the given literals.
     Unknown(IntCst, IntCst),
     /// The variable is known to be absent
     Absent,
diff --git a/model/src/state/domains.rs b/core/src/state/domains.rs
similarity index 89%
rename from model/src/state/domains.rs
rename to core/src/state/domains.rs
index 6ec5c1e4..ca9b6e5f 100644
--- a/model/src/state/domains.rs
+++ b/core/src/state/domains.rs
@@ -1,10 +1,9 @@
-use crate::bounds::{BoundValue, Disjunction, Lit, VarBound};
-use crate::lang::{IntCst, VarRef};
+use crate::literals::{Disjunction, ImplicationGraph};
 use crate::state::cause::{DirectOrigin, Origin};
 use crate::state::event::Event;
 use crate::state::int_domains::IntDomains;
-use crate::state::presence_graph::TwoSatTree;
 use crate::state::{Cause, Explainer, Explanation, InvalidUpdate, OptDomain};
+use crate::*;
 use aries_backtrack::{Backtrack, DecLvl, DecisionLevelClass, EventIndex, ObsTrail};
 use aries_collections::ref_store::RefMap;
 use std::collections::BinaryHeap;
@@ -32,7 +31,7 @@ pub struct Domains {
     /// is true if and only if the variable is present.
     presence: RefMap<VarRef, Lit>,
     /// A graph to encode the relations between presence variables.
-    presence_graph: TwoSatTree,
+    implications: ImplicationGraph,
     /// A queue used internally when building explanations. Only useful to avoid repeated allocations.
     queue: BinaryHeap<InQueueLit>,
 }
@@ -42,7 +41,7 @@ impl Domains {
         let domains = Domains {
             doms: IntDomains::new(),
             presence: Default::default(),
-            presence_graph: Default::default(),
+            implications: Default::default(),
             queue: Default::default(),
         };
         debug_assert!(domains.entails(Lit::TRUE));
@@ -54,13 +53,34 @@ impl Domains {
         self.doms.new_var(lb, ub)
     }
 
-    pub fn new_presence_literal(&mut self, scope: Lit) -> Lit {
-        let lit = self.new_var(0, 1).geq(1);
-        self.presence_graph.add_implication(lit, scope);
-        if self.entails(!scope) {
-            let prop_result = self.set_impl(!lit, DirectOrigin::ImplicationPropagation(!scope));
-            assert_eq!(prop_result, Ok(true));
+    /// Records a direct implication `from => to`
+    ///
+    /// # Assumptions
+    ///
+    /// - `from` and `to` are non-optional
+    /// - Propagating the implication will not create an inconsistencies
+    ///
+    /// The current implementation will panic in those cases, but this check might be done in
+    /// debug-only builds at a later time.
+    #[rustfmt::skip]
+    pub fn add_implication(&mut self, from: Lit, to: Lit) {
+        assert_eq!(self.presence(from.variable()), Lit::TRUE, "Implication only supported between non-optional variables");
+        assert_eq!(self.presence(to.variable()), Lit::TRUE, "Implication only supported between non-optional variables");
+        self.implications.add_implication(from, to);
+        if self.entails(from) {
+            let prop_result = self.set_impl(to, DirectOrigin::ImplicationPropagation(from));
+            assert!(matches!(prop_result, Ok(_)), "Inconsistency on the addition of implies({:?}, {:?}", from, to);
+        }
+        if self.entails(!to) {
+            let prop_result = self.set_impl(!from, DirectOrigin::ImplicationPropagation(!to));
+            assert!(matches!(prop_result, Ok(_)), "Inconsistency on the addition of implies({:?}, {:?}", from, to);
         }
+    }
+
+    #[cfg(test)]
+    fn new_presence_literal(&mut self, scope: Lit) -> Lit {
+        let lit = self.new_var(0, 1).geq(1);
+        self.add_implication(lit, scope);
         lit
     }
 
@@ -82,8 +102,17 @@ impl Domains {
     pub fn only_present_with(&self, a: VarRef, b: VarRef) -> bool {
         let prez_a = self.presence(a);
         let prez_b = self.presence(b);
-        // prez_a => prez_b
-        prez_b == Lit::TRUE || prez_a.entails(prez_b) || self.presence_graph.implies(prez_a, prez_b)
+        self.implications.implies(prez_a, prez_b)
+    }
+
+    /// Returns true if `a` is known to imply `b`
+    pub fn implies(&self, a: Lit, b: Lit) -> bool {
+        self.implications.implies(a, b)
+    }
+
+    /// Returns true if `a` and `b` are known to be exclusive
+    pub fn exclusive(&self, a: Lit, b: Lit) -> bool {
+        self.implies(a, !b) || self.implies(b, !a)
     }
 
     /// Returns true if we know that two variable are always present jointly.
@@ -225,12 +254,11 @@ impl Domains {
         // variable must be optional
         debug_assert_ne!(prez, Lit::TRUE);
         // invariant: optional variable cannot be involved in implications
-        debug_assert_eq!(
-            self.presence_graph
-                .direct_implications_of(Lit::from_parts(affected, new))
-                .next(),
-            None
-        );
+        debug_assert!(self
+            .implications
+            .direct_implications_of(Lit::from_parts(affected, new))
+            .next()
+            .is_none());
 
         let new_bound = Lit::from_parts(affected, new);
 
@@ -283,7 +311,7 @@ impl Domains {
                     let lit = ev.new_literal();
                     // invariant: variables in implications are not optional
                     debug_assert_eq!(self.presence(lit.variable()), Lit::TRUE);
-                    for implied in self.presence_graph.direct_implications_of(lit) {
+                    for implied in self.implications.direct_implications_of(lit) {
                         self.doms.set_bound(
                             implied.affected_bound(),
                             implied.bound_value(),
@@ -415,23 +443,29 @@ impl Domains {
 
         loop {
             for l in explanation.lits.drain(..) {
-                debug_assert!(self.entails(l));
-                // find the location of the event that made it true
-                // if there is no such event, it means that the literal is implied in the initial state and we can ignore it
-                if let Some(loc) = self.implying_event(l) {
-                    match self.trail().decision_level_class(loc) {
-                        DecisionLevelClass::Root => {
-                            // implied at decision level 0, and thus always true, discard it
-                        }
-                        DecisionLevelClass::Current => {
-                            // at the current decision level, add to the queue
-                            self.queue.push(InQueueLit { cause: loc, lit: l })
-                        }
-                        DecisionLevelClass::Intermediate => {
-                            // implied before the current decision level, the negation of the literal will appear in the final clause (1UIP)
-                            result.push(!l)
+                if self.entails(l) {
+                    // find the location of the event that made it true
+                    // if there is no such event, it means that the literal is implied in the initial state and we can ignore it
+                    if let Some(loc) = self.implying_event(l) {
+                        match self.trail().decision_level_class(loc) {
+                            DecisionLevelClass::Root => {
+                                // implied at decision level 0, and thus always true, discard it
+                            }
+                            DecisionLevelClass::Current => {
+                                // at the current decision level, add to the queue
+                                self.queue.push(InQueueLit { cause: loc, lit: l })
+                            }
+                            DecisionLevelClass::Intermediate => {
+                                // implied before the current decision level, the negation of the literal will appear in the final clause (1UIP)
+                                result.push(!l)
+                            }
                         }
                     }
+                } else {
+                    // the event is not entailed, must be part of an eager propagation
+                    // Even if it was not necessary for this propagation to occur, it must be part of
+                    // the clause for correctness
+                    result.push(!l)
                 }
             }
             debug_assert!(explanation.lits.is_empty());
@@ -584,11 +618,9 @@ impl PartialOrd for InQueueLit {
 
 #[cfg(test)]
 mod tests {
-    use crate::bounds::Lit;
-    use crate::lang::VarRef;
     use crate::state::domains::Domains;
-    use crate::state::{Cause, Explainer, Explanation, InferenceCause, InvalidUpdate, OptDomain, Origin};
-    use crate::WriterId;
+    use crate::state::*;
+    use crate::*;
     use aries_backtrack::Backtrack;
     use std::collections::HashSet;
 
diff --git a/model/src/state/event.rs b/core/src/state/event.rs
similarity index 92%
rename from model/src/state/event.rs
rename to core/src/state/event.rs
index 9c58d054..f1575d91 100644
--- a/model/src/state/event.rs
+++ b/core/src/state/event.rs
@@ -1,6 +1,5 @@
-use crate::bounds::{BoundValue, Lit, VarBound};
-use crate::state::cause::Origin;
-use crate::state::int_domains::ValueCause;
+use crate::state::{Origin, ValueCause};
+use crate::*;
 use aries_backtrack::EventIndex;
 
 pub type ChangeIndex = Option<EventIndex>;
diff --git a/model/src/state/explanation.rs b/core/src/state/explanation.rs
similarity index 94%
rename from model/src/state/explanation.rs
rename to core/src/state/explanation.rs
index 85761fdc..7c671480 100644
--- a/model/src/state/explanation.rs
+++ b/core/src/state/explanation.rs
@@ -1,10 +1,10 @@
-use crate::bounds::Lit;
 use crate::state::{Domains, InferenceCause};
+use crate::Lit;
 
 /// Builder for a conjunction of literals that make the explained literal true
 #[derive(Clone, Debug)]
 pub struct Explanation {
-    pub(crate) lits: Vec<Lit>,
+    pub lits: Vec<Lit>,
 }
 impl Explanation {
     pub fn new() -> Self {
diff --git a/model/src/state/int_domains.rs b/core/src/state/int_domains.rs
similarity index 97%
rename from model/src/state/int_domains.rs
rename to core/src/state/int_domains.rs
index d74c0e46..70cc343f 100644
--- a/model/src/state/int_domains.rs
+++ b/core/src/state/int_domains.rs
@@ -1,8 +1,7 @@
-use crate::bounds::{BoundValue, Lit, VarBound};
-use crate::lang::{IntCst, VarRef};
 use crate::state::cause::Origin;
 use crate::state::event::{ChangeIndex, Event};
 use crate::state::InvalidUpdate;
+use crate::*;
 use aries_backtrack::{Backtrack, BacktrackWith, DecLvl, EventIndex, ObsTrail};
 use aries_collections::ref_store::RefVec;
 
@@ -24,14 +23,14 @@ impl ValueCause {
     }
 }
 
-/// Associates every variable to the bounds of its domain.
-/// In addition, it maintains the history of changes that caused the bounds to be in this state,
+/// Associates every variable to the literals of its domain.
+/// In addition, it maintains the history of changes that caused the literals to be in this state,
 /// which enables explanations and backtracking.
 ///
 /// **Invariant:** every domain is non empty. Hence any update that would result in an empty domain
 /// would return an `Error<EmptyDomain>`.
 #[derive(Clone)]
-pub(crate) struct IntDomains {
+pub struct IntDomains {
     /// Associates each lb/ub of each variable to its current value and the event that caused the latest update.
     bounds: RefVec<VarBound, ValueCause>,
     /// All events that updated the bound values.
diff --git a/model/src/bounds/var_bound.rs b/core/src/var_bound.rs
similarity index 79%
rename from model/src/bounds/var_bound.rs
rename to core/src/var_bound.rs
index 137acc60..1eff2a70 100644
--- a/model/src/bounds/var_bound.rs
+++ b/core/src/var_bound.rs
@@ -1,13 +1,12 @@
-use crate::bounds::{BoundValue, Lit, Relation};
-use crate::lang::VarRef;
+use crate::*;
 
 /// Represents the upper or the lower bound of a particular variable.
-/// The type has dense integer values and can by used an index in an array.
+/// The type has dense integer values and can be used as an index in an array.
 ///
 /// It is coded on 32 bits where:
 ///  - the 31 most significant bits represent the variable
 ///  - the least significant bit represents either a lower bound (0) or upper bound (1).
-#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq)]
+#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash)]
 pub struct VarBound(u32);
 
 impl VarBound {
@@ -28,13 +27,13 @@ impl VarBound {
     }
 
     #[inline]
-    pub fn ub(v: VarRef) -> Self {
-        VarBound((u32::from(v) << 1) + 1)
+    pub const fn ub(v: VarRef) -> Self {
+        VarBound((v.to_u32() << 1) + 1)
     }
 
     #[inline]
-    pub fn lb(v: VarRef) -> Self {
-        VarBound(u32::from(v) << 1)
+    pub const fn lb(v: VarRef) -> Self {
+        VarBound(v.to_u32() << 1)
     }
 
     #[inline]
@@ -48,8 +47,7 @@ impl VarBound {
     /// and vice versa.
     ///
     /// ```
-    /// use aries_model::bounds::VarBound;
-    /// use aries_model::lang::VarRef;
+    /// use aries_core::*;
     /// let var = VarRef::from(1u32);
     /// let var_lb = VarBound::lb(var);
     /// let var_ub = VarBound::ub(var);
@@ -57,17 +55,17 @@ impl VarBound {
     /// assert_eq!(var_ub.symmetric_bound(), var_lb);
     /// ```
     #[inline]
-    pub fn symmetric_bound(self) -> Self {
+    pub const fn symmetric_bound(self) -> Self {
         VarBound(self.0 ^ 0x1)
     }
 
     #[inline]
-    pub fn is_lb(self) -> bool {
+    pub const fn is_lb(self) -> bool {
         (self.0 & 0x1) == 0
     }
 
     #[inline]
-    pub fn is_ub(self) -> bool {
+    pub const fn is_ub(self) -> bool {
         (self.0 & 0x1) == 1
     }
 
diff --git a/core/src/variable.rs b/core/src/variable.rs
new file mode 100644
index 00000000..bec49236
--- /dev/null
+++ b/core/src/variable.rs
@@ -0,0 +1,37 @@
+use crate::Lit;
+use aries_collections::create_ref_type;
+use std::hash::Hash;
+
+/// Type representing an integer constant.
+pub type IntCst = i32;
+
+/// Overflow tolerant min value for integer constants.
+/// It is used as a default for the lower bound of integer variable domains
+pub const INT_CST_MIN: IntCst = IntCst::MIN / 2 + 1;
+
+/// Overflow tolerant max value for integer constants.
+/// It is used as a default for the upper bound of integer variable domains
+pub const INT_CST_MAX: IntCst = IntCst::MAX / 2 - 1;
+
+create_ref_type!(VarRef);
+
+impl VarRef {
+    /// A reserved special variable that is always equal to 0. It corresponds to the first representable VarRef.
+    ///
+    /// For efficiency reasons, this special case is not treated separately from the other variables, and it is the responsibility
+    /// of the producers of VarRef to ensure that they only emit this value for variables whose domain is `[0,0]`.
+    pub const ZERO: VarRef = VarRef::from_u32(0);
+
+    pub fn leq(self, i: IntCst) -> Lit {
+        Lit::leq(self, i)
+    }
+    pub fn lt(self, i: IntCst) -> Lit {
+        Lit::lt(self, i)
+    }
+    pub fn geq(self, i: IntCst) -> Lit {
+        Lit::geq(self, i)
+    }
+    pub fn gt(self, i: IntCst) -> Lit {
+        Lit::gt(self, i)
+    }
+}
diff --git a/doc/book.toml b/doc/book.toml
index a4d66599..857acc0a 100644
--- a/doc/book.toml
+++ b/doc/book.toml
@@ -4,3 +4,6 @@ language = "en"
 multilingual = false
 src = "src"
 title = "Aries Documentation"
+
+[output.html]
+mathjax-support = true
\ No newline at end of file
diff --git a/doc/src/SUMMARY.md b/doc/src/SUMMARY.md
index a9f81329..ba0b925c 100644
--- a/doc/src/SUMMARY.md
+++ b/doc/src/SUMMARY.md
@@ -6,13 +6,20 @@
 - [Literals](./literals.md)
 - [Domains]()
 - [Optional Variables](optional-variables.md)
+- [Implication Graph]()
 - [Reasoners]()
 - [Explanations]()
 
+# Optionals
+
+- [Propagation](./optionals/propagation.md)
+- [STN](./optionals/stn.md)
+
 # SAT Reasoning
 
 
 - [Unit Propagation]()
+- [Extension for optionals](./sat/optional-reasoning.md)
 
 # STN / Difference Logic 
 
diff --git a/doc/src/literals.md b/doc/src/literals.md
index 27631b3d..c1c293a3 100644
--- a/doc/src/literals.md
+++ b/doc/src/literals.md
@@ -3,7 +3,7 @@
 
 A literal is a boolean assertion on the bounds of a variable.
 
-Given a variable `var` and an integer value `val, a literal is either:
+Given a variable `var` and an integer value `val`, a literal is either:
 
  - `var <= val`, or
  - `var > val`
diff --git a/doc/src/sat/optional-reasoning.md b/doc/src/sat/optional-reasoning.md
new file mode 100644
index 00000000..cd812c95
--- /dev/null
+++ b/doc/src/sat/optional-reasoning.md
@@ -0,0 +1,101 @@
+# Extension for optionals
+
+Consider the case where be have to optional boolean variables:
+
+ - `a` that is in scope α
+ - `b` that is in scope β
+
+ The boolean disjunction \\( a \vee b \\) is:
+
+ - `true` if `a` is present and true or if `b` is present and true   `(a ∧ α) ∨ (b ∧ β)`
+ - `false` if `a` is present and false and if `b` is  present and false `(¬a ∧ α) ∨ (¬b ∧ β)`
+ - `⊥` (undefined) if `(¬a ∧ ¬β) ∨ (¬b ∧ ¬α) ∨ (¬α ∧ ¬β)`
+
+
+\\[ (a_\alpha \vee b_\beta)_{\alpha\beta} \\]
+
+
+
+# Sufficient conditions for propagation
+
+## Simple case
+
+If both `a` and `b` are present (i.e. we are in the conjunctive scope αβ), then all components of the clause are present and be can 
+
+
+
+!((a && x) || (b && y)) || !( !a && x && !b && y) 
+
+
+
+## WIP
+
+We have a disjunction `¬a ∨ ¬b ∨ c`, a model `M`.
+
+We say that `M ⊨ a` if `a ∈ { true, ⊥}`.
+
+- If `M ⊨ a, b`, can be propagate `c ≠ false` ?
+
+  - if `α ∧ β`, meaning that `a = true`  and `b = true`    (base condition in non-optional reasoning)
+  - if `¬γ`, meaning that `c = ⊥`   (condition where propagation is correct but useless)
+
+
+
+If we know that `(α ∧ β) ∨ ¬γ` holds, we can always propagate because we known that we will always be in one of the two valid propagation conditions.
+Let's reformulate it:
+
+ - `γ ⇒ (α ∧ β)`
+
+
+ ### With validity scope
+
+ Lets call `π` the conjunctive scope of all conjuncts: `π = (α ∧ β ∧ γ)`.
+
+ **Claim:** When unit propagation infers that a literal `a` of the clause is true, it makes this inference under the assumption that the the clause is within its validity scope.
+ In other words, it infers that `π ⇒ a`.
+
+ By construction, the literal `π` is always present.
+ When I learn `π ⇒ a`, when can i deduce that `l ≠ false` (i.e. `l ∈ { ⊥, true}`).
+
+ - when `α ⇒ π`, proof:
+   - assume `¬α` then asserting `l ≠ false` has no effect.
+   - assume `α` then the premise `π` will hold.
+
+
+
+
+ ## Example
+
+ We have an action A with presence `α` and parameter `xa`. It can be supported 
+ 
+  - by action B with presence `β` and parameter `xb`
+    - \\( subB_α  ⇔  [β ∧ (xa_α = xb_β)_{αβ}]_α   \\)
+  - by action C with presence `γ` and parameter `xc`
+
+
+Constraints:
+
+- \\( ¬α ∨ subB_α ∨ subC_α \\)
+  - `¬subB ∧ ¬subC ⇒ ¬α` ✓ (but not the same conditions)
+  - `α ∧ ¬subB ⇒ subC` ✓
+  - `α ∧ ¬subC ⇒ subB` ✓
+- \\( subB_α  ⇔  [β ∧ (xa_α = xb_β)_{αβ}]_α   \\)
+  - `¬subB ∨ β`
+    - `subB ⇒ β` x    |> should be `subB ⇒ (α ⇒ β)`
+    - `¬β ⇒ ¬subB` ✓
+  - `¬subB ∨ (xa=xb)`
+    - `subB ⇒ (xa=xb)` ✓
+    - `(xa≠xb) → ¬subB` ✓
+  - `(xa≠xb) ∨ ¬β ∨ subB`
+    - `(xa=xb) ∧ β ⇒ subB` ✓
+    - `(xa=xb) ∧ ¬subB ⇒ ¬β` x   |> should be `... ⇒ (α ⇒ ¬β)`
+    - `β ∧ ¬subB ⇒ (xa≠xb)` ✓
+- subC ...
+
+Propagation desiderata:
+
+ - if we learn \\( (xa_α = xb_β) ≠ true \\)
+   - infer `subB ≠ true`
+   - infer `subC ≠ false`
+
+- if we learn  
\ No newline at end of file
diff --git a/doc/src/variables.md b/doc/src/variables.md
index 1aba27db..da145f98 100644
--- a/doc/src/variables.md
+++ b/doc/src/variables.md
@@ -12,14 +12,14 @@ A `VarRef` type is the ID of a variable.
 
 ## The Zero variable
 
-Among the possible instances of a `VarRef`, the `VarRef::ZERO` constant is a reserved value.
+Among the possible instances of a `VarRef`, the `VarRef::ZERO` constant is a reserved value that is always equal to `0`.
 
 
 # Creating Variables
 
 ## Manually creating variables.
 
-IT is possible to manually create variables by using `VarRef::from_u32(i)` where is an `u32`. This is typically non what you want to do but allows:
+IT is possible to manually create variables by using `VarRef::from_u32(i)` where `i` an unsigned integer on 32 bits (`u32`). This is typically non what you want to do but allows:
 
  - converting an index to a variable (this what is used internally in the implementation of aries)
  - create variable instances for testing purposes.
@@ -27,7 +27,7 @@ IT is possible to manually create variables by using `VarRef::from_u32(i)` where
 The first variable is reserved as the `VarRef::ZERO` value.
 
 ```rust
-assert_eq!(VarRef::from_u32(0), VarRef::ZERO)
+assert_eq!(VarRef::from_u32(0), VarRef::ZERO);
 ``` 
 
 
diff --git a/examples/jobshop/Cargo.toml b/examples/jobshop/Cargo.toml
index ed90246e..be6943b2 100644
--- a/examples/jobshop/Cargo.toml
+++ b/examples/jobshop/Cargo.toml
@@ -8,6 +8,7 @@ edition = "2018"
 structopt = "0.3"
 anyhow = "1.0"
 aries_backtrack = { path = "../../backtrack" }
+aries_core = { path = "../../core" }
 aries_model = { path = "../../model" }
 aries_solver = { path = "../../solver" }
 aries_tnet = { path = "../../tnet" }
diff --git a/examples/jobshop/main.rs b/examples/jobshop/main.rs
index a4756055..5a16c14a 100644
--- a/examples/jobshop/main.rs
+++ b/examples/jobshop/main.rs
@@ -1,8 +1,8 @@
 use aries_backtrack::{Backtrack, DecLvl};
-use aries_model::bounds::Lit;
+use aries_core::*;
 use aries_model::extensions::{AssignmentExt, Shaped};
 use aries_model::lang::expr::leq;
-use aries_model::lang::{IVar, VarRef};
+use aries_model::lang::IVar;
 use aries_solver::solver::search::activity::{ActivityBrancher, Heuristic};
 use aries_solver::solver::search::{Decision, SearchControl};
 use aries_solver::solver::stats::Stats;
@@ -21,6 +21,12 @@ pub enum Var {
     Prec(usize, usize, usize, usize),
 }
 
+impl std::fmt::Display for Var {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        write!(f, "{:?}", self)
+    }
+}
+
 type Model = aries_model::Model<Var>;
 type Solver = aries_solver::solver::Solver<Var>;
 type ParSolver = aries_solver::parallel_solver::ParSolver<Var>;
diff --git a/examples/sat/Cargo.toml b/examples/sat/Cargo.toml
index d5953eea..a8e98860 100644
--- a/examples/sat/Cargo.toml
+++ b/examples/sat/Cargo.toml
@@ -7,6 +7,7 @@ edition = "2018"
 [dependencies]
 structopt = "0.3"
 anyhow = "1.0"
+aries_core = { path = "../../core" }
 aries_model = { path = "../../model" }
 aries_solver = { path = "../../solver" }
 varisat-dimacs = "0.2.2"
diff --git a/examples/sat/main.rs b/examples/sat/main.rs
index 4c7def0a..bbdef1ff 100644
--- a/examples/sat/main.rs
+++ b/examples/sat/main.rs
@@ -1,7 +1,7 @@
 #![allow(clippy::map_entry)]
 
 use anyhow::*;
-use aries_model::bounds::Lit;
+use aries_core::Lit;
 use aries_model::lang::expr::or;
 use aries_solver::parallel_solver::ParSolver;
 use aries_solver::solver::search::activity::{ActivityBrancher, BranchingParams};
diff --git a/examples/smt/Cargo.toml b/examples/smt/Cargo.toml
index 795d8e7f..f1171298 100644
--- a/examples/smt/Cargo.toml
+++ b/examples/smt/Cargo.toml
@@ -6,7 +6,9 @@ edition = "2018"
 
 [dependencies]
 aries_backtrack = { path = "../../backtrack" }
+aries_core = { path = "../../core" }
 aries_model = { path = "../../model" }
 aries_solver = { path = "../../solver" }
 aries_tnet = { path = "../../tnet" }
+itertools = "0.10.3"
 
diff --git a/model/Cargo.toml b/model/Cargo.toml
index aa755333..2811ce1d 100644
--- a/model/Cargo.toml
+++ b/model/Cargo.toml
@@ -9,26 +9,10 @@ edition = "2018"
 [dependencies]
 anyhow = "1.0.35"
 streaming-iterator = "0.1.5"
+aries_core = { path = "../core" }
 aries_backtrack = { path = "../backtrack" }
 aries_collections = { path = "../collections" }
 aries_utils = { path = "../utils" }
 lazy_static = "1.4.0"
 
 
-[dev-dependencies]
-criterion = "0.3"
-rand = "0.8.3"
-
-
-[[bench]]
-name = "bounds"
-harness = false
-
-[[bench]]
-name = "watches"
-harness = false
-
-
-[[bench]]
-name = "domains"
-harness = false
\ No newline at end of file
diff --git a/model/src/bounds.rs b/model/src/bounds.rs
deleted file mode 100644
index 4c2ab7cd..00000000
--- a/model/src/bounds.rs
+++ /dev/null
@@ -1,11 +0,0 @@
-mod bound;
-mod bound_value;
-mod disjunction;
-mod var_bound;
-mod watches;
-
-pub use bound::*;
-pub use bound_value::*;
-pub use disjunction::*;
-pub use var_bound::VarBound;
-pub use watches::*;
diff --git a/model/src/decomposition.rs b/model/src/decomposition.rs
index e098508e..ecc75089 100644
--- a/model/src/decomposition.rs
+++ b/model/src/decomposition.rs
@@ -1,11 +1,11 @@
-use crate::bounds::Lit;
 use crate::extensions::AssignmentExt;
 use crate::label::Label;
-use crate::lang::expr::{and, leq, opt_leq};
-use crate::lang::normal_form::{NFEq, NFOptEq};
+use crate::lang::expr::{and, leq};
+use crate::lang::normal_form::NFEq;
 use crate::lang::reification::{downcast, BindTarget, BindingCursor, Expr};
 use crate::lang::IVar;
 use crate::Model;
+use aries_core::*;
 use std::sync::Arc;
 
 /// Module to constructs the constraints from the model and store them in a queue.
@@ -95,20 +95,15 @@ impl<Lbl: Label> Decompose<Lbl> for Eq2Leq {
             // decompose `l <=> (a = b)` into `l <=> (a <= b) && (b <= a)`
             let lhs = IVar::new(lhs);
             let rhs = IVar::new(rhs) + rhs_add;
-            let x = model.reify(leq(lhs, rhs));
-            let y = model.reify(leq(rhs, lhs));
-            model.bind(and([x, y]), literal);
-            DecompositionResult::Decomposed
-        } else if let Some(&NFOptEq { lhs, rhs, rhs_add }) = downcast(expr) {
             if model.entails(literal) {
-                let a = IVar::new(lhs);
-                let b = IVar::new(rhs) + rhs_add;
-                model.enforce(opt_leq(a, b));
-                model.enforce(opt_leq(b, a));
-                DecompositionResult::Decomposed
+                model.bind(leq(lhs, rhs), literal);
+                model.bind(leq(rhs, lhs), literal);
             } else {
-                DecompositionResult::Inapplicable
+                let x = model.reify(leq(lhs, rhs));
+                let y = model.reify(leq(rhs, lhs));
+                model.bind(and([x, y]), literal);
             }
+            DecompositionResult::Decomposed
         } else {
             DecompositionResult::Inapplicable
         }
diff --git a/model/src/extensions/assignments.rs b/model/src/extensions/assignments.rs
index 85d38851..5bcb0ee7 100644
--- a/model/src/extensions/assignments.rs
+++ b/model/src/extensions/assignments.rs
@@ -1,9 +1,9 @@
-use crate::bounds::Lit;
 use crate::extensions::SavedAssignment;
-use crate::lang::{Atom, FAtom, IAtom, IVar, IntCst, SAtom, VarRef};
-use crate::state::{FixedDomain, IntDomain, OptDomain};
+use crate::lang::{Atom, FAtom, IAtom, IVar, SAtom};
 use crate::symbols::ContiguousSymbols;
 use crate::symbols::SymId;
+use aries_core::state::{FixedDomain, IntDomain, OptDomain};
+use aries_core::*;
 
 /// Extension methods for an object containing a partial or total assignment to a problem.
 pub trait AssignmentExt {
diff --git a/model/src/extensions/disjunction.rs b/model/src/extensions/disjunction.rs
index a178b459..cda606b2 100644
--- a/model/src/extensions/disjunction.rs
+++ b/model/src/extensions/disjunction.rs
@@ -1,5 +1,5 @@
-use crate::bounds::Lit;
 use crate::extensions::PartialAssignment;
+use aries_core::*;
 
 /// Extension trait that provides convenience methods to query the status of disjunctions.
 pub trait DisjunctionExt<Disj>
diff --git a/model/src/extensions/format.rs b/model/src/extensions/format.rs
index 02584937..7c545702 100644
--- a/model/src/extensions/format.rs
+++ b/model/src/extensions/format.rs
@@ -1,9 +1,9 @@
-use crate::bounds::Lit;
 use crate::label::Label;
-use crate::lang::{Atom, FAtom, IAtom, IVar, Kind, SAtom, VarRef};
+use crate::lang::{Atom, FAtom, IAtom, IVar, Kind, SAtom, Type};
 use crate::symbols::{SymId, SymbolTable};
 use crate::types::TypeId;
 use crate::ModelShape;
+use aries_core::*;
 use aries_utils::input::Sym;
 use aries_utils::Fmt;
 
@@ -17,6 +17,10 @@ where
         self.get_shape().labels.get(var.into())
     }
 
+    fn get_type(&self, var: impl Into<VarRef>) -> Option<Type> {
+        self.get_shape().types.get(var.into()).copied()
+    }
+
     fn get_var(&self, label: &Lbl) -> Option<VarRef> {
         self.get_shape().get_variable(label)
     }
@@ -68,8 +72,25 @@ fn format_impl<Lbl: Label>(ctx: &impl Shaped<Lbl>, atom: Atom, f: &mut std::fmt:
     }
 }
 fn format_impl_bool<Lbl: Label>(ctx: &impl Shaped<Lbl>, b: Lit, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-    format_impl_var(ctx, b.variable(), Kind::Int, f)?;
-    write!(f, " {} {}", b.relation(), b.value())
+    let tpe = ctx.get_type(b.variable());
+    let t = b.variable().geq(1);
+    if b == Lit::TRUE {
+        write!(f, "true")
+    } else if b == Lit::FALSE {
+        write!(f, "false")
+    } else if let Some(Type::Bool) = tpe {
+        if b == t {
+            format_impl_var(ctx, b.variable(), Kind::Bool, f)
+        } else {
+            debug_assert_eq!(b, !t);
+            write!(f, "!")?;
+            format_impl_var(ctx, b.variable(), Kind::Bool, f)
+        }
+    } else {
+        let tpe = tpe.unwrap_or(Type::Int);
+        format_impl_var(ctx, b.variable(), tpe.into(), f)?;
+        write!(f, " {} {}", b.relation(), b.value())
+    }
 }
 
 #[allow(clippy::comparison_chain)]
@@ -120,7 +141,7 @@ fn format_impl_var<Lbl: Label>(
     f: &mut std::fmt::Formatter<'_>,
 ) -> std::fmt::Result {
     if let Some(lbl) = ctx.get_label(v) {
-        write!(f, "{:?}", lbl)
+        write!(f, "{}", lbl)
     } else {
         let prefix = match kind {
             Kind::Bool => "b_",
diff --git a/model/src/extensions/mod.rs b/model/src/extensions/mod.rs
index 9410564d..257e38e0 100644
--- a/model/src/extensions/mod.rs
+++ b/model/src/extensions/mod.rs
@@ -13,10 +13,9 @@ pub use assignments::*;
 pub use disjunction::*;
 pub use format::*;
 
-use crate::bounds::Lit;
-use crate::lang::{IAtom, VarRef};
-use crate::state::{Domains, IntDomain};
-use crate::Model;
+use crate::lang::IAtom;
+use aries_core::state::{Domains, IntDomain};
+use aries_core::*;
 
 pub trait PartialAssignment {
     fn entails(&self, literal: Lit) -> bool;
@@ -39,12 +38,6 @@ impl PartialAssignment for Domains {
 
 pub type SavedAssignment = Domains;
 
-impl SavedAssignment {
-    pub fn from_model<Lbl>(model: &Model<Lbl>) -> SavedAssignment {
-        model.state.clone()
-    }
-}
-
 impl AssignmentExt for SavedAssignment {
     fn entails(&self, literal: Lit) -> bool {
         self.entails(literal)
diff --git a/model/src/label.rs b/model/src/label.rs
index 9e14dc37..83c3bc3d 100644
--- a/model/src/label.rs
+++ b/model/src/label.rs
@@ -1,13 +1,13 @@
-use crate::lang::VarRef;
 use aries_collections::ref_store::RefMap;
+use aries_core::*;
 use std::collections::HashMap;
-use std::fmt::Debug;
+use std::fmt::{Debug, Display};
 use std::hash::Hash;
 use std::sync::Arc;
 
 /// Trait requiring the minimum capabilities for a type to serve as the label of variables.
-pub trait Label: Debug + Clone + Eq + PartialEq + Hash + Send + Sync + 'static {}
-impl<T> Label for T where T: Debug + Clone + Eq + PartialEq + Hash + Send + Sync + 'static {}
+pub trait Label: Display + Debug + Clone + Eq + PartialEq + Hash + Send + Sync + 'static {}
+impl<T> Label for T where T: Display + Debug + Clone + Eq + PartialEq + Hash + Send + Sync + 'static {}
 
 #[derive(Clone)]
 pub struct VariableLabels<Lbl> {
diff --git a/model/src/lang.rs b/model/src/lang.rs
index b3c54c3c..373e4b6e 100644
--- a/model/src/lang.rs
+++ b/model/src/lang.rs
@@ -6,62 +6,40 @@ mod int;
 pub mod normal_form;
 pub mod reification;
 mod sym;
+mod validity_scope;
 mod variables;
 
 pub use atom::Atom;
 pub use boolean::BVar;
-pub use fixed::*;
+pub use fixed::{FAtom, FVar};
 pub use int::{IAtom, IVar};
-pub use sym::{SAtom, SVar};
-pub use variables::Variable;
+pub use validity_scope::*;
 
-use crate::bounds::Lit;
 use crate::types::TypeId;
-use aries_collections::create_ref_type;
-use std::convert::TryFrom;
-use std::hash::Hash;
-
-/// Type representing an integer constant.
-pub type IntCst = i32;
-
-/// Overflow tolerant min value for integer constants.
-/// It is used as a default for the lower bound of integer variable domains
-pub const INT_CST_MIN: IntCst = IntCst::MIN / 2 + 1;
-
-/// Overflow tolerant max value for integer constants.
-/// It is used as a default for the upper bound of integer variable domains
-pub const INT_CST_MAX: IntCst = IntCst::MAX / 2 - 1;
-
-create_ref_type!(VarRef);
-
-impl VarRef {
-    /// A reserved special variable that is always equal to 0. It corresponds to the first representable VarRef.
-    ///
-    /// For efficiency reasons, this special case is not treated separately from the other variables, and it is the responsibility
-    /// of the producers of VarRef to ensure that they only emit this value for variables whose domain is `[0,0]`.
-    pub const ZERO: VarRef = VarRef::from_u32(0);
-
-    pub fn leq(self, i: IntCst) -> Lit {
-        Lit::leq(self, i)
-    }
-    pub fn lt(self, i: IntCst) -> Lit {
-        Lit::lt(self, i)
-    }
-    pub fn geq(self, i: IntCst) -> Lit {
-        Lit::geq(self, i)
-    }
-    pub fn gt(self, i: IntCst) -> Lit {
-        Lit::gt(self, i)
-    }
-}
+use aries_core::IntCst;
+pub use sym::{SAtom, SVar};
+pub use variables::Variable;
 
 #[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq)]
 pub enum Type {
     Sym(TypeId),
     Int,
+    /// A fixed-point numeral, parameterized with its denominator.
+    Fixed(IntCst),
     Bool,
 }
 
+impl From<Type> for Kind {
+    fn from(tpe: Type) -> Self {
+        match tpe {
+            Type::Sym(_) => Kind::Sym,
+            Type::Int => Kind::Int,
+            Type::Fixed(denum) => Kind::Fixed(denum),
+            Type::Bool => Kind::Bool,
+        }
+    }
+}
+
 #[derive(Hash, Ord, PartialOrd, Eq, PartialEq, Copy, Clone, Debug)]
 pub enum Kind {
     Bool,
diff --git a/model/src/lang/atom.rs b/model/src/lang/atom.rs
index 705bc382..eb9a58a5 100644
--- a/model/src/lang/atom.rs
+++ b/model/src/lang/atom.rs
@@ -1,6 +1,6 @@
 use super::*;
-use crate::bounds::Lit;
 use crate::symbols::TypedSym;
+use aries_core::*;
 
 #[derive(Hash, Eq, PartialEq, Copy, Clone, Debug)]
 pub enum Atom {
@@ -80,6 +80,24 @@ impl TryFrom<Atom> for Lit {
     }
 }
 
+impl TryFrom<Atom> for bool {
+    type Error = ConversionError;
+
+    fn try_from(value: Atom) -> Result<Self, Self::Error> {
+        match value {
+            Atom::Bool(Lit::TRUE) => Ok(true),
+            Atom::Bool(Lit::FALSE) => Ok(false),
+            _ => Err(ConversionError::TypeError),
+        }
+    }
+}
+
+impl From<bool> for Atom {
+    fn from(b: bool) -> Self {
+        Atom::Bool(if b { Lit::TRUE } else { Lit::FALSE })
+    }
+}
+
 impl TryFrom<Atom> for IAtom {
     type Error = ConversionError;
 
@@ -128,10 +146,9 @@ impl TryFrom<Atom> for Variable {
 }
 
 use crate::transitive_conversions;
-use std::convert::TryInto;
+use std::convert::{TryFrom, TryInto};
 
 transitive_conversions!(Atom, IAtom, IVar);
 transitive_conversions!(Atom, IAtom, IntCst);
 transitive_conversions!(Atom, SAtom, SVar);
 transitive_conversions!(Atom, SAtom, TypedSym);
-transitive_conversions!(Atom, Lit, bool);
diff --git a/model/src/lang/boolean.rs b/model/src/lang/boolean.rs
index 06dc5e1e..8ad47edf 100644
--- a/model/src/lang/boolean.rs
+++ b/model/src/lang/boolean.rs
@@ -1,5 +1,5 @@
-use crate::bounds::Lit;
-use crate::lang::{IVar, VarRef};
+use crate::lang::IVar;
+use aries_core::*;
 
 /// A boolean variable.
 /// It is a wrapper around an (untyped) discrete variable to provide type safety.
@@ -31,6 +31,12 @@ impl From<BVar> for VarRef {
     }
 }
 
+impl From<BVar> for Lit {
+    fn from(v: BVar) -> Self {
+        v.true_lit()
+    }
+}
+
 impl From<usize> for BVar {
     fn from(i: usize) -> Self {
         BVar(VarRef::from(i))
diff --git a/model/src/lang/expr.rs b/model/src/lang/expr.rs
index b8ca510f..10ae361f 100644
--- a/model/src/lang/expr.rs
+++ b/model/src/lang/expr.rs
@@ -1,7 +1,8 @@
-use crate::bounds::{Disjunction, Lit};
-use crate::lang::normal_form::{NFEq, NFLeq, NFOptEq, NFOptLeq, NormalExpr};
-use crate::lang::reification::ReifiableExpr;
-use crate::lang::{Atom, FAtom, IAtom};
+use crate::lang::normal_form::{NFEq, NFLeq, NormalExpr};
+use crate::lang::reification::{ExprInterface, ReifiableExpr};
+use crate::lang::{Atom, FAtom, IAtom, ValidityScope};
+use aries_core::literals::Disjunction;
+use aries_core::*;
 
 /// Trait denoting the capability of transforming an expression into its normal form.
 ///
@@ -29,6 +30,12 @@ impl Normalize<Never> for Lit {
 #[derive(Eq, PartialEq, Hash, Debug)]
 pub struct Never(());
 
+impl ExprInterface for Never {
+    fn validity_scope(&self, _: &dyn Fn(VarRef) -> Lit) -> ValidityScope {
+        unreachable!()
+    }
+}
+
 pub fn leq(lhs: impl Into<IAtom>, rhs: impl Into<IAtom>) -> Leq {
     Leq(lhs.into(), rhs.into())
 }
@@ -44,10 +51,6 @@ pub fn gt(lhs: impl Into<IAtom>, rhs: impl Into<IAtom>) -> Leq {
     lt(rhs, lhs)
 }
 
-pub fn opt_leq(lhs: impl Into<IAtom>, rhs: impl Into<IAtom>) -> OptLeq {
-    OptLeq(lhs.into(), rhs.into())
-}
-
 pub fn f_leq(lhs: impl Into<FAtom>, rhs: impl Into<FAtom>) -> Leq {
     let lhs = lhs.into();
     let rhs = rhs.into();
@@ -60,12 +63,6 @@ pub fn f_lt(lhs: impl Into<FAtom>, rhs: impl Into<FAtom>) -> Leq {
     assert_eq!(lhs.denom, rhs.denom);
     lt(lhs.num, rhs.num)
 }
-pub fn f_opt_leq(lhs: impl Into<FAtom>, rhs: impl Into<FAtom>) -> OptLeq {
-    let lhs = lhs.into();
-    let rhs = rhs.into();
-    assert_eq!(lhs.denom, rhs.denom);
-    opt_leq(lhs.num, rhs.num)
-}
 
 pub fn eq(lhs: impl Into<Atom>, rhs: impl Into<Atom>) -> Eq {
     let lhs = lhs.into();
@@ -81,13 +78,6 @@ pub fn neq(lhs: impl Into<Atom>, rhs: impl Into<Atom>) -> Neq {
     Neq(lhs, rhs)
 }
 
-pub fn opt_eq(lhs: impl Into<Atom>, rhs: impl Into<Atom>) -> OptEq {
-    let lhs = lhs.into();
-    let rhs = rhs.into();
-    assert_eq!(lhs.kind(), rhs.kind());
-    OptEq(lhs, rhs)
-}
-
 pub fn or(disjuncts: impl Into<Box<[Lit]>>) -> Or {
     Or(disjuncts.into())
 }
@@ -160,28 +150,3 @@ impl Normalize<NFEq> for Neq {
         !NFEq::eq(self.0, self.1)
     }
 }
-
-pub struct OptLeq(IAtom, IAtom);
-
-impl Normalize<NFOptLeq> for OptLeq {
-    fn normalize(&self) -> NormalExpr<NFOptLeq> {
-        let lhs = self.0.var.into();
-        let rhs = self.1.var.into();
-        let rhs_add = self.1.shift - self.0.shift;
-        NormalExpr::Pos(NFOptLeq { lhs, rhs, rhs_add })
-    }
-}
-
-pub struct OptEq(Atom, Atom);
-
-impl Normalize<NFOptEq> for OptEq {
-    fn normalize(&self) -> NormalExpr<NFOptEq> {
-        assert_eq!(self.0.kind(), self.1.kind());
-        let a = self.0.int_view().unwrap();
-        let b = self.1.int_view().unwrap();
-        let lhs = a.var.into();
-        let rhs = b.var.into();
-        let rhs_add = b.shift - a.shift;
-        NormalExpr::Pos(NFOptEq { lhs, rhs, rhs_add })
-    }
-}
diff --git a/model/src/lang/fixed.rs b/model/src/lang/fixed.rs
index 42505d04..64e9ad11 100644
--- a/model/src/lang/fixed.rs
+++ b/model/src/lang/fixed.rs
@@ -1,4 +1,5 @@
-use crate::lang::{ConversionError, IAtom, IVar, IntCst, VarRef};
+use crate::lang::{ConversionError, IAtom, IVar};
+use aries_core::{IntCst, VarRef};
 use std::cmp::Ordering;
 use std::convert::{TryFrom, TryInto};
 
diff --git a/model/src/lang/int.rs b/model/src/lang/int.rs
index b257b7f8..4bbe45be 100644
--- a/model/src/lang/int.rs
+++ b/model/src/lang/int.rs
@@ -1,5 +1,5 @@
-use crate::bounds::Lit;
-use crate::lang::{ConversionError, IntCst, VarRef};
+use crate::lang::ConversionError;
+use aries_core::*;
 use std::cmp::Ordering;
 use std::convert::TryFrom;
 
@@ -12,6 +12,22 @@ impl IVar {
     pub const fn new(dvar: VarRef) -> Self {
         IVar(dvar)
     }
+
+    pub fn leq(self, i: IntCst) -> Lit {
+        Lit::leq(self, i)
+    }
+
+    pub fn geq(self, i: IntCst) -> Lit {
+        Lit::geq(self, i)
+    }
+
+    pub fn lt(self, i: IntCst) -> Lit {
+        Lit::lt(self, i)
+    }
+
+    pub fn gt(self, i: IntCst) -> Lit {
+        Lit::gt(self, i)
+    }
 }
 
 impl From<IVar> for VarRef {
diff --git a/model/src/lang/normal_form.rs b/model/src/lang/normal_form.rs
index 386334be..21d88406 100644
--- a/model/src/lang/normal_form.rs
+++ b/model/src/lang/normal_form.rs
@@ -1,5 +1,6 @@
-use crate::bounds::Lit;
+use crate::lang::reification::ExprInterface;
 use crate::lang::*;
+use aries_core::*;
 use std::cmp::Ordering;
 
 /// Normal form of an expression with component of type `X`.
@@ -11,6 +12,16 @@ pub enum NormalExpr<X> {
     Neg(X),
 }
 
+impl<X: ExprInterface> NormalExpr<X> {
+    pub fn validity_scope(&self, presence: &dyn Fn(VarRef) -> Lit) -> ValidityScope {
+        match self {
+            NormalExpr::Literal(l) => ValidityScope::new([presence(l.variable())], []),
+            NormalExpr::Pos(x) => x.validity_scope(presence),
+            NormalExpr::Neg(x) => x.validity_scope(presence),
+        }
+    }
+}
+
 impl<T> std::ops::Not for NormalExpr<T> {
     type Output = Self;
 
@@ -69,7 +80,7 @@ impl NFLeq {
 
         // Only encode as a LEQ the patterns with two variables.
         // Other are treated either are constant (if provable as so)
-        // or as bounds on a single variable
+        // or as literals on a single variable
         if lhs == rhs {
             // X  <= X + rhs_add   <=>  0 <= rhs_add
             return (0 <= rhs_add).into();
@@ -95,6 +106,12 @@ impl NFLeq {
     }
 }
 
+impl ExprInterface for NFLeq {
+    fn validity_scope(&self, presence: &dyn Fn(VarRef) -> Lit) -> ValidityScope {
+        ValidityScope::new([presence(self.lhs), presence(self.rhs)], [])
+    }
+}
+
 /// lhs = rhs + rhs_add
 #[derive(Debug, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash)]
 pub struct NFEq {
@@ -145,18 +162,8 @@ impl NFEq {
     }
 }
 
-/// (lhs <= rhs + rhs_add) || absent(lhs) || absent(rhs)
-#[derive(Debug, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash)]
-pub struct NFOptLeq {
-    pub lhs: VarRef,
-    pub rhs: VarRef,
-    pub rhs_add: IntCst,
-}
-
-/// (lhs = rhs + rhs_add) || absent(lhs) || absent(rhs)
-#[derive(Debug, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash)]
-pub struct NFOptEq {
-    pub lhs: VarRef,
-    pub rhs: VarRef,
-    pub rhs_add: IntCst,
+impl ExprInterface for NFEq {
+    fn validity_scope(&self, presence: &dyn Fn(VarRef) -> Lit) -> ValidityScope {
+        ValidityScope::new([presence(self.lhs), presence(self.rhs)], [])
+    }
 }
diff --git a/model/src/lang/reification.rs b/model/src/lang/reification.rs
index 33ad9ec9..f8653273 100644
--- a/model/src/lang/reification.rs
+++ b/model/src/lang/reification.rs
@@ -1,6 +1,8 @@
-use crate::bounds::Lit;
 use crate::lang::expr::Normalize;
 use crate::lang::normal_form::NormalExpr;
+use crate::lang::ValidityScope;
+use aries_core::literals::Disjunction;
+use aries_core::*;
 use std::any::{Any, TypeId};
 use std::collections::HashMap;
 use std::fmt::Debug;
@@ -16,8 +18,23 @@ pub trait ReifiableExpr: ExprInterface + PartialEq + Eq + Hash {}
 impl<T: ExprInterface + PartialEq + Eq + Hash> ReifiableExpr for T {}
 
 /// Capabilities of an expression that is read from the model.
-pub trait ExprInterface: Any + Debug + Send + Sync + 'static {}
-impl<T: Any + Debug + Send + Sync + 'static> ExprInterface for T {}
+pub trait ExprInterface: Any + Debug + Send + Sync + 'static {
+    /// Returns a description of the scope in which the expression is valid (i.e. can be evaluated to a value).
+    fn validity_scope(&self, presence: &dyn Fn(VarRef) -> Lit) -> ValidityScope;
+}
+
+impl ExprInterface for Disjunction {
+    fn validity_scope(&self, presence: &dyn Fn(VarRef) -> Lit) -> ValidityScope {
+        ValidityScope::new(
+            self.literals().iter().map(|l| presence(l.variable())),
+            // guard by non optional literal (if one of them is true, the disjunction is defined)
+            self.literals()
+                .iter()
+                .copied()
+                .filter(|l| presence(l.variable()) == Lit::TRUE),
+        )
+    }
+}
 
 /// Dynamically typed interface to an expression in the model.
 ///
@@ -63,15 +80,15 @@ pub struct Reification {
 }
 
 impl Reification {
-    fn add(&mut self, be: WrappedExpr, lit: Lit) {
+    fn intern_raw_expr_as(&mut self, be: WrappedExpr, lit: Lit) {
         assert!(!self.map.contains_key(&be));
         self.binding_events.push((lit, BindTarget::Expr(be.value.clone())));
         self.map.insert(be, lit);
     }
 
     /// If this expression was previously interned, returns the literal it was bound to.
-    pub fn interned<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, expr: Expr) -> Option<Lit> {
-        let nf = expr.normalize();
+    /// TODO: we should accept a borrowed expression.
+    pub fn interned<T: ReifiableExpr>(&mut self, nf: NormalExpr<T>) -> Option<Lit> {
         match nf {
             NormalExpr::Literal(l) => Some(l),
             NormalExpr::Pos(e) => self.map.get(&WrappedExpr::new(e)).copied(),
@@ -80,12 +97,11 @@ impl Reification {
     }
 
     /// Interns the user-facing expression, creating a new literal with the `make_lit` closure if the the expression had no preexisting binding.
-    pub fn intern<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, expr: Expr, make_lit: impl FnOnce() -> Lit) -> Lit {
-        let nf = expr.normalize();
+    pub fn intern_as<T: ReifiableExpr>(&mut self, nf: NormalExpr<T>, lit: Lit) {
         match nf {
-            NormalExpr::Literal(l) => l,
-            NormalExpr::Pos(e) => self.intern_raw_expr(WrappedExpr::new(e), make_lit),
-            NormalExpr::Neg(e) => !self.intern_raw_expr(WrappedExpr::new(e), make_lit),
+            NormalExpr::Literal(_) => panic!("Cannot be interned"),
+            NormalExpr::Pos(e) => self.intern_raw_expr_as(WrappedExpr::new(e), lit),
+            NormalExpr::Neg(e) => self.intern_raw_expr_as(WrappedExpr::new(e), !lit),
         }
     }
 
@@ -99,7 +115,7 @@ impl Reification {
                 if let Some(prev) = self.map.get(&e).copied() {
                     self.bind_literals(prev, literal)
                 } else {
-                    self.add(e, literal);
+                    self.intern_raw_expr_as(e, literal);
                 }
             }
             NormalExpr::Neg(e) => {
@@ -107,7 +123,7 @@ impl Reification {
                 if let Some(prev) = self.map.get(&e).copied() {
                     self.bind_literals(prev, !literal)
                 } else {
-                    self.add(e, !literal);
+                    self.intern_raw_expr_as(e, !literal);
                 }
             }
         }
@@ -118,16 +134,6 @@ impl Reification {
         self.binding_events.push((l1, BindTarget::Literal(l2)));
     }
 
-    fn intern_raw_expr(&mut self, be: WrappedExpr, make_lit: impl FnOnce() -> Lit) -> Lit {
-        if let Some(l) = self.map.get(&be) {
-            *l
-        } else {
-            let l = make_lit();
-            self.add(be, l);
-            l
-        }
-    }
-
     /// Returns the event that this cursors points and advances it to the next one.
     pub fn pop_next_event(&self, cursor: &mut BindingCursor) -> Option<&(Lit, BindTarget)> {
         let ret = self.binding_events.get(cursor.0);
@@ -206,10 +212,9 @@ impl BindingCursor {
 
 #[cfg(test)]
 mod tests {
-    use crate::bounds::Lit;
-    use crate::lang::expr::*;
-    use crate::lang::reification::Reification;
-    use crate::lang::{IVar, VarRef};
+    use super::*;
+    use crate::lang::expr::{eq, geq, leq};
+    use crate::lang::IVar;
 
     const A: IVar = IVar::new(VarRef::from_u32(1));
     const B: IVar = IVar::new(VarRef::from_u32(2));
@@ -228,15 +233,15 @@ mod tests {
         reif.bind(l2, f);
         reif.bind(e1, f);
 
-        assert_eq!(reif.interned(l1), Some(t));
-        assert_eq!(reif.interned(l2), Some(f));
-        assert_eq!(reif.interned(e1), Some(f));
+        assert_eq!(reif.interned(l1.normalize()), Some(t));
+        assert_eq!(reif.interned(l2.normalize()), Some(f));
+        assert_eq!(reif.interned(e1.normalize()), Some(f));
 
         // same as l1
         let l1_prime = geq(B + 3, A);
-        assert_eq!(reif.interned(l1_prime), Some(t));
+        assert_eq!(reif.interned(l1_prime.normalize()), Some(t));
 
         // inverse of l1, should return the opposite literal
-        assert_eq!(reif.interned(!l1), Some(f));
+        assert_eq!(reif.interned((!l1).normalize()), Some(f));
     }
 }
diff --git a/model/src/lang/sym.rs b/model/src/lang/sym.rs
index 093e961e..8a1794a9 100644
--- a/model/src/lang/sym.rs
+++ b/model/src/lang/sym.rs
@@ -1,6 +1,7 @@
-use crate::lang::{ConversionError, IAtom, IVar, IntCst, VarRef};
+use crate::lang::{ConversionError, IAtom, IVar};
 use crate::symbols::{SymId, TypedSym};
 use crate::types::TypeId;
+use aries_core::*;
 use std::convert::TryFrom;
 
 #[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash, Debug)]
diff --git a/model/src/lang/validity_scope.rs b/model/src/lang/validity_scope.rs
new file mode 100644
index 00000000..31dc2c32
--- /dev/null
+++ b/model/src/lang/validity_scope.rs
@@ -0,0 +1,88 @@
+use aries_core::literals::{LitSet, StableLitSet};
+use aries_core::Lit;
+
+/// Represents the scope in which a given expression is valid.
+///
+/// It is composed of:
+/// - `required_presence`: the set of presence literals of all variables that appear in the expression
+/// - `guards`: a set of literals such that if one of them is true, the expression is defined.  
+///
+/// # Example (no guards)
+///
+/// Consider the expression `a <= b`, where `a` and `b` are optional variable.
+/// The value of the expression can only be computed when both `a` and `b` are present,
+/// that is when `presence(a) & presence(b)` holds.
+///
+/// For this example, the scope should be:
+/// `ValidityScope = { required_presence: [presence(a), presence(b)], guards = [] }`
+///
+/// # Example (with guards)
+///
+/// In some cases however, it might be the case that an expression is valid even when some of
+/// its parts are not.
+/// For instance, consider the expression `!presence(a) | (a <= 2)`.
+/// If `!presence(a)` holds, then the expression evaluates to true regardless of the value of
+/// of `(a <= 2)` and in particular when it is undefined due the absence of `a`.
+///
+/// For this example, the scope should be:
+/// `ValidityScope = { required_presence: [presence(a)], guards = [!presence(a)] }`
+///
+/// # Flattening to a conjunction of literals
+///
+/// The [`ValidityScope::flatten`] method allows transforming a `ValidityScope` into a conjunction
+/// of literals that must hold for the expression to be defined.
+#[derive(Debug)]
+pub struct ValidityScope {
+    required_presence: LitSet,
+    guards: Vec<Lit>,
+}
+
+impl ValidityScope {
+    /// An empty scope, corresponding to an always valid expression.
+    pub const EMPTY: ValidityScope = Self {
+        required_presence: LitSet::EMPTY,
+        guards: vec![],
+    };
+
+    pub fn new(required: impl IntoIterator<Item = Lit>, guards: impl IntoIterator<Item = Lit>) -> Self {
+        Self {
+            required_presence: required.into(),
+            guards: guards.into_iter().filter(|&l| l != Lit::FALSE).collect(),
+        }
+    }
+
+    /// Flatten the scope into a conjunction of literals.
+    ///
+    /// This involves:
+    ///  - replacing all presence literal defined as a conjunction of other presence literal by them
+    ///  - removing from the resulting set any literal that is guarded
+    ///  - removing any tautological (i.e. always true) literal from the set
+    ///
+    /// # Parameters
+    ///
+    /// - `flattened`: if the given literal `l` is defined as a conjunction of literals `[l1, l2, ... ln]`,
+    ///   it must return a iterator over them. Returns `None` otherwise.
+    /// - `tautology`: Returns true if the given literal always holds.   
+    pub fn to_conjunction<Lits: IntoIterator<Item = Lit>>(
+        &self,
+        flattened: impl Fn(Lit) -> Option<Lits>,
+        tautology: impl Fn(Lit) -> bool,
+    ) -> StableLitSet {
+        let mut set = LitSet::empty();
+        for l in self.required_presence.literals() {
+            if let Some(flat) = flattened(l) {
+                for l in flat {
+                    if !tautology(l) {
+                        set.insert(l);
+                    }
+                }
+            } else if !tautology(l) {
+                set.insert(l)
+            }
+        }
+        for &guard in &self.guards {
+            set.remove(!guard, &tautology)
+        }
+        set.into_sorted()
+    }
+}
diff --git a/model/src/lang/variables.rs b/model/src/lang/variables.rs
index 8c55e6a2..0b9f5b72 100644
--- a/model/src/lang/variables.rs
+++ b/model/src/lang/variables.rs
@@ -1,5 +1,6 @@
 use crate::lang::variables::Variable::*;
-use crate::lang::{BVar, ConversionError, FVar, IVar, Kind, SVar, VarRef};
+use crate::lang::{BVar, ConversionError, FVar, IVar, Kind, SVar};
+use aries_core::*;
 use std::convert::TryFrom;
 
 /// Contains a variable of any type
diff --git a/model/src/lib.rs b/model/src/lib.rs
index abe1e123..b8da0b86 100644
--- a/model/src/lib.rs
+++ b/model/src/lib.rs
@@ -1,12 +1,10 @@
-pub mod bounds;
+pub use label::Label;
+pub use model::*;
+
 pub mod decomposition;
 pub mod extensions;
 mod label;
 pub mod lang;
 mod model;
-pub mod state;
 pub mod symbols;
 pub mod types;
-
-pub use label::Label;
-pub use model::*;
diff --git a/model/src/model.rs b/model/src/model.rs
index 1daa445f..77224924 100644
--- a/model/src/model.rs
+++ b/model/src/model.rs
@@ -1,15 +1,23 @@
-use crate::bounds::Lit;
+use std::convert::TryFrom;
+use std::sync::Arc;
+
+use aries_backtrack::{Backtrack, DecLvl};
+use aries_collections::ref_store::RefMap;
+use aries_core::WriterId;
+
 use crate::extensions::{AssignmentExt, SavedAssignment, Shaped};
 use crate::label::{Label, VariableLabels};
-use crate::lang::expr::Normalize;
+use crate::lang::expr::{or, Normalize};
 use crate::lang::reification::{ReifiableExpr, Reification};
 use crate::lang::*;
-use crate::state::*;
+use crate::model::scopes::Scopes;
 use crate::symbols::SymbolTable;
 use crate::types::TypeId;
-use aries_backtrack::{Backtrack, DecLvl};
-use aries_collections::ref_store::RefMap;
-use std::sync::Arc;
+use aries_core::literals::StableLitSet;
+use aries_core::state::*;
+use aries_core::*;
+
+mod scopes;
 
 /// Defines the structure of a model: variables names, types, relations, ...
 #[derive(Clone)]
@@ -18,6 +26,7 @@ pub struct ModelShape<Lbl> {
     pub types: RefMap<VarRef, Type>,
     pub expressions: Reification,
     pub labels: VariableLabels<Lbl>,
+    pub conjunctive_scopes: Scopes,
     num_writers: u8,
 }
 
@@ -32,6 +41,7 @@ impl<Lbl: Label> ModelShape<Lbl> {
             types: Default::default(),
             expressions: Default::default(),
             labels: Default::default(),
+            conjunctive_scopes: Default::default(),
             num_writers: 0,
         }
     }
@@ -97,13 +107,98 @@ impl<Lbl: Label> Model<Lbl> {
     }
 
     pub fn new_presence_variable(&mut self, scope: Lit, label: impl Into<Lbl>) -> BVar {
-        let lit = self.state.new_presence_literal(scope);
+        let lit = self.state.new_var(0, 1).geq(1);
+        self.shape.conjunctive_scopes.insert(StableLitSet::from([lit]), lit);
+        self.state.add_implication(lit, scope);
         let var = lit.variable();
         self.shape.set_label(var, label);
         self.shape.set_type(var, Type::Bool);
         BVar::new(var)
     }
 
+    /// Returns a presence literal that is true iff all given presence literal are true.
+    pub fn get_conjunctive_scope(&mut self, presence_variables: &[Lit]) -> Lit {
+        assert!(presence_variables
+            .iter()
+            .all(|l| self.state.presence(l.variable()) == Lit::TRUE));
+        let empty: &[Lit] = &[];
+        let scope = ValidityScope::new(presence_variables.iter().copied(), empty.iter().copied());
+        let scope = scope.to_conjunction(
+            |l| self.shape.conjunctive_scopes.conjuncts(l),
+            |l| self.state.entails(l) && self.state.entailing_level(l) == DecLvl::ROOT,
+        );
+        self.new_conjunctive_presence_variable(scope)
+    }
+
+    /// Returns a literal whose presence is `scope` and that is always true.
+    ///
+    /// THis is functionnaly equivalent to creating a new optional boolean variable
+    /// with domain `[1,1]` with `presence=scope` but will ensure that only one such
+    /// variable is created in this scope.
+    pub fn get_tautology_of_scope(&mut self, scope: Lit) -> Lit {
+        self.shape
+            .conjunctive_scopes
+            .get_tautology_of_scope(scope)
+            .unwrap_or_else(|| {
+                let var = self.state.new_optional_var(1, 1, scope);
+                let lit = var.geq(1);
+                self.shape.set_type(var, Type::Bool);
+                self.shape.conjunctive_scopes.set_tautology_of_scope(scope, lit);
+                lit
+            })
+    }
+
+    fn new_conjunctive_presence_variable(&mut self, set: impl Into<StableLitSet>) -> Lit {
+        let set = set.into();
+        if let Some(l) = self.shape.conjunctive_scopes.get(&set) {
+            // scope already exists, return it immediately
+            return l;
+        }
+
+        // let the scope's set be composed of { v1, v2, ..., vn }
+        // we need to create a new literal `l` such that  `l <=> v1 & v2 & ... & vn`
+
+        // first, try to avoid creating a new literal, by simplifying the conjunctive set
+        let attempt: Option<Lit> = if let Ok([v1]) = <[Lit; 1]>::try_from(&set) {
+            // single literal v1, let l = v1
+            Some(v1)
+        } else if let Ok([v1, v2]) = <[Lit; 2]>::try_from(&set) {
+            // only two literals v1 and v2
+            if self.state.implies(v1, v2) {
+                // v1 => v2, let l = v1
+                Some(v1)
+            } else if self.state.implies(v2, v1) {
+                // v2 => v1, let l = v2
+                Some(v2)
+            } else if self.state.exclusive(v1, v2) {
+                // `v1 & v2` is always false, let l = FALSE
+                Some(Lit::FALSE)
+            } else {
+                None // no simplification found, proceed
+            }
+        } else {
+            None // no simplification found, proceed
+        };
+
+        let l = attempt.unwrap_or_else(|| {
+            // Simplification did not succeed.
+            // create a new literal `l` such that `l <=> (v1 & v2 & ... & vn)`, decomposed to:
+            // - `l => v1`, `l => v2`, ...
+            // - `l | !v1 | !v2 | ... | !vn`
+            let l = self.state.new_var(0, 1).geq(1);
+            self.shape.set_type(l.variable(), Type::Bool);
+            let mut clause = vec![l];
+            for v_i in set.literals() {
+                self.state.add_implication(l, v_i);
+                clause.push(!v_i);
+            }
+            self.enforce(or(clause));
+            l
+        });
+        self.shape.conjunctive_scopes.insert(set, l);
+        l
+    }
+
     fn create_bvar(&mut self, presence: Option<Lit>, label: impl Into<Lbl>) -> BVar {
         let dvar = if let Some(presence) = presence {
             self.state.new_optional_var(0, 1, presence)
@@ -209,29 +304,53 @@ impl<Lbl: Label> Model<Lbl> {
     }
 
     /// Interns the given expression and returns an equivalent literal.
+    /// The returned literal is *optional* and defined such that it is
+    /// present iff the expression is valid (typically meaning that all
+    /// variables involved in the expression are present).
+    ///
     /// If the expression was already interned, the handle to the previously inserted
     /// instance will be returned.
     pub fn reify<T: ReifiableExpr, Expr: Normalize<T>>(&mut self, expr: Expr) -> Lit {
-        // locals to be captured and updated by the closure for literal creation
-        let state = &mut self.state;
-        let mut created: Option<VarRef> = None;
-
-        // intern the expression, creating a new literal if necessary
-        let lit = self.shape.expressions.intern(expr, || {
-            let var = state.new_var(0, 1);
-            // notify caller that a variable was created
-            created = Some(var);
-            BVar::new(var).true_lit()
-        });
-        if let Some(v) = created {
-            // variable was created, give it a type and label
-            self.shape.set_type(v, Type::Bool);
+        let e1 = expr.normalize();
+        let e2 = expr.normalize(); // TODO: avoid this duplicated work (requires update to interned)
+        if let Some(l) = self.shape.expressions.interned(e1) {
+            l
+        } else {
+            let scope = e2.validity_scope(&|var| self.state.presence(var));
+            let scope = scope.to_conjunction(
+                |l| self.shape.conjunctive_scopes.conjuncts(l),
+                |l| self.state.entails(l),
+            );
+            let scope = self.new_conjunctive_presence_variable(scope);
+            let var = self.state.new_optional_var(0, 1, scope);
+            let lit = var.geq(1);
+            self.shape.set_type(var, Type::Bool);
+            self.shape.expressions.intern_as(e2, lit);
+            lit
         }
-        lit
     }
 
-    pub fn enforce<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, b: Expr) {
-        self.shape.expressions.bind(b, Lit::TRUE)
+    /// Enforce the given expression to be true. Similar to posting a constraint in CP sovlers.
+    ///
+    /// Internally, the expression is reified to an optional literal that is true, when the expression
+    /// is valid and absent otherwise.
+    pub fn enforce<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, expr: Expr) {
+        debug_assert_eq!(self.state.current_decision_level(), DecLvl::ROOT);
+        let normalized = expr.normalize();
+
+        // compute the scope in which the expression is valid
+        let scope = normalized.validity_scope(&|var| self.state.presence(var));
+        let scope = scope.to_conjunction(
+            |l| self.shape.conjunctive_scopes.conjuncts(l),
+            |l| self.state.entails(l),
+        );
+        let scope = self.new_conjunctive_presence_variable(scope);
+
+        // retrieve or create an optional variable that is always true in the scope
+        let lit = self.get_tautology_of_scope(scope);
+
+        // bind the expression to the scoped tautology
+        self.shape.expressions.bind(expr, lit);
     }
 
     pub fn enforce_all<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, bools: impl IntoIterator<Item = Expr>) {
@@ -269,19 +388,6 @@ impl<Lbl: Label> Model<Lbl> {
     }
 }
 
-/// Identifies an external writer to the model.
-#[derive(Ord, PartialOrd, PartialEq, Eq, Copy, Clone, Hash, Debug)]
-pub struct WriterId(pub u8);
-impl WriterId {
-    pub fn new(num: impl Into<u8>) -> WriterId {
-        WriterId(num.into())
-    }
-
-    pub fn cause(&self, cause: impl Into<u32>) -> Cause {
-        Cause::inference(*self, cause)
-    }
-}
-
 impl<Lbl: Label> Default for Model<Lbl> {
     fn default() -> Self {
         Self::new()
@@ -320,7 +426,7 @@ impl<Lbl> AssignmentExt for Model<Lbl> {
     }
 
     fn to_owned_assignment(&self) -> SavedAssignment {
-        SavedAssignment::from_model(self)
+        self.state.clone()
     }
 }
 
diff --git a/model/src/model/scopes.rs b/model/src/model/scopes.rs
new file mode 100644
index 00000000..7474cdd8
--- /dev/null
+++ b/model/src/model/scopes.rs
@@ -0,0 +1,81 @@
+use aries_core::literals::StableLitSet;
+use aries_core::*;
+use std::collections::HashMap;
+use std::sync::Arc;
+
+/// A structure to keep track of the conjunctive scopes that have been defined in the problem.
+///
+/// A conjunctive scope is created when we want to refer to a subset of the problem that exists
+/// iff all required scopes are present.
+///
+/// For instance, the expression `a <= b` is defined iff both `a` and `b` are *present*.
+/// It can be said to exist in the conjunctive scope `presence(a) & presence(b)`.  
+#[derive(Clone)]
+pub struct Scopes {
+    conjunctive_scopes: HashMap<Arc<StableLitSet>, Lit>,
+    conjunction_of: HashMap<Lit, Arc<StableLitSet>>,
+    /// Associates for each scope literal (the key) an optional literal that is true in this scope.
+    /// Invariant for any entry `(k, v)` in the map, `k = presence(v)`
+    tautologies: HashMap<Lit, Lit>,
+}
+
+impl Scopes {
+    pub fn new() -> Self {
+        let mut s = Self {
+            conjunctive_scopes: Default::default(),
+            conjunction_of: Default::default(),
+            tautologies: Default::default(),
+        };
+        s.insert(StableLitSet::EMPTY, Lit::TRUE);
+        s
+    }
+
+    /// IF defined, return the literal reprensenting the given conjunction.
+    pub fn get(&self, conjunction: &StableLitSet) -> Option<Lit> {
+        self.conjunctive_scopes.get(conjunction).copied()
+    }
+
+    /// If it already exists, return the literal that is always true in the given
+    /// scope.
+    pub fn get_tautology_of_scope(&self, scope: Lit) -> Option<Lit> {
+        self.tautologies.get(&scope).copied()
+    }
+
+    /// Record the `tautology` literal as the tautology of the given `scope`.
+    ///
+    /// It should be the case that `scope = presence(tautology)` that `tautology`
+    /// is always true when present.
+    pub fn set_tautology_of_scope(&mut self, scope: Lit, tautology: Lit) {
+        assert!(!self.tautologies.contains_key(&scope));
+        self.tautologies.insert(scope, tautology);
+    }
+
+    /// Inserts a new equivalence between a conjunctive scope and a literal.
+    ///
+    /// # Panics
+    ///
+    /// Panics if the scope was already associated to a literal.
+    pub fn insert(&mut self, conjunction: StableLitSet, literal: Lit) {
+        debug_assert!(!self.conjunctive_scopes.contains_key(&conjunction));
+
+        let conjunction = Arc::new(conjunction);
+        self.conjunction_of
+            .entry(literal)
+            .or_insert_with(|| conjunction.clone());
+        // ideally, we want to have each literal pointing to its smallest conjunctive set (relaxable)
+        debug_assert!(self.conjunction_of[&literal].len() <= conjunction.len(), "Sanity check");
+
+        self.conjunctive_scopes.insert(conjunction, literal);
+    }
+
+    /// If the literal was defined as a conjunctive scope, returns the set of literals in the conjunction.
+    pub fn conjuncts(&self, lit: Lit) -> Option<impl IntoIterator<Item = Lit> + '_> {
+        self.conjunction_of.get(&lit).map(|set| set.literals())
+    }
+}
+
+impl Default for Scopes {
+    fn default() -> Self {
+        Self::new()
+    }
+}
diff --git a/model/src/state/presence_graph.rs b/model/src/state/presence_graph.rs
deleted file mode 100644
index f0a9be99..00000000
--- a/model/src/state/presence_graph.rs
+++ /dev/null
@@ -1,51 +0,0 @@
-use crate::bounds::{Lit, Watches};
-use std::collections::VecDeque;
-
-#[derive(Clone, Default)]
-pub struct TwoSatTree {
-    edges: Watches<Lit>,
-    num_edges: usize,
-}
-
-impl TwoSatTree {
-    pub fn add_implication(&mut self, from: Lit, to: Lit) {
-        debug_assert!(
-            from.variable() > to.variable(),
-            "Invariant that should be maintained by OptDomains {:?}  {:?}",
-            from,
-            to
-        );
-        debug_assert!(!self.implies(to, from));
-        self.num_edges += 1;
-        self.edges.add_watch(to, from);
-        self.edges.add_watch(!from, !to);
-        debug_assert!(self.implies(from, to));
-        debug_assert!(self.implies(!to, !from));
-    }
-
-    pub fn implies(&self, x: Lit, y: Lit) -> bool {
-        // TODO: we can be smarter and only go up the tree.
-        //       key observation: by construction, the variable ID of a node is strictly greater than
-        //       the one of its ancestor.
-        let mut num_iter = 0;
-        let mut queue = VecDeque::new();
-        queue.push_back(x);
-        // dfs through implications
-        while let Some(curr) = queue.pop_front() {
-            for next in self.edges.watches_on(curr) {
-                if next.entails(y) {
-                    return true;
-                } else {
-                    queue.push_back(next);
-                }
-            }
-            debug_assert!(num_iter <= self.num_edges, "The structure appear not to be a tree");
-            num_iter += 1;
-        }
-        false
-    }
-
-    pub fn direct_implications_of(&self, lit: Lit) -> impl Iterator<Item = Lit> + '_ {
-        self.edges.watches_on(lit)
-    }
-}
diff --git a/planners/Cargo.toml b/planners/Cargo.toml
index 4a87ac7f..6ea2594d 100644
--- a/planners/Cargo.toml
+++ b/planners/Cargo.toml
@@ -9,6 +9,7 @@ structopt = "0.3"
 anyhow = "1.0"
 aries_collections = { path = "../collections" }
 aries_backtrack = { path = "../backtrack" }
+aries_core = { path = "../core" }
 aries_model = { path = "../model" }
 aries_planning = { path = "../planning" }
 aries_tnet = { path = "../tnet" }
diff --git a/planners/src/encode.rs b/planners/src/encode.rs
index 2441340a..a349792b 100644
--- a/planners/src/encode.rs
+++ b/planners/src/encode.rs
@@ -4,11 +4,10 @@
 use crate::encoding::{conditions, effects, refinements_of, refinements_of_task, TaskRef, HORIZON, ORIGIN};
 use crate::Model;
 use anyhow::{Context, Result};
-use aries_model::bounds::Lit;
+use aries_core::*;
 use aries_model::extensions::{AssignmentExt, Shaped};
 use aries_model::lang::expr::*;
-use aries_model::lang::Variable;
-use aries_model::lang::{FAtom, VarRef};
+use aries_model::lang::{FAtom, Variable};
 use aries_planning::chronicles::constraints::ConstraintType;
 use aries_planning::chronicles::*;
 use env_param::EnvParam;
@@ -269,11 +268,11 @@ fn enforce_refinement(t: TaskRef, supporters: Vec<TaskRef>, model: &mut Model) {
             .only_present_with(s.presence.variable(), t.presence.variable()));
         model.enforce(implies(s.presence, t.presence)); // TODO: can we get rid of this
 
-        model.enforce(opt_eq(s.start, t.start));
-        model.enforce(opt_eq(s.end, t.end));
+        model.enforce(eq(s.start, t.start));
+        model.enforce(eq(s.end, t.end));
         assert_eq!(s.task.len(), t.task.len());
         for (a, b) in s.task.iter().zip(t.task.iter()) {
-            model.enforce(opt_eq(*a, *b))
+            model.enforce(eq(*a, *b))
         }
     }
 }
@@ -417,11 +416,21 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
             supported_by_eff_conjunction.push(model.reify(f_leq(cond.end, eff_ends[eff_id])));
 
             let support_lit = model.reify(and(supported_by_eff_conjunction));
+            debug_assert!({
+                let prez_support = model.presence_literal(support_lit.variable());
+                model.state.implies(prez_cond, prez_support)
+            });
 
             // add this support expression to the support clause
             supported.push(support_lit);
         }
 
+        debug_assert!({
+            let or_reif = model.reify(or(supported.as_slice()));
+            let or_reif_prez = model.presence_literal(or_reif.variable());
+            or_reif_prez == Lit::TRUE
+        });
+
         // enforce necessary conditions for condition' support
         model.enforce(or(supported));
     }
@@ -481,14 +490,14 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
 
     for ch in &pb.chronicles {
         // chronicle finishes before the horizon and has a non negative duration
-        model.enforce(f_opt_leq(ch.chronicle.end, pb.horizon));
-        model.enforce(f_opt_leq(ch.chronicle.start, ch.chronicle.end));
+        model.enforce(f_leq(ch.chronicle.end, pb.horizon));
+        model.enforce(f_leq(ch.chronicle.start, ch.chronicle.end));
 
         // enforce temporal coherence between the chronicle and its subtasks
         for subtask in &ch.chronicle.subtasks {
-            model.enforce(f_opt_leq(subtask.start, subtask.end));
-            model.enforce(f_opt_leq(ch.chronicle.start, subtask.start));
-            model.enforce(f_opt_leq(subtask.end, ch.chronicle.end));
+            model.enforce(f_leq(subtask.start, subtask.end));
+            model.enforce(f_leq(ch.chronicle.start, subtask.start));
+            model.enforce(f_leq(subtask.end, ch.chronicle.end));
         }
     }
     add_decomposition_constraints(pb, &mut model);
diff --git a/planners/src/encoding.rs b/planners/src/encoding.rs
index 0c13330f..886c35b5 100644
--- a/planners/src/encoding.rs
+++ b/planners/src/encoding.rs
@@ -1,6 +1,6 @@
 //! Functions responsible for
 
-use aries_model::bounds::Lit;
+use aries_core::Lit;
 use aries_model::lang::FAtom;
 use aries_planning::chronicles::{ChronicleOrigin, ChronicleTemplate, Condition, Effect, FiniteProblem, Problem, Task};
 
diff --git a/planners/src/forward_search.rs b/planners/src/forward_search.rs
index be188d95..ed5d825f 100644
--- a/planners/src/forward_search.rs
+++ b/planners/src/forward_search.rs
@@ -3,9 +3,9 @@
 use crate::encoding::refinements_of;
 use crate::Model;
 use aries_backtrack::{Backtrack, DecLvl};
-use aries_model::bounds::Lit;
+use aries_core::{Lit, VarRef};
 use aries_model::extensions::{AssignmentExt, Shaped};
-use aries_model::lang::{Atom, IVar, VarRef};
+use aries_model::lang::{Atom, IVar};
 use aries_planning::chronicles::{ChronicleInstance, FiniteProblem, SubTask, VarLabel, VarType};
 use aries_solver::solver::search::{Decision, SearchControl};
 use aries_solver::solver::stats::Stats;
diff --git a/planning/Cargo.toml b/planning/Cargo.toml
index 338ab271..03ee782e 100644
--- a/planning/Cargo.toml
+++ b/planning/Cargo.toml
@@ -13,6 +13,7 @@ itertools = "0.10"
 streaming-iterator = "0.1.5"
 regex = "1"
 aries_collections = { path = "../collections" }
+aries_core = { path = "../core" }
 aries_model = { path = "../model" }
 aries_utils = { path = "../utils" }
 env_param = { path = "../env_param" }
\ No newline at end of file
diff --git a/planning/src/chronicles/concrete.rs b/planning/src/chronicles/concrete.rs
index 39752d72..9ac406bd 100644
--- a/planning/src/chronicles/concrete.rs
+++ b/planning/src/chronicles/concrete.rs
@@ -2,7 +2,7 @@ use core::fmt;
 use std::fmt::Debug;
 
 use crate::chronicles::constraints::Constraint;
-use aries_model::bounds::Lit;
+use aries_core::{Lit, VarRef};
 use aries_model::lang::*;
 
 pub type Sv = Vec<SAtom>;
diff --git a/planning/src/chronicles/constraints.rs b/planning/src/chronicles/constraints.rs
index e8767533..316e29f8 100644
--- a/planning/src/chronicles/constraints.rs
+++ b/planning/src/chronicles/constraints.rs
@@ -1,5 +1,5 @@
 use super::*;
-use aries_model::lang::{IntCst, Type};
+use aries_model::lang::Type;
 
 /// Generic representation of a constraint on a set of variables
 #[derive(Debug, Clone)]
@@ -7,7 +7,9 @@ pub struct Constraint {
     pub variables: Vec<Atom>,
     pub tpe: ConstraintType,
 }
+use aries_core::IntCst;
 use ConstraintType::*;
+
 impl Constraint {
     pub fn lt(a: impl Into<Atom>, b: impl Into<Atom>) -> Constraint {
         Constraint {
diff --git a/planning/src/chronicles/mod.rs b/planning/src/chronicles/mod.rs
index c81e1fcb..04137000 100644
--- a/planning/src/chronicles/mod.rs
+++ b/planning/src/chronicles/mod.rs
@@ -4,18 +4,17 @@ pub mod constraints;
 pub mod preprocessing;
 mod templates;
 
-use aries_model::symbols::{SymId, SymbolTable, TypedSym};
+pub use concrete::*;
 
 use self::constraints::Table;
-use aries_model::lang::{Atom, FAtom, IAtom, IntCst, Type, Variable};
+use aries_core::IntCst;
+use aries_model::extensions::Shaped;
+use aries_model::lang::{Atom, FAtom, IAtom, Type, Variable};
+use aries_model::symbols::{SymId, SymbolTable, TypedSym};
 use aries_model::Model;
-
 use std::fmt::Formatter;
 use std::sync::Arc;
 
-use aries_model::extensions::Shaped;
-pub use concrete::*;
-
 /// Time being represented as a fixed point numeral, this is the denominator of any time numeral.
 /// Having a time scale 100, will allow a resolution of `0.01` for time values.
 pub const TIME_SCALE: IntCst = 10;
@@ -194,6 +193,12 @@ impl std::fmt::Debug for VarLabel {
     }
 }
 
+impl std::fmt::Display for VarLabel {
+    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
+        write!(f, "{:?}", self)
+    }
+}
+
 #[derive(Copy, Clone, Debug, Eq, PartialEq, Hash)]
 pub enum Container {
     Base,
diff --git a/planning/src/classical/mod.rs b/planning/src/classical/mod.rs
index f3a724c4..fe6391ad 100644
--- a/planning/src/classical/mod.rs
+++ b/planning/src/classical/mod.rs
@@ -105,10 +105,7 @@ pub fn from_chronicles(chronicles: &crate::chronicles::Problem) -> Result<Lifted
     let ctx = &chronicles.context;
     for instance in &chronicles.chronicles {
         let ch = &instance.chronicle;
-        anyhow::ensure!(
-            ch.presence == aries_model::bounds::Lit::TRUE,
-            "A chronicle instance is optional",
-        );
+        anyhow::ensure!(ch.presence == aries_core::Lit::TRUE, "A chronicle instance is optional",);
         for eff in &ch.effects {
             anyhow::ensure!(
                 eff.effective_start() == eff.transition_start(),
diff --git a/planning/src/parsing/mod.rs b/planning/src/parsing/mod.rs
index 47174ae4..f1dc3194 100644
--- a/planning/src/parsing/mod.rs
+++ b/planning/src/parsing/mod.rs
@@ -8,7 +8,7 @@ use crate::parsing::pddl::{PddlFeature, TypedSymbol};
 use crate::chronicles::constraints::Constraint;
 use crate::parsing::sexpr::SExpr;
 use anyhow::{Context, Result};
-use aries_model::bounds::Lit;
+use aries_core::*;
 use aries_model::extensions::Shaped;
 use aries_model::lang::*;
 use aries_model::symbols::SymbolTable;
diff --git a/solver/Cargo.toml b/solver/Cargo.toml
index a6d57d27..576dca65 100644
--- a/solver/Cargo.toml
+++ b/solver/Cargo.toml
@@ -17,6 +17,7 @@ cpu_cycles = []
 [dependencies]
 aries_collections = { path = "../collections" }
 aries_backtrack = { path = "../backtrack" }
+aries_core = { path = "../core" }
 aries_model = { path = "../model" }
 env_param = { path = "../env_param" }
 itertools = "0.10.0"
diff --git a/solver/src/clauses/mod.rs b/solver/src/clauses/mod.rs
index b15aee22..72569f40 100644
--- a/solver/src/clauses/mod.rs
+++ b/solver/src/clauses/mod.rs
@@ -1,7 +1,8 @@
 use aries_backtrack::EventIndex;
 use aries_collections::ref_store::{RefMap, RefVec};
 use aries_collections::*;
-use aries_model::bounds::{Disjunction, Lit};
+use aries_core::literals::Disjunction;
+use aries_core::Lit;
 use std::cmp::Ordering::Equal;
 use std::fmt::{Debug, Display, Error, Formatter};
 use std::ops::{Index, IndexMut};
diff --git a/solver/src/lib.rs b/solver/src/lib.rs
index 98e08431..27b6c8f1 100644
--- a/solver/src/lib.rs
+++ b/solver/src/lib.rs
@@ -1,5 +1,12 @@
 extern crate aries_model;
 
+use aries_backtrack::{Backtrack, DecLvl};
+use aries_core::state::{Domains, Explanation, InvalidUpdate};
+use aries_core::*;
+use aries_model::lang::reification::Expr;
+
+use crate::solver::BindingResult;
+
 pub mod clauses;
 pub(crate) mod cpu_time;
 pub mod parallel_solver;
@@ -7,13 +14,6 @@ pub mod signals;
 pub mod solver;
 pub mod theories;
 
-use crate::solver::BindingResult;
-use aries_backtrack::{Backtrack, DecLvl};
-use aries_model::bounds::Lit;
-use aries_model::lang::reification::Expr;
-use aries_model::state::{Domains, Explanation, InvalidUpdate};
-use aries_model::WriterId;
-
 /// A trait that provides the ability to bind an arbitrary expression to a literal.
 pub trait Bind {
     /// When invoke, the module should add constraints to enforce `lit <=> expr`.
diff --git a/solver/src/parallel_solver.rs b/solver/src/parallel_solver.rs
index 31aceccf..b53969db 100644
--- a/solver/src/parallel_solver.rs
+++ b/solver/src/parallel_solver.rs
@@ -1,7 +1,8 @@
 use crate::signals::{InputSignal, InputStream, OutputSignal, SolverOutput, ThreadID};
 use crate::solver::{Exit, Solver};
+use aries_core::IntCst;
 use aries_model::extensions::{SavedAssignment, Shaped};
-use aries_model::lang::{IAtom, IntCst};
+use aries_model::lang::IAtom;
 use aries_model::{Label, ModelShape};
 use std::sync::mpsc::{channel, Receiver, Sender};
 use std::sync::Arc;
diff --git a/solver/src/signals.rs b/solver/src/signals.rs
index 5eeba8f3..c9ebcf5f 100644
--- a/solver/src/signals.rs
+++ b/solver/src/signals.rs
@@ -1,4 +1,4 @@
-use aries_model::bounds::Disjunction;
+use aries_core::literals::Disjunction;
 use aries_model::extensions::SavedAssignment;
 use env_param::EnvParam;
 use std::sync::atomic::{AtomicUsize, Ordering};
diff --git a/solver/src/solver.rs b/solver/src/solver.rs
index 535531d6..19459c9c 100644
--- a/solver/src/solver.rs
+++ b/solver/src/solver.rs
@@ -1,7 +1,18 @@
-pub mod sat_solver;
-pub mod search;
-pub mod stats;
-pub mod theory_solver;
+use std::fmt::Formatter;
+use std::sync::mpsc::Sender;
+use std::sync::Arc;
+use std::time::Instant;
+
+use aries_backtrack::{Backtrack, DecLvl};
+use aries_core::literals::Disjunction;
+use aries_core::state::*;
+use aries_core::*;
+use aries_model::decomposition::Constraints;
+use aries_model::extensions::{AssignmentExt, DisjunctionExt, SavedAssignment, Shaped};
+use aries_model::lang::expr::Normalize;
+use aries_model::lang::reification::{BindTarget, ReifiableExpr};
+use aries_model::lang::IAtom;
+use aries_model::{Label, Model, ModelShape};
 
 use crate::cpu_time::CycleCount;
 use crate::cpu_time::StartCycleCount;
@@ -11,19 +22,11 @@ use crate::solver::search::{default_brancher, Decision, SearchControl};
 use crate::solver::stats::Stats;
 use crate::solver::theory_solver::TheorySolver;
 use crate::{Bind, Contradiction, Theory};
-use aries_backtrack::{Backtrack, DecLvl};
-use aries_model::bounds::{Disjunction, Lit};
-use aries_model::decomposition::Constraints;
-use aries_model::extensions::{AssignmentExt, DisjunctionExt, SavedAssignment, Shaped};
-use aries_model::lang::expr::Normalize;
-use aries_model::lang::reification::{BindTarget, ReifiableExpr};
-use aries_model::lang::{IAtom, IntCst};
-use aries_model::state::{Cause, Domains, Explainer, Explanation, InferenceCause};
-use aries_model::{Label, Model, ModelShape, WriterId};
-use std::fmt::Formatter;
-use std::sync::mpsc::Sender;
-use std::sync::Arc;
-use std::time::Instant;
+
+pub mod sat_solver;
+pub mod search;
+pub mod stats;
+pub mod theory_solver;
 
 /// Result of the `_solve` method.
 enum SolveResult {
@@ -254,7 +257,7 @@ impl<Lbl: Label> Solver<Lbl> {
             }
             match self.brancher.next_decision(&self.stats, &self.model) {
                 Some(Decision::SetLiteral(lit)) => {
-                    // println!("Decision on: {} -- {:?}", self.model.discrete.fmt(lit.variable()), lit);
+                    // println!("Decision: {}", self.model.fmt(lit));
                     self.decide(lit);
                 }
                 Some(Decision::Restart) => {
@@ -332,33 +335,49 @@ impl<Lbl: Label> Solver<Lbl> {
         self.stats.num_decisions += 1;
     }
 
-    /// Determines the appropriate backtrack level for this clause.
-    /// Ideally this should be the earliest level at which the clause is unit.
+    /// Determines the appropriate backtrack level for this clause and returns the literal that
+    /// is asserted at this level.
+    ///
+    /// The common understanding is that it should be the earliest level at which the clause is unit.
+    /// However, the explanation of eager propagation of optional might generate explanations where some
+    /// literals are not violated, those are ignored in determining the asserting level.
     ///
-    /// In the general case, there might not be such level. This means that the two literals
-    /// that became violated the latest, are violated at the same decision level.
-    /// In this case, we select the latest decision level in which the clause is not violated
-    fn backtrack_level_for_clause(&self, clause: &[Lit]) -> Option<DecLvl> {
-        debug_assert_eq!(self.model.state.value_of_clause(clause.iter().copied()), Some(false));
+    /// If there is more that one violated literal at the latest level, then no literal is asserted
+    /// and the bactrack level is set to the ante-last level (might occur with clause sharing).
+    fn backtrack_level_for_clause(&self, clause: &[Lit]) -> Option<(DecLvl, Option<Lit>)> {
+        // debug_assert_eq!(self.model.state.value_of_clause(clause.iter().copied()), Some(false));
+
+        // level of the the two latest set element of the clause
         let mut max = DecLvl::ROOT;
         let mut max_next = DecLvl::ROOT;
+
+        // the latest violated literal, which will be the asserted literal.
+        let mut asserted = None;
+
         for &lit in clause {
-            if let Some(ev) = self.model.state.implying_event(!lit) {
-                let dl = self.model.state.trail().decision_level(ev);
-                if dl > max {
-                    max_next = max;
-                    max = dl;
-                } else if dl > max_next {
-                    max_next = dl;
+            // only consider literals that are violated.
+            // non violated literals might be there because of eager propagation of optionals.
+            if self.model.state.entails(!lit) {
+                if let Some(ev) = self.model.state.implying_event(!lit) {
+                    let dl = self.model.state.trail().decision_level(ev);
+                    if dl > max {
+                        max_next = max;
+                        max = dl;
+                        asserted = Some(lit);
+                    } else if dl > max_next {
+                        max_next = dl;
+                    }
                 }
             }
         }
+
         if max == DecLvl::ROOT {
             None
         } else if max == max_next {
-            Some(max - 1)
+            Some((max - 1, None))
         } else {
-            Some(max_next)
+            assert!(max_next < max);
+            Some((max_next, asserted))
         }
     }
 
@@ -368,8 +387,31 @@ impl<Lbl: Label> Solver<Lbl> {
     /// Returns `false` if the clause is conflicting at the root and thus constitutes a contradiction.
     #[must_use]
     fn add_conflicting_clause_and_backtrack(&mut self, expl: Disjunction) -> bool {
-        // println!("conflict: {:?}", &expl);
-        if let Some(dl) = self.backtrack_level_for_clause(expl.literals()) {
+        // // print the clause before analysis
+        // println!("conflict ({}) :", expl.literals().len());
+        // for &l in expl.literals() {
+        //     if !self.model.state.entails(!l) {
+        //         print!("  > {}", self.model.fmt(l));
+        //     } else {
+        //         print!("    {}", self.model.fmt(l));
+        //     }
+        //     let prez = self.model.state.presence(l.variable());
+        //     let scope: Vec<Lit> = self
+        //         .model
+        //         .state
+        //         .implications
+        //         .direct_implications_of(l.variable().geq(1))
+        //         .collect();
+        //     print!("  / {}   <<<<   ", self.model.fmt(prez));
+        //     print!("  [");
+        //     for prez in scope {
+        //         print!("  & {}", self.model.fmt(prez));
+        //     }
+        //     println!("]");
+        // }
+        // println!();
+
+        if let Some((dl, asserted)) = self.backtrack_level_for_clause(expl.literals()) {
             // backtrack
             self.restore(dl);
             debug_assert_eq!(self.model.state.value_of_clause(&expl), None);
@@ -383,8 +425,13 @@ impl<Lbl: Label> Solver<Lbl> {
                 self.brancher.bump_activity(b.variable(), &self.model);
             }
 
-            // add clause to sat solver
-            self.reasoners.sat.add_forgettable_clause(expl);
+            if let Some(asserted) = asserted {
+                // add clause to sat solver, making sure the asserted literal is set to true
+                self.reasoners.sat.add_learnt_clause(expl, asserted);
+            } else {
+                // no asserted literal, just add a forgettable clause
+                self.reasoners.sat.add_forgettable_clause(expl)
+            }
 
             true
         } else {
diff --git a/solver/src/solver/sat_solver.rs b/solver/src/solver/sat_solver.rs
index d5d6358a..ec2efa3c 100644
--- a/solver/src/solver/sat_solver.rs
+++ b/solver/src/solver/sat_solver.rs
@@ -1,15 +1,17 @@
-use crate::clauses::{Clause, ClauseDb, ClauseId, ClausesParams};
-use crate::solver::BindingResult;
-use crate::BindSplit;
+use itertools::Itertools;
+use smallvec::alloc::collections::VecDeque;
+
 use aries_backtrack::{Backtrack, DecLvl, ObsTrailCursor, Trail};
 use aries_collections::set::RefSet;
-use aries_model::bounds::{Disjunction, Lit, WatchSet, Watches};
+use aries_core::literals::{Disjunction, WatchSet, Watches};
+use aries_core::state::{Domains, Event, Explanation};
+use aries_core::*;
 use aries_model::extensions::DisjunctionExt;
 use aries_model::lang::reification::{downcast, Expr};
-use aries_model::state::{Domains, Event, Explanation};
-use aries_model::WriterId;
-use itertools::Itertools;
-use smallvec::alloc::collections::VecDeque;
+
+use crate::clauses::{Clause, ClauseDb, ClauseId, ClausesParams};
+use crate::solver::BindingResult;
+use crate::BindSplit;
 
 #[derive(Clone)]
 struct ClauseLocks {
@@ -106,6 +108,17 @@ impl Default for Stats {
     }
 }
 
+/// A clause that has been recorded but not propagated yet.
+#[derive(Copy, Clone)]
+struct PendingClause {
+    /// Id of the clause to propagate
+    clause: ClauseId,
+    /// If non-empty, the literal is entailed by the clause at the current level.
+    /// This literal MUST be set to True, with this clause as its cause even if the clause is not unit.
+    /// This situation might happen when the clause was learnt from a conflict involving the eager propagation of optional variables.
+    asserted_literal: Option<Lit>,
+}
+
 #[derive(Clone)]
 pub struct SatSolver {
     clauses: ClauseDb,
@@ -113,7 +126,7 @@ pub struct SatSolver {
     events_stream: ObsTrailCursor<Event>,
     token: WriterId,
     /// Clauses that have been added to the database but not processed and propagated yet
-    pending_clauses: VecDeque<ClauseId>,
+    pending_clauses: VecDeque<PendingClause>,
     /// Clauses that are locked (can't be remove from the database).
     /// A clause is locked if it asserted a literal and thus might be needed for an explanation
     locks: ClauseLocks,
@@ -154,9 +167,25 @@ impl SatSolver {
         self.add_clause_impl(clause.into(), true);
     }
 
+    /// Adds an asserting clause that was learnt.
+    /// On the next propagation, the clause will be propagated and the `asserted` literal set
+    /// to true (even is the clause is not unit).
+    pub fn add_learnt_clause(&mut self, clause: impl Into<Disjunction>, asserted: Lit) {
+        let clause = clause.into();
+        debug_assert!(clause.contains(asserted));
+        let cl_id = self.clauses.add_clause(Clause::new(clause), true);
+        self.pending_clauses.push_back(PendingClause {
+            clause: cl_id,
+            asserted_literal: Some(asserted),
+        });
+    }
+
     fn add_clause_impl(&mut self, clause: Disjunction, learnt: bool) -> ClauseId {
         let cl_id = self.clauses.add_clause(Clause::new(clause), learnt);
-        self.pending_clauses.push_back(cl_id);
+        self.pending_clauses.push_back(PendingClause {
+            clause: cl_id,
+            asserted_literal: None,
+        });
         cl_id
     }
 
@@ -270,11 +299,21 @@ impl SatSolver {
 
     fn propagate_impl(&mut self, model: &mut Domains) -> Result<(), ClauseId> {
         // process all clauses that have been added since last propagation
-        while let Some(cl) = self.pending_clauses.pop_front() {
-            if let Some(conflict) = self.process_arbitrary_clause(cl, model) {
+        while let Some(PendingClause {
+            clause,
+            asserted_literal,
+        }) = self.pending_clauses.pop_front()
+        {
+            if let Some(conflict) = self.process_arbitrary_clause(clause, model) {
                 self.stats.conflicts += 1;
                 return Err(conflict);
             }
+            if let Some(asserted) = asserted_literal {
+                if !model.entails(asserted) {
+                    assert!(!model.entails(!asserted));
+                    self.set_from_unit_propagation(asserted, clause, model);
+                }
+            }
         }
         // grow or shrink database. Placed here to be as close as possible to initial minisat
         // implementation where this appeared in search
@@ -398,9 +437,14 @@ impl SatSolver {
     }
 
     fn set_from_unit_propagation(&mut self, literal: Lit, propagating_clause: ClauseId, model: &mut Domains) {
-        // lock clause to ensure it will not be removed. This is necessary as we might need it to provide an explanation
-        self.lock(propagating_clause);
-        model.set_unchecked(literal, self.token.cause(propagating_clause));
+        // Set the literal to false.
+        // We know that no inconsistency will occur (from the invariants of unit propagation.
+        // However, it might be the case that nothing happens if the literal is already known to be absent.
+        let changed_something = model.set(literal, self.token.cause(propagating_clause)).unwrap();
+        if changed_something {
+            // lock clause to ensure it will not be removed. This is necessary as we might need it to provide an explanation
+            self.lock(propagating_clause);
+        }
     }
 
     fn lock(&mut self, clause: ClauseId) {
@@ -445,13 +489,14 @@ impl SatSolver {
         // bump the activity of any clause use in an explanation
         self.clauses.bump_activity(clause);
         let clause = &self.clauses[clause];
-        debug_assert!(model.unit_clause(clause));
+        // in a normal sat solver, we would expect the clause to be unit,
+        // however it is not necessarily the case with eager propagation of optionals
+        // debug_assert!(model.unit_clause(clause));
         explanation.reserve(clause.len() - 1);
         for l in clause {
             if l.entails(literal) {
                 debug_assert_eq!(model.value(l), None)
             } else {
-                debug_assert_eq!(model.value(l), Some(false));
                 explanation.push(!l);
             }
         }
@@ -566,13 +611,12 @@ impl Backtrack for SatSolver {
 
 #[cfg(test)]
 mod tests {
-    use crate::solver::sat_solver::SatSolver;
+    use super::*;
     use aries_backtrack::Backtrack;
-    use aries_model::bounds::Lit;
+    use aries_core::state::Cause;
     use aries_model::extensions::AssignmentExt;
-    use aries_model::lang::IntCst;
-    use aries_model::state::Cause;
-    use aries_model::WriterId;
+
+    use crate::solver::sat_solver::SatSolver;
 
     type Model = aries_model::Model<&'static str>;
 
diff --git a/solver/src/solver/search.rs b/solver/src/solver/search.rs
index 84fbf307..962d6b07 100644
--- a/solver/src/solver/search.rs
+++ b/solver/src/solver/search.rs
@@ -2,9 +2,8 @@ pub mod activity;
 
 use crate::solver::stats::Stats;
 use aries_backtrack::Backtrack;
-use aries_model::bounds::Lit;
+use aries_core::*;
 use aries_model::extensions::SavedAssignment;
-use aries_model::lang::{IntCst, VarRef};
 use aries_model::{Label, Model};
 
 pub enum Decision {
diff --git a/solver/src/solver/search/activity.rs b/solver/src/solver/search/activity.rs
index c735c3a1..23b106c8 100644
--- a/solver/src/solver/search/activity.rs
+++ b/solver/src/solver/search/activity.rs
@@ -3,10 +3,10 @@ use crate::solver::stats::Stats;
 use aries_backtrack::{Backtrack, DecLvl, ObsTrailCursor, Trail};
 use aries_collections::heap::IdxHeap;
 use aries_collections::ref_store::RefMap;
-use aries_model::bounds::{Lit, Watches};
+use aries_core::literals::Watches;
+use aries_core::state::{Event, IntDomain};
+use aries_core::*;
 use aries_model::extensions::{AssignmentExt, SavedAssignment, Shaped};
-use aries_model::lang::{IntCst, VarRef};
-use aries_model::state::{Event, IntDomain};
 use aries_model::{Label, Model};
 use env_param::EnvParam;
 use itertools::Itertools;
diff --git a/solver/src/solver/theory_solver.rs b/solver/src/solver/theory_solver.rs
index 5ba21b25..e9cb32cc 100644
--- a/solver/src/solver/theory_solver.rs
+++ b/solver/src/solver/theory_solver.rs
@@ -1,9 +1,9 @@
 use crate::solver::BindingResult;
 use crate::{Contradiction, Theory};
 use aries_backtrack::{Backtrack, DecLvl};
-use aries_model::bounds::Lit;
+use aries_core::state::Domains;
+use aries_core::Lit;
 use aries_model::lang::reification::Expr;
-use aries_model::state::Domains;
 
 // TODO: remove this useless layer
 pub struct TheorySolver {
diff --git a/tnet/Cargo.toml b/tnet/Cargo.toml
index 76197046..13152493 100644
--- a/tnet/Cargo.toml
+++ b/tnet/Cargo.toml
@@ -6,10 +6,12 @@ edition = "2018"
 
 [dependencies]
 num-traits = "0.2"
+aries_core = { path = "../core" }
+aries_model = { path = "../model" }
 aries_solver = { path = "../solver" }
 aries_collections = { path = "../collections" }
 aries_backtrack = { path = "../backtrack" }
-aries_model = { path = "../model" }
+
 env_param = { path = "../env_param" }
 
 [dev-dependencies]
diff --git a/tnet/benches/propagation.rs b/tnet/benches/propagation.rs
index cdcc0fe0..dd7e8a58 100644
--- a/tnet/benches/propagation.rs
+++ b/tnet/benches/propagation.rs
@@ -1,4 +1,4 @@
-use aries_model::bounds::{Lit, Relation};
+use aries_core::literals::{Lit, Relation};
 use aries_model::lang::IntCst;
 use aries_tnet::theory::Stn;
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
@@ -107,7 +107,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
             b.iter(|| propagate_bounds(black_box(stn.clone()), black_box(&bounds)))
         });
 
-        // lower bounds only
+        // lower literals only
         let bounds_subset: Vec<_> = bounds
             .iter()
             .copied()
@@ -120,7 +120,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
             });
         }
 
-        // upper bounds only
+        // upper literals only
         let bounds_subset: Vec<_> = bounds
             .iter()
             .copied()
diff --git a/tnet/src/lib.rs b/tnet/src/lib.rs
index dbca7f55..4b18ce98 100644
--- a/tnet/src/lib.rs
+++ b/tnet/src/lib.rs
@@ -1,33 +1,2 @@
-use crate::theory::*;
-
-pub(crate) mod distances;
 pub mod stn;
 pub mod theory;
-
-/// Creates a new edge representing a maximum delay from one timepoint to another.
-///  - constraint: `to - from <= max_delay`
-///  - edge: from ---- (max_delay) ---> to
-pub fn max_delay(from: Timepoint, to: Timepoint, max_delay: W) -> Edge {
-    Edge::new(from, to, max_delay)
-}
-
-/// Creates a new edge representing a minimal delay from one timepoint to another
-///  - constraint `to - from >= min_delay`
-///  - edge: from <--- (-min_delay) ---- to
-pub fn min_delay(from: Timepoint, to: Timepoint, min_delay: W) -> Edge {
-    Edge::new(to, from, -min_delay)
-}
-
-/// Creates a new edge specifying that the first timepoint must be before the second
-///  - constraint: `first <= second`
-///  - edge: first <--- 0 ---- second
-pub fn before_eq(first: Timepoint, second: Timepoint) -> Edge {
-    min_delay(first, second, 0)
-}
-
-/// Creates a new edge specifying that the first timepoint must be strictly before the second
-///  - constraint: `first < second`
-///  - edge: first <--- (epsilon/step) ---- second
-pub fn strictly_before(first: Timepoint, second: Timepoint) -> Edge {
-    min_delay(first, second, 1)
-}
diff --git a/tnet/src/stn.rs b/tnet/src/stn.rs
index 34840060..34ee2a1a 100644
--- a/tnet/src/stn.rs
+++ b/tnet/src/stn.rs
@@ -1,7 +1,8 @@
-use crate::theory::{can_propagate, edge_presence, EdgeId, StnConfig, StnTheory, Timepoint, W};
+use crate::theory::{EdgeId, StnConfig, StnTheory, Timepoint, W};
 use aries_backtrack::Backtrack;
-use aries_model::bounds::{Disjunction, Lit};
-use aries_model::state::{Cause, Domains, Explainer, Explanation, InferenceCause};
+use aries_core::literals::Disjunction;
+use aries_core::state::{Cause, Domains, Explainer, Explanation, InferenceCause};
+use aries_core::Lit;
 use aries_model::Model;
 use aries_solver::{Contradiction, Theory};
 
@@ -16,7 +17,7 @@ impl Stn {
         let stn = StnTheory::new(model.new_write_token(), StnConfig::default());
         Stn { stn, model }
     }
-    pub fn with_config(config: StnConfig) -> Self {
+    pub fn new_with_config(config: StnConfig) -> Self {
         let mut model = Model::new();
         let stn = StnTheory::new(model.new_write_token(), config);
         Stn { stn, model }
@@ -35,51 +36,35 @@ impl Stn {
     }
 
     pub fn add_edge(&mut self, source: Timepoint, target: Timepoint, weight: W) -> EdgeId {
+        let valid_edge = self.get_conjunctive_scope(source, target);
+        let active_edge = self.model.get_tautology_of_scope(valid_edge);
+        debug_assert!(self.model.state.entails(active_edge));
         self.stn
-            .add_reified_edge(Lit::TRUE, source, target, weight, &self.model.state)
-    }
-
-    pub fn add_reified_edge(&mut self, literal: Lit, source: Timepoint, target: Timepoint, weight: W) -> EdgeId {
-        self.stn
-            .add_reified_edge(literal, source, target, weight, &self.model.state)
-    }
-
-    pub fn add_optional_true_edge(
-        &mut self,
-        source: impl Into<Timepoint>,
-        target: impl Into<Timepoint>,
-        weight: W,
-        forward_prop: Lit,
-        backward_prop: Lit,
-        presence: Option<Lit>,
-    ) -> EdgeId {
-        self.stn.add_optional_true_edge(
-            source,
-            target,
-            weight,
-            forward_prop,
-            backward_prop,
-            presence,
-            &self.model.state,
-        )
+            .add_reified_edge(active_edge, source, target, weight, &self.model.state)
     }
 
     pub fn add_inactive_edge(&mut self, source: Timepoint, target: Timepoint, weight: W) -> Lit {
-        let v = self
+        let valid_edge = self.get_conjunctive_scope(source, target);
+        let active_edge = self
             .model
-            .new_bvar(format!("reif({:?} -- {} --> {:?})", source, weight, target));
-        let activation = v.true_lit();
-        self.add_reified_edge(activation, source, target, weight);
-        activation
+            .new_optional_bvar(valid_edge, format!("reif({:?} -- {} --> {:?})", source, weight, target))
+            .true_lit();
+
+        self.stn
+            .add_reified_edge(active_edge, source, target, weight, &self.model.state);
+        active_edge
     }
 
     // add delay between optional variables
-    pub fn add_delay(&mut self, a: Timepoint, b: Timepoint, delay: W) {
-        // edge a <--- -1 --- b
-        let a_to_b = can_propagate(&self.model.state, a, b);
-        let b_to_a = can_propagate(&self.model.state, b, a);
-        let presence = edge_presence(&self.model.state, a, b);
-        self.add_optional_true_edge(b, a, -delay, b_to_a, a_to_b, presence);
+    pub fn add_delay(&mut self, a: impl Into<Timepoint>, b: impl Into<Timepoint>, delay: W) {
+        self.add_edge(b.into(), a.into(), -delay);
+    }
+
+    /// Returns a literal that is true iff both timepoints are present.
+    fn get_conjunctive_scope(&mut self, a: Timepoint, b: Timepoint) -> Lit {
+        let pa = self.model.state.presence(a);
+        let pb = self.model.state.presence(b);
+        self.model.get_conjunctive_scope(&[pa, pb])
     }
 
     pub fn mark_active(&mut self, edge: Lit) {
diff --git a/tnet/src/theory.rs b/tnet/src/theory.rs
index 2e30a443..45e97445 100644
--- a/tnet/src/theory.rs
+++ b/tnet/src/theory.rs
@@ -1,30 +1,35 @@
-use crate::distances::DijkstraState;
+mod contraint_db;
+mod distances;
+mod edges;
+
 use crate::theory::Event::{EdgeActivated, EdgeAdded};
 use aries_backtrack::Backtrack;
 use aries_backtrack::{DecLvl, ObsTrailCursor, Trail};
 use aries_collections::ref_store::{RefMap, RefVec};
 use aries_collections::set::RefSet;
-use aries_model::bounds::{BoundValue, BoundValueAdd, Lit, VarBound, Watches};
-use aries_model::lang::normal_form::{NFLeq, NFOptLeq};
+use aries_core::state::*;
+use aries_core::*;
+use aries_model::lang::normal_form::NFLeq;
 use aries_model::lang::reification::{downcast, Expr};
-use aries_model::lang::{IntCst, VarRef};
-use aries_model::state::Domains;
-use aries_model::state::*;
-use aries_model::WriterId;
 use aries_solver::solver::BindingResult;
 use aries_solver::{Bind, Contradiction, Theory};
+use contraint_db::*;
+use distances::DijkstraState;
+use edges::*;
 use env_param::EnvParam;
-use std::collections::{HashMap, VecDeque};
+use std::collections::VecDeque;
 use std::convert::*;
-use std::hash::Hash;
 use std::marker::PhantomData;
-use std::ops::Index;
-use std::ops::IndexMut;
 use std::str::FromStr;
 
-type ModelEvent = aries_model::state::Event;
+pub use edges::EdgeId;
+
+type ModelEvent = aries_core::state::Event;
 
+/// A temporal reference in an STN, i.e., reference to an absolute time.
 pub type Timepoint = VarRef;
+
+/// The edge weight of an STN, i.e., a fixed duration.
 pub type W = IntCst;
 
 pub static STN_THEORY_PROPAGATION: EnvParam<TheoryPropagationLevel> =
@@ -72,7 +77,7 @@ impl FromStr for TheoryPropagationLevel {
             "edges" => Ok(TheoryPropagationLevel::Edges),
             "full" => Ok(TheoryPropagationLevel::Full),
             x => Err(format!(
-                "Unknown theory propagation level: {}. Valid options: none, bounds, edges, full",
+                "Unknown theory propagation level: {}. Valid options: none, literals, edges, full",
                 x
             )),
         }
@@ -105,415 +110,12 @@ impl Default for StnConfig {
     }
 }
 
-/// A unique identifier for an edge in the STN.
-/// An edge and its negation share the same `base_id` but differ by the `is_negated` property.
-///
-/// For instance, valid edge ids:
-///  -  a - b <= 10
-///    - base_id: 3
-///    - negated: false
-///  - a - b > 10       # negation of the previous one
-///    - base_id: 3     # same
-///    - negated: true  # inverse
-///  - a -b <= 20       # unrelated
-///    - base_id: 4
-///    - negated: false
-#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash, Debug)]
-pub struct EdgeId(u32);
-impl EdgeId {
-    #[inline]
-    fn new(base_id: u32, negated: bool) -> EdgeId {
-        if negated {
-            EdgeId((base_id << 1) + 1)
-        } else {
-            EdgeId(base_id << 1)
-        }
-    }
-
-    #[inline]
-    pub fn base_id(&self) -> u32 {
-        self.0 >> 1
-    }
-
-    #[inline]
-    pub fn is_negated(&self) -> bool {
-        self.0 & 0x1 == 1
-    }
-
-    /// Id of the forward (from source to target) view of this edge
-    fn forward(self) -> DirEdge {
-        DirEdge::forward(self)
-    }
-
-    /// Id of the backward view (from target to source) of this edge
-    fn backward(self) -> DirEdge {
-        DirEdge::backward(self)
-    }
-}
-
-impl std::ops::Not for EdgeId {
-    type Output = Self;
-
-    #[inline]
-    fn not(self) -> Self::Output {
-        EdgeId(self.0 ^ 0x1)
-    }
-}
-
-impl From<EdgeId> for u32 {
-    fn from(e: EdgeId) -> Self {
-        e.0
-    }
-}
-impl From<u32> for EdgeId {
-    fn from(id: u32) -> Self {
-        EdgeId(id)
-    }
-}
-
-impl From<EdgeId> for usize {
-    fn from(e: EdgeId) -> Self {
-        e.0 as usize
-    }
-}
-impl From<usize> for EdgeId {
-    fn from(id: usize) -> Self {
-        EdgeId(id as u32)
-    }
-}
-
-/// An edge in the STN, representing the constraint `target - source <= weight`
-/// An edge can be either in canonical form or in negated form.
-/// Given to edges (tgt - src <= w) and (tgt -src > w) one will be in canonical form and
-/// the other in negated form.
-#[derive(Copy, Clone, Debug, PartialOrd, Ord, PartialEq, Eq, Hash)]
-pub struct Edge {
-    pub source: Timepoint,
-    pub target: Timepoint,
-    pub weight: W,
-}
-
-impl Edge {
-    pub fn new(source: Timepoint, target: Timepoint, weight: W) -> Edge {
-        Edge { source, target, weight }
-    }
-
-    fn is_negated(&self) -> bool {
-        !self.is_canonical()
-    }
-
-    fn is_canonical(&self) -> bool {
-        self.source < self.target || self.source == self.target && self.weight >= 0
-    }
-
-    // not(b - a <= 6)
-    //   = b - a > 6
-    //   = a -b < -6
-    //   = a - b <= -7
-    //
-    // not(a - b <= -7)
-    //   = a - b > -7
-    //   = b - a < 7
-    //   = b - a <= 6
-    fn negated(&self) -> Self {
-        Edge {
-            source: self.target,
-            target: self.source,
-            weight: -self.weight - 1,
-        }
-    }
-}
-
-/// A directional constraint representing the fact that an update on the `source` bound
-/// should be reflected on the `target` bound.
-///
-/// From a classical STN edge `source -- weight --> target` there will be two directional constraints:
-///   - ub(source) = X   implies   ub(target) <= X + weight
-///   - lb(target) = X   implies   lb(source) >= X - weight
-#[derive(Clone, Debug)]
-struct DirConstraint {
-    source: VarBound,
-    target: VarBound,
-    weight: BoundValueAdd,
-    /// Non-empty if the constraint active (participates in propagation)
-    /// If the enabler is Lit::TRUE, then the constraint can be assumed to be always active
-    enabler: Option<Lit>,
-    /// A set of potential enablers for this constraint.
-    /// The edge becomes active once one of its enablers becomes true
-    enablers: Vec<Lit>,
-}
-impl DirConstraint {
-    /// source <= X   =>   target <= X + weight
-    pub fn forward(edge: Edge) -> DirConstraint {
-        DirConstraint {
-            source: VarBound::ub(edge.source),
-            target: VarBound::ub(edge.target),
-            weight: BoundValueAdd::on_ub(edge.weight),
-            enabler: None,
-            enablers: vec![],
-        }
-    }
-
-    /// target >= X   =>   source >= X - weight
-    pub fn backward(edge: Edge) -> DirConstraint {
-        DirConstraint {
-            source: VarBound::lb(edge.target),
-            target: VarBound::lb(edge.source),
-            weight: BoundValueAdd::on_lb(-edge.weight),
-            enabler: None,
-            enablers: vec![],
-        }
-    }
-
-    pub fn as_edge(&self) -> Edge {
-        if self.source.is_ub() {
-            debug_assert!(self.target.is_ub());
-            Edge {
-                source: self.source.variable(),
-                target: self.target.variable(),
-                weight: self.weight.as_ub_add(),
-            }
-        } else {
-            debug_assert!(self.target.is_lb());
-            Edge {
-                source: self.target.variable(),
-                target: self.source.variable(),
-                weight: -self.weight.as_lb_add(),
-            }
-        }
-    }
-}
-
-/// A pair of constraints (a, b) where edge(a) = !edge(b)
-struct ConstraintPair {
-    /// constraint where the edge is in its canonical form
-    base_forward: DirConstraint,
-    base_backward: DirConstraint,
-    /// constraint corresponding to the negation of base
-    negated_forward: DirConstraint,
-    negated_backward: DirConstraint,
-}
-
-impl ConstraintPair {
-    pub fn new_inactives(edge: Edge) -> ConstraintPair {
-        let edge = if edge.is_canonical() { edge } else { edge.negated() };
-        ConstraintPair {
-            base_forward: DirConstraint::forward(edge),
-            base_backward: DirConstraint::backward(edge),
-            negated_forward: DirConstraint::forward(edge.negated()),
-            negated_backward: DirConstraint::backward(edge.negated()),
-        }
-    }
-}
-
-/// Represents an edge together with a particular propagation direction:
-///  - forward (source to target)
-///  - backward (target to source)
-#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Debug)]
-pub(crate) struct DirEdge(u32);
-
-impl DirEdge {
-    /// Forward view of the given edge
-    pub fn forward(e: EdgeId) -> Self {
-        DirEdge(u32::from(e) << 1)
-    }
-
-    /// Backward view of the given edge
-    pub fn backward(e: EdgeId) -> Self {
-        DirEdge((u32::from(e) << 1) + 1)
-    }
-
-    #[allow(unused)]
-    pub fn is_forward(self) -> bool {
-        (u32::from(self) & 0x1) == 0
-    }
-
-    /// The edge underlying this projection
-    pub fn edge(self) -> EdgeId {
-        EdgeId::from(self.0 >> 1)
-    }
-}
-impl From<DirEdge> for usize {
-    fn from(e: DirEdge) -> Self {
-        e.0 as usize
-    }
-}
-impl From<usize> for DirEdge {
-    fn from(u: usize) -> Self {
-        DirEdge(u as u32)
-    }
-}
-impl From<DirEdge> for u32 {
-    fn from(e: DirEdge) -> Self {
-        e.0
-    }
-}
-impl From<u32> for DirEdge {
-    fn from(u: u32) -> Self {
-        DirEdge(u)
-    }
-}
-
-/// Data structures that holds all active and inactive edges in the STN.
-/// Note that some edges might be represented even though they were never inserted if they are the
-/// negation of an inserted edge.
-#[derive(Clone)]
-struct ConstraintDb {
-    /// All directional constraints.
-    ///
-    /// Each time a new edge is create for `DirConstraint` will be added
-    /// - forward view of the canonical edge
-    /// - backward view of the canonical edge
-    /// - forward view of the negated edge
-    /// - backward view of the negated edge
-    constraints: RefVec<DirEdge, DirConstraint>,
-    /// Maps each canonical edge to its base ID.
-    lookup: HashMap<Edge, u32>,
-    /// Associates literals to the edges that should be activated when they become true
-    watches: Watches<DirEdge>,
-    edges: RefVec<VarBound, Vec<EdgeTarget>>,
-    /// Index of the next constraint that has not been returned yet by the `next_new_constraint` method.
-    next_new_constraint: usize,
-}
-
-#[derive(Copy, Clone, Debug)]
-struct EdgeTarget {
-    target: VarBound,
-    weight: BoundValueAdd,
-    /// Literal that is true if and only if the edge must be present in the network.
-    /// Note that handling of optional variables might allow and edge to propagate even it is not known
-    /// to be present yet.
-    presence: Lit,
-}
-
-impl ConstraintDb {
-    pub fn new() -> ConstraintDb {
-        ConstraintDb {
-            constraints: Default::default(),
-            lookup: HashMap::new(),
-            watches: Default::default(),
-            edges: Default::default(),
-            next_new_constraint: 0,
-        }
-    }
-
-    /// A function that acts as a one time iterator over constraints.
-    /// It can be used to check if new constraints have been added since last time this method was called.
-    pub fn next_new_constraint(&mut self) -> Option<&DirConstraint> {
-        if self.next_new_constraint < self.constraints.len() {
-            let out = &self.constraints[self.next_new_constraint.into()];
-            self.next_new_constraint += 1;
-            Some(out)
-        } else {
-            None
-        }
-    }
-
-    /// Record the fact that, when `literal` becomes true, the given edge
-    /// should be made active in both directions.
-    pub fn add_enabler(&mut self, edge: EdgeId, literal: Lit) {
-        self.add_directed_enabler(edge.forward(), literal, Some(literal));
-        self.add_directed_enabler(edge.backward(), literal, Some(literal));
-    }
-
-    /// Record the fact that:
-    ///  - if `propagation_enabler` is true, then propagation of the directed edge should be made active
-    ///  - if the edge is inconsistent with the rest of the network, then the presence literal should be false.
-    pub fn add_directed_enabler(&mut self, edge: DirEdge, propagation_enabler: Lit, presence_literal: Option<Lit>) {
-        self.watches.add_watch(edge, propagation_enabler);
-        let constraint = &mut self.constraints[edge];
-        constraint.enablers.push(propagation_enabler);
-        self.edges.fill_with(constraint.source, Vec::new);
-        if let Some(presence_literal) = presence_literal {
-            self.edges[constraint.source].push(EdgeTarget {
-                target: constraint.target,
-                weight: constraint.weight,
-                presence: presence_literal,
-            });
-        }
-    }
-
-    pub fn potential_out_edges(&self, source: VarBound) -> &[EdgeTarget] {
-        if self.edges.contains(source) {
-            &self.edges[source]
-        } else {
-            &[]
-        }
-    }
-
-    fn find_existing(&self, edge: &Edge) -> Option<EdgeId> {
-        if edge.is_canonical() {
-            self.lookup.get(edge).map(|&id| EdgeId::new(id, false))
-        } else {
-            self.lookup.get(&edge.negated()).map(|&id| EdgeId::new(id, true))
-        }
-    }
-
-    /// Adds a new edge and return a pair (created, edge_id) where:
-    ///  - created is false if NO new edge was inserted (it was merge with an identical edge already in the DB)
-    ///  - edge_id is the id of the edge
-    pub fn push_edge(&mut self, source: Timepoint, target: Timepoint, weight: W) -> (bool, EdgeId) {
-        let edge = Edge::new(source, target, weight);
-        match self.find_existing(&edge) {
-            Some(id) => {
-                // edge already exists in the DB, return its id and say it wasn't created
-                debug_assert_eq!(self[DirEdge::forward(id)].as_edge(), edge);
-                debug_assert_eq!(self[DirEdge::backward(id)].as_edge(), edge);
-                (false, id)
-            }
-            None => {
-                // edge does not exist, record the corresponding pair and return the new id.
-                let pair = ConstraintPair::new_inactives(edge);
-                let base = pair.base_forward.as_edge();
-                let id1 = self.constraints.push(pair.base_forward);
-                let _ = self.constraints.push(pair.base_backward);
-                let id2 = self.constraints.push(pair.negated_forward);
-                let _ = self.constraints.push(pair.negated_backward);
-                self.lookup.insert(base, id1.edge().base_id());
-                debug_assert_eq!(id1.edge().base_id(), id2.edge().base_id());
-                let edge_id = if edge.is_negated() { id2 } else { id1 };
-                (true, edge_id.edge())
-            }
-        }
-    }
-
-    /// Removes the last created ConstraintPair in the DB. Note that this will remove the last edge that was
-    /// pushed and THAT WAS NOT UNIFIED with an existing edge (i.e. edge_push returned : (true, _)).
-    pub fn pop_last(&mut self) {
-        debug_assert_eq!(self.constraints.len() % 4, 0);
-        // remove the four edges (forward and backward) for both the base and negated edge
-        self.constraints.pop();
-        self.constraints.pop();
-        self.constraints.pop();
-        if let Some(c) = self.constraints.pop() {
-            debug_assert!(c.as_edge().is_canonical());
-            self.lookup.remove(&c.as_edge());
-        }
-    }
-
-    pub fn has_edge(&self, id: EdgeId) -> bool {
-        id.base_id() <= self.constraints.len() as u32
-    }
-}
-impl Index<DirEdge> for ConstraintDb {
-    type Output = DirConstraint;
-
-    fn index(&self, index: DirEdge) -> &Self::Output {
-        &self.constraints[index]
-    }
-}
-impl IndexMut<DirEdge> for ConstraintDb {
-    fn index_mut(&mut self, index: DirEdge) -> &mut Self::Output {
-        &mut self.constraints[index]
-    }
-}
-
 type BacktrackLevel = DecLvl;
 
 #[derive(Copy, Clone)]
 enum Event {
     EdgeAdded,
-    EdgeActivated(DirEdge),
+    EdgeActivated(PropagatorId),
     AddedTheoryPropagationCause,
 }
 
@@ -565,13 +167,13 @@ where
 /// (weight on edges) and absolute times (bound on nodes). It is the responsibility
 /// of the caller to ensure that no overflow occurs when adding an absolute and relative time,
 /// either by the choice of an appropriate type (e.g. saturating add) or by the choice of
-/// appropriate initial bounds.
+/// appropriate initial literals.
 #[derive(Clone)]
 pub struct StnTheory {
     pub config: StnConfig,
     constraints: ConstraintDb,
     /// Forward/Backward adjacency list containing active edges.
-    active_propagators: RefVec<VarBound, Vec<Propagator>>,
+    active_propagators: RefVec<VarBound, Vec<InlinedPropagator>>,
     pending_updates: RefSet<VarBound>,
     /// History of changes and made to the STN with all information necessary to undo them.
     trail: Trail<Event>,
@@ -583,7 +185,7 @@ pub struct StnTheory {
     /// When encountering an inconsistency, this vector will be cleared and
     /// a negative cycle will be constructed in it. The explanation returned
     /// will be a slice of this vector to avoid any allocation.
-    explanation: Vec<DirEdge>,
+    explanation: Vec<PropagatorId>,
     theory_propagation_causes: Vec<TheoryPropagationCause>,
     /// Internal data structure used by the `propagate` method to keep track of pending work.
     internal_propagate_queue: VecDeque<VarBound>,
@@ -600,16 +202,16 @@ enum TheoryPropagationCause {
     Path {
         source: VarBound,
         target: VarBound,
-        triggering_edge: DirEdge,
+        triggering_edge: PropagatorId,
     },
-    /// Theory propagation was triggered by the incompatibility of the two bounds with an edge in the graph.
+    /// Theory propagation was triggered by the incompatibility of the two literals with an edge in the graph.
     Bounds { source: Lit, target: Lit },
 }
 
 #[derive(Copy, Clone)]
 pub(crate) enum ModelUpdateCause {
     /// The update was caused by an edge propagation
-    EdgePropagation(DirEdge),
+    EdgePropagation(PropagatorId),
     /// index in the trail of the TheoryPropagationCause
     TheoryPropagation(u32),
 }
@@ -617,7 +219,7 @@ pub(crate) enum ModelUpdateCause {
 impl From<u32> for ModelUpdateCause {
     fn from(enc: u32) -> Self {
         if (enc & 0x1) == 0 {
-            ModelUpdateCause::EdgePropagation(DirEdge::from(enc >> 1))
+            ModelUpdateCause::EdgePropagation(PropagatorId::from(enc >> 1))
         } else {
             ModelUpdateCause::TheoryPropagation(enc >> 1)
         }
@@ -633,17 +235,19 @@ impl From<ModelUpdateCause> for u32 {
     }
 }
 
+/// Contains the id of a propagator as well as its `target` and `weight` fields that
+/// are inlined to facilitate propagation.
 #[derive(Copy, Clone, Debug)]
-struct Propagator {
+struct InlinedPropagator {
     target: VarBound,
     weight: BoundValueAdd,
-    id: DirEdge,
+    id: PropagatorId,
 }
 
 #[derive(Copy, Clone)]
 enum ActivationEvent {
     /// Should activate the given edge, enabled by this literal
-    ToActivate(DirEdge, Lit),
+    ToEnable(PropagatorId, Enabler),
 }
 
 impl StnTheory {
@@ -672,7 +276,7 @@ impl StnTheory {
     }
 
     pub fn reserve_timepoint(&mut self) {
-        // add slots for the propagators of both bounds
+        // add slots for the propagators of both literals
         self.active_propagators.push(Vec::new());
         self.active_propagators.push(Vec::new());
     }
@@ -686,81 +290,73 @@ impl StnTheory {
         domains: &Domains,
     ) -> EdgeId {
         let e = self.add_inactive_constraint(source.into(), target.into(), weight).0;
+        // record all possible propagators
 
-        if domains.entails(literal) {
-            assert_eq!(domains.entailing_level(literal), DecLvl::ROOT);
-            self.mark_active(e, literal);
-        } else if domains.entails(!literal) {
-            assert_eq!(domains.entailing_level(!literal), DecLvl::ROOT);
-            self.mark_active(!e, !literal);
-        } else {
-            self.constraints.add_enabler(e, literal);
-            self.constraints.add_enabler(!e, !literal);
-        }
+        // normal edge:  literal <=> source ---- weight ---> target
+        self.record_propagator(e.forward(), literal, domains);
+        self.record_propagator(e.backward(), literal, domains);
+
+        // negated edge: !literal = source <--- (-weight-1) ---- target
+        self.record_propagator((!e).forward(), !literal, domains);
+        self.record_propagator((!e).backward(), !literal, domains);
 
         e
     }
 
-    /// Adds an edge `source --- weight ---> target` to the network that must hold if
-    /// both `source` and `target` are present.
-    ///
-    /// To control propagation the following literals are provided:
-    ///  - `forward_prop`: true if propagation is allowed from `source` to `target`
-    ///    This is typically equivalent to   `present(source) => present(target)`
-    ///  - `backward_prop`: true if propagation is allowed from `target` to `source`
-    ///    This is typically equivalent to   `present(target) => present(source)`
-    ///  - `presence`: true if both timepoints are present, and thus the edge is active.
-    ///    equivalent to `present(source) and present(target)`. This parameter is optional
-    ///    and is used in theory propagation to deactivate the edge.
-    #[allow(clippy::too_many_arguments)]
-    pub fn add_optional_true_edge(
-        &mut self,
-        source: impl Into<Timepoint>,
-        target: impl Into<Timepoint>,
-        weight: W,
-        forward_prop: Lit,
-        backward_prop: Lit,
-        presence: Option<Lit>,
-        domains: &Domains,
-    ) -> EdgeId {
-        let e = self.add_inactive_constraint(source.into(), target.into(), weight).0;
+    /// Creates and record a new propagator associated with the given [DirEdge], making sure
+    /// to set up the watches to enable it when it becomes active and valid.
+    fn record_propagator(&mut self, prop: PropagatorId, active: Lit, domains: &Domains) {
+        let propagator_source = self.constraints[prop].source.variable();
+        let propagator_target = self.constraints[prop].target.variable();
 
-        self.constraints
-            .add_directed_enabler(e.forward(), forward_prop, presence);
-        if domains.entails(forward_prop) {
-            assert_eq!(domains.entailing_level(forward_prop), DecLvl::ROOT);
-            self.pending_activations
-                .push_back(ActivationEvent::ToActivate(e.forward(), forward_prop));
-        }
-        self.constraints
-            .add_directed_enabler(e.backward(), backward_prop, presence);
-        if domains.entails(backward_prop) {
-            assert_eq!(domains.entailing_level(backward_prop), DecLvl::ROOT);
-            self.pending_activations
-                .push_back(ActivationEvent::ToActivate(e.backward(), backward_prop));
-        }
+        // literal that is true if the edge is within its validity scope (i.e. both timepoints are present)
+        // edge_valid <=> presence(source) & presence(target)
+        let edge_valid = domains.presence(active.variable());
+        debug_assert!(domains.implies(edge_valid, domains.presence(propagator_source)));
+        debug_assert!(domains.implies(edge_valid, domains.presence(propagator_target)));
 
-        e
-    }
+        // the propagator is valid when `presence(target) => edge_valid`.
+        // This is because in this case, the modification to the target's domain are only meaningful if the edge is present.
+        // Once the propagator is valid, it can be propagated as soon as its `active` literal becomes true.
 
-    /// Marks an edge as active and enqueue it for propagation.
-    /// No changes are committed to the network by this function until a call to `propagate_all()`
-    pub fn mark_active(&mut self, edge: EdgeId, enabler: Lit) {
-        debug_assert!(self.constraints.has_edge(edge));
-        self.pending_activations
-            .push_back(ActivationEvent::ToActivate(DirEdge::forward(edge), enabler));
-        self.pending_activations
-            .push_back(ActivationEvent::ToActivate(DirEdge::backward(edge), enabler));
+        // determine a literal that is true iff the propagator is valid
+        let propagator_valid = if domains.implies(domains.presence(propagator_target), edge_valid) {
+            // it is statically known that `presence(target) => edge_valid`,
+            // the propagator is always valid
+            Lit::TRUE
+        } else {
+            // given that `presence(source) & presence(target) <=> edge_valid`, we can infer that the propagator becomes valid
+            // (i.e. `presence(target) => edge_valid` holds) when `presence(source)` becomes true
+            domains.presence(propagator_source)
+        };
+        let enabler = Enabler::new(active, propagator_valid);
+
+        // Add the propagator, with different modalities depending on whether it is currently enabled or not.
+        // Note that we should make sure that when backtracking beyond the current decision level, we should deactivate the edge.
+        if domains.entails(!active) || domains.entails(!edge_valid) {
+            // do nothing as the propagator can never be active/present
+        } else if domains.entails(active) && domains.entails(propagator_valid) {
+            // propagator is always active in the current and following decision levels, enqueue it for activation.
+            self.pending_activations
+                .push_back(ActivationEvent::ToEnable(prop, enabler));
+        } else {
+            if domains.current_decision_level() != DecLvl::ROOT {
+                // FIXME: when backtracking, we should remove this edge (or at least ensure that it is definitely deactivated)
+                println!("WARNING: adding a dynamically enabled edge beyond the root decision level is unsupported.")
+            }
+            self.constraints.add_propagator_enabler(prop, enabler);
+        }
     }
 
-    fn build_contradiction(&self, culprits: &[DirEdge], model: &Domains) -> Contradiction {
+    fn build_contradiction(&self, culprits: &[PropagatorId], model: &Domains) -> Contradiction {
         let mut expl = Explanation::with_capacity(culprits.len());
         for &edge in culprits {
             debug_assert!(self.active(edge));
-            let literal = self.constraints[edge].enabler;
-            let literal = literal.expect("No entailed enabler for this edge");
-            debug_assert!(model.entails(literal));
-            expl.push(literal);
+            let enabler = self.constraints[edge].enabler;
+            let enabler = enabler.expect("No entailed enabler for this edge");
+            debug_assert!(model.entails(enabler.active) && model.entails(enabler.valid));
+            expl.push(enabler.active);
+            expl.push(model.presence(enabler.active.variable()));
         }
         Contradiction::Explanation(expl)
     }
@@ -768,7 +364,7 @@ impl StnTheory {
     fn explain_bound_propagation(
         &self,
         event: Lit,
-        propagator: DirEdge,
+        propagator: PropagatorId,
         model: &Domains,
         out_explanation: &mut Explanation,
     ) {
@@ -777,8 +373,9 @@ impl StnTheory {
         let val = event.bound_value();
         debug_assert_eq!(event.affected_bound(), c.target);
 
-        let literal = self.constraints[propagator].enabler.expect("inactive constraint");
-        out_explanation.push(literal);
+        let enabler = self.constraints[propagator].enabler.expect("inactive constraint");
+        out_explanation.push(enabler.active);
+        out_explanation.push(model.presence(enabler.active.variable()));
 
         let cause = Lit::from_parts(c.source, val - c.weight);
         debug_assert!(model.entails(cause));
@@ -786,7 +383,7 @@ impl StnTheory {
         if self.config.deep_explanation {
             // function that return the stn propagator responsible for this literal being set,
             // of None if it was not set by a bound propagation of the STN.
-            let propagator_of = |lit: Lit, model: &Domains| -> Option<DirEdge> {
+            let propagator_of = |lit: Lit, model: &Domains| -> Option<PropagatorId> {
                 if let Some(event_index) = model.implying_event(lit) {
                     let event = model.get_event(event_index);
                     match event.cause.as_external_inference() {
@@ -806,7 +403,8 @@ impl StnTheory {
             while let Some(propagator) = propagator_of(latest_trigger, model) {
                 let c = &self.constraints[propagator];
                 let propagator_enabler = c.enabler.expect("inactive edge");
-                out_explanation.push(propagator_enabler);
+                out_explanation.push(propagator_enabler.active);
+                out_explanation.push(model.presence(propagator_enabler.active.variable()));
                 debug_assert_eq!(latest_trigger.affected_bound(), c.target);
                 latest_trigger = Lit::from_parts(c.source, latest_trigger.bound_value() - c.weight);
                 debug_assert!(model.entails(latest_trigger));
@@ -833,8 +431,9 @@ impl StnTheory {
                 let path = self.theory_propagation_path(source, target, triggering_edge, model);
 
                 for edge in path {
-                    let literal = self.constraints[edge].enabler.expect("inactive constraint");
-                    out_explanation.push(literal);
+                    let enabler = self.constraints[edge].enabler.expect("inactive constraint");
+                    out_explanation.push(enabler.active);
+                    out_explanation.push(model.presence(enabler.active.variable()));
                 }
             }
             TheoryPropagationCause::Bounds { source, target } => {
@@ -847,7 +446,7 @@ impl StnTheory {
 
     /// Propagates all edges that have been marked as active since the last propagation.
     pub fn propagate_all(&mut self, model: &mut Domains) -> Result<(), Contradiction> {
-        // in first propagation, process each edge once to check if it can be added to the model based on the bounds
+        // in first propagation, process each edge once to check if it can be added to the model based on the literals
         // of its extremities. If it is not the case, make its enablers false.
         // This step is equivalent to "bound theory propagation" but need to be made independently because
         // we do not get events for the initial domain of the variables.
@@ -871,7 +470,7 @@ impl StnTheory {
                             .inference(ModelUpdateCause::TheoryPropagation(cause_index as u32));
                         // make all enablers false
                         for &l in &c.enablers {
-                            model.set(!l, cause)?;
+                            model.set(!l.active, cause)?;
                         }
                     }
                 }
@@ -879,16 +478,18 @@ impl StnTheory {
         }
 
         while self.model_events.num_pending(model.trail()) > 0 || !self.pending_activations.is_empty() {
-            // start by propagating all bounds changes before considering the new edges.
+            // start by propagating all literals changes before considering the new edges.
             // This is necessary because cycle detection on the insertion of a new edge requires
             // a consistent STN and no interference of external bound updates.
             while let Some(ev) = self.model_events.pop(model.trail()).copied() {
                 let literal = ev.new_literal();
-                for edge in self.constraints.watches.watches_on(literal) {
+                for (enabler, edge) in self.constraints.enabled_by(literal) {
                     // mark active
                     debug_assert!(self.constraints.has_edge(edge.edge()));
-                    self.pending_activations
-                        .push_back(ActivationEvent::ToActivate(edge, literal));
+                    if model.entails(enabler.active) && model.entails(enabler.valid) {
+                        self.pending_activations
+                            .push_back(ActivationEvent::ToEnable(edge, enabler));
+                    }
                 }
                 if self.config.theory_propagation.bounds() {
                     self.theory_propagate_bound(literal, model)?;
@@ -905,7 +506,7 @@ impl StnTheory {
                 self.propagate_bound_change(literal, model)?;
             }
             while let Some(event) = self.pending_activations.pop_front() {
-                let ActivationEvent::ToActivate(edge, enabler) = event;
+                let ActivationEvent::ToEnable(edge, enabler) = event;
                 let c = &mut self.constraints[edge];
                 if c.enabler.is_none() {
                     // edge is currently inactive
@@ -925,7 +526,7 @@ impl StnTheory {
                     } else {
                         debug_assert_ne!(c.source, c.target);
 
-                        self.active_propagators[c.source].push(Propagator {
+                        self.active_propagators[c.source].push(InlinedPropagator {
                             target: c.target,
                             weight: c.weight,
                             id: edge,
@@ -1011,7 +612,7 @@ impl StnTheory {
         (id, created)
     }
 
-    fn active(&self, e: DirEdge) -> bool {
+    fn active(&self, e: PropagatorId) -> bool {
         self.constraints[e].enabler.is_some()
     }
 
@@ -1039,7 +640,7 @@ impl StnTheory {
 
     /// Implementation of [Cesta96]
     /// It propagates a **newly_inserted** edge in a **consistent** STN.
-    fn propagate_new_edge(&mut self, new_edge: DirEdge, model: &mut Domains) -> Result<(), Contradiction> {
+    fn propagate_new_edge(&mut self, new_edge: PropagatorId, model: &mut Domains) -> Result<(), Contradiction> {
         let c = &self.constraints[new_edge];
         debug_assert_ne!(c.source, c.target, "This algorithm does not support self loops.");
         let cause = self.identity.inference(ModelUpdateCause::EdgePropagation(new_edge));
@@ -1118,7 +719,8 @@ impl StnTheory {
             curr = c.source;
             // cycle_length += c.edge.weight;
             let trigger = self.constraints[edge].enabler.expect("inactive constraint");
-            expl.push(trigger);
+            expl.push(trigger.active);
+            expl.push(model.presence(trigger.active.variable()));
 
             if curr == vb {
                 // debug_assert!(cycle_length < 0);
@@ -1129,7 +731,7 @@ impl StnTheory {
 
     pub fn print_stats(&self) {
         println!("# nodes: {}", self.num_nodes());
-        println!("# constraints: {}", self.constraints.constraints.len());
+        println!("# propagators: {}", self.constraints.num_propagators());
         println!("# propagations: {}", self.stats.num_propagations);
         println!("# domain updates: {}", self.stats.distance_updates);
     }
@@ -1203,7 +805,7 @@ impl StnTheory {
     /// Then we check if there exist an inactive edge BA where `weight(BA) + dist(AB) < 0`.
     /// For each such edge, we set its enabler to false since its addition would result in a negative cycle.
     #[inline(never)]
-    fn theory_propagate_edge(&mut self, edge: DirEdge, model: &mut Domains) -> Result<(), Contradiction> {
+    fn theory_propagate_edge(&mut self, edge: PropagatorId, model: &mut Domains) -> Result<(), Contradiction> {
         let constraint = &self.constraints[edge];
         let target = constraint.target;
         let source = constraint.source;
@@ -1310,7 +912,7 @@ impl StnTheory {
         &self,
         source: VarBound,
         target: VarBound,
-        through_edge: DirEdge,
+        through_edge: PropagatorId,
         model: &Domains,
         successors: &DijkstraState,
         predecessors: &DijkstraState,
@@ -1419,7 +1021,7 @@ impl StnTheory {
         to: VarBound,
         model: &Domains,
         state: &mut DijkstraState,
-        out: &mut Vec<DirEdge>,
+        out: &mut Vec<PropagatorId>,
     ) {
         state.clear();
         state.enqueue(from, BoundValueAdd::ZERO, None);
@@ -1488,9 +1090,9 @@ impl StnTheory {
         &self,
         source: VarBound,
         target: VarBound,
-        through_edge: DirEdge,
+        through_edge: PropagatorId,
         model: &Domains,
-    ) -> Vec<DirEdge> {
+    ) -> Vec<PropagatorId> {
         let mut path = Vec::with_capacity(8);
 
         let e = &self.constraints[through_edge];
@@ -1567,78 +1169,23 @@ impl Backtrack for StnTheory {
 impl Bind for StnTheory {
     fn bind(&mut self, literal: Lit, expr: &Expr, doms: &mut Domains) -> BindingResult {
         if let Some(&NFLeq { lhs, rhs, rhs_add }) = downcast(expr) {
-            // lhs  <= rhs + rhs_add    <=>   lhs - rhs <= rhs_add
             self.add_reified_edge(literal, rhs, lhs, rhs_add, doms);
             BindingResult::Enforced
-        } else if let Some(&NFOptLeq { lhs, rhs, rhs_add }) = downcast(expr) {
-            if doms.entails(literal) {
-                let a_to_b = can_propagate(doms, lhs, rhs);
-                let b_to_a = can_propagate(doms, rhs, lhs);
-                let presence = edge_presence(doms, lhs, rhs);
-                self.add_optional_true_edge(rhs, lhs, rhs_add, b_to_a, a_to_b, presence, doms);
-                BindingResult::Enforced
-            } else if doms.entails(!literal) {
-                // this constraint is always false, post the opposite
-                // !(lhs <= rhs + rhs_add) <=> lhs > rhs + rhs_add
-                //                         <=> - rhs_add > rhs - lhs
-                //                         <=> rhs -lhs <= -rhs_add -1
-                let delay = -rhs_add - 1;
-                let a = rhs;
-                let b = lhs;
-
-                let a_to_b = can_propagate(doms, a, b);
-                let b_to_a = can_propagate(doms, b, a);
-                let presence = edge_presence(doms, a, b);
-                self.add_optional_true_edge(b, a, delay, b_to_a, a_to_b, presence, doms);
-                BindingResult::Enforced
-            } else {
-                BindingResult::Unsupported
-            }
         } else {
             BindingResult::Unsupported
         }
     }
 }
 
-/// Returns the literal that is true if propagation is allowed from one optional variable to another.
-///
-/// # Panics
-/// If the the two variable are in unrelated scopes.
-pub(crate) fn can_propagate(doms: &Domains, source: Timepoint, target: Timepoint) -> Lit {
-    // lit = (from ---> to)    ,  we can if (lit != false) && p(from) => p(to)
-    if doms.only_present_with(target, source) {
-        // p(target) => p(source)
-        // !p(source) => !p(target)
-        Lit::TRUE
-    } else if doms.only_present_with(source, target) {
-        // p(source) => p(target)
-        doms.presence(source)
-    } else {
-        panic!()
-    }
-}
-
-/// Returns a literal that is true iff both optional variables are present.
-/// Returns None if it was not possible to find such a literal
-pub(crate) fn edge_presence(doms: &Domains, var1: Timepoint, var2: Timepoint) -> Option<Lit> {
-    if doms.only_present_with(var2, var1) {
-        // p(var2) => p(var1)
-        Some(doms.presence(var2))
-    } else if doms.only_present_with(var1, var2) {
-        // p(var1) => p(var2)
-        Some(doms.presence(var1))
-    } else {
-        None
-    }
-}
-
 #[cfg(test)]
 mod tests {
-    use super::*;
-    use crate::stn::Stn;
     use aries_model::extensions::AssignmentExt;
     use aries_model::lang::IVar;
 
+    use crate::stn::Stn;
+
+    use super::*;
+
     #[test]
     fn test_edge_id_conversions() {
         fn check_rountrip(i: u32) {
@@ -1803,9 +1350,7 @@ mod tests {
         let prez_b = stn.model.new_presence_variable(prez_a, "prez_b").true_lit();
         let b = stn.model.new_optional_ivar(0, 10, prez_b, "b");
 
-        let a_implies_b = prez_b;
-        let b_implies_a = Lit::TRUE;
-        stn.add_optional_true_edge(b, a, 0, a_implies_b, b_implies_a, Some(a_implies_b));
+        stn.add_delay(a, b, 0);
 
         stn.propagate_all()?;
         stn.model.state.set_lb(b, 1, Cause::Decision)?;
@@ -1842,7 +1387,7 @@ mod tests {
                 .true_lit();
             let var = stn.model.new_optional_ivar(0, 20, prez, format!("var_{}", i));
             if i > 0 {
-                stn.add_delay(vars[i - 1].1.into(), var.into(), 1);
+                stn.add_delay(vars[i - 1].1, var, 1);
             }
             vars.push((prez, var));
             context = prez;
@@ -1869,7 +1414,7 @@ mod tests {
 
     #[test]
     fn test_theory_propagation_edges_simple() -> Result<(), Contradiction> {
-        let stn = &mut Stn::with_config(StnConfig {
+        let stn = &mut Stn::new_with_config(StnConfig {
             theory_propagation: TheoryPropagationLevel::Edges,
             ..Default::default()
         });
@@ -2047,7 +1592,7 @@ mod tests {
 
     #[test]
     fn test_theory_propagation_edges() -> Result<(), Contradiction> {
-        let stn = &mut Stn::with_config(StnConfig {
+        let stn = &mut Stn::new_with_config(StnConfig {
             theory_propagation: TheoryPropagationLevel::Edges,
             ..Default::default()
         });
@@ -2113,7 +1658,7 @@ mod tests {
 
     #[test]
     fn test_theory_propagation_bounds() -> Result<(), Contradiction> {
-        let stn = &mut Stn::with_config(StnConfig {
+        let stn = &mut Stn::new_with_config(StnConfig {
             theory_propagation: TheoryPropagationLevel::Bounds,
             ..Default::default()
         });
diff --git a/tnet/src/theory/contraint_db.rs b/tnet/src/theory/contraint_db.rs
new file mode 100644
index 00000000..3a7b4a5d
--- /dev/null
+++ b/tnet/src/theory/contraint_db.rs
@@ -0,0 +1,200 @@
+use crate::theory::edges::*;
+use crate::theory::{Timepoint, W};
+use aries_collections::ref_store::RefVec;
+use aries_core::literals::Watches;
+use aries_core::{Lit, VarBound};
+use std::collections::HashMap;
+use std::ops::{Index, IndexMut};
+
+/// Enabling information for a propagator.
+/// A propagator should be enabled iff both literals `active` and `valid` are true.
+#[derive(Clone, Copy, Debug)]
+pub(crate) struct Enabler {
+    /// A literal that is true (but not necessarily present) when the propagator must be active if present
+    pub active: Lit,
+    /// A literal that is true when the propagator is within its validity scope, i.e.,
+    /// when is known to be sound to propagate a change from the source to the target.
+    ///
+    /// In the simplest case, we have `valid = presence(active)` since by construction
+    /// `presence(active)` is true iff both variables of the constraint are present.
+    ///
+    /// `valid` might a more specific literal but always with the constraints that
+    /// `presence(active) => valid`
+    pub valid: Lit,
+}
+
+impl Enabler {
+    pub fn new(active: Lit, valid: Lit) -> Enabler {
+        Enabler { active, valid }
+    }
+}
+
+/// Data structures that holds all active and inactive edges in the STN.
+/// Note that some edges might be represented even though they were never inserted if they are the
+/// negation of an inserted edge.
+#[derive(Clone)]
+pub(crate) struct ConstraintDb {
+    /// All directional constraints.
+    ///
+    /// Each time a new edge is created four `DirConstraint` will be added
+    /// - forward view of the canonical edge
+    /// - backward view of the canonical edge
+    /// - forward view of the negated edge
+    /// - backward view of the negated edge
+    constraints: RefVec<PropagatorId, Propagator>,
+    /// Maps each canonical edge to its base ID.
+    lookup: HashMap<Edge, u32>,
+    /// Associates literals to the edges that should be activated when they become true
+    watches: Watches<(Enabler, PropagatorId)>,
+    edges: RefVec<VarBound, Vec<PropagatorTarget>>,
+    /// Index of the next constraint that has not been returned yet by the `next_new_constraint` method.
+    next_new_constraint: usize,
+}
+
+impl ConstraintDb {
+    pub fn new() -> ConstraintDb {
+        ConstraintDb {
+            constraints: Default::default(),
+            lookup: HashMap::new(),
+            watches: Default::default(),
+            edges: Default::default(),
+            next_new_constraint: 0,
+        }
+    }
+
+    pub fn num_propagators(&self) -> usize {
+        self.constraints.len()
+    }
+
+    /// A function that acts as a one time iterator over constraints.
+    /// It can be used to check if new constraints have been added since last time this method was called.
+    pub fn next_new_constraint(&mut self) -> Option<&Propagator> {
+        if self.next_new_constraint < self.constraints.len() {
+            let out = &self.constraints[self.next_new_constraint.into()];
+            self.next_new_constraint += 1;
+            Some(out)
+        } else {
+            None
+        }
+    }
+
+    /// Record the fact that:
+    ///  - if `enabler` holds (both literals are true), then the propagator should be enabled
+    ///  - if the `propagator` is inconsistent with the rest of the network, then the `enabler.active`
+    ///    literal should be made false.
+    pub fn add_propagator_enabler(&mut self, propagator: PropagatorId, enabler: Enabler) {
+        // watch both the `active` and the `valid` literal.
+        // when notified that one becomes true, we will need to check that the other is true before taking any action
+        self.watches.add_watch((enabler, propagator), enabler.active);
+        self.watches.add_watch((enabler, propagator), enabler.valid);
+        let constraint = &mut self.constraints[propagator];
+        constraint.enablers.push(enabler);
+        self.edges.fill_with(constraint.source, Vec::new);
+
+        self.edges[constraint.source].push(PropagatorTarget {
+            target: constraint.target,
+            weight: constraint.weight,
+            presence: enabler.active,
+        });
+    }
+
+    pub fn potential_out_edges(&self, source: VarBound) -> &[PropagatorTarget] {
+        if self.edges.contains(source) {
+            &self.edges[source]
+        } else {
+            &[]
+        }
+    }
+
+    fn find_existing(&self, edge: &Edge) -> Option<EdgeId> {
+        if edge.is_canonical() {
+            self.lookup.get(edge).map(|&id| EdgeId::new(id, false))
+        } else {
+            self.lookup.get(&edge.negated()).map(|&id| EdgeId::new(id, true))
+        }
+    }
+
+    /// Adds a new edge and return a pair (created, edge_id) where:
+    ///  - created is false if NO new edge was inserted (it was merge with an identical edge already in the DB)
+    ///  - edge_id is the id of the edge
+    pub fn push_edge(&mut self, source: Timepoint, target: Timepoint, weight: W) -> (bool, EdgeId) {
+        let edge = Edge::new(source, target, weight);
+        match self.find_existing(&edge) {
+            Some(id) => {
+                // edge already exists in the DB, return its id and say it wasn't created
+                debug_assert_eq!(self[PropagatorId::forward(id)].as_edge(), edge);
+                debug_assert_eq!(self[PropagatorId::backward(id)].as_edge(), edge);
+                (false, id)
+            }
+            None => {
+                // edge does not exist, record the corresponding pair and return the new id.
+                let pair = ConstraintPair::new_inactives(edge);
+                let base = pair.base_forward.as_edge();
+                let id1 = self.constraints.push(pair.base_forward);
+                let _ = self.constraints.push(pair.base_backward);
+                let id2 = self.constraints.push(pair.negated_forward);
+                let _ = self.constraints.push(pair.negated_backward);
+                self.lookup.insert(base, id1.edge().base_id());
+                debug_assert_eq!(id1.edge().base_id(), id2.edge().base_id());
+                let edge_id = if edge.is_negated() { id2 } else { id1 };
+                (true, edge_id.edge())
+            }
+        }
+    }
+
+    pub fn enabled_by(&self, literal: Lit) -> impl Iterator<Item = (Enabler, PropagatorId)> + '_ {
+        self.watches.watches_on(literal)
+    }
+
+    /// Removes the last created ConstraintPair in the DB. Note that this will remove the last edge that was
+    /// pushed and THAT WAS NOT UNIFIED with an existing edge (i.e. edge_push returned : (true, _)).
+    pub fn pop_last(&mut self) {
+        debug_assert_eq!(self.constraints.len() % 4, 0);
+        // remove the four edges (forward and backward) for both the base and negated edge
+        self.constraints.pop();
+        self.constraints.pop();
+        self.constraints.pop();
+        if let Some(c) = self.constraints.pop() {
+            debug_assert!(c.as_edge().is_canonical());
+            self.lookup.remove(&c.as_edge());
+        }
+    }
+
+    pub fn has_edge(&self, id: EdgeId) -> bool {
+        id.base_id() <= self.constraints.len() as u32
+    }
+}
+impl Index<PropagatorId> for ConstraintDb {
+    type Output = Propagator;
+
+    fn index(&self, index: PropagatorId) -> &Self::Output {
+        &self.constraints[index]
+    }
+}
+impl IndexMut<PropagatorId> for ConstraintDb {
+    fn index_mut(&mut self, index: PropagatorId) -> &mut Self::Output {
+        &mut self.constraints[index]
+    }
+}
+
+/// A pair of constraints (a, b) where edge(a) = !edge(b)
+struct ConstraintPair {
+    /// constraint where the edge is in its canonical form
+    base_forward: Propagator,
+    base_backward: Propagator,
+    /// constraint corresponding to the negation of base
+    negated_forward: Propagator,
+    negated_backward: Propagator,
+}
+
+impl ConstraintPair {
+    pub fn new_inactives(edge: Edge) -> ConstraintPair {
+        let edge = if edge.is_canonical() { edge } else { edge.negated() };
+        ConstraintPair {
+            base_forward: Propagator::forward(edge),
+            base_backward: Propagator::backward(edge),
+            negated_forward: Propagator::forward(edge.negated()),
+            negated_backward: Propagator::backward(edge.negated()),
+        }
+    }
+}
diff --git a/tnet/src/distances.rs b/tnet/src/theory/distances.rs
similarity index 94%
rename from tnet/src/distances.rs
rename to tnet/src/theory/distances.rs
index 44b653a2..404a104b 100644
--- a/tnet/src/distances.rs
+++ b/tnet/src/theory/distances.rs
@@ -1,6 +1,6 @@
-use crate::theory::DirEdge;
+use crate::theory::PropagatorId;
 use aries_collections::ref_store::RefMap;
-use aries_model::bounds::{BoundValueAdd, VarBound};
+use aries_core::*;
 use std::cmp::{Ordering, Reverse};
 use std::collections::BinaryHeap;
 
@@ -34,7 +34,7 @@ pub(crate) struct DijkstraState {
     latest: BoundValueAdd,
     /// Associates each vertex to its distance.
     /// If the node is not an origin, it also indicates the latest edge on the shortest path to this node.
-    pub distances: RefMap<VarBound, (BoundValueAdd, Option<DirEdge>)>,
+    pub distances: RefMap<VarBound, (BoundValueAdd, Option<PropagatorId>)>,
     /// Elements of the queue that have not been extracted yet.
     /// Note that a single node might appear several times in the queue, in which case only
     /// the one with the smallest distance is relevant.
@@ -50,7 +50,7 @@ impl DijkstraState {
 
     /// Add a node to the queue, indicating the distance from the origin and the latest edge
     /// on the path from the origin to this node.
-    pub fn enqueue(&mut self, node: VarBound, dist: BoundValueAdd, incoming_edge: Option<DirEdge>) {
+    pub fn enqueue(&mut self, node: VarBound, dist: BoundValueAdd, incoming_edge: Option<PropagatorId>) {
         let previous_dist = match self.distances.get(node) {
             None => BoundValueAdd::MAX,
             Some((prev, _)) => *prev,
@@ -98,7 +98,7 @@ impl DijkstraState {
     /// Return the predecessor edge from the origin to this node or None if it is an origin.
     ///
     /// **Panics** if the node has no associated distance (i.e. was not reached by the algorithm).
-    pub fn predecessor(&self, node: VarBound) -> Option<DirEdge> {
+    pub fn predecessor(&self, node: VarBound) -> Option<PropagatorId> {
         self.distances[node].1
     }
 
diff --git a/tnet/src/theory/edges.rs b/tnet/src/theory/edges.rs
new file mode 100644
index 00000000..380ee005
--- /dev/null
+++ b/tnet/src/theory/edges.rs
@@ -0,0 +1,240 @@
+use crate::theory::contraint_db::Enabler;
+use crate::theory::{Timepoint, W};
+use aries_core::{BoundValueAdd, Lit, VarBound};
+
+/// A unique identifier for an edge in the STN.
+/// An edge and its negation share the same `base_id` but differ by the `is_negated` property.
+///
+/// For instance, valid edge ids:
+///  -  `a - b <= 10`
+///    - base_id: 3
+///    - negated: false
+///  - `a - b > 10`       # negation of the previous one
+///    - base_id: 3        # same
+///    - negated: true     # inverse
+///  - `a - b <= 20`      # unrelated
+///    - base_id: 4
+///    - negated: false
+#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash, Debug)]
+pub struct EdgeId(u32);
+impl EdgeId {
+    #[inline]
+    pub fn new(base_id: u32, negated: bool) -> EdgeId {
+        if negated {
+            EdgeId((base_id << 1) + 1)
+        } else {
+            EdgeId(base_id << 1)
+        }
+    }
+
+    #[inline]
+    pub fn base_id(&self) -> u32 {
+        self.0 >> 1
+    }
+
+    #[inline]
+    pub fn is_negated(&self) -> bool {
+        self.0 & 0x1 == 1
+    }
+
+    /// Id of the forward (from source to target) view of this edge
+    pub(in crate::theory) fn forward(self) -> PropagatorId {
+        PropagatorId::forward(self)
+    }
+
+    /// Id of the backward view (from target to source) of this edge
+    pub(in crate::theory) fn backward(self) -> PropagatorId {
+        PropagatorId::backward(self)
+    }
+}
+
+impl std::ops::Not for EdgeId {
+    type Output = Self;
+
+    #[inline]
+    fn not(self) -> Self::Output {
+        EdgeId(self.0 ^ 0x1)
+    }
+}
+
+impl From<EdgeId> for u32 {
+    fn from(e: EdgeId) -> Self {
+        e.0
+    }
+}
+impl From<u32> for EdgeId {
+    fn from(id: u32) -> Self {
+        EdgeId(id)
+    }
+}
+
+impl From<EdgeId> for usize {
+    fn from(e: EdgeId) -> Self {
+        e.0 as usize
+    }
+}
+impl From<usize> for EdgeId {
+    fn from(id: usize) -> Self {
+        EdgeId(id as u32)
+    }
+}
+
+/// An edge in the STN, representing the constraint `target - source <= weight`
+/// An edge can be either in canonical form or in negated form.
+/// Given to edges (tgt - src <= w) and (tgt -src > w) one will be in canonical form and
+/// the other in negated form.
+#[derive(Copy, Clone, Debug, PartialOrd, Ord, PartialEq, Eq, Hash)]
+pub struct Edge {
+    pub source: Timepoint,
+    pub target: Timepoint,
+    pub weight: W,
+}
+
+impl Edge {
+    pub fn new(source: Timepoint, target: Timepoint, weight: W) -> Edge {
+        Edge { source, target, weight }
+    }
+
+    pub fn is_negated(&self) -> bool {
+        !self.is_canonical()
+    }
+
+    pub fn is_canonical(&self) -> bool {
+        self.source < self.target || self.source == self.target && self.weight >= 0
+    }
+
+    // not(b - a <= 6)
+    //   = b - a > 6
+    //   = a -b < -6
+    //   = a - b <= -7
+    //
+    // not(a - b <= -7)
+    //   = a - b > -7
+    //   = b - a < 7
+    //   = b - a <= 6
+    pub fn negated(&self) -> Self {
+        Edge {
+            source: self.target,
+            target: self.source,
+            weight: -self.weight - 1,
+        }
+    }
+}
+
+/// A `Propagator` represents the fact that an update on the `source` bound
+/// should be reflected on the `target` bound.
+///
+/// From a classical STN edge `source -- weight --> target` there will be two `Propagator`s:
+///   - ub(source) = X   implies   ub(target) <= X + weight
+///   - lb(target) = X   implies   lb(source) >= X - weight
+#[derive(Clone, Debug)]
+pub(crate) struct Propagator {
+    pub source: VarBound,
+    pub target: VarBound,
+    pub weight: BoundValueAdd,
+    /// Non-empty if the constraint active (participates in propagation)
+    /// If the enabler is Lit::TRUE, then the constraint can be assumed to be always active
+    pub enabler: Option<Enabler>,
+    /// A set of potential enablers for this constraint.
+    /// The edge becomes active once one of its enablers becomes true
+    pub enablers: Vec<Enabler>,
+}
+impl Propagator {
+    /// source <= X   =>   target <= X + weight
+    pub fn forward(edge: Edge) -> Propagator {
+        Propagator {
+            source: VarBound::ub(edge.source),
+            target: VarBound::ub(edge.target),
+            weight: BoundValueAdd::on_ub(edge.weight),
+            enabler: None,
+            enablers: vec![],
+        }
+    }
+
+    /// target >= X   =>   source >= X - weight
+    pub fn backward(edge: Edge) -> Propagator {
+        Propagator {
+            source: VarBound::lb(edge.target),
+            target: VarBound::lb(edge.source),
+            weight: BoundValueAdd::on_lb(-edge.weight),
+            enabler: None,
+            enablers: vec![],
+        }
+    }
+
+    pub fn as_edge(&self) -> Edge {
+        if self.source.is_ub() {
+            debug_assert!(self.target.is_ub());
+            Edge {
+                source: self.source.variable(),
+                target: self.target.variable(),
+                weight: self.weight.as_ub_add(),
+            }
+        } else {
+            debug_assert!(self.target.is_lb());
+            Edge {
+                source: self.target.variable(),
+                target: self.source.variable(),
+                weight: -self.weight.as_lb_add(),
+            }
+        }
+    }
+}
+
+/// Represents an edge together with a particular propagation direction:
+///  - forward (source to target)
+///  - backward (target to source)
+#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Debug)]
+pub(crate) struct PropagatorId(u32);
+
+impl PropagatorId {
+    /// Forward view of the given edge
+    pub fn forward(e: EdgeId) -> Self {
+        PropagatorId(u32::from(e) << 1)
+    }
+
+    /// Backward view of the given edge
+    pub fn backward(e: EdgeId) -> Self {
+        PropagatorId((u32::from(e) << 1) + 1)
+    }
+
+    #[allow(unused)]
+    pub fn is_forward(self) -> bool {
+        (u32::from(self) & 0x1) == 0
+    }
+
+    /// The edge underlying this projection
+    pub fn edge(self) -> EdgeId {
+        EdgeId::from(self.0 >> 1)
+    }
+}
+impl From<PropagatorId> for usize {
+    fn from(e: PropagatorId) -> Self {
+        e.0 as usize
+    }
+}
+impl From<usize> for PropagatorId {
+    fn from(u: usize) -> Self {
+        PropagatorId(u as u32)
+    }
+}
+impl From<PropagatorId> for u32 {
+    fn from(e: PropagatorId) -> Self {
+        e.0
+    }
+}
+impl From<u32> for PropagatorId {
+    fn from(u: u32) -> Self {
+        PropagatorId(u)
+    }
+}
+
+#[derive(Copy, Clone, Debug)]
+pub struct PropagatorTarget {
+    pub target: VarBound,
+    pub weight: BoundValueAdd,
+    /// Literal that is true if and only if the edge must be present in the network.
+    /// Note that handling of optional variables might allow and edge to propagate even it is not known
+    /// to be present yet.
+    pub presence: Lit,
+}


Test patch:
diff --git a/examples/smt/tests/smt.rs b/examples/smt/tests/smt.rs
index 62c4c99e..6f137ad0 100644
--- a/examples/smt/tests/smt.rs
+++ b/examples/smt/tests/smt.rs
@@ -1,13 +1,14 @@
 use aries_backtrack::Backtrack;
-use aries_model::bounds::Lit;
+use aries_core::state::OptDomain;
+use aries_core::Lit;
 use aries_model::extensions::AssignmentExt;
 use aries_model::lang::expr::*;
 use aries_model::lang::IVar;
-use aries_model::state::OptDomain;
-use aries_solver::solver::Solver;
 use aries_tnet::theory::{StnConfig, StnTheory};
+use itertools::Itertools;
 
 type Model = aries_model::Model<String>;
+type Solver = aries_solver::solver::Solver<String>;
 
 #[test]
 fn sat() {
@@ -235,7 +236,7 @@ fn optional_hierarchy() {
     // constraints.push(model.or(&scopes.iter().map(|&lit| BAtom::from(lit)).collect::<Vec<_>>()));
 
     for &sub_var in &vars {
-        model.enforce(opt_eq(i, sub_var));
+        model.enforce(eq(i, sub_var));
     }
 
     let mut solver = Solver::new(model);
@@ -308,3 +309,187 @@ fn optional_hierarchy() {
     // assert_eq!(solver.model.opt_domain_of(vars[1]), Absent);
     // assert_eq!(solver.model.opt_domain_of(vars[2]), Absent);
 }
+
+#[test]
+fn test_reification_optionals() {
+    let mut model = Model::new();
+    let pa = model.new_presence_variable(Lit::TRUE, "pa").true_lit();
+    let a = model.new_optional_ivar(0, 10, pa, "a");
+    let pb = model.new_presence_variable(Lit::TRUE, "pb").true_lit();
+    let b = model.new_optional_ivar(0, 10, pb, "b");
+
+    let a_leq_b = model.reify(leq(a, b));
+    let pab = model.presence_literal(a_leq_b.variable());
+    assert_ne!(pab, Lit::TRUE);
+    assert!(model.state.implies(pab, pa));
+    assert!(model.state.implies(pab, pb));
+
+    // defined if (!pa | pab) <=> (!pa | (pa & pb)) <=> pb
+    let or = model.reify(or([!pa, a_leq_b]));
+    let p_or = model.presence_literal(or.variable());
+    assert!(model.state.implies(p_or, pb));
+    assert!(!model.state.implies(p_or, pa));
+}
+
+/// A test represent the fact when **all** `premices` become true (regardless of the order),
+/// then **all** `expected` literal should be inferred as true.
+#[derive(Debug)]
+struct Test {
+    premices: Vec<Lit>,
+    expected: Vec<Lit>,
+}
+
+impl Test {
+    pub fn new(premices: &[Lit], expected: &[Lit]) -> Test {
+        Test {
+            premices: premices.into(),
+            expected: expected.into(),
+        }
+    }
+}
+
+/// Run all tests on the given solver.
+///
+/// For all possible order of achievement of the premices in a given test, this checks that
+///  - no expected literal becomes true before all premices are entailed.
+///  - all expected literals become true once all premices have been enforced.
+fn run_tests(solver: &mut Solver, tests: &[Test]) {
+    for (test_id, test) in tests.iter().enumerate() {
+        // we have a test, we will test the inferences under all possible realisation order of the premices
+        for premices_order in test.premices.iter().copied().permutations(test.premices.len()) {
+            // new test case, with a unique order of premices
+            let test = Test::new(&premices_order, &test.expected);
+
+            println!(
+                "=============================\n\
+                     Running test {} with premices order: {:?}",
+                test_id,
+                test.premices
+                    .iter()
+                    .map(|l| solver.model.fmt(*l).to_string())
+                    .collect::<Vec<_>>()
+            );
+
+            solver.reset();
+            solver.propagate().unwrap();
+
+            for (decision_id, &decision) in test.premices.iter().enumerate() {
+                solver.save_state();
+                solver.model.print_state();
+
+                // some premices have not been established yet, check that no inference was made
+                for &expected in &test.expected {
+                    assert!(
+                        !solver.model.entails(expected),
+                        "Inferred after only {} established premices: {}",
+                        decision_id,
+                        solver.model.fmt(expected)
+                    );
+                }
+
+                println!("\ndecision: {}\n", solver.model.fmt(decision));
+                solver.decide(decision);
+                let r = solver.propagate();
+                assert_eq!(
+                    r,
+                    Ok(()),
+                    "Failed to propagate the {}th premice: {} ",
+                    decision_id + 1,
+                    solver.model.fmt(decision),
+                );
+            }
+            solver.model.print_state();
+
+            // all premices established, check that all inferences were made
+            for &expected in &test.expected {
+                assert!(
+                    solver.model.entails(expected),
+                    "Not inferred: {}",
+                    solver.model.fmt(expected),
+                );
+            }
+        }
+    }
+}
+
+#[test]
+fn test_leq_propagation() {
+    let mut model = Model::new();
+
+    let a = model.new_ivar(0, 100, "a");
+    let b = model.new_ivar(0, 100, "b");
+    let l = model.reify(leq(a, b));
+
+    let tests = vec![
+        Test::new(&[l, a.geq(4)], &[b.geq(4)]),
+        Test::new(&[l, b.leq(3)], &[a.leq(3)]),
+        Test::new(&[!l, b.geq(4)], &[a.geq(5)]),
+        Test::new(&[!l, a.leq(4)], &[b.leq(3)]),
+        // the two test below should only work in the presence of bounds theory propagation (on by default)
+        Test::new(&[b.lt(5), a.geq(5)], &[!l]),
+        Test::new(&[a.leq(5), b.geq(5)], &[l]),
+    ];
+
+    let mut solver = Solver::new(model);
+    solver.add_theory(|tok| StnTheory::new(tok, StnConfig::default()));
+    run_tests(&mut solver, &tests);
+}
+
+#[test]
+fn test_opt_leq_propagation() {
+    let mut model = Model::new();
+    let pa = model.new_presence_variable(Lit::TRUE, "pa").true_lit();
+    let a = model.new_optional_ivar(0, 100, pa, "a");
+    let pb = model.new_presence_variable(Lit::TRUE, "pb").true_lit();
+    let b = model.new_optional_ivar(0, 100, pb, "b");
+    // the two test below should only work in the presence of bounds theory propagation (on by default)
+    let l = model.reify(leq(a, b));
+    model.shape.labels.insert(l.variable(), "l".to_string());
+    let v = model.presence_literal(l.variable());
+    model.shape.labels.insert(v.variable(), "v".to_string());
+
+    let tests = vec![
+        Test::new(&[v, l, a.geq(4)], &[b.geq(4)]),
+        Test::new(&[pa, l, a.geq(4)], &[b.geq(4)]),
+        Test::new(&[v, l, b.leq(3)], &[a.leq(3)]),
+        Test::new(&[pb, l, b.leq(3)], &[a.leq(3)]),
+        Test::new(&[v, !l, b.geq(4)], &[a.geq(5)]),
+        Test::new(&[pb, !l, b.geq(4)], &[a.geq(5)]),
+        Test::new(&[v, !l, a.leq(4)], &[b.leq(3)]),
+        Test::new(&[pa, !l, a.leq(4)], &[b.leq(3)]),
+        Test::new(&[b.lt(5), a.geq(5)], &[!l]),
+        Test::new(&[a.leq(5), b.geq(5)], &[l]),
+    ];
+
+    let mut solver = Solver::new(model);
+    solver.add_theory(|tok| StnTheory::new(tok, StnConfig::default()));
+    run_tests(&mut solver, &tests);
+}
+
+#[test]
+fn test_opt_leq_eager_propagation() {
+    let mut model = Model::new();
+    let pa = model.new_presence_variable(Lit::TRUE, "pa").true_lit();
+    let a = model.new_optional_ivar(0, 100, pa, "a");
+    let pb = model.new_presence_variable(pa, "pb").true_lit();
+    let b = model.new_optional_ivar(0, 100, pb, "b");
+    // the two test below should only work in the presence of bounds theory propagation (on by default)
+    let l = model.reify(leq(a, b));
+    model.shape.labels.insert(l.variable(), "l".to_string());
+    let v = model.presence_literal(l.variable());
+    debug_assert_eq!(v, pb);
+    // model.shape.labels.insert(v.variable(), "v".to_string());
+
+    let tests = vec![
+        // Test::new(&[l, a.geq(4)], &[b.geq(4)]),
+        Test::new(&[v, l, b.leq(3)], &[a.leq(3)]),
+        Test::new(&[v, !l, b.geq(4)], &[a.geq(5)]),
+        Test::new(&[!l, a.leq(4)], &[b.leq(3)]),
+        Test::new(&[b.lt(5), a.geq(5)], &[!l]),
+        Test::new(&[a.leq(5), b.geq(5)], &[l]),
+    ];
+
+    let mut solver = Solver::new(model);
+    solver.add_theory(|tok| StnTheory::new(tok, StnConfig::default()));
+    run_tests(&mut solver, &tests);
+}