From 4e12fb77272d3e2bd11b805185832d1409518ba8 Mon Sep 17 00:00:00 2001
From: daveringelberg <141670222+daveringelberg@users.noreply.github.com>
Date: Wed, 9 Sep 2026 12:50:07 +0200
Subject: [PATCH 1/8] Add native Hill dose-response fitting

---
 docs/api/tools_index.md                       |  53 +++++-
 .../_perturbation_space.py                    | 153 ++++++++++++++++++
 .../test_perturbation_space_extras.py         | 146 ++++++++++++++++-
 3 files changed, 348 insertions(+), 4 deletions(-)

diff --git a/docs/api/tools_index.md b/docs/api/tools_index.md
index 59aa1a00..f3372f31 100644
--- a/docs/api/tools_index.md
+++ b/docs/api/tools_index.md
@@ -564,4 +564,55 @@ ds_adata = ds.compute(mdata["rna"], target_col="gene_target", metric="edistance"
 similar = ds.nearest_perturbations(ds_adata, "IFNGR2", target_col="gene_target")
 ```
 
-See [perturbation space tutorial](https://pertpy.readthedocs.io/en/latest/tutorials/notebooks/perturbation_space.html).
+### Dose-response curve fitting
+
+`PerturbationSpace.dose_response` calculates a scalar distance from control for every perturbation and dose.
+`PerturbationSpace.fit_dose_response` fits a four-parameter Hill curve to these values and reports the EC50 for this transcriptomic-distance response:
+
+```python
+import pertpy as pt
+import scanpy as sc
+
+adata = pt.dt.srivatsan_2020_sciplex2()
+adata = adata[adata.obs["dose_value"].notna()].copy()
+adata.obs["dose_value"] = adata.obs["dose_value"].astype(float)
+adata.obs["perturbation"] = adata.obs["perturbation"].astype(str)
+adata.obs.loc[adata.obs["dose_value"] == 0, "perturbation"] = "zero_dose"
+sc.pp.normalize_total(adata, target_sum=1e4)
+sc.pp.log1p(adata)
+sc.pp.pca(adata)
+ps = pt.tl.PseudobulkSpace()
+
+responses = ps.dose_response(
+    adata,
+    dose_col="dose_value",
+    reference_key="zero_dose",
+    embedding_key="X_pca",
+)
+fits = ps.fit_dose_response(responses)
+```
+
+This example pools the assigned zero-dose samples as the reference and excludes cells without a dose assignment.
+The existing `control` label in this dataset includes [cells without an assigned sample](https://github.com/sanderlab/scPerturb/blob/master/dataset_processing/scripts/SrivatsanTrapnell2020.py#L55-L56); it is not used as the reference here.
+Doses retain the dataset's recorded units; verify their mapping to physical concentrations before comparing potency between compounds.
+
+The same fitter accepts other scalar assay responses.
+For example, an independently measured and control-normalized viability response can be labelled as inhibition to obtain IC50 rather than EC50 output:
+
+```python
+fits = ps.fit_dose_response(
+    assay_responses,
+    response_col="viability",
+    response_type="inhibition",
+)
+```
+
+`response_type` only determines how the fitted midpoint is interpreted and named; the fitter does not perform biological or control normalization.
+Internal numerical rescaling does not change the output response units.
+EC50 and IC50 are relative midpoints between the fitted `e0` and `emax`.
+Check `r_squared`, the reported standard error and `midpoint_in_range` before interpreting a fit; a midpoint outside the measured positive-dose range is an extrapolation.
+If the uncertainty cannot be estimated numerically, a warning is emitted and the standard error is `NaN`; the fitted parameters remain available for inspection.
+The standard error is a local approximation. Even a small error, high R-squared and an in-range midpoint can be misleading when the doses do not constrain both plateaus; inspect the measured responses and fitted curve.
+The terminology assumes that `dose` represents a concentration, and this first implementation fits a four-parameter curve without automatically selecting fixed-top or fixed-bottom alternatives.
+
+See the [perturbation space tutorial](https://pertpy.readthedocs.io/en/latest/tutorials/notebooks/perturbation_space.html#dose-response) for plotting the measured dose responses and fitted Hill curves together.
diff --git a/src/pertpy/tools/_perturbation_space/_perturbation_space.py b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
index c82d7799..002e3d72 100644
--- a/src/pertpy/tools/_perturbation_space/_perturbation_space.py
+++ b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
@@ -7,6 +7,8 @@
 import numpy as np
 import pandas as pd
 from anndata import AnnData
+from scipy.optimize import curve_fit
+from scipy.special import expit
 from scipy.stats import entropy
 
 from pertpy._logger import logger
@@ -83,6 +85,16 @@ def _vector_distance(u: np.ndarray, v: np.ndarray, metric: str) -> float:
     raise ValueError(f"Unknown metric {metric!r}. Choose from 'euclidean', 'cosine', 'pearson'.")
 
 
+def _four_parameter_logistic(
+    dose: np.ndarray, e0: float, emax: float, log_midpoint: float, hill_coefficient: float
+) -> np.ndarray:
+    """Evaluate a four-parameter Hill curve with a positive Hill coefficient."""
+    log_dose = np.full_like(dose, -np.inf, dtype=float)
+    np.log(dose, out=log_dose, where=dose > 0)
+    fraction = expit(hill_coefficient * (log_dose - log_midpoint))
+    return e0 + (emax - e0) * fraction
+
+
 def _subtract_control_mean(
     matrix: np.ndarray,
     control_mask: np.ndarray,
@@ -648,6 +660,147 @@ def dose_response(
                 result["dose"] = pd.to_numeric(result["dose"])
         return result.sort_values(["perturbation", "dose"]).reset_index(drop=True)
 
+    def fit_dose_response(
+        self,
+        data: pd.DataFrame,
+        *,
+        perturbation_col: str = "perturbation",
+        dose_col: str = "dose",
+        response_col: str = "distance",
+        response_type: Literal["effect", "inhibition"] = "effect",
+    ) -> pd.DataFrame:
+        """Fit a four-parameter Hill curve for each perturbation.
+
+        ``data`` can be the output of :meth:`dose_response` or a table containing another scalar assay response.
+        ``response_type`` names the fitted midpoint according to the meaning of that response: ``"effect"`` returns
+        ``ec50``, while ``"inhibition"`` returns ``ic50``. It does not perform biological or control normalization.
+        Dose values must represent concentrations for the EC50 or IC50 terminology to apply.
+
+        Args:
+            data: Tidy table containing perturbation, dose and response columns.
+            perturbation_col: Column identifying the perturbation.
+            dose_col: Column containing non-negative numeric doses.
+            response_col: Column containing the scalar response to fit.
+            response_type: Whether the response represents an effect or inhibition.
+
+        Returns:
+            One row per perturbation with the fitted zero-dose response (``e0``), asymptotic response (``emax``),
+            Hill coefficient, EC50 or IC50, its approximate standard error, R-squared and whether the midpoint lies
+            within the tested positive-dose range.
+
+        Notes:
+            Responses are rescaled internally for numerical stability; ``e0`` and ``emax`` retain the input units.
+            The standard error uses a local linear approximation. It is NaN, with a warning, if the parameter
+            covariance is non-finite or numerically rank deficient, or there are no residual degrees of freedom.
+            Parameters are retained for inspection. A small standard error, high R-squared or an in-range midpoint
+            does not establish that the doses capture both plateaus or that the Hill model is appropriate.
+
+        Examples:
+            >>> import pertpy as pt
+            >>> import scanpy as sc
+            >>> adata = pt.dt.srivatsan_2020_sciplex2()
+            >>> adata = adata[adata.obs["dose_value"].notna()].copy()
+            >>> adata.obs["dose_value"] = adata.obs["dose_value"].astype(float)
+            >>> adata.obs["perturbation"] = adata.obs["perturbation"].astype(str)
+            >>> adata.obs.loc[adata.obs["dose_value"] == 0, "perturbation"] = "zero_dose"
+            >>> sc.pp.normalize_total(adata, target_sum=1e4)
+            >>> sc.pp.log1p(adata)
+            >>> sc.pp.pca(adata)
+            >>> ps = pt.tl.PseudobulkSpace()
+            >>> responses = ps.dose_response(
+            ...     adata, dose_col="dose_value", reference_key="zero_dose", embedding_key="X_pca"
+            ... )
+            >>> fits = ps.fit_dose_response(responses)
+        """
+        required = {perturbation_col, dose_col, response_col}
+        missing = required.difference(data.columns)
+        if missing:
+            raise ValueError(f"Columns {sorted(missing)} do not exist in the input data.")
+        if response_type not in {"effect", "inhibition"}:
+            raise ValueError("response_type must be either 'effect' or 'inhibition'.")
+
+        fit_data = data[[perturbation_col, dose_col, response_col]].copy()
+        if fit_data[perturbation_col].isna().any():
+            raise ValueError("Perturbation labels must not be missing.")
+        fit_data[dose_col] = pd.to_numeric(fit_data[dose_col], errors="raise")
+        fit_data[response_col] = pd.to_numeric(fit_data[response_col], errors="raise")
+        if (fit_data[dose_col] < 0).any():
+            raise ValueError("Dose values must be non-negative.")
+
+        midpoint_col = "ec50" if response_type == "effect" else "ic50"
+        records: list[dict[str, object]] = []
+        for perturbation, group in fit_data.groupby(perturbation_col, observed=True, sort=True):
+            doses = group[dose_col].to_numpy(dtype=float)
+            responses = group[response_col].to_numpy(dtype=float)
+            if np.unique(doses).size < 4:
+                raise ValueError(f"Perturbation {perturbation!r} needs at least four distinct dose values.")
+            response_offset = float(np.min(responses))
+            response_scale = float(np.ptp(responses))
+            if response_scale == 0:
+                raise ValueError(
+                    f"Perturbation {perturbation!r} has a constant response, so a Hill curve cannot be fit."
+                )
+
+            responses = (responses - response_offset) / response_scale
+            mean_response = (
+                group.groupby(dose_col, sort=True, observed=True)[response_col].mean() - response_offset
+            ) / response_scale
+            e0_guess = float(mean_response.iloc[0])
+            emax_guess = float(mean_response.iloc[-1])
+            halfway = (e0_guess + emax_guess) / 2
+            positive_response = mean_response[mean_response.index > 0]
+            midpoint_guess = float((positive_response - halfway).abs().idxmin())
+
+            parameters, covariance = curve_fit(
+                _four_parameter_logistic,
+                doses,
+                responses,
+                p0=(e0_guess, emax_guess, np.log(midpoint_guess), 1.0),
+                bounds=((-np.inf, -np.inf, -np.inf, np.finfo(float).eps), np.inf),
+                absolute_sigma=True,
+                maxfev=20_000,
+            )
+
+            e0, emax, log_midpoint, hill_coefficient = parameters
+            midpoint = float(np.exp(log_midpoint))
+            fitted = _four_parameter_logistic(doses, *parameters)
+            residual_sum_squares = float(np.sum((responses - fitted) ** 2))
+            total_sum_squares = float(np.sum((responses - responses.mean()) ** 2))
+            positive_doses = doses[doses > 0]
+
+            # Check conditioning before scaling covariance by residual variance, which may be zero.
+            degrees_of_freedom = len(responses) - len(parameters)
+            if (
+                degrees_of_freedom <= 0
+                or not np.isfinite(covariance).all()
+                or np.linalg.matrix_rank(covariance) < len(parameters)
+            ):
+                midpoint_standard_error = np.nan
+                warnings.warn(
+                    f"Cannot estimate the {midpoint_col.upper()} standard error for perturbation {perturbation!r}. "
+                    "Inspect the dose range and fitted curve before interpreting the estimate.",
+                    UserWarning,
+                    stacklevel=2,
+                )
+            else:
+                residual_variance = residual_sum_squares / degrees_of_freedom
+                midpoint_standard_error = midpoint * np.sqrt(float(covariance[2, 2]) * residual_variance)
+
+            records.append(
+                {
+                    perturbation_col: perturbation,
+                    "e0": float(e0 * response_scale + response_offset),
+                    "emax": float(emax * response_scale + response_offset),
+                    "hill_coefficient": float(hill_coefficient),
+                    midpoint_col: midpoint,
+                    f"{midpoint_col}_standard_error": float(midpoint_standard_error),
+                    "r_squared": 1 - residual_sum_squares / total_sum_squares,
+                    "midpoint_in_range": bool(positive_doses.min() < midpoint < positive_doses.max()),
+                }
+            )
+
+        return pd.DataFrame.from_records(records)
+
     def plot_similarity(  # pragma: no cover
         self,
         adata: AnnData,
diff --git a/tests/tools/_perturbation_space/test_perturbation_space_extras.py b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
index 3d1501e8..f85f99c5 100644
--- a/tests/tools/_perturbation_space/test_perturbation_space_extras.py
+++ b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
@@ -66,19 +66,159 @@ def test_evaluate_combinations(rng):
     np.testing.assert_allclose(result.loc["A+B", "distance"], 0.0, atol=1e-6)
 
 
-def test_dose_response(rng):
+def test_dose_response():
+    rng = np.random.default_rng(0)
     groups, doses = [], []
     for pert in ["control", "drug"]:
-        for dose in [0.0] if pert == "control" else [1.0, 10.0, 100.0]:
+        for dose in [0.0] if pert == "control" else [0.1, 1.0, 3.0, 10.0, 30.0, 100.0]:
             groups += [pert] * 15
             doses += [dose] * 15
     groups = np.array(groups)
     doses = np.array(doses, dtype=float)
     X = rng.normal(0, 0.3, (len(groups), 8))
-    X[groups == "drug"] += (doses[groups == "drug"] / 10.0)[:, None]
+    drug_doses = doses[groups == "drug"]
+    X[groups == "drug"] += (5 * drug_doses**1.2 / (10**1.2 + drug_doses**1.2))[:, None]
     adata = AnnData(X, obs=pd.DataFrame({"perturbation": groups, "dose": doses}))
     sc.pp.pca(adata, n_comps=5)
 
     curves = pt.tl.PseudobulkSpace().dose_response(adata, dose_col="dose", metric="euclidean", embedding_key="X_pca")
     drug = curves[curves["perturbation"] == "drug"].sort_values("dose")
     assert drug["distance"].is_monotonic_increasing
+
+    fits = pt.tl.PseudobulkSpace().fit_dose_response(curves)
+    assert fits.loc[0, "ec50"] == pytest.approx(10, rel=0.2)
+    assert fits.loc[0, "r_squared"] > 0.99
+    assert fits.loc[0, "midpoint_in_range"]
+
+
+@pytest.mark.parametrize(
+    ("response_type", "e0", "emax", "midpoint", "midpoint_col"),
+    [("effect", 0.1, 1.8, 3.0, "ec50"), ("inhibition", 1.0, 0.05, 8.0, "ic50")],
+)
+@pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
+def test_fit_dose_response(*, response_type, e0, emax, midpoint, midpoint_col, response_scale):
+    e0, emax = e0 * response_scale, emax * response_scale
+    doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0])
+    hill_coefficient = 1.4
+    fraction = doses**hill_coefficient / (midpoint**hill_coefficient + doses**hill_coefficient)
+    first = pd.DataFrame(
+        {
+            "compound": "drug_a",
+            "concentration": doses,
+            "response": e0 + (emax - e0) * fraction,
+        }
+    )
+    second_midpoint = midpoint * 2
+    second_fraction = doses**hill_coefficient / (second_midpoint**hill_coefficient + doses**hill_coefficient)
+    second = first.assign(compound="drug_b", response=e0 + (emax - e0) * second_fraction)
+    data = pd.concat([first, second], ignore_index=True)
+
+    fits = (
+        pt.tl.PseudobulkSpace()
+        .fit_dose_response(
+            data,
+            perturbation_col="compound",
+            dose_col="concentration",
+            response_col="response",
+            response_type=response_type,
+        )
+        .set_index("compound")
+    )
+
+    assert midpoint_col in fits
+    assert f"{midpoint_col}_standard_error" in fits
+    assert {"e0", "emax", "hill_coefficient", "r_squared", "midpoint_in_range"} <= set(fits)
+    assert fits.loc["drug_a", "e0"] == pytest.approx(e0, rel=1e-6, abs=0)
+    assert fits.loc["drug_a", "emax"] == pytest.approx(emax, rel=1e-6, abs=0)
+    assert fits.loc["drug_a", "hill_coefficient"] == pytest.approx(hill_coefficient)
+    assert fits.loc["drug_a", midpoint_col] == pytest.approx(midpoint)
+    assert fits.loc["drug_b", midpoint_col] == pytest.approx(second_midpoint)
+    assert fits.loc["drug_a", "r_squared"] == pytest.approx(1)
+    assert np.isfinite(fits[f"{midpoint_col}_standard_error"]).all()
+
+
+@pytest.mark.parametrize("response_scale", [1.0, 100.0])
+def test_fit_dose_response_rank_deficient(response_scale):
+    doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
+    responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
+    complete = pd.DataFrame({"perturbation": "complete", "dose": doses, "distance": response_scale * responses})
+    limited = complete.loc[complete["dose"] <= 3].assign(perturbation="limited")
+
+    with pytest.warns(UserWarning, match="Cannot estimate.*'limited'"):
+        fits = pt.tl.PseudobulkSpace().fit_dose_response(pd.concat([complete, limited])).set_index("perturbation")
+
+    assert fits.loc["complete", "ec50"] == pytest.approx(10, rel=0.1)
+    assert np.isfinite(fits.loc["complete", "ec50_standard_error"])
+    # A high R-squared and an in-range midpoint do not expose this poorly determined curve.
+    assert fits.loc["limited", "r_squared"] > 0.98
+    assert fits.loc["limited", "midpoint_in_range"]
+    assert np.isfinite(fits.loc["limited", "ec50"])
+    assert np.isnan(fits.loc["limited", "ec50_standard_error"])
+
+
+@pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
+def test_fit_dose_response_standard_error(response_scale):
+    doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
+    responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
+    responses *= response_scale
+    data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": responses})
+    fit = pt.tl.PseudobulkSpace().fit_dose_response(data).iloc[0]
+
+    # Independent derivatives with respect to EC50 itself, rather than the fitted log(EC50).
+    midpoint, slope = fit["ec50"], fit["hill_coefficient"]
+    fraction = doses**slope / (midpoint**slope + doses**slope)
+    log_ratio = np.zeros_like(doses)
+    np.log(doses / midpoint, out=log_ratio, where=doses > 0)
+    sensitivity = (fit["emax"] - fit["e0"]) * fraction * (1 - fraction)
+    jacobian = np.column_stack((1 - fraction, fraction, -slope * sensitivity / midpoint, sensitivity * log_ratio))
+    residuals = responses - (fit["e0"] + (fit["emax"] - fit["e0"]) * fraction)
+    variance = np.sum(residuals**2) / (len(doses) - 4)
+    expected_error = np.sqrt(np.linalg.inv(jacobian.T @ jacobian)[2, 2] * variance)
+    assert fit["ec50_standard_error"] == pytest.approx(expected_error, rel=1e-4)
+
+
+def test_fit_dose_response_insufficient_dof():
+    doses = np.array([0.0, 1.0, 10.0, 100.0])
+    data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": doses / (10 + doses)})
+    with pytest.warns(UserWarning, match="Cannot estimate.*'drug'"):
+        fit = pt.tl.PseudobulkSpace().fit_dose_response(data).iloc[0]
+    assert fit["ec50"] == pytest.approx(10)
+    assert np.isnan(fit["ec50_standard_error"])
+
+
+@pytest.mark.parametrize(
+    ("data", "kwargs", "match"),
+    [
+        (pd.DataFrame(), {}, "Columns"),
+        (
+            pd.DataFrame(
+                {"perturbation": ["drug", "drug", None, "drug"], "dose": [0, 1, 2, 3], "distance": [0, 1, 2, 3]}
+            ),
+            {},
+            "Perturbation labels",
+        ),
+        (
+            pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [-1, 1, 2, 3], "distance": [0, 1, 2, 3]}),
+            {},
+            "non-negative",
+        ),
+        (
+            pd.DataFrame({"perturbation": ["drug"] * 3, "dose": [0, 1, 2], "distance": [0, 1, 2]}),
+            {},
+            "four distinct",
+        ),
+        (
+            pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [0, 1, 2, 3], "distance": [1, 1, 1, 1]}),
+            {},
+            "constant response",
+        ),
+        (
+            pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [0, 1, 2, 3], "distance": [0, 1, 2, 3]}),
+            {"response_type": "unknown"},
+            "response_type",
+        ),
+    ],
+)
+def test_fit_dose_response_validation(data, kwargs, match):
+    with pytest.raises(ValueError, match=match):
+        pt.tl.PseudobulkSpace().fit_dose_response(data, **kwargs)

From a583b9265c6d853f3ea30c505a90f77d1fe099de Mon Sep 17 00:00:00 2001
From: daveringelberg <141670222+daveringelberg@users.noreply.github.com>
Date: Sun, 13 Sep 2026 08:46:46 +0200
Subject: [PATCH 2/8] Address Hill fitting API and documentation review

---
 docs/api/tools_index.md                       | 52 ++---------
 .../_perturbation_space.py                    | 59 ++++++++-----
 .../test_perturbation_space_extras.py         | 87 ++++++++++++++-----
 3 files changed, 108 insertions(+), 90 deletions(-)

diff --git a/docs/api/tools_index.md b/docs/api/tools_index.md
index f3372f31..07cf3262 100644
--- a/docs/api/tools_index.md
+++ b/docs/api/tools_index.md
@@ -564,55 +564,21 @@ ds_adata = ds.compute(mdata["rna"], target_col="gene_target", metric="edistance"
 similar = ds.nearest_perturbations(ds_adata, "IFNGR2", target_col="gene_target")
 ```
 
+See [perturbation space tutorial](https://pertpy.readthedocs.io/en/latest/tutorials/notebooks/perturbation_space.html).
+
 ### Dose-response curve fitting
 
-`PerturbationSpace.dose_response` calculates a scalar distance from control for every perturbation and dose.
-`PerturbationSpace.fit_dose_response` fits a four-parameter Hill curve to these values and reports the EC50 for this transcriptomic-distance response:
+{meth}`~pertpy.tools.PseudobulkSpace.dose_response` calculates a scalar distance from control for each perturbation and dose.
+{meth}`~pertpy.tools.PseudobulkSpace.fit_dose_response` fits a four-parameter Hill curve and stores the results in AnnData.
+For preprocessed `adata` with `perturbation` and `dose` columns, a `control` group and a PCA representation:
 
 ```python
 import pertpy as pt
-import scanpy as sc
 
-adata = pt.dt.srivatsan_2020_sciplex2()
-adata = adata[adata.obs["dose_value"].notna()].copy()
-adata.obs["dose_value"] = adata.obs["dose_value"].astype(float)
-adata.obs["perturbation"] = adata.obs["perturbation"].astype(str)
-adata.obs.loc[adata.obs["dose_value"] == 0, "perturbation"] = "zero_dose"
-sc.pp.normalize_total(adata, target_sum=1e4)
-sc.pp.log1p(adata)
-sc.pp.pca(adata)
 ps = pt.tl.PseudobulkSpace()
-
-responses = ps.dose_response(
-    adata,
-    dose_col="dose_value",
-    reference_key="zero_dose",
-    embedding_key="X_pca",
-)
-fits = ps.fit_dose_response(responses)
+responses = ps.dose_response(adata, embedding_key="X_pca")
+ps.fit_dose_response(adata, responses)
+fits = adata.uns["dose_response"]["fits"]
 ```
 
-This example pools the assigned zero-dose samples as the reference and excludes cells without a dose assignment.
-The existing `control` label in this dataset includes [cells without an assigned sample](https://github.com/sanderlab/scPerturb/blob/master/dataset_processing/scripts/SrivatsanTrapnell2020.py#L55-L56); it is not used as the reference here.
-Doses retain the dataset's recorded units; verify their mapping to physical concentrations before comparing potency between compounds.
-
-The same fitter accepts other scalar assay responses.
-For example, an independently measured and control-normalized viability response can be labelled as inhibition to obtain IC50 rather than EC50 output:
-
-```python
-fits = ps.fit_dose_response(
-    assay_responses,
-    response_col="viability",
-    response_type="inhibition",
-)
-```
-
-`response_type` only determines how the fitted midpoint is interpreted and named; the fitter does not perform biological or control normalization.
-Internal numerical rescaling does not change the output response units.
-EC50 and IC50 are relative midpoints between the fitted `e0` and `emax`.
-Check `r_squared`, the reported standard error and `midpoint_in_range` before interpreting a fit; a midpoint outside the measured positive-dose range is an extrapolation.
-If the uncertainty cannot be estimated numerically, a warning is emitted and the standard error is `NaN`; the fitted parameters remain available for inspection.
-The standard error is a local approximation. Even a small error, high R-squared and an in-range midpoint can be misleading when the doses do not constrain both plateaus; inspect the measured responses and fitted curve.
-The terminology assumes that `dose` represents a concentration, and this first implementation fits a four-parameter curve without automatically selecting fixed-top or fixed-bottom alternatives.
-
-See the [perturbation space tutorial](https://pertpy.readthedocs.io/en/latest/tutorials/notebooks/perturbation_space.html#dose-response) for plotting the measured dose responses and fitted Hill curves together.
+For an inhibitory assay response, use `response_col` to select the measurement column and `response_type="inhibition"` to report IC50 instead of EC50.
diff --git a/src/pertpy/tools/_perturbation_space/_perturbation_space.py b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
index 002e3d72..1a3405f8 100644
--- a/src/pertpy/tools/_perturbation_space/_perturbation_space.py
+++ b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
@@ -88,7 +88,7 @@ def _vector_distance(u: np.ndarray, v: np.ndarray, metric: str) -> float:
 def _four_parameter_logistic(
     dose: np.ndarray, e0: float, emax: float, log_midpoint: float, hill_coefficient: float
 ) -> np.ndarray:
-    """Evaluate a four-parameter Hill curve with a positive Hill coefficient."""
+    """Calculate predicted responses at the given doses from four-parameter Hill curve parameters."""
     log_dose = np.full_like(dose, -np.inf, dtype=float)
     np.log(dose, out=log_dose, where=dose > 0)
     fraction = expit(hill_coefficient * (log_dose - log_midpoint))
@@ -662,55 +662,60 @@ def dose_response(
 
     def fit_dose_response(
         self,
+        adata: AnnData,
         data: pd.DataFrame,
         *,
         perturbation_col: str = "perturbation",
         dose_col: str = "dose",
         response_col: str = "distance",
         response_type: Literal["effect", "inhibition"] = "effect",
-    ) -> pd.DataFrame:
+        key_added: str = "dose_response",
+    ) -> None:
         """Fit a four-parameter Hill curve for each perturbation.
 
         ``data`` can be the output of :meth:`dose_response` or a table containing another scalar assay response.
-        ``response_type`` names the fitted midpoint according to the meaning of that response: ``"effect"`` returns
-        ``ec50``, while ``"inhibition"`` returns ``ic50``. It does not perform biological or control normalization.
+        ``response_type`` names the fitted midpoint according to the meaning of that response: ``"effect"`` reports
+        ``ec50``, while ``"inhibition"`` reports ``ic50``.
+        It does not perform biological or control normalization.
         Dose values must represent concentrations for the EC50 or IC50 terminology to apply.
 
         Args:
+            adata: AnnData to store the fit results in.
             data: Tidy table containing perturbation, dose and response columns.
             perturbation_col: Column identifying the perturbation.
             dose_col: Column containing non-negative numeric doses.
             response_col: Column containing the scalar response to fit.
             response_type: Whether the response represents an effect or inhibition.
+            key_added: Key in `.uns` for the fit results and parameters.
 
         Returns:
-            One row per perturbation with the fitted zero-dose response (``e0``), asymptotic response (``emax``),
-            Hill coefficient, EC50 or IC50, its approximate standard error, R-squared and whether the midpoint lies
-            within the tested positive-dose range.
+            Updates `.uns[key_added]` with a ``fits`` table and a ``params`` dictionary.
+            The table contains one row per perturbation with the fitted zero-dose response (``e0``), asymptotic
+            response (``emax``), Hill coefficient, EC50 or IC50, its approximate standard error, R-squared and
+            whether the midpoint lies within the tested positive-dose range.
 
         Notes:
+            Fits are stored in `.uns` because they describe groups of observations, not individual cells or genes.
+            EC50 and IC50 are relative midpoints between the fitted ``e0`` and ``emax``.
             Responses are rescaled internally for numerical stability; ``e0`` and ``emax`` retain the input units.
             The standard error uses a local linear approximation. It is NaN, with a warning, if the parameter
             covariance is non-finite or numerically rank deficient, or there are no residual degrees of freedom.
-            Parameters are retained for inspection. A small standard error, high R-squared or an in-range midpoint
+            Parameters are retained for inspection.
+            A small standard error, high R-squared or an in-range midpoint
             does not establish that the doses capture both plateaus or that the Hill model is appropriate.
 
         Examples:
+            A minimal example with a simulated response:
+
             >>> import pertpy as pt
-            >>> import scanpy as sc
-            >>> adata = pt.dt.srivatsan_2020_sciplex2()
-            >>> adata = adata[adata.obs["dose_value"].notna()].copy()
-            >>> adata.obs["dose_value"] = adata.obs["dose_value"].astype(float)
-            >>> adata.obs["perturbation"] = adata.obs["perturbation"].astype(str)
-            >>> adata.obs.loc[adata.obs["dose_value"] == 0, "perturbation"] = "zero_dose"
-            >>> sc.pp.normalize_total(adata, target_sum=1e4)
-            >>> sc.pp.log1p(adata)
-            >>> sc.pp.pca(adata)
-            >>> ps = pt.tl.PseudobulkSpace()
-            >>> responses = ps.dose_response(
-            ...     adata, dose_col="dose_value", reference_key="zero_dose", embedding_key="X_pca"
-            ... )
-            >>> fits = ps.fit_dose_response(responses)
+            >>> import pandas as pd
+            >>> from anndata import AnnData
+            >>> adata = AnnData(obs=pd.DataFrame(index=["drug"]))
+            >>> doses = pd.Series([0, 1, 3, 10, 30, 100])
+            >>> responses = pd.DataFrame({"dose": doses, "distance": doses / (10 + doses)})
+            >>> responses["perturbation"] = "drug"
+            >>> pt.tl.PseudobulkSpace().fit_dose_response(adata, responses)
+            >>> fits = adata.uns["dose_response"]["fits"]
         """
         required = {perturbation_col, dose_col, response_col}
         missing = required.difference(data.columns)
@@ -799,7 +804,15 @@ def fit_dose_response(
                 }
             )
 
-        return pd.DataFrame.from_records(records)
+        adata.uns[key_added] = {
+            "fits": pd.DataFrame.from_records(records),
+            "params": {
+                "perturbation_col": perturbation_col,
+                "dose_col": dose_col,
+                "response_col": response_col,
+                "response_type": response_type,
+            },
+        }
 
     def plot_similarity(  # pragma: no cover
         self,
diff --git a/tests/tools/_perturbation_space/test_perturbation_space_extras.py b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
index f85f99c5..e26be0f0 100644
--- a/tests/tools/_perturbation_space/test_perturbation_space_extras.py
+++ b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
@@ -66,8 +66,8 @@ def test_evaluate_combinations(rng):
     np.testing.assert_allclose(result.loc["A+B", "distance"], 0.0, atol=1e-6)
 
 
-def test_dose_response():
-    rng = np.random.default_rng(0)
+@pytest.mark.parametrize("categorical_doses", [False, True])
+def test_dose_response(rng, categorical_doses):
     groups, doses = [], []
     for pert in ["control", "drug"]:
         for dose in [0.0] if pert == "control" else [0.1, 1.0, 3.0, 10.0, 30.0, 100.0]:
@@ -75,18 +75,22 @@ def test_dose_response():
             doses += [dose] * 15
     groups = np.array(groups)
     doses = np.array(doses, dtype=float)
-    X = rng.normal(0, 0.3, (len(groups), 8))
+    # Match background cells across doses so sampling noise does not change the known curve.
+    X = np.tile(rng.normal(0, 0.3, (15, 8)), (len(groups) // 15, 1))
     drug_doses = doses[groups == "drug"]
     X[groups == "drug"] += (5 * drug_doses**1.2 / (10**1.2 + drug_doses**1.2))[:, None]
     adata = AnnData(X, obs=pd.DataFrame({"perturbation": groups, "dose": doses}))
+    if categorical_doses:
+        adata.obs["dose"] = pd.Categorical(adata.obs["dose"].astype(str))
     sc.pp.pca(adata, n_comps=5)
 
     curves = pt.tl.PseudobulkSpace().dose_response(adata, dose_col="dose", metric="euclidean", embedding_key="X_pca")
     drug = curves[curves["perturbation"] == "drug"].sort_values("dose")
     assert drug["distance"].is_monotonic_increasing
 
-    fits = pt.tl.PseudobulkSpace().fit_dose_response(curves)
-    assert fits.loc[0, "ec50"] == pytest.approx(10, rel=0.2)
+    pt.tl.PseudobulkSpace().fit_dose_response(adata, curves)
+    fits = adata.uns["dose_response"]["fits"]
+    assert fits.loc[0, "ec50"] == pytest.approx(10, rel=1e-4)
     assert fits.loc[0, "r_squared"] > 0.99
     assert fits.loc[0, "midpoint_in_range"]
 
@@ -96,7 +100,7 @@ def test_dose_response():
     [("effect", 0.1, 1.8, 3.0, "ec50"), ("inhibition", 1.0, 0.05, 8.0, "ic50")],
 )
 @pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
-def test_fit_dose_response(*, response_type, e0, emax, midpoint, midpoint_col, response_scale):
+def test_fit_dose_response(*, adata, response_type, e0, emax, midpoint, midpoint_col, response_scale):
     e0, emax = e0 * response_scale, emax * response_scale
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0])
     hill_coefficient = 1.4
@@ -113,17 +117,15 @@ def test_fit_dose_response(*, response_type, e0, emax, midpoint, midpoint_col, r
     second = first.assign(compound="drug_b", response=e0 + (emax - e0) * second_fraction)
     data = pd.concat([first, second], ignore_index=True)
 
-    fits = (
-        pt.tl.PseudobulkSpace()
-        .fit_dose_response(
-            data,
-            perturbation_col="compound",
-            dose_col="concentration",
-            response_col="response",
-            response_type=response_type,
-        )
-        .set_index("compound")
+    pt.tl.PseudobulkSpace().fit_dose_response(
+        adata,
+        data,
+        perturbation_col="compound",
+        dose_col="concentration",
+        response_col="response",
+        response_type=response_type,
     )
+    fits = adata.uns["dose_response"]["fits"].set_index("compound")
 
     assert midpoint_col in fits
     assert f"{midpoint_col}_standard_error" in fits
@@ -138,14 +140,15 @@ def test_fit_dose_response(*, response_type, e0, emax, midpoint, midpoint_col, r
 
 
 @pytest.mark.parametrize("response_scale", [1.0, 100.0])
-def test_fit_dose_response_rank_deficient(response_scale):
+def test_fit_dose_response_rank_deficient(adata, response_scale):
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
     responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
     complete = pd.DataFrame({"perturbation": "complete", "dose": doses, "distance": response_scale * responses})
     limited = complete.loc[complete["dose"] <= 3].assign(perturbation="limited")
 
     with pytest.warns(UserWarning, match="Cannot estimate.*'limited'"):
-        fits = pt.tl.PseudobulkSpace().fit_dose_response(pd.concat([complete, limited])).set_index("perturbation")
+        pt.tl.PseudobulkSpace().fit_dose_response(adata, pd.concat([complete, limited]))
+    fits = adata.uns["dose_response"]["fits"].set_index("perturbation")
 
     assert fits.loc["complete", "ec50"] == pytest.approx(10, rel=0.1)
     assert np.isfinite(fits.loc["complete", "ec50_standard_error"])
@@ -157,12 +160,13 @@ def test_fit_dose_response_rank_deficient(response_scale):
 
 
 @pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
-def test_fit_dose_response_standard_error(response_scale):
+def test_fit_dose_response_standard_error(adata, response_scale):
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
     responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
     responses *= response_scale
     data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": responses})
-    fit = pt.tl.PseudobulkSpace().fit_dose_response(data).iloc[0]
+    pt.tl.PseudobulkSpace().fit_dose_response(adata, data)
+    fit = adata.uns["dose_response"]["fits"].iloc[0]
 
     # Independent derivatives with respect to EC50 itself, rather than the fitted log(EC50).
     midpoint, slope = fit["ec50"], fit["hill_coefficient"]
@@ -177,15 +181,49 @@ def test_fit_dose_response_standard_error(response_scale):
     assert fit["ec50_standard_error"] == pytest.approx(expected_error, rel=1e-4)
 
 
-def test_fit_dose_response_insufficient_dof():
+def test_fit_dose_response_insufficient_dof(adata):
     doses = np.array([0.0, 1.0, 10.0, 100.0])
     data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": doses / (10 + doses)})
     with pytest.warns(UserWarning, match="Cannot estimate.*'drug'"):
-        fit = pt.tl.PseudobulkSpace().fit_dose_response(data).iloc[0]
+        pt.tl.PseudobulkSpace().fit_dose_response(adata, data)
+    fit = adata.uns["dose_response"]["fits"].iloc[0]
     assert fit["ec50"] == pytest.approx(10)
     assert np.isnan(fit["ec50_standard_error"])
 
 
+def test_fit_dose_response_storage(adata, tmp_path):
+    doses = np.array([0.0, 0.1, 1.0, 3.0, 10.0, 30.0, 100.0])
+    data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": doses / (10 + doses)})
+    original_obs, original_data = adata.obs.copy(), data.copy()
+    original_x = adata.X.copy()
+    adata.uns["other_result"] = {"value": 1}
+
+    result = pt.tl.PseudobulkSpace().fit_dose_response(adata, data, key_added="hill")
+
+    assert result is None
+    assert "dose_response" not in adata.uns
+    assert adata.uns["other_result"] == {"value": 1}
+    pd.testing.assert_frame_equal(adata.obs, original_obs)
+    pd.testing.assert_frame_equal(data, original_data)
+    np.testing.assert_array_equal(adata.X, original_x)
+    assert adata.uns["hill"]["params"] == {
+        "perturbation_col": "perturbation",
+        "dose_col": "dose",
+        "response_col": "distance",
+        "response_type": "effect",
+    }
+
+    path = tmp_path / "fits.h5ad"
+    adata.write_h5ad(path)
+    restored = sc.read_h5ad(path)
+    pd.testing.assert_frame_equal(restored.uns["hill"]["fits"], adata.uns["hill"]["fits"])
+    assert restored.uns["hill"]["params"] == adata.uns["hill"]["params"]
+
+    pt.tl.PseudobulkSpace().fit_dose_response(adata, data, response_type="inhibition", key_added="hill")
+    assert "ic50" in adata.uns["hill"]["fits"]
+    assert "ec50" not in adata.uns["hill"]["fits"]
+
+
 @pytest.mark.parametrize(
     ("data", "kwargs", "match"),
     [
@@ -219,6 +257,7 @@ def test_fit_dose_response_insufficient_dof():
         ),
     ],
 )
-def test_fit_dose_response_validation(data, kwargs, match):
+def test_fit_dose_response_validation(adata, data, kwargs, match):
     with pytest.raises(ValueError, match=match):
-        pt.tl.PseudobulkSpace().fit_dose_response(data, **kwargs)
+        pt.tl.PseudobulkSpace().fit_dose_response(adata, data, **kwargs)
+    assert "dose_response" not in adata.uns

From c243d08fe44697219b52c83995d2ead2ac9f6ecc Mon Sep 17 00:00:00 2001
From: Lukas Heumos <lukas.heumos@posteo.net>
Date: Wed, 23 Sep 2026 15:21:02 +0200
Subject: [PATCH 3/8] Return Hill fits as a DataFrame and survive unfittable
 perturbations

fit_dose_response only used the AnnData as a sink for .uns, so callers with
assay tables had to build an empty AnnData. It now takes the tidy table and
returns the fits, matching dose_response and evaluate_combinations.

A perturbation with too few doses, a constant response or a non-converging
fit now warns and gets a NaN row instead of aborting the whole screen.
Drop response_type, which only renamed ec50 to ic50, and document that e0
is extrapolated when the reference group is absent from the input.
---
 docs/api/tools_index.md                       |   7 +-
 .../_perturbation_space.py                    | 187 ++++++++----------
 .../test_perturbation_space_extras.py         | 146 +++++---------
 3 files changed, 131 insertions(+), 209 deletions(-)

diff --git a/docs/api/tools_index.md b/docs/api/tools_index.md
index 07cf3262..5cde0cdf 100644
--- a/docs/api/tools_index.md
+++ b/docs/api/tools_index.md
@@ -569,7 +569,7 @@ See [perturbation space tutorial](https://pertpy.readthedocs.io/en/latest/tutori
 ### Dose-response curve fitting
 
 {meth}`~pertpy.tools.PseudobulkSpace.dose_response` calculates a scalar distance from control for each perturbation and dose.
-{meth}`~pertpy.tools.PseudobulkSpace.fit_dose_response` fits a four-parameter Hill curve and stores the results in AnnData.
+{meth}`~pertpy.tools.PseudobulkSpace.fit_dose_response` fits a four-parameter Hill curve per perturbation and returns EC50 estimates.
 For preprocessed `adata` with `perturbation` and `dose` columns, a `control` group and a PCA representation:
 
 ```python
@@ -577,8 +577,7 @@ import pertpy as pt
 
 ps = pt.tl.PseudobulkSpace()
 responses = ps.dose_response(adata, embedding_key="X_pca")
-ps.fit_dose_response(adata, responses)
-fits = adata.uns["dose_response"]["fits"]
+fits = ps.fit_dose_response(responses)
 ```
 
-For an inhibitory assay response, use `response_col` to select the measurement column and `response_type="inhibition"` to report IC50 instead of EC50.
+`fit_dose_response` also accepts any tidy table of scalar assay responses via `response_col`.
diff --git a/src/pertpy/tools/_perturbation_space/_perturbation_space.py b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
index 1a3405f8..1d33eb78 100644
--- a/src/pertpy/tools/_perturbation_space/_perturbation_space.py
+++ b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
@@ -95,6 +95,61 @@ def _four_parameter_logistic(
     return e0 + (emax - e0) * fraction
 
 
+def _fit_hill(doses: np.ndarray, responses: np.ndarray) -> dict[str, float | bool]:
+    """Fit a four-parameter Hill curve to one perturbation's responses."""
+    unique_doses, dose_index = np.unique(doses, return_inverse=True)
+    if unique_doses.size < 4:
+        raise ValueError("at least four distinct dose values are required.")
+    response_offset = float(np.min(responses))
+    response_scale = float(np.ptp(responses))
+    if response_scale == 0:
+        raise ValueError("the response is constant.")
+
+    responses = (responses - response_offset) / response_scale
+    mean_response = np.bincount(dose_index, weights=responses) / np.bincount(dose_index)
+    e0_guess, emax_guess = float(mean_response[0]), float(mean_response[-1])
+    positive = unique_doses > 0
+    halfway_distance = np.abs(mean_response[positive] - (e0_guess + emax_guess) / 2)
+    midpoint_guess = float(unique_doses[positive][np.argmin(halfway_distance)])
+
+    parameters, covariance = curve_fit(
+        _four_parameter_logistic,
+        doses,
+        responses,
+        p0=(e0_guess, emax_guess, np.log(midpoint_guess), 1.0),
+        bounds=((-np.inf, -np.inf, -np.inf, np.finfo(float).eps), np.inf),
+        absolute_sigma=True,
+        maxfev=20_000,
+    )
+    e0, emax, log_midpoint, hill_coefficient = parameters
+    midpoint = float(np.exp(log_midpoint))
+    residual_sum_squares = float(np.sum((responses - _four_parameter_logistic(doses, *parameters)) ** 2))
+    total_sum_squares = float(np.sum((responses - responses.mean()) ** 2))
+    positive_doses = doses[doses > 0]
+
+    # Check conditioning before scaling covariance by residual variance, which may be zero.
+    degrees_of_freedom = len(responses) - len(parameters)
+    if (
+        degrees_of_freedom <= 0
+        or not np.isfinite(covariance).all()
+        or np.linalg.matrix_rank(covariance) < len(parameters)
+    ):
+        midpoint_standard_error = np.nan
+    else:
+        residual_variance = residual_sum_squares / degrees_of_freedom
+        midpoint_standard_error = midpoint * np.sqrt(float(covariance[2, 2]) * residual_variance)
+
+    return {
+        "e0": float(e0 * response_scale + response_offset),
+        "emax": float(emax * response_scale + response_offset),
+        "hill_coefficient": float(hill_coefficient),
+        "ec50": midpoint,
+        "ec50_standard_error": float(midpoint_standard_error),
+        "r_squared": 1 - residual_sum_squares / total_sum_squares,
+        "midpoint_in_range": bool(positive_doses.min() < midpoint < positive_doses.max()),
+    }
+
+
 def _subtract_control_mean(
     matrix: np.ndarray,
     control_mask: np.ndarray,
@@ -662,68 +717,43 @@ def dose_response(
 
     def fit_dose_response(
         self,
-        adata: AnnData,
         data: pd.DataFrame,
         *,
         perturbation_col: str = "perturbation",
         dose_col: str = "dose",
         response_col: str = "distance",
-        response_type: Literal["effect", "inhibition"] = "effect",
-        key_added: str = "dose_response",
-    ) -> None:
+    ) -> pd.DataFrame:
         """Fit a four-parameter Hill curve for each perturbation.
 
         ``data`` can be the output of :meth:`dose_response` or a table containing another scalar assay response.
-        ``response_type`` names the fitted midpoint according to the meaning of that response: ``"effect"`` reports
-        ``ec50``, while ``"inhibition"`` reports ``ic50``.
         It does not perform biological or control normalization.
-        Dose values must represent concentrations for the EC50 or IC50 terminology to apply.
+        Perturbations whose curve cannot be fit are reported with a warning and NaN parameters.
 
         Args:
-            adata: AnnData to store the fit results in.
             data: Tidy table containing perturbation, dose and response columns.
             perturbation_col: Column identifying the perturbation.
             dose_col: Column containing non-negative numeric doses.
             response_col: Column containing the scalar response to fit.
-            response_type: Whether the response represents an effect or inhibition.
-            key_added: Key in `.uns` for the fit results and parameters.
 
         Returns:
-            Updates `.uns[key_added]` with a ``fits`` table and a ``params`` dictionary.
-            The table contains one row per perturbation with the fitted zero-dose response (``e0``), asymptotic
-            response (``emax``), Hill coefficient, EC50 or IC50, its approximate standard error, R-squared and
-            whether the midpoint lies within the tested positive-dose range.
+            DataFrame with one row per perturbation containing the fitted response at zero dose (``e0``), the asymptotic response (``emax``), the Hill coefficient, ``ec50``, its approximate standard error, R-squared and whether ``ec50`` lies within the tested positive-dose range.
 
         Notes:
-            Fits are stored in `.uns` because they describe groups of observations, not individual cells or genes.
-            EC50 and IC50 are relative midpoints between the fitted ``e0`` and ``emax``.
-            Responses are rescaled internally for numerical stability; ``e0`` and ``emax`` retain the input units.
-            The standard error uses a local linear approximation. It is NaN, with a warning, if the parameter
-            covariance is non-finite or numerically rank deficient, or there are no residual degrees of freedom.
-            Parameters are retained for inspection.
-            A small standard error, high R-squared or an in-range midpoint
-            does not establish that the doses capture both plateaus or that the Hill model is appropriate.
+            EC50 is the relative midpoint between the fitted ``e0`` and ``emax``.
+            :meth:`dose_response` does not return the reference group, so ``e0`` is extrapolated below the lowest tested dose unless ``data`` contains zero doses.
+            The standard error uses a local linear approximation and is NaN, with a warning, if the parameters are not identifiable from the data.
+            A small standard error, high R-squared or an in-range EC50 does not establish that the doses capture both plateaus.
 
         Examples:
-            A minimal example with a simulated response:
-
             >>> import pertpy as pt
-            >>> import pandas as pd
-            >>> from anndata import AnnData
-            >>> adata = AnnData(obs=pd.DataFrame(index=["drug"]))
-            >>> doses = pd.Series([0, 1, 3, 10, 30, 100])
-            >>> responses = pd.DataFrame({"dose": doses, "distance": doses / (10 + doses)})
-            >>> responses["perturbation"] = "drug"
-            >>> pt.tl.PseudobulkSpace().fit_dose_response(adata, responses)
-            >>> fits = adata.uns["dose_response"]["fits"]
+            >>> adata = pt.dt.srivatsan_2020_sciplex2()
+            >>> ps = pt.tl.PseudobulkSpace()
+            >>> responses = ps.dose_response(adata, dose_col="dose_value", embedding_key="X_pca")
+            >>> fits = ps.fit_dose_response(responses)
         """
-        required = {perturbation_col, dose_col, response_col}
-        missing = required.difference(data.columns)
+        missing = {perturbation_col, dose_col, response_col}.difference(data.columns)
         if missing:
             raise ValueError(f"Columns {sorted(missing)} do not exist in the input data.")
-        if response_type not in {"effect", "inhibition"}:
-            raise ValueError("response_type must be either 'effect' or 'inhibition'.")
-
         fit_data = data[[perturbation_col, dose_col, response_col]].copy()
         if fit_data[perturbation_col].isna().any():
             raise ValueError("Perturbation labels must not be missing.")
@@ -732,87 +762,30 @@ def fit_dose_response(
         if (fit_data[dose_col] < 0).any():
             raise ValueError("Dose values must be non-negative.")
 
-        midpoint_col = "ec50" if response_type == "effect" else "ic50"
         records: list[dict[str, object]] = []
         for perturbation, group in fit_data.groupby(perturbation_col, observed=True, sort=True):
             doses = group[dose_col].to_numpy(dtype=float)
             responses = group[response_col].to_numpy(dtype=float)
-            if np.unique(doses).size < 4:
-                raise ValueError(f"Perturbation {perturbation!r} needs at least four distinct dose values.")
-            response_offset = float(np.min(responses))
-            response_scale = float(np.ptp(responses))
-            if response_scale == 0:
-                raise ValueError(
-                    f"Perturbation {perturbation!r} has a constant response, so a Hill curve cannot be fit."
+            try:
+                fit = _fit_hill(doses, responses)
+            except (ValueError, RuntimeError) as e:
+                warnings.warn(
+                    f"Cannot fit a Hill curve for perturbation {perturbation!r}: {e}", UserWarning, stacklevel=2
                 )
-
-            responses = (responses - response_offset) / response_scale
-            mean_response = (
-                group.groupby(dose_col, sort=True, observed=True)[response_col].mean() - response_offset
-            ) / response_scale
-            e0_guess = float(mean_response.iloc[0])
-            emax_guess = float(mean_response.iloc[-1])
-            halfway = (e0_guess + emax_guess) / 2
-            positive_response = mean_response[mean_response.index > 0]
-            midpoint_guess = float((positive_response - halfway).abs().idxmin())
-
-            parameters, covariance = curve_fit(
-                _four_parameter_logistic,
-                doses,
-                responses,
-                p0=(e0_guess, emax_guess, np.log(midpoint_guess), 1.0),
-                bounds=((-np.inf, -np.inf, -np.inf, np.finfo(float).eps), np.inf),
-                absolute_sigma=True,
-                maxfev=20_000,
-            )
-
-            e0, emax, log_midpoint, hill_coefficient = parameters
-            midpoint = float(np.exp(log_midpoint))
-            fitted = _four_parameter_logistic(doses, *parameters)
-            residual_sum_squares = float(np.sum((responses - fitted) ** 2))
-            total_sum_squares = float(np.sum((responses - responses.mean()) ** 2))
-            positive_doses = doses[doses > 0]
-
-            # Check conditioning before scaling covariance by residual variance, which may be zero.
-            degrees_of_freedom = len(responses) - len(parameters)
-            if (
-                degrees_of_freedom <= 0
-                or not np.isfinite(covariance).all()
-                or np.linalg.matrix_rank(covariance) < len(parameters)
-            ):
-                midpoint_standard_error = np.nan
+                fit = dict.fromkeys(
+                    ("e0", "emax", "hill_coefficient", "ec50", "ec50_standard_error", "r_squared"), np.nan
+                )
+                fit["midpoint_in_range"] = False
+            if np.isnan(fit["ec50_standard_error"]) and np.isfinite(fit["ec50"]):
                 warnings.warn(
-                    f"Cannot estimate the {midpoint_col.upper()} standard error for perturbation {perturbation!r}. "
+                    f"Cannot estimate the EC50 standard error for perturbation {perturbation!r}. "
                     "Inspect the dose range and fitted curve before interpreting the estimate.",
                     UserWarning,
                     stacklevel=2,
                 )
-            else:
-                residual_variance = residual_sum_squares / degrees_of_freedom
-                midpoint_standard_error = midpoint * np.sqrt(float(covariance[2, 2]) * residual_variance)
-
-            records.append(
-                {
-                    perturbation_col: perturbation,
-                    "e0": float(e0 * response_scale + response_offset),
-                    "emax": float(emax * response_scale + response_offset),
-                    "hill_coefficient": float(hill_coefficient),
-                    midpoint_col: midpoint,
-                    f"{midpoint_col}_standard_error": float(midpoint_standard_error),
-                    "r_squared": 1 - residual_sum_squares / total_sum_squares,
-                    "midpoint_in_range": bool(positive_doses.min() < midpoint < positive_doses.max()),
-                }
-            )
+            records.append({perturbation_col: perturbation, **fit})
 
-        adata.uns[key_added] = {
-            "fits": pd.DataFrame.from_records(records),
-            "params": {
-                "perturbation_col": perturbation_col,
-                "dose_col": dose_col,
-                "response_col": response_col,
-                "response_type": response_type,
-            },
-        }
+        return pd.DataFrame.from_records(records)
 
     def plot_similarity(  # pragma: no cover
         self,
diff --git a/tests/tools/_perturbation_space/test_perturbation_space_extras.py b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
index e26be0f0..1ea592db 100644
--- a/tests/tools/_perturbation_space/test_perturbation_space_extras.py
+++ b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
@@ -88,67 +88,53 @@ def test_dose_response(rng, categorical_doses):
     drug = curves[curves["perturbation"] == "drug"].sort_values("dose")
     assert drug["distance"].is_monotonic_increasing
 
-    pt.tl.PseudobulkSpace().fit_dose_response(adata, curves)
-    fits = adata.uns["dose_response"]["fits"]
+    fits = pt.tl.PseudobulkSpace().fit_dose_response(curves)
     assert fits.loc[0, "ec50"] == pytest.approx(10, rel=1e-4)
     assert fits.loc[0, "r_squared"] > 0.99
     assert fits.loc[0, "midpoint_in_range"]
 
 
-@pytest.mark.parametrize(
-    ("response_type", "e0", "emax", "midpoint", "midpoint_col"),
-    [("effect", 0.1, 1.8, 3.0, "ec50"), ("inhibition", 1.0, 0.05, 8.0, "ic50")],
-)
+@pytest.mark.parametrize(("e0", "emax", "midpoint"), [(0.1, 1.8, 3.0), (1.0, 0.05, 8.0)])
 @pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
-def test_fit_dose_response(*, adata, response_type, e0, emax, midpoint, midpoint_col, response_scale):
+def test_fit_dose_response(e0, emax, midpoint, response_scale):
     e0, emax = e0 * response_scale, emax * response_scale
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0])
     hill_coefficient = 1.4
-    fraction = doses**hill_coefficient / (midpoint**hill_coefficient + doses**hill_coefficient)
-    first = pd.DataFrame(
-        {
-            "compound": "drug_a",
-            "concentration": doses,
-            "response": e0 + (emax - e0) * fraction,
-        }
+
+    def hill(ec50):
+        return e0 + (emax - e0) * doses**hill_coefficient / (ec50**hill_coefficient + doses**hill_coefficient)
+
+    data = pd.concat(
+        [
+            pd.DataFrame({"compound": "drug_a", "concentration": doses, "response": hill(midpoint)}),
+            pd.DataFrame({"compound": "drug_b", "concentration": doses, "response": hill(2 * midpoint)}),
+        ]
     )
-    second_midpoint = midpoint * 2
-    second_fraction = doses**hill_coefficient / (second_midpoint**hill_coefficient + doses**hill_coefficient)
-    second = first.assign(compound="drug_b", response=e0 + (emax - e0) * second_fraction)
-    data = pd.concat([first, second], ignore_index=True)
-
-    pt.tl.PseudobulkSpace().fit_dose_response(
-        adata,
-        data,
-        perturbation_col="compound",
-        dose_col="concentration",
-        response_col="response",
-        response_type=response_type,
+    fits = (
+        pt.tl.PseudobulkSpace()
+        .fit_dose_response(data, perturbation_col="compound", dose_col="concentration", response_col="response")
+        .set_index("compound")
     )
-    fits = adata.uns["dose_response"]["fits"].set_index("compound")
 
-    assert midpoint_col in fits
-    assert f"{midpoint_col}_standard_error" in fits
-    assert {"e0", "emax", "hill_coefficient", "r_squared", "midpoint_in_range"} <= set(fits)
     assert fits.loc["drug_a", "e0"] == pytest.approx(e0, rel=1e-6, abs=0)
     assert fits.loc["drug_a", "emax"] == pytest.approx(emax, rel=1e-6, abs=0)
     assert fits.loc["drug_a", "hill_coefficient"] == pytest.approx(hill_coefficient)
-    assert fits.loc["drug_a", midpoint_col] == pytest.approx(midpoint)
-    assert fits.loc["drug_b", midpoint_col] == pytest.approx(second_midpoint)
+    assert fits.loc["drug_a", "ec50"] == pytest.approx(midpoint)
+    assert fits.loc["drug_b", "ec50"] == pytest.approx(2 * midpoint)
     assert fits.loc["drug_a", "r_squared"] == pytest.approx(1)
-    assert np.isfinite(fits[f"{midpoint_col}_standard_error"]).all()
+    assert np.isfinite(fits["ec50_standard_error"]).all()
 
 
 @pytest.mark.parametrize("response_scale", [1.0, 100.0])
-def test_fit_dose_response_rank_deficient(adata, response_scale):
+def test_fit_dose_response_rank_deficient(response_scale):
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
     responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
     complete = pd.DataFrame({"perturbation": "complete", "dose": doses, "distance": response_scale * responses})
     limited = complete.loc[complete["dose"] <= 3].assign(perturbation="limited")
 
     with pytest.warns(UserWarning, match="Cannot estimate.*'limited'"):
-        pt.tl.PseudobulkSpace().fit_dose_response(adata, pd.concat([complete, limited]))
-    fits = adata.uns["dose_response"]["fits"].set_index("perturbation")
+        fits = pt.tl.PseudobulkSpace().fit_dose_response(pd.concat([complete, limited]))
+    fits = fits.set_index("perturbation")
 
     assert fits.loc["complete", "ec50"] == pytest.approx(10, rel=0.1)
     assert np.isfinite(fits.loc["complete", "ec50_standard_error"])
@@ -160,13 +146,12 @@ def test_fit_dose_response_rank_deficient(adata, response_scale):
 
 
 @pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
-def test_fit_dose_response_standard_error(adata, response_scale):
+def test_fit_dose_response_standard_error(response_scale):
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
     responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
     responses *= response_scale
     data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": responses})
-    pt.tl.PseudobulkSpace().fit_dose_response(adata, data)
-    fit = adata.uns["dose_response"]["fits"].iloc[0]
+    fit = pt.tl.PseudobulkSpace().fit_dose_response(data).iloc[0]
 
     # Independent derivatives with respect to EC50 itself, rather than the fitted log(EC50).
     midpoint, slope = fit["ec50"], fit["hill_coefficient"]
@@ -181,83 +166,48 @@ def test_fit_dose_response_standard_error(adata, response_scale):
     assert fit["ec50_standard_error"] == pytest.approx(expected_error, rel=1e-4)
 
 
-def test_fit_dose_response_insufficient_dof(adata):
+def test_fit_dose_response_insufficient_dof():
     doses = np.array([0.0, 1.0, 10.0, 100.0])
     data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": doses / (10 + doses)})
     with pytest.warns(UserWarning, match="Cannot estimate.*'drug'"):
-        pt.tl.PseudobulkSpace().fit_dose_response(adata, data)
-    fit = adata.uns["dose_response"]["fits"].iloc[0]
+        fit = pt.tl.PseudobulkSpace().fit_dose_response(data).iloc[0]
     assert fit["ec50"] == pytest.approx(10)
     assert np.isnan(fit["ec50_standard_error"])
 
 
-def test_fit_dose_response_storage(adata, tmp_path):
-    doses = np.array([0.0, 0.1, 1.0, 3.0, 10.0, 30.0, 100.0])
-    data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": doses / (10 + doses)})
-    original_obs, original_data = adata.obs.copy(), data.copy()
-    original_x = adata.X.copy()
-    adata.uns["other_result"] = {"value": 1}
-
-    result = pt.tl.PseudobulkSpace().fit_dose_response(adata, data, key_added="hill")
-
-    assert result is None
-    assert "dose_response" not in adata.uns
-    assert adata.uns["other_result"] == {"value": 1}
-    pd.testing.assert_frame_equal(adata.obs, original_obs)
-    pd.testing.assert_frame_equal(data, original_data)
-    np.testing.assert_array_equal(adata.X, original_x)
-    assert adata.uns["hill"]["params"] == {
-        "perturbation_col": "perturbation",
-        "dose_col": "dose",
-        "response_col": "distance",
-        "response_type": "effect",
-    }
-
-    path = tmp_path / "fits.h5ad"
-    adata.write_h5ad(path)
-    restored = sc.read_h5ad(path)
-    pd.testing.assert_frame_equal(restored.uns["hill"]["fits"], adata.uns["hill"]["fits"])
-    assert restored.uns["hill"]["params"] == adata.uns["hill"]["params"]
+@pytest.mark.parametrize(
+    ("doses", "responses", "match"),
+    [([0, 1, 2], [0, 1, 2], "four distinct"), ([0, 1, 2, 3], [1, 1, 1, 1], "constant")],
+)
+def test_fit_dose_response_unfittable(doses, responses, match):
+    good_doses = np.array([0.0, 0.1, 1.0, 3.0, 10.0, 30.0, 100.0])
+    data = pd.concat(
+        [
+            pd.DataFrame({"perturbation": "good", "dose": good_doses, "distance": good_doses / (10 + good_doses)}),
+            pd.DataFrame({"perturbation": "bad", "dose": doses, "distance": responses}),
+        ]
+    )
+    with pytest.warns(UserWarning, match=f"'bad'.*{match}"):
+        fits = pt.tl.PseudobulkSpace().fit_dose_response(data).set_index("perturbation")
 
-    pt.tl.PseudobulkSpace().fit_dose_response(adata, data, response_type="inhibition", key_added="hill")
-    assert "ic50" in adata.uns["hill"]["fits"]
-    assert "ec50" not in adata.uns["hill"]["fits"]
+    assert fits.loc["good", "ec50"] == pytest.approx(10)
+    assert fits.loc["bad"].drop("midpoint_in_range").isna().all()
+    assert not fits.loc["bad", "midpoint_in_range"]
 
 
 @pytest.mark.parametrize(
-    ("data", "kwargs", "match"),
+    ("data", "match"),
     [
-        (pd.DataFrame(), {}, "Columns"),
+        (pd.DataFrame(), "Columns"),
         (
             pd.DataFrame(
                 {"perturbation": ["drug", "drug", None, "drug"], "dose": [0, 1, 2, 3], "distance": [0, 1, 2, 3]}
             ),
-            {},
             "Perturbation labels",
         ),
-        (
-            pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [-1, 1, 2, 3], "distance": [0, 1, 2, 3]}),
-            {},
-            "non-negative",
-        ),
-        (
-            pd.DataFrame({"perturbation": ["drug"] * 3, "dose": [0, 1, 2], "distance": [0, 1, 2]}),
-            {},
-            "four distinct",
-        ),
-        (
-            pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [0, 1, 2, 3], "distance": [1, 1, 1, 1]}),
-            {},
-            "constant response",
-        ),
-        (
-            pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [0, 1, 2, 3], "distance": [0, 1, 2, 3]}),
-            {"response_type": "unknown"},
-            "response_type",
-        ),
+        (pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [-1, 1, 2, 3], "distance": [0, 1, 2, 3]}), "non-negative"),
     ],
 )
-def test_fit_dose_response_validation(adata, data, kwargs, match):
+def test_fit_dose_response_validation(data, match):
     with pytest.raises(ValueError, match=match):
-        pt.tl.PseudobulkSpace().fit_dose_response(adata, data, **kwargs)
-    assert "dose_response" not in adata.uns
+        pt.tl.PseudobulkSpace().fit_dose_response(data)

From b4da166bf2e9ea7eaa2bf306301a1912c60d7973 Mon Sep 17 00:00:00 2001
From: Lukas Heumos <lukas.heumos@posteo.net>
Date: Wed, 23 Sep 2026 15:28:40 +0200
Subject: [PATCH 4/8] Keep dose-response analysis in AnnData

dose_response now returns a perturbation-by-dose AnnData (group means in X,
distance and group-constant obs columns in .obs) instead of a DataFrame,
matching PseudobulkSpace.compute. This is a breaking change to the 1.2.0 API.

fit_dose_response reads perturbation, dose and response from .obs and writes
the fitted response and the per-perturbation Hill parameters as hill_*
columns to .obs, so results stay in AnnData without .uns. Assay data works
as an AnnData whose .obs holds one row per well.
---
 docs/api/tools_index.md                       |  10 +-
 .../_perturbation_space.py                    | 119 ++++++++-------
 .../test_perturbation_space_extras.py         | 136 ++++++++++--------
 3 files changed, 152 insertions(+), 113 deletions(-)

diff --git a/docs/api/tools_index.md b/docs/api/tools_index.md
index 5cde0cdf..c4336a5d 100644
--- a/docs/api/tools_index.md
+++ b/docs/api/tools_index.md
@@ -568,16 +568,16 @@ See [perturbation space tutorial](https://pertpy.readthedocs.io/en/latest/tutori
 
 ### Dose-response curve fitting
 
-{meth}`~pertpy.tools.PseudobulkSpace.dose_response` calculates a scalar distance from control for each perturbation and dose.
-{meth}`~pertpy.tools.PseudobulkSpace.fit_dose_response` fits a four-parameter Hill curve per perturbation and returns EC50 estimates.
+{meth}`~pertpy.tools.PseudobulkSpace.dose_response` returns an AnnData with one observation per perturbation and dose, holding its distance from control in `.obs`.
+{meth}`~pertpy.tools.PseudobulkSpace.fit_dose_response` fits a four-parameter Hill curve per perturbation and stores EC50 and the other curve parameters in `.obs`.
 For preprocessed `adata` with `perturbation` and `dose` columns, a `control` group and a PCA representation:
 
 ```python
 import pertpy as pt
 
 ps = pt.tl.PseudobulkSpace()
-responses = ps.dose_response(adata, embedding_key="X_pca")
-fits = ps.fit_dose_response(responses)
+dose_adata = ps.dose_response(adata, embedding_key="X_pca")
+ps.fit_dose_response(dose_adata)
 ```
 
-`fit_dose_response` also accepts any tidy table of scalar assay responses via `response_col`.
+Assay measurements such as viability can be fit the same way from an AnnData whose `.obs` holds perturbation, dose and response columns.
diff --git a/src/pertpy/tools/_perturbation_space/_perturbation_space.py b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
index 1d33eb78..6dd2dde1 100644
--- a/src/pertpy/tools/_perturbation_space/_perturbation_space.py
+++ b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
@@ -6,6 +6,7 @@
 
 import numpy as np
 import pandas as pd
+import scanpy as sc
 from anndata import AnnData
 from scipy.optimize import curve_fit
 from scipy.special import expit
@@ -142,11 +143,11 @@ def _fit_hill(doses: np.ndarray, responses: np.ndarray) -> dict[str, float | boo
     return {
         "e0": float(e0 * response_scale + response_offset),
         "emax": float(emax * response_scale + response_offset),
-        "hill_coefficient": float(hill_coefficient),
+        "slope": float(hill_coefficient),
         "ec50": midpoint,
-        "ec50_standard_error": float(midpoint_standard_error),
+        "ec50_se": float(midpoint_standard_error),
         "r_squared": 1 - residual_sum_squares / total_sum_squares,
-        "midpoint_in_range": bool(positive_doses.min() < midpoint < positive_doses.max()),
+        "ec50_in_range": bool(positive_doses.min() < midpoint < positive_doses.max()),
     }
 
 
@@ -657,7 +658,7 @@ def dose_response(
         layer_key: str | None = None,
         embedding_key: str | None = None,
         **kwargs,
-    ) -> pd.DataFrame:
+    ) -> AnnData:
         """Quantify the effect size of each perturbation as a function of dose.
 
         For every (perturbation, dose) group the statistical distance to ``reference_key`` is computed in the chosen representation using :class:`~pertpy.tools.Distance`.
@@ -674,13 +675,15 @@ def dose_response(
             kwargs: Passed to :meth:`~pertpy.tools.Distance.onesided_distances`.
 
         Returns:
-            Tidy DataFrame with ``perturbation``, ``dose`` and ``distance`` columns, sorted by perturbation then dose.
+            AnnData with one observation per (perturbation, dose) group other than ``reference_key``, sorted by perturbation then dose.
+            ``X`` holds the group mean of the chosen representation.
+            `.obs` holds the ``distance`` and every `.obs` column that is constant within each group, including ``target_col`` and ``dose_col``.
 
         Examples:
             >>> import pertpy as pt
             >>> adata = pt.dt.srivatsan_2020_sciplex2()
             >>> ps = pt.tl.PseudobulkSpace()
-            >>> curves = ps.dose_response(adata, dose_col="dose_value", embedding_key="X_pca")
+            >>> dose_adata = ps.dose_response(adata, dose_col="dose_value", embedding_key="X_pca")
         """
         for col in (target_col, dose_col):
             if col not in adata.obs:
@@ -702,45 +705,51 @@ def dose_response(
         if isinstance(dists, tuple):
             dists = dists[0]
 
-        records = []
-        for label, value in dists.items():
-            if label == reference_key:
-                continue
-            perturbation, _, dose = str(label).partition(sep)
-            records.append({"perturbation": perturbation, "dose": dose, "distance": float(value)})
-        result = pd.DataFrame.from_records(records)
+        treated = grouped[~is_control].copy()
+        treated.obs["_dose_group"] = treated.obs["_dose_group"].cat.remove_unused_categories()
+        dose_adata = sc.get.aggregate(treated, by="_dose_group", func="mean", layer=layer_key, obsm=embedding_key)
+        dose_adata.X = dose_adata.layers.pop("mean")
+        _carry_constant_obs(dose_adata, cast_frame(treated.obs), "_dose_group")
+        dose_obs = cast_frame(dose_adata.obs)
+        dose_adata.obs["distance"] = dists.reindex(dose_obs["_dose_group"].astype(str)).to_numpy(dtype=float)
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             with contextlib.suppress(ValueError, TypeError):
-                result["dose"] = pd.to_numeric(result["dose"])
-        return result.sort_values(["perturbation", "dose"]).reset_index(drop=True)
+                dose_adata.obs[dose_col] = pd.to_numeric(dose_obs[dose_col].astype(str))
+        dose_adata.obs = cast_frame(dose_adata.obs).drop(columns="_dose_group")
+        dose_adata.obs_names = dose_adata.obs_names.str.replace(sep, "_")
+        order = cast_frame(dose_adata.obs).sort_values([target_col, dose_col]).index
+        return dose_adata[order].copy()
 
     def fit_dose_response(
         self,
-        data: pd.DataFrame,
+        adata: AnnData,
         *,
-        perturbation_col: str = "perturbation",
+        target_col: str = "perturbation",
         dose_col: str = "dose",
         response_col: str = "distance",
-    ) -> pd.DataFrame:
+        key_added: str = "hill",
+    ) -> None:
         """Fit a four-parameter Hill curve for each perturbation.
 
-        ``data`` can be the output of :meth:`dose_response` or a table containing another scalar assay response.
+        ``adata`` holds one observation per perturbation and dose or per replicate, such as the output of :meth:`dose_response` or assay measurements stored in `.obs`.
         It does not perform biological or control normalization.
         Perturbations whose curve cannot be fit are reported with a warning and NaN parameters.
 
         Args:
-            data: Tidy table containing perturbation, dose and response columns.
-            perturbation_col: Column identifying the perturbation.
-            dose_col: Column containing non-negative numeric doses.
-            response_col: Column containing the scalar response to fit.
+            adata: AnnData with perturbation, dose and response columns in `.obs`.
+            target_col: `.obs` column identifying the perturbation.
+            dose_col: `.obs` column containing non-negative numeric doses.
+            response_col: `.obs` column containing the scalar response to fit.
+            key_added: Prefix of the `.obs` columns the results are written to.
 
         Returns:
-            DataFrame with one row per perturbation containing the fitted response at zero dose (``e0``), the asymptotic response (``emax``), the Hill coefficient, ``ec50``, its approximate standard error, R-squared and whether ``ec50`` lies within the tested positive-dose range.
+            Adds ``{key_added}_fitted`` with the fitted response of every observation to `.obs`.
+            Also adds the per-perturbation parameters, repeated across each perturbation's observations: the fitted response at zero dose (``_e0``), the asymptotic response (``_emax``), the Hill slope (``_slope``), ``_ec50``, its approximate standard error (``_ec50_se``), R-squared (``_r_squared``) and whether EC50 lies within the tested positive-dose range (``_ec50_in_range``).
 
         Notes:
             EC50 is the relative midpoint between the fitted ``e0`` and ``emax``.
-            :meth:`dose_response` does not return the reference group, so ``e0`` is extrapolated below the lowest tested dose unless ``data`` contains zero doses.
+            :meth:`dose_response` does not return the reference group, so ``e0`` is extrapolated below the lowest tested dose unless ``adata`` contains zero doses.
             The standard error uses a local linear approximation and is NaN, with a warning, if the parameters are not identifiable from the data.
             A small standard error, high R-squared or an in-range EC50 does not establish that the doses capture both plateaus.
 
@@ -748,44 +757,50 @@ def fit_dose_response(
             >>> import pertpy as pt
             >>> adata = pt.dt.srivatsan_2020_sciplex2()
             >>> ps = pt.tl.PseudobulkSpace()
-            >>> responses = ps.dose_response(adata, dose_col="dose_value", embedding_key="X_pca")
-            >>> fits = ps.fit_dose_response(responses)
+            >>> dose_adata = ps.dose_response(adata, dose_col="dose_value", embedding_key="X_pca")
+            >>> ps.fit_dose_response(dose_adata, dose_col="dose_value")
         """
-        missing = {perturbation_col, dose_col, response_col}.difference(data.columns)
+        obs = cast_frame(adata.obs)
+        missing = {target_col, dose_col, response_col}.difference(obs.columns)
         if missing:
-            raise ValueError(f"Columns {sorted(missing)} do not exist in the input data.")
-        fit_data = data[[perturbation_col, dose_col, response_col]].copy()
-        if fit_data[perturbation_col].isna().any():
+            raise ValueError(f"Columns {sorted(missing)} do not exist in the .obs attribute.")
+        if obs[target_col].isna().any():
             raise ValueError("Perturbation labels must not be missing.")
-        fit_data[dose_col] = pd.to_numeric(fit_data[dose_col], errors="raise")
-        fit_data[response_col] = pd.to_numeric(fit_data[response_col], errors="raise")
-        if (fit_data[dose_col] < 0).any():
+        doses = obs[dose_col].to_numpy(dtype=float)
+        responses = obs[response_col].to_numpy(dtype=float)
+        if (doses < 0).any():
             raise ValueError("Dose values must be non-negative.")
 
-        records: list[dict[str, object]] = []
-        for perturbation, group in fit_data.groupby(perturbation_col, observed=True, sort=True):
-            doses = group[dose_col].to_numpy(dtype=float)
-            responses = group[response_col].to_numpy(dtype=float)
+        labels = obs[target_col].to_numpy()
+        fitted = np.full(len(obs), np.nan)
+        records: dict[object, dict[str, float | bool]] = {}
+        for perturbation in pd.unique(labels):
+            mask = labels == perturbation
             try:
-                fit = _fit_hill(doses, responses)
+                fit = _fit_hill(doses[mask], responses[mask])
             except (ValueError, RuntimeError) as e:
                 warnings.warn(
                     f"Cannot fit a Hill curve for perturbation {perturbation!r}: {e}", UserWarning, stacklevel=2
                 )
-                fit = dict.fromkeys(
-                    ("e0", "emax", "hill_coefficient", "ec50", "ec50_standard_error", "r_squared"), np.nan
+                fit = dict.fromkeys(("e0", "emax", "slope", "ec50", "ec50_se", "r_squared"), np.nan)
+                fit["ec50_in_range"] = False
+            else:
+                fitted[mask] = _four_parameter_logistic(
+                    doses[mask], fit["e0"], fit["emax"], np.log(fit["ec50"]), fit["slope"]
                 )
-                fit["midpoint_in_range"] = False
-            if np.isnan(fit["ec50_standard_error"]) and np.isfinite(fit["ec50"]):
-                warnings.warn(
-                    f"Cannot estimate the EC50 standard error for perturbation {perturbation!r}. "
-                    "Inspect the dose range and fitted curve before interpreting the estimate.",
-                    UserWarning,
-                    stacklevel=2,
-                )
-            records.append({perturbation_col: perturbation, **fit})
-
-        return pd.DataFrame.from_records(records)
+                if np.isnan(fit["ec50_se"]):
+                    warnings.warn(
+                        f"Cannot estimate the EC50 standard error for perturbation {perturbation!r}. "
+                        "Inspect the dose range and fitted curve before interpreting the estimate.",
+                        UserWarning,
+                        stacklevel=2,
+                    )
+            records[perturbation] = fit
+
+        fits = pd.DataFrame.from_dict(records, orient="index")
+        adata.obs[f"{key_added}_fitted"] = fitted
+        for col in fits.columns:
+            adata.obs[f"{key_added}_{col}"] = fits[col].reindex(labels).to_numpy()
 
     def plot_similarity(  # pragma: no cover
         self,
diff --git a/tests/tools/_perturbation_space/test_perturbation_space_extras.py b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
index 1ea592db..139d2364 100644
--- a/tests/tools/_perturbation_space/test_perturbation_space_extras.py
+++ b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
@@ -5,6 +5,7 @@
 from anndata import AnnData
 
 import pertpy as pt
+from pertpy._types import cast_frame
 
 
 @pytest.fixture
@@ -79,19 +80,33 @@ def test_dose_response(rng, categorical_doses):
     X = np.tile(rng.normal(0, 0.3, (15, 8)), (len(groups) // 15, 1))
     drug_doses = doses[groups == "drug"]
     X[groups == "drug"] += (5 * drug_doses**1.2 / (10**1.2 + drug_doses**1.2))[:, None]
-    adata = AnnData(X, obs=pd.DataFrame({"perturbation": groups, "dose": doses}))
+    adata = AnnData(X, obs=pd.DataFrame({"perturbation": groups, "dose": doses, "line": "A549"}))
     if categorical_doses:
         adata.obs["dose"] = pd.Categorical(adata.obs["dose"].astype(str))
     sc.pp.pca(adata, n_comps=5)
 
-    curves = pt.tl.PseudobulkSpace().dose_response(adata, dose_col="dose", metric="euclidean", embedding_key="X_pca")
-    drug = curves[curves["perturbation"] == "drug"].sort_values("dose")
-    assert drug["distance"].is_monotonic_increasing
+    ps = pt.tl.PseudobulkSpace()
+    dose_adata = ps.dose_response(adata, dose_col="dose", metric="euclidean", embedding_key="X_pca")
+    assert dose_adata.shape == (6, 5)
+    assert dose_adata.obs["dose"].tolist() == [0.1, 1.0, 3.0, 10.0, 30.0, 100.0]
+    assert (dose_adata.obs["line"] == "A549").all()
+    assert dose_adata.obs["distance"].is_monotonic_increasing
 
-    fits = pt.tl.PseudobulkSpace().fit_dose_response(curves)
-    assert fits.loc[0, "ec50"] == pytest.approx(10, rel=1e-4)
-    assert fits.loc[0, "r_squared"] > 0.99
-    assert fits.loc[0, "midpoint_in_range"]
+    ps.fit_dose_response(dose_adata)
+    assert dose_adata.obs["hill_ec50"].iloc[0] == pytest.approx(10, rel=1e-4)
+    assert dose_adata.obs["hill_r_squared"].iloc[0] > 0.99
+    assert dose_adata.obs["hill_ec50_in_range"].all()
+    np.testing.assert_allclose(dose_adata.obs["hill_fitted"], dose_adata.obs["distance"], rtol=1e-3)
+
+
+def _assay(data: pd.DataFrame) -> AnnData:
+    return AnnData(obs=data.reset_index(drop=True).rename(index=str))
+
+
+def _fits(adata: AnnData, target_col: str = "perturbation") -> pd.DataFrame:
+    obs = cast_frame(adata.obs)
+    cols = [col for col in obs.columns.astype(str) if col.startswith("hill_") and col != "hill_fitted"]
+    return obs.drop_duplicates(target_col).set_index(target_col)[cols]
 
 
 @pytest.mark.parametrize(("e0", "emax", "midpoint"), [(0.1, 1.8, 3.0), (1.0, 0.05, 8.0)])
@@ -99,30 +114,32 @@ def test_dose_response(rng, categorical_doses):
 def test_fit_dose_response(e0, emax, midpoint, response_scale):
     e0, emax = e0 * response_scale, emax * response_scale
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0])
-    hill_coefficient = 1.4
+    slope = 1.4
 
     def hill(ec50):
-        return e0 + (emax - e0) * doses**hill_coefficient / (ec50**hill_coefficient + doses**hill_coefficient)
-
-    data = pd.concat(
-        [
-            pd.DataFrame({"compound": "drug_a", "concentration": doses, "response": hill(midpoint)}),
-            pd.DataFrame({"compound": "drug_b", "concentration": doses, "response": hill(2 * midpoint)}),
-        ]
+        return e0 + (emax - e0) * doses**slope / (ec50**slope + doses**slope)
+
+    adata = _assay(
+        pd.concat(
+            [
+                pd.DataFrame({"compound": "drug_a", "concentration": doses, "response": hill(midpoint)}),
+                pd.DataFrame({"compound": "drug_b", "concentration": doses, "response": hill(2 * midpoint)}),
+            ]
+        )
     )
-    fits = (
-        pt.tl.PseudobulkSpace()
-        .fit_dose_response(data, perturbation_col="compound", dose_col="concentration", response_col="response")
-        .set_index("compound")
+    pt.tl.PseudobulkSpace().fit_dose_response(
+        adata, target_col="compound", dose_col="concentration", response_col="response"
     )
+    fits = _fits(adata, "compound")
 
-    assert fits.loc["drug_a", "e0"] == pytest.approx(e0, rel=1e-6, abs=0)
-    assert fits.loc["drug_a", "emax"] == pytest.approx(emax, rel=1e-6, abs=0)
-    assert fits.loc["drug_a", "hill_coefficient"] == pytest.approx(hill_coefficient)
-    assert fits.loc["drug_a", "ec50"] == pytest.approx(midpoint)
-    assert fits.loc["drug_b", "ec50"] == pytest.approx(2 * midpoint)
-    assert fits.loc["drug_a", "r_squared"] == pytest.approx(1)
-    assert np.isfinite(fits["ec50_standard_error"]).all()
+    assert fits.loc["drug_a", "hill_e0"] == pytest.approx(e0, rel=1e-6, abs=0)
+    assert fits.loc["drug_a", "hill_emax"] == pytest.approx(emax, rel=1e-6, abs=0)
+    assert fits.loc["drug_a", "hill_slope"] == pytest.approx(slope)
+    assert fits.loc["drug_a", "hill_ec50"] == pytest.approx(midpoint)
+    assert fits.loc["drug_b", "hill_ec50"] == pytest.approx(2 * midpoint)
+    assert fits.loc["drug_a", "hill_r_squared"] == pytest.approx(1)
+    assert np.isfinite(fits["hill_ec50_se"]).all()
+    np.testing.assert_allclose(adata.obs["hill_fitted"], adata.obs["response"], rtol=1e-6)
 
 
 @pytest.mark.parametrize("response_scale", [1.0, 100.0])
@@ -131,18 +148,19 @@ def test_fit_dose_response_rank_deficient(response_scale):
     responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
     complete = pd.DataFrame({"perturbation": "complete", "dose": doses, "distance": response_scale * responses})
     limited = complete.loc[complete["dose"] <= 3].assign(perturbation="limited")
+    adata = _assay(pd.concat([complete, limited]))
 
     with pytest.warns(UserWarning, match="Cannot estimate.*'limited'"):
-        fits = pt.tl.PseudobulkSpace().fit_dose_response(pd.concat([complete, limited]))
-    fits = fits.set_index("perturbation")
+        pt.tl.PseudobulkSpace().fit_dose_response(adata)
+    fits = _fits(adata)
 
-    assert fits.loc["complete", "ec50"] == pytest.approx(10, rel=0.1)
-    assert np.isfinite(fits.loc["complete", "ec50_standard_error"])
+    assert fits.loc["complete", "hill_ec50"] == pytest.approx(10, rel=0.1)
+    assert np.isfinite(fits.loc["complete", "hill_ec50_se"])
     # A high R-squared and an in-range midpoint do not expose this poorly determined curve.
-    assert fits.loc["limited", "r_squared"] > 0.98
-    assert fits.loc["limited", "midpoint_in_range"]
-    assert np.isfinite(fits.loc["limited", "ec50"])
-    assert np.isnan(fits.loc["limited", "ec50_standard_error"])
+    assert fits.loc["limited", "hill_r_squared"] > 0.98
+    assert fits.loc["limited", "hill_ec50_in_range"]
+    assert np.isfinite(fits.loc["limited", "hill_ec50"])
+    assert np.isnan(fits.loc["limited", "hill_ec50_se"])
 
 
 @pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
@@ -150,29 +168,31 @@ def test_fit_dose_response_standard_error(response_scale):
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
     responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
     responses *= response_scale
-    data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": responses})
-    fit = pt.tl.PseudobulkSpace().fit_dose_response(data).iloc[0]
+    adata = _assay(pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": responses}))
+    pt.tl.PseudobulkSpace().fit_dose_response(adata)
+    fit = _fits(adata).iloc[0]
 
     # Independent derivatives with respect to EC50 itself, rather than the fitted log(EC50).
-    midpoint, slope = fit["ec50"], fit["hill_coefficient"]
+    midpoint, slope = fit["hill_ec50"], fit["hill_slope"]
     fraction = doses**slope / (midpoint**slope + doses**slope)
     log_ratio = np.zeros_like(doses)
     np.log(doses / midpoint, out=log_ratio, where=doses > 0)
-    sensitivity = (fit["emax"] - fit["e0"]) * fraction * (1 - fraction)
+    sensitivity = (fit["hill_emax"] - fit["hill_e0"]) * fraction * (1 - fraction)
     jacobian = np.column_stack((1 - fraction, fraction, -slope * sensitivity / midpoint, sensitivity * log_ratio))
-    residuals = responses - (fit["e0"] + (fit["emax"] - fit["e0"]) * fraction)
+    residuals = responses - (fit["hill_e0"] + (fit["hill_emax"] - fit["hill_e0"]) * fraction)
     variance = np.sum(residuals**2) / (len(doses) - 4)
     expected_error = np.sqrt(np.linalg.inv(jacobian.T @ jacobian)[2, 2] * variance)
-    assert fit["ec50_standard_error"] == pytest.approx(expected_error, rel=1e-4)
+    assert fit["hill_ec50_se"] == pytest.approx(expected_error, rel=1e-4)
 
 
 def test_fit_dose_response_insufficient_dof():
     doses = np.array([0.0, 1.0, 10.0, 100.0])
-    data = pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": doses / (10 + doses)})
+    adata = _assay(pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": doses / (10 + doses)}))
     with pytest.warns(UserWarning, match="Cannot estimate.*'drug'"):
-        fit = pt.tl.PseudobulkSpace().fit_dose_response(data).iloc[0]
-    assert fit["ec50"] == pytest.approx(10)
-    assert np.isnan(fit["ec50_standard_error"])
+        pt.tl.PseudobulkSpace().fit_dose_response(adata)
+    fit = _fits(adata).iloc[0]
+    assert fit["hill_ec50"] == pytest.approx(10)
+    assert np.isnan(fit["hill_ec50_se"])
 
 
 @pytest.mark.parametrize(
@@ -181,24 +201,28 @@ def test_fit_dose_response_insufficient_dof():
 )
 def test_fit_dose_response_unfittable(doses, responses, match):
     good_doses = np.array([0.0, 0.1, 1.0, 3.0, 10.0, 30.0, 100.0])
-    data = pd.concat(
-        [
-            pd.DataFrame({"perturbation": "good", "dose": good_doses, "distance": good_doses / (10 + good_doses)}),
-            pd.DataFrame({"perturbation": "bad", "dose": doses, "distance": responses}),
-        ]
+    adata = _assay(
+        pd.concat(
+            [
+                pd.DataFrame({"perturbation": "good", "dose": good_doses, "distance": good_doses / (10 + good_doses)}),
+                pd.DataFrame({"perturbation": "bad", "dose": doses, "distance": responses}),
+            ]
+        )
     )
     with pytest.warns(UserWarning, match=f"'bad'.*{match}"):
-        fits = pt.tl.PseudobulkSpace().fit_dose_response(data).set_index("perturbation")
+        pt.tl.PseudobulkSpace().fit_dose_response(adata)
+    fits = _fits(adata)
 
-    assert fits.loc["good", "ec50"] == pytest.approx(10)
-    assert fits.loc["bad"].drop("midpoint_in_range").isna().all()
-    assert not fits.loc["bad", "midpoint_in_range"]
+    assert fits.loc["good", "hill_ec50"] == pytest.approx(10)
+    assert fits.loc["bad"].drop("hill_ec50_in_range").isna().all()
+    assert not fits.loc["bad", "hill_ec50_in_range"]
+    assert adata.obs.loc[adata.obs["perturbation"] == "bad", "hill_fitted"].isna().all()
 
 
 @pytest.mark.parametrize(
     ("data", "match"),
     [
-        (pd.DataFrame(), "Columns"),
+        (pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [0, 1, 2, 3]}), "Columns"),
         (
             pd.DataFrame(
                 {"perturbation": ["drug", "drug", None, "drug"], "dose": [0, 1, 2, 3], "distance": [0, 1, 2, 3]}
@@ -210,4 +234,4 @@ def test_fit_dose_response_unfittable(doses, responses, match):
 )
 def test_fit_dose_response_validation(data, match):
     with pytest.raises(ValueError, match=match):
-        pt.tl.PseudobulkSpace().fit_dose_response(data)
+        pt.tl.PseudobulkSpace().fit_dose_response(_assay(data))

From 4feca59ed8922e843c9d64959cdaa6f6419c02bd Mon Sep 17 00:00:00 2001
From: Lukas Heumos <lukas.heumos@posteo.net>
Date: Wed, 23 Sep 2026 15:33:49 +0200
Subject: [PATCH 5/8] Trim dose-response tests and drop redundant label
 validation

---
 .../_perturbation_space.py                    |   2 -
 .../test_perturbation_space_extras.py         | 143 +++++-------------
 2 files changed, 34 insertions(+), 111 deletions(-)

diff --git a/src/pertpy/tools/_perturbation_space/_perturbation_space.py b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
index 6dd2dde1..7d5b0d5e 100644
--- a/src/pertpy/tools/_perturbation_space/_perturbation_space.py
+++ b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
@@ -764,8 +764,6 @@ def fit_dose_response(
         missing = {target_col, dose_col, response_col}.difference(obs.columns)
         if missing:
             raise ValueError(f"Columns {sorted(missing)} do not exist in the .obs attribute.")
-        if obs[target_col].isna().any():
-            raise ValueError("Perturbation labels must not be missing.")
         doses = obs[dose_col].to_numpy(dtype=float)
         responses = obs[response_col].to_numpy(dtype=float)
         if (doses < 0).any():
diff --git a/tests/tools/_perturbation_space/test_perturbation_space_extras.py b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
index 139d2364..f89be83e 100644
--- a/tests/tools/_perturbation_space/test_perturbation_space_extras.py
+++ b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
@@ -76,11 +76,10 @@ def test_dose_response(rng, categorical_doses):
             doses += [dose] * 15
     groups = np.array(groups)
     doses = np.array(doses, dtype=float)
-    # Match background cells across doses so sampling noise does not change the known curve.
     X = np.tile(rng.normal(0, 0.3, (15, 8)), (len(groups) // 15, 1))
     drug_doses = doses[groups == "drug"]
     X[groups == "drug"] += (5 * drug_doses**1.2 / (10**1.2 + drug_doses**1.2))[:, None]
-    adata = AnnData(X, obs=pd.DataFrame({"perturbation": groups, "dose": doses, "line": "A549"}))
+    adata = AnnData(X, obs=pd.DataFrame({"perturbation": groups, "dose": doses}))
     if categorical_doses:
         adata.obs["dose"] = pd.Categorical(adata.obs["dose"].astype(str))
     sc.pp.pca(adata, n_comps=5)
@@ -89,90 +88,48 @@ def test_dose_response(rng, categorical_doses):
     dose_adata = ps.dose_response(adata, dose_col="dose", metric="euclidean", embedding_key="X_pca")
     assert dose_adata.shape == (6, 5)
     assert dose_adata.obs["dose"].tolist() == [0.1, 1.0, 3.0, 10.0, 30.0, 100.0]
-    assert (dose_adata.obs["line"] == "A549").all()
     assert dose_adata.obs["distance"].is_monotonic_increasing
 
     ps.fit_dose_response(dose_adata)
     assert dose_adata.obs["hill_ec50"].iloc[0] == pytest.approx(10, rel=1e-4)
-    assert dose_adata.obs["hill_r_squared"].iloc[0] > 0.99
     assert dose_adata.obs["hill_ec50_in_range"].all()
-    np.testing.assert_allclose(dose_adata.obs["hill_fitted"], dose_adata.obs["distance"], rtol=1e-3)
 
 
 def _assay(data: pd.DataFrame) -> AnnData:
     return AnnData(obs=data.reset_index(drop=True).rename(index=str))
 
 
-def _fits(adata: AnnData, target_col: str = "perturbation") -> pd.DataFrame:
-    obs = cast_frame(adata.obs)
-    cols = [col for col in obs.columns.astype(str) if col.startswith("hill_") and col != "hill_fitted"]
-    return obs.drop_duplicates(target_col).set_index(target_col)[cols]
+def _fits(adata: AnnData) -> pd.DataFrame:
+    return cast_frame(adata.obs).drop_duplicates("perturbation").set_index("perturbation")
 
 
-@pytest.mark.parametrize(("e0", "emax", "midpoint"), [(0.1, 1.8, 3.0), (1.0, 0.05, 8.0)])
-@pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
-def test_fit_dose_response(e0, emax, midpoint, response_scale):
-    e0, emax = e0 * response_scale, emax * response_scale
+@pytest.mark.parametrize(("e0", "emax"), [(0.1, 1.8), (1.0, 0.05)])
+def test_fit_dose_response(e0, emax):
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0])
-    slope = 1.4
-
-    def hill(ec50):
-        return e0 + (emax - e0) * doses**slope / (ec50**slope + doses**slope)
-
-    adata = _assay(
-        pd.concat(
-            [
-                pd.DataFrame({"compound": "drug_a", "concentration": doses, "response": hill(midpoint)}),
-                pd.DataFrame({"compound": "drug_b", "concentration": doses, "response": hill(2 * midpoint)}),
-            ]
-        )
-    )
-    pt.tl.PseudobulkSpace().fit_dose_response(
-        adata, target_col="compound", dose_col="concentration", response_col="response"
-    )
-    fits = _fits(adata, "compound")
+    responses = e0 + (emax - e0) * doses**1.4 / (3.0**1.4 + doses**1.4)
+    adata = _assay(pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": responses}))
+    pt.tl.PseudobulkSpace().fit_dose_response(adata)
+    fit = _fits(adata).loc["drug"]
 
-    assert fits.loc["drug_a", "hill_e0"] == pytest.approx(e0, rel=1e-6, abs=0)
-    assert fits.loc["drug_a", "hill_emax"] == pytest.approx(emax, rel=1e-6, abs=0)
-    assert fits.loc["drug_a", "hill_slope"] == pytest.approx(slope)
-    assert fits.loc["drug_a", "hill_ec50"] == pytest.approx(midpoint)
-    assert fits.loc["drug_b", "hill_ec50"] == pytest.approx(2 * midpoint)
-    assert fits.loc["drug_a", "hill_r_squared"] == pytest.approx(1)
-    assert np.isfinite(fits["hill_ec50_se"]).all()
-    np.testing.assert_allclose(adata.obs["hill_fitted"], adata.obs["response"], rtol=1e-6)
+    assert fit["hill_e0"] == pytest.approx(e0)
+    assert fit["hill_emax"] == pytest.approx(emax)
+    assert fit["hill_slope"] == pytest.approx(1.4)
+    assert fit["hill_ec50"] == pytest.approx(3.0)
+    np.testing.assert_allclose(adata.obs["hill_fitted"], responses, rtol=1e-6)
 
 
-@pytest.mark.parametrize("response_scale", [1.0, 100.0])
-def test_fit_dose_response_rank_deficient(response_scale):
+def test_fit_dose_response_standard_error():
     doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
     responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
-    complete = pd.DataFrame({"perturbation": "complete", "dose": doses, "distance": response_scale * responses})
+    complete = pd.DataFrame({"perturbation": "complete", "dose": doses, "distance": responses})
     limited = complete.loc[complete["dose"] <= 3].assign(perturbation="limited")
     adata = _assay(pd.concat([complete, limited]))
 
     with pytest.warns(UserWarning, match="Cannot estimate.*'limited'"):
         pt.tl.PseudobulkSpace().fit_dose_response(adata)
     fits = _fits(adata)
+    fit = fits.loc["complete"]
 
-    assert fits.loc["complete", "hill_ec50"] == pytest.approx(10, rel=0.1)
-    assert np.isfinite(fits.loc["complete", "hill_ec50_se"])
-    # A high R-squared and an in-range midpoint do not expose this poorly determined curve.
-    assert fits.loc["limited", "hill_r_squared"] > 0.98
-    assert fits.loc["limited", "hill_ec50_in_range"]
-    assert np.isfinite(fits.loc["limited", "hill_ec50"])
-    assert np.isnan(fits.loc["limited", "hill_ec50_se"])
-
-
-@pytest.mark.parametrize("response_scale", [1e-8, 1.0, 1e8])
-def test_fit_dose_response_standard_error(response_scale):
-    doses = np.array([0.0, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0])
-    responses = doses / (10 + doses) + np.random.default_rng(2026).normal(0, 0.025, len(doses))
-    responses *= response_scale
-    adata = _assay(pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": responses}))
-    pt.tl.PseudobulkSpace().fit_dose_response(adata)
-    fit = _fits(adata).iloc[0]
-
-    # Independent derivatives with respect to EC50 itself, rather than the fitted log(EC50).
     midpoint, slope = fit["hill_ec50"], fit["hill_slope"]
     fraction = doses**slope / (midpoint**slope + doses**slope)
     log_ratio = np.zeros_like(doses)
@@ -181,57 +138,25 @@ def test_fit_dose_response_standard_error(response_scale):
     jacobian = np.column_stack((1 - fraction, fraction, -slope * sensitivity / midpoint, sensitivity * log_ratio))
     residuals = responses - (fit["hill_e0"] + (fit["hill_emax"] - fit["hill_e0"]) * fraction)
     variance = np.sum(residuals**2) / (len(doses) - 4)
-    expected_error = np.sqrt(np.linalg.inv(jacobian.T @ jacobian)[2, 2] * variance)
-    assert fit["hill_ec50_se"] == pytest.approx(expected_error, rel=1e-4)
+    assert fit["hill_ec50_se"] == pytest.approx(
+        np.sqrt(np.linalg.inv(jacobian.T @ jacobian)[2, 2] * variance), rel=1e-4
+    )
+    assert np.isnan(fits.loc["limited", "hill_ec50_se"])
 
 
-def test_fit_dose_response_insufficient_dof():
-    doses = np.array([0.0, 1.0, 10.0, 100.0])
-    adata = _assay(pd.DataFrame({"perturbation": "drug", "dose": doses, "distance": doses / (10 + doses)}))
-    with pytest.warns(UserWarning, match="Cannot estimate.*'drug'"):
-        pt.tl.PseudobulkSpace().fit_dose_response(adata)
-    fit = _fits(adata).iloc[0]
-    assert fit["hill_ec50"] == pytest.approx(10)
-    assert np.isnan(fit["hill_ec50_se"])
-
-
-@pytest.mark.parametrize(
-    ("doses", "responses", "match"),
-    [([0, 1, 2], [0, 1, 2], "four distinct"), ([0, 1, 2, 3], [1, 1, 1, 1], "constant")],
-)
-def test_fit_dose_response_unfittable(doses, responses, match):
-    good_doses = np.array([0.0, 0.1, 1.0, 3.0, 10.0, 30.0, 100.0])
-    adata = _assay(
-        pd.concat(
-            [
-                pd.DataFrame({"perturbation": "good", "dose": good_doses, "distance": good_doses / (10 + good_doses)}),
-                pd.DataFrame({"perturbation": "bad", "dose": doses, "distance": responses}),
-            ]
-        )
-    )
-    with pytest.warns(UserWarning, match=f"'bad'.*{match}"):
+def test_fit_dose_response_unfittable():
+    doses = np.array([0.0, 0.1, 1.0, 3.0, 10.0, 30.0, 100.0])
+    good = pd.DataFrame({"perturbation": "good", "dose": doses, "distance": doses / (10 + doses)})
+    adata = _assay(pd.concat([good, good.assign(perturbation="bad", distance=1.0)]))
+
+    with pytest.warns(UserWarning, match="'bad'.*constant"):
         pt.tl.PseudobulkSpace().fit_dose_response(adata)
     fits = _fits(adata)
-
     assert fits.loc["good", "hill_ec50"] == pytest.approx(10)
-    assert fits.loc["bad"].drop("hill_ec50_in_range").isna().all()
-    assert not fits.loc["bad", "hill_ec50_in_range"]
-    assert adata.obs.loc[adata.obs["perturbation"] == "bad", "hill_fitted"].isna().all()
-
-
-@pytest.mark.parametrize(
-    ("data", "match"),
-    [
-        (pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [0, 1, 2, 3]}), "Columns"),
-        (
-            pd.DataFrame(
-                {"perturbation": ["drug", "drug", None, "drug"], "dose": [0, 1, 2, 3], "distance": [0, 1, 2, 3]}
-            ),
-            "Perturbation labels",
-        ),
-        (pd.DataFrame({"perturbation": ["drug"] * 4, "dose": [-1, 1, 2, 3], "distance": [0, 1, 2, 3]}), "non-negative"),
-    ],
-)
-def test_fit_dose_response_validation(data, match):
-    with pytest.raises(ValueError, match=match):
-        pt.tl.PseudobulkSpace().fit_dose_response(_assay(data))
+    assert np.isnan(fits.loc["bad", "hill_ec50"])
+
+
+def test_fit_dose_response_negative_dose():
+    adata = _assay(pd.DataFrame({"perturbation": "drug", "dose": [-1, 1, 2, 3], "distance": [0, 1, 2, 3]}))
+    with pytest.raises(ValueError, match="non-negative"):
+        pt.tl.PseudobulkSpace().fit_dose_response(adata)

From ea236c97133926e650bdc4b4c66d727ff21698c1 Mon Sep 17 00:00:00 2001
From: Lukas Heumos <lukas.heumos@posteo.net>
Date: Wed, 23 Sep 2026 15:34:48 +0200
Subject: [PATCH 6/8] Trim dose-response docstrings

---
 .../_perturbation_space.py                    | 19 ++++---------------
 1 file changed, 4 insertions(+), 15 deletions(-)

diff --git a/src/pertpy/tools/_perturbation_space/_perturbation_space.py b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
index 7d5b0d5e..bfa84d0c 100644
--- a/src/pertpy/tools/_perturbation_space/_perturbation_space.py
+++ b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
@@ -675,9 +675,7 @@ def dose_response(
             kwargs: Passed to :meth:`~pertpy.tools.Distance.onesided_distances`.
 
         Returns:
-            AnnData with one observation per (perturbation, dose) group other than ``reference_key``, sorted by perturbation then dose.
-            ``X`` holds the group mean of the chosen representation.
-            `.obs` holds the ``distance`` and every `.obs` column that is constant within each group, including ``target_col`` and ``dose_col``.
+            AnnData with one observation per non-reference (perturbation, dose) group, holding the group mean in ``X`` and the ``distance`` in `.obs`.
 
         Examples:
             >>> import pertpy as pt
@@ -732,26 +730,17 @@ def fit_dose_response(
     ) -> None:
         """Fit a four-parameter Hill curve for each perturbation.
 
-        ``adata`` holds one observation per perturbation and dose or per replicate, such as the output of :meth:`dose_response` or assay measurements stored in `.obs`.
-        It does not perform biological or control normalization.
-        Perturbations whose curve cannot be fit are reported with a warning and NaN parameters.
+        Perturbations whose curve cannot be fit get a warning and NaN parameters.
 
         Args:
-            adata: AnnData with perturbation, dose and response columns in `.obs`.
+            adata: AnnData with one observation per dose or replicate, such as the output of :meth:`dose_response`.
             target_col: `.obs` column identifying the perturbation.
             dose_col: `.obs` column containing non-negative numeric doses.
             response_col: `.obs` column containing the scalar response to fit.
             key_added: Prefix of the `.obs` columns the results are written to.
 
         Returns:
-            Adds ``{key_added}_fitted`` with the fitted response of every observation to `.obs`.
-            Also adds the per-perturbation parameters, repeated across each perturbation's observations: the fitted response at zero dose (``_e0``), the asymptotic response (``_emax``), the Hill slope (``_slope``), ``_ec50``, its approximate standard error (``_ec50_se``), R-squared (``_r_squared``) and whether EC50 lies within the tested positive-dose range (``_ec50_in_range``).
-
-        Notes:
-            EC50 is the relative midpoint between the fitted ``e0`` and ``emax``.
-            :meth:`dose_response` does not return the reference group, so ``e0`` is extrapolated below the lowest tested dose unless ``adata`` contains zero doses.
-            The standard error uses a local linear approximation and is NaN, with a warning, if the parameters are not identifiable from the data.
-            A small standard error, high R-squared or an in-range EC50 does not establish that the doses capture both plateaus.
+            Adds the fitted response and the per-perturbation ``e0``, ``emax``, ``slope``, ``ec50``, ``ec50_se``, ``r_squared`` and ``ec50_in_range`` to `.obs`, prefixed with ``key_added``.
 
         Examples:
             >>> import pertpy as pt

From 4be6bef4a3b1ee64abe422f06c85b3df70628033 Mon Sep 17 00:00:00 2001
From: Lukas Heumos <lukas.heumos@posteo.net>
Date: Wed, 23 Sep 2026 15:36:39 +0200
Subject: [PATCH 7/8] Use pt.ds in the dose-response examples

---
 src/pertpy/tools/_perturbation_space/_perturbation_space.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pertpy/tools/_perturbation_space/_perturbation_space.py b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
index 253a28f2..fbcd6c97 100644
--- a/src/pertpy/tools/_perturbation_space/_perturbation_space.py
+++ b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
@@ -744,7 +744,7 @@ def fit_dose_response(
 
         Examples:
             >>> import pertpy as pt
-            >>> adata = pt.dt.srivatsan_2020_sciplex2()
+            >>> adata = pt.ds.srivatsan_2020_sciplex2()
             >>> ps = pt.tl.PseudobulkSpace()
             >>> dose_adata = ps.dose_response(adata, dose_col="dose_value", embedding_key="X_pca")
             >>> ps.fit_dose_response(dose_adata, dose_col="dose_value")

From 4e088d20796b519132a00c2a665138989c43cd47 Mon Sep 17 00:00:00 2001
From: Lukas Heumos <lukas.heumos@posteo.net>
Date: Wed, 23 Sep 2026 15:48:18 +0200
Subject: [PATCH 8/8] Fit genes, keep expression in X and add
 plot_dose_response

dose_response now stores the mean expression of adata.X (or layer_key) in
X so var matches the input genes; embedding_key only selects the distance
space. fit_dose_response takes a response that is an .obs column or a gene
and prefixes its .obs columns with it, so several responses can be fit
side by side. plot_dose_response replaces hand-written plotting code.
---
 .../docstring_previews/dose_response.png      | Bin 0 -> 40215 bytes
 docs/api/tools_index.md                       |   5 +-
 .../_perturbation_space.py                    | 112 +++++++++++++++---
 .../test_perturbation_space_extras.py         |  32 ++---
 4 files changed, 115 insertions(+), 34 deletions(-)
 create mode 100644 docs/_static/docstring_previews/dose_response.png

diff --git a/docs/_static/docstring_previews/dose_response.png b/docs/_static/docstring_previews/dose_response.png
new file mode 100644
index 0000000000000000000000000000000000000000..84a8e76b97fbd105bd21f2ae40bc1835e6ff5972
GIT binary patch
literal 40215
zcmbTebzD_z*EI|Zf+9!=2BjMTX^?J}4nYv4LlBTo=@Nqm5do1BDQOVt?gj;<JEa>o
z{m$h%_x(NZ^L>AP{(k!$0oi-4>zdb`V~#QA@_(!#g@5VhB@7G<d>LtRWef~#XABHX
zX50(#9T%;UM);qAgM_-n6Ki7!r)PFX7>}Mg*jQLQSeP2zax}8DH?_9nWaZ#s<zT+`
z!ok7DUVx3w^1uIr)!NR4&2IQS8s6oijkJb62F8_V$P4qEXqqVoM(k@D@dql-u`6RP
z?j*ktn>Q6dSXr=`o4z`SeU(}D%L?k;_ZC$`C2^G<pJ|`+yRyHTt`$!Q2)af{U}H*L
z&G7l=CrX;%^7WjLa@)A_m!FX-S|28;#m=WNy?@@bbL3gM^|^Ice%`;QcPH9D-q2Pc
zPC@*?f0BFhKT~@9{r6Ao4O(Z3|NfMb;s5)qWIr)q0wQwq_AHHpHX%9`3N`=4N9{OZ
z%!;0er>G^#exd6!rNGC<-n8(jC^fSqG#X7ju^nSwWm~9GWT6o@aOu*ePyF_D1KRfZ
zI^~W=HQty)8#U+Q-I@#&l*E(osh8M13#%$IY!1=3Yo!gM6pVD)o(~HTH+1wvp+wea
z*E2YcJMQiswnZ@c2M5muvQ_aV?N5|DvE}FIcSLih4Vq(PV_OX7y*fGCk8+V0Z=j~+
zw|f(%=NZkETk)YQPLT7(&s$Vp$F?8B!`p?dTi)MpXllCWb$YBkqU|_jbeWtdt6_kE
zg73bC1#`E^(TawSSOyiXb7L-FoC1%<;C&U9t8kA-YonzVn@!|9D?{Oe>lNWlIVE2k
z_F{Og=d+8)rw8*58^^a-izke?=R1*W$;x6L*o=*(V6E7ANaZ|<S@v?d=_A7v=h-l+
z!D4G_N<qij*H=U~&<mtQo2ochsGirq!yk27i0g<Gj56&_-Tadj+EZd{Naeoqpi6jP
z+w*u+s-XTV(WdF!2~}Q;!J0Q0N%-8+Hemx-s6{&7;gj<@{P}sK(#z9Rgq+7>vgXx!
zQZv0Y1<KRoZAB}%Y|gTHKWU>plWwZvBFpy7+^V*@D7`M&q*h9&(Fe-P#0#<Z1Wal<
z%$pV5=KXmuexZu2C%Tq%t96%-&QA6xEQY?PGz?riniF>0i~Y@=ymXOW5Q9h&7Y~ny
zho^tH!h9fSWibD>nx02QFgY(V3CWY>6K~A(=)>h)jXZr!-cg&^%jnaSql4){wvv_a
z?dm^XD23UJpeN5i&DLyO87e3}JzUmlKH6JtCsos~4Wp3~+F2AAXADh!6hpQcdbGRT
zJmxg<Xrjua^mwaXqu7d)Me_$c$3r5biEG!cH4hh>yPO^@xNNo3KMJfjE+Ad(l@LGv
z1)<Wiw>qLx>B?<uYnwCk`U1Ye)@<wS{CuSIc=>vj94xxa;bL0pZok&S(Rwxi{;2)f
z%^2RC+Y#ZC8?q5f)t(;S#DX&|Ez%b*U$22+3#3R8cH?>V>Xo*ACwF6)jEoHC7YOP`
zU2+lk5=kCEzmH!XSBER${hphcY>(P^le(;KOw|jW9dGX}4}_MK@V&cyqs}D3W8T%H
z^E0=$O)a*3yijS23j+f~Lu2FS<m)T0bDVyD&QE*^D6exEV6jyn-%C$VH!9L8G$Vz8
zKR({>ZqzlOs1jzY+#$w*bx+Z&_QVi6S}Dj;CNz44$C`d&9bz6%iczf(c9u?gq<K#1
zyAVo2yOZ640thY6gAx=fRaQo()*v4~{@PVH%eprgcmDiJ8r)2h3cY92MbQ*=qrN4K
zCS5I;YQ1tlUTbEni92~q*+(to4czSJlrOd2TIoWhz><F7Hzj7RqM$fIr?L>z{()?b
z0^|I>g#<4V{X4JX6io65&Y$DDuPdI+1;3Z_tp1%lTtto`4gtj%xhU3(?JmKfo7@4#
z<F0}6!fwa~Zg+5YV>I0@YG<pCZKiVTcxk1eKmfaBcy@Y{(*NWdQTj@}$DRhH6z6Gw
z`Zf!v*%l(?j5<0x<`4Ae>aEKX7F&&d=ACdS@&1)0Ro)g8OeI{_w-2jgZDqv@XG8yB
zYfk83Hhk2soka-t+hX{Kziy?g>C+fL6#?YR71F|o3QVHVt}<TBUdNlw6e8}QMy#vj
z-oGc545jLXFj4gMtQwjA&zc<fe4v%nP9oI&@#2Al0~h{HF61MyE?kIo?s|V)433EX
z?=BsP*ya9gLs&1`B5K%I$f>lqZ@=QP9JYZazGpeiHPE6<M8u#4Zz1Tk)-Yv*vGG&v
zO1}MqLL{4B7wpDje>V9Q;T^(@7ca`k@`bIA6rbPN=~cAfY`kV)kI$QD*i5iib!3H{
zc#%`I;eZ7`5s_DKg^T@K*)pqf$EVt+VDj$6lMJ<71-KGSLmsUur-hc57Eagd-0IWM
ziXsPZUY*D3$<-4b_d4_F>Fo^++9M@OzuXi`t(c||cb(hpYMj8!Y1sVvq)-uq_3_G-
z`0cgw^?Q(AzkdCSab3`n8P2caRNY(Jq6zrieC%mML&L?f^5T)4+p2-Y71L2=PU~wo
z-)^+hC$PgE!%1l@9MILWunl-vTU#4+ox^8;!i(Db2B-1QuM&8yRp@Er5v}Z}WX$R{
za1QFY2eifP?Cz1!$$yauB)fc_{Y#6qKCBrlJ3FCo{i`vnW~*XOuTyt8tCNR4ntw7f
zbLy;$w>x<!FeLF?O6?bjuU>6<dri0HhiR|#VoKDX{_L0=)MDZd=16sTC?@8Dp7QM^
zNcjXOCbi{T(RR0J!qV@|O!?Ok8ii)RQ49?J$WlSL<U6ivT2-GoqG)JJ7gJ*bn5zgj
z2K2o04H_>sg-|XFJpFX-ZuXmp1;(A){DI-r>lNDy{;S=yH<b%t5R{@*o?O#d9?YkK
z43OlsK3)MSEk(hST-cTCNt(hOM3~9m%Fw9OgoiqF<h|VA*;>46%C%EA>~4D&NC?6@
zJXu=$j6$7r+h0pjh!cQd`w$=BHC69>u-Qy)I#E@bxf9&G=n%Z|NKVc^ZY1xS80=zq
ziS10uqz_)nLcE*&!DSR`Tm>owmB)?}^24pU_D+H{zbDpRuZY%L=O`fg;srL=K@z{t
zt`ES-#RWMnLYoaF>e>CM<^77IwQ`Iga^4*2Q++s`D<j3B5R$!6)=}%!#Kg)52QA^%
z>+vh{|Lnc8k;B^P4R1ZqgOm;)%aNkFNG+?|G&HYctjZRD<UA|XEQQ*Nl!^qeQ+^1e
zqzA7rn2;9x`H9NP$<4qwr@ptdp8p~ld{dp7>~Mcw&x~?xW5Z_P6f&Fhe6)#QP|$3s
z=;<Vc<}?(Yec2Th%D4>f)OjYD_qwR3yFcF7-`|wZdjrfM<&O^wgK*hb*XD`IZkVCX
zR=tf<9>T)L-b^?<_6-#|^lNAkAK0v@5E41w5Q}%;l6gl!sku@C3(N^|3*TtvtSUmv
zL;`VzveDJ?i%iEoTA#hnjwXbTH|lp_Avy&0@vgDX*5R=+bG*`0pRK6)RgFTa{_a&2
zg`Js&b8{i63PQ3AG7NRB)3{3t-K84z2SKRmK0ZF`1;+UBJ^JD?M{-I5he^nk+x=36
zL>v+}HtgL(+m~siLY@O4Y1FkI{=tCJ+S)1+ayD^xqUXPk-85~ALK%GGHYbC;Npt5;
z@}T*_!GXkS*V#%z_h7Y`DDSv4O=k=*TmR`+JDUmQdPpbv9qYwLn@E;_g0GuS{Jo$+
zr?9O}^Tzebobt8lZqc(@SXjn_u?*GBKR*-C+s}8{(mJd<3ihTcnru!tUcSjC9eCpB
z7i-WS$ujD?oLzA+A0xCkXn2`|@8|9?Qita|V@g+xMgXV0Uq<KXlsoVIlF|vx0<;kF
zJ5!C*dRz#Q<tGlq3q(Xj0LkPt3j?JiY@j5NT)U?BJisOqKu%A&ll8dgk>!BTeFS*H
zTC(i<dO=mk#mAq6-QSQLM4=+;Ad4$3bj4|uIi#DD)#V#DBeg<WL4ozd6X*RgC)2UA
zm;U~0-d-omdMm?)fwXdwv9?0vCGS)`Wj6x!G`~OhHz#eg?-q`_fPb^3YEaknV9qFB
ztJ3wQfuUibJT_}9)XE!>!$XC(a9=L>&&AqzB|<I}+W9Fq`f@;rHsl*2QMwtDcHtNJ
z;q1To`2`b7q%oWmq*|sWc$P_zK*js+eqaX`KpgJi#drn(?$Ym>uYs&6X67Kc*?UfF
z+T-rqkLVQR<-6mBF_YE~j}Nwo%3d-#joSL&pnf9WkjYGj$5yf7lPP(So~jTh`}i^8
zqZpnPD8fNEINw3==QI?--U-2W-DYGA7#tiVZWr<MQ=#45+;mtUj{|%c@4CD6D?^oO
zsXx1YXR+61DLw6zkjp&`tzs*qq3=c#7noz0;F`V~cYe;fS(esLK0tV;Z<}%V13G$l
zdElB0)%EMr;dGB=SvvDPPh8<$<dP(VWj(rARxH$hZD|&oMW7GiPQPgs&;{S*RyMsr
zl)gnkDbN5Bpxyo2$46UcIH0Ai4aLF|4EIb-Mm8@{^>X<soNCFDBFn(I0*@m{h{kr7
z&08qt=#RTPyO4n)tCNtDswU?$|9~H#J9qBIXbA@f6!$NXnxUb&Y8n+DUJvil(%zng
zd`nW&pm^i}ZWq!-GJFGrma4tIeP+G-EY)_VsW1IWTFa5uct!M0ZnK}sG9<;v#e_tg
zMdxt{@~kIBmf+_rZWhMIbWoU?^Unk|pa{W{f+}*s+w1fIC-b)~gneKk%iX*By{V6K
zs?hSew^c)lD<Go({S)y~mg~QLm!FC<1loVU{Qv(IrZLW>0)g*q0OgM}pRzEM1rAO5
zNs}%(2()%}r79&!s%M{ii9Yhk^FtGRt-t=g5a7qAl)!|V7N*@T^^BMtD=@jDq@-l*
z3+b-`K%;3-vj5OC=h#DiqD_g@VVulvQX;)%5jF7!H-$C(6UJIC_>>y7cVc29ofHi{
zETJJR4hGV)efe@LoQ8O_0{R{R1KJ+p)laCIQB@5SNBX;^sA;*rg#dbngj>ogG9Ed4
zw9a^GvTxqK%jn(ZfWpE@K~RG}OCr$7_>gfgSPr>dl-KG%m#dt>6#vEDN4oW&3t*~Z
zVrPF3N4jk9Kb>n^68!k!<0NTro;vdzU*8zYcmFIcI_$>}F6bYcpm8jn2@x!TGAn0M
zY%}$cj73u%IwhndeT{Qv&o0?6^`jyoL&6!0H3LfdfunD@-iIrBHe*(vG9un%Q(S^~
zh|(MS(%@}#w+1GZ5<K>9z;T`jVDApiQF3zf$3ZIq(4(G5t4M3<yrGXkDC6$<)wvfE
zjBH*jV#lkTrdYYoZ^x`mAp2nau%4ib=`^;wuRO=1k<gah5S(m0UlJnm?k&5m&E@lk
zG>?Mh#TjXK8BnO{%;N(ml~5|-K{)x)c;rZwDC|Pb*ajf>-F4RtS|6MmKm_%PV!jpE
z(_``HNIK?1J^VyH>4K3p8Ci7)gowNbR~-xL!x*hv4^^7FG@m4M;z+l5JX^`)X4XhB
zQ+|$8zLBm7$yaY~y*jC`zP>%XRlwuWfwgp&Sk12Z@#=7)@b&1D@U4sUo6RB()S@JJ
zYngNn3N;~}a%XE5H<|qw732sVn|iK1Eh+oiUCG%(dAdnfR9H;BLd*lA9nE+Wau2C?
z4vRf+fi2*KQnywzLr`_NlzFeawX>6!i7B&1tOK@_|2Y53W^=cSwCIg#QcO%t1Q`<1
zND>l$X}e_*pZqb%z-dWg9!r_iURt(wL-#RC+4P@FPnl49hx9N|q{-NHJ5<$lB%ph*
z_q{@Y`*vcVJuxvc0>_0RLMyK4$I2eO39xl~q=D?uCMRi+;yXVThQ*;2r+Jee1@R1?
zf0lnou<bG6B%jx>qj|04&3Doo4W2<0DgTKbGq6#4B^Hf#24ub98XpmnawiwfrP%jz
z_wNEdxeXf-#z<T7)K6=~U;1~d2!yPAgs0Sv;bLn=2M33+^b_b}ZwU&<04RxYR)E?L
zCrg3~Rn_){+vw!qHMpp|8E|~;uIn;KioCYAw>t@>GV)k|z4X6~mx%UoeT_Y%U=sDG
zZ2YpEI3vS<FWLrx?>0c#@DCr(UB`3m1Y!Xa?_(%q;Yr~Ek3|w4;<Ff*|3*scL&!&J
zKzH2Z<ZSKeXjsg0Z_H@)8q6Z2BM8B~+n;_-B!Q3bLQLTg)@r-f?rsxkaE?!p)=AA0
zA|smsnm;r%o73i3X&A}SOn{6RA*G6@EHr8-g96&{<-uz{yIEZ8>Jy%S`|lhcICb;!
zZfl7iM|)6oe1QYl1mur|^yllPu>^A@vJ0)Bzn>2JqGwDPpW7b(Qj~i7l<e8FXU~m|
z%`UwSwdJG8d873FdGgb$tGzB?C4Z&6Cn~yb)Z)0zMX3GhO}jDyC?VRC8i*I9E<(S9
zwD&gvU|Uxm82i6$m7T^*P-|cw|MO&a2kn+HyPlvZKH8j6y>D!6yfqUlN+w8k<AzLg
z2xTC}2qZ~zUaLDoLa{}bBL)RSjP8`XrSr<LfOrF|^HrtGZI=5>etVZqC(uDTQ=t>~
zOh+b+>54rl^J&tjnGsFF-C3MOn{Qs_hwX&=-(sMBy|Y45XMUrXopC0vtcO{*;!|Vb
zHO}WPgx|k^Z?iapMcC>RT;G6d)sv$m2u-abw6uPKfrcHQICKNlleaQ1A?Ha&;q+LL
z5wHRvWz}z%P<CJReGBuYAy!28cKqAFdrJqbUUJ)a4HgG$jqVk8Arv=vgqC$>1awfF
zfd8UmCgT8f!Onbz799ggUvMSRE+%5-N@E|`|J}|QPbq(NTHWXGrqTSg{ZoY)M9*$J
zm_8+^aGd$l>!1118hhW8Z*lw%ZtSZjC;`O7YOFk8+>)_TnKy9k2wdkM{;qt|tp<U7
z$+y1taCaG##x|LO6cs;X!*`wE?hf=7SAd@(6)=G$1LrESB_C?{*1<wTl&b=C?WLuq
zISX^NVgF=&g($swl1p~C&96{&dK=!bS()^gQlOiqetT@#lyoZ%-LHf}fWsq&NkYjH
zWW+^Q<DzV&fk#kN*N}8^FutnoN2uo932!e`l#-8H#n&%iyn*t&><uW4@aom)cH*eh
z<sp-VMrb5GkI|P3*|dLru^Rj8M7BJT+YN*QlBW3hN(Nn!F8niBD+@>9-(vGIxd$$4
zYqEw+dIHyt_kPpUmeGsR<&zt_#X82f+&<m>)S7%8BR$t3eS&-M({=Ri^Ol*Tjx3ux
z$iah-2FQb{;cG<QgB32Eio&~BffJxQ+xI%dsTr?y<9FV;)gHx`2F+m~^j9tXWPxwJ
zy)j^op3b7-qA+nUUuNFa1W*V4Sl6x)6I0-~Z+D$MDIl8@5v5NvJ`%Zoc6~+t;)vHf
zOVRqZiJEvd*CTJU$4QtfnfSk#Vx8Z-dGHt?7xz_-S(=4dP&PMV7N&lo54ZO#GV+uB
zupjhMJ4p(4=GeU{rZ^UuXXn~Cv9T3TtxfXpPrcA9UOT#^7tknJ5jiPrGl{`z()E0<
zEdrre0Q~2v$5guQZ_US8e~68h2W|tpnzVo~NB%W6_kov`ET|{kbc17-1f8sW3{U6F
z6@YgrHnz~S)BWm-{c*P+X1_BAU7%4m++Q1$2w5GkjQRHMn{wKHY#f{@zVB*(>FWj^
zcYy@`Zr{$pXZ7dMZVa>342P~Zci$x7OmRFu)BZ+@{j8b#_i`CGgW35qUWbtp8QPf$
zto5%Ffe2{<<zWFB(XVQ`x{bptW93n>u*AT3|5nrML|O^ByaQl%e`tP1la|*nsxff~
zSVNC33|qoxXEECnmYkBZ3ICuh1j09vMXTs1v?&<Sqjq1@bARyg;kl=oq-njab#L(&
zC#tKjC0!%hydo7!Wz-eRpQT;;0RiD69*0-Wx9hKn=2;Bg3!xI0)z+qf%W5_*h$hU2
z0E;cEKbjmUezxi>(azFmhV@H0QRB*tQf^9Uum!2tTm+Ui^gA*N3UqhwyaD#)<KR3T
z9(wo!`q*>l&m)NsxWhYl?rb!XTRYRp6NPGh7m1MZCzt8$u%!&VjQ4DM^7@4)$<NoL
zRk*hox))yIP}R7&7xrExp<{>E%&KJS!n-TfT|U@&z;yd&_D4Ye=>Wo`6I!9togO(T
zjZ|>=Op%k&uHFUi@~dKk2-l0B-t_U#KJf420BY-U(!g%H5B<7~J|LBL9Z<)`|6(lm
zg7=)kyZcFev<{6|b!9drB_k6`01@O~zxUiPv)}OkeGH(=jK(Wm;t3QhZi$K}04;LO
z8!}j7-{<&Ki+{k(>+O#-keh<y3g4yE%zT_B4dp}HuED%N+HrJ`K3|P~^y`P|6VH=f
zHC@+5TBIaEC{~{yGQn=`I07FUp$Ro&sKl10(0qV(mu_gbH5_nF;qpN+@3=HzdLRrH
zAyp+8KO}_{b=4h}FzHZ7y*;@sHT4H`Mm>q|k1h3&QVI)Iw&vSajV5FoBY@VR7h|!@
z0?O55qzK=)UdaHzWpq^A&fZ?9yWsnGA*!_48xR*&9!J_os=<GOQ<Fif4hFsC!8&Zs
zA~i{(9saaqsPLY8z9GVa4(2_R2mu}{0QwbICT3fb2oIu77NnIg?rd$<0c8glNUWjq
z0S_iKNK0tE-tCK01T`kb^^w*=3;4tLt7UA~Qf_85iPX~*R{+lm1P9Q8y0gx3zI1Zm
zzIk(1kPS#lnQ%H6CV755Kwig36K9UAYTPmmHzeBb#kO4pocjLbhX?zSyA{_y`Kg#O
z9-*eGMNW2}+`D&AHtV;R^m%e|@#Ktq+zEMZ-0sKLkOqtN|Ln4BJF{V7x$qK%*?p|w
zDvQ!!CWA|b#HJrM18ndu$f*sptoPmBg`%UQbDE(SsU*sY#e~iZ^b}LregCBV6CGE_
z#nv%SQ+kxL`rce$jMQ5ry~r#@ZiahA{FG3lXZy9RjYdmsfk=-SW`!Hc%VQs)VPtFz
zRDyc>sriR{ot|FJ&MQf~Vk@~4OMdO+_aDlq!`<n69f1emXr}hM%D)HyMOg8VuKxp$
z-K}`q%FEF8BhF97axTIHGifDoOlenk+E0bf>+FR0<ZyQ;fGKAZsKvVl*%v?s@N@o^
zust;sQxlJq&Gvev%J&_Iu@tQ`V}i^^2nGKKsJ;1sj%vQOmh0D^W8s>jry#(I;juJ;
zNP}jmpnkLLUwHRpNr)`S^taB}WmisSqV!~iw3R?VbSiY6zK}D10$B^CT4~Mi4!#8;
zQCLRX^1&O?h59tZ8dG09$(oHTx9@ap50LWX!O5~3EBjse1u0_DJeD^ho@P>_bnOZu
zU0hbvw#5c+_p9|pmBY$l80a+i3vo`aG~%0iK0fDwxoNI&*}i-CF1sK_Zh;1iJcY&R
zsvvi(7TO)s220Q#6zSLRY^aElNEWR>W97Vn>Wt<)qE)}ULz<A5aDiAbq;(e|vw7I8
zWHYCdr&{^)P4d=eTXm1X%O!S$23(8*vM&&HIUX^cO;F$<&d&p-B`zgp52B76jnL+~
zz{cqeQj_$rneqfd!OxdkN(}H&5-%RyA>x08L|{3hRy|P^#;$?Z5DnUyI`HiXTOA^}
zhM^CN9|fO{H<V1P@>QOJQqTz(hRt)iX|4XYJB99bwU+Gr7m}F8VjkEWT>be-J^c<;
zGfEx{^|v6GjYPuBl@GO9gcL2nj|KJYV>N@)ht1;1QE+u6m?{kVetGR?sksxM@ZLCT
zTJTbPRkG3iM|Ix;AKPtdBR%noTwpF8mizH*3bFk{LL65`fgGa{cKtC}3TJQI-5Wv8
zL26qU5BMeR_WyYLS>ULhD7NaPt8jH`H?&(k#T(P=w?8^!)kx5^aIi_1#0f;F)O}|&
zG|8?#C-I>I7Qr9~kQFGUpYDxld!6j+21r7-Wt=@qd$#WSO_w?+2JZoPM~L+7k3WuY
z{HWC-|K}I<U)6G8VBz$IzvVsm$dC#1_e;AP>EVydr1|bPiunyBQwKvSc%O_QYESel
zV3a~Y9{f0D1xyy`R?$FUx-#5Lzgl#H#Wc>om+MUVkUW|^RK@jwwa`jE6>|GMxc_uA
z#$B-&YrZKZ0czpmy(Mbk#25XsB)`FQSRv?_dpI#~k0{;4qyzEeptSZX1YA)^+sAS2
zR<v6_9usZjI|0y&uv;LgemDEgXx8`bfws1`zh#^PJwV{;v(xpnKRqdOVY{yyid`U0
zGg>>dhi|OwehDP8ptSbkiHh5GllyJlzhVXwFH@`nz1!Z3?P}oBos{QS$xGet(ZBy3
zEv78{PD)nH&#zQS^)JtG+!9}ycfd5bl^r&@OTul|HyN%da+{9s9dzGOLym198LnQv
z$`09VZ*6Qj4x^D}G*0nuymxu9^s)_|q)tmsfu~b?+a(jVgu~n)Hjah#CYp@UJrfZU
z)&OC+Fnn!arkho`)QxUcT>O&hlqAHZXufM|e(;TH-XteXR-o-K6#@gvjfj5LCwt^G
z&7u1?D$`bL9rPbRd)vR^a2%IRz3}c^Uvw-{fr%PZFP$!M>WH4%Nl?%BcHk02^)J&O
zBt)1jqj!N*;I%85>ON@lY*tErXYdZVI?15xwV(#eJb&I;RS){;hnSdUq@$fn-$*_;
z8liYs5Yt4-MU0=xxI3pGWsZu_hK5u8-P^~*xzu=)%BgZc&PLM#zxLDxLnM1XKS_J<
zbbYDN$5^-{q-hx>5xE2rkCLr6CMI_UKt#~V<w&Whd|WtpHJ3jQE4F3vo||JKC`6>d
zX3fq%2c;fbJY(M(HHdr?oMw!xa$96b4=U+&OCtGUTZ|@FQ_BVGEG#DYG^@wz@we9|
z^wTxHNU}5w1NiOde}duxylDrN@gYbiBgNJx(6T3|rA0aG{?8giQ5xH!bAVMqVlzF_
z5wDXC6ph4#2bdt+81!dpz@C1h44JrboxHU#%vSsj>FQln-9sVG1}*X9P(<wYLx=TE
z1Y`GSz503TT9C9^24pIS=dm;N===_U2muD6L`5fdN3DwCbN|!b8&Ej`MWayT@5KrC
zZ1xZSzsY0qXE;|SR-vW!4b9P%v;c)5u8UNV=2I8n&~@r}J`^^yuPZmga6--*YB_f(
zVekd-I|#@QrhS}ATZz)Y^|Si(+dEqJJ(PvG#m2p~FUX75V>fitsk=_jvPbugt=w}j
zm7FQA>nb)#nLdAwoGFdn_<Uh&_P%e{R?KuSetpFbr&naPvlz(H0W?c0$XeB|bMU4P
zOIvfEBbGLKiuwBMWR>y9UE)(KOOPz|6JJQylNV%00K){wcoxzuyI@ccf86`R=92PC
zYZ40iS3y_#O2(v@hPOLH@*Y~;cPfhKv+4WsSUmYdPYXiwi`uT-I^OippKu8(&3Bg8
z)Vv9;16I)6UsVV5<OI|r&<<ECWX|WHd0@SeFmUTPSh?1mH(BG-lk6D5d)4V~td7T$
z!EYj6^(@Vcy6K-|75MV7i`;Rvs(r`IgT$>ge9RtnUU;0%iZ$@&%S`CxXF`?p=Q3Wg
zU|2k8x_VdZOTJ!7;RWyC7s#0OIf}^IWxDxICzKDX%$HfMPUq&@^Zlhu6@n~U7;Gsj
z|L(R5@$&rrS`OUMSH{;qQwB=`p$3c!R*lF07)}k)Z$NWMr!QQl68Z$#9VuH6&hJ;X
z_tjjbvlUD~_^s1i8Ri|H@l4<1&WRxkn5$dDg1#G?O%n?7RkmNSa#dKO%cnwJHkzCB
z)^oq+YFIQyCA_mtw8EQFO5^V{!{xqtU!VKQLo%XGx`H$^f7nHGJ_`CaP6w));=oL9
zt5Nh1ITmfV)k3NON%!%vle{U;3&UHBUhVdG=D7FG&M9ON#Qc7UAUD@XNbZ9ujRrFV
z2ADI>1F<-cufHydK2TK`J#(hi4pJ(7R4$fT^>gQ;pPwqR+OaI03-GBVd!8JU=BTPZ
zep~}EAziOp1QYX+R4K<tJBgsU*S4+R#b9D&H+f87-|&;G(be!OIa=<c2Ds?LbDsuz
zKfDUIpPJ;^K6mx(Nm}U;5uBi?><KZ+dsx18*s?EV=jMyt>&k|j%J_o?^aGcqop_t>
z8t|O#vg%Gb2n<zul$%eOhphOn)75^yVoR&>*^&3@d*Ss<-=2=R<@-6ePI#Sl0(sK~
z4aXyW{ae7%kHX#vo$U0c(77T;8O!3aCg@?l>6G)pz04w@t#W^yhMBpeJ)s6?^$dr{
z12f&LrmFLOoEh_Xb(tmm1pim}z0%$$*YG4iiLoL0MQcq`M5~}~;OJn|wf5uvEZ-uh
z2Jf>1m8v@x69?9*|J+jDukzI#PT4kp6$7hQdOnJck;{hFGcR~0=%xh`q8b+$S06x5
zjz=zf#t^-9yM_N$YQ{yaMq*!VCS<&{%1;K-{GmjC`p0aNrJnBt>;5ZCgA));5pnlH
zx7TSs2wp9sXUAf|aF#o+el_WiZ)k4b1`O<SK!ScRCOEb1jAqI|7p#~aZ7Hl|UjCs!
zoQKS66mx~oGuib<w#USWV*>d$@&_5);wO{Rau-NmFWh6d$@T{j;8*f1k#b{yvpnCk
ztmLx_3*1HWu|9Ju6~(1|wQ>#tq6vss=eNJjQb5WbFI7L~P(Q}aa+z?ZvK7(ZLvgLn
z$1MkbTIsv92AEhn^J(@Vv|s59cj)PR3(eIKx4?^@<crsIoo@rqqYD0=mU2o^ko&6M
zEm+*LUG9sE;)|ab%y{3ve~g=h{faw;?Tb2pc!_G`huw5&@lw@mgSbp|n&wAa`<~fq
zs)&h;AJ<#=2H#WqRr#g1)L4qfu7C7{qPmCs&Kh<-cZSaK2dFH1A;vW;&!8#Mf2Dgx
zELSRMHc*-rRH30#dwS^bXTWjt94H2);bAR2iE|on#pfNZ87;M61iiER>{Jj8D+sC`
z#~Y3U;we<}dh6J`e4AY^Tx&x|6}xvU1oVJiBGvJS{5h&zfddD4_`()E80n*1_WGST
zO`cI&oxu_fY83n!F;!f9s__qwKmJ?sfO*4gYPp~Jd+9WgVjPE>;BBSX!$JXei!zRB
zh2%>>7gGa^LI=DKBo@R~^b6!O<&i7b*}V~G3b>_|b9IHm!d7$sick%hzvP7Jd3mFt
z*&NsgL8jRA!V?c}s_*;v$0FGenQ)aTr`wo;hwST!Rn~JlmmbXz;K$jcPz$nh@i3a7
zLKq`4cxu!unlji=EB2+xRPNg*gJq&?*8?*H(YLQW9}*8>0;vm>V_+mf!}Hlv<kOMj
z)SsV&&PtLm=uj4bYY$ZB!6HlY8#itsBFXW7_1PSF&R}hkYlJSV1G<GlaPfij=?)#8
zUAOQ)^j!bJ2G3l4Wv$Od^E^fLUu8}_u^YIOY!VeeDje<zZ!5yIsAD){%QeUqFLzw8
zu|U$##$L2=#ok&UB_fWnjDd|0L1jtTDrO~f8c6@KP*0zve`fscH?AbxyLTLGmG7i!
z`;(-uAfN{}Sk98uyq_8FMhtpQOuv8tq-y{pmdM>*gxm*CTqhJUI0-t~XT_RXbW`F6
z&vXP!jH3P__92HQ$I0X!@#bMaXzr1-_gCRDiBN1AK0o4lZ*N3*?sQVz{f*V4zZzl*
z;wRq7V)o~K+Hg4wlh^&+u9nNZ04m}QSXYdhG@V?*R09s7`o9!jP|#(yhrC=|AA$cM
z<zJYreGQ6Pq-)&dq&{d7V`RmOErOS-yE*pCCrquFr+J`9B)&oRP}5k!|D`j_eQ|YM
zq-zx9+I)gI0ugH>&~ic4qVeEjh#H(PU5KaQ=a2sL=LNVw^o$uOiNfqowo>zj&u)IZ
zVP(aPvT}HO@CHBOyn^4-g~F5PO^MPMxqQ?ZizgOpGn_b5m8&yd1yD6cQ97`4oM`X}
zA~rklPTckP_eZ=dIXqwvpcfDr-gSq2LJR@!K+qrW4(Pdq^j`rSj}dFfd{j}c_i+;8
zFxKaNF2S3&Br(718fu$5LY_3cU;PXC@Ur~1I^u<6;HrB;r$*#*S6A139!MM@xxEJ6
z@f41%a*`xIYuO(ZkM+34!v@s!H3BQ&YO)(gV-De(4MrzvO+0&M_A`W=?7@=}ZA^mr
zWTCu!Zp~sb?VZ6=4N?F#f|1?^>INmm4ge4Jd;@Gal`jh~;^QNFdd>o5TRu1+9kKl3
zZ{A=*6e6yZl9RoWPBiYk%xTwPQ5n{#j)J5-H#^S<=$k%4Pt)RK4otZ`1w}oyizG=T
zUI@32aA_EH3bE=2G010vu)Yay5QBX6KNHoW2(|}8ni=S_D&BfhBCoI=Q5%6T64v;x
zL<}Fp3_h~tt0%5TK#rSMLa1Yas2d=uoo5Fz?cBMpt%UgAi+VyB{Lw$Zo!@c)iLVrZ
z_5}iL11vF7L-RnxA}>eb;vow}tW;%7bjz4EYAQ0)gT)EGew|$A>=W(><$KcIUAHSD
zl_DR|I5ozZ0Ti~1S+6<@fjzjEPQwR=FOXH|TYri-=*>N^U7F$bEt{}}kfE|^I(YL)
z@A}C6wR8)Av@n7)fY|`^W2pIA)(Jl~j~90T{p(wiPyI;43Ze>e_AKwMZ5018a^#x6
z&6D|djb%BS@|}vUpySFHP_HFIpwsdjD|d=gm!iD&KN4DKGXpB3Vfe~zHu9z7{zqn;
zawf5g=p%E(5zPXR#xn*c=za=@5C#!Gc9!yON+|hJB^4ZflN)h<xbF{r^DS+tR5A_W
zlBbKI1#~R097IL@&QLucq$Er0p(ZOWt>2R@qwB9E&Pe2jndLq*xb8`>Uf2<WPT*{#
zrx1LRm5ce)^@}fxWRp@PYO=GvokmzV4%9C6#L|lDr%vM};6T<~3iKx4)teIbUv(W-
z^K>H$3y_|8Vezt8T3o2bARiGGA!R?FMSQX8aJ%UfUNbW@H+ZdLmt8>ZQc}{y*G_I5
z`G^}W_{G(uWKy#riD=UgKHuViAlA`vLo#kY9#^g*4c)O?30WTaS%nub?n2_DfusR&
z)b7vE+XZ4IdUk-bj|1{27L;&*KEfcf4d1B_vP8SgH{C{gvazy;sAfJT<=+Duni!lm
zOZ?!(gG58-;^()~9@A8`eH}N;7#6DP`L=<f+gb|z+BB_5_@oDTYSc^ZXyF_pCK|v2
zn;@}qKC8om&*Zh?qkujlEHcv3k=xnpvjKA!8P)k~m!0&0=y*%xT+i*JCN5+S<->jV
zw;HhE#NMFT1p`|R$;8kOwCMB|c^{q38nVlA)*OO2i9-2T-F<|qi`)%Yich0VUXyGR
z!8JURkx2v#6_?fMb#L(Kg1b;bRyLqD%wG^W7l?=kug8xc+ZbEh&0-uMiW*i=c-6W_
z7P?D_UmDteb~e#eCAkq$GDxWaN`Fn=%kl%-0}Mys-j(<6<9c2v65xCM)0a--dxd%)
zfXo*VV1%|>XmWT4blq64h@fvlRZcZ(k3<qt`C3^ZnAutzS|pNt4h~ce-D<fXbB>Ox
z%ST2>w}J~t-ptI*KUXkLL0MTDVYvAD`E^4Sy82BQNF1}cF1%CC#d##)CUT#MdXkgd
z8IGG3Sg!;iJxD+w`8U6oIWSPWZ@d~x^d|?&2!xTjVNse)6KF|*E?T&3W9WcVo-?w*
zcB8`TXJ2v)#Fgt$yT@QllO9cfN*MA#Nq%ed{;V6I$U5fFfPsl+>p0pZ{@F2l6Y+3^
zbhOm>YdQpFceCfr=aC);SxMefj@)<45gx>P3*PP_z$3X8TbCfoJB9JUC1HYj^bU0H
zV37smVm)WKP`!@Rn2~ayU73Pl!p1(ap<5RQ{^fASdK2P`C~doo;P_<%ot_Rf)tw!O
zwaR^6#L~CK4@w`94`uNtMn*0)D*Q-WB0+Y6-xAy9{Z3B)#40Oa@11-<X_-Aso+5lr
z^IcP-pxL!S4RgMi6;_(Y1i*$pTWaB>b&*p7nH$JOmBb<-^i#o{6^}11e%z~=dwpd4
zg4k>ONBE8MU?EF|9Yx%iT;P0Z;3;=np9egZ0^<i5VDY0zdMR*f#8U~oI&*09KS_&I
z94s3hELvI!a~5|G>pmuRoAB7@0LN*O)mYowSosj(oxp2stlK%7{{k7`r^I`rIq8Xd
zv|Dn<lq=k*#MQ?I^~6M8;e{iE{r!J*h|23*!T!<sq1=0U5dOkAa0_J+`FdosTQTS)
z-h{)i>58`NV#=uFh$W)^_3?npmov+XI6lNGU3bx-^HdThI*&~`43>`%hOQaGQPrL0
zyf_Xlh)p|7amj}tV6mGo8nH5jtx1U|p{7<eHcnkejK$V(nRJ2nW4tC@9JD;H)9S(W
zu!QwYRNLYCJvkmqc0n$pP2RT`$A3y9XZ}2tg@BJnd=VtlogKFsXzSe``1`e76s^9X
zAo!yK>*h;#!OyZLprSl~{`@D{1eHe+8Yv<|b@q0Ic-z^Hk#`GXofqbtP}34{a`Ex-
zW&qJhgcv=4PDn=f^nausYsR~G5qfP3Ug?jOcIwg-(@tMSN0PL_Myr3Xnw>sCp-k>9
zJ12a_%EetABGqzF{R+HUL-6<gx?znb&u*N~;1uFEuqgoUKpghak(cy2er8$GQeaBR
z+jTP>d6?vsVzwELwGQ%?%W`({ClL<;;MD1nSkA3FWQ5q{5#+Kx4{|klI}lxnz&EL*
z=<`iRlx|dE^B|CS51q3Ki;GBzh<G&*=mzZ*Wv;zQJ<Lz5w<@+rjP~G56FED!2k}kL
zGIOY|Dy8lKD~+V+RDY#-WhM@Rgx$>FfddZx^-YOpzALUBhC`sEpM?E1fwXAsOJ}8D
zFoYhcKAjt8YhRSdUmgD*FIaySIDXZyiTS0rerWH1G=keKoM{a)*<`nXuZ8baQhxT6
z{_VJOE9<lI%D($!JjH|y<zrn)GXv6U8`ykJp`J&tzJC4sA_2jkvV(H`r3I3EMGxl{
z0KoRH=>{k?@S^2?2FwaU!YxK(=o4jtvp3W<=zvJgGSd@L)2akG&ojkeOW*zr*Epal
zyMT99W{WJ97z`t|#;OJnz?}GD{B=b=;W$}WW#7sRMtrykaItFuutKzKkU)%m@1&*I
z9lXhI54KcM!_Aob^-0`4xciQ?#3fU+Z+_0df1TIw=G-%v{6rq;hkz4=g_D;ZRe%2U
z5ujV31=-nodK&%!JN;#}K5vV~uzG<ncwBf5d|6pq#RTlcrQl+m1GnW@XrBMN^$_C&
zz^kCMlXWi`I`I{D3+mB^^FaqNHS1?1$t9ACjA|y;_efs``dH9q!@^Y{nhh9c@cdKX
z1^B5MQRh}`(NLW8fd-(#_i|8hvX?v-h-1kn3xvAu=&{VOaukYAeV3J@-CN)np?%8G
zFLV$f7^E&@09h;rqc$)!f|y>?;$1I!gT?S?`V)E(yU$4ZozNRoqi*Y!x}h2hdT+Nb
z65v=I{em;SS~w`~{ua3=U^I{!ByXiWk-Oh^?zH)*CuZhPHn{X?dLyxEm4@9fE^Abv
zsj{vs-e6HWw6Sg8<e4i~lPWrcR|s<v37&`P9Xw!6-cjb7W*}{&Um1DcmrjsH5IhW~
z1tMyFmTG@Lw^)V5LvX9eTkdcZQ83#H{UZc$ajS@rN&y^)Fl9U1r8L~Js(fm5$_-f5
zAi!?>qTti(BXn`KWPjqmh`ak)oj}8bja{1BrPA{2@&&3`iL=QoGG$99J#mHH@J>QA
z8zP2U1+O5G;ZH0hJ;~Dws}=K|92F(FEDO>2Nz4BFZA^dktrJJWiJkp{NeBHvzZAFf
zBKSPiDf7bn{Or`SalqlzhMpeH=p5xUZNJLll^5df;k{sV;Mj|K$;(Sd6#h;8H98f+
zr2SFR<i39BYQtu^T9r&ov+MA?C<zKdET*8cq@8$-r!nbMe&{T*k>I!C*{^<<gb)4l
z`}7cTUx5aOoMq$%gC}iPv0j%yMBO;D1S*GL2l`u~!#_z7UhhC|oSvQr53^+H0woAe
zphY1wT;QcEg-I8%WbQp|pc8a#8GM>)`0vMOx#QV!Cf2$Ta>5xAmw#d<l9(qXmPA0P
zGFD#ekH-BMK_mPF@`5b1(xl}D9aEE!RM=`Iae|+y%h75v?8ziTV1RLc6ckCMT683c
zR`=N>UV>cRO2k$I0E-2{8Zu=BNyd%Mrl3S}0}>gmg2?oX3((rG?87RYEV8Pql^$TJ
zo#;fEyw8y11`Kb(Lk=3J{c4oyc2rZ*u95yPA2U?i@-HcBfDWlKR@05w^k%^-VGOv2
zHhvqlHKe>CGiytI9K>;kO+lN^St5z-wdBepeeB{HH_Pi~lw~h04)JSBtBTL6R99e@
zYW!p;>ng89RzcND<XQ9PVfiDr9Rn3ql?2k#$8^Hv1F3pv&mtF@`uPYGLfhA%UW0d#
z38}p>=%;K}j-2dv*7A;s(QIYyz%13;ZdQX2w98LGUJC2$pcu$+$s>gP1O#j{LUN(K
zSGa=5N%>`>58l)-+UVoyze?4qycy-kmj2T^0LT*<Nn^9Kvum?vw&UsHt(&SOe2o8D
z_110XY2VRb3u|CS0Z9>h@f#xUe26du$VMi8H9S1rwC@`cFg5gm!$CZ?2Qv>+2#sSk
zG&CfrB0_r8iSbdoMNoMGurPcW#s0s1m=`?ONr`$-He!{&=g646;A~lWAtdXO6I&Qp
zeRE`Ed&fgJ<HtRc3V>UiGtEy5j2{^o7z}<l!ovWgz;#eAwq(7hsbdw=m@sD@1mvtg
z9arXCI^M>@a--)XqMqs!DF(;$?Ch+~WDVxH$G#5qo+3bNJ((}m^Vn4d2@%LJK{N~&
zX;e6K08TRp+V0LBsL}JT;t~>ZPlH5mujBHCu=PF8E=o`IQ^5wQ58Mh2i!t}7^5HM}
zf^F!W$=iqHOem3OXJqOpRrwRvA0{HRo+L1(C3U#h<)}6l$(XF;&`pcA$!}Ca2km|^
z7lm*{Kdi=vUmhdmi$8t*hO^DIC_+)+GrII7RsOc9X!XuJtx$Srt-VbJAbJ|cJgmSR
zg81(z&Q5s|L~ad+qI63d8k)uNN<n0F0>OS@tprQiD~OG(zJ#Iw!RN#~2Bk7lRGQ=Y
z+|Mpynl^<j-YB^i-|?@-VXE>=NbkNb(dNbWc&s$aE=t_a6DP@}Lv4Oaqj?fXD#Z=S
zmz$n&A;}$B>)`D{<_}B%d<_IS#DOw2AOHpZDl|nd0wS9<U;r}L*iE`*_iUuQ$rF#(
zIb~_{d4hTINIQk_1YivV(FC$UY;cn<x%{X5QW2mZnAq(lx`(Q|BV2Cq#D;0OE3V?8
zvKSm|h^NNBi{JM%w|UQZBf06Bnb!7p`w?#OWa8mnatm!-^)Kt<AV=MuyE*F+@kudn
zp+3H_$&&;=*Nn+NtPqOMW0-8LL5rO|2O|=TqkJ;gg9i_e7evn(%71m;xyyHfyDssW
z);bpK(dh35gkS+(nFCUX67c>|P}>mx4m!p)Gvu&ySdH>sx)ij91s+b(e<xcf;rq0}
zF9FB<Dl?}n0C*Z+kB8JMd_KI%Md{D{R7j@^ACCTWPy5h)Oo7p665Oiv(2&BOaya~X
z48sk`u8ri-pd$QwFHTS>E)`W%-E8_q9_q_5a4RjT&%JJZSCZiNoy$vwTmNhYQ3V(Y
zflfRpAvgveV-jNGhv2CYCXZ8Lc{_JQbEjWgApHZ}62y2UOkqZXDvvmW#$6UD$HB6V
zjC84y|D7ds5)SqEzs<oBh71xTD-B)?#M1@FG*Mm;0qPVewG|n14mzm%hj&+tE3o~*
zjSF8|zr%LkcUy<%s);sAyACp=JShA~WFT`?K!w{O8};E}fqd|kWY8zksbc7tb1WEr
z$jUuAiyjGxX%xJb$nd@=wV?oE%c@hB2aOjDu{Hx{K!!g8$yf}Af0&kK(yc{g9hNH>
zmS%RZ!e}7UcNBDj+<^?r>ddNFd!B#^Vl4Wd9zUhvbBn^7&wgtE7<v)QxB8bB1M|1S
zqYbqjgIVKeA&Z^sm53c488AdRD41YE7%zPo0YT&wt8%bLqX-1ukMs<ymn|Qc6SHc&
zeg0@pSUed0+ybcV<Q1X{#LEm5a==f5si+GaBP|^rv#_Nm&>oKmC>`Gi1qtTAXptE(
zb6IC1DNcpHF4JElnfd<cOV)d(e_f6D);Cr0Z`LwRzXTV2w7|U(zr3$T*Fkpu1b+um
zlG-rWDS}374z%{kfJlK=?a4QkLVOukrE`(Wh?D$f_AOmBAxP+JG}06L@|@K#VquZr
z+Tq&?$ad1NFv@Q@-lJzHqE=qh+lQ;xKipXy_mepOR6jLYJQ1`=0tOMbas6i5nv1*L
zZ+}XT<0~}ye(T=-9)3-OgeVvE;QvLF)dq!?({DeeMuu|{0S0jc_LSP20MM@YBcl1H
zTPc9Jxcre!%5K_oC$yOkr5hQJJt1tmq9;U^;Z`jktg%z6#+cFIUN_}oV?(YF?$kGS
zk{jnNV3=>vuN}C^Fp%W)zcic;TJV6(ECE$LZ^)vN+~RC#kf;1-`73|)_1#k4REd6Z
zdBXzDKjYu2#D;lO-P&DK>78obJT{=gC|S-aFPZ%y_f0jE44TL3s-rc;1jP7%^W5BY
zL>#T5Y;|D{-$sbbMs+>iM5ot82CRx_90b~}??x#2AH30NUhe08VHxCkEqxlYIE3j(
z=$tDDxyDoEqR2AG?kbZ(#RfVX#y_1}FOE}IUmr=8-H=uKU&fHBYQ5pUh0nLv{qpg1
zOsyzT1_IW5wPb-Q>K-zgzXJ+`unnNV0tw5qMHnW#cK;p%Q}hm*v9Q2qat%UX*mdVF
z@XTpI17TrdZ2-F&Ik*NI;5?`#w6|_~gSDIk$c};ly^2g&+_gIFelAxvtTvgAs=Z-z
z1Z_8lgu)CmHXfM`%o2xU)2T|YJ(E5e8?^2p8&z*egW-raiEfseGAA{=iP4wZ`)jhS
zu2wgITKP=*S?XIi_k(9Xyx-m-ICimO2<Hks3IXhS-(e6GikFP6>@-Zc&HzAh_?w^A
z52^DO9<<9%j>>h0G-X02@`U;GV<@)>u7C43LMUJx@+&gJxRB^eH4Qy7GDe3?NdqIa
zmL>A~Q<e>WLK9_YuvmlC$LPsNii+U#Odhm_)#@~+VPWYseC6*2&wYRr|2{nYiTiLz
zj_TxH@4UO;#`<D08W)iJfkqmdcm_M6f2==fqYr8xITxY%c?F!Z%FFSU*Lu!V*peX>
ztw{7LC0&5Iejps|@qzUTfMH&g2<w5{;4!Vbg^lSs^$xZtGmqAm?2+=Ds+gbv9s-b5
zjM}b3m^OmqDw(ZWcpV(5h%F49GA0z&&T9GYKD{oe+L60mqtFYBM_(L=(RhF(V+PIF
znAMR?f!Lzd^{4e<6<uK1+!<ai-Q8@5Uo(YXORfjhuGt{e7A>jaT_GL@+oXvI^T!G1
z=LT@cmk?qzb=#O)8#o|cxDUsb!3Vm!R4}~t96r7;E8UI*rRr0O5!pa~6xV#;ncRMc
z;}!4GC@Is^pP6d8sAgp3<!P4sRFS!#328e$#i(&s$}>Eljq3&VB%2p8jVjR_^;dc!
z>6AIHBF|A!fWbYDaz_@Bk7hvLV0+3YH5+6hdK#!YcVo5$8CRikQvA=j$`?41<9HqQ
zNYtNw%UAWk5XPtoBTvqqOE>b?*}gr+aaCj&y!phrpd0o!4b1%#A#Z(s?KdVLz*1ZU
zr?(?`;;q6O$6;Jxlq|I5!PFlRB2_aS3#B#MR=Blb3>}(G#6%8$iaUTp`m?pRk%?<4
zsn(5mJ#WDP6_js^FMrOVuwCe%R$KJfziDZ>mh!K_7V-8_n*_QhKZ-OFb@7Vn7EkrP
zAHW`SeA68mb(S}9Ojx$z!w6t;*0}qxS5$}~{gc7L8)5UHYd2wcq!aC|1LFXoAbQW=
zyxM5MpFcd;syy)E0f9;`&TjXvtM+tEojGZ4^GlvpMv{N?*6f|7*Afjgk*6>;@SuJ>
zappxq5=TVw;V2kek$H3q@bOT?S~I;&wgDH<OjEGEa5u~jg2OY(R(3I9j$|*3vnW{L
z6LF{?8)LzR+#7JX=&4+0iGaju1`5`~k#E`*a8<Se9-j1KnfL(hJ?1w(rxNPk&>578
zt924xnU>5U#?ScM6;TzfnxIc;484A~>pxk%W(`~s4mb=v!VH1+qXb~{EPVkB2Z!Y5
zP0iWcB6U-bebkb%jOd&>IXMx-H9k2{gF*gPw=BKt@$6b;I8+-3Ng+wubn!22G>|Cf
zf%mgmSFGb%@d%8nGC(bZv4aHXXC@!=gO>~x^ryrBnxDuUw-g`!<9@mj^Sj#Wl5=@?
z_7+75$wl3W7li3seK>)BQlzuBQO#0s94S=i&Q5Jj60{yaMh=nwM_sUO+C=DMXuodO
z^*VWp00|gpJb>XH7ucP#)6xdcfLW3rUNK!Btf#(**_{NB8caUAv~K_ep9DQ@4&l8)
zoM_Yy1>2(5Eu?@W&!iAN*$H1(7^ZX@t%V7NF-LqD{jyq!(hK@k2BbJoDohj6FRk;W
z59`?^mmp*U<&9>bHzmll19AT!w|SfiggL3g9j6DhH^s$~Q_$bW`)q;^{ECMPpa=Fs
zhSvdE5s~7Asd=p^V9#g39BvGbU<|NR1|v@`I0N8{(DL9Hg@~tstTPvXg@#J)PQv4L
z3?6SkJg!T+vv-V6xsD;LRzOn+^VtJ4uaKutXv18bM5+w1>CgoWd0K7?LDd2mfX)%d
z^EID(x7wQ+VXN{%3huxE-+s?i9=qcAyhqwo$ymAA=Q5*j2KasYdyqT8I#Ivfbn2a!
zu~KJ7Z?Bdx0=Z(~&)9cu7*9f8CNEhy6vNyf4m&kHV9BFXM)F^a^m_1nEyyCsX!ozL
z63iOk-vIRl!-SJ4rx}Ham>4-Ls~kibisXG5!3pFb;K}>HUcXIpB5%q#P5IbnT==un
zNVuSW0|8U#W+{1!Uj0@*uPu%Q$O5xxLrMN<7S50Iu9X+=K1drQv30Rmm#YW}FSr11
zXY0z{r}3VLreNe2qLOxW<e~OD=0M0$4h{|qc%Tl9&2H|ESi?xpf4S|#tmoQuHw0pU
zR|(xPPGGvPJ9W;SvHKPZImcqBhv#^_j5;(cj*Z4tL=>Qf1uZT4^w|niE%J#p&d87!
zazf!D3jWYA-TLxk8x$Zo9Tde_FbRC22I6nEOaaOVGAxdaaLNkNi-^QS16ztn1Praq
zNzmwK{;d3M)@YCq^N4oHt=IVz<9~i|3-#>^8POF!N5U$hs`mqL_a7pOe)IN{{hUL|
z)<skDd@_o9Vv0*IwLViu`YAGS*1C<o93V}2>h-~Oca+!i5U^bHl+D$}L6}2suYHZf
z4?}G53=HaLq!bkLHa2-9287)(3kQ>hnC!&c51!H|pITL1PXSKaxfo*l;|tsj1kg;A
zAnQNkal%KT@Y;3Zey%SR{<Lr15xq)x67FL>2wSbuoPXmG#+WvWm0LUtk1IeP1B8sB
zz&!N1CY|z6pFR<Xmn*lXoLt-!`lfqj>4mBp(S}t5Q|#Segx`oNP+uWUk6N@l<yn#G
zzq)Px;=#ZZc>98{TbEIdfDm?5LAEV`xsoo?Gj{|z!u0au;-dTBkO{&v6JEdGLK7-9
zJ#1p~Ve?^lSB4lZYQ9C#z(V!w12F!MGVqcUuXFB|j}ts92h1f+yU<48K_RZMiOChh
z34^f-Ppo{L;K+1*<_`S$@hb2o0Tp}j6qldCd=+DXtnCa|a^#T{(2)pFO-&)EcIo%7
zZ?PycuKVUqEZ<~C^+;>1x48Sq;zMoBEj*G6#$UHk7j|n=(_|A0J{@t7%Wv#>05jws
zd8nO!7X*4~Z|t6d8#rJ`Nd5#=2*TV0%uIs}3KGdXFuLcv?OP#(b1*j#P{``aT^^pu
z!9n%n>Rj}vvCLpgLl}y?1i7UXL;_G8nbq^W`}$Ngi!55f_#*?Y6Qjc^JVT})hFumR
z>w_A9I^=$5A{BTWJ7UhN{Qy7e;XTeQ*MY(6X>MHQ(%KKmh88)U=i5ru)R~$0_q(cO
zth$8g`yh%YfJFEU#>#cwVX7J)&yp%1jXaeAEKvT%(;sEKtcnZZK>%R+3K3?aB(UVw
z{Tatgbj@kE6p_`Cu>er!D$Y)i5G@KG(Xs^(H2LUe<?1Q`#9AT%+?DalBMZ?F)=&H|
z5b^uN)wO7}3!V#H<Jp@z#b<Gixf#NSkT3NSieCixTKFPdt};A;j01xs_Kl=!7?Pdg
zK$2kE00V=_7Ck*ZLez<9Y~bs!h#uF#%iE({b4GtG5@CF7{n5m5WhC?pvw(b3VPXKe
zDuj`MbQ%hSj8QJ}@O%pxGnyTDTSx3+Ci~AA&qJjO3;*Mm>la&?5-Qg86xZ_d`S^Q+
zHxp7t)D(Q7a`>_g9qTWK2!D-(5q5Mk>?-h3_P=W^fl5PWEcB0U_~W_^UfPMv8Ww6U
zf+!32{<_bH<+hKLoFoEl)Dp-EAOAS)0*)Q{^(lB{0poTBP?A(|8}Qtdu3Lw)(okb&
z?;ymK$=rJGRPkcdX@vCy=})DbXujQ2JAtgHOmP{jU&8Gv`<vQm{5~d+O7vgSv27#;
zTU`MJ)YBYo89=3`!#~pL4M{|a!MqJiU?mjI**<7CCl{t1XPFQK^RN=vbQ)lLkw<MH
zDyf>zOV#vMaOA5Y`wq{O$yc#yd-9*LMxX>qL%Qh^y!G+ngqy_d4~n5Hc%f0oe;m_f
zIb-Q|OO@j>q^P2(k1pql^~D(b(oBaVv>B{6m-RRWG6)p`z&anW#7(XjH89ear`b+{
z!L0nfju<}}VW^`8k4C|nkoUZm7!aehy!>n}4s|k!{yXav)#@YQ96)qtc+x~?uVyyi
zKafuK{}lG#@m#;{`#3^Flw@QyUPw{ON<xdu$f&Fe8HF-3B0G^)MwGWmTC!(SDI*yv
zJ2N9ABg#tl`W=_L@6Y|cAD`dj(;x4LcjWaP*L9ued7Q^_oUe1MPW*c{ZLryASd(kb
z$HN)5+)W^o=ftx|cx%ed4?AASRrsh2eI}^;$%x!uutIcC@!~lWvJBkB_yC{8#PCfF
zw-1MHE)G8b?68BNPV!o(>wwm$)?b~--g7??>z5f0UPb;a#Mw(xL1A0QAmEm!JU40L
zuYtq=39faA?(5qmKL-N{0X%J8)ve%Q0Wq=K@g_(R`mdRZu?BS@m1Zu;8P3JwhThHx
z-r=1T{d5O5i-mQ`*R99juH`UtkPq_&2yjMqBVR1(lxUcEUPVdpj7{<(P)0})#}(AO
zuKoJ;%Mmshjssy=Kk`y>zPSE~)5Mp*>y*rwJ|A!IkOvR+y3CP~TtK>yAYC;6z2(&|
zsNM`dhVs9FPg>E~zbe*vjbDMqG8#BX-_5qWG~ztMIVlCFweMdt1yHzd=I1X?X<W;^
zrcB{FKTE*8@<w03v5LUBbth5}sUG{9nk7^0CX>}8g9|EWI#wj`peDx`*}?6|NF~7W
zZ9KtKB7mw2`^BT<k#&1NuzEO)h9c_*mA5$UiE5|%v7Fuq$So8BJ5j}Jr9_%{TXjAq
zRvzDg_YKV-%zM%?)x|^2d6Ad$Y^-SVD-08Wg%r}U>mJ1Mf^I<$DEn~8`e;cmQ>fHS
z?eg`ULU4j#k8I53{l)M}Lc#y!wGW0Foz{^&etyG(8T7*d6h7g@;i0uAT;F*+u`Fl=
z$7XhNI1W!u+@h)4{JK&P8^R8VfK(5IwphDRoAoH;z?;k_j}iWWM6CnKJD#Hx0>2YH
z^K>o2w3fHEecp0F`PoTp*Xyr63_w|8J{2*VI|+=9Z`(FMWD9NfEGa+&tX!@!k9vY7
zJXF7S9ldt$GOhbxm^D{@k|^y@myY7I{Ho~e%p+_hWpnko(67r!<yE(d6u6*`awjCD
zfGk8{(tz+oen=o2l$WQ){u)1X`{Xnct(n(*{^t93kasD7ysd&ph=bVd@bWH0T*8rb
zGmGT@kW&ke%Mz{OPlwrlx7(j$cNdt$H8Dm|?3>Fl<SEl=>>dF{84b#3h*x!~dM!&U
zp*e@c-Knr?(=WC<t4mfgVScuZCr|vF#v%2tjY{7AHW+gC92`mB4f3j*MxK+fxuPM4
zSTI<MfZ>Xh#J15>zAp<iyh><j8I3D~jux3ZCRV|s!73rm+uNI&wA+7Z8r9t}r%i(^
ziN*yJ*WtD-b9OxZM>wycO1g_Xu0hU9AkimJpDN4zC!3UzulP$_=^h44OuP;99#hIU
z6fTYkGl!)gGclThzMQz%fRL7&H*+aX_pECBAKtFdDIrm)qINqpaK;7&oC?+MUVovX
zV-R}*qWrCBZBpS3k}N7Jss=|GMQdN+RFN7u_U`%<5(jEmx)hE{@SlmH9k9zS*T2^>
zUdvfXu{{IGrNv?oS?mNK1^)tOc~q@eYE_;^-?@Fe<Bb$M8{1JRP0kvx3t0sqeC~eW
zA)&y12Nwk@g5p2_GI=Nsj|*JyQ3z<Khp#OvV~Y|mLc$7Yb>^{I^I!aoCe7R0F@6Y7
zq$dIA21vrdaD(!Qi02Rr^U8NUqQpIi>&DytftPo(8waa9(T5;C7#H37d@Ds@)MQgD
z0B3lc86mr;KM2Owt9Ey3A+6*NBpdwy=4aUaGK!>h@wXpXCVkDjL!D!foTPu1Cr}EM
zu~Ebw4@qzhbjCyZi|#O|`h-w7)bPX41x#M<U1jS3mN89wu1la*+EJr!hT}<QI>tuZ
z78NwHx&lBCQvJoI^es3r;hW;Zub$9~vwQH;Pp`s$9j=R<cN7%s`kYT6iq6RD4~B(-
z)uvW};Qt8go#}2OoHeZ!&VmB|D%({|Ov2+5*l`ug9dX-B7Fg5d+Cs|`6pB&Y4I)xQ
zI<^II4y>}!2M;PxR&)xoZ^`4Aefzt6DB=Z6f*=A~Z%<G5B{w&0LR4q$w~9TBjkZ9)
zj2>Sa?e?&99`&EXoE}p9KJJ&#EMy)9F+aDozix#-^Y-`m-9`wb80L(mq9wgJoLdC6
z1_CInlX5{dorPKKXFCmp71nDjy=9uX1oyR$iv5$|O$rDI_#$Yn<Ka<JQd&h?M({Ec
z9Otf=l~i*~x#=dEW+oky_@tcnGBQc0pBW$EBGpR(Jz(r~23RVn93#aqb|Jd0;htLl
zi%HJFKE8da3mrP6yUTWcM`G8r{xFK?)!7b!7yy+Dpi(RbIqHa`<6cmnxv_FcLtySf
zTPNSYtRcLs`QHfdNh+Ra|5V}QBze$^P2S#%y6W{JrzQ%{R&4cA;GDk*dE_t1j9+>r
zTv(3-ip3?6`ECRiTUGx$y>Q?fpGEDIk_*kLDenUwIx$N<#?1gc!8XYu04R3KYR>b=
zQ;?rDM(jQDba#*VPkmzx-j?er<<}f2nn<~K3GTh96b_!;im$_7XOSc0QMWDRQ8Wv%
zuL&O`M%DEEzxfz$?6&(l$D~OSzk#~&Ky{#@<S2F-R8hpb4ifj(2O$L`9;m40g(M2N
zHZ2*`=@Lowi`UP>|D#7?-PK(bc>jH_zn8I$o!qgQb(ydz$v^G98$@{Y+Cj_)Ff1X>
zEv2a6RhX00Oix`9;4*c>)AM(+!rVE{A9n=n`($mY^vw)^9<W&}&VQz8SR-+N|LUpP
zK9=njf5ik%rFri)s=zE6QvS6tCF&}<GWj}4@_qwd(UaqHd94C)CkrMrc-%o_2HK2)
z$2Pg$MR!;GyR7T<*mM&AwOG-Ml~-f-4ZIGjPwy6)ZdP4wm3IEEh9&>F;><0r9A0~#
z?0?qU4!$YdO4<L&nl?>B{i7Z;bBNOWB@ENpzxE|Vu!>TvUt6b-$nmpzUKi9h%!bRR
zbzAW}z_y3DqJ)NoP!omG6lm#49^yWRD4wR@V?T1Lqpf_Uvrl-MO!)Of_cyT3Ee|gd
z7YA|cp@v1OE9w~NE?5X!b+yf&9hHROcsJqkHl%2(X?`-S60jZ~ZXIl0E$OamVDRV*
zv;WlU@$9V`o6d>z9u)d~OXC9Sc)@T0n3M|}RD1sV?GOOd)cvvBHL&Un`aPH=4}YjJ
zUC{3`(t0?HWoZj9%4!Hzh#e;O6lo-(rm`uUz2*D!v{yw}jLQvH(?+e|`}*3cJ_*1t
zRmXD@Rxq7Fx=C)3buQ9}(Y{Y!aP$1G9TVIGbj#(gNT-{vPgB@}OWO{Njo=v?dtlQe
zaJ4E%&p(LUBiKjIYOB1xcbR@&`z|Y7M+-Prd5a4!01XvYRQM^B!*0$O(p|DLEL-8(
zb2!)8)TZ`m%&KkB_==ua@SmlaT2zhFHE?D%YM6zfc$QE!^lXfU0g7hqew`TGm1&w8
z-Y0~bQ}21|-jiI%6nk(m)_Txl#qY01$>?Xq++F=S({4o?>tSdj-Lr1r{o-bY-HlpT
zKR-3<@q6!PuB+luzANEYwTlKF(n^BD>KlS?-<3y+lkvp@%5aKi1ojs^754o%i~n8-
zhpeTl0N>Z;(pm?tx9Kjfth2CwxkbC_Uwe=C$e+FU`f|<QkWcSBJ!H)T=Sz4pS_YH#
zEd#!5?mrZBzIa^@EW*0M|4OX>2;`2%viXILi+Y!-Lt1fASn=8GDYuAe<i7JE-vCR?
zOE?~(MfDvPSaFH?LclD_i|!_$6T^pWRc>xOlc)OEoHw{XT*k(wE5`CZcXBN6FO%V)
z7pi{@r2zKzu`BAMb}`)W@%DW>>iS0oPKf02WLizkytd!ucYNsBwGum+8V;l#c>jhF
zX%qWrjWqoSKu$+cJdz$T4B5KJJ)VVrlwQ@(l|duYCog9pkt|rdx$DXPMn3cj3CfcY
z0P6}-`}4XPHL1*|V$I`K(<5znN)|3uLBM8p^QG$_ih?qC$j>4LT{Wdvg~fY*O6O7!
zxut687|yqqIv6m#8`!aV>sI=M)3hNcb)r;F74wy;p0034o*l8*#wg?)>;+;F-D6^7
zon_wDhJ=M3(b3_6R{?@rBu@|(%vEZ4c}FsjbZp(JzFG5A|Ab;#N|mhUvjMNg8JgZU
zDwQ|y?p-bqIw~*ou6?#c307>59xVg)Td%FjOrt)w{|XPoV&u;SgR&9jU(~5CS1)7_
z_O><$4($ZPcKdjCm)1Gi)r!LhUr8Gpc#`5NDs(1B^s-yBer_ZpF}w6e78W_Aljw2$
zJ|5F`K-atG+uc^y?Y>K=_x`S9ypBh@_mX#dtK!mQY+#rb-2>bH!D-v8Rr`!jNF3Y$
zf@f9u{vQ*2%kD%!IOYDd)=c{}Fbz@ERrNIH6jm)0PohOZ>|vG8GpXP_dFm86>pmZ7
zlwr?_y}Qxn+GyzepzIS-=!9v|m`r4A2FX?tq9o-*0$!1k=*#y%4Iv7Cvdic3u*b?%
z`Gw}*ue0<nJl?OZq<MS!#cQ0iK<>G2De#|6rXT*e_4vJaq>=#iqKbHd;)N5h;5<62
z#A2~7#f3iUj@CQye%>SaAA~RTG_XloTgK5)*^(7s=X?eo=X=7bjrWm*E#id#6kp~{
ztKdfoPvqS?>q5lP+K5>t6~zRp74PikF4xx%FE2{Ez5UY>cu^j;?#IZP-2)6){U_4)
zE4ZcJ6?F5;*lep{8Q(%fDiZ>wlcn+*_}U;G570I;fxwu6G$eh7vgsFDz{r7vX<xu9
zmQ?<5QiU1Y#%v$Qk2R~mRqpKwFgxZjv8%+r|H8u~iQ;0ozNC(+75ZmylU`z#`lcNx
z2|Fz_)yJ#lk%QRl&QmKy(qa?;Yfs)%g}ema2_5OM8QtxgD=%z>5dO=TFR5<K+f&PF
zcM(%_LJB&%gZLQ3L~R7mNz7%!0x(aGt={Nu92b@k%xDgCStIvawsx)h&;jWX!Ls{q
zR$;LK7f;e))wl~wwrnfAJ+WLPZH=TmeHLdV59{CdJQzQU6*4ExB|#EIoCZttqS2M;
zr+BH|N+Z$n`CKq<!;Y779T512j{SL5cXL0^EgFH03@pR+M@-O{OIhofnGLOED4M0H
zjEDy;k`q+Yp#<9i&XHC(f#QIOhKPA4=rYhLd>}ZTj?O*yrxO%yMtvBqL1?$`SJ3-o
zjf7xSH+-M%aPB`Q$9hDe##s<yMp}G3`hlF|&Zj%P@&mCmv{(L^+6F(u{jlUfsF)=?
z33znM`@a8Fzx!8rA0L#uc*z^>Ig+>mT_Pw2F~$ZDg95G!u|*)=Z^uf!Ut;2k(up65
ztG2!qF}cem76px|gKcX%2&;iYIceYB_pbGa;cV4f?n?o!fqoLwiAI0y9mJ$pGmqho
z&n}2g+`qC%VtvS`68N-WPZKDOSm40@EHoklK0i!wO@zJz{^axPd^^zqCaeM$%0ISJ
z4+MsCB!>Kes8t;Z4BL%r!PR2VeSdHXU!m9is<gBPdnOcJ9B*JWuj^=s(QQ6+UgQ>V
z#dRAuR^tCY?R?eRHF_^MwPUK{P(69DqnGKNMm}z9nvtb#@Zjn-r@HBv*6il`?pzz2
zdgk>X>AzM|dd4DXFT<J!F41>?=$M#n1w9SDrZMEnHN-WoulU0UG0=pT^oYcMkNNps
z0XTSlWnkly*Y$ucB5b3nQ2WKpv5h`DRVGN-xi4P45LFhBdUMhS)}bp_uIvPVc%y&K
z74MaQUij!ay6oMHCQ~jMqtc6Sq^e&eeL3N;aM5AwhL){ImOuORdKE^{1NJUv)8K;B
z*u{$%U9rrcMYI3G@Z$>ah@F6Yz?{sA#GtH49%x!F-#yUh=Vv<Dp8pWL#rkKH$R4ZK
zQsZFR62Pvyonwoj>e_%QAe{UA^T`&C0a2-m_q1tC(cAGRydu+LE!>%0U3up+*IN`I
zOE9Y`lWt@S5z8FYHs;b4FEGzK69%tIUo}3*;x*Zy#5ign!oc)3N=}v43Ne=wK+S-`
z^=72F=!L+rkhs}_Q?+BVb#XR-@f(bp3~<tbD)*+&A?(fZpg*5n_a_}txMULaL+5{z
z@l5tA=}UF3ZFqAUUO)DG{uAyv5lOLi;C5)TmkG93^3=v5_G|iQN-$P7`+d<JVWmUg
zcpbSC^1}Q!0DkFjOxpzfk`=Zh#ISvdAvT!<AZ|A_V8k^VDoGua7O?8H+kTgqA9omP
z?N1$%t<v<;&O9hzq>F=y2<E^EiR8KI*YcNpaFg#UHL4O0X>&$E8_q*r*}k1mI|b@n
z=jA`@g?vKhS6b%tD01yaV+V0cSWs0^hZ{_s%b5;U-Ys-DqP=`a6*2+}mAH9B0K7Lv
zTd*T{!E5Zc%TP99RVTJiCf8lC)#Y^*x;O`~3%E~|VO=Z1LlZ4C58%sD5`U`ty~>Ui
zD|I|@hQVy9bHSGJ15*JGkjk*N8VpA(ePi2NT^N~11N1&DYF%t)PEO1avVG&~0HwZB
zKi8Jw^wDPnQx9JK4F@D7ph2&~q+XP-@?c*WAoc8MK<e$J(E{)PDtim~@gIp6;O7^y
z?l?#^FVRlD8`4V5SG|Q(tWjQsm?=Ec%3&E_&UB>6lE^IrTV*5uCAY|BuT=Gl2<Q6P
z*0K9_za9V8XMPZW6D=a~r(T%)wg%>&r&xv4*)(0Z$oH_we~P0`dwtIEL+jp3`@V|)
zE2tq39X(2q)Mj$*NRb|FsFt9})@(5X{NUIkoyj8uJwwm{p-8A4pCc4L@CU+#v<dmG
zk$AnDxYAbDqU;v=0?5TbIG9kJXlbW@#_S8J{U=(_*>>`Hvu5e_^Y`{pO{DrpG8+C4
zk5t@tgp}i?@dE&m8NeTr;lmrv&_T!cx#+qaQU<8CZ!Mm8ND+N0U|&ugKR3!YXYHay
z+u;d{IsDXt$G=Ch(PN|^H@v-BFWNyh#=^z*!+8N%_U-%k&;HDf++#T~>ZtOV6VLJq
z$h#2t67Mw_d?k-a2~5@kfQVkD)}4H#F|mmzmeo=_wXWT;q2BuZraFHgt^;0gYKELH
zk1~*ANlj*y7V0kv(7={=p!zZaqh;%{t-`*OhmJnaXpfnCQoVG^a3^(_z9#WfAjS1!
z)|q+96MS7>*fjpPgxo?3A)~@Mr&|Sknf}0ZS6{ME{Wm&)-vWrP7@QhGVuW%yB!YvI
z_m@Qo7V!2ArZhUg6(AD6;9{%$I<u@iw?OCK_AVNZ>Fsr^56n)o=Klh}tR~k5Y$Wz~
zCYi(%faPf@#UXk3%TCM3XC4^s5+YFz{g%mn&o_y_<eZIU<EP{)A*V0mBNTEdo}K~e
zCfgD#?m4`ck;+A2KFfv;{Ehs9{{ySXh4Zo}aQxBZ?&FRVeym08p6+FyGq4=AxMv?0
z(rWpY<)~Xick12mfUt;jb(ilZ_YSbadw|HEM|poC`VQLGT3T8mzK+!|9|l@xGJv-f
zvkpjTYf`$E>U_hIp&Ti0CcNEuz~CSn+)5`vvDGaYeRGZ4DsrPB%y`LQ09a-nSeWwa
zDGMi!GI;Eee(dQyS1AbdgvWELtE(w%n9uzc-ukT&nsM3lN3(!*#Hsi4koCK(YVY4m
z7o4EV&@N;g+aiB}<>Yf5S9x>4lsCFh84|uCKA(7OVm&V*hX;{d{IK~|E;;Hk_<<#-
z!sYJNa}wv!W3x%q^Y*TJyIx8cWLQ{oQybI5dHeP&UKXG$o2#0d>)8GMu4X|HfCtPb
zwFj;sNVabBuqTwGPJz3T=f5h(3(!+RF*5z31As{(&9~{|Z;vOR8q~k$->!$oM`>lH
z2?6Vtb7j7Qhycn^ZovOLM*qulHV$?Kxm*GBCj0`9B(%VjzDr-4fn7%^vLj1pM=Iwr
zwzAWzAJ#ETopo3m@gr05hbl0cb9zGmp5$i^btA?b1!;$Uxmw4*b~VEC>Q?<Fd;Y7`
z`y76jg1}h?g8{7AqZxKs-M3z_uW8O$XMH%=PHAW6lfT}}lNo)5L-i%@8jX})boBJm
z4TnLf*x0je_sf?D2~JN<08f;A6-;tPR|n5n@?FiGJgC5-E6w#Jsfi}TwuJPcNI(I|
zLEJg0h*|-AcZW!}kiM3x6TLd>SPI%%lZx{^tAY+GSbiJ24?E`PeOzag+!=5BO%rRc
zOu4nAdgi7D?B4?DPX2hvz`NgbX6CG%;1Jr8=V4R-!n`qvdKvv~9J8n^h?`+KlS+rY
z=c_MAKpH#wv47WDn}?PMHnf<?YAjclcfJ43(Ni0x6I=`L$ZpDWUmf~YB>Ghzw@+S5
zab5R|(<xd+xutesz6iYeB}ku?D%fGu0=G<i&Tm*&OckkT21Q<26yC}H3&&GHVKW7L
zmx|b9($WhiKrCdlPwlQ*dg3l5!J|^AT?`g>rW<>Y=3JH0olB<=DLqZ7!>_+n0@sai
zm%c_9w`@S%H@)Q9o>O=ru1X`HPG_Q{!9UFFrcB*AG4_3({g#)fnt_uMCn`9Kil`Xk
zpol%^MPwd54=;CSwsd@_UXe(rJ-?<|*Tg}C(zWtzI9E&Wg~mE9@+)V`v?@RSJDXtI
zP47{<hH3ZcO%1HwH_ddkw7apM8)3|K@u!o(1kv%z7fCja8Gfl9851fQ^JD0Uj8IQ|
zhMLPyfdATMi?b8kGb6~(sa5ld%77YU1HyjX*8|FBpDuuiOrQd#UuROZ&b+J(2gXe8
z&5%sp<XL`84&;N?5pN9^DoU<cdQ`r8?-93UcTwwW`69dLwVU0R8C?JV^}-$kZ@3B?
z@tGF=!O}9ig**0v56J7_8>{B$7rg!ClS4}{@OW0gvg0g;H=B!lRb$SY@cj7hRKI+o
z9aYuwvUTyMCVUik{9sSOf8po!uw5(cv`&Xx9~BO6Z%qb1Kflwfj*XnLFvuXfG(QE+
zusd|hTi?O8W~_9=x?**F%lo=?`66XVgj(1xgqu>qT0DVjUZZrwlE>HClW&QqUU@`e
zl$0IFS=;({Ff*^Q<yITH^g9(&D5;a~Hxgun1g`F^{<i=1g@n+z6ix41r<3z-%AYJ#
z$#w(np;7_pv`U1ag=vNkas4$vCtoBET{U5O`}tk%>3jD)t|MHMTA<E~UasAT*cC<i
zIo5t<L=%_B{Gdj-?fs=N!H+w>;-9LFjPXCMef?7z4gyVp`Xfo(cs=v>6NZfX_uc6L
zP&rnlwi|TnrC+`?Dl^n>-7mIQw~V}#*;Qp?;=|O`1wDJTmtINDX44{*EaGNP#m0R^
zUnu9{9(=tgX(-j-FF2!b!`Y>hZ8WFVVrsz9G3{YUYYH<{V}+NH{wZ7Sn8X>=gp8Df
zPM8+JL9USzrhk{(sVIAOx*g?pR-8IkA_D1OOsTy5Cv}}C1#IU*w4qQItj3HSGllqa
zm+gJrW^di6WF}kvvsAul2fRsUR~awrj?+<wo=hj$YSFW~sx5ZBQoLlpslTuNC@BjI
z6L_zVjg3u%+DZKyT-=qb>q@1VMf{h1`yI6=^mXQTTg<gO>c8Ev+WpAWae3N>Y5FLK
z6_*+8p9)>MxwOL1y|Q%(%`Y^u`J8LaT_i4?6;a!N9*rdw<IJ^YArkP9s6r1aak68U
zwzPuw7UP?Cy4-JD;@^C!@%P)+^xHeb((bgk+KIngJ3{)!vWZ1YSDRZynkKc-8l~n=
zoWa&}Y(REW49y9^&1c}2v|1N3)&YVjx5Eib%pK7HA?Atbrk#3arMK(I&m%!#1AyBL
zW1otKo^(f=pR;9-8^3&w&*b`2SspK66ppI&pW%S{23EBZL6Sh!cZ7*)H;G90udcro
z{-y3x*5_CiE}UFm^@(T;ME}N#YF(=$Nh%?92o=F9K)dAvv=)e-NHl?Ll0qI769I4l
z!0iPTo6C_M-V_{emOfDS*w3%@H14jqZN;Y5BB5{Ve-w0RcQ#y58oH@)V*BFq>0^(y
zc+Kh_zmFp`1NP6KFO3k=oqn!KY$`%1#8n&$W74F5jnaJ@gn*UM6C>B&H6BF-C)i-b
zM2y4}G=&2zOpH-7=?AnbW9!JvcmITG#q+q%z{-(Eim}$Y@RrW#`f8i?I+iBKwo&qW
zyML`_Q~qKJDh3#3rduaRE~=y;HGzq1Aljpy)_Y~RqZH~JAaI~m1#?TiM|6bKEXbhH
zoe+Z$TpImB1%(x)aDo>~^y+%<2^nK^22n{Sp#GoYH?MJbTdx0uYH^ZgefPn!XzLHb
z)?dcdQ-Z!QV?KjzlOc>27R5h_d0<VkcAgqL1&h{-=NHQMg6a-+3udjXxHU8Y6VQJ=
zW!VQ_diYqDS?a5+cW}ybK-{f<dz}y^`U>d6<o3ZsAtnvP84(dyFl!#r)N9P%xygMi
zD98kA?Y$2w*UfhQhb?XEB(==PtkOue`Bm!?kt`uXom-vfvWh8|;{y!+@>n`^3OEA=
zC>NBcse@v;B(FzK%A}lH$~&q8I&wRGxc&};B1-5+y{k1S9Z(m6SE}?2kNa%FE;|T4
z@CvIR9a_decq9kjulG>Rl72IIA0=4UF#2*?GIWF;KpmcD=_rm~JolS?HfsOezP*&1
zF|)GvOjLC;2Yr)pTw<!lxT*@wjqS|2y0Sd~i^5QY#R{+3cMpWN^sYjw^qD+6;%f}z
zv$*u##J*s)0Lzx7%^knjy&$Q4m8s7zSu|ImaLSGjhPxytMC^vHpJXL0ymyR#1sAk`
z_%R$U8FM%ef)iWfj=Fmrv(#y%XT(62YK_bY6}#RQW|Zy*3v<)O&1YcH(2)EQdFn0h
zs?yTU;3wR^eS2LR^MACv1doL_B7m+$ITdV@>R)BIFEG@ijKyH`a!JOcNF4v)-PvjD
zcT4UM5hFvwexcT}CJiB|OCotx>qD;oF4~;I;PY^z|Ia<F%MlQ`N{vED30BNsa$R02
zGAky8=1Xejt+rep`0dl>U6k18Iy(qq3j>w1>-NAC+M*RUH5owMt6y67^XnkNgDTwN
zY?S^}2uX_06K^=!x2z%dQ$V$d$pt~SFfwN?M|_^*ajla^l8-pA<k2~|u94un#~EE3
zF*Efon6TktThs$9DihyRS#!A@XpIiWuC1>JRcJ1>>aN*dO(P=;oPRqo5KuTv1#~DB
zXh$}7(Z!vPJ|v$bBqG8bSKIH5K}`h^9~^;_N%sS^f@d37ob!N+C5W`aVah~9Mj8_P
z;l-#!h6nVY99pV#h&2aTU}WkZ8TCfw6&TC3hQIRp+EGs5WeZ9-CTHINObYC_Iw<ko
zU9xpjT6ZQO^zF%u&thV%Q2~koA8k_*E8Mh+t7Xhr$aizsmi^O>gJkc(#)uG6iN$Eb
zAmEvM;uk~Nd0DVPV$%ZV$LO`^Ud*{7D`a>F>BMhk%_1-Nsm|qOD?)Uu21+&&zXM8<
z!KAy?J#dvL=^!yWUz(x<B7_gTk85F01+F+d94w09PiNlnVT7B+%R3i)Bd6|q>Q1={
zCp$7SS%?I&1~ui5=*QPrw@%4dp2SQO-lG9TJ^2XU)2#dAJ_~RCc?%t~ogc~!qYnMe
zqFsdx-a=D0Gcl+PfsPXzzhxF@Fa9{7WpY37M(gnvw>)Fpgj^#Lx^O-qbd5%#Ru5Ue
z(a{58hZQryGOPBNfupJ&qXHXPCrnh*t6H{f*=g`UKx5bp(I4SF-V%OOdf>;Q=CObj
zm89XDewX9FmrL+-6{Q66c-D&<^l0qOY23@Ujvfk(M(IZ;uM;}`UzcIXloDEJz)n&|
zq&ScF9{pjm;THF?m%j8VRo{0m&1g~McYJlZU6n|S09_Em$>v8tM+7bBKYX$n^euTg
zKW>nvIrRe3610c`WWZIIFdS%!gqQXT=eLfV9$J9Ff<Wi@jVlxjCjov<O}WJAo8S4H
z`>&16${6c0?uIvqJ9|+ZZ9!IfpjuHN$fQ-I*#MvT0YR@!%-SPEj8Oc+()|G3Op7aP
z4MZhvnhfj3jq}HLn$GR{vnAa{iTZ?GEO-!LCe^P=*|qOA&yo*B{=y5@tTQRZp<?3-
zCHZdg1dsryR#fP+q-C(b>rD$N8hW`&jeE^y2H7F+rlfI~6LX=MJ%UCVaPV?_3y8fT
z^4X0sQN-!Q$Wz<vW9#{{SJr0IH|l(v4F7!%&*O=|Gca{z)MTb$cF^QNeVKjS(c=@d
z-0xP_u**H@C8By;lb@(`wdS`7QF_VqxLMer-n(qda^HELv{8$UL-)%b{}`-J_GOb<
zm)}4@T?QRuW)05RPK$m<#ShBw?KxH+swXCQRnx~kCvkXo0kJI(YD_uXc*87)!*^Lj
z<(>)YR!90cX$2J(y3L}Z>)?4;2?qzkO0;H7kzy*aK8-N%Ix#`h^oue10o}D-BA3q5
zpzXkO>`+|I+J>UKwMSxB?Y(XB`-o0uqRcXa@iOlpvm~g9gf#7T<zH8CXD{<C?8qsI
zb92UaBZCwSj*F~b^||XP#c18%UV)uiZ=_?<!nJ^yp<z4K9g+0aF=AKpSex}%QghM?
zV{fzL^`3X&SQu`66=<2{uKud{bZ!{Q@qlxuO}&12PU5AF<Y9xX+Wl=t8LuyKofc<2
zE?*u)z$a0>z*ka{)LiWEhwTjeKCygAYV>vB4Y*j-*%!h7&~MEf?Dh(#g2fBZ-VaO$
z&NB?}|MBW?HK#P~kg|u(IGrAyJ}h{+hh!tcY2ym_D~(7;(v0%Gh%YjGt116vdWJ#s
z!?D2A4&MA6JGAzidW6);_MKc&y+RV|y{a@1S@uW=-|zaS^A1Z`)k^*C?N!A6##c_q
z_>>$ovy?euSyRGP+=D{O=k4}e=&z>xpqj&Glb5k?p04q*S&Wm&O+OD2=gC(aO5BH6
zo+&+1q<ps{>^YBuPyV{V`HhjEVDx*C=r&?R;aZy?{ytWIs_$&+frR8^ZEP@xSnDDo
zZf2<;K%vCVT)cVq*ziq-Vq{bo&hFn}{XTnHRie8#UJ#=4a$<&g1~Zw;jp6%H4|tJz
zRKmB`ur@rspP^}{(EolY3<=?!!<*0&?{&1;-*5QKs}1h=O@6+$`JMWww)DPh)p%)Y
z%=n$rl}c00N=$beGg+!qSDCXm{}xi0mp|T{4v1o;N;18Nt0|MoI@&RJA-mdSZIt8+
z&f@I}LX^C+{u^l<Mf1#G<*h#R+o9#`Qi(KgFT`)R4_=za`EB+k+pcF4<%PJ}?WK}y
zsikAGFK5k5S988Acft!YypEi^yxx@WO3P{=`o2h(r+%gV&FNGp7uL`mHZbr!+kS=s
z?%NB5_*+Ic{G4N@Ns?<9sGJ#u!vaA5D{=;~woyb3lP|h&yQVNGEPZe>;5kFWkDGZL
zx#@S#zLT1Nta%`w)llgke(-10mZ=#V<HaX7>T%g;+x7qc>SCS>kF-q^6E3Z3^;j?Y
zO>pT&iTUY($?Cfb(={9wFHI}_e;xFu2=9%q9$cfcFUjJzpN9St$>592BUSO=ud9Cg
z$tUu;-MRfFN_DEsy;RlPx$m9#<m-OR2N?9Y#sle_!HN0Z%p$jz$|GrfqFxd<Q&3n*
zo@|=g{q)$iZX)uNH`v@BX%!k>Hom%7@n3egTTpM_sz@9zr*wU5aHZlpu-u$VSC)}R
zuGM|zU(|a2_^fWz>+0`OX)~K>HLsSn^!D|M;fSMx(tYKsRkmfbH~pga0ODj05$n0~
z{`yXV&}`{D28+F2gO?;kL$fPfI9qIEg?9aKleHx7tIE69eGN>XVm4k3T|CHAb@v3p
zmI#!x^5T!qq){w(M>FoehF59&ZL%MYm_st3pbttp$*d}He$YqO+_37+x0WpH9Zwxw
zoKrR@3d!5U4){>+Kk5a{CWNU8dSg@^x7xp%5OI8IepAuOvb$Y4)^m%QY}Iv4JGLth
zh5Q%5S#rI?_uLq}FL3cP3Om0`xx5JzN~cD@warB1{l=s-JsCjXL}5!@?Y0vRwiHQw
z9_}yNIhF@kh~E?o&GuGo;~cO0{`?QU;T8b9C=WXL4u+uNNqnU8cm-HY7{#nS>kqK|
zs9;Q;arn|XGxI|$;EiEsn#j#3l4i>q-B*Af3%JJ7hhWqimQKeuyRpPv_Q(~GRpN{=
z&p70`TKv*w8H)75qUp6A(ftMg%P@#mTNOMgTk&Q4Ijw#_mf(8NIQH0CFZ1cTgkN_%
zQU*k?79=Mbw#R;ZFeX3>9r9^Gu)g-wV4Rk^9kP0M(Z}sReCnU-pD)qccCkFU*krYX
zSzV{&;aAJVEf}uX>#DTou08Yb`43$uBHkJ{XdYeI9P={dKzwca&v<!Se#fMjz(GUa
zo__mSi!VTzD(k*;_esnK-*M)}3J=Wd_Vl2duVrDrv$P>bxl{$XRV71TFJHJm?;V`q
zS<R$$yql(O1!gs@{o0_=wI=|DsQ|y7T93bX^?P$oi>&FJ@Ov7)Hz)**uB5#EOkhU|
z-CFl(ZQ&FG%ltldf>nz9_U{8zy(`oEI3?_k3VSr4*q9UJg+2|72ZWqXKlf3L$QZW+
zb0f=OyugqCF<isKT`CE+3r#Fi&tph3nMz)l)aP4!h*UN)wi)ak%hoYj2`%isy*fMM
z*+Z^%Ca#Ss*HS}D8G}V$N))OI(W>pd(E!q5h~XA8ci`U5ml)AhqF~6(e5#3|qi_4j
zd}hZL?y60$0~rN|{D1KUEeWjmAC3ys2fjHaEm+qt%IJ$>;*1R@<6^pAksX;POv-5I
z-`@|RoF!cmbcx6ubRFf>a|W<WBqOcKupIQ`F$0;@eWYyyX8Xu5fsx&0gO2WGTF$K~
z8M|OIv?IExqSeXhdU>c$G88>^2l^aEC%NOVy5=8sOOI;nIv=gO0>CtZjL1}j!!@HN
z-oxTH)&XofEX;Nm=cvyGN-DKd>F%ZxzP8fS+*xAvpX1`PX2=5j+#KiJ{^i!3NsCO5
zocn^Me(cI+g#{#;@cSW+mqN6!`YwYK+M<(CQS-vMpLoR}nR`=YM^rgP8ink>cb~&Y
z{_qhfzcMi;Fs*1%V&a=3O%g%A9HVM4?kK&JRi@BkGPtL~qqNy=Zd$lhN!Dg8Mf~Mv
zUjMUkPBq8c%pjSkk=1_xu9o330j<cSik?lB%I9-le(S@+!;e9TEMe{}DRp}L;0&X?
z>q;o(WYV;j3%jv3>-B}r6(lSrxeIT`MX&BESC5U_(Uil{P~mfTlPdRR>Gcpr=B>Tu
z*WM2$Y)Bs{-La^gEikC2&W^ZrCjX-eOBA%>e<=i*WB?cL$=cke|6oKf`%1#My7+B~
z0vV>Uq2d>4`zJ2lu=5jcP-)untga|aCs1O`VJ>#Hu7-<NIa;4~tXCT=Icx@s2kDFE
z=CjBr-X^?3bWtNA#-JwJuIOtEPMBqGWP5(dz9HKuhd143OY-zuhvdEP)igYdmfdS)
zWPWJNR|c)i(eL*6-f$Ccwf$V+J*1B>Z<@Gs!k9f^SB0gp@5V{FJF2GFQflE=^RZlf
z;3+8Th9aG*4;3`Wf#_sjTU-zwb_Qv!7Cpg8)YG*<2S~dUngS4Lbrb#{!bZdEwu<`Y
zeKJPe_wF3-)=uE$;iq&_lVkZ}4gE>GjSf7wUR)3nxx};V{=O~AD{n~nRDWxUWHuP9
zIG~b~wKlTilGIb?ca)Pa$@(BOsB}RdAjSfpm3P@-@Nzo%9B7rBIyhR{lw|vq>y1>M
zoN`GJ*|uXb-lmn4id~n1UDn`kuE-MPJyvn~gnQ{u2C=;lj{4pYTNNSQQLXxP>*=F6
zPq3ISwjZq-yHhW{+$hbHi8uEQ?YsuD_ku4j6$D2tm_~MH2SxK^w5-c(OqeFK;KAYB
zk=2F)E*~L!-n{u;)-|ASaCka1AoX_~8H?fD2XSm405f$d4cZMp4rp@_IkhWM*DaOp
zFiS$~SEvgqRh!U+Bu<W`oX%Gm(faOUbzJT#U*+|kw<iRC>NCD=<aF}BQ-9H_o2gmv
za;xhY^WO||*S2eh^a6MGPkS#&yBC{OUH6pPRw@#!(9P{>?Ay_F$hbXPd!jh8QS#{b
zqm!O&xu8{?y3~D_Z?`%Fr$5l=To=M80Cutwe1dJ6gCJBKhf?BZT}v*g$YhT86a-cO
zuf~($JkcF(rFU+9gH2chW<wBF1?EHe1HZzopC~l_O?>Se<}#}451UaI_+@1ufSe0o
zK0}9XrYVQBeWcoJ4@^;81g-m2s*{X7KT3?MGIu>^$;|PE4h?c=Dl7Gu?u;QJyMS;0
z*RnRtYhjB%*Cie~tiN$&Lv0I;$bBrgbkW4KcHZH!#f)8;9LIlwh0O*K=Lu&K{-gK2
z8Q3xOdg4V5F|V=dDuRmuI$j$gfj?OVL}lyVy-hnuokQc;@p7c}h{3ISV=5F_52c4*
zd!d7L7mnV)V4+b1L!Zn-X`zMYzDrxgnHM{*>?3%C8L46Qoh~E{3My}ZzN&}wry*a@
zv07PchB`}X#xt;uuinc<X^~TipQPC@_pkRAKk|Nc%4Lsu;}XFF>Md(^Ds)nwpRqKC
z3s=*c#Qr;qy3C`Ds3O0GDhoeGe-Cle53ZuD=6CMiEkY6GglV1Bdoe>~8ra3HfKokp
z_w_{?yH})S_5>E&$jhsK0LA_o^e<$N07e}WY7x9PLvG)G2eJ`1iD~h}S>2sCww|1p
z9eKBmcs%Vmmlp+6=iTPi=XWVd)(>6!TdHLnOwAwMD*aNc_KID@(a>k8QI4PUBl4Gi
z81x%pwC?H+qBo!=`e>hH{epwzV+&Dyr_v5%_VGJRmn6PAfPzc_dJ&X|@GW9p>9783
z=-Nx6ux4H9KbiVwKoH6XDgcjL<a|aTe!tKC<;&RQer<-w{u=tk^alOF{q2~xoSk)5
zh&<uR?2r>cqo<_iZ@64%Gw~phGR(%Pgx~%YrJ?z)rdD(+HRqQSKoo|9s+In4XRvou
z=gCV1Ek8RN(t!Ioo`0bc4!PDRGD;2(aee{GbEj16p$AR_#wi`h#jck_1`z*j&1Uf}
z>{ObIPJVuzT?-R>^S*<3&<pJ68>X^6Mn|#dW!jEsrHM)plSr=N`YVo17D-!WW^PCi
zY&u!*Et>}@T#?*?e%D-}m-=K`3yuH`arL`(0H{FAni=<AG;jukefeQ*2+ww-9y_$X
znK;_Ll$8WE{^@@|V0Wabk(iATT>vQlzHrpV%cMJqPBf;bib=i6Jf%6KOM<SptD}cM
z8I4Y`Mx+KRk>ZAcu}pls_jogiVNlA0!_3TV??I|}&)@9@=SU2LgFI(BM&3H?rJO9*
z7>8up)fW;EjIV?n)p{v+xRC+r)f({2_z#ONXWmw^Cm28M`mw9{7MV4HwYUZb$9A$H
zR2n&$S)E}0b3uJ#ypG`*nE*^D0#=^H42*$oM+ce|F_<ew3{*{g``9#{4Cygp3h%_i
zfrE0=TeF8;hBEARF>hTZ#26~MCqONqMV~qy%fODw!<4+Px`l>K(=LCX`qqM3m=Ao8
zQ`>~8^;@Y3I=p|XV;LFCny)zG$!3BH?RUsCrUOf{7eq@k?3ci;7?quARzo2_m?>u`
zgo((Mui!|=TpEj5$CwWmLWZ=c+sto6IJ9p)!Z(CaNYv)Xg18I%cO=ZcMNdzm-WnVz
zTd0ygr<|2U2s(zoQa1Wj@H{$;E;hQbeSZjV8w%cHXf?#3>fJgxT`_QdE!cuk{IE*L
zJx|`cqr`K|uDmIVLEUYk{zWp=SV40XcO?pkN;ne}3M1HdqEL1-Ge3-3T5MYj7)n0*
z@yotCqrA~fxB1e0rQb8rQ074EP4#Utc;ESiRjX#7B-(~&kJu$bbbn5L!%PQO4~4hN
zVw7uX{Ij6kLzg|^mMf!1eWEmA=eD%~lK~*Za}F`aDO#6y?xBZAXkh>tML5(``(VYp
z6<HccG{Gf;rlx4$m;i;<1fLB2A9#d+xL|6!j0pT6cBQ$c&)^_CGGJF+i{weFvwMGI
z*YA+C7Jx#DL}w0obgD6N9R^NE4hL4}HSl#sTqQP0Y`2C}*&;ya#uN=OE&0U7>wa3h
zMlZNL*f>AAPd0jPxaZ%~w4YVVj|dES&ENUgcL%(xqGog#V;1$}jxFXc<mbAqxtTck
zZSTVMa#q#?@nf|)Ci1E$^bcn!4jOfQH+$wEqE}(N_3fkH<0_x;&5J)sj);6-7W3`d
z`mlzN!m$ZAzTITgr>Mvs=Q+Ks%wzqj&guEi7caMzetSRkJm=iCtUm7WUBRa9Qoe<$
zJdrN_W>#iKxj)nO9p>6vA7y)8g=6^u227&0{1jm_MJMAhAoWk><;kOVkas0vW!vw)
zFn9^Ges-;5NDk*|C{xu-(cbp8K%JG7yIE(KH2#j$BDmLQPeU2gj|;xN%}=u*2aZ24
zbpA4`ZRgc{%p6M;v*MHPMn|WjrMFW*M+YKb4bYe#9a~1V^vC1Tdhx3xN?1st_(LO0
z|L3o%>KG2$Ex#<+Wxe<AKmTm?EoA<FVu0??$BbWcoRWTp4>2S+^Ef!)m|O_?9_!4M
z^vgfrm~zSW?3HWy``f4-UVE|s{d!=#@CePHuMf(yYHWJ$^9g#{yZ7%4DdK2Shv<by
z=ETVMUAr35r>lD`x7^!zOzhSGkkjvw$A-p$Whnvj>^l~^{O8X_Sy@?a+$Ge-r@G9n
zD{2f}#9tde&SpPmXC6eJ1|vg7?%v3=cH)itmXhIO;-$z+)y||FTnAhIQB+ja!C;mv
zm_(O~nkW%@>%iw{x7=Y1crbW<9p)Cqz+>t=dN_4x%F)QGR^OK}{|BQz(qK-zea{{(
zG(%SWRu~cE)RpD5@TG5gT)?=dB4uMmmYm&Sb*qnWuU2qx^W$uXfhXCz9gMa}a>8r-
zHC$a?A01n%1T|s}#xPu&8ZlIf5IuqoZuIEbs!LTeV*zl;0HcQs58V6a%~jcLRy<O+
z^2_K~S1;e>y7qfWhpG5$XYnQ`#+PzlTpJ_L%0`;w%|FimpAWA)V&`;=eiOp0@<Z2O
z{%Ex09co?77dZ;=)8^+FZh_}^H!3RGJT<!{bFuwZ#^PDDf9o?%tM)yR<sfsSKb4l=
zF6;#KEPJKEm~cc-=R+2L<j4`H3Cl0fkDzt^4<A}DX0T##kC0X;j}VkmK|#R*&qXhh
z5;Gi&pDHW0u3UM6pLmKOOK<SGLtR5723b0B*E^-8^ocGQU2noQu0saJ&(FVIRP=M^
zaJxD8#*L|H7^PO4`tWf6<of$)AN%&s_$)n6Gpp~9v-9v=%~Bs2SO~K_y6px%a~@ZH
zGbrevmccxOD<ok~IhU_b7}LbBqz!5G45Ut1a2hM^yUH~S_K(nfHnKO_?x(YG(=fxc
zx_R^F$s<RK>7K$8AQmF^OTakcN!<uORT{#95Z<_HWmF9)Vlu;npPw9IU}a~YooyMp
zk9nXVTAsQDn*k*yC1H*>bY2c(Cp$xa2m+HZtcw?Gp&iAUnNM|`dv6H`>0Dx;YhWHW
zH`52Yr_#X6(0NzGj<3gWyBisq2#;gu`9g(I60RQXwR!UC8pr0{X6cvB(s?*_?Zo$r
z*XXua>q?8i>cy$@<j?c4s-UyCV-$`k>k%#Cxbx3HCiqaTi;M753gjcRY8_t&ZV)R$
zeNfSG&K%aoy`-N0ZRk~vV6Afcbc7RtgyZ=0Y#vtkTQOs4%q=`T{60#!`$$!sgUEo@
z?HDSFUmhNA3NNjJhn>qh9*$Y!_G1B$giWNKeO0L-ie!1jf(tn#`1ESejNjV-Fu`~*
z-K&N4I!;d6uqn-uvLCn`fkvGFty@V~0&u;d<)HHgh?GRhx-_ij*!>Qjd1oF;D19;2
z+3}RV{*doSLofPSAF|6Nj(N$kRE#r{jfQRRvo-Pu7uf?1y1OOACp*RuL0pp5mXE@(
zxTM5+{)z=m@|*Ek>8@#Z9mZ%k>@WORS2!F6__x~f0zQ8JXU%F_x_PeI;LpY5d{OMP
zI(rR%@?jAX3E;J<gkC|95sveU%{fk)$Yh@2WSY4!aWMYlv;)+0+k}NHGcC!dytKRz
z!(l!}<mh*5`|JE!1YSUTldYS&?}3y$7=dA-p%l!O?GZr;z3J~?gLW&wVz{0|1p19>
z=p`P7GfEuaMMTsYtY-<EE_uvPK8v(swD5QQ_P@KKk(-f8MZ8f`ypgFQMJkdj(k%6|
z83N2%yMk6-;gN>N5DOg9dvj&9d7{eYZg)>l{dXyJN62Xgrk+D{wp|>oOZ0!((I>%#
zMca8(D7Kqw2}hNHP9E};H(fV{#GQVfw;yUfjrX^vDE`pnt$i)%a%RD^qWHsyAJ`F%
zNNDEwo2JJP%q7OhSEEwof8f&o3<O*uDXG)g*v{zV(1-{MGc!&wl9Ot%x}#do&W#+6
zShbSBtN$$DM%71=WiGg89Odx4ZZ<y)nFgPtfLe5NjUB4$vi>U#6K3OwFJ4R_*?XEL
zng9FqP7Ev&`k^QpowX>4Pfk9LjQ3abMEe-bM^2DprQMX_=(A_f!uEHL<@agBhN->p
z@ZjTuyvGF&WCL+K|9xawH9i~vE}!O=8g|)Uub3n*e{Hp!9X}AX>40$m+%0NNb8PIs
z?!{w@ivFgXa%}Rxbk41E%h~E$XkNtBcb}&+t-8?son+3|pVmv?lFgyVx8CDh4-P&!
zr`wSt2Yjs>*-Oq#Nav3srwOxaJbD2sZwm3*$Qpe30k2_Ak;0*r#Ic<+GP<yDqK*`F
zv~l;|k2kK*OD-87f~QA*ZADh$G%LA&N0lZVf;9H-HrHP<Ap=wbx3&R|-z|#djezWe
zY0jHLYTf#-mM$NTLgr{{l2s+!!X~U6kLwvg63NtYzMvm}sqOBb*H<xGR#t}FmuW+@
zIL!S4O<!di>a+oHFn(K{eQo_1?Q_BTH%tQp|JS*EP=ahkGN3X6OTrv&Y3!_f#H>p3
zGos#l9-b78-cTg&gb*30pvK#F;J`T`zs;Bf=^?O&2{q)Mmca&?3fRw0+2oEloDGeE
zdvX$HGrt*p+437@qdCrJXK0K7(#1jYYmUtQI*-gjl1+InU-PNQ_T!*~OXVx@iS$Us
zf&J!hAhzZh|ILXxGe|jp<NTLs8ncAX9+B&$s%j8IfJ|TGjg>6j)_cV<95fJ64QL6@
zQ=Pt9pe+!wJEKGH0Pp{G`v#U)_e4tS`n6XWBH{{e+ZO*~9xQ83=(=Qr?d<Hbz)nOk
z-7X?>A|xbafm_Us+~EwY2!$-KfhnKcU<R<(Hd~j(BqZ(6)Gk|W10hROM@Q%4EZlMG
zU`|x^4V(EQ9xcM2mXXnuNXg4XRU%#`U@2a~O)9k@S|vckt{D0!?&EBNUlIXi><8*T
zwm<D!S7Xp#Z6Hlj2`1UU=LwEfxix>!)<IfxyX<!if73RMxG(shuh|``|NQ-<=q{RF
ze}AOzlx_#zpRYZfZ!+Zk{cCkKo9q8;@>fi!PaB4{Zw_tQ^^H%p2AeYE&~j(|%|T^-
zLY#wY&bODwW3gKI#y;T@O8i<qsukRvtEm>`?5N*eGZQ|{^B)&+Y5jNV>U1z8&MB#>
zY1`OjI8Xd^Dlmo;ac5j@0!kk^ORJ(74F>luH;5!2a1p>^!P{o@fd?n?=ufc{GpwLb
zA3vTzRe1$GD42q(vOh>kPge($nu_^PkTjlgdbvr)Q5`MNMocQT|NZ?0$*DuZyLx=@
z9tR`OK8-+5&H2yGuWg=!LNyL#!RcjBZ|?y7AL;;pbd0r9sh<*`!fm|UWHuxv&hBxb
zy4B9is<413XT`m)E~D}O+GB_Vv51$+Fu1}Q=EP&yRkX4V^kdZSkw{cyVeq65BG9L@
zvN3L`?S{XjZf(SVegW`+tXPm$!mbE6k%bEnP7NI$oRB+;i;MelZs&h@1<KI?+J8KF
zh)$M^7`^<$wnvE^dgimpbIBe_LMpt^uJ?FY)mxiyN8JD_>9;_AlIyNl6rN4}Xsg+X
zlQ_pgy{FlIHM&&xfxXek?6}r0qb&tsavb$r#m74Wfg&Yy0E2r`(2JA9nVX~|I0MX(
zA^j`NbBdP+aEQFOJFcr6MLh9kgR7#XjPXJE#{|-d1&yv3C`9x;p0^~JAQDLS*&m2V
zNYulY*U|Dvb5>FWrW?}*A3!<*3-Pzps2p?i^bvh?^Wq@0b(zlwQdjq5cEa9@>NN+a
z>qbU;^+vNHCvo;am*48+^~&tvJ2KOV6SJL>Df1~J^J;*z2Y&I6H28IGc#6VU@iWGz
z5#Wb}X3>&uJ9Zq02B#X{Q1)QcBtf2mvBn|jyKGTBlVRkIMK%tOV;KjpaR9<dD!_#S
z)GqAv5;x$uR`Ebu8?}Qv{2#h)v#$V1Ay`|qF@WJjD8A3&7~JPJ^MX^-3Mi4o^ou=t
zm5fyEa$ZSCy{PN<M4pw1WO($lrZusqwYvDM_+|qnGe|ssWToJz^@AB!trE0+JCK;9
zp`{fI5r$moZ=CN*pcmZ?4%WteBynUFtZZz~c|OSSVfv$~ub+xKgPcb@cd9V?etRN@
z{dfv1Y^M8;8lR%_KRaq5mAryyq#W9is+WWg#&%v_xevjc!uh;mL$WY?a(bU6mF{ZB
zR8BnK7Bj1QGj!69+9e&elViU+LNP$H5Qk^rY-#tlO+Ps8%|N133$7ZE`wINb6hcW+
z1w1zK$&<s#Je^VUOF%-Days)UbnhoXCHoCOUfo|&m^^sq`Fa0X)dz`~LR1kYRhxMk
z*Xj&jDdD#|zJwp`)%yMI0Hsq=%{W^M?c53HQZBn?2krm2MawkQafGCrm*AQ$C?s^}
z%1}if`_41n_Qz2U!!{)mu4xLP2o4A{oBCnAuLY3_-$nRHg;vgsCtV%>N~j_Bc!P?c
z@zG>PN<XC4`c*${voFb*Vzd>iNmE?sDbzQdQnoRs?bsXjONTcYoVdA!=Dr`T*_lvY
z_D$#Uf9${F%ya6`<?d#r?>HbH-b*vGDN@;vi$X3HV)je@pP%7HHX`{+C@(KpI$#+o
z<MbRUt@Gqi3TiaKEJ6sSC_kc&O{$_(-oHPLlr0H45*Z~zT55`)QLJR4GH8R92)QVn
zD~XsWP4Ydc=aXs(O;79aVNX~Cr$SM=4Q{QO!=R{DBQR>9hJywB#pz!?VeS>*@^Qdy
z#L`kB<tdWcERZ-*7P<WY+{BaJ#sMpTuQhwEq^z6;%d4=vcjK|i-@vmf%X$1#?#O$F
z=>I&kEnBqkeACF89c_w~01*+~undqHV&G-})%+T~HJ@i@k_S3I*rE;2<I6~xxsv&b
z7${u`iuJBEH$YG`cy$QrckkVc$79EuH81jG4vhhU!wh$<q@prD7lFf<ATT&_$v#14
zMzlt0!*BQW(t8Y48{ou7%&DMTdMA6t|1;3x1o4`{Cphn?8_TC(voRSGzXNKtRE*0R
z>8TbUXH((durHMOCjSkY2W7Es|Hejq&GUc8$?<Ce#g9a^uL<2G|4Q+Qs(k8U1JC~t
D0x@tb

literal 0
HcmV?d00001

diff --git a/docs/api/tools_index.md b/docs/api/tools_index.md
index 0947beb8..aec4577c 100644
--- a/docs/api/tools_index.md
+++ b/docs/api/tools_index.md
@@ -556,8 +556,8 @@ See [perturbation space tutorial](https://pertpy.readthedocs.io/en/latest/tutori
 
 ### Dose-response curve fitting
 
-{meth}`~pertpy.tools.PseudobulkSpace.dose_response` returns an AnnData with one observation per perturbation and dose, holding its distance from control in `.obs`.
-{meth}`~pertpy.tools.PseudobulkSpace.fit_dose_response` fits a four-parameter Hill curve per perturbation and stores EC50 and the other curve parameters in `.obs`.
+{meth}`~pertpy.tools.PseudobulkSpace.dose_response` returns an AnnData with one observation per perturbation and dose, holding the mean expression in `X` and the distance from control in `.obs`.
+{meth}`~pertpy.tools.PseudobulkSpace.fit_dose_response` fits a four-parameter Hill curve per perturbation to the distance or any gene and stores EC50 and the other curve parameters in `.obs`.
 For preprocessed `adata` with `perturbation` and `dose` columns, a `control` group and a PCA representation:
 
 ```python
@@ -566,6 +566,7 @@ import pertpy as pt
 ps = pt.tl.PseudobulkSpace()
 dose_adata = ps.dose_response(adata, embedding_key="X_pca")
 ps.fit_dose_response(dose_adata)
+ps.plot_dose_response(dose_adata)
 ```
 
 Assay measurements such as viability can be fit the same way from an AnnData whose `.obs` holds perturbation, dose and response columns.
diff --git a/src/pertpy/tools/_perturbation_space/_perturbation_space.py b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
index fbcd6c97..69abda0f 100644
--- a/src/pertpy/tools/_perturbation_space/_perturbation_space.py
+++ b/src/pertpy/tools/_perturbation_space/_perturbation_space.py
@@ -12,12 +12,15 @@
 from scipy.special import expit
 from scipy.stats import entropy
 
+from pertpy._doc import _doc_params, doc_common_plot_args
 from pertpy._logger import logger
 from pertpy._types import cast_dense, cast_frame
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Iterable, Sequence
 
+    from matplotlib.figure import Figure
+
     from pertpy._types import RandomStateLike
     from pertpy.tools._distances._distances import Metric
 
@@ -670,12 +673,12 @@ def dose_response(
             dose_col: `.obs` column with the (numeric) dose.
             reference_key: Control perturbation all doses are compared against.
             metric: Distance metric passed to :class:`~pertpy.tools.Distance`.
-            layer_key: Layer to compute distances from.
+            layer_key: Layer to compute distances and mean expression from.
             embedding_key: `.obsm` embedding to compute distances from.
             kwargs: Passed to :meth:`~pertpy.tools.Distance.onesided_distances`.
 
         Returns:
-            AnnData with one observation per non-reference (perturbation, dose) group, holding the group mean in ``X`` and the ``distance`` in `.obs`.
+            AnnData with one observation per non-reference (perturbation, dose) group, holding the mean expression in ``X`` and the ``distance`` in `.obs`.
 
         Examples:
             >>> import pertpy as pt
@@ -705,7 +708,7 @@ def dose_response(
 
         treated = grouped[~is_control].copy()
         treated.obs["_dose_group"] = treated.obs["_dose_group"].cat.remove_unused_categories()
-        dose_adata = sc.get.aggregate(treated, by="_dose_group", func="mean", layer=layer_key, obsm=embedding_key)
+        dose_adata = sc.get.aggregate(treated, by="_dose_group", func="mean", layer=layer_key)
         dose_adata.X = dose_adata.layers.pop("mean")
         _carry_constant_obs(dose_adata, cast_frame(treated.obs), "_dose_group")
         dose_obs = cast_frame(dose_adata.obs)
@@ -722,11 +725,11 @@ def dose_response(
     def fit_dose_response(
         self,
         adata: AnnData,
+        response: str = "distance",
         *,
         target_col: str = "perturbation",
         dose_col: str = "dose",
-        response_col: str = "distance",
-        key_added: str = "hill",
+        key_added: str | None = None,
     ) -> None:
         """Fit a four-parameter Hill curve for each perturbation.
 
@@ -734,10 +737,10 @@ def fit_dose_response(
 
         Args:
             adata: AnnData with one observation per dose or replicate, such as the output of :meth:`dose_response`.
+            response: `.obs` column or gene in ``var_names`` to fit.
             target_col: `.obs` column identifying the perturbation.
             dose_col: `.obs` column containing non-negative numeric doses.
-            response_col: `.obs` column containing the scalar response to fit.
-            key_added: Prefix of the `.obs` columns the results are written to.
+            key_added: Prefix of the `.obs` columns the results are written to. Defaults to ``response``.
 
         Returns:
             Adds the fitted response and the per-perturbation ``e0``, ``emax``, ``slope``, ``ec50``, ``ec50_se``, ``r_squared`` and ``ec50_in_range`` to `.obs`, prefixed with ``key_added``.
@@ -748,18 +751,16 @@ def fit_dose_response(
             >>> ps = pt.tl.PseudobulkSpace()
             >>> dose_adata = ps.dose_response(adata, dose_col="dose_value", embedding_key="X_pca")
             >>> ps.fit_dose_response(dose_adata, dose_col="dose_value")
+            >>> ps.fit_dose_response(dose_adata, "CDKN1A", dose_col="dose_value")
         """
-        obs = cast_frame(adata.obs)
-        missing = {target_col, dose_col, response_col}.difference(obs.columns)
-        if missing:
-            raise ValueError(f"Columns {sorted(missing)} do not exist in the .obs attribute.")
-        doses = obs[dose_col].to_numpy(dtype=float)
-        responses = obs[response_col].to_numpy(dtype=float)
+        data = sc.get.obs_df(adata, keys=[target_col, dose_col, response])
+        doses = data[dose_col].to_numpy(dtype=float)
+        responses = data[response].to_numpy(dtype=float)
         if (doses < 0).any():
             raise ValueError("Dose values must be non-negative.")
 
-        labels = obs[target_col].to_numpy()
-        fitted = np.full(len(obs), np.nan)
+        labels = data[target_col].to_numpy()
+        fitted = np.full(len(data), np.nan)
         records: dict[object, dict[str, float | bool]] = {}
         for perturbation in pd.unique(labels):
             mask = labels == perturbation
@@ -784,10 +785,87 @@ def fit_dose_response(
                     )
             records[perturbation] = fit
 
+        prefix = response if key_added is None else key_added
         fits = pd.DataFrame.from_dict(records, orient="index")
-        adata.obs[f"{key_added}_fitted"] = fitted
+        adata.obs[f"{prefix}_fitted"] = fitted
         for col in fits.columns:
-            adata.obs[f"{key_added}_{col}"] = fits[col].reindex(labels).to_numpy()
+            adata.obs[f"{prefix}_{col}"] = fits[col].reindex(labels).to_numpy()
+
+    @_doc_params(common_plot_args=doc_common_plot_args)
+    def plot_dose_response(  # pragma: no cover # noqa: D417
+        self,
+        adata: AnnData,
+        response: str = "distance",
+        *,
+        target_col: str = "perturbation",
+        dose_col: str = "dose",
+        key_added: str | None = None,
+        perturbations: Sequence[str] | None = None,
+        ncols: int = 4,
+        return_fig: bool = False,
+    ) -> Figure | None:
+        """Plot the measured responses and fitted Hill curves of each perturbation.
+
+        Args:
+            adata: AnnData with one observation per dose or replicate, such as the output of :meth:`dose_response`.
+            response: `.obs` column or gene in ``var_names`` to plot.
+            target_col: `.obs` column identifying the perturbation.
+            dose_col: `.obs` column containing the doses.
+            key_added: Prefix passed to :meth:`fit_dose_response`. Defaults to ``response``.
+            perturbations: Perturbations to plot. Defaults to all.
+            ncols: Number of panels per row.
+            {common_plot_args}
+
+        Returns:
+            If `return_fig` is `True`, returns the figure, otherwise `None`.
+
+        Examples:
+            >>> import pertpy as pt
+            >>> adata = pt.ds.srivatsan_2020_sciplex2()
+            >>> ps = pt.tl.PseudobulkSpace()
+            >>> dose_adata = ps.dose_response(adata, dose_col="dose_value", embedding_key="X_pca")
+            >>> ps.fit_dose_response(dose_adata, dose_col="dose_value")
+            >>> ps.plot_dose_response(dose_adata, dose_col="dose_value")
+
+        Preview:
+            .. image:: /_static/docstring_previews/dose_response.png
+        """
+        import matplotlib.pyplot as plt
+
+        prefix = response if key_added is None else key_added
+        data = sc.get.obs_df(adata, keys=[target_col, dose_col, response])
+        obs = cast_frame(adata.obs)
+        labels = list(pd.unique(data[target_col]) if perturbations is None else perturbations)
+        ncols = min(ncols, len(labels))
+        nrows = -(-len(labels) // ncols)
+        fig, axes = plt.subplots(nrows, ncols, figsize=(3.5 * ncols, 3 * nrows), squeeze=False, layout="constrained")
+        for ax, label in zip(axes.flat, labels, strict=False):
+            mask = (data[target_col] == label).to_numpy()
+            doses = data.loc[mask, dose_col].to_numpy(dtype=float)
+            ax.scatter(doses, data.loc[mask, response], zorder=3)
+            if f"{prefix}_ec50" in obs and np.isfinite((fit := obs.loc[mask].iloc[0])[f"{prefix}_ec50"]):
+                grid = np.geomspace(doses[doses > 0].min(), doses.max(), 200)
+                curve = _four_parameter_logistic(
+                    grid,
+                    fit[f"{prefix}_e0"],
+                    fit[f"{prefix}_emax"],
+                    np.log(fit[f"{prefix}_ec50"]),
+                    fit[f"{prefix}_slope"],
+                )
+                ax.plot(grid, curve, color="tab:orange")
+                if fit[f"{prefix}_ec50_in_range"]:
+                    ax.axvline(fit[f"{prefix}_ec50"], color="0.5", linestyle=":")
+            lowest = doses[doses > 0].min()
+            ax.set_xscale("symlog", linthresh=lowest)
+            ax.set_xlim(0 if (doses == 0).any() else lowest / 2, doses.max() * 2)
+            ax.set(title=str(label), xlabel=dose_col, ylabel=response)
+        for ax in axes.flat[len(labels) :]:
+            ax.set_visible(False)
+
+        if return_fig:
+            return fig
+        plt.show()
+        return None
 
     def plot_similarity(  # pragma: no cover
         self,
diff --git a/tests/tools/_perturbation_space/test_perturbation_space_extras.py b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
index f89be83e..db084706 100644
--- a/tests/tools/_perturbation_space/test_perturbation_space_extras.py
+++ b/tests/tools/_perturbation_space/test_perturbation_space_extras.py
@@ -86,13 +86,15 @@ def test_dose_response(rng, categorical_doses):
 
     ps = pt.tl.PseudobulkSpace()
     dose_adata = ps.dose_response(adata, dose_col="dose", metric="euclidean", embedding_key="X_pca")
-    assert dose_adata.shape == (6, 5)
+    assert dose_adata.shape == (6, 8)
     assert dose_adata.obs["dose"].tolist() == [0.1, 1.0, 3.0, 10.0, 30.0, 100.0]
     assert dose_adata.obs["distance"].is_monotonic_increasing
 
     ps.fit_dose_response(dose_adata)
-    assert dose_adata.obs["hill_ec50"].iloc[0] == pytest.approx(10, rel=1e-4)
-    assert dose_adata.obs["hill_ec50_in_range"].all()
+    ps.fit_dose_response(dose_adata, dose_adata.var_names[0])
+    assert dose_adata.obs["distance_ec50"].iloc[0] == pytest.approx(10, rel=1e-4)
+    assert dose_adata.obs["distance_ec50_in_range"].all()
+    assert dose_adata.obs[f"{dose_adata.var_names[0]}_ec50"].iloc[0] == pytest.approx(10, rel=1e-4)
 
 
 def _assay(data: pd.DataFrame) -> AnnData:
@@ -111,11 +113,11 @@ def test_fit_dose_response(e0, emax):
     pt.tl.PseudobulkSpace().fit_dose_response(adata)
     fit = _fits(adata).loc["drug"]
 
-    assert fit["hill_e0"] == pytest.approx(e0)
-    assert fit["hill_emax"] == pytest.approx(emax)
-    assert fit["hill_slope"] == pytest.approx(1.4)
-    assert fit["hill_ec50"] == pytest.approx(3.0)
-    np.testing.assert_allclose(adata.obs["hill_fitted"], responses, rtol=1e-6)
+    assert fit["distance_e0"] == pytest.approx(e0)
+    assert fit["distance_emax"] == pytest.approx(emax)
+    assert fit["distance_slope"] == pytest.approx(1.4)
+    assert fit["distance_ec50"] == pytest.approx(3.0)
+    np.testing.assert_allclose(adata.obs["distance_fitted"], responses, rtol=1e-6)
 
 
 def test_fit_dose_response_standard_error():
@@ -130,18 +132,18 @@ def test_fit_dose_response_standard_error():
     fits = _fits(adata)
     fit = fits.loc["complete"]
 
-    midpoint, slope = fit["hill_ec50"], fit["hill_slope"]
+    midpoint, slope = fit["distance_ec50"], fit["distance_slope"]
     fraction = doses**slope / (midpoint**slope + doses**slope)
     log_ratio = np.zeros_like(doses)
     np.log(doses / midpoint, out=log_ratio, where=doses > 0)
-    sensitivity = (fit["hill_emax"] - fit["hill_e0"]) * fraction * (1 - fraction)
+    sensitivity = (fit["distance_emax"] - fit["distance_e0"]) * fraction * (1 - fraction)
     jacobian = np.column_stack((1 - fraction, fraction, -slope * sensitivity / midpoint, sensitivity * log_ratio))
-    residuals = responses - (fit["hill_e0"] + (fit["hill_emax"] - fit["hill_e0"]) * fraction)
+    residuals = responses - (fit["distance_e0"] + (fit["distance_emax"] - fit["distance_e0"]) * fraction)
     variance = np.sum(residuals**2) / (len(doses) - 4)
-    assert fit["hill_ec50_se"] == pytest.approx(
+    assert fit["distance_ec50_se"] == pytest.approx(
         np.sqrt(np.linalg.inv(jacobian.T @ jacobian)[2, 2] * variance), rel=1e-4
     )
-    assert np.isnan(fits.loc["limited", "hill_ec50_se"])
+    assert np.isnan(fits.loc["limited", "distance_ec50_se"])
 
 
 def test_fit_dose_response_unfittable():
@@ -152,8 +154,8 @@ def test_fit_dose_response_unfittable():
     with pytest.warns(UserWarning, match="'bad'.*constant"):
         pt.tl.PseudobulkSpace().fit_dose_response(adata)
     fits = _fits(adata)
-    assert fits.loc["good", "hill_ec50"] == pytest.approx(10)
-    assert np.isnan(fits.loc["bad", "hill_ec50"])
+    assert fits.loc["good", "distance_ec50"] == pytest.approx(10)
+    assert np.isnan(fits.loc["bad", "distance_ec50"])
 
 
 def test_fit_dose_response_negative_dose():