diff --git a/foundry/__init__.py b/foundry/__init__.py
index 44b1806..407b8a2 100644
--- a/foundry/__init__.py
+++ b/foundry/__init__.py
@@ -1 +1 @@
-__version__ = '0.2.6'
+__version__ = '0.2.7'
diff --git a/foundry/uplift/__init__.py b/foundry/uplift/__init__.py
new file mode 100644
index 0000000..4d03b2e
--- /dev/null
+++ b/foundry/uplift/__init__.py
@@ -0,0 +1 @@
+from .meta_learners import SLearner, TLearner, XLearner
diff --git a/foundry/uplift/meta_learners/__init__.py b/foundry/uplift/meta_learners/__init__.py
new file mode 100644
index 0000000..c224da8
--- /dev/null
+++ b/foundry/uplift/meta_learners/__init__.py
@@ -0,0 +1,3 @@
+from .t_learner import TLearner
+from .s_learner import SLearner
+from .x_learner import XLearner
diff --git a/foundry/uplift/meta_learners/s_learner.py b/foundry/uplift/meta_learners/s_learner.py
new file mode 100644
index 0000000..121f914
--- /dev/null
+++ b/foundry/uplift/meta_learners/s_learner.py
@@ -0,0 +1,114 @@
+import warnings
+from typing import Tuple, Union, overload, Literal, Optional
+
+import numpy as np
+import pandas as pd
+from sklearn.base import BaseEstimator, clone
+
+from foundry.util import SliceDict, to_1d, safe_predict
+from ..util import get_qini_curve, get_cumulative_gain_score
+
+
+class SLearner(BaseEstimator):
+    """
+    An S-learner. Please note that current implementation assumes randomized treatment/control!
+
+    :param estimator: Any instance that supports the sklearn API (fit/predict and can call ``clone()`` on it).
+    :param include_interaction: Whether to include X * treatment interaction terms.
+    """
+    estimator_: Optional[BaseEstimator] = None
+
+    def __init__(self, estimator: BaseEstimator, include_interaction: bool = True) -> None:
+        self.estimator = estimator
+        self.include_interaction = include_interaction
+
+    def fit(self, X: Union[pd.DataFrame, np.ndarray], y: SliceDict, **fit_kwargs) -> "SLearner":
+        y_arr, treatment_ind = self._normalize_y(y)
+
+        X_aug = self._augment_with_treatment(X, treatment_ind)
+
+        self.estimator_ = clone(self.estimator).fit(X=X_aug, y=y_arr, **fit_kwargs)
+        return self
+
+
+    @overload
+    def predict(self, X: Union[pd.DataFrame, np.ndarray], return_components: Literal[False] = ..., **predict_kwargs) -> np.ndarray: ...
+    @overload
+    def predict(self, X: Union[pd.DataFrame, np.ndarray], return_components: Literal[True], **predict_kwargs) -> Tuple[np.ndarray, np.ndarray]: ...
+
+    def predict(self, X, return_components=False, **predict_kwargs):
+        X_t = self._augment_with_treatment(X, np.ones(len(X), dtype=bool))
+        X_c = self._augment_with_treatment(X, np.zeros(len(X), dtype=bool))
+
+        yhat_t = safe_predict(self.estimator_, X=X_t, **predict_kwargs)
+        yhat_c = safe_predict(self.estimator_, X=X_c, **predict_kwargs)
+
+        if return_components:
+            return yhat_t, yhat_c
+
+        return yhat_t - yhat_c
+
+    def score(
+        self,
+        X: Union[pd.DataFrame, np.ndarray],
+        y: SliceDict,
+        sample_weight: Union[np.ndarray, None] = None,
+        method: str = 'qini',
+        normalize: bool = True,
+        **kwargs,
+    ) -> float:
+        y_arr, treatment_ind = self._normalize_y(y)
+        if sample_weight is not None:
+            raise NotImplementedError
+        pred = self.predict(X=X)
+        if method == 'cumulative_gain':
+            return get_cumulative_gain_score(
+                y_true=y_arr,
+                treatment=treatment_ind,
+                score=pred,
+                **kwargs,
+            )
+        qini = get_qini_curve(
+            y_true=y_arr,
+            treatment=treatment_ind,
+            score=pred,
+            normalize=normalize,
+            **kwargs,
+        )
+        random_area = np.linspace(0, qini[-1], qini.shape[0]).sum()
+        return (np.nansum(qini) - random_area) / qini.shape[0]
+
+    def _augment_with_treatment(
+        self,
+        X: Union[pd.DataFrame, np.ndarray],
+        treatment_ind: np.ndarray,
+    ) -> Union[pd.DataFrame, np.ndarray]:
+        treatment_col = treatment_ind.astype(int)
+
+        if isinstance(X, pd.DataFrame):
+            X_aug = X.copy()
+            X_aug["treatment"] = treatment_col
+
+            if self.include_interaction:
+                for col in X.columns:
+                    X_aug[f"{col}_x_treatment"] = X[col] * treatment_col
+
+            return X_aug
+
+        else:
+            treatment_col = treatment_col.reshape(-1, 1)
+
+            if self.include_interaction:
+                interaction_terms = X * treatment_col
+                return np.hstack([X, treatment_col, interaction_terms])
+
+            return np.hstack([X, treatment_col])
+
+    @staticmethod
+    def _normalize_y(y: SliceDict) -> Tuple[np.ndarray, np.ndarray]:
+        y = y.copy()
+        y_arr = to_1d(np.asanyarray(y.pop("value")))
+        treatment_ind = to_1d(np.asanyarray(y.pop("is_treatment")).astype(bool))
+        if len(y.keys()):
+            warnings.warn(f"Unused keys in ``y``: {set(y)}")
+        return y_arr, treatment_ind
diff --git a/foundry/uplift/meta_learners/t_learner.py b/foundry/uplift/meta_learners/t_learner.py
new file mode 100644
index 0000000..4e50f46
--- /dev/null
+++ b/foundry/uplift/meta_learners/t_learner.py
@@ -0,0 +1,80 @@
+import warnings
+from typing import Optional, Tuple, Union, overload, Literal
+
+import numpy as np
+import pandas as pd
+from sklearn.base import BaseEstimator, clone
+
+from foundry.util import SliceDict, to_1d, safe_predict
+from ..util import get_qini_curve, get_cumulative_gain_score
+
+
+class TLearner(BaseEstimator):
+    """
+    A T-learner. Please note that current implementation assumes randomized treatment/control!
+
+    :param estimator: Any instance that supports the sklearn API (fit/predict and can call ``clone()`` on it).
+    """
+    treatment_est_: Optional[BaseEstimator] = None
+    control_est_: Optional[BaseEstimator] = None
+
+    def __init__(self, estimator: BaseEstimator) -> None:
+        self.estimator = estimator
+
+    def fit(self, X: Union[pd.DataFrame, np.ndarray], y: SliceDict, **fit_kwargs) -> "TLearner":
+        y_arr, treatment_ind = self._normalize_y(y)
+        self.treatment_est_ = clone(self.estimator).fit(X=X[treatment_ind], y=y_arr[treatment_ind], **fit_kwargs)
+        self.control_est_ = clone(self.estimator).fit(X[~treatment_ind], y_arr[~treatment_ind], **fit_kwargs)
+
+        return self
+
+    @overload
+    def predict(self, X: Union[pd.DataFrame, np.ndarray], return_components: Literal[False] = ..., **predict_kwargs) -> np.ndarray: ...
+    @overload
+    def predict(self, X: Union[pd.DataFrame, np.ndarray], return_components: Literal[True], **predict_kwargs) -> Tuple[np.ndarray, np.ndarray]: ...
+
+    def predict(self, X, return_components=False, **predict_kwargs):
+        yhat_t = safe_predict(self.treatment_est_, X=X, **predict_kwargs)
+        yhat_c = safe_predict(self.control_est_, X=X, **predict_kwargs)
+        if return_components:
+            return yhat_t, yhat_c
+        return yhat_t - yhat_c
+
+    def score(
+        self,
+        X: Union[pd.DataFrame, np.ndarray],
+        y: SliceDict,
+        sample_weight: Optional[np.ndarray] = None,
+        method: str = 'qini',
+        normalize: bool = True,
+        **kwargs,
+    ) -> float:
+        y_arr, treatment_ind = self._normalize_y(y)
+        if sample_weight is not None:
+            raise NotImplementedError
+        pred = self.predict(X=X)
+        if method == 'cumulative_gain':
+            return get_cumulative_gain_score(
+                y_true=y_arr,
+                treatment=treatment_ind,
+                score=pred,
+                **kwargs,
+            )
+        qini = get_qini_curve(
+            y_true=y_arr,
+            treatment=treatment_ind,
+            score=pred,
+            normalize=normalize,
+            **kwargs
+        )
+        random_area = np.linspace(0, qini[-1], qini.shape[0]).sum()
+        return (np.nansum(qini) - random_area) / qini.shape[0]
+
+    @staticmethod
+    def _normalize_y(y: SliceDict) -> Tuple[np.ndarray, np.ndarray]:
+        y = y.copy()
+        y_arr = to_1d(np.asanyarray(y.pop('value')))
+        treatment_ind = to_1d(np.asanyarray(y.pop('is_treatment')).astype(bool))
+        if len(y.keys()):
+            warnings.warn(f"Unused keys in ``y``: {set(y)}")
+        return y_arr, treatment_ind
diff --git a/foundry/uplift/meta_learners/x_learner.py b/foundry/uplift/meta_learners/x_learner.py
new file mode 100644
index 0000000..89e19b5
--- /dev/null
+++ b/foundry/uplift/meta_learners/x_learner.py
@@ -0,0 +1,133 @@
+import warnings
+from typing import Any, Dict, Optional, Tuple, Union, overload, Literal
+
+import numpy as np
+import pandas as pd
+from sklearn.base import BaseEstimator, clone
+
+from foundry.util import SliceDict, to_1d, safe_predict
+from ..util import get_qini_curve, get_cumulative_gain_score
+
+
+class XLearner(BaseEstimator):
+    """
+    An X-learner. Please note that current implementation assumes randomized treatment/control!
+    Adapted from https://matheusfacure.github.io/python-causality-handbook/21-Meta-Learners.html 
+    and the original paper by Kunzel et al. (2019): https://arxiv.org/abs/1706.03461
+
+    :param first_stage_estimator: Any instance that supports the sklearn API (fit/predict and can call ``clone()`` on it).
+    :param second_stage_estimator: Any instance that supports the sklearn API (fit/predict and can call ``clone()`` on it).
+    :param propensity_estimator: Any instance that supports the sklearn API (fit/predict and can call ``clone()`` on it).
+    :param first_stage_fit_params: Optional dict of kwargs passed to ``fit()`` for the first-stage models.
+    :param second_stage_fit_params: Optional dict of kwargs passed to ``fit()`` for the second-stage models.
+    :param propensity_fit_params: Optional dict of kwargs passed to ``fit()`` for the propensity model.
+    """
+    first_treatment_est_: Optional[BaseEstimator] = None
+    first_control_est_: Optional[BaseEstimator] = None
+    second_treatment_est_: Optional[BaseEstimator] = None
+    second_control_est_: Optional[BaseEstimator] = None
+    propensity_est_: Optional[BaseEstimator] = None
+
+    def __init__(
+        self,
+        first_stage_estimator: BaseEstimator,
+        second_stage_estimator: BaseEstimator,
+        propensity_estimator: BaseEstimator,
+        first_stage_fit_params: Optional[Dict[str, Any]] = None,
+        second_stage_fit_params: Optional[Dict[str, Any]] = None,
+        propensity_fit_params: Optional[Dict[str, Any]] = None,
+    ) -> None:
+        self.first_stage_estimator = first_stage_estimator
+        self.second_stage_estimator = second_stage_estimator
+        self.propensity_estimator = propensity_estimator
+        self.first_stage_fit_params = first_stage_fit_params
+        self.second_stage_fit_params = second_stage_fit_params
+        self.propensity_fit_params = propensity_fit_params
+
+    def fit(self, X: Union[pd.DataFrame, np.ndarray], y: SliceDict) -> "XLearner":
+        y_arr, treatment_ind = self._normalize_y(y)
+
+        _first_stage_fit_params = self.first_stage_fit_params or {}
+        _second_stage_fit_params = self.second_stage_fit_params or {}
+        _propensity_fit_params = self.propensity_fit_params or {}
+
+        self.first_control_est_ = clone(self.first_stage_estimator).fit(
+            X[~treatment_ind], y_arr[~treatment_ind], **_first_stage_fit_params
+        )
+        self.first_treatment_est_ = clone(self.first_stage_estimator).fit(
+            X[treatment_ind], y_arr[treatment_ind], **_first_stage_fit_params
+        )
+
+        self.propensity_est_ = clone(self.propensity_estimator).fit(
+            X, treatment_ind, **_propensity_fit_params
+        )
+
+        imputed_te = np.where(
+            treatment_ind,
+            y_arr - safe_predict(self.first_control_est_, X),
+            safe_predict(self.first_treatment_est_, X) - y_arr,
+        )
+
+        self.second_control_est_ = clone(self.second_stage_estimator).fit(
+            X[~treatment_ind], imputed_te[~treatment_ind], **_second_stage_fit_params
+        )
+        self.second_treatment_est_ = clone(self.second_stage_estimator).fit(
+            X[treatment_ind], imputed_te[treatment_ind], **_second_stage_fit_params
+        )
+
+        return self
+
+    @overload
+    def predict(self, X: Union[pd.DataFrame, np.ndarray], return_components: Literal[False] = ..., **predict_kwargs) -> np.ndarray: ...
+    @overload
+    def predict(self, X: Union[pd.DataFrame, np.ndarray], return_components: Literal[True], **predict_kwargs) -> Tuple[np.ndarray, np.ndarray]: ...
+
+    def predict(self, X, return_components=False, **predict_kwargs):
+        p_treatment = safe_predict(self.propensity_est_, X)
+        p_control = 1 - p_treatment
+
+        tau0 = safe_predict(self.second_control_est_,   X, **predict_kwargs)
+        tau1 = safe_predict(self.second_treatment_est_, X, **predict_kwargs)
+
+        if return_components:
+            return tau0, tau1
+        return p_treatment * tau0 + p_control * tau1
+
+    def score(
+        self,
+        X: Union[pd.DataFrame, np.ndarray],
+        y: SliceDict,
+        sample_weight: Optional[np.ndarray] = None,
+        method: str = 'qini',
+        normalize: bool = True,
+        **kwargs,
+    ) -> float:
+        y_arr, treatment_ind = self._normalize_y(y)
+        if sample_weight is not None:
+            raise NotImplementedError
+        pred = self.predict(X=X)
+        if method == 'cumulative_gain':
+            return get_cumulative_gain_score(
+                y_true=y_arr,
+                treatment=treatment_ind,
+                score=pred,
+                **kwargs,
+            )
+        qini = get_qini_curve(
+            y_true=y_arr,
+            treatment=treatment_ind,
+            score=pred,
+            normalize=normalize,
+            **kwargs,
+        )
+        random_area = np.linspace(0, qini[-1], qini.shape[0]).sum()
+        return (np.nansum(qini) - random_area) / qini.shape[0]
+
+    @staticmethod
+    def _normalize_y(y: SliceDict) -> Tuple[np.ndarray, np.ndarray]:
+        y = y.copy()
+        y_arr = to_1d(np.asanyarray(y.pop('value')))
+        treatment_ind = to_1d(np.asanyarray(y.pop('is_treatment')).astype(bool))
+        if len(y.keys()):
+            warnings.warn(f"Unused keys in ``y``: {set(y)}")
+        return y_arr, treatment_ind
diff --git a/foundry/uplift/util.py b/foundry/uplift/util.py
new file mode 100644
index 0000000..6a6b980
--- /dev/null
+++ b/foundry/uplift/util.py
@@ -0,0 +1,255 @@
+import numpy as np
+import pandas as pd
+from sklearn.pipeline import Pipeline
+
+
+def get_qini_curve(y_true: np.ndarray,
+                   treatment: np.ndarray,
+                   score: np.ndarray,
+                   min_n_per: int = 1,
+                   normalize: bool = True) -> np.ndarray:
+    """
+    Adapted from https://www.uplift-modeling.com/en/latest/_modules/sklift/metrics/metrics.html#qini_curve
+
+    :param y_true: The true values
+    :param treatment: A treatment indicator (boolean).
+    :param score: The uplift score predicted for each record.
+    :param min_n_per: Minimum number of treatment and control records. For example, ``min_n_per=2`` means no qini
+     calculations until both (1) treatment has at least 2 records **and** (2) control has at least 2 records.
+    :param normalize: Whether to normalize to 0-1.
+    :return: np.ndarray with
+    """
+
+    y_true = np.asarray(y_true)
+    score = np.asarray(score)
+    treatment = np.asarray(treatment, dtype='bool')
+
+    desc_score_indices = np.argsort(score, kind="mergesort")[::-1]
+    y_true = y_true[desc_score_indices]
+    treatment = treatment[desc_score_indices]
+    uplift = score[desc_score_indices]
+
+    distinct_value_indices = np.where(np.diff(uplift))[0]
+    threshold_indices = np.concatenate([distinct_value_indices, [uplift.size - 1]])
+
+    cumu_num_trmnt = np.cumsum(treatment)[threshold_indices]
+    y_trmnt = np.cumsum(np.where(treatment, y_true, 0))[threshold_indices]
+
+    cumu_num_all = threshold_indices + 1
+    cumu_num_ctrl = cumu_num_all - cumu_num_trmnt
+    y_ctrl = np.cumsum(np.where(~treatment, y_true, 0))[threshold_indices]
+
+    mask = (cumu_num_trmnt >= min_n_per) & (cumu_num_ctrl >= min_n_per)
+    ratio = cumu_num_trmnt[mask] / cumu_num_ctrl[mask]
+    curve_values = y_trmnt[mask] - y_ctrl[mask] * ratio
+    # TODO
+    # if num_all.size == 0 or curve_values[0] != 0 or num_all[0] != 0:
+    #     # Add an extra threshold position if necessary
+    #     # to make sure that the curve starts at (0, 0)
+    #     curve_values = np.r_[0, curve_values]
+
+    out = np.full(len(cumu_num_all), np.nan)
+    out[mask] = curve_values
+    if normalize:
+        out /= out[-1]
+
+    return out
+
+
+def qini_scorer(estimator, X, y, normalize=True) -> float:
+    if isinstance(estimator, Pipeline):
+        X_transformed = estimator[:-1].transform(X)
+        return estimator[-1].score(X_transformed, y, method='qini', normalize=normalize)
+    return estimator.score(X, y, method='qini', normalize=normalize)
+
+
+# ── Elasticity helpers ────────────────────────────────────────────────────────
+# Adapted from https://matheusfacure.github.io/python-causality-handbook/21-Meta-Learners.html
+
+def elast(data, y: str, t: str) -> float:
+    """
+    OLS slope of y on t -- equivalent to the ATE under randomization.
+    Matches the @curry elast from the causality handbook.
+    """
+    t_vals = np.asarray(data[t], dtype=float)
+    y_vals = np.asarray(data[y], dtype=float)
+    num = np.sum((t_vals - t_vals.mean()) * (y_vals - y_vals.mean()))
+    den = np.sum((t_vals - t_vals.mean()) ** 2)
+    return float(num / den) if den != 0 else np.nan
+
+
+def elast_ci(data, y: str, t: str, z: float = 1.96) -> np.ndarray:
+    """95% confidence interval around the elasticity estimate."""
+    n = len(data)
+    t_bar  = np.asarray(data[t], dtype=float).mean()
+    beta1  = elast(data, y, t)
+    beta0  = np.asarray(data[y], dtype=float).mean() - beta1 * t_bar
+    e      = np.asarray(data[y], dtype=float) - (beta0 + beta1 * np.asarray(data[t], dtype=float))
+    se     = np.sqrt(((1 / (n - 2)) * np.sum(e ** 2)) /
+                     np.sum((np.asarray(data[t], dtype=float) - t_bar) ** 2))
+    return np.array([beta1 - z * se, beta1 + z * se])
+
+
+# ── Cumulative gain / elasticity curves ───────────────────────────────────────────
+# Adapted from https://matheusfacure.github.io/python-causality-handbook/21-Meta-Learners.html
+# Code here: https://github.com/matheusfacure/python-causality-handbook/blob/master/causal-inference-for-the-brave-and-true/nb21.py
+
+def cumulative_gain(dataset, prediction: str, y: str, t: str,
+                    min_periods: int = 30, steps: int = 100) -> np.ndarray:
+    """
+    Cumulative gain curve -- matches the causality handbook implementation.
+    Returns a 1-D array of length ~steps suitable for plt.plot().
+    """
+    size       = dataset.shape[0]
+    ordered_df = dataset.sort_values(prediction, ascending=False).reset_index(drop=True)
+    n_rows     = list(range(min_periods, size, size // steps)) + [size]
+    return np.array([0.0] + [elast(ordered_df.head(rows), y, t) * (rows / size)
+                             for rows in n_rows])
+
+
+def cumulative_gain_ci(dataset, prediction: str, y: str, t: str,
+                       min_periods: int = 30, steps: int = 100) -> np.ndarray:
+    """
+    Cumulative gain curve with 95% CI bands.
+    Returns an (N, 2) array of [lower, upper] at each step.
+    """
+    size       = dataset.shape[0]
+    ordered_df = dataset.sort_values(prediction, ascending=False).reset_index(drop=True)
+    n_rows     = list(range(min_periods, size, size // steps)) + [size]
+    return np.array([[0.0, 0.0]] + [elast_ci(ordered_df.head(rows), y, t) * (rows / size)
+                                    for rows in n_rows])
+
+
+def cumulative_elast_curve_ci(dataset, prediction: str, y: str, t: str,
+                               min_periods: int = 30, steps: int = 100) -> np.ndarray:
+    """
+    Cumulative elasticity curve with 95% CI (not multiplied by rows/size).
+    Returns an (N, 2) array of [lower, upper] at each step.
+    """
+    size       = dataset.shape[0]
+    ordered_df = dataset.sort_values(prediction, ascending=False).reset_index(drop=True)
+    n_rows     = list(range(min_periods, size, size // steps)) + [size]
+    return np.array([elast_ci(ordered_df.head(rows), y, t) for rows in n_rows])
+
+
+# ── Cumulative gain scalar score ──────────────────────────────────────────────
+
+def get_cumulative_gain_score(y_true: np.ndarray,
+                               treatment: np.ndarray,
+                               score: np.ndarray,
+                               steps: int = 100,
+                               min_periods: int = 30) -> float:
+    """
+    Compute the area between the cumulative gain curve and the random baseline
+    diagonal as a scalar model quality score.
+
+    The cumulative gain curve sorts observations by predicted CATE descending,
+    then at each population fraction k computes elast(top k%) * k -- the
+    expected gain from treating only the top k% of players. The random baseline
+    is a straight line from (0, 0) to (1, ATE). A model that ranks the most
+    treatment-responsive players first arcs above the baseline early, yielding
+    a positive score. A model with no ranking ability tracks the diagonal and
+    scores near zero.
+
+    :param y_true: Binary outcome array (0/1).
+    :param treatment: Binary treatment indicator (0/1 or bool).
+    :param score: Predicted CATE or uplift score for each observation. Higher
+        values are assumed to indicate higher predicted treatment responsiveness.
+    :param steps: Number of evaluation points along the population fraction
+        axis. More steps give a smoother curve and more precise AUC estimate.
+        Default is 100.
+    :param min_periods: Minimum number of observations required before
+        computing elasticity at a given population fraction. Avoids unstable
+        estimates at very small sample sizes. Default is 30.
+    :return: float -- area between the cumulative gain curve and the random
+        baseline. Positive = model beats random targeting. Returns np.nan if
+        fewer than 2 finite evaluation points exist.
+    """
+    n     = len(y_true)
+    order = np.argsort(score)[::-1]
+    y_s   = y_true[order]
+    t_s   = treatment[order]
+
+    n_rows = list(range(min_periods, n, max(1, n // steps))) + [n]
+
+    sorted_df  = pd.DataFrame({'y': y_s, 't': t_s})
+    full_df    = pd.DataFrame({'y': y_true, 't': treatment})
+
+    gains    = np.array([elast(sorted_df.head(k), 'y', 't') * (k / n) for k in n_rows])
+    xs       = np.array([k / n for k in n_rows])
+    baseline = elast(full_df, 'y', 't')
+    baseline_curve = xs * baseline
+
+    finite = np.isfinite(gains) & np.isfinite(baseline_curve)
+    if finite.sum() < 2:
+        return np.nan
+
+    return float(
+        np.trapezoid(gains[finite], xs[finite]) -
+        np.trapezoid(baseline_curve[finite], xs[finite])
+    )
+
+
+def cumulative_gain_scorer(estimator, X, y) -> float:
+    """
+    Scorer callable with signature ``(estimator, X, y)`` for use with
+    ``sklearn.model_selection.cross_validate`` or ``GridSearchCV``.
+
+    :param estimator: A fitted estimator with a ``score(X, y, method=...)`` method.
+    :param X: Feature matrix passed to ``estimator.score``.
+    :param y: SliceDict with keys ``'value'`` (outcome) and ``'is_treatment'``
+        (treatment indicator), passed to ``estimator.score``.
+    :return: float -- area between the cumulative gain curve and the random
+        baseline. Positive values indicate the model beats random targeting.
+
+    Example (cross_validate)
+    ------------------------
+    ::
+
+        from foundry.uplift.util import cumulative_gain_scorer
+        from sklearn.model_selection import StratifiedKFold, cross_validate
+
+        skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
+        splits = list(skf.split(df, t))
+
+        scores_dict = cross_validate(
+            model,
+            X=df,
+            y=SliceDict(value=df['outcome'], is_treatment=df['treatment']),
+            cv=splits,
+            scoring=cumulative_gain_scorer,
+            return_train_score=True,
+        )
+
+        print(scores_dict['test_score'])
+        print(scores_dict['train_score'])
+
+    Example (GridSearchCV)
+    ----------------------
+    ::
+
+        from foundry.uplift.util import cumulative_gain_scorer
+        from sklearn.model_selection import GridSearchCV
+
+        gs = GridSearchCV(
+            estimator=model,
+            param_grid={'tlearner__estimator__max_depth': [2, 3, 4]},
+            scoring=cumulative_gain_scorer,
+            cv=splits,
+            refit=True,
+        )
+
+        gs.fit(
+            df,
+            SliceDict(value=df['outcome'], is_treatment=df['treatment']),
+        )
+
+        print(f"Best params: {gs.best_params_}")
+        print(f"Best score:  {gs.best_score_:.4f}")
+    """
+    # unwrap pipeline if needed to reach the learner's score method
+    if isinstance(estimator, Pipeline):
+        # transform X through all steps except the last
+        X_transformed = estimator[:-1].transform(X)
+        return estimator[-1].score(X_transformed, y, method='cumulative_gain')
+    return estimator.score(X, y, method='cumulative_gain')
diff --git a/foundry/util.py b/foundry/util.py
index bbd9abb..9a90b7a 100644
--- a/foundry/util.py
+++ b/foundry/util.py
@@ -116,6 +116,22 @@ class FitFailedException(RuntimeError):
     pass
 
 
+def safe_predict(estimator, *args, **kwargs) -> np.ndarray:
+    if hasattr(estimator, 'predict_proba'):
+        try:
+            out = estimator.predict_proba(*args, **kwargs)
+        except NotImplementedError:
+            out = None
+
+        if out is not None:
+            if out.shape[1] == 2:
+                out = out[:, 1]
+            elif out.shape[1] > 2:
+                raise NotImplementedError("Multi-class predict_proba not supported.")
+            return out
+    return estimator.predict(*args, **kwargs)
+
+
 class SliceDict(dict):
     """
     Adapted from https://github.com/skorch-dev/skorch/blob/baf0580/skorch/helper.py#L20
@@ -135,8 +151,23 @@ def __init__(self, **kwargs):
         else:
             self._len = lengths[0]
 
+        # sklearn checks if it should use pandas indexing by checking if there's an iloc attribute
+        if self.is_pandas:
+            self.__dict__['iloc'] = True
+        else:
+            self.__dict__.pop('iloc', None)
+
         super().__init__(**kwargs)
 
+    @property
+    def is_pandas(self) -> bool:
+        is_pandas = [hasattr(v, 'iloc') for v in self.values() if hasattr(v, 'shape')]
+        any_pandas = any(is_pandas)
+        all_pandas = all(is_pandas)
+        if any_pandas and not all_pandas:
+            raise ValueError("Currenlty SliceDict does not support a mix of pandas and non-pandas")
+        return any_pandas
+
     @staticmethod
     def _standardize_val(val):
         return np.asarray(val) if is_array(val) and not hasattr(val, 'shape') else val
@@ -151,7 +182,8 @@ def __getitem__(self, sl: Union[int, str, slice]) -> Union['SliceDict', ArrayTyp
             )
         if isinstance(sl, str):
             return super(SliceDict, self).__getitem__(sl)
-        return SliceDict(**{k: (v[sl] if hasattr(v, 'shape') else v) for k, v in self.items()})
+        cls = type(self)
+        return cls(**{k: (v[sl] if hasattr(v, 'shape') else v) for k, v in self.items()})
 
     def __setitem__(self, key: str, value: ArrayType):
         value = self._standardize_val(value)
@@ -177,6 +209,18 @@ def __setitem__(self, key: str, value: ArrayType):
 
         super().__setitem__(key, value)
 
+        # sklearn checks if it should use pandas indexing by checking if there's an iloc attribute
+        if self.is_pandas:
+            self.__dict__['iloc'] = True
+        else:
+            self.__dict__.pop('iloc', None)
+
+    def take(self, indices, axis: int = 0, **kwargs) -> 'SliceDict':
+        if axis:
+            raise ValueError("Only axis=0 is supported")
+        cls = type(self)
+        return cls(**{k: (v.take(indices, axis, **kwargs) if hasattr(v, 'shape') else v) for k, v in self.items()})
+
     def update(self, kwargs: dict):
         for key, value in kwargs.items():
             self.__setitem__(key, value)