From a9b926c69b3e5441eae7abf3f611aa8916cbeac5 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Wed, 15 Mar 2023 11:36:31 -0400
Subject: [PATCH 01/14] add expr

---
 motile/expressions.py | 163 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 163 insertions(+)
 create mode 100644 motile/expressions.py

diff --git a/motile/expressions.py b/motile/expressions.py
new file mode 100644
index 0000000..a42ad48
--- /dev/null
+++ b/motile/expressions.py
@@ -0,0 +1,163 @@
+from __future__ import annotations
+
+import ast
+from numbers import Number
+from typing import Any, Sequence
+
+
+class Expr(ast.AST):
+    # comparisons
+    @staticmethod
+    def _cast(obj: Any) -> Expr:
+        """Cast object into an Expression."""
+        return obj if isinstance(obj, Expr) else Constant(obj)
+
+    def __lt__(self, other: Expr | float) -> Compare:
+        return Compare(self, [ast.Lt()], [other])
+
+    def __le__(self, other: Expr | float) -> Compare:
+        return Compare(self, [ast.LtE()], [other])
+
+    def __eq__(self, other: Expr | float) -> Compare:  # type: ignore
+        return Compare(self, [ast.Eq()], [other])
+
+    def __ne__(self, other: Expr | float) -> Compare:  # type: ignore
+        return Compare(self, [ast.NotEq()], [other])
+
+    def __gt__(self, other: Expr | float) -> Compare:
+        return Compare(self, [ast.Gt()], [other])
+
+    def __ge__(self, other: Expr | float) -> Compare:
+        return Compare(self, [ast.GtE()], [other])
+
+    # binary operators
+    # (note that __and__ and __or__ are reserved for boolean operators.)
+
+    def __add__(self, other: Expr) -> BinOp:
+        return BinOp(self, ast.Add(), other)
+
+    def __radd__(self, other: Expr) -> BinOp:
+        return BinOp(self, ast.Add(), other)
+
+    def __sub__(self, other: Expr) -> BinOp:
+        return BinOp(self, ast.Sub(), other)
+
+    def __rmul__(self, other: float) -> BinOp:
+        return BinOp(other, ast.Mult(), self)
+
+    def __mul__(self, other: float) -> BinOp:
+        return BinOp(self, ast.Mult(), other)
+
+    def __truediv__(self, other: float) -> BinOp:
+        return BinOp(self, ast.Div(), other)
+
+    # unary operators
+
+    def __neg__(self) -> UnaryOp:
+        return UnaryOp(ast.USub(), self)
+
+    def __pos__(self) -> UnaryOp:
+        # usually a no-op
+        return UnaryOp(ast.UAdd(), self)
+
+
+class Compare(Expr, ast.Compare):
+    """A comparison of two or more values.
+
+    `left` is the first value in the comparison, `ops` the list of operators,
+    and `comparators` the list of values after the first element in the
+    comparison.
+    """
+
+    def __init__(
+        self,
+        left: Expr,
+        ops: Sequence[ast.cmpop],
+        comparators: Sequence[Expr | float],
+        **kwargs: Any,
+    ) -> None:
+        super().__init__(
+            Expr._cast(left),
+            ops,
+            [Expr._cast(c) for c in comparators],
+            **kwargs,
+        )
+
+
+class BinOp(Expr, ast.BinOp):
+    """A binary operation (like addition or division).
+
+    `op` is the operator, and `left` and `right` are any expression nodes.
+    """
+
+    def __init__(
+        self,
+        left: T | Expr,
+        op: ast.operator,
+        right: T | Expr,
+        **k: Any,
+    ) -> None:
+        super().__init__(Expr._cast(left), op, Expr._cast(right), **k)
+
+
+class UnaryOp(Expr, ast.UnaryOp):
+    """A unary operation.
+
+    `op` is the operator, and `operand` any expression node.
+    """
+
+    def __init__(self, op: ast.unaryop, operand: Expr, **kwargs: Any) -> None:
+        super().__init__(op, Expr._cast(operand), **kwargs)
+
+
+class Constant(Expr, ast.Constant):
+    """A constant value.
+
+    The `value` attribute contains the Python object it represents.
+    types supported: NoneType, str, bytes, bool, int, float
+    """
+
+    def __init__(self, value: Number, kind: str | None = None, **kwargs: Any) -> None:
+        if not isinstance(value, Number):
+            raise TypeError("Constants must be numbers")
+        super().__init__(value, kind, **kwargs)
+
+
+class Index(Expr, ast.Name):
+    """A solution index.
+
+    `id` holds the index as a string (becuase ast.Name requires a string).
+    """
+
+    def __init__(self, index: int) -> None:
+        self.index = index
+        super().__init__(str(index), ctx=ast.Load())
+
+
+from ilpy import LinearConstraint, Relation
+
+op_map: dict[ast.cmpop, Relation] = {
+    ast.LtE(): Relation.LessEqual,
+    ast.Eq(): Relation.Equal,
+    ast.Gt(): Relation.GreaterEqual,
+}
+
+
+def to_constraint(expr: Expr) -> LinearConstraint:
+    constraint = LinearConstraint()
+
+    if isinstance(expr, Compare):
+        if len(expr.ops) != 1:
+            raise ValueError("Only single comparisons are supported")
+
+        try:
+            constraint.set_relation(op_map[expr.ops[0]])
+        except KeyError as e:
+            raise ValueError(f"Unsupported comparison operator: {expr.ops[0]}") from e
+
+        right = expr.comparators[0]
+        if not isinstance(right, Constant):
+            raise ValueError("Right side of comparison must be a constant")
+        right_value = right.value
+
+        for 

From 5008fb32b468f5737d2d01d448e0e7e9640b7d9f Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Wed, 15 Mar 2023 14:30:10 -0400
Subject: [PATCH 02/14] add expression

---
 motile/constraints/select_edge_nodes.py | 26 +++----
 motile/expressions.py                   | 91 ++++++++++++++++++++-----
 motile/variables/variable.py            |  5 ++
 3 files changed, 94 insertions(+), 28 deletions(-)

diff --git a/motile/constraints/select_edge_nodes.py b/motile/constraints/select_edge_nodes.py
index a320cb6..4009e14 100644
--- a/motile/constraints/select_edge_nodes.py
+++ b/motile/constraints/select_edge_nodes.py
@@ -31,17 +31,19 @@ def instantiate(self, solver: Solver) -> list[ilpy.LinearConstraint]:
         constraints = []
         for edge in solver.graph.edges:
             u, v = edge
-
-            ind_e = edge_indicators[edge]
-            ind_u = node_indicators[u]
-            ind_v = node_indicators[v]
-
-            constraint = ilpy.LinearConstraint()
-            constraint.set_coefficient(ind_e, 2)
-            constraint.set_coefficient(ind_u, -1)
-            constraint.set_coefficient(ind_v, -1)
-            constraint.set_relation(ilpy.Relation.LessEqual)
-            constraint.set_value(0)
-            constraints.append(constraint)
+            ind_e = edge_indicators.expr(edge)
+            ind_u = node_indicators.expr(u)
+            ind_v = node_indicators.expr(v)
+
+            expr = 2 * ind_e - ind_u + ind_u - ind_v <= 0
+            constraints.append(expr.constraint())
+
+            # constraint = ilpy.LinearConstraint()
+            # constraint.set_coefficient(ind_e, 2)
+            # constraint.set_coefficient(ind_u, -1)
+            # constraint.set_coefficient(ind_v, -1)
+            # constraint.set_relation(ilpy.Relation.LessEqual)
+            # constraint.set_value(0)
+            # constraints.append(constraint)
 
         return constraints
diff --git a/motile/expressions.py b/motile/expressions.py
index a42ad48..15a09a3 100644
--- a/motile/expressions.py
+++ b/motile/expressions.py
@@ -4,8 +4,14 @@
 from numbers import Number
 from typing import Any, Sequence
 
+from ilpy import LinearConstraint, Relation
+
 
 class Expr(ast.AST):
+    def constraint(self) -> LinearConstraint:
+        """Return a linear constraint from this expression."""
+        return
+
     # comparisons
     @staticmethod
     def _cast(obj: Any) -> Expr:
@@ -134,30 +140,83 @@ def __init__(self, index: int) -> None:
         super().__init__(str(index), ctx=ast.Load())
 
 
-from ilpy import LinearConstraint, Relation
-
-op_map: dict[ast.cmpop, Relation] = {
-    ast.LtE(): Relation.LessEqual,
-    ast.Eq(): Relation.Equal,
-    ast.Gt(): Relation.GreaterEqual,
+op_map: dict[type[ast.cmpop], Relation] = {
+    ast.LtE: Relation.LessEqual,
+    ast.Eq: Relation.Equal,
+    ast.Gt: Relation.GreaterEqual,
 }
 
 
-def to_constraint(expr: Expr) -> LinearConstraint:
+def to_constraint(expr: ast.expr) -> LinearConstraint:
     constraint = LinearConstraint()
 
-    if isinstance(expr, Compare):
-        if len(expr.ops) != 1:
+    seen_compare = False
+    for sub in ast.walk(expr):
+        if not isinstance(sub, Compare):
+            continue
+        if seen_compare:
             raise ValueError("Only single comparisons are supported")
 
+        op_type = type(sub.ops[0])
         try:
-            constraint.set_relation(op_map[expr.ops[0]])
+            constraint.set_relation(op_map[op_type])
         except KeyError as e:
-            raise ValueError(f"Unsupported comparison operator: {expr.ops[0]}") from e
+            raise ValueError(f"Unsupported comparison operator: {op_type}") from e
 
-        right = expr.comparators[0]
-        if not isinstance(right, Constant):
-            raise ValueError("Right side of comparison must be a constant")
-        right_value = right.value
+        seen_compare = True
 
-        for 
+    for index, coeff in get_coefficients(expr).items():
+        if index is None:
+            constraint.set_value(-coeff)
+        else:
+            constraint.set_coefficient(index, coeff)
+    return constraint
+
+
+def get_coefficients(
+    expr: ast.expr, coeffs: dict[int | None, float] | None = None, scale: int = 1
+) -> dict[int | None, float]:
+    if coeffs is None:
+        coeffs = {}
+
+    if isinstance(expr, Compare):
+        if len(expr.ops) != 1:
+            raise ValueError("Only single comparisons are supported")
+        get_coefficients(expr.left, coeffs)
+        get_coefficients(expr.comparators[0], coeffs, -1)
+
+    elif isinstance(expr, BinOp):
+        if isinstance(expr.op, (ast.Mult, ast.Div)):
+            if isinstance(expr.right, Constant):
+                e = expr.left
+                v = expr.right.value
+            elif isinstance(expr.left, Constant):
+                e = expr.right
+                v = expr.left.value
+            else:
+                raise ValueError("Multiplication must be by a constant")
+            assert isinstance(e, Index)
+            scale *= 1 / v if isinstance(expr.op, ast.Div) else v
+            get_coefficients(e, coeffs, scale)
+
+        else:
+            get_coefficients(expr.left, coeffs, scale)
+            if isinstance(expr.op, (ast.USub, ast.Sub)):
+                scale = -scale
+            get_coefficients(expr.right, coeffs, scale)
+    elif isinstance(expr, UnaryOp):
+        if isinstance(expr.op, ast.USub):
+            scale = -scale
+        get_coefficients(expr.operand, coeffs, scale)
+    elif isinstance(expr, Constant):
+        coeffs[None] = expr.value * scale
+    elif isinstance(expr, Index):
+        coeffs.setdefault(expr.index, 0)
+        coeffs[expr.index] += scale
+    else:
+        raise ValueError("Unsupported expression")
+
+    return coeffs
+
+
+# -u + 2*e - (v + 4 - u)
diff --git a/motile/variables/variable.py b/motile/variables/variable.py
index 7ba74da..bd7fb0c 100644
--- a/motile/variables/variable.py
+++ b/motile/variables/variable.py
@@ -5,6 +5,8 @@
 
 import ilpy
 
+from motile.expressions import Index
+
 if TYPE_CHECKING:
     from motile.solver import Solver
 
@@ -122,3 +124,6 @@ def keys(self) -> Iterable[Hashable]:
 
     def values(self) -> Iterable[int]:
         return self._index_map.values()
+
+    def expr(self, key: Hashable) -> Index:
+        return Index(self._index_map[key])

From 0b717abd9c5c31b99d499b629f4f0aefc8d2eb23 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Wed, 15 Mar 2023 15:21:36 -0400
Subject: [PATCH 03/14] wip

---
 motile/constraints/max_children.py      | 19 ++++++++++-------
 motile/constraints/select_edge_nodes.py | 28 ++++++++++++++-----------
 motile/expressions.py                   |  9 ++++----
 3 files changed, 32 insertions(+), 24 deletions(-)

diff --git a/motile/constraints/max_children.py b/motile/constraints/max_children.py
index df1a3bb..0202d95 100644
--- a/motile/constraints/max_children.py
+++ b/motile/constraints/max_children.py
@@ -6,6 +6,7 @@
 
 from ..variables import EdgeSelected
 from .constraint import Constraint
+from motile.expressions import Constant, Expr
 
 if TYPE_CHECKING:
     from motile.solver import Solver
@@ -35,16 +36,20 @@ def instantiate(self, solver: Solver) -> list[ilpy.LinearConstraint]:
 
         constraints = []
         for node in solver.graph.nodes:
-            constraint = ilpy.LinearConstraint()
+            
+            n_edges = sum(edge_indicators.expr(e) for e in solver.graph.next_edges[node], Constant(0))
+            expr = (n_edges <= self.max_children)
+            constraints.append(expr.constraint())
 
-            # all outgoing edges
-            for edge in solver.graph.next_edges[node]:
-                constraint.set_coefficient(edge_indicators[edge], 1)
+            # constraint = ilpy.LinearConstraint()
+            # # all outgoing edges
+            # for edge in solver.graph.next_edges[node]:
+            #     constraint.set_coefficient(edge_indicators[edge], 1)
 
             # relation, value
-            constraint.set_relation(ilpy.Relation.LessEqual)
+            # constraint.set_relation(ilpy.Relation.LessEqual)
 
-            constraint.set_value(self.max_children)
-            constraints.append(constraint)
+            # constraint.set_value(self.max_children)
+            # constraints.append(constraint)
 
         return constraints
diff --git a/motile/constraints/select_edge_nodes.py b/motile/constraints/select_edge_nodes.py
index 4009e14..88f68fa 100644
--- a/motile/constraints/select_edge_nodes.py
+++ b/motile/constraints/select_edge_nodes.py
@@ -29,21 +29,25 @@ def instantiate(self, solver: Solver) -> list[ilpy.LinearConstraint]:
         edge_indicators = solver.get_variables(EdgeSelected)
 
         constraints = []
-        for edge in solver.graph.edges:
-            u, v = edge
-            ind_e = edge_indicators.expr(edge)
+        for (u, v) in solver.graph.edges:
+            ind_e = edge_indicators.expr((u, v))
+            ind_u = node_indicators.expr(u)
             ind_u = node_indicators.expr(u)
             ind_v = node_indicators.expr(v)
 
-            expr = 2 * ind_e - ind_u + ind_u - ind_v <= 0
-            constraints.append(expr.constraint())
+            # 2 x_e - x_u - x_v \leq 0
+            constraint = ilpy.LinearConstraint()
+            constraint.set_coefficient(ind_e, 2)
+            constraint.set_coefficient(ind_u, -1)
+            constraint.set_coefficient(ind_v, -1)
+            constraint.set_relation(ilpy.Relation.LessEqual)
+            constraint.set_value(0)
+            constraints.append(constraint)
+
+
+            expression = 2 * ind_e - ind_u - ind_v <= 0
+            constraints.append(expression.constraint())
 
-            # constraint = ilpy.LinearConstraint()
-            # constraint.set_coefficient(ind_e, 2)
-            # constraint.set_coefficient(ind_u, -1)
-            # constraint.set_coefficient(ind_v, -1)
-            # constraint.set_relation(ilpy.Relation.LessEqual)
-            # constraint.set_value(0)
-            # constraints.append(constraint)
 
         return constraints
+
diff --git a/motile/expressions.py b/motile/expressions.py
index 15a09a3..f659336 100644
--- a/motile/expressions.py
+++ b/motile/expressions.py
@@ -1,7 +1,6 @@
 from __future__ import annotations
 
 import ast
-from numbers import Number
 from typing import Any, Sequence
 
 from ilpy import LinearConstraint, Relation
@@ -10,7 +9,7 @@
 class Expr(ast.AST):
     def constraint(self) -> LinearConstraint:
         """Return a linear constraint from this expression."""
-        return
+        return to_constraint(self)
 
     # comparisons
     @staticmethod
@@ -123,8 +122,8 @@ class Constant(Expr, ast.Constant):
     types supported: NoneType, str, bytes, bool, int, float
     """
 
-    def __init__(self, value: Number, kind: str | None = None, **kwargs: Any) -> None:
-        if not isinstance(value, Number):
+    def __init__(self, value: float, kind: str | None = None, **kwargs: Any) -> None:
+        if not isinstance(value, (float, int)):
             raise TypeError("Constants must be numbers")
         super().__init__(value, kind, **kwargs)
 
@@ -147,7 +146,7 @@ def __init__(self, index: int) -> None:
 }
 
 
-def to_constraint(expr: ast.expr) -> LinearConstraint:
+def to_constraint(expr: Expr) -> LinearConstraint:
     constraint = LinearConstraint()
 
     seen_compare = False

From 94ca52f10345e85aaf48bf49190351de5b9a86cd Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Thu, 16 Mar 2023 20:31:57 -0400
Subject: [PATCH 04/14] wip

---
 motile/constraints/max_children.py      |  10 +-
 motile/constraints/select_edge_nodes.py |  25 +--
 motile/expressions.py                   | 221 ------------------------
 motile/solver.py                        |  15 +-
 motile/variables/node_appear.py         |  51 +-----
 motile/variables/variable.py            |   9 +-
 6 files changed, 40 insertions(+), 291 deletions(-)
 delete mode 100644 motile/expressions.py

diff --git a/motile/constraints/max_children.py b/motile/constraints/max_children.py
index 0202d95..f780d5a 100644
--- a/motile/constraints/max_children.py
+++ b/motile/constraints/max_children.py
@@ -3,10 +3,10 @@
 from typing import TYPE_CHECKING
 
 import ilpy
+from ilpy.expressions import Constant
 
 from ..variables import EdgeSelected
 from .constraint import Constraint
-from motile.expressions import Constant, Expr
 
 if TYPE_CHECKING:
     from motile.solver import Solver
@@ -36,9 +36,11 @@ def instantiate(self, solver: Solver) -> list[ilpy.LinearConstraint]:
 
         constraints = []
         for node in solver.graph.nodes:
-            
-            n_edges = sum(edge_indicators.expr(e) for e in solver.graph.next_edges[node], Constant(0))
-            expr = (n_edges <= self.max_children)
+            n_edges = sum(
+                [edge_indicators.expr(e) for e in solver.graph.next_edges[node]],
+                Constant(0),
+            )
+            expr = n_edges <= self.max_children
             constraints.append(expr.constraint())
 
             # constraint = ilpy.LinearConstraint()
diff --git a/motile/constraints/select_edge_nodes.py b/motile/constraints/select_edge_nodes.py
index 88f68fa..696a1bd 100644
--- a/motile/constraints/select_edge_nodes.py
+++ b/motile/constraints/select_edge_nodes.py
@@ -29,25 +29,12 @@ def instantiate(self, solver: Solver) -> list[ilpy.LinearConstraint]:
         edge_indicators = solver.get_variables(EdgeSelected)
 
         constraints = []
-        for (u, v) in solver.graph.edges:
-            ind_e = edge_indicators.expr((u, v))
-            ind_u = node_indicators.expr(u)
-            ind_u = node_indicators.expr(u)
-            ind_v = node_indicators.expr(v)
-
-            # 2 x_e - x_u - x_v \leq 0
-            constraint = ilpy.LinearConstraint()
-            constraint.set_coefficient(ind_e, 2)
-            constraint.set_coefficient(ind_u, -1)
-            constraint.set_coefficient(ind_v, -1)
-            constraint.set_relation(ilpy.Relation.LessEqual)
-            constraint.set_value(0)
-            constraints.append(constraint)
-
-
-            expression = 2 * ind_e - ind_u - ind_v <= 0
-            constraints.append(expression.constraint())
+        for u, v in solver.graph.edges:
+            x_e = edge_indicators.expr((u, v), "e")
+            x_u = node_indicators.expr(u, "u")
+            x_v = node_indicators.expr(v, "v")
 
+            expression = 2 * x_e - x_u - x_v <= 0
+            constraints.append(expression.constraint())
 
         return constraints
-
diff --git a/motile/expressions.py b/motile/expressions.py
deleted file mode 100644
index f659336..0000000
--- a/motile/expressions.py
+++ /dev/null
@@ -1,221 +0,0 @@
-from __future__ import annotations
-
-import ast
-from typing import Any, Sequence
-
-from ilpy import LinearConstraint, Relation
-
-
-class Expr(ast.AST):
-    def constraint(self) -> LinearConstraint:
-        """Return a linear constraint from this expression."""
-        return to_constraint(self)
-
-    # comparisons
-    @staticmethod
-    def _cast(obj: Any) -> Expr:
-        """Cast object into an Expression."""
-        return obj if isinstance(obj, Expr) else Constant(obj)
-
-    def __lt__(self, other: Expr | float) -> Compare:
-        return Compare(self, [ast.Lt()], [other])
-
-    def __le__(self, other: Expr | float) -> Compare:
-        return Compare(self, [ast.LtE()], [other])
-
-    def __eq__(self, other: Expr | float) -> Compare:  # type: ignore
-        return Compare(self, [ast.Eq()], [other])
-
-    def __ne__(self, other: Expr | float) -> Compare:  # type: ignore
-        return Compare(self, [ast.NotEq()], [other])
-
-    def __gt__(self, other: Expr | float) -> Compare:
-        return Compare(self, [ast.Gt()], [other])
-
-    def __ge__(self, other: Expr | float) -> Compare:
-        return Compare(self, [ast.GtE()], [other])
-
-    # binary operators
-    # (note that __and__ and __or__ are reserved for boolean operators.)
-
-    def __add__(self, other: Expr) -> BinOp:
-        return BinOp(self, ast.Add(), other)
-
-    def __radd__(self, other: Expr) -> BinOp:
-        return BinOp(self, ast.Add(), other)
-
-    def __sub__(self, other: Expr) -> BinOp:
-        return BinOp(self, ast.Sub(), other)
-
-    def __rmul__(self, other: float) -> BinOp:
-        return BinOp(other, ast.Mult(), self)
-
-    def __mul__(self, other: float) -> BinOp:
-        return BinOp(self, ast.Mult(), other)
-
-    def __truediv__(self, other: float) -> BinOp:
-        return BinOp(self, ast.Div(), other)
-
-    # unary operators
-
-    def __neg__(self) -> UnaryOp:
-        return UnaryOp(ast.USub(), self)
-
-    def __pos__(self) -> UnaryOp:
-        # usually a no-op
-        return UnaryOp(ast.UAdd(), self)
-
-
-class Compare(Expr, ast.Compare):
-    """A comparison of two or more values.
-
-    `left` is the first value in the comparison, `ops` the list of operators,
-    and `comparators` the list of values after the first element in the
-    comparison.
-    """
-
-    def __init__(
-        self,
-        left: Expr,
-        ops: Sequence[ast.cmpop],
-        comparators: Sequence[Expr | float],
-        **kwargs: Any,
-    ) -> None:
-        super().__init__(
-            Expr._cast(left),
-            ops,
-            [Expr._cast(c) for c in comparators],
-            **kwargs,
-        )
-
-
-class BinOp(Expr, ast.BinOp):
-    """A binary operation (like addition or division).
-
-    `op` is the operator, and `left` and `right` are any expression nodes.
-    """
-
-    def __init__(
-        self,
-        left: T | Expr,
-        op: ast.operator,
-        right: T | Expr,
-        **k: Any,
-    ) -> None:
-        super().__init__(Expr._cast(left), op, Expr._cast(right), **k)
-
-
-class UnaryOp(Expr, ast.UnaryOp):
-    """A unary operation.
-
-    `op` is the operator, and `operand` any expression node.
-    """
-
-    def __init__(self, op: ast.unaryop, operand: Expr, **kwargs: Any) -> None:
-        super().__init__(op, Expr._cast(operand), **kwargs)
-
-
-class Constant(Expr, ast.Constant):
-    """A constant value.
-
-    The `value` attribute contains the Python object it represents.
-    types supported: NoneType, str, bytes, bool, int, float
-    """
-
-    def __init__(self, value: float, kind: str | None = None, **kwargs: Any) -> None:
-        if not isinstance(value, (float, int)):
-            raise TypeError("Constants must be numbers")
-        super().__init__(value, kind, **kwargs)
-
-
-class Index(Expr, ast.Name):
-    """A solution index.
-
-    `id` holds the index as a string (becuase ast.Name requires a string).
-    """
-
-    def __init__(self, index: int) -> None:
-        self.index = index
-        super().__init__(str(index), ctx=ast.Load())
-
-
-op_map: dict[type[ast.cmpop], Relation] = {
-    ast.LtE: Relation.LessEqual,
-    ast.Eq: Relation.Equal,
-    ast.Gt: Relation.GreaterEqual,
-}
-
-
-def to_constraint(expr: Expr) -> LinearConstraint:
-    constraint = LinearConstraint()
-
-    seen_compare = False
-    for sub in ast.walk(expr):
-        if not isinstance(sub, Compare):
-            continue
-        if seen_compare:
-            raise ValueError("Only single comparisons are supported")
-
-        op_type = type(sub.ops[0])
-        try:
-            constraint.set_relation(op_map[op_type])
-        except KeyError as e:
-            raise ValueError(f"Unsupported comparison operator: {op_type}") from e
-
-        seen_compare = True
-
-    for index, coeff in get_coefficients(expr).items():
-        if index is None:
-            constraint.set_value(-coeff)
-        else:
-            constraint.set_coefficient(index, coeff)
-    return constraint
-
-
-def get_coefficients(
-    expr: ast.expr, coeffs: dict[int | None, float] | None = None, scale: int = 1
-) -> dict[int | None, float]:
-    if coeffs is None:
-        coeffs = {}
-
-    if isinstance(expr, Compare):
-        if len(expr.ops) != 1:
-            raise ValueError("Only single comparisons are supported")
-        get_coefficients(expr.left, coeffs)
-        get_coefficients(expr.comparators[0], coeffs, -1)
-
-    elif isinstance(expr, BinOp):
-        if isinstance(expr.op, (ast.Mult, ast.Div)):
-            if isinstance(expr.right, Constant):
-                e = expr.left
-                v = expr.right.value
-            elif isinstance(expr.left, Constant):
-                e = expr.right
-                v = expr.left.value
-            else:
-                raise ValueError("Multiplication must be by a constant")
-            assert isinstance(e, Index)
-            scale *= 1 / v if isinstance(expr.op, ast.Div) else v
-            get_coefficients(e, coeffs, scale)
-
-        else:
-            get_coefficients(expr.left, coeffs, scale)
-            if isinstance(expr.op, (ast.USub, ast.Sub)):
-                scale = -scale
-            get_coefficients(expr.right, coeffs, scale)
-    elif isinstance(expr, UnaryOp):
-        if isinstance(expr.op, ast.USub):
-            scale = -scale
-        get_coefficients(expr.operand, coeffs, scale)
-    elif isinstance(expr, Constant):
-        coeffs[None] = expr.value * scale
-    elif isinstance(expr, Index):
-        coeffs.setdefault(expr.index, 0)
-        coeffs[expr.index] += scale
-    else:
-        raise ValueError("Unsupported expression")
-
-    return coeffs
-
-
-# -u + 2*e - (v + 4 - u)
diff --git a/motile/solver.py b/motile/solver.py
index 1aef073..dcbd643 100644
--- a/motile/solver.py
+++ b/motile/solver.py
@@ -1,7 +1,7 @@
 from __future__ import annotations
 
 import logging
-from typing import TYPE_CHECKING, TypeVar, cast
+from typing import TYPE_CHECKING, Hashable, TypeVar, cast
 
 import ilpy
 import numpy as np
@@ -156,6 +156,19 @@ def solve(self, timeout: float = 0.0, num_threads: int = 1) -> ilpy.Solution:
 
         return self.solution
 
+    def var(self, cls: type[V], name: str = ""):
+        from ilpy import expressions
+
+        class GlobalVariable:
+            def __init__(self, vars: Variable, name: str = ""):
+                self.vars = vars
+                self.name = name
+
+            def expr(self, key: Hashable, name: str = "") -> expressions.Variable:
+                return expressions.Variable(name or self.name, index=self.vars[key])
+
+        return GlobalVariable(self.get_variables(cls), name or cls.__name__)
+
     def get_variables(self, cls: type[V]) -> V:
         """Get variables by their class name.
 
diff --git a/motile/variables/node_appear.py b/motile/variables/node_appear.py
index 9c953c3..6586902 100644
--- a/motile/variables/node_appear.py
+++ b/motile/variables/node_appear.py
@@ -46,20 +46,13 @@ def instantiate_constraints(solver: Solver) -> list[ilpy.LinearConstraint]:
         constraints = []
         for node in solver.graph.nodes:
             prev_edges = solver.graph.prev_edges[node]
-            num_prev_edges = len(prev_edges)
-
-            if num_prev_edges == 0:
+            selected = node_indicators.expr(node, "n")
+            appear = appear_indicators.expr(node, "a")
 
+            if not prev_edges:
                 # special case: no incoming edges, appear indicator is equal to
                 # selection indicator
-                constraint = ilpy.LinearConstraint()
-                constraint.set_coefficient(node_indicators[node], 1.0)
-                constraint.set_coefficient(appear_indicators[node], -1.0)
-                constraint.set_relation(ilpy.Relation.Equal)
-                constraint.set_value(0.0)
-
-                constraints.append(constraint)
-
+                constraints.append((selected == appear).constraint())
                 continue
 
             # Ensure that the following holds:
@@ -72,36 +65,10 @@ def instantiate_constraints(solver: Solver) -> list[ilpy.LinearConstraint]:
             # let s = num_prev * selected - sum(prev_selected)
             # (1) s - appear <= num_prev - 1
             # (2) s - appear * num_prev >= 0
-
-            constraint1 = ilpy.LinearConstraint()
-            constraint2 = ilpy.LinearConstraint()
-
-            # set s for both constraints:
-
-            # num_prev * selected
-            constraint1.set_coefficient(node_indicators[node], num_prev_edges)
-            constraint2.set_coefficient(node_indicators[node], num_prev_edges)
-
-            # - sum(prev_selected)
-            for prev_edge in prev_edges:
-                constraint1.set_coefficient(edge_indicators[prev_edge], -1.0)
-                constraint2.set_coefficient(edge_indicators[prev_edge], -1.0)
-
-            # constraint specific parts:
-
-            # - appear
-            constraint1.set_coefficient(appear_indicators[node], -1.0)
-
-            # - appear * num_prev
-            constraint2.set_coefficient(appear_indicators[node], -num_prev_edges)
-
-            constraint1.set_relation(ilpy.Relation.LessEqual)
-            constraint2.set_relation(ilpy.Relation.GreaterEqual)
-
-            constraint1.set_value(num_prev_edges - 1)
-            constraint2.set_value(0)
-
-            constraints.append(constraint1)
-            constraints.append(constraint2)
+            num_prev = len(prev_edges)
+            s = num_prev * selected - sum(edge_indicators.expr(e) for e in prev_edges)
+            expr1 = s - appear <= num_prev - 1
+            expr2 = s - appear >= 0
+            constraints.extend([expr1.constraint(), expr2.constraint()])
 
         return constraints
diff --git a/motile/variables/variable.py b/motile/variables/variable.py
index bd7fb0c..e8f3f2f 100644
--- a/motile/variables/variable.py
+++ b/motile/variables/variable.py
@@ -4,8 +4,7 @@
 from typing import TYPE_CHECKING, Hashable, Iterable, Sequence
 
 import ilpy
-
-from motile.expressions import Index
+from ilpy import expressions
 
 if TYPE_CHECKING:
     from motile.solver import Solver
@@ -125,5 +124,7 @@ def keys(self) -> Iterable[Hashable]:
     def values(self) -> Iterable[int]:
         return self._index_map.values()
 
-    def expr(self, key: Hashable) -> Index:
-        return Index(self._index_map[key])
+    def expr(self, key: Hashable, name: str = "") -> expressions.Variable:
+        if not name:
+            name = f"{type(self).__name__}({key})"
+        return expressions.Variable(name, index=self._index_map[key])

From 6d0c8c6147e08d879d67d9124ecefbf3f61e4c47 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Fri, 17 Mar 2023 10:02:37 -0400
Subject: [PATCH 05/14] 2 more constraints

---
 motile/constraints/max_children.py | 13 +------------
 motile/constraints/max_parents.py  | 16 ++++++----------
 2 files changed, 7 insertions(+), 22 deletions(-)

diff --git a/motile/constraints/max_children.py b/motile/constraints/max_children.py
index f780d5a..cd764be 100644
--- a/motile/constraints/max_children.py
+++ b/motile/constraints/max_children.py
@@ -37,21 +37,10 @@ def instantiate(self, solver: Solver) -> list[ilpy.LinearConstraint]:
         constraints = []
         for node in solver.graph.nodes:
             n_edges = sum(
-                [edge_indicators.expr(e) for e in solver.graph.next_edges[node]],
+                (edge_indicators.expr(e) for e in solver.graph.next_edges[node]),
                 Constant(0),
             )
             expr = n_edges <= self.max_children
             constraints.append(expr.constraint())
 
-            # constraint = ilpy.LinearConstraint()
-            # # all outgoing edges
-            # for edge in solver.graph.next_edges[node]:
-            #     constraint.set_coefficient(edge_indicators[edge], 1)
-
-            # relation, value
-            # constraint.set_relation(ilpy.Relation.LessEqual)
-
-            # constraint.set_value(self.max_children)
-            # constraints.append(constraint)
-
         return constraints
diff --git a/motile/constraints/max_parents.py b/motile/constraints/max_parents.py
index 38893f8..5bf6f8f 100644
--- a/motile/constraints/max_parents.py
+++ b/motile/constraints/max_parents.py
@@ -3,6 +3,7 @@
 from typing import TYPE_CHECKING
 
 import ilpy
+from ilpy.expressions import Constant
 
 from ..variables import EdgeSelected
 from .constraint import Constraint
@@ -35,16 +36,11 @@ def instantiate(self, solver: Solver) -> list[ilpy.LinearConstraint]:
 
         constraints = []
         for node in solver.graph.nodes:
-            constraint = ilpy.LinearConstraint()
-
             # all incoming edges
-            for edge in solver.graph.prev_edges[node]:
-                constraint.set_coefficient(edge_indicators[edge], 1)
-
-            # relation, value
-            constraint.set_relation(ilpy.Relation.LessEqual)
-
-            constraint.set_value(self.max_parents)
-            constraints.append(constraint)
+            s = sum(
+                (edge_indicators.expr(e) for e in solver.graph.prev_edges[node]),
+                start=Constant(0),
+            )
+            constraints.append((s <= self.max_parents).constraint())
 
         return constraints

From 93c2cc1b6caf3a9e1e56b30a0200a075aa200366 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Fri, 17 Mar 2023 15:33:53 -0400
Subject: [PATCH 06/14] yield expressions

---
 motile/constraints/constraint.py   |  5 ++++-
 motile/constraints/max_children.py |  9 ++++-----
 motile/constraints/max_parents.py  |  7 +++----
 motile/solver.py                   |  5 +++++
 motile/variables/variable.py       | 10 ++++++----
 5 files changed, 22 insertions(+), 14 deletions(-)

diff --git a/motile/constraints/constraint.py b/motile/constraints/constraint.py
index 721cf54..c4ef1e5 100644
--- a/motile/constraints/constraint.py
+++ b/motile/constraints/constraint.py
@@ -5,13 +5,16 @@
 
 if TYPE_CHECKING:
     import ilpy
+    from ilpy.expressions import Expression
 
     from motile.solver import Solver
 
 
 class Constraint(ABC):
     @abstractmethod
-    def instantiate(self, solver: Solver) -> Iterable[ilpy.LinearConstraint]:
+    def instantiate(
+        self, solver: Solver
+    ) -> Iterable[ilpy.LinearConstraint | Expression]:
         """Create and return specific linear constraints for the given solver.
 
         Args:
diff --git a/motile/constraints/max_children.py b/motile/constraints/max_children.py
index 413ba0c..f68285f 100644
--- a/motile/constraints/max_children.py
+++ b/motile/constraints/max_children.py
@@ -2,8 +2,7 @@
 
 from typing import TYPE_CHECKING, Iterable
 
-import ilpy
-from ilpy.expressions import Constant
+from ilpy.expressions import Constant, Expression
 
 from ..variables import EdgeSelected
 from .constraint import Constraint
@@ -31,7 +30,7 @@ class MaxChildren(Constraint):
     def __init__(self, max_children: int) -> None:
         self.max_children = max_children
 
-    def instantiate(self, solver: Solver) -> Iterable[ilpy.LinearConstraint]:
+    def instantiate(self, solver: Solver) -> Iterable[Expression]:
         edge_indicators = solver.get_variables(EdgeSelected)
 
         for node in solver.graph.nodes:
@@ -39,5 +38,5 @@ def instantiate(self, solver: Solver) -> Iterable[ilpy.LinearConstraint]:
                 (edge_indicators.expr(e) for e in solver.graph.next_edges[node]),
                 Constant(0),
             )
-            expr = n_edges <= self.max_children
-            yield expr.constraint()
+
+            yield n_edges <= self.max_children
diff --git a/motile/constraints/max_parents.py b/motile/constraints/max_parents.py
index 96da71c..eb582a4 100644
--- a/motile/constraints/max_parents.py
+++ b/motile/constraints/max_parents.py
@@ -2,8 +2,7 @@
 
 from typing import TYPE_CHECKING, Iterable
 
-import ilpy
-from ilpy.expressions import Constant
+from ilpy.expressions import Constant, Expression
 
 from ..variables import EdgeSelected
 from .constraint import Constraint
@@ -31,7 +30,7 @@ class MaxParents(Constraint):
     def __init__(self, max_parents: int) -> None:
         self.max_parents = max_parents
 
-    def instantiate(self, solver: Solver) -> Iterable[ilpy.LinearConstraint]:
+    def instantiate(self, solver: Solver) -> Iterable[Expression]:
         edge_indicators = solver.get_variables(EdgeSelected)
 
         for node in solver.graph.nodes:
@@ -40,4 +39,4 @@ def instantiate(self, solver: Solver) -> Iterable[ilpy.LinearConstraint]:
                 (edge_indicators.expr(e) for e in solver.graph.prev_edges[node]),
                 start=Constant(0),
             )
-            yield (s <= self.max_parents).constraint()
+            yield s <= self.max_parents
diff --git a/motile/solver.py b/motile/solver.py
index 75fd64b..15a21cb 100644
--- a/motile/solver.py
+++ b/motile/solver.py
@@ -5,6 +5,7 @@
 
 import ilpy
 import numpy as np
+from ilpy.expressions import Expression
 
 from .constraints import SelectEdgeNodes
 from .constraints.constraint import Constraint
@@ -108,6 +109,8 @@ def add_constraints(self, constraints: Constraint) -> None:
         logger.info("Adding %s constraints...", type(constraints).__name__)
 
         for constraint in constraints.instantiate(self):
+            if isinstance(constraint, Expression):
+                constraint = constraint.constraint()
             self.constraints.add(constraint)
 
     def solve(self, timeout: float = 0.0, num_threads: int = 1) -> ilpy.Solution:
@@ -248,6 +251,8 @@ def _add_variables(self, cls: type[V]) -> None:
             self.variable_types[index] = cls.variable_type
 
         for constraint in cls.instantiate_constraints(self):
+            if isinstance(constraint, Expression):
+                constraint = constraint.constraint()
             self.constraints.add(constraint)
 
         self.features.resize(num_variables=self.num_variables)
diff --git a/motile/variables/variable.py b/motile/variables/variable.py
index 75811b9..c375fba 100644
--- a/motile/variables/variable.py
+++ b/motile/variables/variable.py
@@ -90,7 +90,9 @@ def instantiate(solver):
         pass
 
     @staticmethod
-    def instantiate_constraints(solver: Solver) -> Iterable[ilpy.LinearConstraint]:
+    def instantiate_constraints(
+        solver: Solver,
+    ) -> Iterable[ilpy.LinearConstraint | expressions.Expression]:
         """Add linear constraints to the solver to ensure that these variables
         are coupled to other variables of the solver.
 
@@ -101,9 +103,9 @@ def instantiate_constraints(solver: Solver) -> Iterable[ilpy.LinearConstraint]:
 
         Returns:
 
-            A iterable of :class:`ilpy.LinearConstraint`. See
-            :class:`motile.constraints.Constraint` for how to create linear
-            constraints.
+            A iterable of :class:`ilpy.LinearConstraint` or
+            :class:`ilpy.expressions.Expression.` See
+            :class:`motile.constraints.Constraint` for how to create linear constraints.
         """
         return []
 

From 447cdfee4705df4819c9b2b5112b9e55de5f5d39 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Fri, 17 Mar 2023 16:05:55 -0400
Subject: [PATCH 07/14] cleanup

---
 motile/constraints/select_edge_nodes.py |  9 ++++-----
 motile/variables/node_appear.py         | 13 +++++++------
 2 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/motile/constraints/select_edge_nodes.py b/motile/constraints/select_edge_nodes.py
index 7e27dea..68755c9 100644
--- a/motile/constraints/select_edge_nodes.py
+++ b/motile/constraints/select_edge_nodes.py
@@ -2,12 +2,12 @@
 
 from typing import TYPE_CHECKING, Iterable
 
-import ilpy
-
 from ..variables import EdgeSelected, NodeSelected
 from .constraint import Constraint
 
 if TYPE_CHECKING:
+    from ilpy.expressions import Expression
+
     from motile.solver import Solver
 
 
@@ -24,7 +24,7 @@ class SelectEdgeNodes(Constraint):
     This constraint will be added by default to any :class:`Solver` instance.
     """
 
-    def instantiate(self, solver: Solver) -> Iterable[ilpy.LinearConstraint]:
+    def instantiate(self, solver: Solver) -> Iterable[Expression]:
         node_indicators = solver.get_variables(NodeSelected)
         edge_indicators = solver.get_variables(EdgeSelected)
 
@@ -32,5 +32,4 @@ def instantiate(self, solver: Solver) -> Iterable[ilpy.LinearConstraint]:
             nodes = list(solver.graph.nodes_of(edge))
             x_e = edge_indicators.expr(edge, "e")
 
-            expr = len(nodes) * x_e - sum(node_indicators.expr(n) for n in nodes) <= 0
-            yield expr.constraint()
+            yield len(nodes) * x_e - sum(node_indicators.expr(n) for n in nodes) <= 0
diff --git a/motile/variables/node_appear.py b/motile/variables/node_appear.py
index 16bdbd9..30e6343 100644
--- a/motile/variables/node_appear.py
+++ b/motile/variables/node_appear.py
@@ -2,13 +2,13 @@
 
 from typing import TYPE_CHECKING, Collection, Iterable
 
-import ilpy
-
 from .edge_selected import EdgeSelected
 from .node_selected import NodeSelected
 from .variable import Variable
 
 if TYPE_CHECKING:
+    from ilpy.expressions import Expression
+
     from motile._types import NodeId
     from motile.solver import Solver
 
@@ -39,7 +39,7 @@ def instantiate(solver: Solver) -> Collection[NodeId]:
         return solver.graph.nodes
 
     @staticmethod
-    def instantiate_constraints(solver: Solver) -> Iterable[ilpy.LinearConstraint]:
+    def instantiate_constraints(solver: Solver) -> Iterable[Expression]:
         appear_indicators = solver.get_variables(NodeAppear)
         node_indicators = solver.get_variables(NodeSelected)
         edge_indicators = solver.get_variables(EdgeSelected)
@@ -52,7 +52,7 @@ def instantiate_constraints(solver: Solver) -> Iterable[ilpy.LinearConstraint]:
             if not prev_edges:
                 # special case: no incoming edges, appear indicator is equal to
                 # selection indicator
-                yield (selected == appear).constraint()
+                yield selected == appear
                 continue
 
             # Ensure that the following holds:
@@ -65,8 +65,9 @@ def instantiate_constraints(solver: Solver) -> Iterable[ilpy.LinearConstraint]:
             # let s = num_prev * selected - sum(prev_selected)
             # (1) s - appear <= num_prev - 1
             # (2) s - appear * num_prev >= 0
+
             num_prev = len(prev_edges)
             s = num_prev * selected - sum(edge_indicators.expr(e) for e in prev_edges)
 
-            yield (s - appear <= num_prev - 1).constraint()
-            yield (s - appear >= 0).constraint()
+            yield s - appear <= num_prev - 1
+            yield s - appear >= 0

From 3fa0f655486745ace9ae321e527ef0d5068ba996 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Wed, 3 May 2023 10:46:06 -0400
Subject: [PATCH 08/14] fix test

---
 motile/solver.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/motile/solver.py b/motile/solver.py
index ec45bd3..cec014c 100644
--- a/motile/solver.py
+++ b/motile/solver.py
@@ -110,7 +110,7 @@ def add_constraints(self, constraints: Constraint) -> None:
 
         for constraint in constraints.instantiate(self):
             if isinstance(constraint, Expression):
-                constraint = constraint.constraint()
+                constraint = constraint.as_constraint()
             self.constraints.add(constraint)
 
     def solve(self, timeout: float = 0.0, num_threads: int = 1) -> ilpy.Solution:
@@ -254,7 +254,7 @@ def _add_variables(self, cls: type[V]) -> None:
 
         for constraint in cls.instantiate_constraints(self):
             if isinstance(constraint, Expression):
-                constraint = constraint.constraint()
+                constraint = constraint.as_constraint()
             self.constraints.add(constraint)
 
         self.features.resize(num_variables=self.num_variables)

From 744dbe258a1935f84ec7ef0f42a7d086d7db4378 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Sat, 6 May 2023 12:54:14 -0400
Subject: [PATCH 09/14] indexing motile variables  returns ilpy.Variable

---
 motile/constraints/max_children.py      |  3 +--
 motile/constraints/max_parents.py       |  2 +-
 motile/constraints/select_edge_nodes.py |  5 ++---
 motile/solver.py                        |  5 -----
 motile/variables/node_appear.py         |  6 +++---
 motile/variables/variable.py            | 15 +++++++--------
 6 files changed, 14 insertions(+), 22 deletions(-)

diff --git a/motile/constraints/max_children.py b/motile/constraints/max_children.py
index f68285f..e22d4be 100644
--- a/motile/constraints/max_children.py
+++ b/motile/constraints/max_children.py
@@ -35,8 +35,7 @@ def instantiate(self, solver: Solver) -> Iterable[Expression]:
 
         for node in solver.graph.nodes:
             n_edges = sum(
-                (edge_indicators.expr(e) for e in solver.graph.next_edges[node]),
-                Constant(0),
+                (edge_indicators[e] for e in solver.graph.next_edges[node]), Constant(0)
             )
 
             yield n_edges <= self.max_children
diff --git a/motile/constraints/max_parents.py b/motile/constraints/max_parents.py
index eb582a4..d26b13e 100644
--- a/motile/constraints/max_parents.py
+++ b/motile/constraints/max_parents.py
@@ -36,7 +36,7 @@ def instantiate(self, solver: Solver) -> Iterable[Expression]:
         for node in solver.graph.nodes:
             # all incoming edges
             s = sum(
-                (edge_indicators.expr(e) for e in solver.graph.prev_edges[node]),
+                (edge_indicators[e] for e in solver.graph.prev_edges[node]),
                 start=Constant(0),
             )
             yield s <= self.max_parents
diff --git a/motile/constraints/select_edge_nodes.py b/motile/constraints/select_edge_nodes.py
index 68755c9..a9c6488 100644
--- a/motile/constraints/select_edge_nodes.py
+++ b/motile/constraints/select_edge_nodes.py
@@ -30,6 +30,5 @@ def instantiate(self, solver: Solver) -> Iterable[Expression]:
 
         for edge in solver.graph.edges:
             nodes = list(solver.graph.nodes_of(edge))
-            x_e = edge_indicators.expr(edge, "e")
-
-            yield len(nodes) * x_e - sum(node_indicators.expr(n) for n in nodes) <= 0
+            x_e = edge_indicators[edge]
+            yield len(nodes) * x_e - sum(node_indicators[n] for n in nodes) <= 0
diff --git a/motile/solver.py b/motile/solver.py
index cec014c..fccc6d7 100644
--- a/motile/solver.py
+++ b/motile/solver.py
@@ -5,7 +5,6 @@
 
 import ilpy
 import numpy as np
-from ilpy.expressions import Expression
 
 from .constraints import SelectEdgeNodes
 from .constraints.constraint import Constraint
@@ -109,8 +108,6 @@ def add_constraints(self, constraints: Constraint) -> None:
         logger.info("Adding %s constraints...", type(constraints).__name__)
 
         for constraint in constraints.instantiate(self):
-            if isinstance(constraint, Expression):
-                constraint = constraint.as_constraint()
             self.constraints.add(constraint)
 
     def solve(self, timeout: float = 0.0, num_threads: int = 1) -> ilpy.Solution:
@@ -253,8 +250,6 @@ def _add_variables(self, cls: type[V]) -> None:
             self.variable_types[index] = cls.variable_type
 
         for constraint in cls.instantiate_constraints(self):
-            if isinstance(constraint, Expression):
-                constraint = constraint.as_constraint()
             self.constraints.add(constraint)
 
         self.features.resize(num_variables=self.num_variables)
diff --git a/motile/variables/node_appear.py b/motile/variables/node_appear.py
index 0c44459..f513114 100644
--- a/motile/variables/node_appear.py
+++ b/motile/variables/node_appear.py
@@ -46,8 +46,8 @@ def instantiate_constraints(solver: Solver) -> Iterable[ilpy.Expression]:
 
         for node in solver.graph.nodes:
             prev_edges = solver.graph.prev_edges[node]
-            selected = node_indicators.expr(node, "n")
-            appear = appear_indicators.expr(node, "a")
+            selected = node_indicators[node]
+            appear = appear_indicators[node]
 
             if not prev_edges:
                 # special case: no incoming edges, appear indicator is equal to
@@ -67,7 +67,7 @@ def instantiate_constraints(solver: Solver) -> Iterable[ilpy.Expression]:
             # (2) s - appear * num_prev >= 0
 
             num_prev = len(prev_edges)
-            s = num_prev * selected - sum(edge_indicators.expr(e) for e in prev_edges)
+            s = num_prev * selected - sum(edge_indicators[e] for e in prev_edges)
 
             yield s - appear <= num_prev - 1
             yield s - appear >= 0
diff --git a/motile/variables/variable.py b/motile/variables/variable.py
index 4ee10e8..c6c2c98 100644
--- a/motile/variables/variable.py
+++ b/motile/variables/variable.py
@@ -20,7 +20,7 @@
 _KT = TypeVar("_KT", bound=Hashable)
 
 
-class Variable(ABC, Mapping[_KT, int]):
+class Variable(ABC, Mapping[_KT, ilpy.Variable]):
     """Base class for solver variables.
 
     New variables can be introduced by inheriting from this base class and
@@ -38,6 +38,9 @@ class Variable(ABC, Mapping[_KT, int]):
 
         solution = solver.solve()
 
+        # here `node_selected` is an instance of a Variable subclass
+        # specifically, it will be an instance of NodeSelected, which
+        # maps node Ids to variables in the solver.
         node_selected = solver.get_variables(NodeSelected)
 
         for node in graph.nodes:
@@ -127,8 +130,9 @@ def __repr__(self) -> str:
             rs.append(r)
         return "\n".join(rs)
 
-    def __getitem__(self, key: _KT) -> int:
-        return self._index_map[key]
+    def __getitem__(self, key: _KT) -> ilpy.Variable:
+        name = f"{type(self).__name__}({key})"
+        return ilpy.Variable(name, index=self._index_map[key])
 
     def __iter__(self) -> Iterator[_KT]:
         return iter(self._index_map)
@@ -136,11 +140,6 @@ def __iter__(self) -> Iterator[_KT]:
     def __len__(self) -> int:
         return len(self._index_map)
 
-    def expr(self, key: _KT, name: str = "") -> ilpy.Variable:
-        if not name:
-            name = f"{type(self).__name__}({key})"
-        return ilpy.Variable(name, index=self._index_map[key])
-
     # All of these methods are provided by subclassing typing.Mapping
     # __contains__
     # keys

From 01ae3bb16bb7f391bf4c9627728c73682dc3acf2 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Sat, 6 May 2023 13:42:23 -0400
Subject: [PATCH 10/14] add ilpy to pre-commit typing

---
 .pre-commit-config.yaml  | 2 ++
 motile/costs/features.py | 8 +++++++-
 motile/solver.py         | 4 +++-
 3 files changed, 12 insertions(+), 2 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 321f9da..fb62ced 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -20,3 +20,5 @@ repos:
     hooks:
       - id: mypy
         files: "^motile/"
+        additional_dependencies:
+          - ilpy
diff --git a/motile/costs/features.py b/motile/costs/features.py
index f61f5e4..dcefaea 100644
--- a/motile/costs/features.py
+++ b/motile/costs/features.py
@@ -1,7 +1,12 @@
 from __future__ import annotations
 
+from typing import TYPE_CHECKING
+
 import numpy as np
 
+if TYPE_CHECKING:
+    import ilpy
+
 
 class Features:
     def __init__(self) -> None:
@@ -33,10 +38,11 @@ def _increase_features(self, num_features: int) -> None:
         self._values = np.hstack((self._values, new_features))
 
     def add_feature(
-        self, variable_index: int, feature_index: int, value: float
+        self, variable_index: int | ilpy.Variable, feature_index: int, value: float
     ) -> None:
         num_variables, num_features = self._values.shape
 
+        variable_index = int(variable_index)
         if variable_index >= num_variables or feature_index >= num_features:
             self.resize(
                 max(variable_index + 1, num_variables),
diff --git a/motile/solver.py b/motile/solver.py
index fccc6d7..43c76cf 100644
--- a/motile/solver.py
+++ b/motile/solver.py
@@ -180,7 +180,9 @@ def get_variables(self, cls: type[V]) -> V:
             self._add_variables(cls)
         return cast("V", self.variables[cls])
 
-    def add_variable_cost(self, index: int, value: float, weight: Weight) -> None:
+    def add_variable_cost(
+        self, index: int | ilpy.Variable, value: float, weight: Weight
+    ) -> None:
         """Add costs for an individual variable.
 
         To be used within implementations of :class:`motile.costs.Costs`.

From 38e1b234b0553e29c7aa4a64076fdf148eaae8a1 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Sat, 6 May 2023 14:24:01 -0400
Subject: [PATCH 11/14] bump ilpy

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 7381183..1793742 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,7 +17,7 @@ authors = [
     { name = 'Florian Jug', email = 'florian.jug@fht.org' },
 ]
 dynamic = ["version"]
-dependencies = ['networkx', 'ilpy>=0.3.0', 'numpy', 'structsvm']
+dependencies = ['networkx', 'ilpy>=0.3.1', 'numpy', 'structsvm']
 
 [project.optional-dependencies]
 dev = ["pre-commit", "pytest", "pytest-cov", "ruff", "twine", "build"]

From 5c0c7d393dbf37bf7df54f6f901e77eb03318cb8 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Sat, 6 May 2023 14:25:19 -0400
Subject: [PATCH 12/14] bump pre-commit dep

---
 .pre-commit-config.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index fb62ced..d7cdb60 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -21,4 +21,4 @@ repos:
       - id: mypy
         files: "^motile/"
         additional_dependencies:
-          - ilpy
+          - ilpy>=0.3.1

From 4a515d6942e56859af972339132cdd725442c632 Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Sat, 6 May 2023 14:31:57 -0400
Subject: [PATCH 13/14] remove ilpy from pre-commit

---
 .pre-commit-config.yaml | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index d7cdb60..b8a2f53 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -20,5 +20,6 @@ repos:
     hooks:
       - id: mypy
         files: "^motile/"
-        additional_dependencies:
-          - ilpy>=0.3.1
+        # this won't work cause pre-commit can't install ilpy on CI
+        # additional_dependencies:
+          # - ilpy>=0.3.1
\ No newline at end of file

From 937b08591266a8bd09960e24d9f537242e9694bb Mon Sep 17 00:00:00 2001
From: Talley Lambert <talley.lambert@gmail.com>
Date: Sat, 6 May 2023 14:34:07 -0400
Subject: [PATCH 14/14] undo changes to pre-commit

---
 .pre-commit-config.yaml | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index b8a2f53..321f9da 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -20,6 +20,3 @@ repos:
     hooks:
       - id: mypy
         files: "^motile/"
-        # this won't work cause pre-commit can't install ilpy on CI
-        # additional_dependencies:
-          # - ilpy>=0.3.1
\ No newline at end of file