Repository: matthewwardrop/formulaic

Issue:
Null rows that are dropped are not consistent across left and right hand sides
e.g. 
```
from formulaic import Formula
import pandas

Formula("y ~ x").get_model_matrix(pandas.DataFrame({'y': [1,2,3,4,None], 'x': [None, None, 1,2,3]}))
```

outputs:

```
(     y
 0  1.0
 1  2.0
 2  3.0
 3  4.0,
    Intercept    x
 2        1.0  1.0
 3        1.0  2.0
 4        1.0  3.0)
```

Discussion:

Patch:
diff --git a/formulaic/formula.py b/formulaic/formula.py
index 1c952c4..8a0f5bf 100644
--- a/formulaic/formula.py
+++ b/formulaic/formula.py
@@ -3,7 +3,7 @@ import inspect
 from .errors import FormulaInvalidError, FormulaMaterializerInvalidError
 from .materializers.base import FormulaMaterializer
 from .parser import FormulaParser
-from .parser.types import Term
+from .parser.types import Structured, Term
 from .utils.calculus import differentiate_term
 
 
@@ -49,7 +49,9 @@ class Formula:
 
     @classmethod
     def __check_terms(cls, terms, depth=0):
-        if depth == 0 and isinstance(terms, tuple):
+        if isinstance(terms, Structured):
+            terms._map(cls.__check_terms)
+        elif depth == 0 and isinstance(terms, tuple):
             for termset in terms:
                 cls.__check_terms(termset, depth=depth + 1)
         else:
@@ -81,12 +83,32 @@ class Formula:
             [differentiate_term(term, vars, use_sympy=use_sympy) for term in self.terms]
         )
 
+    def __getattr__(self, attr):
+        if attr in self.terms._to_dict(recurse=False):
+            return Formula(self.terms[attr])
+        raise AttributeError(f"This formula has no substructures keyed by '{attr}'.")
+
+    def __getitem__(self, item):
+        if (
+            isinstance(self.terms, tuple)
+            or isinstance(self.terms, Structured)
+            and self.terms._has_root
+            and isinstance(self.terms.root, tuple)
+        ):
+            return Formula(self.terms[item])
+        raise KeyError(
+            f"This formula does not have any sub-parts indexable via `{repr(item)}`."
+        )
+
     def __str__(self):
+        terms = self.terms
         if isinstance(self.terms, tuple):
-            return " ~ ".join(
-                " + ".join(str(term) for term in terms) for terms in self._terms
+            terms = Structured(terms)
+        if isinstance(terms, Structured):
+            return str(
+                terms._map(lambda terms: " + ".join(str(term) for term in terms))
             )
-        return " + ".join(str(term) for term in self.terms)
+        return " + ".join(str(term) for term in terms)
 
     def __eq__(self, other):
         if isinstance(other, Formula):
diff --git a/formulaic/materializers/base.py b/formulaic/materializers/base.py
index c502ebc..3b01531 100644
--- a/formulaic/materializers/base.py
+++ b/formulaic/materializers/base.py
@@ -27,13 +27,12 @@ from formulaic.errors import (
 )
 from formulaic.materializers.types.factor_values import FactorValuesMetadata
 from formulaic.model_matrix import ModelMatrix
+from formulaic.parser.types import Factor, Structured, Term
+from formulaic.transforms import TRANSFORMS
 from formulaic.utils.cast import as_columns
 from formulaic.utils.layered_mapping import LayeredMapping
 from formulaic.utils.stateful_transforms import stateful_eval
 
-from formulaic.parser.types import Factor, Term
-from formulaic.transforms import TRANSFORMS
-
 from .types import EvaluatedFactor, FactorValues, ScopedFactor, ScopedTerm
 
 if TYPE_CHECKING:
@@ -127,64 +126,82 @@ class FormulaMaterializer(metaclass=FormulaMaterializerMeta):
     def get_model_matrix(
         self, spec, ensure_full_rank=True, na_action="drop", output=None
     ):
-        from formulaic.formula import Formula
         from formulaic.model_spec import ModelSpec
 
-        # Prepare ModelSpec instance
+        # Prepare arguments
         if output is None:
             output = self.REGISTER_OUTPUTS[0]
-
-        if isinstance(spec, Formula):
-            spec = ModelSpec(
-                formula=spec,
-                materializer=self,
-                ensure_full_rank=ensure_full_rank,
-                na_action=na_action,
-                output=output,
-            )
-        if not isinstance(spec, ModelSpec):
-            spec = ModelSpec(
-                formula=Formula(spec),
-                materializer=self,
-                ensure_full_rank=ensure_full_rank,
-                na_action=na_action,
-                output=output,
-            )
-
         if output not in self.REGISTER_OUTPUTS:
             raise FormulaMaterializationError(
                 f"Nominated output {repr(output)} is invalid. Available output types are: {set(self.REGISTER_OUTPUTS)}."
             )
 
-        # Step 0: Check whether formula separators are in play, and if so, recurse.
-        if isinstance(spec.formula.terms, tuple):
-            return tuple(
-                self.get_model_matrix(
-                    Formula(terms),
-                    ensure_full_rank=ensure_full_rank,
-                    na_action=na_action,
-                    output=output,
-                )
-                for terms in spec.formula.terms
+        # Prepare (potentially structured) ModelSpec(s)
+        model_specs = self._prepare_model_specs(
+            spec, ensure_full_rank=ensure_full_rank, na_action=na_action, output=output
+        )
+        if not isinstance(model_specs, Structured):
+            model_specs = Structured(model_specs)
+
+        # Step 0: Pool all factors and transform state, ensuring consistency
+        # across model specs
+        factors = set()
+        transform_state = {}
+        model_specs._map(
+            lambda ms: factors.update(
+                itertools.chain(*(term.factors for term in ms.formula.terms))
+            )
+        )
+        model_specs._map(
+            lambda ms: transform_state.update(ms.transform_state)
+        )  # TODO: Enforce consistency
+
+        # Step 1: Evaluate all factors and cache the results, keeping track of
+        # which rows need dropping (if `self.config.na_action == 'drop'`).
+        drop_rows = set()
+        dummy_model_spec = ModelSpec(
+            formula=[],
+            ensure_full_rank=ensure_full_rank,
+            na_action=na_action,
+            output=output,
+            transform_state=transform_state,
+        )
+        for factor in factors:
+            self._evaluate_factor(factor, dummy_model_spec, drop_rows)
+        drop_rows = sorted(drop_rows)
+
+        # Step 2: Update the structured model specs with the information from
+        # the shared transform state pool.
+        model_specs._map(
+            lambda ms: ms.transform_state.update(
+                {
+                    factor.expr: transform_state[factor.expr]
+                    for term in ms.formula.terms
+                    for factor in term.factors
+                    if factor.expr in transform_state
+                }
             )
+        )
 
-        drop_rows = (
-            set()
-        )  # Keep track of which rows to drop if `self.config.na_action == 'drop'`.
+        # Step 3: Build the model matrices using the shared factor cache, and
+        # by recursing over the structured model matrices.
+        model_matrices = model_specs._map(
+            lambda model_spec: self._build_model_matrix(model_spec, drop_rows=drop_rows)
+        )
+        model_matrices._self_mapped_attrs = {"model_spec"}
 
-        # Step 1: Evaluate all factors
-        for term in spec.formula.terms:
-            for factor in term.factors:
-                self._evaluate_factor(factor, spec, drop_rows)
+        if len(model_matrices) == 1 and model_matrices._has_root:
+            return model_matrices.root
+        return model_matrices
 
-        drop_rows = sorted(drop_rows)
+    def _build_model_matrix(self, spec: ModelSpec, drop_rows):
 
-        # Step 2: Determine strategy to maintain structural full-rankness of output matrix
+        # Step 1: Determine strategy to maintain structural full-rankness of output matrix
         scoped_terms_for_terms = self._get_scoped_terms(
             spec.formula.terms, ensure_full_rank=spec.ensure_full_rank
         )
 
-        # Step 3: Generate the columns which will be collated into the full matrix
+        # Step 2: Generate the columns which will be collated into the full matrix
         cols = []
         for term, scoped_terms in scoped_terms_for_terms:
             scoped_cols = OrderedDict()
@@ -213,7 +230,7 @@ class FormulaMaterializer(metaclass=FormulaMaterializerMeta):
                     )
             cols.append((term, scoped_terms, scoped_cols))
 
-        # Step 4: Populate remaining model spec fields
+        # Step 3: Populate remaining model spec fields
         spec.materializer = self
         if spec.structure:
             cols = self._enforce_structure(cols, spec, drop_rows)
@@ -227,7 +244,7 @@ class FormulaMaterializer(metaclass=FormulaMaterializerMeta):
                 for term, scoped_terms, scoped_cols in cols
             ]
 
-        # Step 5: Collate factors into one ModelMatrix
+        # Step 4: Collate factors into one ModelMatrix
         return ModelMatrix(
             self._combine_columns(
                 [
@@ -241,6 +258,47 @@ class FormulaMaterializer(metaclass=FormulaMaterializerMeta):
             spec=spec,
         )
 
+    # Methods related to input preparation
+
+    def _prepare_model_specs(
+        self,
+        spec,
+        ensure_full_rank,
+        na_action,
+        output,
+    ) -> Union[ModelSpec, Structured[ModelSpec]]:
+        from formulaic.formula import Formula
+        from formulaic.model_spec import ModelSpec
+
+        if isinstance(spec, Structured):
+            return spec._map(
+                functools.partial(
+                    self._prepare_model_specs,
+                    ensure_full_rank=ensure_full_rank,
+                    na_action=na_action,
+                    output=output,
+                ),
+                recurse=False,
+            )
+        elif isinstance(spec, ModelSpec):
+            return spec
+
+        formula = Formula.from_spec(spec)
+        if isinstance(formula.terms, Structured):
+            return self._prepare_model_specs(
+                formula.terms,
+                ensure_full_rank=ensure_full_rank,
+                na_action=na_action,
+                output=output,
+            )
+        return ModelSpec(
+            formula=formula,
+            materializer=self,
+            ensure_full_rank=ensure_full_rank,
+            na_action=na_action,
+            output=output,
+        )
+
     # Methods related to ensuring out matrices are structurally full-rank
 
     def _get_scoped_terms(self, terms, ensure_full_rank=True):
diff --git a/formulaic/model_spec.py b/formulaic/model_spec.py
index 524609a..8cdb204 100644
--- a/formulaic/model_spec.py
+++ b/formulaic/model_spec.py
@@ -23,8 +23,8 @@ class ModelSpec:
         self.materializer = materializer
         self.na_action = NAAction(na_action)
         self.output = output
-        self.transform_state = transform_state or {}
-        self.encoder_state = encoder_state or {}
+        self.transform_state = transform_state if transform_state is not None else {}
+        self.encoder_state = encoder_state if encoder_state is not None else {}
 
     @property
     def materializer(self):
diff --git a/formulaic/parser/algos/tokens_to_ast.py b/formulaic/parser/algos/tokens_to_ast.py
index 800246c..f4e77cd 100644
--- a/formulaic/parser/algos/tokens_to_ast.py
+++ b/formulaic/parser/algos/tokens_to_ast.py
@@ -87,7 +87,11 @@ def tokens_to_ast(
                     if operator_stack
                     else len(output_queue)
                 )
-                operators = operator_resolver.resolve(token, max_prefix_arity)
+                operators = operator_resolver.resolve(
+                    token,
+                    max_prefix_arity=max_prefix_arity,
+                    context=[s.operator for s in operator_stack],
+                )
 
                 for operator in operators:
 
diff --git a/formulaic/parser/parser.py b/formulaic/parser/parser.py
index 26594f1..d2966ee 100644
--- a/formulaic/parser/parser.py
+++ b/formulaic/parser/parser.py
@@ -2,13 +2,18 @@ import ast
 import itertools
 import functools
 import re
-from typing import Any, Iterable, Union
-
-from formulaic.errors import FormulaParsingError
+from typing import Any, List, Iterable, Union
 
 from .algos.tokens_to_ast import tokens_to_ast
 from .algos.tokenize import tokenize
-from .types import ASTNode, Factor, Operator, OperatorResolver, Term, Token
+from .types import (
+    ASTNode,
+    Operator,
+    OperatorResolver,
+    Structured,
+    Term,
+    Token,
+)
 from .utils import exc_for_token
 
 
@@ -39,6 +44,8 @@ class FormulaParser:
             it will default to `DefaultOperatorResolver`.
     """
 
+    ZERO_PATTERN = re.compile(r"(?:^|(?<=\W))0(?=\W|$)")
+
     def __init__(self, operator_resolver: OperatorResolver = None):
         self.operator_resolver = operator_resolver or DefaultOperatorResolver()
 
@@ -55,50 +62,26 @@ class FormulaParser:
                 (formulas can still omit this intercept in the usual manner:
                 adding a '-1' or '+0' term).
         """
-        tokens = list(tokenize(formula))
 
-        # Insert "1" or "1 + " to beginning of RHS formula
-        one = Token(token="1", kind="value")
-        plus = Token(token="+", kind="operator")
-        minus = Token(token="-", kind="operator")
+        # Transform formula to add intercepts and replace 0 with -1. We do this
+        # as a string transform to reduce the complexity of the code, and also
+        # to avoid the ambiguity in the AST around intentionally unary vs.
+        # incidentally unary operations (e.g. "+0" vs. "x + (+0)").
 
+        # Substitute "0" -> "-1"
+        formula = self.ZERO_PATTERN.sub("-1", formula)
+
+        # Insert ones
         if include_intercept:
-            if len(tokens) == 0:
-                tokens = [one]
-            else:
-                try:
-                    tilde_index = len(tokens) - 1 - tokens[::-1].index("~")
-                    if tilde_index == len(tokens) - 1:
-                        tokens.append(one)
-                    else:
-                        tokens.insert(tilde_index + 1, one)
-                        tokens.insert(tilde_index + 2, plus)
-                except ValueError:
-                    tokens.insert(0, one)
-                    tokens.insert(1, plus)
-
-        # Replace all "0"s with "-1"
-        zero_index = -1
-        try:
-            while True:
-                zero_index = tokens.index("0", zero_index + 1)
-                if zero_index - 1 < 0 or tokens[zero_index - 1] == "~":
-                    tokens.pop(zero_index)
-                    zero_index -= 1
-                    continue
-                elif tokens[zero_index - 1] == "+":
-                    tokens[zero_index - 1] = minus
-                elif tokens[zero_index - 1] == "-":
-                    tokens[zero_index - 1] = plus
-                else:
-                    raise FormulaParsingError(
-                        f"Unrecognised use of `0` at index: {tokens[zero_index-1].source_start}."
-                    )
-                tokens[zero_index] = one
-        except ValueError:
-            pass
+            rhs_offset = formula.index("~") + 1 if "~" in formula else 0
+            formula = formula[:rhs_offset] + "|".join(
+                [
+                    " 1" if len(rhs_part.strip()) == 0 else f" 1 +{rhs_part}"
+                    for rhs_part in formula[rhs_offset:].split("|")
+                ]
+            )
 
-        return tokens
+        return tokenize(formula)
 
     def get_ast(self, formula: str, *, include_intercept: bool = True) -> ASTNode:
         """
@@ -133,7 +116,9 @@ class FormulaParser:
         terms = self.get_ast(formula, include_intercept=include_intercept).to_terms()
 
         if sort:
-            if isinstance(terms, tuple):
+            if isinstance(terms, Structured):
+                terms = terms._map(sorted)
+            elif isinstance(terms, tuple):
                 terms = tuple(sorted(ts) for ts in terms)
             else:
                 terms = sorted(terms)
@@ -154,7 +139,7 @@ class DefaultOperatorResolver(OperatorResolver):
 
     @property
     def operators(self):
-        def formula_separator_expansion(lhs, rhs):
+        def formula_part_expansion(lhs, rhs):
             terms = (lhs.to_terms(), rhs.to_terms())
 
             out = []
@@ -170,17 +155,6 @@ class DefaultOperatorResolver(OperatorResolver):
             common = functools.reduce(lambda x, y: x * y, terms)
             return terms.union({common * term for term in nested.to_terms()})
 
-        def unary_negation(arg):
-            # TODO: FormulaParser().get_terms('a * ( - b)') Should return `a`
-            terms = arg.to_terms()
-            if len(terms) > 1 or list(terms)[0] != "0":
-                raise FormulaParsingError(
-                    "Unary negation is only implemented for '0', where it is substituted for '1'."
-                )
-            return {
-                Term(factors=[Factor("1", eval_method="literal")])
-            }  # pragma: no cover; All zero handling is currently done in the token pre-processor.
-
         def power(arg, power):
             if (
                 not isinstance(power, Token)
@@ -201,7 +175,10 @@ class DefaultOperatorResolver(OperatorResolver):
                 arity=2,
                 precedence=-100,
                 associativity=None,
-                to_terms=formula_separator_expansion,
+                to_terms=lambda lhs, rhs: Structured(
+                    lhs=lhs.to_terms(), rhs=rhs.to_terms()
+                ),
+                accepts_context=lambda context: len(context) == 0,
             ),
             Operator(
                 "~",
@@ -209,7 +186,18 @@ class DefaultOperatorResolver(OperatorResolver):
                 precedence=-100,
                 associativity=None,
                 fixity="prefix",
-                to_terms=lambda expr: (expr.to_terms(),),
+                to_terms=lambda expr: expr.to_terms(),
+                accepts_context=lambda context: len(context) == 0,
+            ),
+            Operator(
+                "|",
+                arity=2,
+                precedence=-50,
+                associativity=None,
+                to_terms=formula_part_expansion,
+                accepts_context=lambda context: all(
+                    isinstance(c, Operator) and c.symbol in "~|" for c in context
+                ),
             ),
             Operator(
                 "+",
@@ -243,7 +231,7 @@ class DefaultOperatorResolver(OperatorResolver):
                 precedence=100,
                 associativity="right",
                 fixity="prefix",
-                to_terms=unary_negation,
+                to_terms=lambda arg: set(),
             ),
             Operator(
                 "*",
@@ -279,22 +267,30 @@ class DefaultOperatorResolver(OperatorResolver):
             ),
         ]
 
-    def resolve(self, token: Token, max_prefix_arity) -> Iterable[Operator]:
+    def resolve(
+        self, token: Token, max_prefix_arity: int, context: List[Union[Token, Operator]]
+    ) -> Iterable[Operator]:
         if token.token in self.operator_table:
-            return super().resolve(token, max_prefix_arity)
+            return super().resolve(token, max_prefix_arity, context)
 
         symbol = token.token
 
-        # Apply R-like transformations to operator
-        symbol = re.sub(
-            r"[+\-]*\-[+\-]*", "-", symbol
-        )  # Any sequence of '+' and '-' -> '-'
-        symbol = re.sub(r"[+]{2,}", "+", symbol)  # multiple sequential '+' -> '+'
+        # Keep track the number of "+" and "-" characters; if an odd number "-"
+        # than "-", else "+"
+        while True:
+            m = re.search(r"[+\-]{2,}", symbol)
+            if not m:
+                break
+            symbol = (
+                symbol[: m.start(0)] + "-"
+                if len(m.group(0).replace("+", "")) % 2
+                else "+" + symbol[m.end(0) :]
+            )
 
         if symbol in self.operator_table:
-            return [self._resolve(token, symbol, max_prefix_arity)]
+            return [self._resolve(token, symbol, max_prefix_arity, context)]
 
         return [
-            self._resolve(token, sym, max_prefix_arity if i == 0 else 0)
+            self._resolve(token, sym, max_prefix_arity if i == 0 else 0, context)
             for i, sym in enumerate(symbol)
         ]
diff --git a/formulaic/parser/types/__init__.py b/formulaic/parser/types/__init__.py
index 1d8d97d..be08d8e 100644
--- a/formulaic/parser/types/__init__.py
+++ b/formulaic/parser/types/__init__.py
@@ -2,8 +2,17 @@ from .ast_node import ASTNode
 from .factor import Factor
 from .operator import Operator
 from .operator_resolver import OperatorResolver
+from .structured import Structured
 from .term import Term
 from .token import Token
 
 
-__all__ = ["ASTNode", "Factor", "Operator", "OperatorResolver", "Term", "Token"]
+__all__ = [
+    "ASTNode",
+    "Factor",
+    "Operator",
+    "OperatorResolver",
+    "Structured",
+    "Term",
+    "Token",
+]
diff --git a/formulaic/parser/types/operator.py b/formulaic/parser/types/operator.py
index 0225376..fd16981 100644
--- a/formulaic/parser/types/operator.py
+++ b/formulaic/parser/types/operator.py
@@ -1,8 +1,11 @@
+from __future__ import annotations
+
 from enum import Enum
 from numbers import Number
-from typing import Callable, Iterable, Union
+from typing import Callable, List, Iterable, Union
 
 from .term import Term
+from .token import Token
 
 
 class Operator:
@@ -27,6 +30,9 @@ class Operator:
             if 'postfix', the operator comes after its arguments.
         to_terms: A callable that maps the arguments pass to the operator to
             an iterable of `Term` instances.
+        accepts_context: A callable that will receive a list of Operator and
+            Token instances that describe the context in which the operator
+            would be applied if this callable returns `True`.
     """
 
     class Associativity(Enum):
@@ -48,6 +54,7 @@ class Operator:
         associativity: Union[str, Associativity] = "none",
         fixity: Union[str, Fixity] = "infix",
         to_terms: Callable[..., Iterable[Term]] = None,
+        accepts_context: Callable[[List[Union[Token, Operator]]], bool] = None,
     ):
         self.symbol = symbol
         self.arity = arity
@@ -55,6 +62,7 @@ class Operator:
         self.associativity = associativity
         self.fixity = fixity
         self._to_terms = to_terms
+        self._accepts_context = accepts_context
 
     @property
     def associativity(self):
@@ -77,5 +85,19 @@ class Operator:
             raise RuntimeError(f"`to_terms` is not implemented for '{self.symbol}'.")
         return self._to_terms(*args)
 
+    def accepts_context(self, context: List[Union[Token, Operator]]):
+        if self._accepts_context:
+            # We only need to pass on tokens and operators with precedence less
+            # than or equal to ourselves, since all other operators will be
+            # evaluated before us.
+            return self._accepts_context(
+                [
+                    c
+                    for c in context
+                    if isinstance(c, Token) or c.precedence <= self.precedence
+                ]
+            )
+        return True
+
     def __repr__(self):
         return self.symbol
diff --git a/formulaic/parser/types/operator_resolver.py b/formulaic/parser/types/operator_resolver.py
index 7244900..2024e42 100644
--- a/formulaic/parser/types/operator_resolver.py
+++ b/formulaic/parser/types/operator_resolver.py
@@ -1,6 +1,6 @@
 import abc
 from collections import defaultdict
-from typing import List
+from typing import List, Union
 
 from ..utils import exc_for_token
 from .operator import Operator
@@ -42,7 +42,9 @@ class OperatorResolver(metaclass=abc.ABCMeta):
         """
         ...  # pragma: no cover
 
-    def resolve(self, token: Token, max_prefix_arity: int) -> List[Operator]:
+    def resolve(
+        self, token: Token, max_prefix_arity: int, context: List[Union[Token, Operator]]
+    ) -> List[Operator]:
         """
         Return a list of operators to apply for a given token in the AST
         generation.
@@ -52,10 +54,19 @@ class OperatorResolver(metaclass=abc.ABCMeta):
                 be resolved.
             max_prefix_arity: The number operator unclaimed tokens preceding the
                 operator in the formula string.
+            context: The current list of operators into which the operator to be
+                resolved will be placed. This will be a list of `Operator`
+                instances or tokens (tokens are return for grouping operators).
         """
-        return [self._resolve(token, token.token, max_prefix_arity)]
+        return [self._resolve(token, token.token, max_prefix_arity, context)]
 
-    def _resolve(self, token: Token, symbol: str, max_prefix_arity: int) -> Operator:
+    def _resolve(
+        self,
+        token: Token,
+        symbol: str,
+        max_prefix_arity: int,
+        context: List[Union[Token, Operator]],
+    ) -> Operator:
         """
         The default operator resolving logic.
         """
@@ -64,10 +75,13 @@ class OperatorResolver(metaclass=abc.ABCMeta):
         candidates = [
             candidate
             for candidate in self.operator_table[symbol]
-            if max_prefix_arity == 0
-            and candidate.fixity is Operator.Fixity.PREFIX
-            or max_prefix_arity > 0
-            and candidate.fixity is not Operator.Fixity.PREFIX
+            if (
+                max_prefix_arity == 0
+                and candidate.fixity is Operator.Fixity.PREFIX
+                or max_prefix_arity > 0
+                and candidate.fixity is not Operator.Fixity.PREFIX
+            )
+            and candidate.accepts_context(context)
         ]
         if not candidates:
             raise exc_for_token(token, f"Operator `{symbol}` is incorrectly used.")
diff --git a/formulaic/parser/types/structured.py b/formulaic/parser/types/structured.py
new file mode 100644
index 0000000..57d32fc
--- /dev/null
+++ b/formulaic/parser/types/structured.py
@@ -0,0 +1,183 @@
+from __future__ import annotations
+
+from typing import (
+    Any,
+    Callable,
+    Dict,
+    Generator,
+    Generic,
+    Optional,
+    Iterable,
+    TypeVar,
+    Union,
+)
+
+
+ItemType = TypeVar("ItemType")
+_MISSING = object()
+
+
+class Structured(Generic[ItemType]):
+    """
+    Represents an arbitrarily nested structure.
+
+    In the context of Formulaic, `Structured` is used to represent the structure
+    of parsed formulae (e.g. separating the `lhs` and `rhs` of a formula), and
+    the derived the model specs/matrices. Using `Structured` allows Formulaic to
+    abstract out concerns regarding structured data, and focus on getting the
+    containers right for singular use-cases.
+
+    `Structured` instances have an (optional) root node (accessible via `.root`
+    and `[None]`), and arbitrarily many named structured children. A
+    `Structured` instance can be constructed using:
+    ```
+    Structured(<root object (optional)>, key=<child object>, key2=<another object>)
+    ```
+    The structure is then explored using <instance>.<key>.
+
+    `Structured` instances also have a notion of "mapped attributes", which when
+    accessed will recursively request the attribute on all children and return
+    a new `Structured` instance with the result (and the same structure).
+
+    Notes:
+        - tuples have special meaning in `Structured` instances, and are treated
+            part of the structure, and so when maps are applied, they tuples are
+            iterated over (vs. lists which would be acted upon as a whole).
+    """
+
+    def __init__(
+        self,
+        root: Any = _MISSING,
+        *,
+        _mapped_attrs: Iterable[str] = None,
+        _metadata: Dict[str, Any] = None,
+        **structure,
+    ):
+        if any(key.startswith("_") for key in structure):
+            raise ValueError(
+                "Substructure keys cannot start with an underscore. "
+                f"The invalid keys are: {set(key for key in structure if key.startswith('_'))}."
+            )
+        if root is not _MISSING:
+            self.root = root
+        self._structure = structure
+        self._mapped_attrs = set(_mapped_attrs or ())
+        self._metadata = _metadata
+
+    @property
+    def _has_root(self) -> bool:
+        "Whether this instance of `Structured` has a root node."
+        return hasattr(self, "root")
+
+    def _map(
+        self, func: Callable[[ItemType], Any], recurse: bool = True
+    ) -> Structured[Any]:
+        """
+        Map a callable object onto all the structured objects, returning a
+        `Structured` instance with identical structure, where the original
+        objects are replaced with the output of `func`.
+
+        Args:
+            func: The callable to apply to all objects contained in the
+                `Structured` instance.
+            recurse: Whether to recursively map, or only map one level deep (the
+                objects directly referenced by this `StructuredInstance`).
+                When `True`, if objects within this structure are `Structured`
+                instances also, then the map will be applied only on the leaf
+                nodes (otherwise `func` will received `Structured` instances).
+                (default: True).
+
+        Returns:
+            A `Structured` instance with the same structure as this instance,
+            but with all objects transformed under `func`.
+        """
+
+        def apply_func(obj):
+            if recurse and isinstance(obj, Structured):
+                return obj._map(func, recurse=True)
+            if isinstance(obj, tuple):
+                return tuple(func(o) for o in obj)
+            return func(obj)
+
+        return Structured[ItemType](
+            apply_func(self.root) if self._has_root else _MISSING,
+            **{key: apply_func(obj) for key, obj in self._structure.items()},
+        )
+
+    def _to_dict(self, recurse: bool = True) -> Dict[Optional[str], Any]:
+        """
+        Generate a dictionary representation of this structure.
+
+        The root node, if present, will have a `None` key assigned to it.
+
+        Args:
+            recurse: Whether to recursively convert any nested `Structured`
+                instances into dictionaries also. If `False`, any nested
+                `Structured` instances will be surfaced in the generated
+                dictionary.
+
+        Returns:
+            The dictionary representation of this `Structured` instance.
+        """
+
+        def do_recursion(obj):
+            if recurse and isinstance(obj, Structured):
+                return obj._to_dict()
+            return obj
+
+        items = {}
+        if self._has_root:
+            items[None] = do_recursion(self.root)
+        items.update(
+            {key: do_recursion(value) for key, value in self._structure.items()}
+        )
+        return items
+
+    def __dir__(self):
+        return super().__dir__() + list(self._structure)
+
+    def __getattr__(self, attr):
+        if not attr.startswith("_") and attr in self._structure:
+            return self._structure[attr]
+        if attr in self._mapped_attrs:
+            return self._map(lambda x: getattr(x, attr))
+        raise AttributeError(attr)
+
+    def __getitem__(self, key):
+        if key is None and self._has_root:
+            return self.root
+        if isinstance(key, str) and not key.startswith("_") and key in self._structure:
+            return self._structure[key]
+        if isinstance(key, int) and self._has_root and isinstance(self.root, tuple):
+            return self.root[key]
+        raise KeyError(key)
+
+    def __iter__(self) -> Generator[Union[ItemType, Structured[ItemType]]]:
+        if self._has_root:
+            yield self.root
+        yield from self._structure.values()
+
+    def __len__(self) -> int:
+        return int(self._has_root) + len(self._structure)
+
+    def __str__(self):
+        return self.__repr__(to_str=str)
+
+    def __repr__(self, to_str=repr):
+        import textwrap
+
+        out = []
+        for key, value in self._to_dict(recurse=False).items():
+            if not key:
+                out.append("root:")
+            else:
+                out.append(f".{key}")
+            if isinstance(value, tuple):
+                for i, obj in enumerate(value):
+                    out.append(f"    [{i}]:")
+                    out.append(textwrap.indent(to_str(obj), "        "))
+            else:
+                out.append(textwrap.indent(to_str(value), "    "))
+            if not key and len(self) > 1:
+                out.append("\nStructured children:\n")
+        return "\n".join(out)


Test patch:
diff --git a/tests/parser/test_parser.py b/tests/parser/test_parser.py
index e727a23..17c4010 100644
--- a/tests/parser/test_parser.py
+++ b/tests/parser/test_parser.py
@@ -1,8 +1,10 @@
+import re
+
 import pytest
 
 from formulaic.errors import FormulaParsingError, FormulaSyntaxError
 from formulaic.parser import FormulaParser, DefaultOperatorResolver
-from formulaic.parser.types import Token
+from formulaic.parser.types import Structured, Token
 
 
 FORMULA_TO_TOKENS = {
@@ -32,15 +34,23 @@ FORMULA_TO_TERMS = {
     "(a - a) + b": ["1", "b"],
     "a + (b - a)": ["1", "a", "b"],
     # Check that "0" -> "-1" substitution works as expected
+    "+0": [],
+    "(+0)": ["1"],
     "0 + 0": [],
     "0 + 0 + 1": ["1"],
     "0 + 0 + 1 + 0": [],
     "0 - 0": ["1"],
-    "0 ~ 0": ([],),
+    "0 ~ 0": {"lhs": [], "rhs": []},
     # Formula separators
-    "~ a + b": (["1", "a", "b"],),
-    "a ~ b + c": (["a"], ["1", "b", "c"]),
-    "a ~ b ~ c": (["a"], ["b"], ["1", "c"]),
+    "~ a + b": ["1", "a", "b"],
+    "a ~ b + c": {"lhs": ["a"], "rhs": ["1", "b", "c"]},
+    # Formula parts
+    "a | b": (["1", "a"], ["1", "b"]),
+    "a | b | c": (["1", "a"], ["1", "b"], ["1", "c"]),
+    "a | b ~ c | d": {
+        "lhs": (["a"], ["b"]),
+        "rhs": (["1", "c"], ["1", "d"]),
+    },
     # Products
     "a:b": ["1", "a:b"],
     "a * b": ["1", "a", "a:b", "b"],
@@ -55,20 +65,24 @@ FORMULA_TO_TERMS = {
     # Unary operations
     "+1": ["1"],
     "-0": ["1"],
+    "+x": ["1", "x"],
+    "-x": ["1"],
 }
 
 PARSER = FormulaParser()
 
 
 class TestFormulaParser:
-    @pytest.mark.parametrize("formula,tokens", FORMULA_TO_TOKENS.items())
-    def test_get_tokens(self, formula, tokens):
-        assert PARSER.get_tokens(formula) == tokens
+    # @pytest.mark.parametrize("formula,tokens", FORMULA_TO_TOKENS.items())
+    # def test_get_tokens(self, formula, tokens):
+    #     assert PARSER.get_tokens(formula) == tokens
 
     @pytest.mark.parametrize("formula,terms", FORMULA_TO_TERMS.items())
     def test_to_terms(self, formula, terms):
         generated_terms = PARSER.get_terms(formula)
-        if isinstance(generated_terms, tuple):
+        if isinstance(generated_terms, Structured):
+            comp = generated_terms._map(sorted)._to_dict()
+        elif isinstance(generated_terms, tuple):
             comp = tuple(
                 sorted([str(term) for term in group]) for group in generated_terms
             )
@@ -76,9 +90,13 @@ class TestFormulaParser:
             comp = sorted([str(term) for term in generated_terms])
         assert comp == terms
 
-    def test_invalid_unary_negation(self):
+    def test_invalid_formula_separation(selF):
+        with pytest.raises(FormulaParsingError):
+            PARSER.get_terms("a ~ b ~ c")
+
+    def test_invalid_part_separation(selF):
         with pytest.raises(FormulaParsingError):
-            PARSER.get_terms("(-10)")
+            PARSER.get_terms("(a | b)")
 
     def test_invalid_use_of_zero(self):
         with pytest.raises(FormulaParsingError):
@@ -104,21 +122,30 @@ class TestDefaultOperatorResolver:
         return DefaultOperatorResolver()
 
     def test_resolve(self, resolver):
-        assert len(resolver.resolve(Token("+++++"), 1)) == 1
-        assert resolver.resolve(Token("+++++"), 1)[0].symbol == "+"
-        assert resolver.resolve(Token("+++++"), 1)[0].arity == 2
 
-        assert len(resolver.resolve(Token("+++-+"), 1)) == 1
-        assert resolver.resolve(Token("+++-+"), 1)[0].symbol == "-"
-        assert resolver.resolve(Token("+++-+"), 1)[0].arity == 2
+        assert len(resolver.resolve(Token("+++++"), 1, [])) == 1
+        assert resolver.resolve(Token("+++++"), 1, [])[0].symbol == "+"
+        assert resolver.resolve(Token("+++++"), 1, [])[0].arity == 2
+
+        assert len(resolver.resolve(Token("+++-+"), 1, [])) == 1
+        assert resolver.resolve(Token("+++-+"), 1, [])[0].symbol == "-"
+        assert resolver.resolve(Token("+++-+"), 1, [])[0].arity == 2
 
-        assert len(resolver.resolve(Token("*+++-+"), 1)) == 2
-        assert resolver.resolve(Token("*+++-+"), 1)[0].symbol == "*"
-        assert resolver.resolve(Token("*+++-+"), 1)[0].arity == 2
-        assert resolver.resolve(Token("*+++-+"), 1)[1].symbol == "-"
-        assert resolver.resolve(Token("*+++-+"), 1)[1].arity == 1
+        assert len(resolver.resolve(Token("*+++-+"), 1, [])) == 2
+        assert resolver.resolve(Token("*+++-+"), 1, [])[0].symbol == "*"
+        assert resolver.resolve(Token("*+++-+"), 1, [])[0].arity == 2
+        assert resolver.resolve(Token("*+++-+"), 1, [])[1].symbol == "-"
+        assert resolver.resolve(Token("*+++-+"), 1, [])[1].arity == 1
 
         with pytest.raises(
             FormulaSyntaxError, match="Operator `/` is incorrectly used."
         ):
-            resolver.resolve(Token("*/"), 2)
+            resolver.resolve(Token("*/"), 2, [])
+
+    def test_accepts_context(self, resolver):
+        tilde_operator = resolver.resolve(Token("~"), 1, [])[0]
+
+        with pytest.raises(
+            FormulaSyntaxError, match=re.escape("Operator `~` is incorrectly used.")
+        ):
+            resolver.resolve(Token("~"), 1, [tilde_operator])
diff --git a/tests/parser/types/test_operator_resolver.py b/tests/parser/types/test_operator_resolver.py
index 223f6d2..0335f67 100644
--- a/tests/parser/types/test_operator_resolver.py
+++ b/tests/parser/types/test_operator_resolver.py
@@ -26,17 +26,17 @@ class TestOperatorResolver:
         return DummyOperatorResolver()
 
     def test_resolve(self, resolver):
-        assert resolver.resolve(Token("+"), 1)[0] is OPERATOR_PLUS
-        assert resolver.resolve(Token("-"), 0)[0] is OPERATOR_UNARY_MINUS
+        assert resolver.resolve(Token("+"), 1, [])[0] is OPERATOR_PLUS
+        assert resolver.resolve(Token("-"), 0, [])[0] is OPERATOR_UNARY_MINUS
 
         with pytest.raises(FormulaSyntaxError):
-            resolver.resolve(Token("@"), 0)
+            resolver.resolve(Token("@"), 0, [])
 
         with pytest.raises(FormulaSyntaxError):
-            resolver.resolve(Token("+"), 0)
+            resolver.resolve(Token("+"), 0, [])
 
         with pytest.raises(FormulaSyntaxError):
-            resolver.resolve(Token("-"), 1)
+            resolver.resolve(Token("-"), 1, [])
 
         with pytest.raises(FormulaParsingError, match="Ambiguous operator `:`"):
-            resolver.resolve(Token(":"), 1)
+            resolver.resolve(Token(":"), 1, [])
diff --git a/tests/parser/types/test_structured.py b/tests/parser/types/test_structured.py
new file mode 100644
index 0000000..7cbee94
--- /dev/null
+++ b/tests/parser/types/test_structured.py
@@ -0,0 +1,103 @@
+import re
+
+import pytest
+
+from formulaic.parser.types import Structured
+
+
+class TestStructured:
+    def test_constructor(self):
+        obj = object()
+        assert Structured(obj).root is obj
+        assert Structured(obj, key="asd").root is obj
+        assert len(Structured(obj, key="asd")) == 2
+        assert {
+            d for d in dir(Structured(obj, key="asd")) if not d.startswith("_")
+        } == {"root", "key"}
+
+        with pytest.raises(
+            ValueError,
+            match=re.escape(
+                "Substructure keys cannot start with an underscore. The invalid keys are: {'_invalid'}."
+            ),
+        ):
+            Structured(_invalid=True)
+
+    def test_access_structure(self):
+
+        s = Structured("Hello", key="asd")
+        assert s.root == "Hello"
+        assert s[None] == "Hello"
+        assert s.key == "asd"
+        assert s["key"] == "asd"
+
+        assert Structured(("1", "2"))[0] == "1"
+
+        assert Structured()._has_root is False
+        with pytest.raises(AttributeError, match="root"):
+            Structured().root
+        with pytest.raises(KeyError, match="root"):
+            Structured()["root"]
+
+    def test_mapped_attrs(self):
+        class O:
+            def __init__(self, o):
+                self.o = o
+
+        assert (
+            Structured(O(1), a=O(2), b=(O(3), O(4)), _mapped_attrs={"o"}).o._to_dict()
+            == Structured(1, a=2, b=(3, 4))._to_dict()
+        )
+
+    def test__map(self):
+        assert Structured("Hi", a="Hello", b="Greetings")._map(len)._to_dict() == {
+            None: 2,
+            "a": 5,
+            "b": 9,
+        }
+        assert Structured(("Hi", "Dave"), a=["Response", "not", "forthcoming"])._map(
+            len
+        )._to_dict() == {
+            None: (2, 4),
+            "a": 3,
+        }
+        assert Structured(("Hi", Structured("Hi!")), a=Structured(hello="world"))._map(
+            len
+        )._to_dict() == {
+            None: (2, 1),
+            "a": {
+                "hello": 5,
+            },
+        }
+
+    def test_iteration(self):
+        assert list(Structured()) == []
+        assert list(Structured("a")) == ["a"]
+        assert list(Structured(b="b", c="c")) == ["b", "c"]
+        assert list(Structured(c="c", b="b")) == ["c", "b"]
+        assert list(Structured("a", b="b", c="c")) == ["a", "b", "c"]
+        assert list(Structured("a", b="b", c=["c"])) == ["a", "b", ["c"]]
+        assert list(Structured("a", b="b", c=("c",))) == ["a", "b", ("c",)]
+
+    def test_repr(self):
+        assert repr(Structured("a")) == ("root:\n" "    'a'")
+        assert repr(Structured("a", b="b")) == (
+            "root:\n" "    'a'\n\n" "Structured children:\n\n" ".b\n" "    'b'"
+        )
+        assert repr(Structured(("a",), b=("b", "c"))) == (
+            "root:\n"
+            "    [0]:\n"
+            "        'a'\n\n"
+            "Structured children:\n\n"
+            ".b\n"
+            "    [0]:\n"
+            "        'b'\n"
+            "    [1]:\n"
+            "        'c'"
+        )
+
+        # Verify that string representations correctly avoid use of `repr`
+        assert str(Structured("a")) == ("root:\n" "    a")
+        assert str(Structured("a", b="b")) == (
+            "root:\n" "    a\n\n" "Structured children:\n\n" ".b\n" "    b"
+        )
diff --git a/tests/test_formula.py b/tests/test_formula.py
index 17336b7..44446c8 100644
--- a/tests/test_formula.py
+++ b/tests/test_formula.py
@@ -1,9 +1,11 @@
-import pytest
+import re
 
 import pandas
+import pytest
 
 from formulaic import Formula
 from formulaic.errors import FormulaInvalidError, FormulaMaterializerInvalidError
+from formulaic.parser.types import Structured
 
 
 class TestFormula:
@@ -18,7 +20,7 @@ class TestFormula:
 
     @pytest.fixture
     def formula_exprs(self):
-        return Formula("a ~ b ~ c")
+        return Formula("a ~ b")
 
     @pytest.fixture
     def data(self):
@@ -42,7 +44,7 @@ class TestFormula:
         assert Formula.from_spec(f) is f
         assert Formula.from_spec(["a"]) == f
 
-    def test_terms(self, formula_expr, formula_exprs):
+    def test_terms(self, formula_expr):
         assert [str(t) for t in formula_expr.terms] == [
             "1",
             "a",
@@ -53,18 +55,31 @@ class TestFormula:
             "b:c",
             "a:b:c",
         ]
-        assert tuple([str(t) for t in tg] for tg in formula_exprs.terms) == (
-            ["a"],
-            ["b"],
-            ["1", "c"],
-        )
 
     def test_get_model_matrix(self, formula_expr, formula_exprs, data):
         mm_expr = formula_expr.get_model_matrix(data)
         mm_exprs = formula_exprs.get_model_matrix(data, materializer="pandas")
 
         assert mm_expr.shape == (3, 8)
-        assert isinstance(mm_exprs, tuple) and len(mm_exprs) == 3
+        assert isinstance(mm_exprs, Structured) and len(mm_exprs) == 2
+
+    def test_structured(self, formula_exprs):
+        assert formula_exprs.lhs.terms == ["a"]
+        assert formula_exprs.rhs.terms == ["1", "b"]
+        assert Formula("a | b")[0].terms == ["1", "a"]
+
+        with pytest.raises(
+            AttributeError,
+            match=re.escape("This formula has no substructures keyed by 'missing'."),
+        ):
+            formula_exprs.missing
+        with pytest.raises(
+            KeyError,
+            match=re.escape(
+                "This formula does not have any sub-parts indexable via `0`."
+            ),
+        ):
+            formula_exprs[0]
 
     def test_differentiate(self):
         f = Formula("a + b + log(c) - 1")
@@ -74,7 +89,10 @@ class TestFormula:
 
     def test_repr(self, formula_expr, formula_exprs):
         assert repr(formula_expr) == "1 + a + b + c + a:b + a:c + b:c + a:b:c"
-        assert repr(formula_exprs) == "a ~ b ~ 1 + c"
+        assert repr(formula_exprs) == (".lhs\n" "    a\n" ".rhs\n" "    1 + b")
+        assert repr(Formula("a | b")) == (
+            "root:\n" "    [0]:\n" "        1 + a\n" "    [1]:\n" "        1 + b"
+        )
 
     def test_equality(self):
         assert Formula("a + b") == Formula("a+b")
diff --git a/tests/transforms/test_encode_categorical.py b/tests/transforms/test_encode_categorical.py
index a2818a8..36084f8 100644
--- a/tests/transforms/test_encode_categorical.py
+++ b/tests/transforms/test_encode_categorical.py
@@ -160,7 +160,6 @@ def test_encode_categorical_sparse():
 
 
 def _compare_factor_values(a, b, comp=lambda x, y: numpy.allclose(x, y)):
-    print(a, b)
     assert type(a) is type(b)
     if isinstance(a, dict):
         assert sorted(a) == sorted(b)