Repository: plaans/aries (language: python)

Issue:
Add structured labeled on variables
We should allow labels﻿ in a user provided type to allow the definition of informed branching heuristics.

- [x] Theories should on require the `Domains` information and not need to create new varibles (move to an earlier decompose step)
- [x] Remove the `Label` type.



Discussion / hints:

Fix patch:
diff --git a/Cargo.lock b/Cargo.lock
index 78e2438a..c9a8c26a 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -78,6 +78,7 @@ dependencies = [
  "aries_collections",
  "aries_utils",
  "criterion",
+ "lazy_static",
  "rand",
  "streaming-iterator",
 ]
@@ -130,12 +131,10 @@ dependencies = [
 name = "aries_smt"
 version = "0.1.0"
 dependencies = [
- "anyhow",
  "aries_backtrack",
  "aries_model",
  "aries_solver",
  "aries_tnet",
- "structopt",
 ]
 
 [[package]]
diff --git a/examples/jobshop/main.rs b/examples/jobshop/main.rs
index 3de5bd47..a4756055 100644
--- a/examples/jobshop/main.rs
+++ b/examples/jobshop/main.rs
@@ -1,5 +1,29 @@
-use aries_model::extensions::{AssignmentExt, ExpressionFactoryExt};
+use aries_backtrack::{Backtrack, DecLvl};
+use aries_model::bounds::Lit;
+use aries_model::extensions::{AssignmentExt, Shaped};
+use aries_model::lang::expr::leq;
+use aries_model::lang::{IVar, VarRef};
+use aries_solver::solver::search::activity::{ActivityBrancher, Heuristic};
+use aries_solver::solver::search::{Decision, SearchControl};
+use aries_solver::solver::stats::Stats;
+use aries_tnet::theory::{StnConfig, StnTheory};
 use std::fmt::Write;
+use std::fs;
+use std::str::FromStr;
+use structopt::StructOpt;
+
+#[derive(Copy, Clone, Debug, Eq, PartialEq, Hash)]
+pub enum Var {
+    /// Variable representing the makespan (constrained to be after the end of tasks
+    Makespan,
+    /// Variable representing the start time of (job_number, task_number_in_job)
+    Start(usize, usize),
+    Prec(usize, usize, usize, usize),
+}
+
+type Model = aries_model::Model<Var>;
+type Solver = aries_solver::solver::Solver<Var>;
+type ParSolver = aries_solver::parallel_solver::ParSolver<Var>;
 
 #[derive(Clone, Debug)]
 struct JobShop {
@@ -10,12 +34,6 @@ struct JobShop {
 }
 
 impl JobShop {
-    pub fn op_id(&self, job: usize, op: usize) -> usize {
-        job * self.num_machines + op
-    }
-    pub fn tvar(&self, job: usize, op: usize) -> TVar {
-        TVar(self.op_id(job, op) + 2)
-    }
     pub fn duration(&self, job: usize, op: usize) -> i32 {
         self.times[job * self.num_machines + op]
     }
@@ -59,30 +77,6 @@ impl JobShop {
     }
 }
 
-#[derive(Copy, Clone, Debug, Ord, PartialOrd, PartialEq, Eq, Hash)]
-struct TVar(usize);
-
-impl From<TVar> for usize {
-    fn from(t: TVar) -> Self {
-        t.0
-    }
-}
-
-use aries_model::lang::IVar;
-use aries_solver::solver::Solver;
-
-use aries_backtrack::{Backtrack, DecLvl};
-use aries_model::bounds::Lit;
-use aries_model::Model;
-use aries_solver::parallel_solver::ParSolver;
-use aries_solver::solver::search::{Decision, SearchControl};
-use aries_solver::solver::stats::Stats;
-use aries_tnet::theory::{StnConfig, StnTheory};
-use std::collections::HashMap;
-use std::fs;
-use std::str::FromStr;
-use structopt::StructOpt;
-
 #[derive(Debug, StructOpt)]
 #[structopt(name = "jobshop")]
 struct Opt {
@@ -138,14 +132,14 @@ fn main() {
     let lower_bound = (opt.lower_bound).max(pb.makespan_lower_bound() as u32);
     println!("Initial lower bound: {}", lower_bound);
 
-    let (model, makespan, var_map) = encode(&pb, lower_bound, opt.upper_bound);
+    let model = encode(&pb, lower_bound, opt.upper_bound);
+    let makespan: IVar = IVar::new(model.shape.get_variable(&Var::Makespan).unwrap());
 
     let mut solver = Solver::new(model);
     solver.add_theory(|tok| StnTheory::new(tok, StnConfig::default()));
 
     let est_brancher = EstBrancher {
         pb: pb.clone(),
-        var_map: var_map.clone(),
         saved: DecLvl::ROOT,
     };
     let mut solver = get_solver(solver, opt.search, est_brancher);
@@ -154,7 +148,7 @@ fn main() {
 
     if let Some((optimum, solution)) = result {
         println!("Found optimal solution with makespan: {}", optimum);
-        assert_eq!(solution.lower_bound(makespan), optimum);
+        assert_eq!(solution.var_domain(makespan).lb, optimum);
 
         // Format the solution in resource order : each machine is given an ordered list of tasks to process.
         let mut formatted_solution = String::new();
@@ -163,9 +157,9 @@ fn main() {
             let mut tasks = Vec::new();
             for j in 0..pb.num_jobs {
                 let op = pb.op_with_machine(j, m);
-                let task = pb.tvar(j, op);
-                let start_var = var_map[&task];
-                let start_time = solution.state.bounds(start_var.into()).0;
+                let task = Var::Start(j, op);
+                let start_var = solver.get_int_var(&task).unwrap();
+                let start_time = solution.var_domain(start_var).lb;
                 tasks.push(((j, op), start_time));
             }
             // sort task by their start time
@@ -179,6 +173,7 @@ fn main() {
         println!("\n=== Solution (resource order) ===");
         print!("{}", formatted_solution);
         println!("=================================\n");
+
         if let Some(output) = &opt.output {
             // write solution to file
             std::fs::write(output, formatted_solution).unwrap();
@@ -229,26 +224,24 @@ fn parse(input: &str) -> JobShop {
     }
 }
 
-fn encode(pb: &JobShop, lower_bound: u32, upper_bound: u32) -> (Model, IVar, HashMap<TVar, IVar>) {
+fn encode(pb: &JobShop, lower_bound: u32, upper_bound: u32) -> Model {
+    let start = |model: &Model, j: usize, t: usize| IVar::new(model.shape.get_variable(&Var::Start(j, t)).unwrap());
+    let end = |model: &Model, j: usize, t: usize| start(model, j, t) + pb.duration(j, t);
+
     let lower_bound = lower_bound as i32;
     let upper_bound = upper_bound as i32;
     let mut m = Model::new();
-    let mut hmap: HashMap<TVar, IVar> = HashMap::new();
-    let mut constraints = Vec::new();
 
-    let makespan_variable = m.new_ivar(lower_bound, upper_bound, "makespan");
+    let makespan_variable = m.new_ivar(lower_bound, upper_bound, Var::Makespan);
     for j in 0..pb.num_jobs {
         for i in 0..pb.num_machines {
-            let tji = pb.tvar(j, i);
-            let task_start = m.new_ivar(0, upper_bound, format!("start({}, {})", j, i));
-            hmap.insert(tji, task_start);
+            let task_start = m.new_ivar(0, upper_bound, Var::Start(j, i));
 
             let left_on_job: i32 = (i..pb.num_machines).map(|t| pb.duration(j, t)).sum();
-            constraints.push(m.leq(task_start + left_on_job, makespan_variable));
+            m.enforce(leq(task_start + left_on_job, makespan_variable));
 
             if i > 0 {
-                let end_of_previous = hmap[&pb.tvar(j, i - 1)] + pb.duration(j, i - 1);
-                constraints.push(m.leq(end_of_previous, task_start));
+                m.enforce(leq(end(&m, j, i - 1), task_start));
             }
         }
     }
@@ -258,28 +251,40 @@ fn encode(pb: &JobShop, lower_bound: u32, upper_bound: u32) -> (Model, IVar, Has
                 let i1 = pb.op_with_machine(j1, machine);
                 let i2 = pb.op_with_machine(j2, machine);
 
-                let tji1 = hmap[&pb.tvar(j1, i1)];
-                let tji2 = hmap[&pb.tvar(j2, i2)];
-                let o1 = m.leq(tji1 + pb.duration(j1, i1), tji2);
-                let o2 = m.leq(tji2 + pb.duration(j2, i2), tji1);
-                constraints.push(m.or2(o1, o2));
+                // variable that is true if (j1, i1) comes first and false otherwise.
+                // in any case, setting a value to it enforces that the two tasks do not overlap
+                let prec = m.new_bvar(Var::Prec(j1, i1, j2, i2));
+                m.bind(leq(end(&m, j1, i1), start(&m, j2, i2)), prec.true_lit());
+                m.bind(leq(end(&m, j2, i2), start(&m, j1, i1)), prec.false_lit());
             }
         }
     }
-    m.enforce_all(&constraints);
-    (m, makespan_variable, hmap)
+    m
+}
+
+struct ResourceOrderingFirst;
+impl Heuristic<Var> for ResourceOrderingFirst {
+    fn decision_stage(&self, _var: VarRef, label: Option<&Var>, _model: &aries_model::Model<Var>) -> u8 {
+        match label {
+            Some(&Var::Prec(_, _, _, _)) => 0, // a reification of (a <= b), decide in the first stage
+            Some(&Var::Makespan) | Some(&Var::Start(_, _)) => 1, // delay decisions on the temporal variable to the second stage
+            _ => 2,
+        }
+    }
 }
 
 /// Builds a solver for the given strategy.
 fn get_solver(base: Solver, strategy: SearchStrategy, est_brancher: EstBrancher) -> ParSolver {
     let base_solver = Box::new(base);
+    let make_act = |s: &mut Solver| s.set_brancher(ActivityBrancher::new_with_heuristic(ResourceOrderingFirst));
+    let make_est = |s: &mut Solver| s.set_brancher(est_brancher.clone());
     match strategy {
-        SearchStrategy::Activity => ParSolver::new(base_solver, 1, |_, _| {}),
-        SearchStrategy::Est => ParSolver::new(base_solver, 1, |_, s| s.set_brancher(est_brancher.clone())),
-        SearchStrategy::Parallel => ParSolver::new(base_solver, 2, |id, s| {
-            if id == 1 {
-                s.set_brancher(est_brancher.clone())
-            }
+        SearchStrategy::Activity => ParSolver::new(base_solver, 1, |_, s| make_act(s)),
+        SearchStrategy::Est => ParSolver::new(base_solver, 1, |_, s| make_est(s)),
+        SearchStrategy::Parallel => ParSolver::new(base_solver, 2, |id, s| match id {
+            0 => make_act(s),
+            1 => make_est(s),
+            _ => unreachable!(),
         }),
     }
 }
@@ -287,15 +292,14 @@ fn get_solver(base: Solver, strategy: SearchStrategy, est_brancher: EstBrancher)
 #[derive(Clone)]
 struct EstBrancher {
     pb: JobShop,
-    var_map: HashMap<TVar, IVar>,
     saved: DecLvl,
 }
 
-impl SearchControl for EstBrancher {
+impl SearchControl<Var> for EstBrancher {
     fn next_decision(&mut self, _stats: &Stats, model: &Model) -> Option<Decision> {
         let active_in_job = |j: usize| {
             for t in 0..self.pb.num_machines {
-                let v = self.var_map[&self.pb.tvar(j, t)];
+                let v = model.shape.get_variable(&Var::Start(j, t)).unwrap();
                 let (lb, ub) = model.domain_of(v);
                 if lb < ub {
                     return Some((v, lb, ub));
@@ -313,7 +317,7 @@ impl SearchControl for EstBrancher {
         best.map(|(var, est, _)| Decision::SetLiteral(Lit::leq(var, est)))
     }
 
-    fn clone_to_box(&self) -> Box<dyn SearchControl + Send> {
+    fn clone_to_box(&self) -> Box<dyn SearchControl<Var> + Send> {
         Box::new(self.clone())
     }
 }
diff --git a/examples/sat/main.rs b/examples/sat/main.rs
index 22bcd41d..70650a72 100644
--- a/examples/sat/main.rs
+++ b/examples/sat/main.rs
@@ -2,9 +2,7 @@
 
 use anyhow::*;
 use aries_model::bounds::Lit;
-use aries_model::extensions::ExpressionFactoryExt;
-use aries_model::lang::BAtom;
-use aries_model::Model;
+use aries_model::lang::expr::or;
 use aries_solver::parallel_solver::ParSolver;
 use aries_solver::solver::search::activity::{ActivityBrancher, BranchingParams};
 use aries_solver::solver::Solver;
@@ -14,6 +12,8 @@ use std::io::Read;
 use std::path::{Path, PathBuf};
 use structopt::StructOpt;
 
+type Model = aries_model::Model<String>;
+
 #[derive(Debug, StructOpt)]
 #[structopt(name = "minisat")]
 struct Opt {
@@ -108,7 +108,7 @@ fn solve_multi_threads(model: Model, opt: &Opt, num_threads: usize) -> Result<()
     let search_params = search_params();
 
     let mut par_solver = ParSolver::new(solver, num_threads, |id, solver| {
-        solver.set_brancher(ActivityBrancher::with_params(search_params[id].clone()))
+        solver.set_brancher(ActivityBrancher::new_with_params(search_params[id].clone()))
     });
 
     match par_solver.solve()? {
@@ -136,9 +136,8 @@ fn solve_multi_threads(model: Model, opt: &Opt, num_threads: usize) -> Result<()
 pub fn load(cnf: varisat_formula::CnfFormula) -> Result<Model> {
     let mut var_bindings = HashMap::new();
     let mut model = Model::new();
-    let mut clauses = Vec::new();
 
-    let mut lits: Vec<BAtom> = Vec::new();
+    let mut lits: Vec<Lit> = Vec::new();
     for clause in cnf.iter() {
         lits.clear();
         for &lit in clause {
@@ -151,11 +150,11 @@ pub fn load(cnf: varisat_formula::CnfFormula) -> Result<Model> {
                 model_var
             };
             let lit: Lit = if lit.is_positive() { var.into() } else { !var };
-            lits.push(lit.into());
+            lits.push(lit);
         }
-        clauses.push(model.or(&lits));
+        model.enforce(or(lits.as_slice()));
     }
-    model.enforce_all(&clauses);
+
     Ok(model)
 }
 
diff --git a/examples/smt/Cargo.toml b/examples/smt/Cargo.toml
index 546fb70b..795d8e7f 100644
--- a/examples/smt/Cargo.toml
+++ b/examples/smt/Cargo.toml
@@ -5,8 +5,6 @@ authors = ["Arthur Bit-Monnot <abitmonnot@laas.fr>"]
 edition = "2018"
 
 [dependencies]
-structopt = "0.3"
-anyhow = "1.0"
 aries_backtrack = { path = "../../backtrack" }
 aries_model = { path = "../../model" }
 aries_solver = { path = "../../solver" }
diff --git a/model/Cargo.toml b/model/Cargo.toml
index c140f2c3..aa755333 100644
--- a/model/Cargo.toml
+++ b/model/Cargo.toml
@@ -12,6 +12,7 @@ streaming-iterator = "0.1.5"
 aries_backtrack = { path = "../backtrack" }
 aries_collections = { path = "../collections" }
 aries_utils = { path = "../utils" }
+lazy_static = "1.4.0"
 
 
 [dev-dependencies]
diff --git a/model/src/bindings.rs b/model/src/bindings.rs
deleted file mode 100644
index d1af08f0..00000000
--- a/model/src/bindings.rs
+++ /dev/null
@@ -1,52 +0,0 @@
-use crate::bounds::Lit;
-use crate::expressions::ExprHandle;
-use aries_collections::ref_store::RefMap;
-
-#[derive(Copy, Clone, Debug)]
-pub enum BindTarget {
-    Expr(ExprHandle),
-    Literal(Lit),
-}
-
-#[derive(Copy, Clone)]
-pub struct BindingCursor(usize);
-
-impl BindingCursor {
-    pub fn first() -> Self {
-        BindingCursor(0)
-    }
-}
-
-/// A structure to keep track of bindings between a literal and an expression or between to literals.
-#[derive(Clone, Default)]
-pub struct Bindings {
-    /// If the expression was reified, associates the handle with the equivalent literal
-    bindings: RefMap<ExprHandle, Lit>,
-    /// All binding events in chronological order. This is intended to easily process
-    /// binding events and detect whether new events have been added.
-    binding_events: Vec<(Lit, BindTarget)>,
-}
-
-impl Bindings {
-    pub fn as_lit(&self, expr_id: ExprHandle) -> Option<Lit> {
-        self.bindings.get(expr_id).copied()
-    }
-
-    pub fn bind(&mut self, expr_id: ExprHandle, literal: Lit) {
-        assert!(!self.bindings.contains(expr_id));
-        self.bindings.insert(expr_id, literal);
-        self.binding_events.push((literal, BindTarget::Expr(expr_id)));
-    }
-
-    pub fn bind_literals(&mut self, l1: Lit, l2: Lit) {
-        self.binding_events.push((l1, BindTarget::Literal(l2)))
-    }
-
-    pub fn pop_next_event(&self, cursor: &mut BindingCursor) -> Option<(Lit, BindTarget)> {
-        let ret = self.binding_events.get(cursor.0).copied();
-        if ret.is_some() {
-            cursor.0 += 1;
-        }
-        ret
-    }
-}
diff --git a/model/src/bounds/bound.rs b/model/src/bounds/bound.rs
index 94490d1a..4c64adf6 100644
--- a/model/src/bounds/bound.rs
+++ b/model/src/bounds/bound.rs
@@ -1,9 +1,10 @@
 use crate::bounds::var_bound::VarBound;
 use crate::bounds::BoundValue;
-use crate::lang::BVar;
+use crate::lang::{BVar, ConversionError};
 use crate::lang::{IntCst, VarRef};
 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
@@ -21,7 +22,7 @@ use std::mem::transmute;
 /// use aries_model::Model;
 /// use aries_model::lang::VarRef;
 /// use aries_model::bounds::{Lit, Relation};
-/// let mut model = Model::new();
+/// 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();
@@ -51,7 +52,7 @@ use std::mem::transmute;
 /// ```
 /// use aries_model::Model;
 /// use aries_model::bounds::Lit;
-/// let mut model = Model::new();
+/// 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)];
@@ -229,6 +230,28 @@ impl From<BVar> for Lit {
     }
 }
 
+impl From<bool> for Lit {
+    fn from(b: bool) -> Self {
+        if b {
+            Lit::TRUE
+        } else {
+            Lit::FALSE
+        }
+    }
+}
+
+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/model/src/bounds/disjunction.rs b/model/src/bounds/disjunction.rs
index 1991fad8..130a3545 100644
--- a/model/src/bounds/disjunction.rs
+++ b/model/src/bounds/disjunction.rs
@@ -4,7 +4,7 @@ use std::borrow::Borrow;
 /// A set of literals representing a disjunction.
 /// A `Disjunction` maintains the invariant that there are not duplicated literals (a literal that entails another one).
 /// Implementation maintains the literals sorted.
-#[derive(Debug, PartialEq, Clone)]
+#[derive(Debug, PartialEq, Clone, Eq, Hash)]
 pub struct Disjunction {
     literals: Vec<Lit>,
 }
diff --git a/model/src/bounds/watches.rs b/model/src/bounds/watches.rs
index 956b4a28..2ff40465 100644
--- a/model/src/bounds/watches.rs
+++ b/model/src/bounds/watches.rs
@@ -174,7 +174,7 @@ mod tests {
 
     #[test]
     fn test_watches() {
-        let mut model = Model::new();
+        let mut model: Model<&'static str> = Model::new();
         let a = model.new_ivar(0, 10, "a");
         let b = model.new_ivar(0, 10, "b");
 
diff --git a/model/src/decomposition.rs b/model/src/decomposition.rs
new file mode 100644
index 00000000..e098508e
--- /dev/null
+++ b/model/src/decomposition.rs
@@ -0,0 +1,116 @@
+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::reification::{downcast, BindTarget, BindingCursor, Expr};
+use crate::lang::IVar;
+use crate::Model;
+use std::sync::Arc;
+
+/// Module to constructs the constraints from the model and store them in a queue.
+#[derive(Clone)]
+pub struct Constraints<Lbl> {
+    /// Object in charge of transforming high-level constraints into simpler ones that the solver can handle
+    decomposer: Arc<dyn Decompose<Lbl>>,
+    /// A queue of constraints that has been decomposed from the model's bindings.
+    constraints: Vec<(Lit, BindTarget)>,
+    /// Cursor into the model's bindings
+    next_to_decompose: BindingCursor,
+    /// Cursor into the constraints queue, that can be pulled by a solver to import constraints
+    next_constraint: usize,
+}
+
+impl<Lbl> Constraints<Lbl> {
+    pub fn new<T: Decompose<Lbl> + 'static>(decomposer: T) -> Self {
+        Constraints {
+            decomposer: Arc::new(decomposer),
+            constraints: vec![],
+            next_to_decompose: BindingCursor::first(),
+            next_constraint: Default::default(),
+        }
+    }
+
+    /// Pulls all bindings in the model, apply decomposition and load maximally decomposed constraints into our queue.
+    ///
+    /// Not that in the process, the subconstraints resulting from the decomposition are added to the model.
+    pub fn decompose_all(&mut self, model: &mut Model<Lbl>) {
+        while let Some((lit, target)) = model
+            .shape
+            .expressions
+            .pop_next_event(&mut self.next_to_decompose)
+            .cloned()
+        {
+            match target {
+                BindTarget::Expr(e) => {
+                    match self.decomposer.decompose(lit, e.as_ref(), model) {
+                        DecompositionResult::Decomposed => {
+                            // constraint was decomposed and subconstraints were added to the model's expressions.
+                            // do nothing as we will encounter them later in the loop
+                        }
+                        DecompositionResult::Inapplicable => {
+                            // could not be further decomposed, add it as is to the constraints
+                            self.constraints.push((lit, BindTarget::Expr(e)));
+                        }
+                    }
+                }
+                BindTarget::Literal(l2) => {
+                    // nothing to decompose, transfer it immediately
+                    self.constraints.push((lit, BindTarget::Literal(l2)))
+                }
+            }
+        }
+    }
+
+    /// Removes and returns the next constraint in the queue.
+    pub fn pop_next_constraint(&mut self) -> Option<&(Lit, BindTarget)> {
+        let ret = self.constraints.get(self.next_constraint);
+        if ret.is_some() {
+            self.next_constraint += 1;
+        }
+        ret
+    }
+}
+
+impl<Lbl: Label> Default for Constraints<Lbl> {
+    fn default() -> Self {
+        Self::new(Eq2Leq)
+    }
+}
+
+pub enum DecompositionResult {
+    Decomposed,
+    Inapplicable,
+}
+
+pub trait Decompose<Lbl>: Send + Sync {
+    fn decompose(&self, binding: Lit, expression: &Expr, model: &mut Model<Lbl>) -> DecompositionResult;
+}
+
+pub struct Eq2Leq;
+
+impl<Lbl: Label> Decompose<Lbl> for Eq2Leq {
+    fn decompose(&self, literal: Lit, expr: &Expr, model: &mut Model<Lbl>) -> DecompositionResult {
+        if let Some(&NFEq { lhs, rhs, rhs_add }) = downcast(expr) {
+            // 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
+            } else {
+                DecompositionResult::Inapplicable
+            }
+        } else {
+            DecompositionResult::Inapplicable
+        }
+    }
+}
diff --git a/model/src/expressions.rs b/model/src/expressions.rs
deleted file mode 100644
index 068fb3be..00000000
--- a/model/src/expressions.rs
+++ /dev/null
@@ -1,63 +0,0 @@
-use crate::lang::{BExpr, Expr};
-use aries_collections::ref_store::RefVec;
-use std::collections::HashMap;
-use std::sync::Arc;
-
-/// Datastructure that interns expressions and gives them a handle that can be used to retrieve the full expression.
-///
-/// Boolean expressions can optionally be associated with a literal that is constrained to have the same value.
-#[derive(Default, Clone)]
-pub struct Expressions {
-    interned: HashMap<Arc<Expr>, ExprHandle>,
-    expressions: RefVec<ExprHandle, Arc<Expr>>,
-}
-#[derive(Eq, PartialEq)]
-pub enum NExpr {
-    Pos(Arc<Expr>),
-    Neg(Arc<Expr>),
-}
-
-// Identifier of an expression which can be retrieved with [Expressions::get]
-aries_collections::create_ref_type!(ExprHandle);
-
-impl Expressions {
-    pub fn is_interned(&self, expr: &Expr) -> bool {
-        self.interned.contains_key(expr)
-    }
-
-    pub fn handle_of(&self, expr: &Expr) -> Option<ExprHandle> {
-        self.interned.get(expr).copied()
-    }
-
-    pub fn get(&self, expr_id: ExprHandle) -> Arc<Expr> {
-        self.expressions[expr_id].clone()
-    }
-
-    pub fn get_ref(&self, expr_id: ExprHandle) -> &Expr {
-        &self.expressions[expr_id]
-    }
-
-    /// Interns the given expression and returns the corresponding handle.
-    /// If the expression was already interned, the handle to the previously inserted
-    /// instance will be returned.
-    pub fn intern(&mut self, expr: Expr) -> ExprHandle {
-        if let Some(handle) = self.interned.get(&expr) {
-            *handle
-        } else {
-            let expr = Arc::new(expr);
-            let handle = self.expressions.push(expr.clone());
-            self.interned.insert(expr, handle);
-            handle
-        }
-    }
-
-    pub fn expr_of(&self, atom: impl Into<BExpr>) -> NExpr {
-        let atom = atom.into();
-        let e = self.get(atom.expr);
-        if atom.negated {
-            NExpr::Neg(e)
-        } else {
-            NExpr::Pos(e)
-        }
-    }
-}
diff --git a/model/src/extensions/assignments.rs b/model/src/extensions/assignments.rs
index 7a8fbb1f..611d9ccc 100644
--- a/model/src/extensions/assignments.rs
+++ b/model/src/extensions/assignments.rs
@@ -1,18 +1,16 @@
 use crate::bounds::Lit;
 use crate::extensions::SavedAssignment;
-use crate::lang::{Atom, BAtom, BExpr, IAtom, IVar, IntCst, SAtom, VarRef};
+use crate::lang::{Atom, IAtom, IVar, IntCst, SAtom, VarRef};
 use crate::state::{IntDomain, OptDomain};
+use crate::symbols::ContiguousSymbols;
 use crate::symbols::SymId;
-use crate::symbols::{ContiguousSymbols, SymbolTable};
 
+/// Extension methods for an object containing a partial or total assignment to a problem.
 pub trait AssignmentExt {
-    fn symbols(&self) -> &SymbolTable;
-
     fn entails(&self, literal: Lit) -> bool;
 
-    fn literal_of_expr(&self, expr: BExpr) -> Option<Lit>;
+    fn var_domain(&self, var: impl Into<IAtom>) -> IntDomain;
 
-    fn var_domain(&self, var: impl Into<VarRef>) -> IntDomain;
     fn presence_literal(&self, variable: VarRef) -> Lit;
 
     fn value_of_literal(&self, literal: Lit) -> Option<bool> {
@@ -38,19 +36,15 @@ pub trait AssignmentExt {
 
     fn domain_of(&self, atom: impl Into<IAtom>) -> (IntCst, IntCst) {
         let atom = atom.into();
-        let base = atom
-            .var
-            .map(|v| self.var_domain(v))
-            .unwrap_or_else(|| IntDomain::new(0, 0));
+        let base = self.var_domain(atom.var);
         (base.lb + atom.shift, base.ub + atom.shift)
     }
 
     /// Returns the domain of an optional integer expression.
     fn opt_domain_of(&self, atom: impl Into<IAtom>) -> OptDomain {
         let atom = atom.into();
-        let var = atom.var.map(VarRef::from).unwrap_or(VarRef::ZERO);
         let (lb, ub) = self.domain_of(atom);
-        let prez = self.presence_literal(var);
+        let prez = self.presence_literal(atom.var.into());
         match self.value_of_literal(prez) {
             Some(true) => OptDomain::Present(lb, ub),
             Some(false) => OptDomain::Absent,
@@ -84,13 +78,8 @@ pub trait AssignmentExt {
     /// Return None otherwise meaning the value con be
     ///  - either true or false
     ///  - neither true nor false (empty domain)
-    fn boolean_value_of(&self, batom: impl Into<BAtom>) -> Option<bool> {
-        let batom = batom.into();
-        match batom {
-            BAtom::Cst(value) => Some(value),
-            BAtom::Literal(b) => self.value_of_literal(b),
-            BAtom::Expr(e) => self.literal_of_expr(e).and_then(|l| self.value_of_literal(l)),
-        }
+    fn boolean_value_of(&self, bool_atom: impl Into<Lit>) -> Option<bool> {
+        self.value_of_literal(bool_atom.into())
     }
 
     /// Return an integer view of the domain of any kind of atom.
diff --git a/model/src/extensions/expression_factory.rs b/model/src/extensions/expression_factory.rs
deleted file mode 100644
index 08218e91..00000000
--- a/model/src/extensions/expression_factory.rs
+++ /dev/null
@@ -1,281 +0,0 @@
-use crate::bounds::Lit;
-use crate::lang::{Atom, BAtom, BExpr, Expr, Fun, IAtom, SAtom, VarRef};
-use std::cmp::Ordering;
-
-/// Provides extension methods to allow the construction of expressions.
-///
-/// The implementer of this trait only requires the capability of interning expressions.
-pub trait ExpressionFactoryExt {
-    fn intern_bool(&mut self, expr: Expr) -> BExpr;
-    fn presence_literal(&self, variable: VarRef) -> Lit;
-
-    // ======= Convenience methods to create expressions ========
-
-    fn or(&mut self, disjuncts: &[BAtom]) -> BAtom {
-        self.or_from_iter(disjuncts.iter().copied())
-    }
-
-    fn or_from_iter(&mut self, disjuncts: impl IntoIterator<Item = BAtom>) -> BAtom {
-        let mut or: Vec<BAtom> = disjuncts.into_iter().collect();
-        or.sort_by(BAtom::lexical_cmp);
-        or.dedup();
-        let e = Expr::new(Fun::Or, or.iter().copied().map(Atom::from).collect());
-        self.intern_bool(e).into()
-    }
-
-    fn and(&mut self, conjuncts: &[BAtom]) -> BAtom {
-        self.and_from_iter(conjuncts.iter().copied())
-    }
-
-    fn and_from_iter(&mut self, conjuncts: impl Iterator<Item = BAtom>) -> BAtom {
-        !self.or_from_iter(conjuncts.map(|b| !b))
-    }
-
-    fn and2(&mut self, a: BAtom, b: BAtom) -> BAtom {
-        self.and(&[a, b])
-    }
-    fn or2(&mut self, a: BAtom, b: BAtom) -> BAtom {
-        let and = Expr::new2(Fun::Or, a, b);
-        self.intern_bool(and).into()
-    }
-
-    fn leq<A: Into<IAtom>, B: Into<IAtom>>(&mut self, a: A, b: B) -> BAtom {
-        let mut a = a.into();
-        let mut b = b.into();
-
-        // normalize, transfer the shift from right to left
-        a.shift -= b.shift;
-        b.shift = 0;
-
-        let x = a.shift;
-        // we are in the form va + X <= vb
-
-        // 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
-        match (a.var, b.var) {
-            (None, None) => {
-                // X <= 0
-                return BAtom::Cst(x <= 0);
-            }
-            (Some(va), Some(vb)) if va == vb => {
-                // va +X <= va   <=>  X <= 0
-                return BAtom::Cst(x <= 0);
-            }
-            (Some(va), None) => {
-                // va + X <= 0   <=> va <= -X
-                return Lit::leq(va, -x).into();
-            }
-            (None, Some(vb)) => {
-                // X <= vb   <=>  vb >= X
-                return Lit::geq(vb, x).into();
-            }
-            (_, _) => {
-                // general, form, continue
-            }
-        }
-
-        // maintain the invariant that left side of the LEQ has a small lexical order
-        match a.lexical_cmp(&b) {
-            Ordering::Less => {
-                let leq = Expr::new2(Fun::Leq, a, b);
-                self.intern_bool(leq).into()
-            }
-            Ordering::Equal => true.into(),
-            Ordering::Greater => {
-                // swap the order by making !(b + 1 <= a)
-                // normalize, transfer the shift from right to left
-                b.shift -= a.shift;
-                a.shift = 0;
-
-                let leq = Expr::new2(Fun::Leq, b + 1, a);
-                (!self.intern_bool(leq)).into()
-            }
-        }
-    }
-
-    fn geq<A: Into<IAtom>, B: Into<IAtom>>(&mut self, a: A, b: B) -> BAtom {
-        self.leq(b, a)
-    }
-
-    fn lt<A: Into<IAtom>, B: Into<IAtom>>(&mut self, a: A, b: B) -> BAtom {
-        let a = a.into();
-        let b = b.into();
-        self.leq(a + 1, b)
-    }
-
-    fn gt<A: Into<IAtom>, B: Into<IAtom>>(&mut self, a: A, b: B) -> BAtom {
-        self.lt(b, a)
-    }
-
-    fn eq<A: Into<Atom>, B: Into<Atom>>(&mut self, a: A, b: B) -> BAtom {
-        let a = a.into();
-        let b = b.into();
-        if a == b {
-            BAtom::Cst(true)
-        } else if a.kind() != b.kind() {
-            BAtom::Cst(false)
-        } else {
-            use Atom::*;
-            match (a, b) {
-                (Bool(_a), Bool(_b)) => todo!(),
-                (Int(a), Int(b)) => self.int_eq(a, b),
-                (Sym(a), Sym(b)) => self.sym_eq(a, b),
-                _ => unreachable!(), // guarded by kind comparison
-            }
-        }
-    }
-
-    fn int_eq<A: Into<IAtom>, B: Into<IAtom>>(&mut self, a: A, b: B) -> BAtom {
-        let mut a = a.into();
-        let mut b = b.into();
-
-        // normalize, transfer the shift from right to left
-        a.shift -= b.shift;
-        b.shift = 0;
-
-        match a.lexical_cmp(&b) {
-            Ordering::Less => {
-                let eq = Expr::new2(Fun::Eq, a, b);
-                self.intern_bool(eq).into()
-            }
-            Ordering::Equal => true.into(),
-            Ordering::Greater => {
-                // normalize, transfer the shift from right to left
-                b.shift -= a.shift;
-                a.shift = 0;
-                let eq = Expr::new2(Fun::Eq, b, a);
-                self.intern_bool(eq).into()
-            }
-        }
-    }
-
-    fn sym_eq<A: Into<SAtom>, B: Into<SAtom>>(&mut self, a: A, b: B) -> BAtom {
-        self.int_eq(a.into().int_view(), b.into().int_view())
-    }
-
-    fn neq<A: Into<Atom>, B: Into<Atom>>(&mut self, a: A, b: B) -> BAtom {
-        !self.eq(a, b)
-    }
-
-    fn implies<A: Into<BAtom>, B: Into<BAtom>>(&mut self, a: A, b: B) -> BAtom {
-        let a = a.into();
-        let b = b.into();
-        let implication = Expr::new2(Fun::Or, !a, b);
-        self.intern_bool(implication).into()
-    }
-
-    // =========== Optionals ===============
-
-    /// Specifies that two optional variables must be equal if present.
-    fn opt_eq<A: Into<Atom>, B: Into<Atom>>(&mut self, a: A, b: B) -> BAtom {
-        let a = a.into();
-        let b = b.into();
-        if a == b {
-            BAtom::Cst(true)
-        } else if a.kind() != b.kind() {
-            BAtom::Cst(false)
-        } else {
-            use Atom::*;
-            match (a, b) {
-                (Bool(_a), Bool(_b)) => todo!(),
-                (Int(a), Int(b)) => self.opt_int_eq(a, b),
-                (Sym(a), Sym(b)) => self.opt_sym_eq(a, b),
-                _ => unreachable!(), // guarded by kind comparison
-            }
-        }
-    }
-
-    /// Specifies that two optional int variables must be equal if present.
-    fn opt_int_eq<A: Into<IAtom>, B: Into<IAtom>>(&mut self, a: A, b: B) -> BAtom {
-        let mut a = a.into();
-        let mut b = b.into();
-
-        // normalize, transfer the shift from right to left
-        a.shift -= b.shift;
-        b.shift = 0;
-
-        match a.lexical_cmp(&b) {
-            Ordering::Less => {
-                let eq = Expr::new2(Fun::OptEq, a, b);
-                self.intern_bool(eq).into()
-            }
-            Ordering::Equal => true.into(),
-            Ordering::Greater => {
-                // normalize, transfer the shift from right to left
-                b.shift -= a.shift;
-                a.shift = 0;
-                let eq = Expr::new2(Fun::OptEq, b, a);
-                self.intern_bool(eq).into()
-            }
-        }
-    }
-
-    /// Specifies that two optional symbolic variables must be equal if present.
-    fn opt_sym_eq<A: Into<SAtom>, B: Into<SAtom>>(&mut self, a: A, b: B) -> BAtom {
-        self.opt_int_eq(a.into().int_view(), b.into().int_view())
-    }
-
-    /// Specifies that, if the two variables are present, then a <= b
-    fn opt_leq<A: Into<IAtom>, B: Into<IAtom>>(&mut self, a: A, b: B) -> BAtom {
-        let mut a = a.into();
-        let mut b = b.into();
-
-        // normalize, transfer the shift from right to left
-        a.shift -= b.shift;
-        b.shift = 0;
-
-        let x = a.shift;
-        // we are in the form va + X <= vb
-
-        // 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
-        match (a.var, b.var) {
-            (None, None) => {
-                // X <= 0
-                return BAtom::Cst(x <= 0);
-            }
-            (Some(va), Some(vb)) if va == vb => {
-                // va +X <= va   <=>  X <= 0
-                return if x <= 0 {
-                    // this is always true
-                    BAtom::Cst(true)
-                } else {
-                    // base expression is always violated, so only valid
-                    // if the variable is absent
-                    (!self.presence_literal(va.into())).into()
-                };
-            }
-            (Some(va), None) => {
-                // va + X <= 0   <=> va <= -X
-                return Lit::leq(va, -x).into();
-            }
-            (None, Some(vb)) => {
-                // X <= vb   <=>  vb >= X
-                return Lit::geq(vb, x).into();
-            }
-            (_, _) => {
-                // general, form, continue
-            }
-        }
-
-        // maintain the invariant that left side of the LEQ has a small lexical order
-        match a.lexical_cmp(&b) {
-            Ordering::Less => {
-                let leq = Expr::new2(Fun::OptLeq, a, b);
-                self.intern_bool(leq).into()
-            }
-            Ordering::Equal => true.into(),
-            Ordering::Greater => {
-                // swap the order by making !(b + 1 <= a)
-                // normalize, transfer the shift from right to left
-                b.shift -= a.shift;
-                a.shift = 0;
-
-                let leq = Expr::new2(Fun::OptLeq, b + 1, a);
-                (!self.intern_bool(leq)).into()
-            }
-        }
-    }
-}
diff --git a/model/src/extensions/format.rs b/model/src/extensions/format.rs
index 5e260667..c6a30dbe 100644
--- a/model/src/extensions/format.rs
+++ b/model/src/extensions/format.rs
@@ -1,16 +1,28 @@
-use crate::expressions::ExprHandle;
-use crate::lang::{Atom, BAtom, BExpr, Expr, IAtom, Kind, SAtom, VarRef};
+use crate::bounds::Lit;
+use crate::label::Label;
+use crate::lang::{Atom, IAtom, IVar, Kind, SAtom, VarRef};
 use crate::symbols::{SymId, SymbolTable};
 use crate::types::TypeId;
 use crate::ModelShape;
 use aries_utils::input::Sym;
 use aries_utils::Fmt;
 
-pub trait Shaped {
-    fn get_shape(&self) -> &ModelShape;
+pub trait Shaped<Lbl>
+where
+    Lbl: Label,
+{
+    fn get_shape(&self) -> &ModelShape<Lbl>;
 
-    fn get_label(&self, var: impl Into<VarRef>) -> Option<&str> {
-        self.get_shape().labels.get(var.into()).map(|s| s.as_str())
+    fn get_label(&self, var: impl Into<VarRef>) -> Option<&Lbl> {
+        self.get_shape().labels.get(var.into())
+    }
+
+    fn get_var(&self, label: &Lbl) -> Option<VarRef> {
+        self.get_shape().get_variable(label)
+    }
+
+    fn get_int_var(&self, label: &Lbl) -> Option<IVar> {
+        self.get_var(label).map(IVar::new)
     }
 
     fn get_symbol(&self, sym: SymId) -> &Sym {
@@ -21,10 +33,6 @@ pub trait Shaped {
         self.get_shape().symbols.type_of(sym)
     }
 
-    fn get_expr(&self, expr: ExprHandle) -> &Expr {
-        self.get_shape().expressions.get_ref(expr)
-    }
-
     fn get_symbol_table(&self) -> &SymbolTable {
         &self.get_shape().symbols
     }
@@ -40,59 +48,34 @@ pub trait Shaped {
 /// ```
 /// use aries_model::Model;
 /// use aries_model::extensions::fmt;
-/// let mut i = Model::default();
+/// let mut i = Model::<&'static str>::default();
 /// let x = i.new_ivar(0, 10, "X");
 /// let y = x + 10;
 /// println!("x: {}", fmt(x, &i));
 /// println!("y: {}", fmt(y, &i));
 /// ```
-pub fn fmt(atom: impl Into<Atom>, ctx: &impl Shaped) -> impl std::fmt::Display + '_ {
+pub fn fmt<Lbl: Label>(atom: impl Into<Atom>, ctx: &impl Shaped<Lbl>) -> impl std::fmt::Display + '_ {
     let atom = atom.into();
     Fmt(move |f| format_impl(ctx, atom, f))
 }
 
-fn format_impl(ctx: &impl Shaped, atom: Atom, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+fn format_impl<Lbl: Label>(ctx: &impl Shaped<Lbl>, atom: Atom, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
     match atom {
         Atom::Bool(b) => format_impl_bool(ctx, b, f),
         Atom::Int(i) => format_impl_int(ctx, i, f),
         Atom::Sym(s) => format_impl_sym(ctx, s, f),
     }
 }
-fn format_impl_bool(ctx: &impl Shaped, atom: BAtom, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-    match atom {
-        BAtom::Cst(b) => write!(f, "{}", b),
-        BAtom::Literal(b) => {
-            format_impl_var(ctx, b.variable(), Kind::Int, f)?;
-            write!(f, " {} {}", b.relation(), b.value())
-        }
-        BAtom::Expr(BExpr { expr, negated }) => {
-            if negated {
-                write!(f, "(not ")?;
-            }
-            format_impl_expr(ctx, expr, f)?;
-            if negated {
-                write!(f, ")")?;
-            }
-            Ok(())
-        }
-    }
-}
-
-fn format_impl_expr(ctx: &impl Shaped, expr: ExprHandle, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-    let expr = ctx.get_expr(expr);
-    write!(f, "({}", expr.fun)?;
-    for arg in &expr.args {
-        write!(f, " ")?;
-        format_impl(ctx, *arg, f)?;
-    }
-    write!(f, ")")
+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())
 }
 
 #[allow(clippy::comparison_chain)]
-fn format_impl_int(ctx: &impl Shaped, i: IAtom, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+fn format_impl_int<Lbl: Label>(ctx: &impl Shaped<Lbl>, i: IAtom, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
     match i.var {
-        None => write!(f, "{}", i.shift),
-        Some(v) => {
+        IVar::ZERO => write!(f, "{}", i.shift),
+        v => {
             if i.shift > 0 {
                 write!(f, "(+ ")?;
             } else if i.shift < 0 {
@@ -107,16 +90,25 @@ fn format_impl_int(ctx: &impl Shaped, i: IAtom, f: &mut std::fmt::Formatter<'_>)
     }
 }
 
-fn format_impl_sym(ctx: &impl Shaped, atom: SAtom, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+fn format_impl_sym<Lbl: Label>(
+    ctx: &impl Shaped<Lbl>,
+    atom: SAtom,
+    f: &mut std::fmt::Formatter<'_>,
+) -> std::fmt::Result {
     match atom {
         SAtom::Var(v) => format_impl_var(ctx, v.var, Kind::Sym, f),
         SAtom::Cst(s) => write!(f, "{}", ctx.get_symbol(s.sym)),
     }
 }
 
-fn format_impl_var(ctx: &impl Shaped, v: VarRef, kind: Kind, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+fn format_impl_var<Lbl: Label>(
+    ctx: &impl Shaped<Lbl>,
+    v: VarRef,
+    kind: Kind,
+    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 78816877..9410564d 100644
--- a/model/src/extensions/mod.rs
+++ b/model/src/extensions/mod.rs
@@ -3,21 +3,19 @@
 //!
 //! - [DisjunctionExt] allows querying the value of a disjunction,
 //! whether it is currently unit, ...
-//! - [ExpressionFactoryExt] provides method to create expressions in a given [Model]
 //! - [AssignmentExt] provides methods to query the value of expressions.
 
 mod assignments;
 mod disjunction;
-mod expression_factory;
 mod format;
 
 pub use assignments::*;
 pub use disjunction::*;
-pub use expression_factory::*;
 pub use format::*;
 
 use crate::bounds::Lit;
-use crate::state::Domains;
+use crate::lang::{IAtom, VarRef};
+use crate::state::{Domains, IntDomain};
 use crate::Model;
 
 pub trait PartialAssignment {
@@ -39,12 +37,33 @@ impl PartialAssignment for Domains {
     }
 }
 
-// TODO: this is correct but wasteful
-//       also, it should be moved to state
-pub type SavedAssignment = Model;
+pub type SavedAssignment = Domains;
 
 impl SavedAssignment {
-    pub fn from_model(model: &Model) -> SavedAssignment {
-        model.clone()
+    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)
+    }
+
+    fn var_domain(&self, int: impl Into<IAtom>) -> IntDomain {
+        let int = int.into();
+        let (lb, ub) = self.bounds(int.var.into());
+        IntDomain {
+            lb: lb + int.shift,
+            ub: ub + int.shift,
+        }
+    }
+
+    fn presence_literal(&self, variable: VarRef) -> Lit {
+        self.presence(variable)
+    }
+
+    fn to_owned_assignment(&self) -> SavedAssignment {
+        todo!()
     }
 }
diff --git a/model/src/label.rs b/model/src/label.rs
index 430acae5..9e14dc37 100644
--- a/model/src/label.rs
+++ b/model/src/label.rs
@@ -1,91 +1,54 @@
-use aries_utils::input::Sym;
+use crate::lang::VarRef;
+use aries_collections::ref_store::RefMap;
+use std::collections::HashMap;
+use std::fmt::Debug;
+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 {}
 
-/// An optional user facing label to an object in the model.
-/// In essence this is just an `Option<String>` with the string guaranteed to be non empty.
-/// The added value comes from the various automatic conversions common label types.
-/// TODO: refactor
 #[derive(Clone)]
-pub struct Label {
-    pub(crate) lbl: Option<String>,
+pub struct VariableLabels<Lbl> {
+    labels: RefMap<VarRef, Arc<Lbl>>,
+    labeled_variables: HashMap<Arc<Lbl>, Vec<VarRef>>,
 }
-impl Label {
-    pub fn new(str: String) -> Label {
-        Label { lbl: Some(str) }
-    }
-
-    pub fn empty() -> Label {
-        Label { lbl: None }
-    }
 
-    pub fn get(&self) -> Option<&str> {
-        self.lbl.as_deref()
-    }
-}
-impl From<Label> for Option<String> {
-    fn from(lbl: Label) -> Self {
-        lbl.lbl
-    }
-}
-impl From<String> for Label {
-    fn from(str: String) -> Self {
-        if str.is_empty() {
-            Label::empty()
-        } else {
-            Label::new(str)
+impl<Lbl> VariableLabels<Lbl> {
+    pub fn new() -> Self {
+        Self {
+            labels: Default::default(),
+            labeled_variables: Default::default(),
         }
     }
-}
-impl From<&String> for Label {
-    fn from(str: &String) -> Self {
-        if str.is_empty() {
-            Label::empty()
-        } else {
-            Label::new(str.clone())
-        }
+    pub fn get(&self, var: VarRef) -> Option<&Lbl> {
+        self.labels.get(var).map(|l| l.as_ref())
     }
-}
-impl From<&str> for Label {
-    fn from(str: &str) -> Self {
-        if str.is_empty() {
-            Label::empty()
-        } else {
-            Label::new(str.into())
-        }
-    }
-}
 
-impl From<Sym> for Label {
-    fn from(sym: Sym) -> Self {
-        Label::from(sym.as_str())
-    }
-}
-impl From<&Sym> for Label {
-    fn from(sym: &Sym) -> Self {
-        Label::from(sym.as_str())
+    pub fn insert(&mut self, var: VarRef, label: impl Into<Arc<Lbl>>)
+    where
+        Lbl: Label,
+    {
+        let label = label.into();
+        self.labels.insert(var, label.clone());
+        let vars = self
+            .labeled_variables
+            .entry(label)
+            .or_insert_with(|| Vec::with_capacity(1));
+        vars.push(var);
     }
-}
 
-impl From<Option<String>> for Label {
-    fn from(lbl: Option<String>) -> Self {
-        match lbl {
-            Some(lbl) => Label::from(lbl),
-            None => Label::empty(),
-        }
+    pub fn variables_with_label(&self, label: &Lbl) -> &[VarRef]
+    where
+        Lbl: Label,
+    {
+        self.labeled_variables.get(label).map(|v| v.as_ref()).unwrap_or(&[])
     }
 }
-impl<'a> From<&'a Option<String>> for Label {
-    fn from(lbl: &'a Option<String>) -> Self {
-        match lbl {
-            Some(lbl) => Label::from(lbl.as_str()),
-            None => Label::empty(),
-        }
-    }
-}
-impl<'a> From<Option<&'a str>> for Label {
-    fn from(lbl: Option<&'a str>) -> Self {
-        match lbl {
-            Some(lbl) => Label::from(lbl),
-            None => Label::empty(),
-        }
+
+impl<Lbl> Default for VariableLabels<Lbl> {
+    fn default() -> Self {
+        Self::new()
     }
 }
diff --git a/model/src/lang.rs b/model/src/lang.rs
index dc6589a6..9588b042 100644
--- a/model/src/lang.rs
+++ b/model/src/lang.rs
@@ -1,7 +1,9 @@
 mod atom;
 mod boolean;
-mod expr;
+pub mod expr;
 mod int;
+pub mod normal_form;
+pub mod reification;
 mod sym;
 mod variables;
 
@@ -45,8 +47,7 @@ impl VarRef {
 }
 
 pub use atom::Atom;
-pub use boolean::{BAtom, BExpr, BVar};
-pub use expr::{Expr, Fun};
+pub use boolean::BVar;
 pub use int::{IAtom, IVar};
 
 use crate::bounds::Lit;
@@ -73,7 +74,7 @@ pub enum ConversionError {
     TypeError,
     NotConstant,
     NotVariable,
-    NotBound,
+    NotLiteral,
     NotExpression,
     /// This conversion occurs when trying to convert an expression into a variable and that,
     /// there is a variable but its value is modified. For instance, this would
@@ -89,7 +90,7 @@ impl std::fmt::Display for ConversionError {
             ConversionError::NotVariable => write!(f, "not a variable"),
             ConversionError::NotPure => write!(f, "not a pure"),
             ConversionError::NotExpression => write!(f, "not an expression"),
-            ConversionError::NotBound => write!(f, "not a bound"),
+            ConversionError::NotLiteral => write!(f, "not a bound"),
         }
     }
 }
@@ -133,8 +134,7 @@ macro_rules! transitive_conversions {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::extensions::ExpressionFactoryExt;
-    use crate::Model;
+    type Model = crate::Model<&'static str>;
 
     fn check(m: &Model, x: impl Into<Atom>, result: &str) {
         assert_eq!(m.fmt(x).to_string(), result);
@@ -159,13 +159,13 @@ mod tests {
         let x = x + 1;
         check(&m, x, "b");
 
-        let x = m.leq(a + 1, 6);
-        check(&m, x, "(<= (+ a 1) 6)");
-
-        let x = m.eq(a - 3, b);
-        check(&m, x, "(= (- a 3) b)");
-
-        let x = m.implies(true, x);
-        check(&m, x, "(or false (= (- a 3) b))")
+        // let x = m.leq(a + 1, 6);
+        // check(&m, x, "(<= (+ a 1) 6)");
+        //
+        // let x = m.eq(a - 3, b);
+        // check(&m, x, "(= (- a 3) b)");
+        //
+        // let x = m.implies(true, x);
+        // check(&m, x, "(or false (= (- a 3) b))")
     }
 }
diff --git a/model/src/lang/atom.rs b/model/src/lang/atom.rs
index 21623f28..016534a7 100644
--- a/model/src/lang/atom.rs
+++ b/model/src/lang/atom.rs
@@ -4,7 +4,7 @@ use crate::symbols::TypedSym;
 
 #[derive(Hash, Eq, PartialEq, Copy, Clone, Debug)]
 pub enum Atom {
-    Bool(BAtom),
+    Bool(Lit),
     Int(IAtom),
     Sym(SAtom),
 }
@@ -22,25 +22,17 @@ impl Atom {
     /// This might fail in the case of a negated boolean or of a complex boolean expression.
     pub fn int_view(self) -> Option<IAtom> {
         match self {
-            Atom::Bool(b) => match b {
-                BAtom::Cst(x) => {
-                    if x {
-                        Some(1.into())
-                    } else {
-                        Some(0.into())
-                    }
-                }
-                BAtom::Literal(_) => None,
-                BAtom::Expr(_) => None,
-            },
+            Atom::Bool(Lit::TRUE) => Some(1.into()),
+            Atom::Bool(Lit::FALSE) => Some(0.into()),
+            Atom::Bool(_) => None,
             Atom::Int(i) => Some(i),
             Atom::Sym(s) => Some(s.int_view()),
         }
     }
 }
 
-impl From<BAtom> for Atom {
-    fn from(b: BAtom) -> Self {
+impl From<Lit> for Atom {
+    fn from(b: Lit) -> Self {
         Atom::Bool(b)
     }
 }
@@ -67,7 +59,7 @@ impl From<Variable> for Atom {
     }
 }
 
-impl TryFrom<Atom> for BAtom {
+impl TryFrom<Atom> for Lit {
     type Error = ConversionError;
 
     fn try_from(value: Atom) -> Result<Self, Self::Error> {
@@ -77,6 +69,7 @@ impl TryFrom<Atom> for BAtom {
         }
     }
 }
+
 impl TryFrom<Atom> for IAtom {
     type Error = ConversionError;
 
@@ -113,10 +106,9 @@ impl TryFrom<Atom> for Variable {
 
 use crate::transitive_conversions;
 use std::convert::TryInto;
-transitive_conversions!(Atom, BAtom, Lit);
-transitive_conversions!(Atom, BAtom, BExpr);
-transitive_conversions!(Atom, BAtom, bool);
+
 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 3824dbe9..06dc5e1e 100644
--- a/model/src/lang/boolean.rs
+++ b/model/src/lang/boolean.rs
@@ -1,9 +1,5 @@
-use crate::expressions::ExprHandle;
-
 use crate::bounds::Lit;
-use crate::lang::{ConversionError, IVar, VarRef};
-use std::cmp::Ordering;
-use std::convert::TryFrom;
+use crate::lang::{IVar, VarRef};
 
 /// A boolean variable.
 /// It is a wrapper around an (untyped) discrete variable to provide type safety.
@@ -60,119 +56,3 @@ impl std::ops::Not for BVar {
         self.false_lit()
     }
 }
-
-// equivalent to lit
-#[derive(Hash, Ord, PartialOrd, Eq, PartialEq, Copy, Clone, Debug)]
-pub enum BAtom {
-    Cst(bool),
-    Literal(Lit),
-    Expr(BExpr),
-}
-
-#[derive(Hash, Ord, PartialOrd, Eq, PartialEq, Copy, Clone, Debug)]
-pub struct BExpr {
-    pub expr: ExprHandle,
-    pub negated: bool,
-}
-
-impl std::ops::Not for BExpr {
-    type Output = Self;
-
-    fn not(self) -> Self::Output {
-        BExpr {
-            expr: self.expr,
-            negated: !self.negated,
-        }
-    }
-}
-
-impl BAtom {
-    pub fn lexical_cmp(&self, other: &BAtom) -> Ordering {
-        match self {
-            BAtom::Cst(c1) => match other {
-                BAtom::Cst(c2) => c1.cmp(c2),
-                _ => Ordering::Less,
-            },
-            BAtom::Literal(x) => match other {
-                BAtom::Literal(y) => x.lexical_cmp(y),
-                BAtom::Cst(_) => Ordering::Greater,
-                BAtom::Expr(_) => Ordering::Less,
-            },
-
-            BAtom::Expr(e) => match other {
-                BAtom::Expr(e2) => e.cmp(e2),
-                BAtom::Cst(_) => Ordering::Greater,
-                BAtom::Literal(_) => Ordering::Greater,
-            },
-        }
-    }
-}
-
-impl std::ops::Not for BAtom {
-    type Output = BAtom;
-
-    fn not(self) -> Self::Output {
-        match self {
-            BAtom::Cst(x) => BAtom::Cst(!x),
-            BAtom::Literal(b) => BAtom::Literal(!b),
-            BAtom::Expr(e) => BAtom::Expr(!e),
-        }
-    }
-}
-
-impl From<bool> for BAtom {
-    fn from(value: bool) -> Self {
-        BAtom::Cst(value)
-    }
-}
-
-impl From<Lit> for BAtom {
-    fn from(bnd: Lit) -> Self {
-        BAtom::Literal(bnd)
-    }
-}
-
-impl From<BVar> for BAtom {
-    fn from(var: BVar) -> Self {
-        BAtom::Literal(var.into())
-    }
-}
-
-impl From<BExpr> for BAtom {
-    fn from(e: BExpr) -> Self {
-        BAtom::Expr(e)
-    }
-}
-
-impl TryFrom<BAtom> for bool {
-    type Error = ConversionError;
-
-    fn try_from(value: BAtom) -> Result<Self, Self::Error> {
-        match value {
-            BAtom::Cst(b) => Ok(b),
-            _ => Err(ConversionError::NotConstant),
-        }
-    }
-}
-
-impl TryFrom<BAtom> for Lit {
-    type Error = ConversionError;
-
-    fn try_from(value: BAtom) -> Result<Self, Self::Error> {
-        match value {
-            BAtom::Literal(b) => Ok(b),
-            _ => Err(ConversionError::NotBound),
-        }
-    }
-}
-
-impl TryFrom<BAtom> for BExpr {
-    type Error = ConversionError;
-
-    fn try_from(value: BAtom) -> Result<Self, Self::Error> {
-        match value {
-            BAtom::Expr(b) => Ok(b),
-            _ => Err(ConversionError::NotExpression),
-        }
-    }
-}
diff --git a/model/src/lang/expr.rs b/model/src/lang/expr.rs
index 1373e473..3de36afd 100644
--- a/model/src/lang/expr.rs
+++ b/model/src/lang/expr.rs
@@ -1,53 +1,156 @@
-use crate::lang::Atom;
+use crate::bounds::{Disjunction, Lit};
+use crate::lang::normal_form::{NFEq, NFLeq, NFOptEq, NFOptLeq, NormalExpr};
+use crate::lang::reification::ReifiableExpr;
+use crate::lang::{Atom, IAtom};
 
-pub type Args = Vec<Atom>;
+/// Trait denoting the capability of transforming an expression into its normal form.
+///
+/// This transformation is necessary to transform user defined constraints into constraints
+/// that can be efficiently reified.
+///
+/// # Example transformations
+///
+/// ```text
+/// - (a = b)   =>  (a = b)
+/// - (a != b)  =>  not(a = b)
+/// - (b = a)   =>  (a = b)     # sorted by lexical order
+/// ```
+/// After normalization, a single literal is necessary to reify the 3 expressions.
+pub trait Normalize<X: ReifiableExpr> {
+    fn normalize(&self) -> NormalExpr<X>;
+}
+impl Normalize<Never> for Lit {
+    fn normalize(&self) -> NormalExpr<Never> {
+        NormalExpr::Literal(*self)
+    }
+}
+
+/// A type that can never be created.
+#[derive(Eq, PartialEq, Hash, Debug)]
+pub struct Never(());
+
+pub fn leq(lhs: impl Into<IAtom>, rhs: impl Into<IAtom>) -> Leq {
+    Leq(lhs.into(), rhs.into())
+}
+
+pub fn lt(lhs: impl Into<IAtom>, rhs: impl Into<IAtom>) -> Leq {
+    leq(lhs.into(), rhs.into() - 1)
+}
+
+pub fn geq(lhs: impl Into<IAtom>, rhs: impl Into<IAtom>) -> Leq {
+    leq(rhs, lhs)
+}
+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 eq(lhs: impl Into<Atom>, rhs: impl Into<Atom>) -> Eq {
+    Eq(lhs.into(), rhs.into())
+}
+
+pub fn neq(lhs: impl Into<Atom>, rhs: impl Into<Atom>) -> Neq {
+    Neq(lhs.into(), rhs.into())
+}
+
+pub fn opt_eq(lhs: impl Into<Atom>, rhs: impl Into<Atom>) -> OptEq {
+    OptEq(lhs.into(), rhs.into())
+}
+
+pub fn or(disjuncts: impl Into<Box<[Lit]>>) -> Or {
+    Or(disjuncts.into())
+}
+pub fn and(disjuncts: impl Into<Box<[Lit]>>) -> And {
+    And(disjuncts.into())
+}
+pub fn implies(a: impl Into<Lit>, b: impl Into<Lit>) -> Or {
+    or([!a.into(), b.into()])
+}
+
+pub struct Or(Box<[Lit]>);
 
-#[derive(Hash, Ord, PartialOrd, Eq, PartialEq, Copy, Clone)]
-pub enum Fun {
-    Or,
-    Eq,
-    Leq,
-    Max,
-    /// Equality between two optional variables: if both are present, then they must be equal.
-    OptEq,
-    /// Inequality that must hold if both (optional) variables are present.
-    OptLeq,
+impl Normalize<Disjunction> for Or {
+    fn normalize(&self) -> NormalExpr<Disjunction> {
+        let vec = self.0.iter().copied().collect();
+        if let Some(disj) = Disjunction::new_non_tautological(vec) {
+            NormalExpr::Pos(disj)
+        } else {
+            NormalExpr::Literal(Lit::TRUE)
+        }
+    }
 }
 
-impl std::fmt::Display for Fun {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        write!(
-            f,
-            "{}",
-            match self {
-                Fun::Or => "or",
-                Fun::Eq => "=",
-                Fun::Leq => "<=",
-                Fun::Max => "max",
-                Fun::OptEq => "opt_eq",
-                Fun::OptLeq => "opt_leq",
-            }
-        )
+pub struct And(Box<[Lit]>);
+
+impl Normalize<Disjunction> for And {
+    fn normalize(&self) -> NormalExpr<Disjunction> {
+        // (and a b c)  <=>  !(or !a !b !c)
+        let vec = self.0.iter().copied().map(|l| !l).collect();
+        if let Some(disj) = Disjunction::new_non_tautological(vec) {
+            NormalExpr::Neg(disj)
+        } else {
+            NormalExpr::Literal(Lit::FALSE)
+        }
     }
 }
-impl std::fmt::Debug for Fun {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        write!(f, "{}", self)
+
+#[derive(Copy, Clone, Debug)]
+pub struct Leq(IAtom, IAtom);
+
+impl std::ops::Not for Leq {
+    type Output = Leq;
+
+    fn not(self) -> Self::Output {
+        gt(self.0, self.1)
     }
 }
 
-#[derive(Hash, Eq, PartialEq, Clone, Debug)]
-pub struct Expr {
-    pub fun: Fun,
-    pub args: Args,
+impl Normalize<NFLeq> for Leq {
+    fn normalize(&self) -> NormalExpr<NFLeq> {
+        NFLeq::leq(self.0, self.1)
+    }
 }
-impl Expr {
-    pub fn new(fun: Fun, args: Args) -> Expr {
-        Expr { fun, args }
+
+#[derive(Copy, Clone, Debug)]
+pub struct Eq(Atom, Atom);
+
+impl Normalize<NFEq> for Eq {
+    fn normalize(&self) -> NormalExpr<NFEq> {
+        NFEq::eq(self.0, self.1)
     }
+}
+
+pub struct Neq(Atom, Atom);
+
+impl Normalize<NFEq> for Neq {
+    fn normalize(&self) -> NormalExpr<NFEq> {
+        !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);
 
-    pub fn new2(fun: Fun, arg1: impl Into<Atom>, arg2: impl Into<Atom>) -> Expr {
-        let args = vec![arg1.into(), arg2.into()];
-        Expr { fun, args }
+impl Normalize<NFOptEq> for OptEq {
+    fn normalize(&self) -> NormalExpr<NFOptEq> {
+        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/int.rs b/model/src/lang/int.rs
index ef3bf348..b257b7f8 100644
--- a/model/src/lang/int.rs
+++ b/model/src/lang/int.rs
@@ -7,7 +7,9 @@ use std::convert::TryFrom;
 pub struct IVar(VarRef);
 
 impl IVar {
-    pub fn new(dvar: VarRef) -> Self {
+    pub const ZERO: IVar = IVar(VarRef::ZERO);
+
+    pub const fn new(dvar: VarRef) -> Self {
         IVar(dvar)
     }
 }
@@ -18,34 +20,36 @@ impl From<IVar> for VarRef {
     }
 }
 
-// var + cst
+/// An int-valued atom `(variable + constant)`
+/// It can be used to represent a constant value by using [IVar::ZERO] as the variable.
 #[derive(Hash, Eq, PartialEq, Copy, Clone, Debug)]
 pub struct IAtom {
-    pub var: Option<IVar>,
+    pub var: IVar,
     pub shift: IntCst,
 }
 impl IAtom {
-    pub fn new(var: Option<IVar>, shift: IntCst) -> IAtom {
+    pub const ZERO: IAtom = IAtom {
+        var: IVar::ZERO,
+        shift: 0,
+    };
+    pub fn new(var: IVar, shift: IntCst) -> IAtom {
         IAtom { var, shift }
     }
 
     /// A total order between the names of the atoms, not on their expected values.
     pub fn lexical_cmp(&self, other: &IAtom) -> Ordering {
-        match (self.var, other.var) {
-            (Some(v1), Some(v2)) if v1 != v2 => v1.cmp(&v2),
-            (Some(_), None) => Ordering::Greater,
-            (None, Some(_)) => Ordering::Less,
-            _ => self.shift.cmp(&other.shift),
-        }
+        self.var.cmp(&other.var).then(self.shift.cmp(&other.shift))
     }
 
     /// Returns a literal representing whether this atom is lesser than the given value.
     pub fn lt_lit(self, value: IntCst) -> Lit {
         let rhs = value - self.shift;
-        match self.var {
-            Some(v) => VarRef::from(v).lt(rhs),
-            None if 0 < rhs => Lit::TRUE,
-            _ => Lit::FALSE,
+        if self.var != IVar::ZERO {
+            VarRef::from(self.var).lt(rhs)
+        } else if 0 < rhs {
+            Lit::TRUE
+        } else {
+            Lit::FALSE
         }
     }
 }
@@ -65,12 +69,17 @@ impl PartialOrd for IAtom {
 
 impl From<IVar> for IAtom {
     fn from(v: IVar) -> Self {
-        IAtom::new(Some(v), 0)
+        IAtom::new(v, 0)
+    }
+}
+impl From<VarRef> for IAtom {
+    fn from(v: VarRef) -> Self {
+        IAtom::new(IVar::new(v), 0)
     }
 }
 impl From<IntCst> for IAtom {
     fn from(i: i32) -> Self {
-        IAtom::new(None, i)
+        IAtom::new(IVar::ZERO, i)
     }
 }
 
@@ -78,15 +87,10 @@ impl TryFrom<IAtom> for IVar {
     type Error = ConversionError;
 
     fn try_from(value: IAtom) -> Result<Self, Self::Error> {
-        match value.var {
-            None => Err(ConversionError::NotVariable),
-            Some(v) => {
-                if value.shift == 0 {
-                    Ok(v)
-                } else {
-                    Err(ConversionError::NotPure)
-                }
-            }
+        if value.shift == 0 {
+            Ok(value.var)
+        } else {
+            Err(ConversionError::NotPure)
         }
     }
 }
@@ -96,8 +100,8 @@ impl TryFrom<IAtom> for IntCst {
 
     fn try_from(value: IAtom) -> Result<Self, Self::Error> {
         match value.var {
-            None => Ok(value.shift),
-            Some(_) => Err(ConversionError::NotConstant),
+            IVar::ZERO => Ok(value.shift),
+            _ => Err(ConversionError::NotConstant),
         }
     }
 }
@@ -113,7 +117,7 @@ impl std::ops::Add<IntCst> for IVar {
     type Output = IAtom;
 
     fn add(self, rhs: IntCst) -> Self::Output {
-        IAtom::new(Some(self), rhs)
+        IAtom::new(self, rhs)
     }
 }
 impl std::ops::Sub<IntCst> for IAtom {
@@ -127,6 +131,6 @@ impl std::ops::Sub<IntCst> for IVar {
     type Output = IAtom;
 
     fn sub(self, rhs: IntCst) -> Self::Output {
-        IAtom::new(Some(self), -rhs)
+        IAtom::new(self, -rhs)
     }
 }
diff --git a/model/src/lang/normal_form.rs b/model/src/lang/normal_form.rs
new file mode 100644
index 00000000..b7738a13
--- /dev/null
+++ b/model/src/lang/normal_form.rs
@@ -0,0 +1,158 @@
+use crate::bounds::Lit;
+use crate::lang::*;
+use std::cmp::Ordering;
+
+/// Normal form of an expression with component of type `X`.
+/// It can be either: a literal, an expression of type `X`
+/// or the negation of an expression of type `X`.
+pub enum NormalExpr<X> {
+    Literal(Lit),
+    Pos(X),
+    Neg(X),
+}
+
+impl<T> std::ops::Not for NormalExpr<T> {
+    type Output = Self;
+
+    fn not(self) -> Self::Output {
+        match self {
+            NormalExpr::Literal(l) => NormalExpr::Literal(!l),
+            NormalExpr::Pos(e) => NormalExpr::Neg(e),
+            NormalExpr::Neg(e) => NormalExpr::Pos(e),
+        }
+    }
+}
+
+impl<T> From<Lit> for NormalExpr<T> {
+    fn from(l: Lit) -> Self {
+        Self::Literal(l)
+    }
+}
+impl<T> From<bool> for NormalExpr<T> {
+    fn from(l: bool) -> Self {
+        Self::Literal(l.into())
+    }
+}
+
+/// Canonical representation of an inequality:  `lhs <= rhs + rhs_add`
+#[derive(Debug, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash)]
+pub struct NFLeq {
+    pub lhs: VarRef,
+    pub rhs: VarRef,
+    pub rhs_add: IntCst,
+}
+
+impl NFLeq {
+    pub fn new(lhs: VarRef, rhs: VarRef, rhs_add: IntCst) -> Self {
+        debug_assert!(lhs < rhs, "Violated lexical order invariant");
+        debug_assert_ne!(lhs, VarRef::ZERO);
+        debug_assert_ne!(rhs, VarRef::ZERO);
+        NFLeq { lhs, rhs, rhs_add }
+    }
+    pub fn geq<A: Into<IAtom>, B: Into<IAtom>>(a: A, b: B) -> NormalExpr<Self> {
+        Self::leq(b, a)
+    }
+    pub fn gt<A: Into<IAtom>, B: Into<IAtom>>(a: A, b: B) -> NormalExpr<Self> {
+        Self::lt(b, a)
+    }
+    pub fn lt<A: Into<IAtom>, B: Into<IAtom>>(a: A, b: B) -> NormalExpr<Self> {
+        Self::leq(a.into() + 1, b)
+    }
+    pub fn leq<A: Into<IAtom>, B: Into<IAtom>>(lhs: A, rhs: B) -> NormalExpr<Self> {
+        let lhs = lhs.into();
+        let rhs = rhs.into();
+
+        // normalize, transfer the shift from right to left
+        let rhs_add = rhs.shift - lhs.shift;
+        let lhs: VarRef = lhs.var.into();
+        let rhs: VarRef = rhs.var.into();
+
+        // 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
+        if lhs == rhs {
+            // X  <= X + rhs_add   <=>  0 <= rhs_add
+            return (0 <= rhs_add).into();
+        }
+        if rhs == VarRef::ZERO {
+            // lhs  <= rhs_add
+            return Lit::leq(lhs, rhs_add).into();
+        }
+        if lhs == VarRef::ZERO {
+            // 0 <= rhs + rhs_add   <=>  -rhs_add <= rhs
+            return Lit::geq(rhs, -rhs_add).into();
+        }
+
+        // maintain the invariant that left side of the LEQ has a small lexical order
+        match lhs.cmp(&rhs) {
+            Ordering::Less => NormalExpr::Pos(Self::new(lhs, rhs, rhs_add)),
+            Ordering::Equal => unreachable!(),
+            Ordering::Greater => {
+                // (lhs <= rhs + rhs_add) <=> (not (rhs <= lhs - rhs_add -1)
+                NormalExpr::Neg(Self::new(rhs, lhs, -rhs_add - 1))
+            }
+        }
+    }
+}
+
+/// lhs = rhs + rhs_add
+#[derive(Debug, Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Hash)]
+pub struct NFEq {
+    pub lhs: VarRef,
+    pub rhs: VarRef,
+    pub rhs_add: IntCst,
+}
+
+impl NFEq {
+    pub fn new(lhs: VarRef, rhs: VarRef, rhs_add: IntCst) -> Self {
+        debug_assert!(lhs < rhs, "Invariant violated");
+        NFEq { lhs, rhs, rhs_add }
+    }
+
+    /// Builds the normal form for this equality.
+    pub fn eq<A: Into<Atom>, B: Into<Atom>>(a: A, b: B) -> NormalExpr<NFEq> {
+        let a = a.into();
+        let b = b.into();
+        if a == b {
+            Lit::TRUE.into()
+        } else if a.kind() != b.kind() {
+            Lit::FALSE.into()
+        } else {
+            use Atom::*;
+            match (a, b) {
+                (Bool(_a), Bool(_b)) => todo!(),
+                (Int(a), Int(b)) => Self::int_eq(a, b),
+                (Sym(a), Sym(b)) => Self::int_eq(a.int_view(), b.int_view()),
+                _ => unreachable!(), // guarded by kind comparison
+            }
+        }
+    }
+
+    pub fn int_eq(a: IAtom, b: IAtom) -> NormalExpr<NFEq> {
+        let b_add = b.shift - a.shift;
+        let a: VarRef = a.var.into();
+        let b: VarRef = b.var.into();
+
+        match a.cmp(&b) {
+            Ordering::Less => NormalExpr::Pos(Self::new(a, b, b_add)),
+            Ordering::Equal => (b_add == 0).into(),
+            Ordering::Greater => NormalExpr::Pos(Self::new(b, a, -b_add)),
+        }
+    }
+}
+
+/// (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,
+}
diff --git a/model/src/lang/reification.rs b/model/src/lang/reification.rs
new file mode 100644
index 00000000..33ad9ec9
--- /dev/null
+++ b/model/src/lang/reification.rs
@@ -0,0 +1,242 @@
+use crate::bounds::Lit;
+use crate::lang::expr::Normalize;
+use crate::lang::normal_form::NormalExpr;
+use std::any::{Any, TypeId};
+use std::collections::HashMap;
+use std::fmt::Debug;
+use std::hash::{Hash, Hasher};
+use std::sync::Arc;
+
+/// Automatically derived trait for the capabilities of an expression that
+/// requires to be inserted into the model (enforced or reified).
+///
+/// It is typically required that the normal form of any user-facing expression
+/// be reifiable.
+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 {}
+
+/// Dynamically typed interface to an expression in the model.
+///
+/// To transform this dynamically typed expression to its original type,
+/// one can use the [`downcast`].
+pub type Expr = dyn ExprInterface;
+
+/// If the (dynamically typed) referenced `expr` is of type `T`, then casts
+/// it and returns it. Returns None if the expr is not of type `T`.
+///
+/// Source: this implementation is taken from the one of `core::any::downcast_ref` which
+/// cannot be directly called as a method on an `&Expr`.
+pub fn downcast<T: Any + 'static>(expr: &Expr) -> Option<&T> {
+    let tpe = expr.type_id();
+    if tpe == TypeId::of::<T>() {
+        // CONTEXT: a reference to trait-object is composed of two pointers.
+        //  - the first one is a pointer to the original data (here of type T).
+        //  - the second one is a pointer to the vtable (here of the trait ExprInterface)
+        // In the following, we cast the first pointer to a pointer on T.
+
+        // SAFETY: just checked whether we are pointing to the correct type, and we can rely on
+        // that check for memory safety because we have implemented Any for all types; no other
+        // impls can exist as they would conflict with our impl.
+        unsafe { Some(&*(expr as *const Expr as *const T)) }
+    } else {
+        None
+    }
+}
+
+/// A structure to keep track of all binding events that result from reifying and enforcing expressions.
+///
+/// A correspondence between canonical expression and the literal that they have reified to is maintained.
+/// When two literals are bound to the same expression, we record a binding event between the two literals.
+///
+/// The structure maintains a list of binding events that can be iterated with a [BindingCursor]
+#[derive(Default, Clone)]
+pub struct Reification {
+    /// Associates each canonical atom to a single literal.
+    map: HashMap<WrappedExpr, Lit>,
+    /// All binding events in chronological order. This is intended to easily process
+    /// binding events and detect whether new events have been added.
+    binding_events: Vec<(Lit, BindTarget)>,
+}
+
+impl Reification {
+    fn add(&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();
+        match nf {
+            NormalExpr::Literal(l) => Some(l),
+            NormalExpr::Pos(e) => self.map.get(&WrappedExpr::new(e)).copied(),
+            NormalExpr::Neg(e) => self.map.get(&WrappedExpr::new(e)).map(|&l| !l),
+        }
+    }
+
+    /// 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();
+        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),
+        }
+    }
+
+    /// Binds the user-provided expression `expr` to the given `literal`.
+    pub fn bind<X: Normalize<T>, T: ReifiableExpr>(&mut self, expr: X, literal: Lit) {
+        let nf = expr.normalize();
+        match nf {
+            NormalExpr::Literal(l) => self.bind_literals(l, literal),
+            NormalExpr::Pos(e) => {
+                let e = WrappedExpr::new(e);
+                if let Some(prev) = self.map.get(&e).copied() {
+                    self.bind_literals(prev, literal)
+                } else {
+                    self.add(e, literal);
+                }
+            }
+            NormalExpr::Neg(e) => {
+                let e = WrappedExpr::new(e);
+                if let Some(prev) = self.map.get(&e).copied() {
+                    self.bind_literals(prev, !literal)
+                } else {
+                    self.add(e, !literal);
+                }
+            }
+        }
+    }
+
+    /// Impose the constraint that the two literals are equal.
+    pub fn bind_literals(&mut self, l1: Lit, l2: Lit) {
+        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);
+        if ret.is_some() {
+            cursor.0 += 1;
+        }
+        ret
+    }
+}
+
+/// A wrapper around a dynamically typed expression to allow usage as a key in hash map.
+#[derive(Clone)]
+struct WrappedExpr {
+    /// The expression that we are interested in.
+    value: Arc<Expr>,
+    /// A lambda function to allow comparing two expressions (the first one being of the same type as `value`
+    eq_any: Arc<dyn Fn(&Expr, &Expr) -> bool + Send + Sync>,
+    /// A precomputed hash of `value`.
+    hash: u64,
+}
+
+impl WrappedExpr {
+    /// Creates a wrapper around this expression that provides [PartialEq] and [Hash] implementations.
+    pub fn new<E: ReifiableExpr>(x: E) -> Self {
+        let mut hasher = std::collections::hash_map::DefaultHasher::new();
+        x.hash(&mut hasher);
+        let hash = hasher.finish();
+        let value = Arc::new(x);
+        let eq_any = Arc::new(|x: &Expr, y: &Expr| {
+            let x = downcast::<E>(x).expect("wrong type");
+            if let Some(y) = downcast::<E>(y) {
+                // x and y are both of type E, we can compare them.
+                return x == y;
+            }
+            false // objects are of different types
+        });
+
+        WrappedExpr { value, eq_any, hash }
+    }
+}
+
+impl Hash for WrappedExpr {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self.value.type_id().hash(state);
+        self.hash.hash(state);
+    }
+}
+
+impl PartialEq for WrappedExpr {
+    fn eq(&self, other: &Self) -> bool {
+        // use the equality lambda built in the object to compare to others.
+        self.eq_any.as_ref()(self.value.as_ref(), other.value.as_ref())
+    }
+}
+impl Eq for WrappedExpr {}
+
+/// Target of the binding of literal: either an expression or another literal.
+#[derive(Clone, Debug)]
+pub enum BindTarget {
+    Expr(Arc<Expr>),
+    Literal(Lit),
+}
+
+/// A cursor into a sequence of bindings.
+///
+/// The next bidding can be popped with [Reification::pop_next_event].
+#[derive(Copy, Clone)]
+pub struct BindingCursor(usize);
+
+impl BindingCursor {
+    /// Creates a new cursor on the first event.
+    pub fn first() -> Self {
+        BindingCursor(0)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::bounds::Lit;
+    use crate::lang::expr::*;
+    use crate::lang::reification::Reification;
+    use crate::lang::{IVar, VarRef};
+
+    const A: IVar = IVar::new(VarRef::from_u32(1));
+    const B: IVar = IVar::new(VarRef::from_u32(2));
+    const C: IVar = IVar::new(VarRef::from_u32(3));
+
+    #[test]
+    fn test_reif() {
+        let t = Lit::TRUE;
+        let f = Lit::FALSE;
+        let l1 = leq(A, B + 3);
+        let l2 = leq(A, C);
+        let e1 = eq(A, B + 3);
+
+        let mut reif = Reification::default();
+        reif.bind(l1, t);
+        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));
+
+        // same as l1
+        let l1_prime = geq(B + 3, A);
+        assert_eq!(reif.interned(l1_prime), Some(t));
+
+        // inverse of l1, should return the opposite literal
+        assert_eq!(reif.interned(!l1), Some(f));
+    }
+}
diff --git a/model/src/lang/sym.rs b/model/src/lang/sym.rs
index c9f631e4..093e961e 100644
--- a/model/src/lang/sym.rs
+++ b/model/src/lang/sym.rs
@@ -41,8 +41,8 @@ impl SAtom {
 
     pub fn int_view(self) -> IAtom {
         match self {
-            SAtom::Var(v) => IAtom::new(Some(IVar::new(v.var)), 0),
-            SAtom::Cst(s) => IAtom::new(None, usize::from(s.sym) as IntCst),
+            SAtom::Var(v) => IAtom::new(IVar::new(v.var), 0),
+            SAtom::Cst(s) => IAtom::new(IVar::ZERO, usize::from(s.sym) as IntCst),
         }
     }
 }
diff --git a/model/src/lib.rs b/model/src/lib.rs
index d01cb909..abe1e123 100644
--- a/model/src/lib.rs
+++ b/model/src/lib.rs
@@ -1,10 +1,5 @@
-// reexport the Label type
-pub use label::Label;
-pub use model::*;
-
-pub mod bindings;
 pub mod bounds;
-pub mod expressions;
+pub mod decomposition;
 pub mod extensions;
 mod label;
 pub mod lang;
@@ -12,3 +7,6 @@ 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 e61f4640..d013dfff 100644
--- a/model/src/model.rs
+++ b/model/src/model.rs
@@ -1,43 +1,39 @@
-use crate::bindings::Bindings;
 use crate::bounds::Lit;
-use crate::expressions::*;
-use crate::extensions::{AssignmentExt, ExpressionFactoryExt, SavedAssignment, Shaped};
+use crate::extensions::{AssignmentExt, SavedAssignment, Shaped};
+use crate::label::{Label, VariableLabels};
+use crate::lang::expr::Normalize;
+use crate::lang::reification::{ReifiableExpr, Reification};
 use crate::lang::*;
 use crate::state::*;
 use crate::symbols::SymbolTable;
 use crate::types::TypeId;
-use crate::Label;
 use aries_backtrack::{Backtrack, DecLvl};
 use aries_collections::ref_store::RefMap;
 use std::sync::Arc;
 
 /// Defines the structure of a model: variables names, types, relations, ...
 #[derive(Clone)]
-pub struct ModelShape {
+pub struct ModelShape<Lbl> {
     pub symbols: Arc<SymbolTable>,
     pub types: RefMap<VarRef, Type>,
-    pub expressions: Expressions,
-    pub bindings: Bindings,
-    pub labels: RefMap<VarRef, String>,
+    pub expressions: Reification,
+    pub labels: VariableLabels<Lbl>,
     num_writers: u8,
 }
 
-impl ModelShape {
+impl<Lbl: Label> ModelShape<Lbl> {
     pub fn new() -> Self {
         Self::new_with_symbols(Arc::new(SymbolTable::empty()))
     }
 
     pub fn new_with_symbols(symbols: Arc<SymbolTable>) -> Self {
-        let mut m = ModelShape {
+        ModelShape {
             symbols,
             types: Default::default(),
             expressions: Default::default(),
-            bindings: Default::default(),
             labels: Default::default(),
             num_writers: 0,
-        };
-        m.set_label(VarRef::ZERO, "ZERO");
-        m
+        }
     }
 
     pub fn new_write_token(&mut self) -> WriterId {
@@ -45,46 +41,38 @@ impl ModelShape {
         WriterId(self.num_writers - 1)
     }
 
-    fn set_label(&mut self, var: VarRef, l: impl Into<Label>) {
-        if let Some(str) = l.into().lbl {
-            self.labels.insert(var, str)
+    fn set_label(&mut self, var: VarRef, l: impl Into<Lbl>) {
+        self.labels.insert(var, l.into())
+    }
+    pub fn get_variable(&self, label: &Lbl) -> Option<VarRef> {
+        match *self.labels.variables_with_label(label) {
+            [] => None,
+            [var] => Some(var),
+            _ => panic!("More than one variable with label: {:?}", label),
         }
     }
     fn set_type(&mut self, var: VarRef, typ: Type) {
         self.types.insert(var, typ);
     }
-
-    pub fn intern_boolean_expression(&mut self, e: Expr) -> BExpr {
-        let handle = self.expressions.intern(e);
-        BExpr {
-            expr: handle,
-            negated: false,
-        }
-    }
-
-    // ======= Expression reification =====
-
-    pub fn interned_expr(&self, handle: ExprHandle) -> Option<Lit> {
-        self.bindings.as_lit(handle)
-    }
 }
 
-impl Default for ModelShape {
+impl<Lbl: Label> Default for ModelShape<Lbl> {
     fn default() -> Self {
         Self::new()
     }
 }
 
-///
+/// Description problem, composed of its shape (variable declaration, composed
+/// expressions, ...) and its state (the currently admissible values for each variable).
 #[derive(Clone)]
-pub struct Model {
+pub struct Model<Lbl> {
     /// Structure of the model and metadata of its various components.
-    pub shape: ModelShape,
+    pub shape: ModelShape<Lbl>,
     /// Domain of all variables, defining the current state of the Model.
     pub state: Domains,
 }
 
-impl Model {
+impl<Lbl: Label> Model<Lbl> {
     pub fn new() -> Self {
         Self::new_with_symbols(Arc::new(SymbolTable::empty()))
     }
@@ -100,15 +88,15 @@ impl Model {
         self.shape.new_write_token()
     }
 
-    pub fn new_bvar<L: Into<Label>>(&mut self, label: L) -> BVar {
+    pub fn new_bvar(&mut self, label: impl Into<Lbl>) -> BVar {
         self.create_bvar(None, label)
     }
 
-    pub fn new_optional_bvar<L: Into<Label>>(&mut self, presence: Lit, label: L) -> BVar {
+    pub fn new_optional_bvar(&mut self, presence: Lit, label: impl Into<Lbl>) -> BVar {
         self.create_bvar(Some(presence), label)
     }
 
-    pub fn new_presence_variable(&mut self, scope: Lit, label: impl Into<Label>) -> BVar {
+    pub fn new_presence_variable(&mut self, scope: Lit, label: impl Into<Lbl>) -> BVar {
         let lit = self.state.new_presence_literal(scope);
         let var = lit.variable();
         self.shape.set_label(var, label);
@@ -116,7 +104,7 @@ impl Model {
         BVar::new(var)
     }
 
-    fn create_bvar(&mut self, presence: Option<Lit>, label: impl Into<Label>) -> BVar {
+    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)
         } else {
@@ -127,15 +115,15 @@ impl Model {
         BVar::new(dvar)
     }
 
-    pub fn new_ivar(&mut self, lb: IntCst, ub: IntCst, label: impl Into<Label>) -> IVar {
+    pub fn new_ivar(&mut self, lb: IntCst, ub: IntCst, label: impl Into<Lbl>) -> IVar {
         self.create_ivar(lb, ub, None, label)
     }
 
-    pub fn new_optional_ivar(&mut self, lb: IntCst, ub: IntCst, presence: Lit, label: impl Into<Label>) -> IVar {
+    pub fn new_optional_ivar(&mut self, lb: IntCst, ub: IntCst, presence: Lit, label: impl Into<Lbl>) -> IVar {
         self.create_ivar(lb, ub, Some(presence), label)
     }
 
-    fn create_ivar(&mut self, lb: IntCst, ub: IntCst, presence: Option<Lit>, label: impl Into<Label>) -> IVar {
+    fn create_ivar(&mut self, lb: IntCst, ub: IntCst, presence: Option<Lit>, label: impl Into<Lbl>) -> IVar {
         let dvar = if let Some(presence) = presence {
             self.state.new_optional_var(lb, ub, presence)
         } else {
@@ -146,15 +134,15 @@ impl Model {
         IVar::new(dvar)
     }
 
-    pub fn new_sym_var(&mut self, tpe: TypeId, label: impl Into<Label>) -> SVar {
+    pub fn new_sym_var(&mut self, tpe: TypeId, label: impl Into<Lbl>) -> SVar {
         self.create_sym_var(tpe, None, label)
     }
 
-    pub fn new_optional_sym_var(&mut self, tpe: TypeId, presence: impl Into<Lit>, label: impl Into<Label>) -> SVar {
+    pub fn new_optional_sym_var(&mut self, tpe: TypeId, presence: impl Into<Lit>, label: impl Into<Lbl>) -> SVar {
         self.create_sym_var(tpe, Some(presence.into()), label)
     }
 
-    fn create_sym_var(&mut self, tpe: TypeId, presence: Option<Lit>, label: impl Into<Label>) -> SVar {
+    fn create_sym_var(&mut self, tpe: TypeId, presence: Option<Lit>, label: impl Into<Lbl>) -> SVar {
         let instances = self.shape.symbols.instances_of_type(tpe);
         if let Some((lb, ub)) = instances.bounds() {
             let lb = usize::from(lb) as IntCst;
@@ -204,72 +192,46 @@ impl Model {
         }
     }
 
-    fn create_binding(&mut self, expr: ExprHandle) -> Lit {
-        if let Some(lit) = self.shape.interned_expr(expr) {
-            lit
-        } else {
-            let name = format!("{}", self.fmt(BAtom::Expr(BExpr { expr, negated: false })));
-            let var = self.new_bvar(name);
-            let lit = var.true_lit();
-            self.record_binding(expr, lit);
-            lit
+    /// Interns the given expression and returns an equivalent literal.
+    /// 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);
         }
+        lit
     }
 
-    fn record_binding(&mut self, handle: ExprHandle, literal: Lit) {
-        assert!(self.shape.bindings.as_lit(handle).is_none());
-        self.shape.bindings.bind(handle, literal);
+    pub fn enforce<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, b: Expr) {
+        self.shape.expressions.bind(b, Lit::TRUE)
     }
 
-    pub fn reify(&mut self, b: BAtom) -> Lit {
-        match b {
-            BAtom::Cst(true) => Lit::TRUE,
-            BAtom::Cst(false) => Lit::FALSE,
-            BAtom::Literal(b) => b,
-            BAtom::Expr(e) => {
-                let BExpr { expr: handle, negated } = e;
-                let lit = self.create_binding(handle);
-                if negated {
-                    !lit
-                } else {
-                    lit
-                }
-            }
+    pub fn enforce_all<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, bools: impl IntoIterator<Item = Expr>) {
+        for b in bools {
+            self.enforce(b);
         }
     }
 
-    fn make_literals_eq(&mut self, a: Lit, b: Lit) {
-        self.shape.bindings.bind_literals(a, b);
-    }
-
-    pub fn enforce(&mut self, b: BAtom) {
-        self.bind(b, Lit::TRUE);
-    }
-
-    pub fn enforce_all(&mut self, bools: &[BAtom]) {
-        for &b in bools {
-            self.enforce(b);
-        }
+    /// Record that `b <=> literal`
+    pub fn bind<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, expr: Expr, literal: Lit) {
+        self.shape.expressions.bind(expr, literal);
     }
 
     /// Record that `b <=> literal`
-    pub fn bind(&mut self, b: BAtom, literal: Lit) {
-        match b {
-            BAtom::Cst(true) => self.make_literals_eq(literal, Lit::TRUE),
-            BAtom::Cst(false) => self.make_literals_eq(literal, Lit::FALSE),
-            BAtom::Literal(b) => self.make_literals_eq(literal, b),
-            BAtom::Expr(e) => {
-                let BExpr { expr: handle, negated } = e;
-                if let Some(prev) = self.shape.interned_expr(handle) {
-                    let prev = if negated { !prev } else { prev };
-                    self.make_literals_eq(prev, literal);
-                } else if negated {
-                    self.record_binding(handle, !literal);
-                } else {
-                    self.record_binding(handle, literal);
-                }
-            }
-        }
+    pub fn bind_literals(&mut self, l1: Lit, l2: Lit) {
+        self.shape.expressions.bind_literals(l1, l2);
     }
 
     // =========== Formatting ==============
@@ -278,6 +240,17 @@ impl Model {
         let atom = atom.into();
         crate::extensions::fmt(atom, self)
     }
+
+    pub fn print_state(&self) {
+        for v in self.state.variables() {
+            print!("{:?} <- {:?}", v, self.state.domain(v));
+            if let Some(lbl) = self.get_label(v) {
+                println!("    {:?}", lbl);
+            } else {
+                println!()
+            }
+        }
+    }
 }
 
 /// Identifies an external writer to the model.
@@ -293,15 +266,13 @@ impl WriterId {
     }
 }
 
-/// Provides write access to a model, making sure the built-in `WriterId` is always set.
-
-impl Default for Model {
+impl<Lbl: Label> Default for Model<Lbl> {
     fn default() -> Self {
         Self::new()
     }
 }
 
-impl Backtrack for Model {
+impl<Lbl> Backtrack for Model<Lbl> {
     fn save_state(&mut self) -> DecLvl {
         self.state.save_state()
     }
@@ -319,41 +290,13 @@ impl Backtrack for Model {
     }
 }
 
-impl ExpressionFactoryExt for Model {
-    fn intern_bool(&mut self, expr: Expr) -> BExpr {
-        self.shape.intern_boolean_expression(expr)
-    }
-
-    fn presence_literal(&self, variable: VarRef) -> Lit {
-        self.state.presence(variable)
-    }
-}
-
-impl AssignmentExt for Model {
-    fn symbols(&self) -> &SymbolTable {
-        &self.shape.symbols
-    }
-
+impl<Lbl> AssignmentExt for Model<Lbl> {
     fn entails(&self, literal: Lit) -> bool {
         self.state.entails(literal)
     }
 
-    fn literal_of_expr(&self, expr: BExpr) -> Option<Lit> {
-        match self.shape.bindings.as_lit(expr.expr) {
-            Some(l) => {
-                if expr.negated {
-                    Some(!l)
-                } else {
-                    Some(l)
-                }
-            }
-            None => None,
-        }
-    }
-
-    fn var_domain(&self, var: impl Into<VarRef>) -> IntDomain {
-        let (lb, ub) = self.state.bounds(var.into());
-        IntDomain { lb, ub }
+    fn var_domain(&self, var: impl Into<IAtom>) -> IntDomain {
+        self.state.var_domain(var)
     }
 
     fn presence_literal(&self, variable: VarRef) -> Lit {
@@ -365,8 +308,8 @@ impl AssignmentExt for Model {
     }
 }
 
-impl Shaped for Model {
-    fn get_shape(&self) -> &ModelShape {
+impl<Lbl: Label> Shaped<Lbl> for Model<Lbl> {
+    fn get_shape(&self) -> &ModelShape<Lbl> {
         &self.shape
     }
 }
diff --git a/model/src/state/domain.rs b/model/src/state/domain.rs
index f3aecd59..10db715e 100644
--- a/model/src/state/domain.rs
+++ b/model/src/state/domain.rs
@@ -1,4 +1,5 @@
 use crate::lang::IntCst;
+use std::fmt::Formatter;
 
 #[derive(Clone, Ord, PartialOrd, Eq, PartialEq, Debug)]
 pub struct IntDomain {
@@ -11,7 +12,7 @@ impl IntDomain {
     }
 
     pub fn size(&self) -> i64 {
-        (self.ub as i64) - (self.lb as i64)
+        (self.ub as i64) - (self.lb as i64) + 1
     }
 
     pub fn is_bound(&self) -> bool {
@@ -36,7 +37,7 @@ impl std::fmt::Display for IntDomain {
 }
 
 /// Represents the domain of an optional variable
-#[derive(Eq, PartialEq, Copy, Clone, Debug)]
+#[derive(Eq, PartialEq, Copy, Clone)]
 pub enum OptDomain {
     /// The variable is necessarily present and must take a value in the given bounds.
     Present(IntCst, IntCst),
@@ -45,3 +46,15 @@ pub enum OptDomain {
     /// The variable is known to be absent
     Absent,
 }
+
+impl std::fmt::Debug for OptDomain {
+    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
+        match *self {
+            OptDomain::Present(lb, ub) if lb == ub => write!(f, "[{}]", lb),
+            OptDomain::Present(lb, ub) => write!(f, "[{}, {}]", lb, ub),
+            OptDomain::Unknown(lb, ub) if lb == ub => write!(f, "?[{}]", lb),
+            OptDomain::Unknown(lb, ub) => write!(f, "?[{}, {}]", lb, ub),
+            OptDomain::Absent => write!(f, "_"),
+        }
+    }
+}
diff --git a/planners/src/bin/lcp.rs b/planners/src/bin/lcp.rs
index 35298297..4279454f 100644
--- a/planners/src/bin/lcp.rs
+++ b/planners/src/bin/lcp.rs
@@ -1,18 +1,17 @@
+use crate::Strat::{Activity, Forward};
 use anyhow::*;
 use aries_model::extensions::SavedAssignment;
 use aries_planners::encode::{encode, populate_with_task_network, populate_with_template_instances};
 use aries_planners::fmt::{format_hddl_plan, format_partial_plan, format_pddl_plan};
+use aries_planners::forward_search::ForwardSearcher;
+use aries_planners::Solver;
+use aries_planning::chronicles::analysis::hierarchical_is_non_recursive;
 use aries_planning::chronicles::*;
 use aries_planning::parsing::pddl::{find_domain_of, parse_pddl_domain, parse_pddl_problem, PddlFeature};
 use aries_planning::parsing::pddl_to_chronicles;
-use aries_solver::solver::Solver;
+use aries_solver::parallel_solver::ParSolver;
 use aries_tnet::theory::{StnConfig, StnTheory, TheoryPropagationLevel};
 use aries_utils::input::Input;
-
-use crate::Strat::{Activity, Forward};
-use aries_planners::forward_search::ForwardSearcher;
-use aries_planning::chronicles::analysis::hierarchical_is_non_recursive;
-use aries_solver::parallel_solver::ParSolver;
 use std::fs::File;
 use std::io::Write;
 use std::path::PathBuf;
diff --git a/planners/src/encode.rs b/planners/src/encode.rs
index cb9be7a0..750746b9 100644
--- a/planners/src/encode.rs
+++ b/planners/src/encode.rs
@@ -2,12 +2,13 @@
 //! into a combinatorial problem from Aries core.
 
 use crate::encoding::{conditions, effects, refinements_of, refinements_of_task, TaskRef, HORIZON, ORIGIN};
+use crate::Model;
 use anyhow::*;
 use aries_model::bounds::Lit;
-use aries_model::extensions::{AssignmentExt, ExpressionFactoryExt};
-use aries_model::lang::{BAtom, VarRef};
+use aries_model::extensions::AssignmentExt;
+use aries_model::lang::expr::*;
+use aries_model::lang::VarRef;
 use aries_model::lang::{IAtom, Variable};
-use aries_model::Model;
 use aries_planning::chronicles::constraints::ConstraintType;
 use aries_planning::chronicles::{
     ChronicleInstance, ChronicleOrigin, ChronicleTemplate, FiniteProblem, InvalidSubstitution, Problem, Sub,
@@ -187,7 +188,7 @@ pub fn populate_with_task_network(pb: &mut FiniteProblem, spec: &Problem, max_de
                     // We can thus unify the presence, start, end and parameters of   subtask/task pair.
                     // Unification is a best effort and might not succeed due to syntactical difference.
                     // We ignore any failed unification and let normal instantiation run its course.
-                    let _ = sub.add_bool_expr_unification(template.chronicle.presence.into(), task.scope.into());
+                    let _ = sub.add_bool_expr_unification(template.chronicle.presence, task.scope);
                     let _ = sub.add_int_expr_unification(template.chronicle.start, task.start);
                     let _ = sub.add_int_expr_unification(template.chronicle.end, task.end);
 
@@ -221,7 +222,7 @@ pub fn populate_with_task_network(pb: &mut FiniteProblem, spec: &Problem, max_de
     Ok(())
 }
 
-fn add_decomposition_constraints(pb: &FiniteProblem, model: &mut Model, constraints: &mut Vec<BAtom>) {
+fn add_decomposition_constraints(pb: &FiniteProblem, model: &mut Model) {
     for (instance_id, chronicle) in pb.chronicles.iter().enumerate() {
         for (task_id, task) in chronicle.chronicle.subtasks.iter().enumerate() {
             let subtask = TaskRef {
@@ -231,25 +232,25 @@ fn add_decomposition_constraints(pb: &FiniteProblem, model: &mut Model, constrai
                 task: &task.task_name,
             };
             let refiners = refinements_of(instance_id, task_id, pb);
-            enforce_refinement(subtask, refiners, model, constraints);
+            enforce_refinement(subtask, refiners, model);
         }
     }
 }
 
-fn enforce_refinement(t: TaskRef, supporters: Vec<TaskRef>, model: &mut Model, constraints: &mut Vec<BAtom>) {
+fn enforce_refinement(t: TaskRef, supporters: Vec<TaskRef>, model: &mut Model) {
     // if t is present then at least one supporter is present
-    let mut clause: Vec<BAtom> = Vec::with_capacity(supporters.len() + 1);
-    clause.push((!t.presence).into());
+    let mut clause: Vec<Lit> = Vec::with_capacity(supporters.len() + 1);
+    clause.push(!t.presence);
     for s in &supporters {
-        clause.push(s.presence.into());
+        clause.push(s.presence);
     }
-    constraints.push(model.or(&clause));
+    model.enforce(or(clause));
 
     // if a supporter is present, then all others are absent
     for (i, s1) in supporters.iter().enumerate() {
         for (j, s2) in supporters.iter().enumerate() {
             if i != j {
-                constraints.push(model.implies(s1.presence, !s2.presence));
+                model.enforce(implies(s1.presence, !s2.presence));
             }
         }
     }
@@ -260,23 +261,18 @@ fn enforce_refinement(t: TaskRef, supporters: Vec<TaskRef>, model: &mut Model, c
         assert!(model
             .state
             .only_present_with(s.presence.variable(), t.presence.variable()));
-        constraints.push(model.implies(s.presence, t.presence)); // TODO: can we get rid of this
+        model.enforce(implies(s.presence, t.presence)); // TODO: can we get rid of this
 
-        constraints.push(model.opt_eq(s.start, t.start));
-        constraints.push(model.opt_eq(s.end, t.end));
+        model.enforce(opt_eq(s.start, t.start));
+        model.enforce(opt_eq(s.end, t.end));
         assert_eq!(s.task.len(), t.task.len());
         for (a, b) in s.task.iter().zip(t.task.iter()) {
-            constraints.push(model.opt_eq(*a, *b))
+            model.enforce(opt_eq(*a, *b))
         }
     }
 }
 
-fn add_symmetry_breaking(
-    pb: &FiniteProblem,
-    model: &mut Model,
-    constraints: &mut Vec<BAtom>,
-    tpe: SymmetryBreakingType,
-) {
+fn add_symmetry_breaking(pb: &FiniteProblem, model: &mut Model, tpe: SymmetryBreakingType) {
     match tpe {
         SymmetryBreakingType::None => {}
         SymmetryBreakingType::Simple => {
@@ -292,8 +288,8 @@ fn add_symmetry_breaking(
             for (instance1, template_id1, generation_id1) in chronicles() {
                 for (instance2, template_id2, generation_id2) in chronicles() {
                     if template_id1 == template_id2 && generation_id1 < generation_id2 {
-                        constraints.push(model.implies(instance1.chronicle.presence, instance2.chronicle.presence));
-                        constraints.push(model.leq(instance1.chronicle.start, instance2.chronicle.start))
+                        model.enforce(implies(instance1.chronicle.presence, instance2.chronicle.presence));
+                        model.enforce(leq(instance1.chronicle.start, instance2.chronicle.start));
                     }
                 }
             }
@@ -305,23 +301,20 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
     let mut model = pb.model.clone();
     let symmetry_breaking_tpe = SYMMETRY_BREAKING.get();
 
-    // the set of constraints that should be enforced
-    let mut constraints: Vec<BAtom> = Vec::new();
-
     let effs: Vec<_> = effects(pb).collect();
     let conds: Vec<_> = conditions(pb).collect();
     let eff_ends: Vec<_> = effs.iter().map(|_| model.new_ivar(ORIGIN, HORIZON, "")).collect();
 
     // for each condition, make sure the end is after the start
     for &(_prez_cond, cond) in &conds {
-        constraints.push(model.leq(cond.start, cond.end));
+        model.enforce(leq(cond.start, cond.end));
     }
 
     // for each effect, make sure the three time points are ordered
     for ieff in 0..effs.len() {
         let (_prez_eff, eff) = effs[ieff];
-        constraints.push(model.leq(eff.persistence_start, eff_ends[ieff]));
-        constraints.push(model.leq(eff.transition_start, eff.persistence_start))
+        model.enforce(leq(eff.persistence_start, eff_ends[ieff]));
+        model.enforce(leq(eff.transition_start, eff.persistence_start));
     }
 
     // are two state variables unifiable?
@@ -339,7 +332,7 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
     };
 
     // for each pair of effects, enforce coherence constraints
-    let mut clause: Vec<BAtom> = Vec::with_capacity(32);
+    let mut clause: Vec<Lit> = Vec::with_capacity(32);
     for (i, &(p1, e1)) in effs.iter().enumerate() {
         for j in i + 1..effs.len() {
             let &(p2, e2) = &effs[j];
@@ -350,8 +343,8 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
             }
 
             clause.clear();
-            clause.push((!p1).into());
-            clause.push((!p2).into());
+            clause.push(!p1);
+            clause.push(!p2);
             assert_eq!(e1.state_var.len(), e2.state_var.len());
             for idx in 0..e1.state_var.len() {
                 let a = e1.state_var[idx];
@@ -359,23 +352,23 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
                 // enforce different : a < b || a > b
                 // if they are the same variable, there is nothing we can do to separate them
                 if a != b {
-                    clause.push(model.neq(a, b));
+                    clause.push(model.reify(neq(a, b)));
                 }
             }
 
-            clause.push(model.leq(eff_ends[j], e1.transition_start));
-            clause.push(model.leq(eff_ends[i], e2.transition_start));
+            clause.push(model.reify(leq(eff_ends[j], e1.transition_start)));
+            clause.push(model.reify(leq(eff_ends[i], e2.transition_start)));
 
             // add coherence constraint
-            constraints.push(model.or(&clause));
+            model.enforce(or(clause.as_slice()));
         }
     }
 
     // support constraints
-    for (cond_id, &(prez_cond, cond)) in conds.iter().enumerate() {
-        let mut supported: Vec<BAtom> = Vec::with_capacity(128);
+    for (_cond_id, &(prez_cond, cond)) in conds.iter().enumerate() {
+        let mut supported: Vec<Lit> = Vec::with_capacity(128);
         // no need to support if the condition is not present
-        supported.push((!prez_cond).into());
+        supported.push(!prez_cond);
 
         for (eff_id, &(prez_eff, eff)) in effs.iter().enumerate() {
             // quick check that the condition and effect are not trivially incompatible
@@ -386,9 +379,9 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
                 continue;
             }
             // vector to store the AND clause
-            let mut supported_by_eff_conjunction: Vec<BAtom> = Vec::with_capacity(32);
+            let mut supported_by_eff_conjunction: Vec<Lit> = Vec::with_capacity(32);
             // support only possible if the effect is present
-            supported_by_eff_conjunction.push(prez_eff.into());
+            supported_by_eff_conjunction.push(prez_eff);
 
             assert_eq!(cond.state_var.len(), eff.state_var.len());
             // same state variable
@@ -396,29 +389,25 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
                 let a = cond.state_var[idx];
                 let b = eff.state_var[idx];
 
-                supported_by_eff_conjunction.push(model.eq(a, b));
+                supported_by_eff_conjunction.push(model.reify(eq(a, b)));
             }
             // same value
             let condition_value = cond.value;
             let effect_value = eff.value;
-            supported_by_eff_conjunction.push(model.eq(condition_value, effect_value));
+            supported_by_eff_conjunction.push(model.reify(eq(condition_value, effect_value)));
 
             // effect's persistence contains condition
-            supported_by_eff_conjunction.push(model.leq(eff.persistence_start, cond.start));
-            supported_by_eff_conjunction.push(model.leq(cond.end, eff_ends[eff_id]));
+            supported_by_eff_conjunction.push(model.reify(leq(eff.persistence_start, cond.start)));
+            supported_by_eff_conjunction.push(model.reify(leq(cond.end, eff_ends[eff_id])));
 
-            let support_lit = model
-                .new_optional_bvar(prez_cond, format!("support_{}_by_{}", cond_id, eff_id))
-                .true_lit();
-            let support_expr = model.and(&supported_by_eff_conjunction);
-            model.bind(support_expr, support_lit);
+            let support_lit = model.reify(and(supported_by_eff_conjunction));
 
             // add this support expression to the support clause
-            supported.push(support_lit.into());
+            supported.push(support_lit);
         }
 
         // enforce necessary conditions for condition' support
-        constraints.push(model.or(&supported));
+        model.enforce(or(supported));
     }
 
     // chronicle constraints
@@ -426,25 +415,26 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
         for constraint in &instance.chronicle.constraints {
             match constraint.tpe {
                 ConstraintType::InTable { table_id } => {
-                    let mut supported_by_a_line: Vec<BAtom> = Vec::with_capacity(256);
-                    supported_by_a_line.push((!instance.chronicle.presence).into());
+                    let mut supported_by_a_line: Vec<Lit> = Vec::with_capacity(256);
+                    supported_by_a_line.push(!instance.chronicle.presence);
                     let vars = &constraint.variables;
                     for values in pb.tables[table_id as usize].lines() {
                         assert_eq!(vars.len(), values.len());
                         let mut supported_by_this_line = Vec::with_capacity(16);
                         for (&var, &val) in vars.iter().zip(values.iter()) {
-                            // TODO: using 2 LEQ might avoid the need for reification
-                            supported_by_this_line.push(model.eq(var, val));
+                            let var = var.int_view().unwrap();
+                            supported_by_this_line.push(model.reify(leq(var, val)));
+                            supported_by_this_line.push(model.reify(geq(var, val)));
                         }
-                        supported_by_a_line.push(model.and(&supported_by_this_line));
+                        supported_by_a_line.push(model.reify(and(supported_by_this_line)));
                     }
-                    constraints.push(model.or(&supported_by_a_line));
+                    model.enforce(or(supported_by_a_line));
                 }
                 ConstraintType::Lt => match constraint.variables.as_slice() {
                     &[a, b] => {
                         let a: IAtom = a.try_into()?;
                         let b: IAtom = b.try_into()?;
-                        constraints.push(model.lt(a, b))
+                        model.enforce(lt(a, b));
                     }
                     x => bail!("Invalid variable pattern for LT constraint: {:?}", x),
                 },
@@ -455,7 +445,7 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
                             constraint.variables.len()
                         );
                     }
-                    constraints.push(model.eq(constraint.variables[0], constraint.variables[1]));
+                    model.enforce(eq(constraint.variables[0], constraint.variables[1]));
                 }
                 ConstraintType::Neq => {
                     if constraint.variables.len() != 2 {
@@ -464,7 +454,7 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
                             constraint.variables.len()
                         );
                     }
-                    constraints.push(model.neq(constraint.variables[0], constraint.variables[1]));
+                    model.enforce(neq(constraint.variables[0], constraint.variables[1]));
                 }
             }
         }
@@ -472,20 +462,18 @@ pub fn encode(pb: &FiniteProblem) -> anyhow::Result<Model> {
 
     for ch in &pb.chronicles {
         // chronicle finishes before the horizon and has a non negative duration
-        constraints.push(model.opt_leq(ch.chronicle.end, pb.horizon));
-        constraints.push(model.opt_leq(ch.chronicle.start, ch.chronicle.end));
+        model.enforce(opt_leq(ch.chronicle.end, pb.horizon));
+        model.enforce(opt_leq(ch.chronicle.start, ch.chronicle.end));
 
         // enforce temporal coherence between the chronicle and its subtasks
         for subtask in &ch.chronicle.subtasks {
-            constraints.push(model.opt_leq(subtask.start, subtask.end));
-            constraints.push(model.opt_leq(ch.chronicle.start, subtask.start));
-            constraints.push(model.opt_leq(subtask.end, ch.chronicle.end));
+            model.enforce(opt_leq(subtask.start, subtask.end));
+            model.enforce(opt_leq(ch.chronicle.start, subtask.start));
+            model.enforce(opt_leq(subtask.end, ch.chronicle.end));
         }
     }
-    add_decomposition_constraints(pb, &mut model, &mut constraints);
-    add_symmetry_breaking(pb, &mut model, &mut constraints, symmetry_breaking_tpe);
-
-    model.enforce_all(&constraints);
+    add_decomposition_constraints(pb, &mut model);
+    add_symmetry_breaking(pb, &mut model, symmetry_breaking_tpe);
 
     Ok(model)
 }
diff --git a/planners/src/fmt.rs b/planners/src/fmt.rs
index 88384f1a..3c5854c9 100644
--- a/planners/src/fmt.rs
+++ b/planners/src/fmt.rs
@@ -3,10 +3,10 @@
 use anyhow::*;
 use std::fmt::Write;
 
+use crate::Model;
 use aries_model::extensions::{AssignmentExt, SavedAssignment, Shaped};
 use aries_model::lang::SAtom;
 use aries_model::symbols::SymId;
-use aries_model::Model;
 use aries_planning::chronicles::{ChronicleInstance, ChronicleKind, ChronicleOrigin, FiniteProblem, SubTask};
 
 pub fn format_partial_symbol(x: &SAtom, ass: &Model, out: &mut String) {
@@ -142,20 +142,20 @@ pub fn format_pddl_plan(problem: &FiniteProblem, ass: &SavedAssignment) -> Resul
     let mut out = String::new();
     let mut plan = Vec::new();
     for ch in &problem.chronicles {
-        if ass.boolean_value_of(ch.chronicle.presence) != Some(true) {
+        if ass.value(ch.chronicle.presence) != Some(true) {
             continue;
         }
         if ch.origin == ChronicleOrigin::Original {
             continue;
         }
-        let start = ass.domain_of(ch.chronicle.start).0;
+        let start = ass.var_domain(ch.chronicle.start).lb;
         let name: Vec<SymId> = ch
             .chronicle
             .name
             .iter()
             .map(|satom| ass.sym_domain_of(*satom).into_singleton().unwrap())
             .collect();
-        let name = ass.symbols().format(&name);
+        let name = problem.model.shape.symbols.format(&name);
         plan.push((start, name));
     }
 
@@ -172,14 +172,14 @@ pub fn format_hddl_plan(problem: &FiniteProblem, ass: &SavedAssignment) -> Resul
     writeln!(f, "==>")?;
     let fmt1 = |x: &SAtom| -> String {
         let sym = ass.sym_domain_of(*x).into_singleton().unwrap();
-        ass.symbols().symbol(sym).to_string()
+        problem.model.shape.symbols.symbol(sym).to_string()
     };
     let fmt = |name: &[SAtom]| -> String {
         let syms: Vec<_> = name
             .iter()
             .map(|x| ass.sym_domain_of(*x).into_singleton().unwrap())
             .collect();
-        ass.symbols().format(&syms)
+        problem.model.shape.symbols.format(&syms)
     };
     let mut chronicles: Vec<_> = problem
         .chronicles
diff --git a/planners/src/forward_search.rs b/planners/src/forward_search.rs
index 33ef5b02..c0e498b6 100644
--- a/planners/src/forward_search.rs
+++ b/planners/src/forward_search.rs
@@ -1,11 +1,11 @@
 //! A search controller that mimics forward search for HTN planning.
 
 use crate::encoding::refinements_of;
+use crate::{Model, Var};
 use aries_backtrack::{Backtrack, DecLvl};
 use aries_model::bounds::Lit;
 use aries_model::extensions::AssignmentExt;
 use aries_model::lang::{Atom, IVar, VarRef};
-use aries_model::Model;
 use aries_planning::chronicles::{ChronicleInstance, FiniteProblem, SubTask};
 use aries_solver::solver::search::{Decision, SearchControl};
 use aries_solver::solver::stats::Stats;
@@ -122,7 +122,7 @@ impl ForwardSearcher {
     }
 }
 
-impl SearchControl for ForwardSearcher {
+impl SearchControl<Var> for ForwardSearcher {
     fn next_decision(&mut self, _stats: &Stats, model: &Model) -> Option<Decision> {
         let xx = earliest_pending_chronicle(&self.problem, model);
         let yy = earliest_pending_task(&self.problem, model);
@@ -153,7 +153,7 @@ impl SearchControl for ForwardSearcher {
         res.map(Decision::SetLiteral)
     }
 
-    fn clone_to_box(&self) -> Box<dyn SearchControl + Send> {
+    fn clone_to_box(&self) -> Box<dyn SearchControl<Var> + Send> {
         Box::new(self.clone())
     }
 }
diff --git a/planners/src/lib.rs b/planners/src/lib.rs
index a4a27101..530f6e2d 100644
--- a/planners/src/lib.rs
+++ b/planners/src/lib.rs
@@ -2,3 +2,10 @@ pub mod encode;
 pub mod encoding;
 pub mod fmt;
 pub mod forward_search;
+
+/// Label of a variable
+pub type Var = String;
+
+pub type Model = aries_model::Model<Var>;
+pub type Solver = aries_solver::solver::Solver<Var>;
+pub type ParSolver = aries_solver::parallel_solver::ParSolver<Var>;
diff --git a/planning/src/chronicles/concrete.rs b/planning/src/chronicles/concrete.rs
index 29e41463..ebbb92dc 100644
--- a/planning/src/chronicles/concrete.rs
+++ b/planning/src/chronicles/concrete.rs
@@ -25,24 +25,14 @@ pub trait Substitution {
 
     fn sub(&self, atom: Atom) -> Atom {
         match atom {
-            Atom::Bool(b) => self.bsub(b).into(),
+            Atom::Bool(b) => self.sub_lit(b).into(),
             Atom::Int(i) => self.isub(i).into(),
             Atom::Sym(s) => self.ssub(s).into(),
         }
     }
 
     fn isub(&self, i: IAtom) -> IAtom {
-        match i.var {
-            Some(x) => IAtom::new(Some(self.sub_ivar(x)), i.shift),
-            None => IAtom::new(None, i.shift),
-        }
-    }
-    fn bsub(&self, b: BAtom) -> BAtom {
-        match b {
-            BAtom::Cst(b) => BAtom::Cst(b),
-            BAtom::Literal(b) => BAtom::Literal(self.sub_lit(b)),
-            BAtom::Expr(_) => panic!("UNSUPPORTED substitution in an expression"),
-        }
+        IAtom::new(self.sub_ivar(i.var), i.shift)
     }
 
     fn ssub(&self, s: SAtom) -> SAtom {
@@ -113,33 +103,29 @@ impl Sub {
     }
     pub fn add_int_expr_unification(&mut self, param: IAtom, instance: IAtom) -> Result<(), InvalidSubstitution> {
         match (param, instance) {
-            (
-                IAtom {
-                    var: Some(x),
-                    shift: dx,
-                },
-                IAtom {
-                    var: Some(y),
-                    shift: dy,
-                },
-            ) if dx == dy => self.add(x.into(), y.into()),
-            (IAtom { var: None, shift: dx }, IAtom { var: None, shift: dy }) if dx == dy => Ok(()),
+            (IAtom { var: x, shift: dx }, IAtom { var: y, shift: dy }) if dx == dy => {
+                if x == y {
+                    Ok(())
+                } else {
+                    self.add(x.into(), y.into())
+                }
+            }
             _ => Err(InvalidSubstitution::IncompatibleStructures(
                 param.into(),
                 instance.into(),
             )),
         }
     }
-    pub fn add_bool_expr_unification(&mut self, param: BAtom, instance: BAtom) -> Result<(), InvalidSubstitution> {
-        match (param, instance) {
-            (BAtom::Literal(l), BAtom::Literal(l2)) if l.relation() == l2.relation() && l.value() == l2.value() => {
-                self.add_untyped(l.variable(), l2.variable())
-            }
-            (BAtom::Cst(a), BAtom::Cst(b)) if a == b => Ok(()),
-            _ => Err(InvalidSubstitution::IncompatibleStructures(
+    pub fn add_bool_expr_unification(&mut self, param: Lit, instance: Lit) -> Result<(), InvalidSubstitution> {
+        if param == instance {
+            Ok(())
+        } else if param.relation() == instance.relation() && param.value() == instance.value() {
+            self.add_untyped(param.variable(), instance.variable())
+        } else {
+            Err(InvalidSubstitution::IncompatibleStructures(
                 param.into(),
                 instance.into(),
-            )),
+            ))
         }
     }
 
diff --git a/planning/src/chronicles/mod.rs b/planning/src/chronicles/mod.rs
index d495cc6a..fd982fea 100644
--- a/planning/src/chronicles/mod.rs
+++ b/planning/src/chronicles/mod.rs
@@ -47,7 +47,7 @@ impl StateFun {
 
 #[derive(Clone)]
 pub struct Ctx {
-    pub model: Model,
+    pub model: Model<String>,
     pub state_functions: Vec<StateFun>,
     origin: IAtom,
     horizon: IAtom,
@@ -169,7 +169,7 @@ pub struct Problem {
 
 #[derive(Clone)]
 pub struct FiniteProblem {
-    pub model: Model,
+    pub model: Model<String>,
     pub origin: IAtom,
     pub horizon: IAtom,
     pub chronicles: Vec<ChronicleInstance>,
diff --git a/planning/src/parsing/mod.rs b/planning/src/parsing/mod.rs
index 0ce04270..385b2ef4 100644
--- a/planning/src/parsing/mod.rs
+++ b/planning/src/parsing/mod.rs
@@ -563,7 +563,7 @@ fn read_conjunction_impl(e: &SExpr, t: &impl Fn(&sexpr::SAtom) -> Result<SAtom>,
         let TermLoc(t, _) = read_term(to_negate, &t)?;
         let negated = match t {
             Term::Binding(sv, value) => {
-                if let Ok(value) = BAtom::try_from(value) {
+                if let Ok(value) = Lit::try_from(value) {
                     Term::Binding(sv, Atom::from(!value))
                 } else {
                     return Err(to_negate.invalid("Could not apply 'not' to this expression").into());
diff --git a/solver/src/lib.rs b/solver/src/lib.rs
index d4c0545a..98e08431 100644
--- a/solver/src/lib.rs
+++ b/solver/src/lib.rs
@@ -1,3 +1,5 @@
+extern crate aries_model;
+
 pub mod clauses;
 pub(crate) mod cpu_time;
 pub mod parallel_solver;
@@ -7,36 +9,33 @@ pub mod theories;
 
 use crate::solver::BindingResult;
 use aries_backtrack::{Backtrack, DecLvl};
-use aries_model::{Model, WriterId};
-
 use aries_model::bounds::Lit;
-use aries_model::expressions::ExprHandle;
-use aries_model::lang::Expr;
+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`.
     ///
     /// The return value should provide feedback on whether it succeeded or failed to do so.
-    fn bind(&mut self, literal: Lit, expr: ExprHandle, i: &mut Model) -> BindingResult;
+    fn bind(&mut self, literal: Lit, expr: &Expr, doms: &mut Domains) -> BindingResult;
 }
 
 /// A convenience trait that when implemented  will allow deriving the [Bind] trait.
 pub trait BindSplit {
-    fn enforce_true(&mut self, expr: &Expr, model: &mut Model) -> BindingResult;
-    fn enforce_false(&mut self, expr: &Expr, model: &mut Model) -> BindingResult;
-    fn enforce_eq(&mut self, literal: Lit, expr: &Expr, model: &mut Model) -> BindingResult;
+    fn enforce_true(&mut self, expr: &Expr, doms: &mut Domains) -> BindingResult;
+    fn enforce_false(&mut self, expr: &Expr, doms: &mut Domains) -> BindingResult;
+    fn enforce_eq(&mut self, literal: Lit, expr: &Expr, doms: &mut Domains) -> BindingResult;
 }
 
 impl<T: BindSplit> Bind for T {
-    fn bind(&mut self, literal: Lit, expr: ExprHandle, i: &mut Model) -> BindingResult {
-        debug_assert_eq!(i.state.current_decision_level(), DecLvl::ROOT);
-        let expr = i.shape.expressions.get(expr);
-        match i.state.value(literal) {
-            Some(true) => self.enforce_true(&expr, i),
-            Some(false) => self.enforce_false(&expr, i),
-            None => self.enforce_eq(literal, &expr, i),
+    fn bind(&mut self, literal: Lit, expr: &Expr, doms: &mut Domains) -> BindingResult {
+        debug_assert_eq!(doms.current_decision_level(), DecLvl::ROOT);
+        match doms.value(literal) {
+            Some(true) => self.enforce_true(expr, doms),
+            Some(false) => self.enforce_false(expr, doms),
+            None => self.enforce_eq(literal, expr, doms),
         }
     }
 }
diff --git a/solver/src/parallel_solver.rs b/solver/src/parallel_solver.rs
index 47af6bc5..31aceccf 100644
--- a/solver/src/parallel_solver.rs
+++ b/solver/src/parallel_solver.rs
@@ -1,24 +1,26 @@
 use crate::signals::{InputSignal, InputStream, OutputSignal, SolverOutput, ThreadID};
 use crate::solver::{Exit, Solver};
-use aries_model::extensions::SavedAssignment;
+use aries_model::extensions::{SavedAssignment, Shaped};
 use aries_model::lang::{IAtom, IntCst};
+use aries_model::{Label, ModelShape};
 use std::sync::mpsc::{channel, Receiver, Sender};
 use std::sync::Arc;
 use std::thread;
 
-pub struct ParSolver {
-    solvers: Vec<Worker>,
+pub struct ParSolver<Lbl> {
+    base_model: ModelShape<Lbl>,
+    solvers: Vec<Worker<Lbl>>,
 }
 
 /// A worker is a solver that is either running on another thread or available.
 /// If it is running we only have its input stream to commonuicate with it.
-enum Worker {
+enum Worker<Lbl> {
     Running(InputStream),
-    Idle(Box<Solver>),
+    Idle(Box<Solver<Lbl>>),
 }
-impl Worker {
+impl<Lbl: Label> Worker<Lbl> {
     /// Marks the solver a running and return it if it wasn't previously running.
-    pub fn extract(&mut self) -> Option<Box<Solver>> {
+    pub fn extract(&mut self) -> Option<Box<Solver<Lbl>>> {
         let stream = match self {
             Worker::Running(_) => return None,
             Worker::Idle(solver) => solver.input_stream(),
@@ -35,19 +37,20 @@ impl Worker {
 
 /// Result of running a computation with a result of type `O` on a solver.
 /// The solver it self is provided as a part of the result.
-struct WorkerResult<O> {
+struct WorkerResult<O, Lbl> {
     id: ThreadID,
     output: Result<O, Exit>,
-    solver: Box<Solver>,
+    solver: Box<Solver<Lbl>>,
 }
 
-impl ParSolver {
+impl<Lbl: Label> ParSolver<Lbl> {
     /// Creates a new parallel solver.
     ///
     /// All solvers will be based on a clone of `base_solver`, on which the provided `adapt` function
     /// will be called to allow its customisation.
-    pub fn new(mut base_solver: Box<Solver>, num_workers: usize, adapt: impl Fn(usize, &mut Solver)) -> Self {
+    pub fn new(mut base_solver: Box<Solver<Lbl>>, num_workers: usize, adapt: impl Fn(usize, &mut Solver<Lbl>)) -> Self {
         let mut solver = ParSolver {
+            base_model: base_solver.model.shape.clone(),
             solvers: Vec::with_capacity(num_workers),
         };
         for i in 0..(num_workers - 1) {
@@ -93,13 +96,13 @@ impl ParSolver {
     fn race_solvers<O, F>(&mut self, run: F) -> Result<O, Exit>
     where
         O: Send + 'static,
-        F: Fn(&mut Solver) -> Result<O, Exit> + Send + 'static + Copy,
+        F: Fn(&mut Solver<Lbl>) -> Result<O, Exit> + Send + 'static + Copy,
     {
         let solvers_output = self.plug_solvers_output();
         let (result_snd, result_rcv) = channel();
 
         // lambda used to start a thread and run a solver on it.
-        let spawn = |id: usize, mut solver: Box<Solver>, result_snd: Sender<WorkerResult<O>>| {
+        let spawn = |id: usize, mut solver: Box<Solver<Lbl>>, result_snd: Sender<WorkerResult<O, Lbl>>| {
             thread::spawn(move || {
                 let output = run(&mut solver);
                 let answer = WorkerResult { id, output, solver };
@@ -173,3 +176,9 @@ impl ParSolver {
         }
     }
 }
+
+impl<Lbl: Label> Shaped<Lbl> for ParSolver<Lbl> {
+    fn get_shape(&self) -> &ModelShape<Lbl> {
+        &self.base_model
+    }
+}
diff --git a/solver/src/solver.rs b/solver/src/solver.rs
index 353b0a3d..535531d6 100644
--- a/solver/src/solver.rs
+++ b/solver/src/solver.rs
@@ -12,12 +12,14 @@ use crate::solver::stats::Stats;
 use crate::solver::theory_solver::TheorySolver;
 use crate::{Bind, Contradiction, Theory};
 use aries_backtrack::{Backtrack, DecLvl};
-use aries_model::bindings::{BindTarget, BindingCursor};
 use aries_model::bounds::{Disjunction, Lit};
-use aries_model::extensions::{AssignmentExt, DisjunctionExt, SavedAssignment};
-use aries_model::lang::{BAtom, IAtom, IntCst};
+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::{Model, WriterId};
+use aries_model::{Label, Model, ModelShape, WriterId};
 use std::fmt::Formatter;
 use std::sync::mpsc::Sender;
 use std::sync::Arc;
@@ -82,38 +84,38 @@ impl std::fmt::Display for Exit {
 }
 impl std::error::Error for Exit {}
 
-pub struct Solver {
-    pub model: Model,
-    pub brancher: Box<dyn SearchControl + Send>,
+pub struct Solver<Lbl> {
+    pub model: Model<Lbl>,
+    constraints: Constraints<Lbl>,
+    pub brancher: Box<dyn SearchControl<Lbl> + Send>,
     reasoners: Reasoners,
     decision_level: DecLvl,
     pub stats: Stats,
     /// A data structure with the various communication channels
     /// need to receive/sent updates and commands.
     sync: Synchro,
-    next_binding: BindingCursor,
     /// A queue of literals that we know to be tautologies but that have not been propagated yet.
     /// Invariant: if the queue is non-empty, we are at root level.
     pending_tautologies: Vec<Lit>,
 }
-impl Solver {
-    pub fn new(mut model: Model) -> Solver {
+impl<Lbl: Label> Solver<Lbl> {
+    pub fn new(mut model: Model<Lbl>) -> Solver<Lbl> {
         let sat_id = model.shape.new_write_token();
         let sat = SatSolver::new(sat_id);
 
         Solver {
             model,
+            constraints: Constraints::default(),
             brancher: default_brancher(),
             reasoners: Reasoners::new(sat, sat_id),
             decision_level: DecLvl::ROOT,
             stats: Default::default(),
             sync: Synchro::new(),
-            next_binding: BindingCursor::first(),
             pending_tautologies: vec![],
         }
     }
 
-    pub fn set_brancher(&mut self, brancher: impl SearchControl + 'static + Send) {
+    pub fn set_brancher(&mut self, brancher: impl SearchControl<Lbl> + 'static + Send) {
         self.brancher = Box::new(brancher)
     }
 
@@ -143,28 +145,30 @@ impl Solver {
         self.pending_tautologies.push(lit);
     }
 
-    pub fn enforce(&mut self, bool_expr: impl Into<BAtom>) {
+    pub fn enforce<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, bool_expr: Expr) {
         assert_eq!(self.decision_level, DecLvl::ROOT);
-        self.model.enforce(bool_expr.into());
-        self.process_bindings();
+        self.model.enforce(bool_expr);
+        self.post_constraints();
     }
-
-    pub fn enforce_all(&mut self, bool_exprs: &[BAtom]) {
+    pub fn enforce_all<Expr: Normalize<T>, T: ReifiableExpr>(&mut self, bools: impl IntoIterator<Item = Expr>) {
         assert_eq!(self.decision_level, DecLvl::ROOT);
-        self.model.enforce_all(bool_exprs);
-        self.process_bindings();
+        self.model.enforce_all(bools);
+        self.post_constraints();
     }
 
     // TODO: we should clean the call places: it should be invoked as early as possible but after all reasoners are added
-    fn process_bindings(&mut self) {
+    pub fn post_constraints(&mut self) {
+        self.constraints.decompose_all(&mut self.model);
+
         use BindingResult::*;
         let start_time = Instant::now();
         let start_cycles = StartCycleCount::now();
 
-        while let Some((llit, expr)) = self.model.shape.bindings.pop_next_event(&mut self.next_binding) {
+        while let Some((llit, expr)) = self.constraints.pop_next_constraint() {
+            let llit = *llit;
             assert_eq!(self.model.current_decision_level(), DecLvl::ROOT);
             match expr {
-                BindTarget::Literal(rlit) => {
+                &BindTarget::Literal(rlit) => {
                     if self.model.entails(llit) {
                         self.set_tautology(rlit);
                     } else if self.model.entails(!llit) {
@@ -181,26 +185,23 @@ impl Solver {
                     }
                 }
                 BindTarget::Expr(expr) => {
+                    let expr = expr.clone();
                     // while let Some(binding) = reader.pop(&queue).copied() {
                     let mut supported = false;
 
                     // expr <=> lit_of_expr
-                    match self.reasoners.sat.bind(llit, expr, &mut self.model) {
-                        Enforced | Refined => supported = true,
+                    match self.reasoners.sat.bind(llit, expr.as_ref(), &mut self.model.state) {
+                        Enforced => supported = true,
                         Unsupported => {}
                     }
                     for theory in &mut self.reasoners.theories {
-                        match theory.bind(llit, expr, &mut self.model) {
-                            Enforced | Refined => supported = true,
+                        match theory.bind(llit, expr.as_ref(), &mut self.model.state) {
+                            Enforced => supported = true,
                             Unsupported => {}
                         }
                     }
                     if !supported {
-                        panic!(
-                            "Unsupported binding: {:?}  {:?}",
-                            self.model.shape.expressions.get_ref(expr),
-                            self.reasoners.theories.len()
-                        );
+                        panic!("Unsupported binding: {:?} <=> {:?}", llit, expr);
                     }
                 }
             }
@@ -212,9 +213,9 @@ impl Solver {
     /// Searches for the first satisfying assignment, returning none if the search
     /// space was exhausted without encountering a solution.
     pub fn solve(&mut self) -> Result<Option<Arc<SavedAssignment>>, Exit> {
-        self.process_bindings();
+        self.post_constraints();
         match self._solve()? {
-            SolveResult::AtSolution => Ok(Some(Arc::new(self.model.clone()))),
+            SolveResult::AtSolution => Ok(Some(Arc::new(self.model.state.clone()))),
             SolveResult::ExternalSolution(s) => Ok(Some(s)),
             SolveResult::Unsat => Ok(None),
         }
@@ -280,7 +281,7 @@ impl Solver {
         objective: impl Into<IAtom>,
         mut on_new_solution: impl FnMut(IntCst, &SavedAssignment),
     ) -> Result<Option<(IntCst, Arc<SavedAssignment>)>, Exit> {
-        self.process_bindings();
+        self.post_constraints();
         let objective = objective.into();
         // best solution found so far
         let mut best = None;
@@ -288,10 +289,10 @@ impl Solver {
             let sol = match self._solve()? {
                 SolveResult::AtSolution => {
                     // solver stopped at a solution, this is necessarily an improvement on the best solution found so far
-                    let sol = Arc::new(self.model.clone());
+                    let sol = Arc::new(self.model.state.clone());
                     // notify other solvers that we have found a new solution
                     self.sync.notify_solution_found(sol.clone());
-                    let lb = sol.domain_of(objective).0;
+                    let lb = sol.var_domain(objective).lb;
                     on_new_solution(lb, &sol);
                     sol
                 }
@@ -300,7 +301,7 @@ impl Solver {
             };
 
             // determine whether the solution found is an improvement on the previous one (might not be the case if sent by another solver)
-            let lb = sol.domain_of(objective).0;
+            let lb = sol.var_domain(objective).lb;
             let is_improvement = match best {
                 None => true,
                 Some((previous_best, _)) => lb < previous_best,
@@ -428,7 +429,7 @@ impl Solver {
     /// - `Err(clause)`: if a conflict was found. In this case, `clause` is a conflicting cause in the current
     ///   decision level that   
     pub fn propagate(&mut self) -> Result<(), Disjunction> {
-        self.process_bindings();
+        self.post_constraints();
         let global_start = StartCycleCount::now();
         while let Some(lit) = self.pending_tautologies.pop() {
             debug_assert_eq!(self.current_decision_level(), DecLvl::ROOT);
@@ -512,7 +513,7 @@ impl Solver {
     }
 }
 
-impl Backtrack for Solver {
+impl<Lbl> Backtrack for Solver<Lbl> {
     fn save_state(&mut self) -> DecLvl {
         self.decision_level += 1;
         let n = self.decision_level;
@@ -557,23 +558,28 @@ impl Backtrack for Solver {
     }
 }
 
-impl Clone for Solver {
+impl<Lbl: Label> Clone for Solver<Lbl> {
     fn clone(&self) -> Self {
         Solver {
             model: self.model.clone(),
+            constraints: self.constraints.clone(),
             brancher: self.brancher.clone_to_box(),
             reasoners: self.reasoners.clone(),
             decision_level: self.decision_level,
             stats: self.stats.clone(),
             sync: self.sync.clone(),
-            next_binding: self.next_binding,
             pending_tautologies: self.pending_tautologies.clone(),
         }
     }
 }
 
+impl<Lbl: Label> Shaped<Lbl> for Solver<Lbl> {
+    fn get_shape(&self) -> &ModelShape<Lbl> {
+        self.model.get_shape()
+    }
+}
+
 pub enum BindingResult {
     Enforced,
     Unsupported,
-    Refined,
 }
diff --git a/solver/src/solver/sat_solver.rs b/solver/src/solver/sat_solver.rs
index fa218636..73a8c5ff 100644
--- a/solver/src/solver/sat_solver.rs
+++ b/solver/src/solver/sat_solver.rs
@@ -5,12 +5,11 @@ use aries_backtrack::{Backtrack, DecLvl, ObsTrailCursor, Trail};
 use aries_collections::set::RefSet;
 use aries_model::bounds::{Disjunction, Lit, WatchSet, Watches};
 use aries_model::extensions::DisjunctionExt;
-use aries_model::lang::*;
+use aries_model::lang::reification::{downcast, Expr};
 use aries_model::state::{Domains, Event, Explanation};
-use aries_model::{Model, WriterId};
+use aries_model::WriterId;
 use itertools::Itertools;
 use smallvec::alloc::collections::VecDeque;
-use std::convert::TryFrom;
 
 #[derive(Clone)]
 struct ClauseLocks {
@@ -499,69 +498,51 @@ impl SatSolver {
 }
 
 impl BindSplit for SatSolver {
-    fn enforce_true(&mut self, e: &Expr, model: &mut Model) -> BindingResult {
-        match e.fun {
-            Fun::Or => {
-                let mut lits = Vec::with_capacity(e.args.len());
-                for &a in &e.args {
-                    let a = BAtom::try_from(a).expect("not a boolean");
-                    let lit = model.reify(a);
-                    lits.push(lit);
-                }
-                if let Some(clause) = Disjunction::new_non_tautological(lits) {
-                    self.add_clause(clause);
-                }
-                BindingResult::Refined
-            }
-            _ => BindingResult::Unsupported,
+    fn enforce_true(&mut self, expr: &Expr, _: &mut Domains) -> BindingResult {
+        if let Some(disjunction) = downcast::<Disjunction>(expr) {
+            self.add_clause(disjunction);
+            BindingResult::Enforced
+        } else {
+            BindingResult::Unsupported
         }
     }
 
-    fn enforce_false(&mut self, expr: &Expr, model: &mut Model) -> BindingResult {
-        match expr.fun {
-            Fun::Or => {
-                // (not (or a b ...))
-                //enforce the equivalent (and (not a) (not b) ....)
-                for &a in &expr.args {
-                    let a = BAtom::try_from(a).expect("not a boolean");
-                    let lit = model.reify(a);
-                    self.add_clause([!lit]);
-                }
-                BindingResult::Refined
+    fn enforce_false(&mut self, expr: &Expr, _: &mut Domains) -> BindingResult {
+        // (not (or a b ...))
+        //enforce the equivalent (and (not a) (not b) ....)
+        if let Some(disjunction) = downcast::<Disjunction>(expr) {
+            for l in disjunction {
+                self.add_clause([!l]);
             }
-            _ => BindingResult::Unsupported,
+            BindingResult::Enforced
+        } else {
+            BindingResult::Unsupported
         }
     }
 
-    fn enforce_eq(&mut self, reif: Lit, e: &Expr, model: &mut Model) -> BindingResult {
-        match e.fun {
-            Fun::Or => {
-                // l  <=>  (or a b ...)
-                // first, transform all Atoms into literals
-                let mut disjuncts = Vec::with_capacity(e.args.len());
-                for &a in &e.args {
-                    let a = BAtom::try_from(a).expect("not a boolean");
-                    let lit = model.reify(a);
-                    disjuncts.push(lit);
-                }
-                let mut clause = Vec::with_capacity(disjuncts.len() + 1);
-                // make l => (or a b ...)    <=>   (or (not l) a b ...)
-                clause.push(!reif);
-                disjuncts.iter().for_each(|l| clause.push(*l));
+    fn enforce_eq(&mut self, literal: Lit, expr: &Expr, _: &mut Domains) -> BindingResult {
+        if let Some(disjunction) = downcast::<Disjunction>(expr) {
+            // l  <=>  (or a b ...)
+            // first, transform all Atoms into literals
+
+            let mut clause = Vec::with_capacity(disjunction.len() + 1);
+            // make l => (or a b ...)    <=>   (or (not l) a b ...)
+            clause.push(!literal);
+            disjunction.into_iter().for_each(|l| clause.push(l));
+            if let Some(clause) = Disjunction::new_non_tautological(clause) {
+                self.add_clause(clause);
+            }
+            // make (or a b ...) => l    <=> (and (a => l) (b => l) ...)
+            for disjunct in disjunction {
+                // enforce a => l
+                let clause = vec![!disjunct, literal];
                 if let Some(clause) = Disjunction::new_non_tautological(clause) {
                     self.add_clause(clause);
                 }
-                // make (or a b ...) => l    <=> (and (a => l) (b => l) ...)
-                for disjunct in disjuncts {
-                    // enforce a => l
-                    let clause = vec![!disjunct, reif];
-                    if let Some(clause) = Disjunction::new_non_tautological(clause) {
-                        self.add_clause(clause);
-                    }
-                }
-                BindingResult::Refined
             }
-            _ => BindingResult::Unsupported,
+            BindingResult::Enforced
+        } else {
+            BindingResult::Unsupported
         }
     }
 }
@@ -589,7 +570,9 @@ mod tests {
     use aries_model::extensions::AssignmentExt;
     use aries_model::lang::IntCst;
     use aries_model::state::Cause;
-    use aries_model::{Model, WriterId};
+    use aries_model::WriterId;
+
+    type Model = aries_model::Model<&'static str>;
 
     #[test]
     fn test_propagation_simple() {
diff --git a/solver/src/solver/search.rs b/solver/src/solver/search.rs
index bd3f78f0..84fbf307 100644
--- a/solver/src/solver/search.rs
+++ b/solver/src/solver/search.rs
@@ -5,19 +5,19 @@ use aries_backtrack::Backtrack;
 use aries_model::bounds::Lit;
 use aries_model::extensions::SavedAssignment;
 use aries_model::lang::{IntCst, VarRef};
-use aries_model::Model;
+use aries_model::{Label, Model};
 
 pub enum Decision {
     SetLiteral(Lit),
     Restart,
 }
 
-pub fn default_brancher() -> Box<dyn SearchControl + Send> {
+pub fn default_brancher<Lbl: Label>() -> Box<dyn SearchControl<Lbl> + Send> {
     Box::new(activity::ActivityBrancher::new())
 }
 
 #[allow(unused_variables)]
-pub trait SearchControl: Backtrack {
+pub trait SearchControl<Lbl>: Backtrack {
     /// Select the next decision to make while maintaining the invariant that every non bound variable remains in the queue.
     ///
     /// This invariant allows to invoke this function at the decision level preceding the one of the decision that will be returned.
@@ -25,9 +25,9 @@ pub trait SearchControl: Backtrack {
     /// to the level preceding the decision to be made.
     ///
     /// Returns `None` if no decision is left to be made.
-    fn next_decision(&mut self, stats: &Stats, model: &Model) -> Option<Decision>;
+    fn next_decision(&mut self, stats: &Stats, model: &Model<Lbl>) -> Option<Decision>;
 
-    fn import_vars(&mut self, model: &Model) {}
+    fn import_vars(&mut self, model: &Model<Lbl>) {}
 
     fn set_default_value(&mut self, var: VarRef, val: IntCst) {}
 
@@ -36,9 +36,9 @@ pub trait SearchControl: Backtrack {
 
     /// Increase the activity of the variable and perform an reordering in the queue.
     /// The activity is then used to select the next variable.
-    fn bump_activity(&mut self, bvar: VarRef, model: &Model) {}
+    fn bump_activity(&mut self, bvar: VarRef, model: &Model<Lbl>) {}
 
     fn decay_activities(&mut self) {}
 
-    fn clone_to_box(&self) -> Box<dyn SearchControl + Send>;
+    fn clone_to_box(&self) -> Box<dyn SearchControl<Lbl> + Send>;
 }
diff --git a/solver/src/solver/search/activity.rs b/solver/src/solver/search/activity.rs
index 577035df..c735c3a1 100644
--- a/solver/src/solver/search/activity.rs
+++ b/solver/src/solver/search/activity.rs
@@ -4,12 +4,13 @@ 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_model::extensions::{AssignmentExt, SavedAssignment};
+use aries_model::extensions::{AssignmentExt, SavedAssignment, Shaped};
 use aries_model::lang::{IntCst, VarRef};
 use aries_model::state::{Event, IntDomain};
-use aries_model::Model;
+use aries_model::{Label, Model};
 use env_param::EnvParam;
 use itertools::Itertools;
+use std::sync::Arc;
 
 pub static PREFER_MIN_VALUE: EnvParam<bool> = EnvParam::new("ARIES_SMT_PREFER_MIN_VALUE", "true");
 pub static INITIALLY_ALLOWED_CONFLICTS: EnvParam<u64> = EnvParam::new("ARIES_SMT_INITIALLY_ALLOWED_CONFLICT", "100");
@@ -34,10 +35,28 @@ impl Default for BranchingParams {
     }
 }
 
+pub trait Heuristic<Lbl>: Send + Sync + 'static {
+    /// Specifies at which decision stage the given variable should be handled.
+    ///
+    /// The brancher only starts branching on the variables of stage N when
+    /// all variables of stage (N - 1) are set.
+    fn decision_stage(&self, var: VarRef, label: Option<&Lbl>, model: &Model<Lbl>) -> u8;
+}
+
+/// Default branching heuristic that puts all variables in the same decision stage.
+pub struct DefaultHeuristic;
+
+impl<L> Heuristic<L> for DefaultHeuristic {
+    fn decision_stage(&self, _: VarRef, _: Option<&L>, _: &Model<L>) -> u8 {
+        0
+    }
+}
+
 /// A branching scheme that first select variables that were recently involved in conflicts.
 #[derive(Clone)]
-pub struct ActivityBrancher {
+pub struct ActivityBrancher<Lbl> {
     pub params: BranchingParams,
+    heuristic: Arc<dyn Heuristic<Lbl>>,
     heap: VarSelect,
     default_assignment: DefaultValues,
     conflicts_at_last_restart: u64,
@@ -56,14 +75,23 @@ struct DefaultValues {
     values: RefMap<VarRef, IntCst>,
 }
 
-impl ActivityBrancher {
+impl<Lbl: Label> ActivityBrancher<Lbl> {
     pub fn new() -> Self {
-        Self::with_params(Default::default())
+        Self::new_with(Default::default(), DefaultHeuristic)
     }
 
-    pub fn with_params(params: BranchingParams) -> Self {
+    pub fn new_with_heuristic(h: impl Heuristic<Lbl>) -> Self {
+        Self::new_with(BranchingParams::default(), h)
+    }
+
+    pub fn new_with_params(params: BranchingParams) -> Self {
+        Self::new_with(params, DefaultHeuristic)
+    }
+
+    pub fn new_with(params: BranchingParams, h: impl Heuristic<Lbl>) -> Self {
         ActivityBrancher {
             params,
+            heuristic: Arc::new(h),
             heap: VarSelect::new(Default::default()),
             default_assignment: DefaultValues::default(),
             conflicts_at_last_restart: 0,
@@ -73,16 +101,12 @@ impl ActivityBrancher {
         }
     }
 
-    /// TODO: we should clarify and document this variable priority
-    fn priority(&self, variable: VarRef, model: &Model) -> u8 {
-        if model.var_domain(variable).size() <= 1 {
-            0
-        } else {
-            1
-        }
+    fn priority(&self, variable: VarRef, model: &Model<Lbl>) -> u8 {
+        self.heuristic
+            .decision_stage(variable, model.get_label(variable), model)
     }
 
-    pub fn import_vars(&mut self, model: &Model) {
+    pub fn import_vars(&mut self, model: &Model<Lbl>) {
         let mut count = 0;
         // go through the model's variables and declare any newly declared variable
         // TODO: use `advance_by` when it is stabilized. The current usage of `dropping` is very expensive
@@ -118,7 +142,7 @@ impl ActivityBrancher {
     /// to the level preceding the decision to be made.
     ///
     /// Returns `None` if no decision is left to be made.
-    pub fn next_decision(&mut self, stats: &Stats, model: &Model) -> Option<Decision> {
+    pub fn next_decision(&mut self, stats: &Stats, model: &Model<Lbl>) -> Option<Decision> {
         self.import_vars(model);
 
         let mut popper = self.heap.extractor();
@@ -196,7 +220,7 @@ impl ActivityBrancher {
     /// Increase the activity of the variable and perform an reordering in the queue.
     /// If the variable is optional, the activity of the presence variable is increased as well.
     /// The activity is then used to select the next variable.
-    pub fn bump_activity(&mut self, bvar: VarRef, model: &Model) {
+    pub fn bump_activity(&mut self, bvar: VarRef, model: &Model<Lbl>) {
         self.heap.var_bump_activity(bvar);
         match model.state.presence(bvar).variable() {
             VarRef::ZERO => {}
@@ -205,7 +229,7 @@ impl ActivityBrancher {
     }
 }
 
-impl Default for ActivityBrancher {
+impl<Lbl: Label> Default for ActivityBrancher<Lbl> {
     fn default() -> Self {
         Self::new()
     }
@@ -393,7 +417,7 @@ impl<'a> Popper<'a> {
     }
 }
 
-impl Backtrack for ActivityBrancher {
+impl<Lbl> Backtrack for ActivityBrancher<Lbl> {
     fn save_state(&mut self) -> DecLvl {
         self.heap.save_state()
     }
@@ -407,12 +431,12 @@ impl Backtrack for ActivityBrancher {
     }
 }
 
-impl SearchControl for ActivityBrancher {
-    fn next_decision(&mut self, stats: &Stats, model: &Model) -> Option<Decision> {
+impl<Lbl: Label> SearchControl<Lbl> for ActivityBrancher<Lbl> {
+    fn next_decision(&mut self, stats: &Stats, model: &Model<Lbl>) -> Option<Decision> {
         self.next_decision(stats, model)
     }
 
-    fn import_vars(&mut self, model: &Model) {
+    fn import_vars(&mut self, model: &Model<Lbl>) {
         self.import_vars(model)
     }
 
@@ -430,13 +454,13 @@ impl SearchControl for ActivityBrancher {
             .unwrap_or(true);
         if USE_LNS.get() && is_improvement {
             self.default_assignment.objective_found = Some(objective);
-            for (var, val) in assignment.state.bound_variables() {
+            for (var, val) in assignment.bound_variables() {
                 self.set_default_value(var, val);
             }
         }
     }
 
-    fn bump_activity(&mut self, bvar: VarRef, model: &Model) {
+    fn bump_activity(&mut self, bvar: VarRef, model: &Model<Lbl>) {
         self.bump_activity(bvar, model)
     }
 
@@ -444,7 +468,7 @@ impl SearchControl for ActivityBrancher {
         self.heap.decay_activities()
     }
 
-    fn clone_to_box(&self) -> Box<dyn SearchControl + Send> {
+    fn clone_to_box(&self) -> Box<dyn SearchControl<Lbl> + Send> {
         Box::new(self.clone())
     }
 }
diff --git a/solver/src/solver/theory_solver.rs b/solver/src/solver/theory_solver.rs
index 23c48c6b..5ba21b25 100644
--- a/solver/src/solver/theory_solver.rs
+++ b/solver/src/solver/theory_solver.rs
@@ -2,12 +2,10 @@ use crate::solver::BindingResult;
 use crate::{Contradiction, Theory};
 use aries_backtrack::{Backtrack, DecLvl};
 use aries_model::bounds::Lit;
-use aries_model::expressions::ExprHandle;
+use aries_model::lang::reification::Expr;
 use aries_model::state::Domains;
-use aries_model::Model;
 
 // TODO: remove this useless layer
-
 pub struct TheorySolver {
     pub theory: Box<dyn Theory>,
 }
@@ -17,12 +15,12 @@ impl TheorySolver {
         TheorySolver { theory }
     }
 
-    pub fn bind(&mut self, lit: Lit, expr: ExprHandle, model: &mut Model) -> BindingResult {
-        self.theory.bind(lit, expr, model)
+    pub fn bind(&mut self, lit: Lit, expr: &Expr, domains: &mut Domains) -> BindingResult {
+        self.theory.bind(lit, expr, domains)
     }
 
-    pub fn process(&mut self, model: &mut Domains) -> Result<(), Contradiction> {
-        self.theory.propagate(model)
+    pub fn process(&mut self, domains: &mut Domains) -> Result<(), Contradiction> {
+        self.theory.propagate(domains)
     }
 
     pub fn print_stats(&self) {
diff --git a/tnet/src/stn.rs b/tnet/src/stn.rs
index fd88bdc5..34840060 100644
--- a/tnet/src/stn.rs
+++ b/tnet/src/stn.rs
@@ -8,7 +8,7 @@ use aries_solver::{Contradiction, Theory};
 #[derive(Clone)]
 pub struct Stn {
     pub(crate) stn: StnTheory,
-    pub model: Model,
+    pub model: Model<String>,
 }
 impl Stn {
     pub fn new() -> Self {
@@ -36,11 +36,12 @@ impl Stn {
 
     pub fn add_edge(&mut self, source: Timepoint, target: Timepoint, weight: W) -> EdgeId {
         self.stn
-            .add_reified_edge(Lit::TRUE, source, target, weight, &self.model)
+            .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)
+        self.stn
+            .add_reified_edge(literal, source, target, weight, &self.model.state)
     }
 
     pub fn add_optional_true_edge(
@@ -59,7 +60,7 @@ impl Stn {
             forward_prop,
             backward_prop,
             presence,
-            &self.model,
+            &self.model.state,
         )
     }
 
diff --git a/tnet/src/theory.rs b/tnet/src/theory.rs
index 8b707117..2e30a443 100644
--- a/tnet/src/theory.rs
+++ b/tnet/src/theory.rs
@@ -5,12 +5,12 @@ 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::expressions::ExprHandle;
-use aries_model::extensions::{AssignmentExt, ExpressionFactoryExt, Shaped};
-use aries_model::lang::{Fun, IAtom, IntCst, VarRef};
+use aries_model::lang::normal_form::{NFLeq, NFOptLeq};
+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::{Model, WriterId};
+use aries_model::WriterId;
 use aries_solver::solver::BindingResult;
 use aries_solver::{Bind, Contradiction, Theory};
 use env_param::EnvParam;
@@ -683,15 +683,15 @@ impl StnTheory {
         source: impl Into<Timepoint>,
         target: impl Into<Timepoint>,
         weight: W,
-        model: &Model,
+        domains: &Domains,
     ) -> EdgeId {
         let e = self.add_inactive_constraint(source.into(), target.into(), weight).0;
 
-        if model.state.entails(literal) {
-            assert_eq!(model.state.entailing_level(literal), DecLvl::ROOT);
+        if domains.entails(literal) {
+            assert_eq!(domains.entailing_level(literal), DecLvl::ROOT);
             self.mark_active(e, literal);
-        } else if model.entails(!literal) {
-            assert_eq!(model.state.entailing_level(!literal), DecLvl::ROOT);
+        } 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);
@@ -721,21 +721,21 @@ impl StnTheory {
         forward_prop: Lit,
         backward_prop: Lit,
         presence: Option<Lit>,
-        model: &Model,
+        domains: &Domains,
     ) -> EdgeId {
         let e = self.add_inactive_constraint(source.into(), target.into(), weight).0;
 
         self.constraints
             .add_directed_enabler(e.forward(), forward_prop, presence);
-        if model.entails(forward_prop) {
-            assert_eq!(model.state.entailing_level(forward_prop), DecLvl::ROOT);
+        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 model.entails(backward_prop) {
-            assert_eq!(model.state.entailing_level(backward_prop), DecLvl::ROOT);
+        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));
         }
@@ -1565,84 +1565,37 @@ impl Backtrack for StnTheory {
 }
 
 impl Bind for StnTheory {
-    fn bind(&mut self, literal: Lit, expr: ExprHandle, model: &mut Model) -> BindingResult {
-        let expr = model.get_expr(expr);
-
-        // function that transforms the parameters into two `IAtom`s, panicking if it is not possible
-        let args_as_two_integers = || {
-            assert_eq!(expr.args.len(), 2);
-            let a = IAtom::try_from(expr.args[0]).expect("type error");
-            let b = IAtom::try_from(expr.args[1]).expect("type error");
-            (a, b)
-        };
-        // function that extracts the variable inside an IAtom, panicking if it is not possible
-        let var_in = |a: IAtom| match a.var {
-            Some(v) => v,
-            None => panic!("leq with no variable on the left side"),
-        };
-
-        match expr.fun {
-            Fun::Leq => {
-                let (a, b) = args_as_two_integers();
-                let va = var_in(a);
-                let vb = var_in(b);
-
-                // va + da <= vb + db    <=>   va - vb <= db - da
-                self.add_reified_edge(literal, vb, va, b.shift - a.shift, model);
-
+    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
-            }
-            Fun::Eq => {
-                let (a, b) = args_as_two_integers();
-                let x = model.leq(a, b);
-                let y = model.leq(b, a);
-                let and_x_y = model.and2(x, y);
-                model.bind(and_x_y, literal);
-                BindingResult::Refined
-            }
-            Fun::OptEq if model.entails(literal) => {
-                let (a, b) = args_as_two_integers();
-                let leq_a_b = model.opt_leq(a, b);
-                model.enforce(leq_a_b);
-                let leq_b_a = model.opt_leq(b, a);
-                model.enforce(leq_b_a);
-                BindingResult::Refined
-            }
-            Fun::OptLeq if model.entails(literal) => {
-                let (a, b) = args_as_two_integers();
-                let va = var_in(a);
-                let vb = var_in(b);
-
-                // va + da <= vb + db    <=>   va - vb <= db - da
-                let delay = b.shift - a.shift;
-                let a = va.into();
-                let b = vb.into();
-
-                let a_to_b = can_propagate(&model.state, a, b);
-                let b_to_a = can_propagate(&model.state, b, a);
-                let presence = edge_presence(&model.state, a, b);
-                self.add_optional_true_edge(b, a, delay, b_to_a, a_to_b, presence, model);
-                BindingResult::Enforced
-            }
-            Fun::OptLeq if model.entails(!literal) => {
+            } else if doms.entails(!literal) {
                 // this constraint is always false, post the opposite
-                let (a, b) = args_as_two_integers();
-                let va = var_in(a);
-                let vb = var_in(b);
-
-                // va + da <= vb + db    <=>   va - vb <= db - da
-                let delay = a.shift - b.shift - 1;
-                let a = vb.into();
-                let b = va.into();
-
-                let a_to_b = can_propagate(&model.state, a, b);
-                let b_to_a = can_propagate(&model.state, b, a);
-                let presence = edge_presence(&model.state, a, b);
-                self.add_optional_true_edge(b, a, delay, b_to_a, a_to_b, presence, model);
+                // !(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
             }
-
-            _ => BindingResult::Unsupported,
+        } else {
+            BindingResult::Unsupported
         }
     }
 }
@@ -1683,6 +1636,7 @@ pub(crate) fn edge_presence(doms: &Domains, var1: Timepoint, var2: Timepoint) ->
 mod tests {
     use super::*;
     use crate::stn::Stn;
+    use aries_model::extensions::AssignmentExt;
     use aries_model::lang::IVar;
 
     #[test]
diff --git a/utils/src/input.rs b/utils/src/input.rs
index 525b924b..82a872d1 100644
--- a/utils/src/input.rs
+++ b/utils/src/input.rs
@@ -370,6 +370,11 @@ impl From<&Sym> for Sym {
         s.clone()
     }
 }
+impl From<&Sym> for String {
+    fn from(s: &Sym) -> Self {
+        s.symbol.clone()
+    }
+}
 
 impl std::fmt::Display for Sym {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {


Test patch:
diff --git a/examples/smt/tests/smt.rs b/examples/smt/tests/smt.rs
index 7c7d04d7..62c4c99e 100644
--- a/examples/smt/tests/smt.rs
+++ b/examples/smt/tests/smt.rs
@@ -1,12 +1,14 @@
 use aries_backtrack::Backtrack;
 use aries_model::bounds::Lit;
-use aries_model::extensions::{AssignmentExt, ExpressionFactoryExt};
-use aries_model::lang::{BAtom, IVar};
+use aries_model::extensions::AssignmentExt;
+use aries_model::lang::expr::*;
+use aries_model::lang::IVar;
 use aries_model::state::OptDomain;
-use aries_model::Model;
 use aries_solver::solver::Solver;
 use aries_tnet::theory::{StnConfig, StnTheory};
 
+type Model = aries_model::Model<String>;
+
 #[test]
 fn sat() {
     let mut model = Model::new();
@@ -17,9 +19,8 @@ fn sat() {
     solver.enforce(a);
     assert!(solver.solve().unwrap().is_some());
     assert_eq!(solver.model.boolean_value_of(a), Some(true));
-    let c = solver.model.implies(a, b);
     solver.reset();
-    solver.enforce(c);
+    solver.enforce(implies(a, b));
     assert!(solver.solve().unwrap().is_some());
     assert_eq!(solver.model.boolean_value_of(a), Some(true));
     assert_eq!(solver.model.boolean_value_of(b), Some(true));
@@ -36,12 +37,12 @@ fn diff_logic() {
     let b = model.new_ivar(0, 10, "b");
     let c = model.new_ivar(0, 10, "c");
 
-    let constraints = vec![model.lt(a, b), model.lt(b, c), model.lt(c, a)];
+    let constraints = vec![lt(a, b), lt(b, c), lt(c, a)];
 
     let mut solver = Solver::new(model);
 
     solver.add_theory(|tok| StnTheory::new(tok, StnConfig::default()));
-    solver.enforce_all(&constraints);
+    solver.enforce_all(constraints);
     assert!(solver.solve().unwrap().is_none());
 }
 
@@ -52,20 +53,21 @@ fn minimize() {
     let b = model.new_ivar(0, 10, "b");
     let c = model.new_ivar(0, 10, "c");
 
-    let x = model.geq(b, 6);
-    let y = model.geq(b, 8);
+    model.enforce(lt(a, b));
+    model.enforce(lt(b, c));
+    model.enforce(lt(a, c));
+    let x = model.reify(geq(b, 6));
+    let y = model.reify(geq(b, 8));
+    model.enforce(or([x, y]));
 
-    let constraints = vec![model.lt(a, b), model.lt(b, c), model.lt(a, c), model.or2(x, y)];
     let mut solver = Solver::new(model);
-
     solver.add_theory(|tok| StnTheory::new(tok, StnConfig::default()));
-    solver.enforce_all(&constraints);
+
     assert!(solver.solve().unwrap().is_some());
     match solver.minimize(c).unwrap() {
         None => panic!(),
         Some((val, _)) => assert_eq!(val, 7),
     }
-    solver.print_stats()
 }
 
 #[test]
@@ -73,15 +75,14 @@ fn minimize_small() {
     let mut model = Model::new();
     let a = model.new_ivar(0, 10, "a");
 
-    let x = model.geq(a, 6);
-    let y = model.geq(a, 8);
+    let x = model.reify(geq(a, 6));
+    let y = model.reify(geq(a, 8));
 
-    let constraints = vec![model.or2(x, y)];
+    model.enforce(or([x, y]));
 
     let mut solver = Solver::new(model);
 
     solver.add_theory(|tok| StnTheory::new(tok, StnConfig::default()));
-    solver.enforce_all(&constraints);
     assert!(solver.solve().unwrap().is_some());
     match solver.minimize(a).unwrap() {
         None => panic!(),
@@ -103,26 +104,26 @@ fn int_bounds() {
     let h = model.new_ivar(0, 10, "h");
 
     let constraints = vec![
-        model.leq(a, 8),
-        model.leq(2, a),
-        model.lt(1, b),
-        model.lt(b, 9),
-        model.geq(c, 2),
-        model.geq(8, c),
-        model.gt(d, 1),
-        model.gt(9, d),
-        !model.gt(e, 8),
-        !model.gt(2, e),
-        !model.geq(1, f),
-        !model.geq(f, 9),
-        !model.lt(g, 2),
-        !model.lt(8, g),
-        !model.leq(h, 1),
-        !model.leq(9, h),
+        leq(a, 8),
+        leq(2, a),
+        lt(1, b),
+        lt(b, 9),
+        geq(c, 2),
+        geq(8, c),
+        gt(d, 1),
+        gt(9, d),
+        !gt(e, 8),
+        !gt(2, e),
+        !geq(1, f),
+        !geq(f, 9),
+        !lt(g, 2),
+        !lt(8, g),
+        !leq(h, 1),
+        !leq(9, h),
     ];
 
     let mut solver = Solver::new(model);
-    solver.enforce_all(&constraints);
+    solver.enforce_all(constraints);
     assert!(solver.propagate_and_backtrack_to_consistent());
     let check_dom = |v, lb, ub| {
         assert_eq!(solver.model.domain_of(v), (lb, ub));
@@ -162,8 +163,10 @@ fn bools_as_ints() {
     assert_eq!(solver.model.boolean_value_of(a), None);
     assert_eq!(solver.model.domain_of(ia), (0, 1));
 
-    let constraints: Vec<BAtom> = vec![a.into(), (!b).into(), solver.model.geq(ic, 1), solver.model.leq(id, 0)];
-    solver.enforce_all(&constraints);
+    solver.enforce(a.true_lit());
+    solver.enforce(b.false_lit());
+    solver.enforce(geq(ic, 1));
+    solver.enforce(leq(id, 0));
 
     solver.propagate().unwrap();
     assert_eq!(solver.model.boolean_value_of(a), Some(true));
@@ -191,22 +194,19 @@ fn ints_and_bools() {
     assert_eq!(solver.model.domain_of(ia), (0, 1));
     assert_eq!(solver.model.boolean_value_of(a), None);
 
-    let constraint = solver.model.leq(i, ia);
-    solver.enforce(constraint);
+    solver.enforce(leq(i, ia));
     assert!(solver.propagate_and_backtrack_to_consistent());
     assert_eq!(solver.model.domain_of(i), (-10, 1));
     assert_eq!(solver.model.domain_of(ia), (0, 1));
     assert_eq!(solver.model.boolean_value_of(a), None);
 
-    let constraint = solver.model.gt(ia, i);
-    solver.enforce(constraint);
+    solver.enforce(gt(ia, i));
     assert!(solver.propagate_and_backtrack_to_consistent());
     assert_eq!(solver.model.domain_of(i), (-10, 0));
     assert_eq!(solver.model.domain_of(ia), (0, 1));
     assert_eq!(solver.model.boolean_value_of(a), None);
 
-    let constraint = solver.model.geq(i, 0);
-    solver.enforce(constraint);
+    solver.enforce(geq(i, 0));
     assert!(solver.propagate_and_backtrack_to_consistent());
     assert_eq!(solver.model.domain_of(i), (0, 0));
     assert_eq!(solver.model.domain_of(ia), (1, 1));
@@ -231,24 +231,21 @@ fn optional_hierarchy() {
         .map(|(i, (lb, ub))| model.new_optional_ivar(*lb, *ub, scopes[i], format!("i_{}", i)))
         .collect();
 
-    let mut constraints = Vec::with_capacity(32);
-
     // at least one must be present
     // constraints.push(model.or(&scopes.iter().map(|&lit| BAtom::from(lit)).collect::<Vec<_>>()));
 
     for &sub_var in &vars {
-        constraints.push(model.opt_eq(i, sub_var));
+        model.enforce(opt_eq(i, sub_var));
     }
 
     let mut solver = Solver::new(model);
     solver.add_theory(|tok| StnTheory::new(tok, StnConfig::default()));
 
-    // solver.model.state.print();
+    // solver.model.print_state();
 
-    solver.enforce_all(&constraints);
     assert!(solver.propagate_and_backtrack_to_consistent());
 
-    // solver.model.state.print();
+    // solver.model.print_state();
 
     assert_eq!(solver.model.opt_domain_of(i), Unknown(-10, 10));
     assert_eq!(solver.model.opt_domain_of(vars[0]), Unknown(0, 8));