diff --git a/docs/core_concepts/outputs_and_inspection.md b/docs/core_concepts/outputs_and_inspection.md
index b8b4d80..87c2889 100644
--- a/docs/core_concepts/outputs_and_inspection.md
+++ b/docs/core_concepts/outputs_and_inspection.md
@@ -94,7 +94,9 @@ Two parameters control the physical layout of the stacked output.
 
 `output_format`
 : One of `"dense"`, `"sparse"`, or `"auto"`. `"auto"` picks sparse when it saves memory (for
-example wide one-hot blocks) and dense otherwise. Default `"dense"`.
+example wide one-hot blocks) and dense otherwise. The choice is made once at `fit`, from the
+training output, and stored in `output_format_`, so every `transform` returns the same container,
+even for a single row. Default `"dense"`.
 
 `dtype`
 : The floating-point precision of the output, for example `numpy.float32` to halve memory.
diff --git a/docs/tutorials/adaptive_resolution.md b/docs/tutorials/adaptive_resolution.md
index 9159e01..6f7870a 100644
--- a/docs/tutorials/adaptive_resolution.md
+++ b/docs/tutorials/adaptive_resolution.md
@@ -57,15 +57,16 @@ for name, y in [("simple", simple), ("wiggly", wiggly)]:
 ```
 
 ```text
-simple   -> selected width 15
-wiggly   -> selected width 15
+simple   -> selected width 20
+wiggly   -> selected width 20
 ```
 
 Both widths land inside the `[5, 20]` window without you having to guess a number up front.
-The two happen to match here because the underlying CART selector's split count is governed
-more by its own tree depth and minimum-samples settings than by how wiggly the signal looks;
-with noisier or smaller data, or a narrower window, the two searches can land on different
-widths. The bound is what you control directly, the exact count inside it is data-driven.
+Here both reach the upper bound: with 3000 samples the CART selector finds more informative
+splits than the window admits, so it keeps the most informative ones. Its split count is governed
+more by its own tree depth and minimum-samples settings than by how wiggly the signal looks; with
+smaller data or a wider window the two searches can land on different widths. The bound is what
+you control directly, the exact count inside it is data-driven.
 
 ```{note}
 Fitting a target-aware transformer directly like this, outside a `Pipeline`, normally emits a
diff --git a/pretab/compose/config.py b/pretab/compose/config.py
index b5c7dc7..a9c5ff2 100644
--- a/pretab/compose/config.py
+++ b/pretab/compose/config.py
@@ -175,10 +175,11 @@ def imputation_plan(self, *, is_numerical: bool) -> dict:
         booleans and the imputer ``strategy``. When ``missing_policy`` is ``None``
         the explicit ``*_imputation`` / ``add_missing_indicator`` parameters stay
         authoritative (historical behaviour); otherwise ``missing_policy`` decides.
+        ``add_indicator`` appends a missing indicator built on the raw input next
+        to the representation (see :func:`~pretab.compose.factory.create_transformer`).
         ``add_missing_indicator=True`` with imputation disabled for this column
         kind routes through the standalone ``MissingStateIndicator`` (via
-        ``separate_state``) instead of the imputer's own indicator, since
-        ``SimpleImputer.add_indicator`` only takes effect when the imputer runs.
+        ``separate_state``) instead, preserving its one-column-per-feature output.
         """
         strategy = (
             (self.numerical_imputation or "median")
diff --git a/pretab/compose/factory.py b/pretab/compose/factory.py
index 8ced4d6..9c14479 100644
--- a/pretab/compose/factory.py
+++ b/pretab/compose/factory.py
@@ -10,8 +10,9 @@
 
 import warnings
 
+import numpy as np
 from sklearn.compose import ColumnTransformer
-from sklearn.impute import SimpleImputer
+from sklearn.impute import MissingIndicator, SimpleImputer
 from sklearn.pipeline import FeatureUnion, Pipeline
 from sklearn.preprocessing import MinMaxScaler, StandardScaler
 
@@ -49,6 +50,20 @@
 _KNOT_SPLINE_METHODS = _BMI_SPLINE_METHODS | frozenset({"cubicspline", "naturalspline"})
 
 
+class _PositiveMinMaxScaler(MinMaxScaler):
+    """MinMaxScaler whose output is floored at ``feature_range[0]``.
+
+    Rescales into ``(1e-3, 1)`` ahead of Box-Cox, which requires strictly
+    positive input. Fitted on the training data, a plain MinMaxScaler maps an
+    unseen value below the training minimum to ``<= 0``; flooring only the lower
+    end makes such values transform like the training minimum while every other
+    output (including values above the training maximum) is unchanged.
+    """
+
+    def transform(self, X):
+        return np.maximum(super().transform(X), self.feature_range[0])
+
+
 def _filter_kwargs(allowed, kwargs):
     """Keep only the ``allowed`` keyword arguments that are present in ``kwargs``."""
     return {key: kwargs[key] for key in allowed if key in kwargs}
@@ -154,7 +169,7 @@ def get_numerical_transformer_steps(
     placement = _placement_kwargs(spec, kwargs)
 
     if method == "box-cox":
-        steps.append(("scale_positive", MinMaxScaler(feature_range=(1e-3, 1))))
+        steps.append(("scale_positive", _PositiveMinMaxScaler(feature_range=(1e-3, 1))))
         steps.append(("boxcox", cls(method="box-cox", **filtered)))
     elif method == "yeo-johnson":
         steps.append(("yeojohnson", cls(method="yeo-johnson", **filtered)))
@@ -263,7 +278,6 @@ def create_transformer(method: str, *, is_numerical: bool, config: PreprocessorC
             target_aware=config.target_aware,
             add_imputer=plan["add_imputer"],
             imputer_strategy=plan["strategy"],
-            add_missing_indicator=plan["add_indicator"],
             output_dim=config.output_dim,
             adaptive=config.adaptive,
             min_output_dim=config.min_output_dim if config.adaptive else None,
@@ -293,7 +307,6 @@ def create_transformer(method: str, *, is_numerical: bool, config: PreprocessorC
             method,
             add_imputer=plan["add_imputer"],
             imputer_strategy=plan["strategy"],
-            add_missing_indicator=plan["add_indicator"],
             **constructor_kwargs,
         )
 
@@ -302,9 +315,44 @@ def create_transformer(method: str, *, is_numerical: bool, config: PreprocessorC
         # Emit a dedicated ``__missing`` column (built on the raw input) alongside
         # the imputed representation, so the indicator never enters the basis.
         return FeatureUnion([("representation", pipeline), ("missing", MissingStateIndicator())])
+    if plan["add_indicator"]:
+        # The imputer's missing indicator, built on the raw input next to the
+        # representation instead of being fed through the scaler and basis as a
+        # second feature. It keeps SimpleImputer(add_indicator=True)'s semantics
+        # (a column only for features with missing values at fit) and names.
+        indicator = MissingIndicator(error_on_new=False)
+        return FeatureUnion(
+            [("representation", pipeline), ("missing", indicator)],
+            verbose_feature_names_out=False,
+        )
     return pipeline
 
 
+def _step_names(blocks) -> list[str]:
+    """Return one valid, unique ColumnTransformer step name per ``(kind, feature)`` block.
+
+    A step is named ``f"{kind}_{feature}"`` whenever scikit-learn accepts that
+    name. A label that would put scikit-learn's ``__`` parameter separator into it
+    (one starting with ``_`` or containing ``__``) gets a positional
+    ``f"{kind}_col{i}"`` name instead, made unique against every other step. The
+    public block / feature names never use these internal names: they are derived
+    from the step's kind and its column label (see :mod:`pretab.compose.inspection`).
+    """
+    natural = [f"{kind}_{feature}" for kind, feature in blocks]
+    taken = {name for name in natural if "__" not in name}
+    names = []
+    for position, name in enumerate(natural):
+        if "__" in name:
+            kind = blocks[position][0]
+            name, suffix = f"{kind}_col{position}", 0
+            while name in taken:
+                suffix += 1
+                name = f"{kind}_col{position}_{suffix}"
+            taken.add(name)
+        names.append(name)
+    return names
+
+
 def build_column_transformer(
     config: PreprocessorConfig,
     numerical_features,
@@ -314,19 +362,22 @@ def build_column_transformer(
 ) -> ColumnTransformer:
     """Assemble the per-column pipelines into the final ColumnTransformer.
 
-    Numerical features are prefixed ``num_`` and categorical features ``cat_`` to
-    match the transformer names the Preprocessor exposes; untransformed columns
-    pass through via ``remainder="passthrough"``.
+    Numerical steps are prefixed ``num_`` and categorical steps ``cat_`` to match
+    the block names the Preprocessor exposes (see :func:`_step_names` for labels
+    that cannot be embedded in a step name); untransformed columns pass through via
+    ``remainder="passthrough"``. Each step selects its column by the label's
+    string form: scikit-learn reads an integer selector as a *position*, so the
+    Preprocessor fits and transforms the ColumnTransformer on a frame whose labels
+    are strings (see :func:`~pretab.compose.feature_detection.with_string_labels`).
     """
+    blocks = [("num", feature) for feature in numerical_features]
+    blocks += [("cat", feature) for feature in categorical_features]
     transformers = []
-    for feature in numerical_features:
-        method = config.method_for(feature, is_numerical=True)
-        pipeline = create_transformer(method, is_numerical=True, config=config)
-        transformers.append((f"num_{feature}", pipeline, [feature]))
-    for feature in categorical_features:
-        method = config.method_for(feature, is_numerical=False)
-        pipeline = create_transformer(method, is_numerical=False, config=config)
-        transformers.append((f"cat_{feature}", pipeline, [feature]))
+    for (kind, feature), name in zip(blocks, _step_names(blocks), strict=True):
+        is_numerical = kind == "num"
+        method = config.method_for(feature, is_numerical=is_numerical)
+        pipeline = create_transformer(method, is_numerical=is_numerical, config=config)
+        transformers.append((name, pipeline, [str(feature)]))
     return ColumnTransformer(
         transformers=transformers,
         remainder="passthrough",
diff --git a/pretab/compose/feature_detection.py b/pretab/compose/feature_detection.py
index ebdb1cc..5f1cec3 100644
--- a/pretab/compose/feature_detection.py
+++ b/pretab/compose/feature_detection.py
@@ -5,12 +5,14 @@
 sequence of delegations rather than inlining the classification heuristic.
 """
 
+from collections import Counter
+
 import numpy as np
 import pandas as pd
 
 from ..exceptions import PretabDataError, invalid_param_error
 
-__all__ = ["detect_column_types", "to_dataframe"]
+__all__ = ["bool_columns_as_object", "detect_column_types", "to_dataframe", "with_string_labels"]
 
 
 def to_dataframe(X, *, copy: bool = False) -> pd.DataFrame:
@@ -41,6 +43,58 @@ def to_dataframe(X, *, copy: bool = False) -> pd.DataFrame:
     return X
 
 
+def with_string_labels(X: pd.DataFrame) -> pd.DataFrame:
+    """Return ``X`` with every column label converted to its ``str`` form.
+
+    scikit-learn's :class:`~sklearn.compose.ColumnTransformer` treats an integer
+    column selector as a *position*, so routing a column by an integer label (as
+    in ``pd.DataFrame(array)`` or ``read_csv(header=None)`` after reordering or
+    dropping a column) would select the wrong column. The Preprocessor therefore
+    fits and transforms its ColumnTransformer on string labels. ``X`` itself is
+    returned when every label is already a string; otherwise a shallow copy with
+    relabelled columns, so the caller's frame is never modified.
+
+    Raises
+    ------
+    PretabDataError
+        If two labels share a string form (e.g. ``1`` and ``"1"``), which would
+        make the columns indistinguishable.
+    """
+    if all(isinstance(label, str) for label in X.columns):
+        return X
+    labels = [str(label) for label in X.columns]
+    collisions = sorted(label for label, count in Counter(labels).items() if count > 1)
+    if collisions:
+        raise PretabDataError(
+            f"Column labels must stay unique when converted to strings; {collisions} would collide.\n"
+            "Fix: rename the columns so their string forms are unique."
+        )
+    relabelled = X.copy(deep=False)
+    relabelled.columns = pd.Index(labels)
+    return relabelled
+
+
+def bool_columns_as_object(X: pd.DataFrame) -> pd.DataFrame:
+    """Return ``X`` with its boolean columns cast to ``object``.
+
+    Boolean columns are categorical (see :func:`detect_column_types`), but the
+    categorical pipeline starts with a :class:`~sklearn.impute.SimpleImputer`,
+    which rejects the ``bool`` dtype. As ``object`` columns of ``True`` /
+    ``False`` they are encoded like any other binary categorical column; a
+    missing value of pandas' nullable ``boolean`` dtype becomes ``NaN``, the
+    missing marker the imputer recognizes. ``X`` itself is returned when it has
+    no boolean column, so the caller's frame is never modified.
+    """
+    bool_columns = [label for label, dtype in X.dtypes.items() if dtype.kind == "b"]
+    if not bool_columns:
+        return X
+    cast = X.copy(deep=False)
+    for label in bool_columns:
+        column = X[label].astype(object)
+        cast[label] = column.where(column.notna(), np.nan)
+    return cast
+
+
 def detect_column_types(X, *, cat_cutoff, treat_all_integers_as_numerical, estimator_name="Preprocessor"):
     """Classify each column of ``X`` as numerical or categorical.
 
diff --git a/pretab/compose/inspection.py b/pretab/compose/inspection.py
index 7f3b9b3..a5d8efa 100644
--- a/pretab/compose/inspection.py
+++ b/pretab/compose/inspection.py
@@ -7,8 +7,11 @@
 :func:`build_transformer_summary` renders that metadata as an aligned table.
 """
 
+import copy
+from typing import Any
+
 import numpy as np
-from sklearn.pipeline import FeatureUnion
+from sklearn.pipeline import FeatureUnion, Pipeline
 
 from ..core.logging import get_logger
 from ..core.representation import FeatureLineage
@@ -16,23 +19,48 @@
 logger = get_logger(__name__)
 
 __all__ = [
+    "block_name",
     "build_feature_info",
     "build_feature_lineage",
     "build_transformer_summary",
     "clean_feature_names",
+    "feature_names_out",
     "get_output_slices",
 ]
 
 
+def _block(step_name, columns):
+    """Return ``(kind, feature)`` for a per-column ``num_*`` / ``cat_*`` step, else ``None``."""
+    if step_name == "remainder" or len(columns) != 1:
+        return None
+    kind, sep, _ = str(step_name).partition("_")
+    if not sep or kind not in ("num", "cat"):
+        return None
+    return kind, str(columns[0])
+
+
+def block_name(step_name, columns):
+    """Return the public block name (``"num_<label>"`` / ``"cat_<label>"``) of a step.
+
+    The name is built from the step's kind and its column label rather than read
+    from the step name itself, which differs for labels scikit-learn cannot embed
+    in a step name (see :func:`pretab.compose.factory._step_names`). Steps that are
+    not per-column representation blocks keep their own name.
+    """
+    block = _block(step_name, columns)
+    return step_name if block is None else f"{block[0]}_{block[1]}"
+
+
 def get_output_slices(column_transformer):
     """Return ordered ``(name, start, width)`` spans for each output block.
 
     Reads widths from ``output_indices_`` — the fitted index map that
     ``ColumnTransformer`` already maintains — so no second transform is needed.
+    ``name`` is the public block name (see :func:`block_name`).
     """
     indices = column_transformer.output_indices_
     slices = []
-    for name, transformer, _columns in column_transformer.transformers_:
+    for name, transformer, columns in column_transformer.transformers_:
         if transformer == "drop":
             continue
         span = indices.get(name)
@@ -41,22 +69,78 @@ def get_output_slices(column_transformer):
         width = span.stop - span.start
         if width == 0:
             continue
-        slices.append((name, span.start, width))
+        slices.append((block_name(name, columns), span.start, width))
     return slices
 
 
-def clean_feature_names(column_transformer, names):
+def feature_names_out(column_transformer, input_features=None):
+    """Return the cleaned output feature names of a fitted ColumnTransformer.
+
+    ``input_features`` -- one name per input column, in input order -- replaces the
+    column labels the ColumnTransformer was fitted on (e.g. the synthetic
+    ``feature_i`` labels of an array input), following scikit-learn's contract
+    for an estimator fitted without feature names. Each step's names are rebuilt
+    from its fitted transformer with the new labels, so names that embed a label
+    anywhere (e.g. polynomial terms) are renamed consistently.
+    """
+    fitted_labels = getattr(column_transformer, "feature_names_in_", None)
+    if input_features is None or fitted_labels is None:
+        raw_names = column_transformer.get_feature_names_out(input_features)
+        return clean_feature_names(column_transformer, [str(name) for name in raw_names])
+
+    rename = {str(label): str(name) for label, name in zip(fitted_labels, input_features, strict=True)}
+    raw_names = []
+    for name, transformer, columns in column_transformer.transformers_:
+        if transformer == "drop" or len(columns) == 0:
+            continue
+        labels = [rename[str(column)] for column in columns]
+        if transformer == "passthrough":
+            inner = labels
+        else:
+            inner = _without_feature_names(transformer).get_feature_names_out(labels)
+        raw_names.extend(f"{name}__{inner_name}" for inner_name in inner)
+    return clean_feature_names(column_transformer, raw_names, rename=rename)
+
+
+def _without_feature_names(estimator) -> Any:
+    """Return a shallow copy of a fitted step that accepts any input feature names.
+
+    The ColumnTransformer hands each step a DataFrame column, so scikit-learn
+    estimators inside it record ``feature_names_in_`` and only accept those exact
+    names in ``get_feature_names_out``. Dropping the attribute on a shallow copy
+    (and on copies of any nested pipeline / union steps) makes them fall back to
+    the length check scikit-learn applies to estimators fitted without feature
+    names, while the fitted step itself is left untouched.
+    """
+    if not hasattr(estimator, "get_params"):
+        return estimator
+    stripped = copy.copy(estimator)
+    stripped.__dict__.pop("feature_names_in_", None)
+    if isinstance(stripped, Pipeline):
+        stripped.steps = [(name, _without_feature_names(step)) for name, step in stripped.steps]
+    elif isinstance(stripped, FeatureUnion):
+        stripped.transformer_list = [
+            (name, _without_feature_names(transformer)) for name, transformer in stripped.transformer_list
+        ]
+    return stripped
+
+
+def clean_feature_names(column_transformer, names, *, rename=None):
     """Collapse the per-feature name that sklearn's ColumnTransformer duplicates.
 
-    Each per-column step is named ``f"{kind}_{feature}"`` (see ``compose/factory.py``),
-    and every PreTab transformer's own ``get_feature_names_out`` already bakes the
-    input feature name into each output column, so sklearn's default
-    ``f"{step}__{inner}"`` naming doubles it, e.g. ``"num_age__age_bs0"``. This
-    collapses that back to ``"num_age_bs0"``, leaving passthrough/remainder columns
-    and any name it cannot confidently match unchanged.
+    Each per-column step is a ``num_*`` / ``cat_*`` block for one feature (see
+    ``compose/factory.py``), and every PreTab transformer's own
+    ``get_feature_names_out`` already bakes the input feature name into each output
+    column, so sklearn's default ``f"{step}__{inner}"`` naming doubles it, e.g.
+    ``"num_age__age_bs0"``. This collapses that back to ``"num_age_bs0"``. Other
+    names of a block keep the public block name as their prefix, and
+    passthrough/remainder columns are left unchanged. ``rename`` maps fitted
+    column labels to the names ``names`` were built from (see
+    :func:`feature_names_out`).
     """
-    step_to_feature = {
-        name: columns[0]
+    rename = rename or {}
+    step_to_columns = {
+        name: [rename.get(str(column), column) for column in columns]
         for name, _transformer, columns in column_transformer.transformers_
         if name != "remainder" and len(columns) == 1
     }
@@ -64,15 +148,22 @@ def clean_feature_names(column_transformer, names):
     for raw in names:
         raw = str(raw)
         step_name, sep, inner_name = raw.partition("__")
-        feature = step_to_feature.get(step_name)
-        if not sep or feature is None:
+        columns = step_to_columns.get(step_name)
+        if not sep or columns is None:
             cleaned.append(raw)
             continue
-        if inner_name == feature or inner_name.startswith(f"{feature}_"):
+        block = _block(step_name, columns)
+        if block is not None:
+            kind_prefix, feature = block
+            public_name = f"{kind_prefix}_{feature}"
+        else:
+            feature = str(columns[0])
             kind_prefix = step_name[: -(len(feature) + 1)] if step_name.endswith(f"_{feature}") else ""
+            public_name = step_name
+        if inner_name == feature or inner_name.startswith(f"{feature}_"):
             cleaned.append(f"{kind_prefix}_{inner_name}" if kind_prefix else inner_name)
         else:
-            cleaned.append(raw)
+            cleaned.append(f"{public_name}__{inner_name}")
     return cleaned
 
 
@@ -318,12 +409,12 @@ def build_feature_lineage(column_transformer):
 
         separate_state = _separate_state_branches(transformer)
         if separate_state is not None:
-            representation_pipeline, _missing_indicator = separate_state
-            union_names = [str(value) for value in transformer.get_feature_names_out(list(columns))]
-            representation_width = sum(value.startswith("representation__") for value in union_names)
-            missing_width = sum(value.startswith("missing__") for value in union_names)
-            if representation_width + missing_width != width:
-                representation_width = width - missing_width
+            representation_pipeline, missing_indicator = separate_state
+            # The missing branch comes last in the union; read its width from the
+            # fitted indicator, which may emit no column (MissingIndicator only
+            # marks features that had missing values at fit).
+            missing_width = len(missing_indicator.get_feature_names_out([str(column) for column in columns]))
+            representation_width = width - missing_width
 
             family, component, uses_target, is_interaction = _resolve_block_representation(
                 representation_pipeline, columns
diff --git a/pretab/compose/output.py b/pretab/compose/output.py
index 5a0223c..b670f48 100644
--- a/pretab/compose/output.py
+++ b/pretab/compose/output.py
@@ -22,6 +22,7 @@
     "compute_output_report",
     "format_output",
     "resolve_embedding_dimensions",
+    "resolve_output_format",
     "to_dataframe_output",
     "validate_embedding_request",
 ]
@@ -62,6 +63,26 @@ def validate_embedding_request(embeddings, *, expected: bool, output_kind: str =
         raise IncompatibleParamsError(_EMBEDDINGS_DICT_ONLY)
 
 
+def _density(array) -> float:
+    """Fraction of non-zero entries in a dense or sparse array (0.0 when empty)."""
+    size = int(np.prod(array.shape, dtype=np.int64))
+    nonzero = int(array.count_nonzero()) if sp.issparse(array) else int(np.count_nonzero(array))
+    return float(nonzero) / size if size else 0.0
+
+
+def resolve_output_format(array, output_format, *, threshold=_SPARSE_AUTO_THRESHOLD) -> str:
+    """Resolve ``output_format`` to ``"dense"`` or ``"sparse"`` for ``array``.
+
+    ``"auto"`` picks ``"sparse"`` when the density of ``array`` is below
+    ``threshold``. The Preprocessor calls this once at ``fit`` on the training
+    output, so the container is fixed per fit rather than chosen per batch (a
+    sparse single row would break a dense-only downstream estimator).
+    """
+    if output_format == "auto":
+        return "sparse" if _density(array) < threshold else "dense"
+    return output_format
+
+
 def compute_output_report(array, output_format, *, threshold=_SPARSE_AUTO_THRESHOLD):
     """Resolve the concrete output format and build the memory report.
 
@@ -71,8 +92,9 @@ def compute_output_report(array, output_format, *, threshold=_SPARSE_AUTO_THRESH
         The stacked output. Sparse inputs are inspected through their shape,
         dtype, and stored values without converting them to a dense array.
     output_format : {"auto", "dense", "sparse"}
-        Requested format. ``"auto"`` picks ``"sparse"`` when the density is below
-        ``threshold``.
+        Requested format. ``"auto"`` picks ``"sparse"`` when the density of this
+        ``array`` is below ``threshold``; pass the format resolved at fit (see
+        :func:`resolve_output_format`) to keep the container stable across calls.
     threshold : float, default=0.3
         Density cut-off for the ``"auto"`` decision.
 
@@ -84,16 +106,9 @@ def compute_output_report(array, output_format, *, threshold=_SPARSE_AUTO_THRESH
         ``memory_saved_bytes``.
     """
     size = int(np.prod(array.shape, dtype=np.int64))
-    nonzero = int(array.count_nonzero()) if sp.issparse(array) else int(np.count_nonzero(array))
-    density = float(nonzero) / size if size else 0.0
+    density = _density(array)
     dense_bytes = size * int(array.dtype.itemsize)
-
-    if output_format == "sparse":
-        use_sparse = True
-    elif output_format == "auto":
-        use_sparse = density < threshold
-    else:  # "dense"
-        use_sparse = False
+    use_sparse = resolve_output_format(array, output_format, threshold=threshold) == "sparse"
 
     if use_sparse:
         csr = array.tocsr(copy=False) if sp.issparse(array) else sp.csr_matrix(array)
diff --git a/pretab/compose/serialize.py b/pretab/compose/serialize.py
index c4b213f..025fbef 100644
--- a/pretab/compose/serialize.py
+++ b/pretab/compose/serialize.py
@@ -52,7 +52,18 @@
 )
 
 
+# Builtin types that are valid numpy dtype specifiers, encoded as that dtype.
+_BUILTIN_DTYPES = (bool, int, float, complex)
+
+
 # --- helpers -------------------------------------------------------------
+def _map_elements(data, depth: int, func):
+    """Apply ``func`` to the elements of a ``depth``-level nested list (an array's ``tolist()``)."""
+    if depth == 0:
+        return func(data)
+    return [_map_elements(item, depth - 1, func) for item in data]
+
+
 def _qualname(obj) -> str:
     cls = type(obj)
     return f"{cls.__module__}:{cls.__qualname__}"
@@ -89,10 +100,24 @@ def _encode(obj):
     if obj is UNSET:
         return {"__unset__": True}
     if isinstance(obj, np.ndarray):
-        return {"__ndarray__": {"dtype": obj.dtype.str, "shape": list(obj.shape), "data": obj.tolist()}}
+        data = obj.tolist()
+        if obj.dtype.kind == "O":
+            # tolist() keeps the elements of an object array as they are (e.g.
+            # numpy scalars such as np.bool_), so encode each one.
+            data = _map_elements(data, obj.ndim, _encode)
+        return {"__ndarray__": {"dtype": obj.dtype.str, "shape": list(obj.shape), "data": data}}
     if isinstance(obj, np.dtype):
         return {"__npdtype__": obj.str}
     if isinstance(obj, type):
+        if obj in _BUILTIN_DTYPES:
+            # A dtype-like builtin (e.g. ``Preprocessor(dtype=float)``) is stored as
+            # the numpy dtype it denotes; ``builtins`` is never importable from a spec.
+            return {"__npdtype__": np.dtype(obj).str}
+        if not _module_allowed(obj.__module__):
+            raise PretabSerializationError(
+                f"Cannot serialize the type {obj.__module__}.{obj.__qualname__}: a spec can only reference "
+                f"types from {sorted(_ALLOWED_TOP_LEVEL)}."
+            )
         return {"__type__": f"{obj.__module__}:{obj.__qualname__}"}
     if isinstance(obj, slice):
         return {"__slice__": [obj.start, obj.stop, obj.step]}
@@ -138,6 +163,16 @@ def _json_safe(obj):
 # --- decoding ------------------------------------------------------------
 def _decode_ndarray(payload: dict) -> np.ndarray:
     dtype = np.dtype(payload["dtype"])
+    if dtype.kind == "O":
+        # Fill element by element so a decoded element that is itself a sequence
+        # stays a single element instead of adding a dimension.
+        shape = payload["shape"]
+        elements: list = []
+        _map_elements(payload["data"], len(shape), lambda item: elements.append(_decode(item)))
+        arr = np.empty(len(elements), dtype=object)
+        for index, element in enumerate(elements):
+            arr[index] = element
+        return arr.reshape(shape)
     arr = np.array(payload["data"], dtype=dtype)
     return arr.reshape(payload["shape"])
 
diff --git a/pretab/core/knots.py b/pretab/core/knots.py
index 3a5b7b4..3bf81c2 100644
--- a/pretab/core/knots.py
+++ b/pretab/core/knots.py
@@ -23,6 +23,7 @@
     "quantile_knots",
     "select_knots",
     "spanning_knots",
+    "supplement_interior_knots",
     "uniform_knots",
 ]
 
@@ -142,3 +143,43 @@ def select_knots(knots: np.ndarray, count: int) -> np.ndarray:
         return knots
     idx = np.linspace(0, len(knots) - 1, count).round().astype(int)
     return knots[idx]
+
+
+def supplement_interior_knots(x: np.ndarray, knots: np.ndarray, count: int) -> np.ndarray:
+    """Return unique knots strictly inside the range of ``x``, topped up to ``count``.
+
+    Knots equal to (or outside) ``min(x)`` / ``max(x)`` and repeated knots give a
+    spline basis function with zero-width support -- a dead column -- so they are
+    dropped first. Every remaining knot is kept; when fewer than ``count`` remain,
+    the shortfall is filled with evenly spread quantile and uniform candidates
+    that are themselves strictly interior and not already present. On heavily tied
+    data the quantile candidates coincide with the tied values (often the range
+    boundary), so the uniform candidates guarantee the top-up for any feature with
+    a positive range. Never down-samples: callers trim an overfull set themselves.
+
+    Parameters
+    ----------
+    x : ndarray
+        Values of a single feature (finite).
+    knots : ndarray
+        Candidate interior knots, e.g. quantile knots or selector locations.
+    count : int
+        Minimum number of interior knots to return.
+
+    Returns
+    -------
+    ndarray
+        Sorted, unique knots strictly inside ``(min(x), max(x))``.
+    """
+    x = np.asarray(x, dtype=float)
+    x_min, x_max = x.min(), x.max()
+    knots = np.asarray(knots, dtype=float)
+    knots = np.unique(knots[(knots > x_min) & (knots < x_max)])
+    missing = count - len(knots)
+    if missing <= 0:
+        return knots
+
+    candidates = np.unique(np.concatenate([quantile_knots(x, count), uniform_knots(x, count)]))
+    candidates = candidates[(candidates > x_min) & (candidates < x_max)]
+    candidates = np.setdiff1d(candidates, knots)
+    return np.unique(np.concatenate([knots, select_knots(candidates, missing)]))
diff --git a/pretab/core/logging.py b/pretab/core/logging.py
index 358d46f..cb75f4c 100644
--- a/pretab/core/logging.py
+++ b/pretab/core/logging.py
@@ -34,9 +34,29 @@ def get_logger(name: str = "pretab") -> logging.Logger:
     return logging.getLogger(name)
 
 
+# Marks the handler :func:`configure_logging` attaches, so later calls can tell
+# PreTab's own handler apart from one owned by the host application.
+_OWN_HANDLER_FLAG = "_pretab_own_handler"
+
+
+def _is_host_handler(handler: logging.Handler) -> bool:
+    """Whether ``handler`` belongs to the host (neither a ``NullHandler`` nor PreTab's own)."""
+    return not isinstance(handler, logging.NullHandler) and not getattr(handler, _OWN_HANDLER_FLAG, False)
+
+
 def _has_real_handler(logger: logging.Logger) -> bool:
-    """Whether ``logger`` already carries a handler other than ``NullHandler``."""
-    return any(not isinstance(handler, logging.NullHandler) for handler in logger.handlers)
+    """Whether ``logger`` itself carries a host handler (one PreTab did not attach)."""
+    return any(_is_host_handler(handler) for handler in logger.handlers)
+
+
+def _ancestor_has_real_handler(logger: logging.Logger) -> bool:
+    """Whether a host handler on an ancestor receives ``logger``'s records through propagation."""
+    current = logger
+    while current.propagate and current.parent is not None:
+        current = current.parent
+        if _has_real_handler(current):
+            return True
+    return False
 
 
 def set_verbosity(level: int = 1) -> None:
@@ -57,8 +77,15 @@ def configure_logging(level: int = 1, handler: "logging.Handler | None" = None)
 
     Sets the ``"pretab"`` logger level and attaches a stream handler so PreTab's
     messages become visible. This never touches the root logger, and it is a
-    no-op when a real (non-``NullHandler``) handler is already attached -- so an
-    embedding host (such as DeepTab) that owns handler/level policy always wins.
+    no-op when the host has attached a real (non-``NullHandler``) handler to the
+    ``"pretab"`` logger -- so an embedding host (such as DeepTab) that owns
+    handler/level policy always wins.
+
+    Repeated calls update the level (the handler PreTab attached earlier is
+    reused, never duplicated). When a host handler on an ancestor logger -- such
+    as a root handler from :func:`logging.basicConfig` -- already receives
+    PreTab's records, no handler is attached (and one attached earlier is
+    removed), so every message is printed once.
 
     Parameters
     ----------
@@ -71,7 +98,15 @@ def configure_logging(level: int = 1, handler: "logging.Handler | None" = None)
     if _has_real_handler(_LOGGER):
         return
     set_verbosity(level)
+    own_handlers = [existing for existing in _LOGGER.handlers if getattr(existing, _OWN_HANDLER_FLAG, False)]
+    if _ancestor_has_real_handler(_LOGGER):
+        for existing in own_handlers:
+            _LOGGER.removeHandler(existing)
+        return
+    if own_handlers:
+        return
     if handler is None:
         handler = logging.StreamHandler()
         handler.setFormatter(logging.Formatter("%(name)s: %(message)s"))
+    setattr(handler, _OWN_HANDLER_FLAG, True)
     _LOGGER.addHandler(handler)
diff --git a/pretab/core/selectors.py b/pretab/core/selectors.py
index f4cef30..ffe86b7 100644
--- a/pretab/core/selectors.py
+++ b/pretab/core/selectors.py
@@ -20,16 +20,20 @@
 """
 
 from abc import ABC, abstractmethod
-from typing import Literal
+from typing import Literal, cast
 
 import numpy as np
 from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor
 
 from ..exceptions import IncompatibleParamsError, OptionalDependencyError
-from .knots import quantile_knots, select_knots
+from .knots import quantile_knots, select_knots, uniform_knots
 
 Task = Literal["regression", "classification"]
 
+# Magnitude bound of LightGBM's "split at zero" sentinel threshold (its
+# ``kZeroThreshold`` is 1e-35, reported as the float32 value ~1.0000000180e-35).
+_LIGHTGBM_ZERO_THRESHOLD = 1e-34
+
 
 class BaseLocationSelector(ABC):
     """Abstract base class for count-based, target-aware location selectors.
@@ -97,11 +101,11 @@ def select(
         y_valid = y[valid_mask]
 
         if len(x_valid) < self.min_samples_floor:
-            return quantile_knots(x_valid, min_count)
+            return np.array(self._supplement([], x_valid, min_count))
 
         points, context = self._ordered_candidates(x_valid, y_valid, task)
         if len(points) == 0:
-            return quantile_knots(x_valid, min_count)
+            return np.array(self._supplement([], x_valid, min_count))
 
         locations = self._enforce_spacing(points, x_valid)
 
@@ -116,10 +120,11 @@ def select(
     def _ordered_candidates(self, x_valid: np.ndarray, y_valid: np.ndarray, task: Task) -> tuple[list[float], object]:
         """Fit a model and return candidate locations plus trimming context.
 
-        The candidates must be returned in the selector's preferred order (the
-        order :meth:`_enforce_spacing` should honour): location order for a
-        single tree, gain-descending order for a boosted ensemble. ``context`` is
-        an opaque object passed straight through to :meth:`_trim_over_max`.
+        The candidates must be returned in the selector's preferred order -- the
+        order :meth:`_enforce_spacing` honours, which keeps the earlier of two
+        candidates that are too close: impurity-decrease order for a single tree,
+        gain-descending order for a boosted ensemble. ``context`` is an opaque
+        object passed straight through to :meth:`_trim_over_max`.
         """
         raise NotImplementedError
 
@@ -131,9 +136,10 @@ def _trim_over_max(self, points: list[float], context: object, max_count: int) -
     def _enforce_spacing(self, split_points: list[float], x: np.ndarray) -> list[float]:
         """Drop locations closer than ``min_location_spacing`` of the range.
 
-        Compares each candidate against every already-kept location so the filter
-        is order-independent and honours both ascending (CART) and gain-descending
-        (LightGBM) ordering.
+        Compares each candidate against every already-kept location, so of two
+        candidates that are too close the one earlier in ``split_points`` -- the
+        more important one, given the importance-descending order both selectors
+        produce -- is kept.
         """
         if len(split_points) <= 1:
             return split_points
@@ -149,22 +155,44 @@ def _enforce_spacing(self, split_points: list[float], x: np.ndarray) -> list[flo
         return spaced
 
     def _supplement(self, existing: list[float], x: np.ndarray, target_count: int) -> list[float]:
-        """Top up an under-filled location set with quantile locations.
+        """Top up an under-filled location set to ``target_count`` distinct locations.
 
         Keeps every existing (selector-found) location and fills only the
-        shortfall with quantile candidates, rather than truncating the union
-        (which would preferentially drop the largest existing values).
+        shortfall, rather than truncating the union (which would preferentially
+        drop the largest existing values). Candidates lie strictly inside the
+        range of ``x`` and keep ``min_location_spacing`` from every kept location:
+        evenly spread quantile locations first, then uniform ones. On a tied or
+        discrete feature the quantiles collapse onto a few values (often the range
+        boundary), and the uniform candidates still provide distinct interior
+        locations for any feature with a positive range. Only when the requested
+        count is too dense for the spacing are the remaining uniform locations
+        added without it.
         """
         missing = target_count - len(existing)
         if missing <= 0:
             return existing
 
-        existing_set = set(existing)
-        candidates = [c for c in quantile_knots(x, target_count).tolist() if c not in existing_set]
-        combined = sorted(existing_set | set(candidates[:missing]))
-        if len(combined) > target_count:
-            combined = select_knots(np.array(combined), target_count).tolist()
-        return combined
+        x = np.asarray(x, dtype=float).ravel()
+        x_min, x_max = float(x.min()), float(x.max())
+        if x_max <= x_min:
+            # A zero-range feature has no interior; keep the historical repeated
+            # location so the requested count (and output width) still holds.
+            return sorted([*existing, *quantile_knots(x, missing).tolist()])
+
+        min_distance = self.min_location_spacing * (x_max - x_min)
+        kept = sorted(float(location) for location in existing)
+        for candidates in (quantile_knots(x, target_count), uniform_knots(x, target_count)):
+            eligible = []
+            for candidate in np.unique(candidates[(candidates > x_min) & (candidates < x_max)]):
+                if all(abs(candidate - other) >= min_distance for other in (*kept, *eligible)):
+                    eligible.append(float(candidate))
+            kept = sorted(kept + select_knots(np.array(eligible), missing).tolist())
+            missing = target_count - len(kept)
+            if missing <= 0:
+                return kept
+
+        leftovers = np.setdiff1d(uniform_knots(x, target_count), kept)
+        return sorted(kept + select_knots(leftovers, missing).tolist())
 
 
 class CARTLocationSelector(BaseLocationSelector):
@@ -172,8 +200,9 @@ class CARTLocationSelector(BaseLocationSelector):
 
     A ``DecisionTreeRegressor`` or ``DecisionTreeClassifier`` is fitted to the
     feature against the target, and its split thresholds become the candidate
-    locations. Candidates are spaced out, and if there are too many they are
-    ranked by weighted impurity decrease so the most informative splits are kept.
+    locations. Candidates are ranked by weighted impurity decrease, so the most
+    informative splits are kept both when two candidates are too close to each
+    other and when there are too many.
 
     Parameters
     ----------
@@ -222,62 +251,50 @@ def _ordered_candidates(self, x_valid: np.ndarray, y_valid: np.ndarray, task: Ta
             )
 
         tree.fit(x_valid, y_valid)
-        split_points = self._extract_split_points(tree, x_valid)
-        return split_points, tree
+        importance = self._split_importance(tree, x_valid)
+        # Most informative split first, so that when two candidates are closer
+        # than ``min_location_spacing`` the spacing filter keeps the stronger one.
+        return self._rank(importance, importance), importance
 
     def _trim_over_max(self, points: list[float], context: object, max_count: int) -> list[float]:
-        return self._select_top_locations(points, context, max_count)
-
-    def _extract_split_points(self, tree, x: np.ndarray) -> list[float]:
-        """Collect in-range split thresholds from a fitted decision tree."""
-        tree_structure = tree.tree_
-        split_points = []
-
-        x_min, x_max = float(x.min()), float(x.max())
-
-        for node_id in range(tree_structure.node_count):
-            is_split = tree_structure.children_left[node_id] != tree_structure.children_right[node_id]
-            if not is_split:
-                continue
-            if tree_structure.feature[node_id] != 0:
-                continue
-            threshold = tree_structure.threshold[node_id]
-            if x_min < threshold < x_max:
-                split_points.append(threshold)
+        importance = cast(dict[float, float], context)
+        return sorted(self._rank(points, importance)[:max_count])
 
-        return sorted(set(split_points))
+    @staticmethod
+    def _rank(points, importance: dict[float, float]) -> list[float]:
+        """Order ``points`` by decreasing impurity decrease (ties by location)."""
+        return sorted(points, key=lambda point: (-importance[point], point))
 
-    def _select_top_locations(self, candidates: list[float], tree, max_count: int) -> list[float]:
-        """Keep the locations whose splits reduce impurity the most."""
-        if len(candidates) <= max_count:
-            return candidates
+    def _split_importance(self, tree, x: np.ndarray) -> dict[float, float]:
+        """Map each in-range split threshold to its weighted impurity decrease.
 
+        In a single-feature tree every threshold occurs at most once: after a split
+        at ``t`` no descendant holds samples on both sides of ``t``.
+        """
         tree_structure = tree.tree_
-        split_importance = {}
+        x_min, x_max = float(x.min()), float(x.max())
+        split_importance: dict[float, float] = {}
 
         for node_id in range(tree_structure.node_count):
-            is_split = tree_structure.children_left[node_id] != tree_structure.children_right[node_id]
-            if not is_split:
+            left_child = tree_structure.children_left[node_id]
+            right_child = tree_structure.children_right[node_id]
+            if left_child == right_child or tree_structure.feature[node_id] != 0:
                 continue
-
-            threshold = tree_structure.threshold[node_id]
-            if threshold not in candidates:
+            threshold = float(tree_structure.threshold[node_id])
+            if not x_min < threshold < x_max:
                 continue
 
             n_samples = tree_structure.n_node_samples[node_id]
             impurity = tree_structure.impurity[node_id]
-
-            left_child = tree_structure.children_left[node_id]
-            right_child = tree_structure.children_right[node_id]
             n_left = tree_structure.n_node_samples[left_child]
             n_right = tree_structure.n_node_samples[right_child]
             impurity_left = tree_structure.impurity[left_child]
             impurity_right = tree_structure.impurity[right_child]
+            split_importance[threshold] = float(
+                n_samples * impurity - (n_left * impurity_left + n_right * impurity_right)
+            )
 
-            split_importance[threshold] = n_samples * impurity - (n_left * impurity_left + n_right * impurity_right)
-
-        top = sorted(split_importance, key=lambda k: split_importance[k], reverse=True)[:max_count]
-        return sorted(top)
+        return split_importance
 
 
 class LightGBMLocationSelector(BaseLocationSelector):
@@ -338,9 +355,21 @@ def _import_lightgbm():
     def _ordered_candidates(self, x_valid: np.ndarray, y_valid: np.ndarray, task: Task) -> tuple[list[float], object]:
         lgb = self._import_lightgbm()
 
+        if task == "regression":
+            objective = {"objective": "regression", "metric": "rmse"}
+        else:
+            # LightGBM needs integer class codes: its binary objective treats every
+            # label > 0 as positive and it cannot read string labels.
+            classes, y_valid = np.unique(y_valid, return_inverse=True)
+            if len(classes) < 2:
+                return [], None
+            if len(classes) == 2:
+                objective = {"objective": "binary", "metric": "binary_logloss"}
+            else:
+                objective = {"objective": "multiclass", "metric": "multi_logloss", "num_class": len(classes)}
+
         params = {
-            "objective": "regression" if task == "regression" else "binary",
-            "metric": "rmse" if task == "regression" else "binary_logloss",
+            **objective,
             "num_leaves": 2**self.max_depth,
             "max_depth": self.max_depth,
             "learning_rate": self.learning_rate,
@@ -370,16 +399,39 @@ def _trim_over_max(self, points: list[float], context: object, max_count: int) -
 
     def _extract_split_points_with_gains(self, model, x: np.ndarray) -> dict[float, float]:
         """Collect split thresholds and their cumulative gains from a model."""
+        x = np.asarray(x, dtype=float).ravel()
         x_min, x_max = float(x.min()), float(x.max())
         split_importance: dict[float, float] = {}
 
+        # LightGBM reports a split at zero as the sentinel threshold +/-1e-35
+        # rather than as a bin midpoint. Map it to the midpoint between the
+        # values the split separates, the location CART would report.
+        zero_splits = {
+            1.0: self._midpoint(x[x <= 0], x[x > 0]),
+            -1.0: self._midpoint(x[x < 0], x[x >= 0]),
+        }
+
         model_dict = model.dump_model()
         for tree_info in model_dict["tree_info"]:
-            self._traverse_tree(tree_info["tree_structure"], split_importance, x_min, x_max)
+            self._traverse_tree(tree_info["tree_structure"], split_importance, x_min, x_max, zero_splits)
 
         return split_importance
 
-    def _traverse_tree(self, node: dict, split_importance: dict, x_min: float, x_max: float):
+    @staticmethod
+    def _midpoint(left: np.ndarray, right: np.ndarray) -> float | None:
+        """Midpoint between the largest ``left`` and the smallest ``right`` value."""
+        if left.size == 0 or right.size == 0:
+            return None
+        return float((left.max() + right.min()) / 2.0)
+
+    def _traverse_tree(
+        self,
+        node: dict,
+        split_importance: dict,
+        x_min: float,
+        x_max: float,
+        zero_splits: dict[float, float | None] | None = None,
+    ):
         """Recursively accumulate split gains for the single feature."""
         if "split_feature" not in node:
             return
@@ -387,10 +439,12 @@ def _traverse_tree(self, node: dict, split_importance: dict, x_min: float, x_max
         if node["split_feature"] == 0:
             threshold = node["threshold"]
             gain = node.get("split_gain", 0.0)
-            if x_min < threshold < x_max:
+            if zero_splits is not None and abs(threshold) <= _LIGHTGBM_ZERO_THRESHOLD:
+                threshold = zero_splits[1.0 if threshold >= 0 else -1.0]
+            if threshold is not None and x_min < threshold < x_max:
                 split_importance[threshold] = split_importance.get(threshold, 0.0) + gain
 
         if "left_child" in node:
-            self._traverse_tree(node["left_child"], split_importance, x_min, x_max)
+            self._traverse_tree(node["left_child"], split_importance, x_min, x_max, zero_splits)
         if "right_child" in node:
-            self._traverse_tree(node["right_child"], split_importance, x_min, x_max)
+            self._traverse_tree(node["right_child"], split_importance, x_min, x_max, zero_splits)
diff --git a/pretab/core/supervised.py b/pretab/core/supervised.py
index f6aca3b..a6f0531 100644
--- a/pretab/core/supervised.py
+++ b/pretab/core/supervised.py
@@ -87,6 +87,25 @@ def warn_target_leakage(estimator, y) -> None:
     )
 
 
+def _as_2d(X):
+    """Return ``X`` as 2D input for the wrapped transformer.
+
+    A pandas DataFrame is passed through unchanged, so column names and per-column
+    dtypes reach the wrapped estimator (e.g. a Preprocessor detecting numerical
+    vs categorical columns, or a name-based ColumnTransformer); any other input
+    is converted to a NumPy array, with 1D input reshaped to a single column.
+    """
+    if hasattr(X, "iloc") and getattr(X, "ndim", None) == 2:
+        return X
+    X = np.asarray(X)
+    return X.reshape(-1, 1) if X.ndim == 1 else X
+
+
+def _take_rows(X, indices):
+    """Select rows of a DataFrame or array by position."""
+    return X.iloc[indices] if hasattr(X, "iloc") else X[indices]
+
+
 class CrossFittedTransformer(RepresentationSpecMixin, TransformerMixin, BaseEstimator):
     """Cross-fit a supervised transformer to remove target leakage on training data.
 
@@ -98,7 +117,10 @@ class CrossFittedTransformer(RepresentationSpecMixin, TransformerMixin, BaseEsti
     Parameters
     ----------
     transformer : estimator
-        A supervised (target-aware) PreTab transformer to cross-fit.
+        A supervised (target-aware) PreTab transformer to cross-fit, or an
+        estimator such as a :class:`~pretab.Preprocessor` or ``Pipeline`` built from
+        them. A pandas DataFrame ``X`` is passed to it unchanged (each fold is a row
+        subset), so column names and dtypes are preserved.
     n_folds : int, default=5
         Number of cross-fitting folds. Must be at least 2.
     task : {"regression", "classification"}, default="regression"
@@ -141,9 +163,7 @@ def _fit_full(self, X, y):
             raise InvalidParamError(f"n_folds must be an integer >= 2; got {self.n_folds!r}.")
         if self.task not in ("regression", "classification"):
             raise InvalidParamError(f"task must be 'regression' or 'classification'; got {self.task!r}.")
-        X_arr = np.asarray(X)
-        if X_arr.ndim == 1:
-            X_arr = X_arr.reshape(-1, 1)
+        X_arr = _as_2d(X)
         y_arr = np.asarray(y).ravel()
         if len(X_arr) != len(y_arr):
             raise PretabDataError(f"X and y must have same length. Got {len(X_arr)} and {len(y_arr)}")
@@ -165,10 +185,7 @@ def fit(self, X, y=None):
     def transform(self, X):
         """Transform ``X`` using the transformer fit on all training data."""
         check_is_fitted(self, "estimator_")
-        X_arr = np.asarray(X)
-        if X_arr.ndim == 1:
-            X_arr = X_arr.reshape(-1, 1)
-        return self.estimator_.transform(X_arr)
+        return self.estimator_.transform(_as_2d(X))
 
     def fit_transform(self, X, y=None):
         """Fit and return leakage-free out-of-fold features for the training data."""
@@ -180,8 +197,8 @@ def fit_transform(self, X, y=None):
         try:
             for train_idx, test_idx in splitter.split(X_arr, y_arr):
                 fold = cast(TransformerLike, clone(self.transformer))
-                fold.fit(X_arr[train_idx], y_arr[train_idx])
-                fold_out = np.asarray(fold.transform(X_arr[test_idx]))
+                fold.fit(_take_rows(X_arr, train_idx), y_arr[train_idx])
+                fold_out = np.asarray(fold.transform(_take_rows(X_arr, test_idx)))
                 if fold_out.shape[1] != width:
                     raise IncompatibleParamsError(
                         "Cross-fitting requires a fixed output width across folds; expected "
diff --git a/pretab/encoding/numerical/ple.py b/pretab/encoding/numerical/ple.py
index 81a639c..cca1e86 100644
--- a/pretab/encoding/numerical/ple.py
+++ b/pretab/encoding/numerical/ple.py
@@ -276,7 +276,7 @@ def _apply_piecewise_linear_vectorized(
         if len(thresholds) == 0:
             lower, upper = edges[0], edges[-1]
             width = upper - lower
-            if width > 1e-10:
+            if width > 0:
                 values = np.clip((feature - lower) / width, 0.0, 1.0)
             else:
                 values = np.full(n_samples, 0.5)
@@ -297,7 +297,7 @@ def _apply_piecewise_linear_vectorized(
             upper_edge = edges[bin_idx + 1]
             bin_width = upper_edge - lower_edge
 
-            if bin_width > 1e-10:
+            if bin_width > 0:
                 ple_encoded[mask, bin_idx] = np.clip((values - lower_edge) / bin_width, 0.0, 1.0)
             else:
                 ple_encoded[mask, bin_idx] = 0.5
diff --git a/pretab/expansion/spline/base.py b/pretab/expansion/spline/base.py
index ab365bf..a6cfa61 100644
--- a/pretab/expansion/spline/base.py
+++ b/pretab/expansion/spline/base.py
@@ -23,9 +23,8 @@
 from ...core.knots import (
     basis_to_knots,
     generate_internal_knots,
-    quantile_knots,
     select_knots,
-    uniform_knots,
+    supplement_interior_knots,
 )
 from ...core.parameters import UNSET, validate_placement
 from ...core.policy import RepresentationPolicy, resolve_out_of_range
@@ -195,30 +194,22 @@ def _generate_knots(self, x: np.ndarray, n_knots: int, strategy: str) -> np.ndar
     def _adjust_internal_knots(
         self, x: np.ndarray, internal_knots: np.ndarray, min_knots: int, max_knots: int
     ) -> np.ndarray:
-        """Clip, deduplicate and rebalance internal knots to the allowed count."""
-        internal_knots = np.clip(internal_knots, x.min(), x.max())
-        internal_knots = np.unique(np.sort(internal_knots))
-        if len(internal_knots) < min_knots:
-            internal_knots = self._supplement_knots(x, internal_knots, min_knots)
+        """Restrict, deduplicate and rebalance internal knots to the allowed count.
+
+        Interior knots must lie strictly inside the data range and be unique: a
+        knot on the boundary or a repeated knot collapses a basis function's
+        support to a point, giving a dead column (and, at the boundary, all-zero
+        B/M-spline rows at ``x_max``). Such knots are dropped and the set is
+        topped up to ``min_knots`` before an overfull set is trimmed.
+        """
+        internal_knots = self._supplement_knots(x, internal_knots, min_knots)
         if len(internal_knots) > max_knots:
             internal_knots = select_knots(internal_knots, max_knots)
         return internal_knots
 
     def _supplement_knots(self, x: np.ndarray, internal_knots: np.ndarray, target_count: int) -> np.ndarray:
-        """Add quantile and uniform candidates until ``target_count`` knots exist."""
-        if target_count <= len(internal_knots):
-            return internal_knots
-
-        candidates = [internal_knots]
-        if target_count > 0:
-            candidates.append(quantile_knots(x, target_count))
-            candidates.append(uniform_knots(x, target_count))
-
-        combined = np.unique(np.concatenate(candidates))
-        combined = np.sort(combined)
-        if len(combined) < target_count:
-            combined = uniform_knots(x, target_count)
-        return select_knots(np.asarray(combined), target_count)
+        """Keep the strictly interior, unique knots and top them up to ``target_count``."""
+        return supplement_interior_knots(x, internal_knots, target_count)
 
     def _column_knots(
         self,
@@ -248,16 +239,18 @@ def _column_knots(
                 )
             internal_knots = self._adjust_internal_knots(x_valid, np.asarray(self.knot_locations), min_knots, max_knots)
         elif selector is not None:
-            selected = selector.get_knot_locations(x_valid.reshape(-1, 1), y_valid, task=self.task)
+            # Search exactly the window this feature needs, so the selector's own
+            # importance ranking picks the knots instead of a positional trim.
+            min_knots = min(min_knots, max_knots)
+            selected = selector.get_knot_locations(
+                x_valid.reshape(-1, 1), y_valid, task=self.task, min_knots=min_knots, max_knots=max_knots
+            )
             internal_knots = self._adjust_internal_knots(x_valid, np.asarray(selected), min_knots, max_knots)
         else:
             n_internal = self._basis_to_knots(n_basis)
             internal_knots = self._generate_knots(x_valid, n_internal, strategy)
             internal_knots = self._adjust_internal_knots(x_valid, internal_knots, min_knots, max_knots)
 
-        internal_knots = np.clip(internal_knots, x_min, x_max)
-        internal_knots = np.unique(np.sort(internal_knots))
-
         boundary_left = np.repeat(x_min, self.degree + 1)
         boundary_right = np.repeat(x_max, self.degree + 1)
         return np.concatenate([boundary_left, internal_knots, boundary_right])
diff --git a/pretab/expansion/spline/i_spline.py b/pretab/expansion/spline/i_spline.py
index 85f374c..b95c609 100644
--- a/pretab/expansion/spline/i_spline.py
+++ b/pretab/expansion/spline/i_spline.py
@@ -75,45 +75,45 @@ def __init__(
     def _feature_suffix(self) -> str:
         return "is"
 
-    def _ispline_basis(self, x: np.ndarray, knots: np.ndarray, basis_idx: int) -> np.ndarray:
+    def _ispline_columns(self, x: np.ndarray, knots: np.ndarray) -> np.ndarray:
         """
-        Compute a single I-spline basis function.
-
-        The M-spline is evaluated on a fine grid, integrated with the trapezoidal
-        rule to build a cumulative integral, then linearly interpolated at the
-        requested points and normalized by the full-range integral.
+        Evaluate every I-spline basis function at ``x`` in closed form.
+
+        ``I_k(x)`` is the integral of the normalized M-spline ``M_k`` from the left
+        boundary, i.e. the antiderivative of the B-spline ``B_k`` divided by its
+        full-range integral ``(t_{k+p+1} - t_k) / (p + 1)``. The antiderivative of
+        a B-spline is itself a B-spline of degree ``p + 1``, so the basis is exact
+        for any knot spacing (a fixed quadrature grid cannot resolve basis
+        functions whose support is narrower than the grid step). Values below /
+        above the knot range are 0 / 1. A degenerate basis function with
+        zero-width support (only possible for a knot vector with a repeated
+        knot, e.g. one fitted before boundary knots were de-duplicated)
+        integrates to the step function at that knot, its limit as the support
+        shrinks.
         """
-        x_min_knot = knots[0]
-        x_max_knot = knots[-1]
+        n_coef = len(knots) - self.degree - 1
+        lower, upper = knots[0], knots[-1]
+        antiderivative = BSpline(knots, np.eye(n_coef), self.degree).antiderivative()
+        total = (knots[self.degree + 1 : self.degree + 1 + n_coef] - knots[:n_coef]) / (self.degree + 1)
+
+        x_clipped = np.clip(x, lower, upper)
+        values = antiderivative(x_clipped) - antiderivative(lower)
+        # Every basis function is fully integrated at the right boundary; pin it
+        # instead of evaluating there, where a repeated boundary knot would put
+        # the evaluation on a zero-width interval.
+        values = np.where(x_clipped[:, None] >= upper, total, values)
+        live = total > 0
+        values[:, live] = values[:, live] / total[live]
+        values[:, ~live] = (x_clipped[:, None] >= knots[:n_coef][~live]).astype(float)
+        return np.clip(values, 0.0, 1.0)
+
+    def _ispline_basis(self, x: np.ndarray, knots: np.ndarray, basis_idx: int) -> np.ndarray:
+        """Compute a single I-spline basis function (see :meth:`_ispline_columns`)."""
         n_coef = len(knots) - self.degree - 1
         if basis_idx >= n_coef:
             return np.zeros(len(x))
-
-        grid = np.linspace(x_min_knot, x_max_knot, 200)
-        coef = np.zeros(n_coef)
-        coef[basis_idx] = 1.0
-        spline = BSpline(knots, coef, self.degree, extrapolate=False)
-        m_values = np.nan_to_num(spline(grid), nan=0.0)
-
-        knot_span = knots[basis_idx + self.degree + 1] - knots[basis_idx]
-        if knot_span > 1e-10:
-            m_values = m_values * (self.degree + 1) / knot_span
-
-        cumulative = np.zeros(len(grid))
-        for i in range(1, len(grid)):
-            dx = grid[i] - grid[i - 1]
-            cumulative[i] = cumulative[i - 1] + 0.5 * (m_values[i - 1] + m_values[i]) * dx
-
-        ispline_values = np.interp(x, grid, cumulative, left=0.0, right=cumulative[-1])
-
-        max_integral = cumulative[-1]
-        if max_integral > 1e-10:
-            ispline_values = ispline_values / max_integral
-        return ispline_values
+        return self._ispline_columns(np.asarray(x, dtype=float), knots)[:, basis_idx]
 
     def _design_matrix(self, x: np.ndarray, knots: np.ndarray) -> np.ndarray:
-        n_coef = len(knots) - self.degree - 1
-        design = np.zeros((len(x), n_coef))
-        for i in range(n_coef):
-            design[:, i] = self._ispline_basis(x, knots, i)
+        design = self._ispline_columns(np.asarray(x, dtype=float), knots)
         return np.nan_to_num(design, nan=0.0)
diff --git a/pretab/expansion/spline/m_spline.py b/pretab/expansion/spline/m_spline.py
index f98682e..beac100 100644
--- a/pretab/expansion/spline/m_spline.py
+++ b/pretab/expansion/spline/m_spline.py
@@ -74,7 +74,14 @@ def _feature_suffix(self) -> str:
         return "ms"
 
     def _mspline_basis(self, x: np.ndarray, knots: np.ndarray, basis_idx: int) -> np.ndarray:
-        """Compute a single M-spline basis function."""
+        """Compute a single M-spline basis function.
+
+        ``M_k = (p + 1) / (t_{k+p+1} - t_k) * B_k``, which integrates to one. The
+        values scale like ``1 / span`` and are left unclipped, so the basis does
+        not depend on the units of the feature. Only a basis function whose
+        support is degenerate relative to the knot range (zero width, possible
+        only for a knot vector with a repeated knot) is returned as zeros.
+        """
         n_coef = len(knots) - self.degree - 1
         coef = np.zeros(n_coef)
         coef[basis_idx] = 1.0
@@ -82,17 +89,13 @@ def _mspline_basis(self, x: np.ndarray, knots: np.ndarray, basis_idx: int) -> np
         values = np.nan_to_num(spline(x), nan=0.0)
 
         knot_span = knots[basis_idx + self.degree + 1] - knots[basis_idx]
-        if knot_span > 1e-6:
-            values = values * (self.degree + 1) / knot_span
-            values = np.clip(values, 0.0, 1e6)
-        else:
-            values = np.zeros(len(x))
-        return np.maximum(values, 0.0)
+        if knot_span <= np.finfo(float).eps * (knots[-1] - knots[0]):
+            return np.zeros(len(x))
+        return np.maximum(values * (self.degree + 1) / knot_span, 0.0)
 
     def _design_matrix(self, x: np.ndarray, knots: np.ndarray) -> np.ndarray:
         n_coef = len(knots) - self.degree - 1
         design = np.zeros((len(x), n_coef))
         for i in range(n_coef):
             design[:, i] = self._mspline_basis(x, knots, i)
-        design = np.nan_to_num(design, nan=0.0, posinf=0.0, neginf=0.0)
-        return np.clip(design, -1e3, 1e3)
+        return np.nan_to_num(design, nan=0.0)
diff --git a/pretab/expansion/spline/mixins.py b/pretab/expansion/spline/mixins.py
index b325181..f288697 100644
--- a/pretab/expansion/spline/mixins.py
+++ b/pretab/expansion/spline/mixins.py
@@ -13,7 +13,7 @@
 import numpy as np
 
 from ...core.base import BasePreTabTransformer
-from ...core.knots import generate_internal_knots, select_knots, spanning_knots
+from ...core.knots import generate_internal_knots, select_knots, spanning_knots, supplement_interior_knots
 from ...exceptions import IncompatibleParamsError, PretabDataError
 
 
@@ -84,13 +84,19 @@ def _place_spanning_knots(self, x, y, n_basis, strategy, selector, task, min_int
         window before bracketing.
         """
         x, y = self._finite_column(x, y)
+        x_min, x_max = x.min(), x.max()
         if selector is not None:
             interior = self._place_interior_knots(
                 x, y, n_basis - 2, strategy, selector, task, min_interior, max_interior
             )
-            x_min, x_max = x.min(), x.max()
             return np.concatenate([[x_min], interior, [x_max]])
-        return spanning_knots(x, n_basis, strategy)
+        knots = spanning_knots(x, n_basis, strategy)
+        if len(knots) <= 2:
+            return knots
+        # On tied data quantile knots repeat or land on the endpoints; keep the
+        # endpoints once and make the interior unique and strictly inside.
+        interior = supplement_interior_knots(x, knots[1:-1], n_basis - 2)
+        return np.concatenate([[x_min], interior, [x_max]])
 
     def _place_interior_knots(self, x, y, n_interior, strategy, selector, task, min_interior=None, max_interior=None):
         """Return the interior knots (endpoints excluded) for one feature.
@@ -106,47 +112,43 @@ def _place_interior_knots(self, x, y, n_interior, strategy, selector, task, min_
         On the adaptive selector path ``min_interior`` / ``max_interior`` clamp the
         data-driven count into that window instead. Without a ``selector``,
         ``n_interior`` knots are placed with
-        :func:`pretab.core.knots.generate_internal_knots`.
+        :func:`pretab.core.knots.generate_internal_knots` and made unique and
+        strictly interior (see :func:`pretab.core.knots.supplement_interior_knots`).
         """
         x, y = self._finite_column(x, y)
         if selector is not None:
             if y is None:
                 raise IncompatibleParamsError("A knot selector requires y during fit for target-aware knot placement.")
-            selected = np.asarray(selector.get_knot_locations(x.reshape(-1, 1), y, task=task), dtype=float)
-            x_min, x_max = x.min(), x.max()
-            selected = np.unique(selected[(selected > x_min) & (selected < x_max)])
             if min_interior is None and max_interior is None:
                 # Fixed (non-adaptive) selector path: force exactly ``n_interior``
                 # interior knots so the width stays ``output_dim``.
                 min_interior = max_interior = n_interior
+            # Search exactly this window, so the selector's own importance ranking
+            # picks the knots instead of a positional trim afterwards.
+            selected = selector.get_knot_locations(
+                x.reshape(-1, 1), y, task=task, min_knots=min_interior, max_knots=max_interior
+            )
+            selected = np.asarray(selected, dtype=float)
+            x_min, x_max = x.min(), x.max()
+            selected = np.unique(selected[(selected > x_min) & (selected < x_max)])
             selected = self._clamp_interior_knots(x, selected, min_interior, max_interior, strategy)
             return selected
-        return generate_internal_knots(x, n_interior, strategy)
+        return supplement_interior_knots(x, generate_internal_knots(x, n_interior, strategy), n_interior)
 
     def _clamp_interior_knots(self, x, knots, min_count, max_count, strategy):
         """Clamp a data-driven set of interior knots into ``[min_count, max_count]``.
 
         Down-samples with :func:`pretab.core.knots.select_knots` when there are too
-        many knots and supplements with quantile / uniform interior candidates when
-        there are too few. Endpoints are never added -- the result stays strictly
-        interior.
+        many knots and keeps every knot while topping up with quantile / uniform
+        interior candidates when there are too few. Endpoints are never added --
+        the result stays strictly interior.
         """
         x = np.asarray(x)
         knots = np.unique(np.sort(np.asarray(knots, dtype=float)))
         if max_count is not None and len(knots) > max_count:
             knots = select_knots(knots, max_count)
         if min_count is not None and len(knots) < min_count:
-            x_min, x_max = x.min(), x.max()
-            candidates = [
-                knots,
-                generate_internal_knots(x, min_count, "quantile"),
-                generate_internal_knots(x, min_count, "uniform"),
-            ]
-            combined = np.unique(np.concatenate(candidates))
-            combined = combined[(combined > x_min) & (combined < x_max)]
-            if len(combined) > min_count:
-                combined = select_knots(combined, min_count)
-            knots = combined
+            knots = supplement_interior_knots(x, knots, min_count)
         return knots
 
     def _adaptive_interior_bounds(self, output_dim, selector, *, floor, offset):
diff --git a/pretab/placement/adapters.py b/pretab/placement/adapters.py
index 15a0e65..0f2e816 100644
--- a/pretab/placement/adapters.py
+++ b/pretab/placement/adapters.py
@@ -33,9 +33,9 @@
     "SplinePlacementAdapter",
 ]
 
-# The spline knot selectors have always searched a fixed basis-function window,
-# independent of the requested output_dim (the transformer clamps to output_dim
-# afterwards). These reproduce ``CART/LightGBMKnotSelector``'s defaults.
+# Default basis-function search window of a standalone adapter, reproducing the
+# historical ``CART/LightGBMKnotSelector`` defaults. The spline transformers do
+# not rely on it: they pass their own resolved knot window per call.
 _SPLINE_MIN_BASIS = 3
 _SPLINE_MAX_BASIS = 15
 # Historical default seed used by the spline knot selectors when random_state is
@@ -101,14 +101,24 @@ def get_knot_locations(
         X: np.ndarray,
         y: np.ndarray | None = None,
         task: Task | None = "regression",
+        *,
+        min_knots: int | None = None,
+        max_knots: int | None = None,
     ) -> np.ndarray:
-        """Return sorted internal knot locations for a single feature."""
+        """Return sorted internal knot locations for a single feature.
+
+        ``min_knots`` / ``max_knots`` override the adapter's search window for this
+        call. The spline transformers pass the interior-knot window resolved from
+        their own ``output_dim`` (or adaptive bounds), so that the selector's
+        importance ranking -- not a positional down-sample afterwards -- decides
+        which knots are kept.
+        """
         seed = self.random_state if self.random_state is not None else _SPLINE_DEFAULT_SEED
         strategy = create_placement_strategy(
             target_aware=True,
             placement_strategy=self.placement_strategy,
-            min_count=self.min_knots,
-            max_count=self.max_knots,
+            min_count=self.min_knots if min_knots is None else min_knots,
+            max_count=self.max_knots if max_knots is None else max_knots,
             task=task,
             random_state=seed,
         )
diff --git a/pretab/preprocessor.py b/pretab/preprocessor.py
index a8193b8..5be94d2 100644
--- a/pretab/preprocessor.py
+++ b/pretab/preprocessor.py
@@ -1,5 +1,6 @@
 import hashlib
 import json
+import logging
 import os
 import time
 import warnings
@@ -13,18 +14,25 @@
 
 from .compose.config import PreprocessorConfig
 from .compose.factory import build_column_transformer
-from .compose.feature_detection import detect_column_types, to_dataframe
+from .compose.feature_detection import (
+    bool_columns_as_object,
+    detect_column_types,
+    to_dataframe,
+    with_string_labels,
+)
 from .compose.inspection import (
+    block_name,
     build_feature_info,
     build_feature_lineage,
     build_transformer_summary,
-    clean_feature_names,
+    feature_names_out,
     get_output_slices,
 )
 from .compose.output import (
     compute_output_report,
     format_output,
     resolve_embedding_dimensions,
+    resolve_output_format,
     to_dataframe_output,
     validate_embedding_request,
 )
@@ -93,6 +101,14 @@ def _preset_numerical_method(task) -> str:
 }
 
 
+def _spec_json(spec: dict, **kwargs) -> str:
+    """Serialize a spec to JSON text, raising a typed error for unsupported state."""
+    try:
+        return json.dumps(spec, **kwargs)
+    except (TypeError, ValueError) as exc:
+        raise PretabSerializationError(f"The fitted state cannot be written as JSON: {exc}") from exc
+
+
 def _method_summary(features, *, is_numerical: bool, config: PreprocessorConfig) -> str:
     """Summarize the resolved method(s) actually used across ``features``.
 
@@ -215,10 +231,12 @@ class Preprocessor(TransformerMixin, BaseEstimator):
         Strategy for the ``SimpleImputer`` that runs *before* every categorical method. ``None``
         disables imputation for categorical columns.
     add_missing_indicator : bool, default=False
-        If True, append a binary missing-value indicator column for each imputed feature (via the
-        imputer's ``add_indicator``; a standalone :class:`~pretab.transformers.MissingStateIndicator`
-        is used instead when imputation is disabled for that column kind). Applies to both
-        numerical and categorical pipelines.
+        If True, append a binary (0/1) missing-value indicator column for each imputed feature
+        that has missing values at ``fit`` (scikit-learn's ``MissingIndicator``, computed on the raw
+        input next to the representation, so it never passes through the scaler or basis; a
+        standalone :class:`~pretab.transformers.MissingStateIndicator` is used instead when
+        imputation is disabled for that column kind). Applies to both numerical and categorical
+        pipelines.
     missing_policy : {"error", "propagate", "impute", "impute_with_indicator", "separate_state"} or None, default=None
         High-level missing-value strategy. ``None`` (default) keeps the explicit
         ``numerical_imputation`` / ``categorical_imputation`` / ``add_missing_indicator``
@@ -269,10 +287,12 @@ class Preprocessor(TransformerMixin, BaseEstimator):
         Container used for the transformed output. ``"dense"`` (the default, for
         backward compatibility) returns NumPy arrays; ``"sparse"`` returns SciPy
         CSR matrices (a single stacked CSR when ``return_array=True``, otherwise CSR
-        blocks in the output dict); ``"auto"`` selects ``"sparse"`` when the output
-        density falls below ``0.3`` and ``"dense"`` otherwise. Ignored when
+        blocks in the output dict); ``"auto"`` selects ``"sparse"`` when the density of
+        the training output falls below ``0.3`` and ``"dense"`` otherwise. The choice is
+        made once at ``fit`` (stored in ``output_format_``), so every ``transform``
+        returns the same container regardless of the batch. Ignored when
         :meth:`set_output` requests a pandas or polars DataFrame. Every ``transform``
-        records the resolved choice and its memory footprint in ``output_report_``.
+        records the format and its memory footprint in ``output_report_``.
     dtype : numpy dtype or None, default=None
         Optional dtype to cast the transformed output to (e.g. ``numpy.float32`` to
         halve memory). ``None`` keeps the native ``float64`` output.
@@ -314,6 +334,9 @@ class Preprocessor(TransformerMixin, BaseEstimator):
     output_dims\_ : dict
         Per-feature expanded output-column counts, keyed by input feature name.
         The values sum to ``total_output_dim_``.
+    output_format\_ : {"dense", "sparse"}
+        The output container resolved at ``fit`` from ``output_format`` (``"auto"``
+        is decided from the training output's density).
     output_report\_ : dict
         Memory report for the most recent ``transform``, with keys ``format``
         (``"dense"`` or ``"sparse"``), ``shape``, ``density``, ``dense_bytes``,
@@ -514,6 +537,7 @@ def fit(self, X, y=None, embeddings=None):
             verbose=resolved["verbose"],
         )
 
+        fitted_on_array = isinstance(X, np.ndarray)
         X = to_dataframe(X)
 
         if self.feature_preprocessing:
@@ -580,8 +604,16 @@ def fit(self, X, y=None, embeddings=None):
             categorical_features,
             sparse_threshold=sparse_threshold,
         )
-        self.column_transformer_.fit(X, y)
+        # ColumnTransformer.fit runs fit_transform internally; keep the training
+        # output to resolve output_format="auto" once, from the training density.
+        training_output = self.column_transformer_.fit_transform(self._column_transformer_input(X), y)
+        self.output_format_ = resolve_output_format(training_output, self.output_format)
         self.n_features_in_ = X.shape[1]
+        # scikit-learn convention: feature names are recorded only for named input.
+        if not fitted_on_array and all(isinstance(label, str) for label in X.columns):
+            self.feature_names_in_ = np.asarray(X.columns, dtype=object)
+        elif hasattr(self, "feature_names_in_"):
+            del self.feature_names_in_
 
         self._enforce_output_budget(X.shape[0])
 
@@ -636,6 +668,14 @@ def transform(self, X, embeddings=None, return_array=None):
 
         check_is_fitted(self)
 
+        # Array columns are matched by position, so a different width can never be
+        # routed correctly; an extra column would silently shift or drop features.
+        if isinstance(X, np.ndarray) and X.ndim == 2 and X.shape[1] != self.n_features_in_:
+            raise PretabDataError(
+                f"X has {X.shape[1]} features, but {type(self).__name__} is expecting "
+                f"{self.n_features_in_} features as input."
+            )
+
         X = to_dataframe(X, copy=True)
 
         if self.missing_policy == "error":
@@ -647,13 +687,15 @@ def transform(self, X, embeddings=None, return_array=None):
         output_kind = container if container in ("pandas", "polars") else ("array" if resolved_return_array else "dict")
         validate_embedding_request(embeddings, expected=self.embeddings_, output_kind=output_kind)
 
-        transformed_X = self.column_transformer_.transform(X)
+        transformed_X = self.column_transformer_.transform(self._column_transformer_input(X))
         if not sp.issparse(transformed_X):
             transformed_X = np.asarray(transformed_X)
         if self.dtype is not None:
             transformed_X = transformed_X.astype(self.dtype, copy=False)
 
-        fmt, self.output_report_ = compute_output_report(transformed_X, self.output_format)
+        # Spec files written before output_format_ existed keep the per-call choice.
+        output_format = getattr(self, "output_format_", self.output_format)
+        fmt, self.output_report_ = compute_output_report(transformed_X, output_format)
 
         if container in ("pandas", "polars"):
             return to_dataframe_output(transformed_X, self.get_feature_names_out(), container)
@@ -754,7 +796,13 @@ def get_feature_names_out(self, input_features=None):
         Parameters
         ----------
         input_features : array-like of str or None, default=None
-            Input feature names. Passed through to the underlying column transformer.
+            Input feature names, following scikit-learn's convention. When the
+            preprocessor was fitted on named (string-labelled) columns they must equal
+            ``feature_names_in_``. When it was fitted without feature names (a NumPy
+            array, whose columns are named ``feature_0``, ``feature_1``, ...), any
+            ``n_features_in_`` names are accepted and used in place of the fitted
+            ones, so names propagate through a ``Pipeline`` whose upstream steps
+            output arrays.
 
         Returns
         -------
@@ -763,8 +811,19 @@ def get_feature_names_out(self, input_features=None):
         """
 
         check_is_fitted(self)
-        raw_names = self.column_transformer_.get_feature_names_out(input_features)
-        return np.array(clean_feature_names(self.column_transformer_, raw_names))
+        if input_features is not None:
+            input_features = np.asarray(input_features, dtype=object)
+            if len(input_features) != self.n_features_in_:
+                raise ValueError(
+                    f"input_features should have length equal to number of features "
+                    f"({self.n_features_in_}), got {len(input_features)}"
+                )
+            fitted_names = getattr(self, "feature_names_in_", None)
+            if fitted_names is not None:
+                if not np.array_equal(fitted_names, input_features):
+                    raise ValueError("input_features is not equal to feature_names_in_")
+                input_features = None
+        return np.array(feature_names_out(self.column_transformer_, input_features))
 
     def get_feature_lineage(self):
         """Return per-output-column provenance for the fitted preprocessor.
@@ -819,9 +878,29 @@ def output_dims_(self) -> dict:
                     feature = full.split("__", 1)[-1] if "__" in full else full
                     dims[feature] = dims.get(feature, 0) + 1
             else:
-                dims[columns[0]] = width
+                dims[self._input_label(columns[0])] = width
         return dims
 
+    @staticmethod
+    def _column_transformer_input(X):
+        """Adapt a frame to what the internal ColumnTransformer expects.
+
+        Column labels become strings (scikit-learn reads integer selectors as
+        positions) and boolean columns become ``object`` columns (the categorical
+        imputer rejects ``bool``). The caller's frame is never modified.
+        """
+        return with_string_labels(bool_columns_as_object(X))
+
+    def _input_label(self, column):
+        """Map a ColumnTransformer column back to the label it has in the input.
+
+        The ColumnTransformer is fitted on the string form of every label (see
+        :func:`~pretab.compose.feature_detection.with_string_labels`), so a
+        non-string label such as ``1`` appears there as ``"1"``.
+        """
+        labels = {str(label): label for label in (*self.numerical_features_, *self.categorical_features_)}
+        return labels.get(str(column), column)
+
     def _output_itemsize(self) -> int:
         """Bytes per element of the dense transformed array.
 
@@ -967,9 +1046,12 @@ def get_feature_info(self, verbose=True):
             embeddings=self.embeddings_,
             embedding_dimensions=self.embedding_dimensions_,
         )
+        numerical_feature_info = {self._input_label(key): info for key, info in numerical_feature_info.items()}
+        categorical_feature_info = {self._input_label(key): info for key, info in categorical_feature_info.items()}
 
         if verbose:
-            configure_logging(1)
+            # Render at INFO without lowering a more verbose level set earlier.
+            configure_logging(2 if logger.isEnabledFor(logging.DEBUG) else 1)
             for line in build_transformer_summary(
                 numerical_feature_info,
                 categorical_feature_info,
@@ -981,7 +1063,8 @@ def get_feature_info(self, verbose=True):
 
     def _log_internal_decisions(self):
         """Log fitted internal decisions (bins / knots / centers) at DEBUG."""
-        for name, transformer, _columns in self.column_transformer_.transformers_:
+        for step_name, transformer, columns in self.column_transformer_.transformers_:
+            name = block_name(step_name, columns)
             last_step = transformer.steps[-1][1] if hasattr(transformer, "steps") else transformer
             for attr in (
                 "thresholds_",
@@ -1018,9 +1101,12 @@ def to_spec(self, path=None) -> dict:
         """
         check_is_fitted(self)
         spec = preprocessor_to_spec(self)
+        # Serialize before opening the file, so a failure never truncates an
+        # existing spec at ``path``.
+        text = _spec_json(spec, indent=2)
         if path is not None:
             with open(path, "w", encoding="utf-8") as handle:
-                json.dump(spec, handle, indent=2)
+                handle.write(text)
         return spec
 
     @classmethod
@@ -1080,7 +1166,7 @@ def fingerprint_(self) -> str:
         the same fitted state and configuration produce the same fingerprint.
         """
         check_is_fitted(self)
-        canonical = json.dumps(self._canonical_spec(), sort_keys=True, separators=(",", ":"), ensure_ascii=True)
+        canonical = _spec_json(self._canonical_spec(), sort_keys=True, separators=(",", ":"), ensure_ascii=True)
         return hashlib.sha256(canonical.encode("utf-8")).hexdigest()
 
     def reproducibility_report(self) -> dict:
diff --git a/tests/compose/test_feature_detection.py b/tests/compose/test_feature_detection.py
index 670e8c4..c04986d 100644
--- a/tests/compose/test_feature_detection.py
+++ b/tests/compose/test_feature_detection.py
@@ -4,7 +4,12 @@
 import pandas as pd
 import pytest
 
-from pretab.compose.feature_detection import detect_column_types, to_dataframe
+from pretab.compose.feature_detection import (
+    bool_columns_as_object,
+    detect_column_types,
+    to_dataframe,
+    with_string_labels,
+)
 from pretab.exceptions import InvalidParamError, PretabDataError
 
 
@@ -105,3 +110,28 @@ def test_invalid_cat_cutoff_type_raises():
     df = pd.DataFrame({"x": [1, 2, 3]})
     with pytest.raises(InvalidParamError):
         detect_column_types(df, cat_cutoff="bad", treat_all_integers_as_numerical=False)
+
+
+def test_bool_columns_as_object_casts_only_boolean_columns():
+    frame = pd.DataFrame({"x": [0.5, 1.5], "flag": [True, False], "nullable": pd.array([True, None], dtype="boolean")})
+    cast = bool_columns_as_object(frame)
+    assert cast["x"].dtype == np.float64
+    assert cast["flag"].dtype == object and cast["flag"].tolist() == [True, False]
+    assert cast["nullable"].dtype == object
+    assert cast["nullable"][0] is True and np.isnan(cast["nullable"][1])
+    # The caller's frame keeps its dtypes.
+    assert frame["flag"].dtype == bool and frame["nullable"].dtype == "boolean"
+
+
+def test_bool_columns_as_object_returns_frame_without_bool_columns_unchanged():
+    frame = pd.DataFrame({"x": [0.5, 1.5]})
+    assert bool_columns_as_object(frame) is frame
+
+
+def test_with_string_labels_relabels_without_touching_the_input():
+    frame = pd.DataFrame({1: [1.0], "b": [2.0]})
+    relabelled = with_string_labels(frame)
+    assert list(relabelled.columns) == ["1", "b"]
+    assert list(frame.columns) == [1, "b"]
+    string_frame = pd.DataFrame({"a": [1.0]})
+    assert with_string_labels(string_frame) is string_frame
diff --git a/tests/core/test_cross_fitted.py b/tests/core/test_cross_fitted.py
index 723d971..f39ddda 100644
--- a/tests/core/test_cross_fitted.py
+++ b/tests/core/test_cross_fitted.py
@@ -6,6 +6,7 @@
 """
 
 import warnings
+from typing import cast
 
 import numpy as np
 import pytest
@@ -104,3 +105,67 @@ def test_invalid_task_rejected(data):
     X, y = data
     with pytest.raises(InvalidParamError, match="task"):
         CrossFittedTransformer(PLETransformer(), task="classificaton").fit_transform(X, y)
+
+
+@pytest.mark.parametrize("seed", range(5))
+def test_cross_fitting_a_discrete_feature_keeps_a_fixed_width(seed):
+    """Regression guard for issue #57: tied features gave fold-dependent widths."""
+    from pretab.transformers import RBFExpansionTransformer
+
+    rng = np.random.default_rng(seed)
+    x = rng.integers(0, 6, size=200).astype(float).reshape(-1, 1)
+    y = np.sin(x[:, 0]) + rng.normal(0, 0.5, size=200)
+    out = CrossFittedTransformer(
+        RBFExpansionTransformer(output_dim=10, target_aware=True), n_folds=5, random_state=0
+    ).fit_transform(x, y)
+    assert out.shape == (200, 10)
+
+
+# --- DataFrame input is passed through (issue #58) --------------------------------
+
+
+@pytest.fixture
+def mixed_frame():
+    import pandas as pd
+
+    rng = np.random.default_rng(0)
+    X = pd.DataFrame(
+        {"num": rng.normal(size=300), "city": rng.choice(["a", "b", "c"], size=300)},
+        index=rng.permutation(1000)[:300],
+    )
+    return X, X["num"].to_numpy() ** 2 + rng.normal(scale=0.1, size=300)
+
+
+def test_wrapped_preprocessor_keeps_dataframe_column_types(mixed_frame):
+    """Regression guard for issue #58: np.asarray turned the frame into an object
+    array, so the wrapped Preprocessor treated every column as categorical."""
+    from pretab import Preprocessor
+
+    X, y = mixed_frame
+    direct = Preprocessor(output_dim=6, random_state=0).fit(X, y)
+    cross_fitted = CrossFittedTransformer(Preprocessor(output_dim=6, random_state=0), n_folds=3, random_state=0)
+    out = cross_fitted.fit_transform(X, y)
+
+    wrapped = cast(Preprocessor, cross_fitted.estimator_)
+    assert wrapped.numerical_features_ == ["num"]
+    assert wrapped.categorical_features_ == ["city"]
+    assert out.shape == (300, direct.total_output_dim_)
+    assert list(cross_fitted.get_feature_names_out()) == list(direct.get_feature_names_out())
+    assert len(np.unique(out[:, 0])) > 1  # a supervised encoding, not an "unseen category" code
+    assert cross_fitted.transform(X).shape == out.shape
+
+
+def test_wrapped_column_transformer_selects_by_column_name(mixed_frame):
+    from sklearn.compose import ColumnTransformer
+
+    X, y = mixed_frame
+    wrapped = ColumnTransformer([("ple", PLETransformer(output_dim=4), ["num"])])
+    out = CrossFittedTransformer(wrapped, n_folds=3, random_state=0).fit_transform(X, y)
+    assert out.shape == (300, 4)
+
+
+def test_one_dimensional_input_is_still_a_single_column():
+    rng = np.random.default_rng(0)
+    x = rng.normal(size=200)
+    out = CrossFittedTransformer(PLETransformer(output_dim=4), n_folds=3, random_state=0).fit_transform(x, x**2)
+    assert out.shape == (200, 4)
diff --git a/tests/core/test_feature_map_selector.py b/tests/core/test_feature_map_selector.py
index 2b886c6..20d4ac5 100644
--- a/tests/core/test_feature_map_selector.py
+++ b/tests/core/test_feature_map_selector.py
@@ -100,3 +100,17 @@ def test_lightgbm_selector_places_centers(Cls, data):
     t = Cls(output_dim=5, target_aware=True, task="regression", placement_strategy="lightgbm").fit(X, y)
     assert all(len(c) == 5 for c in t.centers_)
     assert all(np.all(np.diff(c) > 0) for c in t.centers_)
+
+
+@pytest.mark.parametrize("strategy", ["cart", "lightgbm"])
+def test_target_aware_centers_keep_output_dim_on_a_discrete_feature(strategy):
+    """Regression guard for issue #57: tied features returned fewer, duplicate centers."""
+    if strategy == "lightgbm":
+        pytest.importorskip("lightgbm")
+    rng = np.random.default_rng(0)
+    x = rng.integers(0, 5, size=300).astype(float).reshape(-1, 1)
+    y = (x[:, 0] >= 1) + rng.normal(0, 0.1, size=300)
+    transformer = RBFExpansionTransformer(output_dim=10, target_aware=True, placement_strategy=strategy).fit(x, y)
+    centers = transformer.centers_[0]
+    assert transformer.transform(x).shape[1] == 10
+    assert len(np.unique(centers)) == 10
diff --git a/tests/core/test_location_selectors.py b/tests/core/test_location_selectors.py
index e1e6e7f..423969b 100644
--- a/tests/core/test_location_selectors.py
+++ b/tests/core/test_location_selectors.py
@@ -1,3 +1,5 @@
+from typing import cast
+
 import numpy as np
 import pytest
 
@@ -160,3 +162,134 @@ def test_supplement_preserves_high_end_split_end_to_end():
     locations = CARTLocationSelector().select(xs.reshape(-1, 1), ys, task="regression", min_count=6, max_count=6)
 
     assert locations.max() > 5.0, f"expected a high-end location to survive supplementing, got {locations}"
+
+
+def _single_step(seed):
+    rng = np.random.default_rng(seed)
+    x = rng.uniform(0, 10, 1000)
+    step = rng.uniform(2, 8)
+    y = np.where(x > step, 1.0, 0.0) + 0.5 * rng.normal(size=1000)
+    return x.reshape(-1, 1), y
+
+
+def test_cart_spacing_keeps_the_dominant_split_over_a_weak_neighbour():
+    """Regression guard for issue #70: candidates were spaced in ascending location
+    order, so a weak split just below the root split evicted it."""
+    from sklearn.tree import DecisionTreeRegressor
+
+    X, y = _single_step(34)
+    root = DecisionTreeRegressor(max_depth=1).fit(X, y).tree_.threshold[0]
+    locations = CARTLocationSelector().select(X, y, task="regression", min_count=6, max_count=6)
+    assert np.isclose(locations, root).any(), f"root split {root:.3f} missing from {locations}"
+
+
+@pytest.mark.parametrize("seed", range(10))
+def test_cart_keeps_the_root_split_of_a_single_step_target(seed):
+    from sklearn.tree import DecisionTreeRegressor
+
+    X, y = _single_step(seed)
+    root = DecisionTreeRegressor(max_depth=1).fit(X, y).tree_.threshold[0]
+    locations = CARTLocationSelector().select(X, y, task="regression", min_count=6, max_count=6)
+    assert np.isclose(locations, root).any()
+
+
+def test_cart_candidates_are_ordered_by_impurity_decrease(data):
+    X, y = data
+    selector = CARTLocationSelector()
+    candidates, context = selector._ordered_candidates(X, y, "regression")
+    importance = cast("dict[float, float]", context)
+    gains = [importance[c] for c in candidates]
+    assert gains == sorted(gains, reverse=True)
+
+
+# --- tied / discrete features (issue #57) ------------------------------------
+
+
+def _discrete(levels, n=300, seed=0):
+    rng = np.random.default_rng(seed)
+    x = rng.integers(0, levels, size=n).astype(float)
+    return x.reshape(-1, 1), (x >= 2) + rng.normal(0, 0.1, size=n)
+
+
+def _assert_distinct_interior(locations, X, count):
+    assert len(locations) == count
+    assert len(np.unique(locations)) == count
+    assert locations.min() > X.min() and locations.max() < X.max()
+
+
+@pytest.mark.parametrize("selector_cls", [CARTLocationSelector, LightGBMLocationSelector])
+@pytest.mark.parametrize("levels", [2, 3, 5])
+def test_selector_tops_up_tied_features_with_distinct_interior_locations(selector_cls, levels):
+    """Regression guard for issue #57: the quantile top-up collapsed onto the tied
+    values and the range boundary, returning too few / duplicate locations."""
+    if selector_cls is LightGBMLocationSelector:
+        pytest.importorskip("lightgbm")
+    X, y = _discrete(levels)
+    locations = selector_cls().select(X, y, task="regression", min_count=8, max_count=8)
+    _assert_distinct_interior(locations, X, 8)
+
+
+def test_supplement_on_tied_data_respects_spacing_and_count():
+    x = np.repeat([0.0, 1.0, 2.0], 50)
+    selector = CARTLocationSelector()
+    result = selector._supplement([0.5], x, 6)
+    assert len(result) == 6
+    assert 0.5 in result
+    assert min(np.diff(result)) >= selector.min_location_spacing * 2.0
+    assert min(result) > 0.0 and max(result) < 2.0
+
+
+def test_supplement_fills_dense_requests_beyond_the_spacing():
+    x = np.linspace(0.0, 1.0, 1000)
+    result = CARTLocationSelector()._supplement([], x, 150)
+    assert len(result) == 150
+    assert len(np.unique(result)) == 150
+
+
+def test_small_sample_fallback_on_tied_data_is_distinct_and_interior():
+    X = np.array([[0.0], [0.0], [0.0], [1.0], [1.0]])
+    locations = CARTLocationSelector().select(X, np.arange(5.0), min_count=3, max_count=3)
+    _assert_distinct_interior(locations, X, 3)
+
+
+def test_constant_feature_fallback_keeps_the_requested_count():
+    X = np.full((40, 1), 3.0)
+    locations = CARTLocationSelector().select(X, np.arange(40.0), min_count=4, max_count=4)
+    np.testing.assert_array_equal(locations, np.full(4, 3.0))
+
+
+def test_lightgbm_zero_split_sentinel_maps_to_the_bin_midpoint():
+    pytest.importorskip("lightgbm")
+    rng = np.random.default_rng(0)
+    x = rng.integers(0, 5, size=300).astype(float)
+    y = (x >= 1) + rng.normal(0, 0.1, size=300)
+    locations = LightGBMLocationSelector().select(x.reshape(-1, 1), y, min_count=4, max_count=4)
+    # The 0-vs-positive split is reported by LightGBM as threshold 1e-35.
+    assert 0.5 in locations
+    assert (np.abs(locations) > 1e-30).all()
+
+
+@pytest.mark.parametrize(
+    "make_target",
+    [
+        pytest.param(lambda x: np.where(x < 0.9, 1, 2), id="labels-1-2"),
+        pytest.param(lambda x: np.digitize(x, [0.1, 0.9]), id="three-classes"),
+        pytest.param(lambda x: np.where(x < 0.9, "no", "yes"), id="string-labels"),
+    ],
+)
+def test_lightgbm_classification_places_a_location_at_the_class_boundary(make_target):
+    """Regression guard for issue #61: the binary objective on raw labels lost the
+    supervision for non-0/1, multiclass and string targets."""
+    pytest.importorskip("lightgbm")
+    x = np.random.default_rng(0).uniform(0, 1, size=1000)
+    locations = LightGBMLocationSelector().select(
+        x.reshape(-1, 1), make_target(x), task="classification", min_count=2, max_count=2
+    )
+    assert np.abs(locations - 0.9).min() < 0.01, locations
+
+
+def test_lightgbm_classification_with_a_single_class_falls_back():
+    pytest.importorskip("lightgbm")
+    x = np.linspace(0, 1, 200).reshape(-1, 1)
+    locations = LightGBMLocationSelector().select(x, np.ones(200), task="classification", min_count=2, max_count=2)
+    assert len(locations) == 2
diff --git a/tests/core/test_ple_selector.py b/tests/core/test_ple_selector.py
index 498b9ac..497e09c 100644
--- a/tests/core/test_ple_selector.py
+++ b/tests/core/test_ple_selector.py
@@ -95,3 +95,15 @@ def test_lightgbm_selector_places_thresholds(data):
     t = PLETransformer(output_dim=5, task="regression", placement_strategy="lightgbm").fit(X, y)
     assert t.n_bins_per_feature_ == [5]
     assert np.all(np.diff(t.thresholds_[0]) > 0)
+
+
+@pytest.mark.parametrize("levels", [2, 3, 4, 5, 6])
+def test_ple_keeps_output_dim_and_no_constant_columns_on_discrete_features(levels):
+    """Regression guard for issue #57: thresholds on the range boundary produced a
+    constant column and a narrower-than-requested output."""
+    rng = np.random.default_rng(levels)
+    x = rng.integers(0, levels, size=400).astype(float).reshape(-1, 1)
+    y = np.sin(x[:, 0]) + rng.normal(0, 0.3, size=400)
+    out = PLETransformer(output_dim=7).fit(x, y).transform(x)
+    assert out.shape[1] == 7
+    assert (out.std(axis=0) > 0).all()
diff --git a/tests/expansion/spline/test_spline_expansions.py b/tests/expansion/spline/test_spline_expansions.py
index a55d9ca..3ab0b88 100644
--- a/tests/expansion/spline/test_spline_expansions.py
+++ b/tests/expansion/spline/test_spline_expansions.py
@@ -3,8 +3,10 @@
 
 from pretab.transformers import (
     BSplineTransformer,
+    CubicRegressionSplineTransformer,
     ISplineTransformer,
     MSplineTransformer,
+    NaturalCubicSplineTransformer,
 )
 
 
@@ -116,6 +118,27 @@ def test_mspline_handles_nan():
     assert np.isfinite(Xt).all()
 
 
+@pytest.mark.parametrize("scale", [1.0, 1e-2, 1e-3, 1e-7])
+def test_mspline_basis_integrates_to_one_at_any_feature_scale(scale):
+    """Regression guard for issue #63: absolute clipping thresholds flattened the
+    basis on small-scale features and zeroed it on tiny-range ones."""
+    from scipy.integrate import trapezoid
+
+    X = np.random.default_rng(0).uniform(0, 1, (2000, 1)) * scale
+    transformer = MSplineTransformer(output_dim=8, placement_strategy="uniform").fit(X)
+    grid = np.linspace(X.min(), X.max(), 200_001)
+    integrals = trapezoid(transformer.transform(grid.reshape(-1, 1)), grid, axis=0)
+    np.testing.assert_allclose(integrals, 1.0, atol=1e-3)
+
+
+def test_mspline_basis_is_scale_equivariant():
+    X = np.random.default_rng(0).lognormal(0, 2, (2000, 1))
+    unit = MSplineTransformer(output_dim=8).fit(X).transform(X)
+    scaled = MSplineTransformer(output_dim=8).fit(X * 1e-4).transform(X * 1e-4)
+    # M-spline values scale like 1 / span, so rescaling x by c rescales M by 1 / c.
+    np.testing.assert_allclose(scaled * 1e-4, unit, rtol=1e-8)
+
+
 def test_ispline_monotonic_increasing():
     X = np.linspace(0, 10, 200).reshape(-1, 1)
     transformer = ISplineTransformer(output_dim=8, include_bias=False)
@@ -132,6 +155,42 @@ def test_ispline_bounded_unit_interval():
     assert np.all(Xt <= 1.0 + 1e-9)
 
 
+def _exact_ispline(x, knots, degree):
+    """Reference I-spline: the normalized antiderivative of each B-spline basis function."""
+    from scipy.interpolate import BSpline
+
+    n_coef = len(knots) - degree - 1
+    columns = []
+    for i in range(n_coef):
+        antiderivative = BSpline(knots, np.eye(n_coef)[i], degree).antiderivative()
+        lower, upper = antiderivative(knots[0]), antiderivative(knots[-1])
+        columns.append((antiderivative(x) - lower) / (upper - lower))
+    return np.column_stack(columns)
+
+
+@pytest.mark.parametrize("output_dim", [6, 10, 30])
+def test_ispline_matches_exact_integral_on_tight_knots(output_dim):
+    """Regression guard for issue #53: knot spans narrower than a fixed quadrature
+    grid step produced all-zero and misplaced columns."""
+    X = np.random.default_rng(0).lognormal(0, 2, (2000, 1))
+    transformer = ISplineTransformer(output_dim=output_dim).fit(X)
+    Xt = transformer.transform(X)
+
+    exact = _exact_ispline(X[:, 0], transformer.knots_[0], transformer.degree)
+    np.testing.assert_allclose(Xt, exact, atol=1e-10)
+    assert (Xt.max(axis=0) > 0).all()
+    # Every I-spline row is ordered I_0 >= I_1 >= ... >= I_{K-1}.
+    assert (np.diff(Xt, axis=1) <= 1e-12).all()
+
+
+def test_ispline_reaches_zero_and_one_at_the_range_boundaries():
+    X = np.random.default_rng(1).exponential(size=(500, 1))
+    transformer = ISplineTransformer(output_dim=8).fit(X)
+    boundary = transformer.transform(np.array([[X.min()], [X.max()]]))
+    np.testing.assert_allclose(boundary[0], 0.0, atol=1e-12)
+    np.testing.assert_allclose(boundary[1], 1.0, atol=1e-12)
+
+
 def test_ispline_shape_multi_feature():
     rng = np.random.RandomState(2)
     X = rng.uniform(0, 5, size=(100, 2))
@@ -154,3 +213,43 @@ def test_spline_penalty_matrix_symmetric(data):
     P = transformer.get_penalty_matrix()
     assert P.shape[0] == P.shape[1]
     assert np.allclose(P, P.T, atol=1e-9)
+
+
+# --- target-aware splines search their own knot window (issue #69) ------------
+
+
+def _step_data(seed=1):
+    rng = np.random.default_rng(seed)
+    x = rng.uniform(0, 10, 1000)
+    y = np.where(x > 5.0, 1.0, 0.0) + 0.3 * rng.normal(size=1000)
+    return x.reshape(-1, 1), y
+
+
+@pytest.mark.parametrize("cls", [BSplineTransformer, MSplineTransformer, ISplineTransformer])
+@pytest.mark.parametrize("output_dim", [5, 6, 7])
+def test_target_aware_bmi_spline_keeps_the_dominant_split(cls, output_dim):
+    """Regression guard for issue #69: selector knots were trimmed by position, so
+    the split where the target changes was routinely discarded."""
+    X, y = _step_data()
+    transformer = cls(output_dim=output_dim, target_aware=True, placement_strategy="cart").fit(X, y)
+    interior = transformer.knots_[0][transformer.degree + 1 : -(transformer.degree + 1)]
+    assert len(interior) == output_dim - transformer.degree - 1
+    assert np.abs(interior - 5.0).min() < 0.1, interior
+
+
+@pytest.mark.parametrize("cls", [CubicRegressionSplineTransformer, NaturalCubicSplineTransformer])
+def test_target_aware_cubic_families_keep_the_dominant_split(cls):
+    X, y = _step_data()
+    knots = cls(output_dim=4, target_aware=True, placement_strategy="cart").fit(X, y).knots_[0]
+    assert np.abs(knots - 5.0).min() < 0.1, knots
+
+
+@pytest.mark.parametrize("cls", [BSplineTransformer, CubicRegressionSplineTransformer, NaturalCubicSplineTransformer])
+def test_adaptive_spline_width_is_not_capped_at_the_legacy_window(cls):
+    rng = np.random.default_rng(0)
+    x = rng.uniform(0, 10, 3000)
+    y = np.sin(2 * x) * 3 + rng.normal(0, 0.3, 3000)
+    transformer = cls(
+        adaptive=True, min_output_dim=5, max_output_dim=40, target_aware=True, placement_strategy="cart"
+    ).fit(x.reshape(-1, 1), y)
+    assert 15 < transformer.total_output_dim_ <= 40
diff --git a/tests/expansion/spline/test_tied_knot_placement.py b/tests/expansion/spline/test_tied_knot_placement.py
new file mode 100644
index 0000000..1345bc2
--- /dev/null
+++ b/tests/expansion/spline/test_tied_knot_placement.py
@@ -0,0 +1,118 @@
+"""Knot placement on tied data (top-coded, zero-inflated, low-cardinality features).
+
+Regression guards for issue #56: quantile knots on tied data coincide with each
+other and with the range boundary, which collapsed basis functions to a point
+(dead columns) and, for B/M-splines, left every row at ``x_max`` all zero.
+"""
+
+import numpy as np
+import pytest
+
+from pretab.core.knots import supplement_interior_knots
+from pretab.transformers import (
+    BSplineTransformer,
+    CubicRegressionSplineTransformer,
+    ISplineTransformer,
+    MSplineTransformer,
+    NaturalCubicSplineTransformer,
+    TensorProductSplineTransformer,
+)
+
+
+@pytest.fixture
+def top_coded():
+    # ~40% of the rows sit exactly on the cap of 100.
+    return np.minimum(np.random.default_rng(0).uniform(0, 160, (1000, 1)), 100.0)
+
+
+@pytest.fixture
+def zero_inflated():
+    rng = np.random.default_rng(1)
+    return np.where(rng.random((1000, 1)) < 0.6, 0.0, rng.exponential(5.0, (1000, 1)))
+
+
+def _interior(knots, degree):
+    return knots[degree + 1 : len(knots) - degree - 1]
+
+
+def _dense_grid(X, n=2001):
+    return np.linspace(X.min(), X.max(), n).reshape(-1, 1)
+
+
+@pytest.mark.parametrize("cls", [BSplineTransformer, MSplineTransformer, ISplineTransformer])
+@pytest.mark.parametrize("data", ["top_coded", "zero_inflated"])
+def test_bmi_interior_knots_are_unique_and_strictly_inside(cls, data, request):
+    X = request.getfixturevalue(data)
+    transformer = cls(output_dim=6).fit(X)
+    knots = transformer.knots_[0]
+    interior = _interior(knots, transformer.degree)
+
+    assert len(interior) == 6 - transformer.degree - 1
+    assert len(np.unique(interior)) == len(interior)
+    assert (interior > X.min()).all() and (interior < X.max()).all()
+    # Boundary multiplicity stays exactly degree + 1.
+    assert (knots == X.min()).sum() == transformer.degree + 1
+    assert (knots == X.max()).sum() == transformer.degree + 1
+    assert (np.abs(transformer.transform(_dense_grid(X))).max(axis=0) > 0).all()
+
+
+def test_bspline_is_a_partition_of_unity_at_the_cap(top_coded):
+    transformer = BSplineTransformer(output_dim=6).fit(top_coded)
+    basis = transformer.transform(top_coded)
+    np.testing.assert_allclose(basis.sum(axis=1), 1.0, atol=1e-9)
+    # Values on and beyond the cap (clipped onto x_max) keep a full basis row.
+    np.testing.assert_allclose(transformer.transform(np.array([[100.0], [150.0]])).sum(axis=1), 1.0, atol=1e-9)
+
+
+def test_ispline_equals_one_at_the_cap(top_coded):
+    transformer = ISplineTransformer(output_dim=6).fit(top_coded)
+    np.testing.assert_allclose(transformer.transform(np.array([[100.0]])), 1.0, atol=1e-12)
+
+
+@pytest.mark.parametrize("cls", [NaturalCubicSplineTransformer, CubicRegressionSplineTransformer])
+def test_cubic_families_place_unique_knots_on_zero_inflated_data(cls, zero_inflated):
+    transformer = cls(output_dim=6, placement_strategy="quantile").fit(zero_inflated)
+    knots = transformer.knots_[0]
+    assert len(np.unique(knots)) == len(knots)
+    basis = transformer.transform(_dense_grid(zero_inflated))
+    assert np.linalg.matrix_rank(basis) == basis.shape[1]
+    assert (np.abs(basis).max(axis=0) > 0).all()
+
+
+def test_natural_spline_keeps_both_endpoints_once(zero_inflated):
+    knots = NaturalCubicSplineTransformer(output_dim=6, placement_strategy="quantile").fit(zero_inflated).knots_[0]
+    assert len(knots) == 7
+    assert knots[0] == zero_inflated.min() and knots[-1] == zero_inflated.max()
+    assert (knots[1:-1] > knots[0]).all() and (knots[1:-1] < knots[-1]).all()
+
+
+def test_tensor_product_has_no_degenerate_columns_on_zero_inflated_marginal(zero_inflated):
+    X = np.hstack([zero_inflated, np.random.default_rng(2).normal(size=(len(zero_inflated), 1))])
+    transformer = TensorProductSplineTransformer(output_dim=8, placement_strategy="quantile").fit(X)
+    g0 = np.linspace(X[:, 0].min(), X[:, 0].max(), 200)
+    g1 = np.linspace(X[:, 1].min(), X[:, 1].max(), 200)
+    grid = np.array(np.meshgrid(g0, g1)).reshape(2, -1).T
+    assert (np.abs(transformer.transform(grid)).max(axis=0) > 0).all()
+
+
+@pytest.mark.parametrize("levels", [[0.0, 1.0], [0.0, 1.0, 2.0], [1.0, 2.0, 3.0, 4.0]])
+def test_low_cardinality_feature_keeps_a_full_rank_bspline_basis(levels):
+    X = np.random.default_rng(3).choice(levels, size=(300, 1))
+    transformer = BSplineTransformer(output_dim=6).fit(X)
+    basis = transformer.transform(_dense_grid(X))
+    assert np.linalg.matrix_rank(basis) == 6
+
+
+def test_supplement_interior_knots_drops_boundary_and_duplicate_knots():
+    x = np.array([0.0, 0.0, 0.0, 1.0, 2.0, 10.0, 10.0])
+    knots = supplement_interior_knots(x, np.array([0.0, 2.0, 2.0, 10.0, 12.0]), 3)
+    assert len(knots) == 3
+    assert 2.0 in knots  # an existing interior knot is always kept
+    assert (knots > 0.0).all() and (knots < 10.0).all()
+    assert len(np.unique(knots)) == 3
+
+
+def test_supplement_interior_knots_never_down_samples():
+    x = np.linspace(0, 1, 50)
+    knots = np.array([0.1, 0.2, 0.3, 0.4])
+    np.testing.assert_array_equal(supplement_interior_knots(x, knots, 2), knots)
diff --git a/tests/integration/test_column_labels.py b/tests/integration/test_column_labels.py
new file mode 100644
index 0000000..1c9829f
--- /dev/null
+++ b/tests/integration/test_column_labels.py
@@ -0,0 +1,116 @@
+"""Column labels the ColumnTransformer cannot take verbatim.
+
+Regression guards for issue #55 (integer labels were read as *positions* by
+scikit-learn's ColumnTransformer, silently swapping columns) and issue #71 (labels
+starting with ``_`` or containing ``__`` were embedded into estimator names that
+scikit-learn rejects).
+"""
+
+from typing import cast
+
+import numpy as np
+import pandas as pd
+import pytest
+from sklearn.compose import ColumnTransformer
+from sklearn.impute import SimpleImputer
+from sklearn.pipeline import make_pipeline
+
+from pretab import Preprocessor
+from pretab.exceptions import PretabDataError
+
+
+@pytest.fixture
+def int_labelled():
+    # Labels are a permutation of the positions: label 1 sits at position 0.
+    return pd.DataFrame({1: [10.0, 20, 30, 40, 50, 60], 0: [1.0, 2, 3, 4, 5, 6]})
+
+
+def _passthrough(**kwargs):
+    return Preprocessor(numerical_method="none", scaling="none", output_structure="blocks", **kwargs)
+
+
+def _blocks(out):
+    return cast("dict[str, np.ndarray]", out)
+
+
+def test_integer_labels_route_each_column_to_its_own_block(int_labelled):
+    out = _blocks(_passthrough().fit_transform(int_labelled))
+    np.testing.assert_array_equal(out["num_1"].ravel(), int_labelled[1].to_numpy())
+    np.testing.assert_array_equal(out["num_0"].ravel(), int_labelled[0].to_numpy())
+
+
+def test_feature_preprocessing_with_integer_keys_targets_the_labelled_column(int_labelled):
+    pre = Preprocessor(feature_preprocessing={1: "minmax", 0: "none"}, scaling="none", output_structure="blocks")
+    out = _blocks(pre.fit_transform(int_labelled))
+    np.testing.assert_allclose(out["num_1"].ravel(), np.linspace(0.0, 1.0, 6))
+    np.testing.assert_array_equal(out["num_0"].ravel(), int_labelled[0].to_numpy())
+
+
+def test_non_contiguous_integer_labels_fit_and_transform():
+    frame = pd.DataFrame(np.random.default_rng(0).normal(size=(20, 3))).drop(columns=[0])
+    pre = Preprocessor(numerical_method="minmax").fit(frame)
+    assert pre.get_feature_names_out().tolist() == ["num_1", "num_2"]
+    assert np.asarray(pre.transform(frame)).shape == (20, 2)
+
+
+def test_integer_labels_are_matched_by_label_at_transform(int_labelled):
+    pre = _passthrough(output_format="dense").fit(int_labelled)
+    reordered = int_labelled[[0, 1]]
+    np.testing.assert_array_equal(
+        _blocks(pre.transform(reordered))["num_1"], _blocks(pre.transform(int_labelled))["num_1"]
+    )
+
+
+def test_integer_label_metadata_uses_the_input_labels(int_labelled):
+    pre = Preprocessor(numerical_method="minmax").fit(int_labelled)
+    assert pre.output_dims_ == {1: 1, 0: 1}
+    numerical, _, _ = pre.get_feature_info(verbose=False)
+    assert list(numerical) == [1, 0]
+    assert [record.source_features for record in pre.get_feature_lineage()] == [("1",), ("0",)]
+
+
+def test_labels_colliding_as_strings_are_rejected():
+    frame = pd.DataFrame({1: [1.0, 2.0, 3.0], "1": [4.0, 5.0, 6.0]})
+    with pytest.raises(PretabDataError, match="unique when converted to strings"):
+        Preprocessor(numerical_method="minmax").fit(frame)
+
+
+def test_fit_does_not_modify_the_callers_column_labels(int_labelled):
+    Preprocessor(numerical_method="minmax").fit(int_labelled)
+    assert list(int_labelled.columns) == [1, 0]
+
+
+@pytest.mark.parametrize("label", ["_score", "income__log", "__index_level_0__", "trailing__"])
+def test_numerical_labels_with_underscores_fit(label):
+    frame = pd.DataFrame({label: np.random.default_rng(0).normal(size=50)})
+    pre = Preprocessor(numerical_method="minmax").fit(frame)
+    assert pre.get_feature_names_out().tolist() == [f"num_{label}"]
+    assert list(_blocks(pre.transform(frame, return_array=False))) == [f"num_{label}"]
+    assert pre.output_dims_ == {label: 1}
+
+
+def test_categorical_and_expanded_labels_with_underscores_fit():
+    rng = np.random.default_rng(0)
+    frame = pd.DataFrame({"_city": rng.choice(["a", "b"], 50), "x__y": rng.normal(size=50)})
+    pre = Preprocessor(categorical_method="one-hot", output_dim=3).fit(frame, rng.normal(size=50))
+    names = pre.get_feature_names_out().tolist()
+    assert names == ["num_x__y_ple0", "num_x__y_ple1", "num_x__y_ple2", "cat__city_a", "cat__city_b"]
+    assert [record.output_feature for record in pre.get_feature_lineage()] == names
+    assert pre.output_dims_ == {"x__y": 3, "_city": 2}
+
+
+def test_underscore_labels_do_not_collide_with_other_steps():
+    rng = np.random.default_rng(0)
+    # "_col0" would map to the fallback name of the first step without collision handling.
+    frame = pd.DataFrame({"_a": rng.normal(size=30), "col0": rng.normal(size=30)})
+    pre = Preprocessor(numerical_method="minmax").fit(frame)
+    assert pre.get_feature_names_out().tolist() == ["num__a", "num_col0"]
+
+
+def test_preprocessor_chains_after_a_pandas_column_transformer():
+    rng = np.random.default_rng(0)
+    X = pd.DataFrame({"age": rng.normal(40, 10, 50), "income": rng.normal(50, 10, 50)})
+    upstream = ColumnTransformer([("imp", SimpleImputer(), ["age", "income"])]).set_output(transform="pandas")
+    pipe = make_pipeline(upstream, Preprocessor(numerical_method="minmax")).fit(X, rng.normal(size=50))
+    assert cast(Preprocessor, pipe[-1]).get_feature_names_out().tolist() == ["num_imp__age", "num_imp__income"]
+    assert np.asarray(pipe.transform(X)).shape == (50, 2)
diff --git a/tests/integration/test_missing_policy.py b/tests/integration/test_missing_policy.py
index c7ec70e..ac728a3 100644
--- a/tests/integration/test_missing_policy.py
+++ b/tests/integration/test_missing_policy.py
@@ -262,3 +262,57 @@ def test_add_missing_indicator_no_longer_requires_imputation_enabled(frame_with_
     assert out.size > 0
     assert np.isnan(out).any()
     assert p.total_output_dim_ > 0
+
+
+# --- imputer indicator stays outside the representation (issue #62) -----------
+
+
+def _indicator_values(pre, frame, name):
+    names = list(pre.get_feature_names_out())
+    return np.asarray(pre.transform(frame, return_array=True), dtype=float)[:, names.index(name)]
+
+
+@pytest.mark.parametrize("options", [{"add_missing_indicator": True}, {"missing_policy": "impute_with_indicator"}])
+def test_imputer_indicator_is_a_raw_binary_column(frame_with_nan, y, options):
+    """Regression guard for issue #62: the indicator was standardized / expanded
+    by the representation like a second input feature."""
+    pre = Preprocessor(numerical_method="standardization", **options).fit(frame_with_nan, y)
+    values = _indicator_values(pre, frame_with_nan, "num_a__missingindicator_a")
+    np.testing.assert_array_equal(values, frame_with_nan["a"].isna().astype(float))
+
+
+@pytest.mark.parametrize("method", ["bspline", "rbf", "pspline", "fourier", "cubicspline", "polynomial"])
+def test_imputer_indicator_adds_exactly_one_column_per_feature_with_missing(frame_with_nan, y, method):
+    plain = Preprocessor(numerical_method=method, output_dim=5, missing_policy="impute").fit(frame_with_nan, y)
+    pre = Preprocessor(numerical_method=method, output_dim=5, missing_policy="impute_with_indicator").fit(
+        frame_with_nan, y
+    )
+    names = list(pre.get_feature_names_out())
+    indicators = [name for name in names if "missingindicator" in name]
+    assert indicators == ["num_a__missingindicator_a", "num_b__missingindicator_b"]
+    assert [name for name in names if name not in indicators] == list(plain.get_feature_names_out())
+    assert pre.output_dims_ == {key: width + 1 for key, width in plain.output_dims_.items()}
+
+
+@pytest.mark.parametrize("method", ["bspline", "rbf", "pspline", "fourier", "cubicspline"])
+def test_imputer_indicator_lineage_and_feature_info(frame_with_nan, y, method):
+    pre = Preprocessor(numerical_method=method, output_dim=5, add_missing_indicator=True).fit(frame_with_nan, y)
+    lineage = pre.get_feature_lineage()
+    assert [record.output_feature for record in lineage] == list(pre.get_feature_names_out())
+    missing = [record for record in lineage if record.family == "missing_state"]
+    assert [record.output_feature for record in missing] == ["num_a__missingindicator_a", "num_b__missingindicator_b"]
+    numerical, _, _ = pre.get_feature_info(verbose=False)
+    assert {feature: info["dimension"] for feature, info in numerical.items()} == pre.output_dims_
+
+
+def test_imputer_indicator_on_categorical_integer_codes(y):
+    frame = pd.DataFrame({"c": ["x", "y", np.nan, "x", "y", "x"]})
+    pre = Preprocessor(categorical_method="int", add_missing_indicator=True).fit(frame, y)
+    assert list(pre.get_feature_names_out()) == ["cat_c", "cat_c__missingindicator_c"]
+    np.testing.assert_array_equal(_indicator_values(pre, frame, "cat_c__missingindicator_c"), frame["c"].isna())
+
+
+def test_imputer_indicator_only_marks_features_with_missing_values_at_fit(clean_frame, frame_with_nan, y):
+    pre = Preprocessor(numerical_method="minmax", add_missing_indicator=True).fit(clean_frame, y)
+    assert not any("missing" in name for name in pre.get_feature_names_out())
+    assert np.asarray(pre.transform(frame_with_nan, return_array=True)).shape == (6, 2)
diff --git a/tests/integration/test_output_format.py b/tests/integration/test_output_format.py
index abe2613..b16ccd3 100644
--- a/tests/integration/test_output_format.py
+++ b/tests/integration/test_output_format.py
@@ -5,6 +5,8 @@
 wrapping.
 """
 
+from typing import cast
+
 import numpy as np
 import pandas as pd
 import pytest
@@ -241,3 +243,36 @@ def test_invalid_output_format_raises(frame, y):
 
     with pytest.raises(InvalidParamError):
         _bspline(output_format="nope").fit(frame, y)
+
+
+# --- "auto" is resolved once per fit (issue #73) ----------------------------------
+
+
+def test_auto_format_is_fixed_at_fit_for_dense_training_output():
+    """Regression guard for issue #73: low-density single rows came back sparse."""
+    from sklearn.datasets import load_diabetes
+
+    X, y = cast("tuple[pd.DataFrame, pd.Series]", load_diabetes(return_X_y=True, as_frame=True))
+    pre = Preprocessor(output_format="auto", random_state=0).fit(X, y)
+    assert pre.output_format_ == "dense"
+    for i in range(len(X)):
+        assert isinstance(pre.transform(X.iloc[[i]]), np.ndarray)
+    # The report still describes each call's own density.
+    pre.transform(X.iloc[[0]])
+    assert pre.output_report_["format"] == "dense"
+    assert pre.output_report_["shape"] == (1, pre.total_output_dim_)
+
+
+def test_auto_format_is_fixed_at_fit_for_sparse_training_output():
+    rng = np.random.default_rng(0)
+    frame = pd.DataFrame({"city": rng.choice([f"c{i}" for i in range(20)], 200)})
+    pre = Preprocessor(categorical_method="one-hot", output_format="auto").fit(frame)
+    assert pre.output_format_ == "sparse"
+    assert sp.issparse(pre.transform(frame.iloc[:1]))
+    assert sp.issparse(pre.transform(frame))
+
+
+@pytest.mark.parametrize("output_format", ["dense", "sparse"])
+def test_explicit_formats_are_stored_as_resolved(frame, y, output_format):
+    pre = Preprocessor(numerical_method="minmax", output_format=output_format).fit(frame, y)
+    assert pre.output_format_ == output_format
diff --git a/tests/integration/test_preprocessor.py b/tests/integration/test_preprocessor.py
index df9d224..62148ef 100644
--- a/tests/integration/test_preprocessor.py
+++ b/tests/integration/test_preprocessor.py
@@ -377,3 +377,165 @@ def test_unknown_feature_preprocessing_key_raises(sample_data):
     X, y = sample_data
     with pytest.raises(InvalidParamError, match="feature_preprocessing"):
         Preprocessor(feature_preprocessing={"nnum1": "minmax"}).fit(X, y)
+
+
+# --- boolean columns (issue #54) ------------------------------------------------
+
+
+@pytest.fixture
+def bool_frame():
+    return pd.DataFrame({"x": [0.1, 0.4, 0.2, 0.9, 0.5, 0.7], "flag": [True, False, True, False, True, True]})
+
+
+@pytest.mark.parametrize(
+    ("options", "expected"),
+    [
+        ({}, ["cat_flag"]),
+        ({"categorical_method": "one-hot"}, ["cat_flag_False", "cat_flag_True"]),
+        ({"categorical_method": "none"}, ["cat_flag"]),
+        ({"preset": "expanded"}, ["cat_flag_False", "cat_flag_True"]),
+        ({"missing_policy": "impute_with_indicator"}, ["cat_flag"]),
+        ({"missing_policy": "separate_state"}, ["cat_flag__representation__flag", "cat_flag__missing__flag__missing"]),
+    ],
+)
+def test_bool_column_is_encoded_as_a_binary_categorical(bool_frame, options, expected):
+    """Regression guard for issue #54: SimpleImputer rejected the bool dtype."""
+    pre = Preprocessor(**options).fit(bool_frame, np.arange(6.0))
+    names = pre.get_feature_names_out().tolist()
+    assert names[-len(expected) :] == expected
+    assert pre.categorical_features_ == ["flag"]
+    assert np.asarray(pre.transform(bool_frame)).shape == (6, len(names))
+
+
+def test_bool_column_integer_codes_follow_the_values(bool_frame):
+    pre = Preprocessor(categorical_method="int").fit(bool_frame, np.arange(6.0))
+    codes = np.asarray(pre.transform(bool_frame))[:, -1]
+    np.testing.assert_array_equal(codes, np.where(bool_frame["flag"], 2.0, 1.0))
+
+
+def test_nullable_boolean_column_with_missing_values_is_imputed(bool_frame):
+    frame = bool_frame.astype({"flag": "boolean"})
+    frame.loc[2, "flag"] = pd.NA
+    pre = Preprocessor(categorical_method="one-hot").fit(frame, np.arange(6.0))
+    one_hot = np.asarray(pre.transform(frame))[:, -2:]
+    np.testing.assert_array_equal(one_hot[2], [0.0, 1.0])  # imputed with the most frequent value, True
+
+
+def test_bool_numpy_array_is_supported():
+    X = np.array([[True], [False], [True], [True]])
+    pre = Preprocessor().fit(X)
+    assert pre.categorical_features_ == ["feature_0"]
+    assert np.asarray(pre.transform(X)).shape == (4, 1)
+
+
+# --- box-cox on unseen low values (issue #59) -------------------------------------
+
+
+@pytest.fixture
+def boxcox_fitted():
+    X = pd.DataFrame({"x": np.random.default_rng(0).normal(10, 2, 200)})
+    return X, Preprocessor(numerical_method="box-cox").fit(X)
+
+
+def test_boxcox_transforms_values_below_the_training_minimum(boxcox_fitted):
+    """Regression guard for issue #59: the positivity scaler mapped them to <= 0."""
+    X, pre = boxcox_fitted
+    low = pre.transform(pd.DataFrame({"x": [X.x.min() - 5.0, X.x.min() - 0.5]}))
+    at_min = pre.transform(pd.DataFrame({"x": [X.x.min()]}))
+    assert np.isfinite(low).all()
+    np.testing.assert_allclose(low, np.repeat(at_min, 2, axis=0))
+
+
+def test_boxcox_still_extrapolates_above_the_training_maximum(boxcox_fitted):
+    X, pre = boxcox_fitted
+    out = pre.transform(pd.DataFrame({"x": [X.x.max(), X.x.max() + 3.0]})).ravel()
+    assert np.isfinite(out).all()
+    assert out[1] > out[0]
+
+
+def test_boxcox_survives_cross_validation():
+    from sklearn.datasets import make_regression
+    from sklearn.linear_model import Ridge
+    from sklearn.model_selection import cross_val_score
+    from sklearn.pipeline import make_pipeline
+
+    data = make_regression(n_samples=500, n_features=3, noise=1.0, random_state=0)
+    X, y = data[0], data[1]
+    scores = cross_val_score(make_pipeline(Preprocessor(numerical_method="box-cox"), Ridge()), X, y, cv=5)
+    assert np.isfinite(scores).all()
+
+
+# --- ndarray width check at transform (issue #72) ---------------------------------
+
+
+@pytest.mark.parametrize("width", [2, 4])
+def test_transform_rejects_an_ndarray_of_a_different_width(width):
+    """Regression guard for issue #72: an extra array column silently shifted features."""
+    rng = np.random.default_rng(0)
+    X, y = rng.normal(size=(60, 3)), rng.normal(size=60)
+    pre = Preprocessor(random_state=0).fit(X, y)
+    with pytest.raises(ValueError, match=f"X has {width} features, but Preprocessor is expecting 3"):
+        pre.transform(rng.normal(size=(5, width)))
+
+
+def test_transform_accepts_an_ndarray_of_the_fitted_width():
+    rng = np.random.default_rng(0)
+    X, y = rng.normal(size=(60, 3)), rng.normal(size=60)
+    pre = Preprocessor(random_state=0).fit(X, y)
+    assert np.asarray(pre.transform(X)).shape[0] == 60
+
+
+# --- get_feature_names_out(input_features) (issue #65) ----------------------------
+
+
+@pytest.fixture
+def named_numeric():
+    rng = np.random.default_rng(0)
+    return pd.DataFrame({"age": rng.normal(40, 10, 60), "income": rng.normal(50, 10, 60)}), rng.normal(size=60)
+
+
+@pytest.mark.parametrize(
+    "options",
+    [
+        {"numerical_method": "minmax"},
+        {"numerical_method": "polynomial"},
+        {"numerical_method": "bspline", "output_dim": 5, "missing_policy": "separate_state"},
+        {"numerical_method": "ple", "output_dim": 3, "add_missing_indicator": True},
+    ],
+)
+def test_array_fit_uses_the_given_input_features(named_numeric, options):
+    """Regression guard for issue #65: input_features had to equal the synthetic names."""
+    X, y = named_numeric
+    X = X.copy()
+    X.iloc[:3, 0] = np.nan
+    from_array = Preprocessor(**options).fit(X.to_numpy(), y).get_feature_names_out(["age", "income"])
+    from_frame = Preprocessor(**options).fit(X, y).get_feature_names_out()
+    assert from_array.tolist() == from_frame.tolist()
+
+
+def test_array_fit_rejects_input_features_of_the_wrong_length(named_numeric):
+    X, y = named_numeric
+    pre = Preprocessor(numerical_method="minmax").fit(X.to_numpy(), y)
+    with pytest.raises(ValueError, match="length equal to number of features"):
+        pre.get_feature_names_out(["age"])
+
+
+def test_named_fit_records_and_validates_feature_names_in(named_numeric):
+    X, y = named_numeric
+    pre = Preprocessor(numerical_method="minmax").fit(X, y)
+    assert pre.feature_names_in_.tolist() == ["age", "income"]
+    assert pre.get_feature_names_out(["age", "income"]).tolist() == ["num_age", "num_income"]
+    with pytest.raises(ValueError, match="not equal to feature_names_in_"):
+        pre.get_feature_names_out(["a", "b"])
+    assert not hasattr(pre.fit(X.to_numpy(), y), "feature_names_in_")
+
+
+def test_pipeline_feature_names_propagate_through_an_array_step(named_numeric):
+    from sklearn.impute import SimpleImputer
+    from sklearn.pipeline import make_pipeline
+
+    X, y = named_numeric
+    pipe = make_pipeline(SimpleImputer(), Preprocessor(numerical_method="minmax")).fit(X, y)
+    names = pipe.get_feature_names_out()
+    assert names is not None
+    assert names.tolist() == ["num_age", "num_income"]
diff --git a/tests/integration/test_serialization.py b/tests/integration/test_serialization.py
index 493000c..77530fb 100644
--- a/tests/integration/test_serialization.py
+++ b/tests/integration/test_serialization.py
@@ -274,3 +274,67 @@ def test_from_spec_refuses_dataclass_not_in_allowlist():
 def test_from_spec_rejects_bad_source_type():
     with pytest.raises(PretabSerializationError):
         Preprocessor.from_spec(12345)
+
+
+# --- object arrays and builtin dtypes (issue #66) ---------------------------------
+
+
+@pytest.fixture
+def numpy_scalar_frame():
+    rng = np.random.default_rng(0)
+    n = 80
+    income = np.where(rng.random(n) < 0.1, np.nan, rng.normal(50, 10, n))
+    frame = pd.DataFrame({"income": income, "age": rng.normal(40, 10, n)})
+    # An object column of np.bool_ / None, as produced by a row-wise apply.
+    frame["high_income"] = frame.apply(lambda r: r.income > 50 if pd.notna(r.income) else None, axis=1)
+    return frame, rng.normal(size=n)
+
+
+def test_object_arrays_of_numpy_scalars_serialize(numpy_scalar_frame, tmp_path):
+    """Regression guard for issue #66: numpy scalars in object arrays leaked into the spec."""
+    frame, y = numpy_scalar_frame
+    pre = Preprocessor(random_state=0).fit(frame, y)
+    json.dumps(pre.to_spec())
+    assert len(pre.fingerprint_) == 64
+    path = tmp_path / "spec.json"
+    pre.to_spec(path)
+    np.testing.assert_array_equal(Preprocessor.from_spec(path).transform(frame), pre.transform(frame))
+
+
+@pytest.mark.parametrize("dtype", [float, int, np.float32])
+def test_builtin_and_numpy_dtypes_round_trip(dtype):
+    frame = pd.DataFrame({"age": np.linspace(20.0, 60.0, 30)})
+    pre = Preprocessor(numerical_method="minmax", dtype=dtype).fit(frame)
+    loaded = Preprocessor.from_spec(json.loads(json.dumps(pre.to_spec())))
+    np.testing.assert_array_equal(loaded.transform(frame), pre.transform(frame))
+    assert loaded.fingerprint_ == pre.fingerprint_
+
+
+def test_failed_save_leaves_an_existing_spec_untouched(tmp_path):
+    frame = pd.DataFrame({"age": np.linspace(20.0, 60.0, 30)})
+    pre = Preprocessor(numerical_method="minmax").fit(frame)
+    path = tmp_path / "spec.json"
+    pre.to_spec(path)
+    saved = path.read_text(encoding="utf-8")
+
+    vars(pre)["unsupported_"] = lambda value: value
+    with pytest.raises(PretabSerializationError):
+        pre.to_spec(path)
+    assert path.read_text(encoding="utf-8") == saved
+
+
+def test_types_a_spec_cannot_load_are_rejected_at_save_time():
+    from pretab.compose.serialize import _encode
+
+    with pytest.raises(PretabSerializationError, match="can only reference types"):
+        _encode(list)
+
+
+def test_object_arrays_written_before_element_encoding_still_load():
+    from pretab.compose.serialize import _decode
+
+    payload = {"__ndarray__": {"dtype": "|O", "shape": [2, 2], "data": [["a", None], [True, 1.5]]}}
+    decoded = _decode(payload)
+    assert isinstance(decoded, np.ndarray)
+    assert decoded.dtype == object
+    assert decoded.tolist() == [["a", None], [True, 1.5]]
diff --git a/tests/integration/test_verbosity.py b/tests/integration/test_verbosity.py
index 3c44e1e..8b44f78 100644
--- a/tests/integration/test_verbosity.py
+++ b/tests/integration/test_verbosity.py
@@ -6,6 +6,7 @@
 :class:`~pretab.PretabWarning`.
 """
 
+import contextlib
 import io
 import logging
 import warnings
@@ -186,9 +187,30 @@ def test_set_verbosity_sets_logger_level():
     assert logger.level == logging.INFO
 
 
+@contextlib.contextmanager
+def _root_handlers(*handlers):
+    """Run with exactly ``handlers`` on the root logger (none: a plain Python process).
+
+    pytest adds its capture handlers to the root logger when a test body starts,
+    so this is used inside the test body; the original handlers are restored after.
+    """
+    root = logging.getLogger()
+    saved = root.handlers[:]
+    root.handlers = list(handlers)
+    try:
+        yield
+    finally:
+        root.handlers = saved
+
+
+def _own_handlers():
+    return [h for h in logging.getLogger("pretab").handlers if not isinstance(h, logging.NullHandler)]
+
+
 def test_configure_logging_attaches_stream_handler_when_none():
     logger = logging.getLogger("pretab")
-    configure_logging(1)
+    with _root_handlers():
+        configure_logging(1)
     assert any(isinstance(h, logging.StreamHandler) and not isinstance(h, logging.NullHandler) for h in logger.handlers)
     assert logger.level == logging.INFO
 
@@ -205,6 +227,55 @@ def test_configure_logging_respects_existing_handler():
     assert logger.level == logging.CRITICAL
 
 
+def test_configure_logging_raises_the_level_on_later_calls():
+    """Regression guard for issue #67: the first call pinned the verbosity."""
+    with _root_handlers():
+        configure_logging(1)
+        configure_logging(2)
+    assert logging.getLogger("pretab").level == logging.DEBUG
+    assert len(_own_handlers()) == 1  # the handler is reused, never duplicated
+
+
+def test_fit_verbose_2_after_verbose_1_logs_the_feature_table(sample_data, capsys):
+    X, y = sample_data
+    with _root_handlers():
+        Preprocessor(numerical_method="ple", verbose=1).fit(X, y)
+        capsys.readouterr()
+        Preprocessor(numerical_method="ple", verbose=2).fit(X, y)
+    err = capsys.readouterr().err
+    assert "fit complete" in err
+    assert "pipeline" in err
+
+
+def test_get_feature_info_does_not_lower_a_debug_level(sample_data):
+    X, y = sample_data
+    pre = Preprocessor(numerical_method="ple").fit(X, y)
+    with _root_handlers():
+        configure_logging(2)
+        pre.get_feature_info(verbose=True)
+    assert logging.getLogger("pretab").level == logging.DEBUG
+
+
+def test_configure_logging_defers_to_a_root_handler(sample_data, capsys):
+    """Regression guard for issue #67: a host root handler printed every line twice."""
+    X, y = sample_data
+    stream = io.StringIO()  # a host root handler, as logging.basicConfig installs
+    with _root_handlers(logging.StreamHandler(stream)):
+        Preprocessor(numerical_method="ple", verbose=1).fit(X, y)
+        assert _own_handlers() == []
+    assert stream.getvalue().count("fit complete") == 1
+    assert "fit complete" not in capsys.readouterr().err
+
+
+def test_configure_logging_hands_over_to_a_host_configured_later():
+    with _root_handlers():
+        configure_logging(1)
+    assert len(_own_handlers()) == 1
+    with _root_handlers(logging.StreamHandler(io.StringIO())):
+        configure_logging(1)
+    assert _own_handlers() == []
+
+
 # --------------------------------------------------------------------------- #
 # warning categories (PretabWarning family)
 # --------------------------------------------------------------------------- #
diff --git a/tests/placement/test_spline_placement_adapter.py b/tests/placement/test_spline_placement_adapter.py
index 5d24fa1..52a1350 100644
--- a/tests/placement/test_spline_placement_adapter.py
+++ b/tests/placement/test_spline_placement_adapter.py
@@ -99,3 +99,18 @@ def test_lightgbm_adapter_reproducible(data):
     a = SplinePlacementAdapter(placement_strategy="lightgbm", degree=3).get_knot_locations(X, y)
     b = SplinePlacementAdapter(placement_strategy="lightgbm", degree=3).get_knot_locations(X, y)
     np.testing.assert_array_equal(a, b)
+
+
+def _step_data(seed=1):
+    rng = np.random.default_rng(seed)
+    x = rng.uniform(0, 10, 1000)
+    y = np.where(x > 5.0, 1.0, 0.0) + 0.3 * rng.normal(size=1000)
+    return x.reshape(-1, 1), y
+
+
+def test_get_knot_locations_window_override():
+    X, y = _step_data()
+    adapter = SplinePlacementAdapter(placement_strategy="cart", degree=3)
+    knots = adapter.get_knot_locations(X, y, min_knots=1, max_knots=1)
+    assert len(knots) == 1
+    assert abs(knots[0] - 5.0) < 0.1