diff --git a/docs/source/ENVI.rst b/docs/source/ENVI.rst
index 8933a3e..88ab98c 100644
--- a/docs/source/ENVI.rst
+++ b/docs/source/ENVI.rst
@@ -1,7 +1,7 @@
 ENVI
 =========
 
-.. automodule:: scenvi.ENVI
+.. automodule:: scenvi._envi
    :members:
    :undoc-members:
    :show-inheritance:
diff --git a/scenvi/__init__.py b/scenvi/__init__.py
index 9c18a73..6f5fca0 100644
--- a/scenvi/__init__.py
+++ b/scenvi/__init__.py
@@ -1,2 +1,28 @@
-from scenvi.ENVI import ENVI  # noqa: F401
+"""scENVI — ENVI and COVET.
+
+``ENVI`` is resolved lazily via ``__getattr__`` (PEP 562), so importing scenvi does
+not import jax, flax, optax, clu or tensorflow_probability. COVET is pure
+numpy/sklearn/scanpy and uses none of them, and resolving ENVI on first use keeps a
+breakage anywhere in that stack from taking ``compute_covet`` down with it — which
+is what #9 was, and what the tensorflow_probability pin does today.
+
+``from scenvi import ENVI`` behaves exactly as before.
+"""
+
 from scenvi.utils import compute_covet  # noqa: F401
+
+__all__ = ["ENVI", "compute_covet"]
+
+
+def __getattr__(name):
+    """Resolve ``ENVI`` on first access, so importing scenvi stays free of jax."""
+    if name == "ENVI":
+        from scenvi._envi import ENVI
+
+        return ENVI
+    raise AttributeError(f"module {__name__!r} has no attribute {name!r}")
+
+
+def __dir__():
+    """Keep ``ENVI`` discoverable despite the lazy import."""
+    return sorted(__all__)
diff --git a/scenvi/ENVI.py b/scenvi/_envi.py
similarity index 99%
rename from scenvi/ENVI.py
rename to scenvi/_envi.py
index 51c7f37..3413798 100644
--- a/scenvi/ENVI.py
+++ b/scenvi/_envi.py
@@ -21,7 +21,8 @@
     log_zinb_pdf,
 )
 
-from scenvi.utils import CVAE, Metrics, TrainState, compute_covet, niche_cell_type
+from scenvi._nn import CVAE, Metrics, TrainState
+from scenvi.utils import compute_covet, niche_cell_type
 
 
 class ENVI:
diff --git a/scenvi/_nn.py b/scenvi/_nn.py
new file mode 100644
index 0000000..c3c9264
--- /dev/null
+++ b/scenvi/_nn.py
@@ -0,0 +1,142 @@
+"""Flax/CLU model components for the ENVI CVAE.
+
+Kept apart from ``utils.py`` so that COVET — which is pure numpy/sklearn/scanpy
+— can be imported without jax, flax, clu or tensorflow_probability. Importing
+this module requires the ``envi`` extra.
+"""
+
+import jax
+import jax.numpy as jnp
+from clu import metrics
+from flax import linen as nn
+from flax import struct
+from flax.training import train_state
+from jax import random
+
+class FeedForward(nn.Module):
+    """
+    :meta private:
+    """
+
+    n_layers: int
+    n_neurons: int
+    n_output: int
+
+    @nn.compact
+    def __call__(self, x):
+        """
+        :meta private:
+        """
+
+        n_layers = self.n_layers
+        n_neurons = self.n_neurons
+        n_output = self.n_output
+
+        x = nn.Dense(
+            features=n_neurons,
+            dtype=jnp.float32,
+            kernel_init=nn.initializers.glorot_uniform(),
+            bias_init=nn.initializers.zeros_init(),
+        )(x)
+        x = nn.leaky_relu(x)
+        x = nn.LayerNorm(dtype=jnp.float32)(x)
+
+        for _ in range(n_layers - 1):
+
+            x = nn.Dense(
+                features=n_neurons,
+                dtype=jnp.float32,
+                kernel_init=nn.initializers.glorot_uniform(),
+                bias_init=nn.initializers.zeros_init(),
+            )(x)
+            x = nn.leaky_relu(x) + x
+            x = nn.LayerNorm(dtype=jnp.float32)(x)
+
+        output = nn.Dense(
+            features=n_output,
+            dtype=jnp.float32,
+            kernel_init=nn.initializers.glorot_uniform(),
+            bias_init=nn.initializers.zeros_init(),
+        )(x)
+
+        return output
+
+
+class CVAE(nn.Module):
+    """
+    :meta private:
+    """
+
+    n_layers: int
+    n_neurons: int
+    n_latent: int
+    n_output_exp: int
+    n_output_cov: int
+
+    def setup(self):
+        """
+        :meta private:
+        """
+
+        n_layers = self.n_layers
+        n_neurons = self.n_neurons
+        n_latent = self.n_latent
+        n_output_exp = self.n_output_exp
+        n_output_cov = self.n_output_cov
+
+        self.encoder = FeedForward(
+            n_layers=n_layers, n_neurons=n_neurons, n_output=n_latent * 2
+        )
+
+        self.decoder_exp = FeedForward(
+            n_layers=n_layers, n_neurons=n_neurons, n_output=n_output_exp
+        )
+
+        self.decoder_cov = FeedForward(
+            n_layers=n_layers, n_neurons=n_neurons, n_output=n_output_cov
+        )
+
+    def __call__(self, x, mode="spatial", key=random.key(0)):
+        """
+        :meta private:
+        """
+
+        conf_const = 0 if mode == "spatial" else 1
+        conf_neurons = jax.nn.one_hot(
+            conf_const * jnp.ones(x.shape[0], dtype=jnp.int8), 2, dtype=jnp.float32
+        )
+
+        x_conf = jnp.concatenate([x, conf_neurons], axis=-1)
+        enc_mu, enc_logstd = jnp.split(self.encoder(x_conf), 2, axis=-1)
+
+        key, subkey = random.split(key)
+        z = enc_mu + random.normal(key=subkey, shape=enc_logstd.shape) * jnp.exp(
+            enc_logstd
+        )
+        z_conf = jnp.concatenate([z, conf_neurons], axis=-1)
+
+        dec_exp = self.decoder_exp(z_conf)
+
+        if mode == "spatial":
+            dec_cov = self.decoder_cov(z)
+            return (enc_mu, enc_logstd, dec_exp, dec_cov)
+        return (enc_mu, enc_logstd, dec_exp)
+
+
+@struct.dataclass
+class Metrics(metrics.Collection):
+    """
+    :meta private:
+    """
+
+    enc_loss: metrics.Average
+    dec_loss: metrics.Average
+    enc_corr: metrics.Average
+
+
+class TrainState(train_state.TrainState):
+    """
+    :meta private:
+    """
+
+    metrics: Metrics
diff --git a/scenvi/utils.py b/scenvi/utils.py
index 8b04f56..3d98877 100644
--- a/scenvi/utils.py
+++ b/scenvi/utils.py
@@ -1,148 +1,22 @@
+"""COVET: covariance-environment niche representation (Haviv et al., Nat Biotechnol 2024).
+
+COVET's maths is numpy + sklearn only — ``sklearn.neighbors`` for the spatial
+kNN, ``np.matmul`` for the shifted covariance, ``np.linalg.eigh`` for the matrix
+square root — so this module imports nothing from the deep-learning stack and
+``compute_covet`` is usable without it. The ENVI CVAE components live in
+``scenvi/_nn.py``.
+"""
+
 import warnings
 
-import jax
-import jax.numpy as jnp
 import numpy as np
 import pandas as pd
 import scanpy as sc
+import scipy.sparse
 import sklearn.neighbors
-from clu import metrics
-from flax import linen as nn
-from flax import struct
-from flax.training import train_state
-from jax import random
-
 from sklearn.preprocessing import OneHotEncoder
-from tqdm import tqdm   
-import scipy.sparse
-
-class FeedForward(nn.Module):
-    """
-    :meta private:
-    """
-
-    n_layers: int
-    n_neurons: int
-    n_output: int
-
-    @nn.compact
-    def __call__(self, x):
-        """
-        :meta private:
-        """
-
-        n_layers = self.n_layers
-        n_neurons = self.n_neurons
-        n_output = self.n_output
-
-        x = nn.Dense(
-            features=n_neurons,
-            dtype=jnp.float32,
-            kernel_init=nn.initializers.glorot_uniform(),
-            bias_init=nn.initializers.zeros_init(),
-        )(x)
-        x = nn.leaky_relu(x)
-        x = nn.LayerNorm(dtype=jnp.float32)(x)
-
-        for _ in range(n_layers - 1):
-
-            x = nn.Dense(
-                features=n_neurons,
-                dtype=jnp.float32,
-                kernel_init=nn.initializers.glorot_uniform(),
-                bias_init=nn.initializers.zeros_init(),
-            )(x)
-            x = nn.leaky_relu(x) + x
-            x = nn.LayerNorm(dtype=jnp.float32)(x)
-
-        output = nn.Dense(
-            features=n_output,
-            dtype=jnp.float32,
-            kernel_init=nn.initializers.glorot_uniform(),
-            bias_init=nn.initializers.zeros_init(),
-        )(x)
-
-        return output
-
-
-class CVAE(nn.Module):
-    """
-    :meta private:
-    """
-
-    n_layers: int
-    n_neurons: int
-    n_latent: int
-    n_output_exp: int
-    n_output_cov: int
-
-    def setup(self):
-        """
-        :meta private:
-        """
-
-        n_layers = self.n_layers
-        n_neurons = self.n_neurons
-        n_latent = self.n_latent
-        n_output_exp = self.n_output_exp
-        n_output_cov = self.n_output_cov
-
-        self.encoder = FeedForward(
-            n_layers=n_layers, n_neurons=n_neurons, n_output=n_latent * 2
-        )
-
-        self.decoder_exp = FeedForward(
-            n_layers=n_layers, n_neurons=n_neurons, n_output=n_output_exp
-        )
-
-        self.decoder_cov = FeedForward(
-            n_layers=n_layers, n_neurons=n_neurons, n_output=n_output_cov
-        )
-
-    def __call__(self, x, mode="spatial", key=random.key(0)):
-        """
-        :meta private:
-        """
-
-        conf_const = 0 if mode == "spatial" else 1
-        conf_neurons = jax.nn.one_hot(
-            conf_const * jnp.ones(x.shape[0], dtype=jnp.int8), 2, dtype=jnp.float32
-        )
-
-        x_conf = jnp.concatenate([x, conf_neurons], axis=-1)
-        enc_mu, enc_logstd = jnp.split(self.encoder(x_conf), 2, axis=-1)
-
-        key, subkey = random.split(key)
-        z = enc_mu + random.normal(key=subkey, shape=enc_logstd.shape) * jnp.exp(
-            enc_logstd
-        )
-        z_conf = jnp.concatenate([z, conf_neurons], axis=-1)
-
-        dec_exp = self.decoder_exp(z_conf)
-
-        if mode == "spatial":
-            dec_cov = self.decoder_cov(z)
-            return (enc_mu, enc_logstd, dec_exp, dec_cov)
-        return (enc_mu, enc_logstd, dec_exp)
-
-
-@struct.dataclass
-class Metrics(metrics.Collection):
-    """
-    :meta private:
-    """
-
-    enc_loss: metrics.Average
-    dec_loss: metrics.Average
-    enc_corr: metrics.Average
-
-
-class TrainState(train_state.TrainState):
-    """
-    :meta private:
-    """
+from tqdm import tqdm
 
-    metrics: Metrics
 
 def batch_matrix_sqrt(Mats):
     """
diff --git a/tests/test_covet.py b/tests/test_covet.py
new file mode 100644
index 0000000..ae829e3
--- /dev/null
+++ b/tests/test_covet.py
@@ -0,0 +1,86 @@
+"""COVET does not depend on ENVI's deep-learning stack.
+
+COVET is pure numpy/sklearn/scanpy. Keeping it that way is what lets
+``compute_covet`` survive a breakage in jax/flax/optax/clu/tensorflow_probability
+instead of being taken down with it, so the separation is asserted here rather
+than left to convention.
+"""
+
+import subprocess
+import sys
+import textwrap
+
+import anndata
+import numpy as np
+import pytest
+
+import scenvi
+
+#: Modules the COVET path must never pull in.
+DEEP_LEARNING_STACK = ("jax", "flax", "optax", "clu", "tensorflow_probability")
+
+
+@pytest.fixture
+def spatial_data():
+    rng = np.random.default_rng(0)
+    return anndata.AnnData(
+        X=rng.uniform(low=0, high=100, size=(64, 8)),
+        obsm={"spatial": rng.normal(size=(64, 2))},
+    )
+
+
+def test_compute_covet_runs(spatial_data):
+    covet, covet_sqrt, cov_genes = scenvi.compute_covet(spatial_data, k=6, g=8, batch_key=-1)
+
+    n_cells, n_genes = spatial_data.shape
+    assert covet.shape == (n_cells, n_genes, n_genes)
+    assert covet_sqrt.shape == (n_cells, n_genes, n_genes)
+    assert len(cov_genes) == n_genes
+    assert np.isfinite(covet).all()
+    assert np.isfinite(covet_sqrt).all()
+
+
+def test_covet_matrices_are_symmetric_psd(spatial_data):
+    covet, covet_sqrt, _ = scenvi.compute_covet(spatial_data, k=6, g=8, batch_key=-1)
+
+    np.testing.assert_allclose(covet, covet.transpose(0, 2, 1), rtol=1e-6, atol=1e-8)
+    np.testing.assert_allclose(covet_sqrt, covet_sqrt.transpose(0, 2, 1), rtol=1e-6, atol=1e-8)
+    # The square root is the defining property, and it is what ENVI's OT loss consumes.
+    np.testing.assert_allclose(covet_sqrt @ covet_sqrt, covet, rtol=1e-5, atol=1e-6)
+    assert np.linalg.eigvalsh(covet).min() > -1e-8
+
+
+def test_covet_path_does_not_import_the_deep_learning_stack():
+    """Running COVET end to end must leave jax and friends unimported.
+
+    In a subprocess, because the check is on ``sys.modules`` and this session has
+    almost certainly imported ENVI already.
+    """
+    script = textwrap.dedent(
+        f"""
+        import sys
+
+        import anndata
+        import numpy as np
+
+        import scenvi
+
+        rng = np.random.default_rng(0)
+        spatial_data = anndata.AnnData(
+            X=rng.uniform(low=0, high=100, size=(64, 8)),
+            obsm={{"spatial": rng.normal(size=(64, 2))}},
+        )
+        scenvi.compute_covet(spatial_data, k=6, g=8, batch_key=-1)
+
+        imported = [name for name in {DEEP_LEARNING_STACK!r} if name in sys.modules]
+        assert not imported, "the COVET path imported " + ", ".join(imported)
+        """
+    )
+
+    subprocess.run([sys.executable, "-c", script], check=True)
+
+
+def test_envi_attribute_is_stable_across_repeated_access():
+    """``scenvi.ENVI`` must be the class every time, never the module of the same name."""
+    assert isinstance(scenvi.ENVI, type)
+    assert scenvi.ENVI is scenvi.ENVI