diff --git a/CHANGELOG.md b/CHANGELOG.md
index 83d69f5ef..a8d63b376 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -47,12 +47,14 @@ breaking release).
 
 - **The aliases `GradientOptimizer`, `MomentumOptimizer` and `AdamOptimizer` are removed**; use `GradientMethod`, `MomentumMethod` and `Adam`. These are no longer exported: `OptimizerMethod`, `GradientState`, `MomentumState`, `AdamState`, `optimization_step!`, `optimize_for_one_epoch!` and `Optimizer`. `TrainingOptimizer` (re-exported from GeometricOptimizers), `train!` and `AdamWithEuclideanDecay` are newly exported, and the package re-exports `GradientMethod`, `MomentumMethod`, `Adam`, `AdamOptimizerWithDecay`, `DecayingStatic`, `Cayley` and `Geodesic` from GeometricOptimizers. Methods carry no element type in GeometricOptimizers 0.9: `Adam(T)` becomes `Adam()`, and `MomentumMethod(α)` becomes `MomentumMethod(; α)`. `check(::Optimizer)` is removed.
 
-- **The three 3-tensor functors of `Dense` and `Linear` and the HDF5 `save` and `load` methods on `NeuralNetwork` are removed**, because AbstractNeuralNetworks 0.9 provides them. The HDF5 extension and the HDF5 weak dependency are gone; with HDF5 loaded, AbstractNeuralNetworks' extension provides `save` and `load`. Its `load` has no `backend` keyword: load on the CPU, then call `changebackend(backend, nn)`. Aqua's `piracies` check runs again with 0 piracies (was 3), and the ambiguity count is 0 (was 1). This closes issue #145 and known issue B7.
+- **The three 3-tensor functors of `Dense` and `Linear` and the HDF5 `save` and `load` methods on `NeuralNetwork` are removed**, because AbstractNeuralNetworks 0.9 provides them. The HDF5 extension and the HDF5 weak dependency are gone; with HDF5 loaded, AbstractNeuralNetworks' extension provides `save` and `load`. Its `load` has no `backend` keyword: load on the CPU, then call `changebackend(backend, nn)`. Aqua's `piracies` check runs again, with `Batch` given as this package's own, which keeps the two deliberate methods that take it — `train!` on a network and the `NetworkLoss` fallback — out of its report; the ambiguity count is 0 (was 1). This closes issue #145 and known issue B7.
 
 - **`Batch(::Int, ::Nothing, ::Int)` is removed.** It existed only to raise `"Cannot provide prediction window alone. Need sequence length!"`, and `Batch(batch_size, sequence_length, prediction_window)` is unchanged, so a caller of that form gets a `MethodError` where it got the error.
 
 ### Changed (breaking)
 
+- **The data loader is `ReducedComplexityModeling`'s, and `DataLoader`, `Batch` and `train!` are re-exported from it.** `parentmodule` of all three is that package, this package's method of `train!` extends its generic, and `using GeometricMachineLearning` alone still gives the three names. The loader's routines and their reference documentation move with the code: `GeometricMachineLearning.number_of_batches`, `.convert_input_and_batch_indices_to_array`, `.accuracy`, `.assign_output_estimate`, `.data_tensors_from_geometric_solution`, `.map_to_type` and `.map_to_new_backend` are `ReducedComplexityModeling`'s, and the docstrings of `DataLoader`, `Batch`, `number_of_batches`, `convert_input_and_batch_indices_to_array`, `accuracy` and `assign_output_estimate` are in that package's manual, and so are those of `data_tensors_from_geometric_solution`, `map_to_type` and `map_to_new_backend`, which ReducedComplexityModeling 0.4.1 adds. That release also declares these seven routines `public`, so this package's `[compat]` entry for it is `"0.4.1"`. `accuracy`, `number_of_batches`, `map_to_type` and `map_to_new_backend` leave this package's namespace with them; `assign_output_estimate`, `convert_input_and_batch_indices_to_array` and `data_tensors_from_geometric_solution` keep a binding here, because this package's own code calls them. `ParametricDataLoader` is public there and not here, so `using GeometricMachineLearning: ParametricDataLoader` becomes `using ReducedComplexityModeling: ParametricDataLoader`. A package that names one of these routines declares `ReducedComplexityModeling` in its own `[deps]`, as this one now does; `GeometricSolutions` leaves `[deps]` and `[compat]` with the code that read its types.
+
 - **`DataLoader`'s `autoencoder` defaults to `false` in every constructor that takes it, except the copying constructor, where `autoencoder = nothing` inherits the loader's kind, and a loader given an input and an output is a set of samples.** A matrix or a vector is a time series unless `autoencoder = true` is passed: that keyword is what reads a matrix's second axis as the parameter index, and the matrix and the vector method defaulted to `true`, so migrate `DataLoader(matrix)` to `DataLoader(matrix; autoencoder = true)` where the columns are samples — `scripts/reproduction/symplectic_autoencoders/training.jl` and the tutorials do. `DataLoader(input, output)` builds a `:RegularData` loader whatever its kind, so its target is its output and never the input shifted by the sequence length: in 0.8.0 a `DataLoader(x, y)` of ten pairs drew the samples 1 to 9 and, where the input and the output differed in array type, returned one minibatch holding the whole data set. The count is `cld(n · n_params, batch_size)` with `n` the start positions of one parameter index, and a loader and a `Batch` that leave no start position raise an `ArgumentError` naming the time steps, the sequence length and the prediction window, where 0.8.0 asserted and an epoch could be empty. `Batch(batch_size, dl)` takes its two lengths from the loader, so a classifier's minibatch holds the whole patch sequence and the whole one-hot output of each of its samples. On a loader whose output's time axis differs from its input's, such as a classifier's with one output column, `Batch(batch_size)` raises an `ArgumentError` that names `Batch(batch_size, dl)`, where 0.8.0's minibatch held the whole input sequence; migrate it to `Batch(batch_size, dl)`. That batch is a `Batch{:Transformer}`, so `train!` on a `NeuralNetworkIntegrator` with it needs an explicit `loss`. `suppress_info` is a keyword of every constructor that takes a data set, which makes `DataLoader(solution, true)` a `MethodError`.
 
 - **`QPT` and `QPTOAT` are gone; the `(q, p)` methods dispatch on AbstractNeuralNetworks' `NamedTupleOfArrays` and `ArrayOrNamedTuple`.** `GeometricMachineLearning.QPT{T}` becomes `NamedTupleOfArrays{T, (:q, :p)}` and `GeometricMachineLearning.QPTOAT{T}` becomes `ArrayOrNamedTuple{T}`. Both names are public in AbstractNeuralNetworks from 0.9.2, so its `[compat]` floor is `"0.9.2"`. `NamedTupleOfArrays` ties the element types of `q` and `p` together but not their array types, so `Dense`, `ResNetLayer`, `MultiHeadAttention`, `PositionalEncoding`, `TransformerLoss`, `AutoEncoderLoss`, `ReducedLoss`, `ZygotePullback` and `DataLoader(dl, T)` take `(q, p)` data whose `q` is a view and whose `p` is an `Array`, which `QPT` did not admit. The signatures that took `QPTOAT` — `AutoEncoderLoss`, `ReducedLoss`, `HNNLoss`, `ZygotePullback`, `SymbolicPullback(arch)` and `DataLoader(dl, T)` — admit a named tuple of arrays with any keys, as AbstractNeuralNetworks' own methods do. Migrate `using GeometricMachineLearning: QPT, QPTOAT` to `using AbstractNeuralNetworks: NamedTupleOfArrays, ArrayOrNamedTuple`, `QPT` to `NamedTupleOfArrays{<:Any, (:q, :p)}` and `QPTOAT` to `ArrayOrNamedTuple`.
@@ -69,11 +71,11 @@ breaking release).
 
 - **`ResNetLayer` encodes its bias as a `Symbol` type parameter, the two volume-preserving feedforward layers are one struct, and two fields and a keyword are renamed.** `ResNetLayer`'s third type parameter is `:bias` or `:no_bias`, as in the volume-preserving feedforward layers, where it was `true` or `false`; its keyword `use_bias::Bool` stays. `VolumePreservingFeedForwardLayer{M, N, bias, AT, Triangle}` is now the one struct of the two volume-preserving feedforward layers, with `Triangle` `:lower` or `:upper`, where it was their abstract supertype; `VolumePreservingLowerLayer{M, N, bias, AT}` and `VolumePreservingUpperLayer{M, N, bias, AT}` are its aliases and keep their constructors. Their `initialparameters` takes the bias initializer as the keyword `init_bias`, as `ResNetLayer` and the gradient layers do, where it took `init_bias!`. `GSympNet` stores its activation in the field `activation`, as `LASympNet` does, where it was `act`. `StandardTransformerIntegrator`'s field `upsacling_activation` is spelled `upscaling_activation`. `GradientLayer` stores `second_dim` as an `Int`, and the six SympNet layer constructors (`GradientLayerQ`/`P`, `ActivationLayerQ`/`P`, `LinearLayerQ`/`P`) convert the dimension to an `Int`, so `GradientLayerQ(Int32(4), Int32(8), tanh) === GradientLayerQ(4, 8, tanh)` and `ActivationLayerQ(Int32(4), tanh) === ActivationLayerQ(4, tanh)`.
 
-- **The `[compat]` bounds move to the versions the training API needs**: GeometricOptimizers `"0.9"`, AbstractNeuralNetworks `"0.9.2"`, NeuralNetworkParameters `"0.5"`, SymbolicNeuralNetworks `"0.9"`, GeometricBase `"0.15"`, GeometricEquations `"0.21.5"` and GeometricSolutions `"0.6.6"`. HDF5 leaves `[compat]` and `[weakdeps]`.
+- **The `[compat]` bounds move to the versions the training API needs**: GeometricOptimizers `"0.9"`, AbstractNeuralNetworks `"0.9.2"`, NeuralNetworkParameters `"0.5"`, GeometricBase `"0.15"` and GeometricEquations `"0.21.5"`. HDF5 leaves `[compat]` and `[weakdeps]`.
 
 ### Fixed
 
-- **`accuracy` returns `float` of the element type of the loader's input.** `size(dl.output, 3)` is an `Int` and the loader's label type an `Integer`, so the result was a `Float64` for every loader: a `Float32` loader of six correctly labelled samples gave `1.0` where it gives `one(Float32)`, and a `Float64` model over a `Float32` loader's data gave a `Float64` measure where it now reads `Float32`. The difference is counted in an integer and halved in floating point, so a `Bool` one-hot target gives the same value as an `Int64` one, and the fraction is computed in `F = float(T)`, where `T` is the element type of the loader's input — `n = size(dl.output, 3)`, `errors = sum(abs.(dl.output - tensor_of_maximum_elements)) / 2`, `F = float(T)`, `F(n - errors) / F(n)` — so the value of the fraction does not move, an integer input gives a `Float64`, and a sample that carries no label still counts as half an error, because its one-hot column differs from the prediction in one entry where a wrong label differs in two. The test is `test/data_loader/accuracy.jl`.
+- **`accuracy` returns `float` of the element type of the loader's input.** `size(dl.output, 3)` is an `Int` and the loader's label type an `Integer`, so the result was a `Float64` for every loader: a `Float32` loader of six correctly labelled samples gave `1.0` where it gives `one(Float32)`, and a `Float64` model over a `Float32` loader's data gave a `Float64` measure where it now reads `Float32`. The difference is counted in an integer and halved in floating point, so a `Bool` one-hot target gives the same value as an `Int64` one, and the fraction is computed in `F = float(T)`, where `T` is the element type of the loader's input — `n = size(dl.output, 3)`, `errors = sum(abs.(dl.output - tensor_of_maximum_elements)) / 2`, `F = float(T)`, `F(n - errors) / F(n)` — so the value of the fraction does not move, an integer input gives a `Float64`, and a sample that carries no label still counts as half an error, because its one-hot column differs from the prediction in one entry where a wrong label differs in two. The test is `ReducedComplexityModeling`'s `test/data_loader/accuracy.jl`.
 
 - **`map_to_cpu` returns a network whose backend is the host's.** The copy rebuilt every parameter leaf as a host `Array` but kept the network's backend and its fourth type parameter, so a network that came from a device held host data behind a device backend and its next forward pass sent that data to the device it came from. It is now `NeuralNetwork{AT, MT, typeof(ps), typeof(backend)}(nn.architecture, nn.model, ps, backend)` with `backend = nn.backend isa CPU ? nn.backend : CPU()`: a device network comes back on `CPU()`, a network already on a host backend keeps the one it was handed, so `CPU(; static = true)` comes back static, and the fourth type parameter is the backend's type. `test/map_to_cpu.jl` pins both halves on a `JLArray` — a `JLBackend` network, and `CPU()` and `CPU(; static = true)`.
 
@@ -117,7 +119,7 @@ breaking release).
 
 - **`ForcingLayer`, `ForcingLayerQ`, `ForcingLayerP`, `ForcingLayerQP` and `ForcedSympNet` add a learned non-conservative term to a symplectic integrator.** The five names are `public` and not exported. A forcing layer changes `p` and leaves `q` alone, by adding the output of a `Chain` of its own to `p`; the `Q`, `P` or `QP` of its type says what that network reads — `q`, `p`, or both — and `parameters` fixes the length of the flat vector of system parameters it takes beside the state. `return_parameters = true` makes it return `(output, μ)` instead of the output, so that the layer after it in a `Chain` takes the system parameters too. `dim` is the full phase-space dimension, as for `GradientLayerQ`, and the forcing has `dim ÷ 2` rows; an odd dimension and a forcing type other than `:Q`, `:P` and `:QP` are `ArgumentError`s at construction. `ForcedSympNet(d; upscaling_dimension = 2d, n_layers = 2, n_forcing_layers = 2, activation = tanh, init_upper = true, forcing_type = :P)` is the chain of a `GSympNet` with a forcing layer in front of every layer that changes `p`: it holds `n_layers` gradient layers and one forcing layer per step that changes `p`, it takes no parameters of the system, and it answers `dim` as the `SympNet`s do. A forcing layer holds a `Chain`, so its parameters are a set within the set of the network, which `train!`, `ZygotePullback` and `optimization_step!` carry as they carry any nested set. `scripts/reproduction/damped_oscillator_forced_sympnet.jl` trains one on the damped harmonic oscillator through `TrainingOptimizer` and `train!`. The code is ported from pull request #207 by Benedikt Brantner (relates to #207); there a forcing layer raised "This shouldn't be the case!" for any other tuple, which is a `MethodError` now, the forcing type of the architecture was accepted unchecked and failed at the first call, and the script trained through `Optimizer(AdamOptimizer(), nn)`.
 
-- **`apply_parametric`, `ParametricLoss` and `ParametricDataLoader` train a network whose forward pass takes the parameters of the system beside its input.** The three names are `public` and not exported, so a caller qualifies them or names them in a `using GeometricMachineLearning: …`. `apply_parametric(model, x, μ, ps)` applies a `Chain` to one sample — a vector or a `(q, p)` pair of vectors — with one `μ`, or to a batch — a `2d × 1 × k` tensor or a `(q, p)` pair of `d × 1 × k` tensors — with a vector of `k` values of `μ`, sample by sample; the output has the shape of the input. The `Chain` receives the pair `(x, μ)`, and the internal `concatenate_array_with_parameters` puts the flattened `μ` below the state in the element type of the state, side by side for a matrix and a vector of `μ`, repeated in every column for a matrix and one `μ`, and repeated in every column of every slice for a 3-tensor and one `μ`. `ParametricLoss()(model, ps, input, output, μ)` is the relative norm of the error of `apply_parametric`, and `ZygotePullback` takes the triple `(input, output, μ)` with one `μ` or a vector of them. `ParametricDataLoader(data, parameters)` takes a `2d × nt × m` array or a `(q, p)` pair of `d × nt × m` arrays of one array type with a vector of `m` named tuples, and `ParametricDataLoader(ensemble_solution)` takes an `EnsembleSolution` of `(q, p)` states with the parameters of its problem. A `Batch(b)` draws each pair of a start position and a trajectory once per epoch, as for `DataLoader`, and a minibatch is the state, the state one time step later and the `μ` of each sample. `train!` takes a `ParametricDataLoader` and trains on `ParametricLoss()` when no `loss` is given. A `Batch` with a sequence length or a prediction window other than one, a parameter vector whose length is not the number of trajectories, no trajectory, too few time steps for the batch, a `μ` whose element type is not the data's, and a `μ` whose keys, their order or the sizes of its leaves differ from the first `μ`'s raise an `ArgumentError`; a `(q, p)` pair of two array types is a `MethodError`, as for `DataLoader` (K20). The code is ported from pull request #207 by Benedikt Brantner (relates to #207); there `concatenate_array_with_parameters` promoted a `Float32` state to `Float64`, rejected a matrix of more than one column with one `μ`, and the loader and the batch checked their counts with `@assert`.
+- **`apply_parametric`, `ParametricLoss` and `ParametricDataLoader` train a network whose forward pass takes the parameters of the system beside its input.** `apply_parametric` and `ParametricLoss` are `public` and not exported, so a caller qualifies them or names them in a `using GeometricMachineLearning: …`; `ParametricDataLoader` is `ReducedComplexityModeling`'s and is public there, not here, so a caller takes it from that package or writes `GeometricMachineLearning.ParametricDataLoader`. `apply_parametric(model, x, μ, ps)` applies a `Chain` to one sample — a vector or a `(q, p)` pair of vectors — with one `μ`, or to a batch — a `2d × 1 × k` tensor or a `(q, p)` pair of `d × 1 × k` tensors — with a vector of `k` values of `μ`, sample by sample; the output has the shape of the input. The `Chain` receives the pair `(x, μ)`, and the internal `concatenate_array_with_parameters` puts the flattened `μ` below the state in the element type of the state, side by side for a matrix and a vector of `μ`, repeated in every column for a matrix and one `μ`, and repeated in every column of every slice for a 3-tensor and one `μ`. `ParametricLoss()(model, ps, input, output, μ)` is the relative norm of the error of `apply_parametric`, and `ZygotePullback` takes the triple `(input, output, μ)` with one `μ` or a vector of them. `ParametricDataLoader(data, parameters)` takes a `2d × nt × m` array or a `(q, p)` pair of `d × nt × m` arrays of one array type with a vector of `m` named tuples, and `ParametricDataLoader(ensemble_solution)` takes an `EnsembleSolution` of `(q, p)` states with the parameters of its problem. A `Batch(b)` draws each pair of a start position and a trajectory once per epoch, as for `DataLoader`, and a minibatch is the state, the state one time step later and the `μ` of each sample. `train!` takes a `ParametricDataLoader` and trains on `ParametricLoss()` when no `loss` is given. A `Batch` with a sequence length or a prediction window other than one, a parameter vector whose length is not the number of trajectories, no trajectory, too few time steps for the batch, a `μ` whose element type is not the data's, and a `μ` whose keys, their order or the sizes of its leaves differ from the first `μ`'s raise an `ArgumentError`; a `(q, p)` pair of two array types is a `MethodError`, as for `DataLoader` (K20). The code is ported from pull request #207 by Benedikt Brantner (relates to #207); there `concatenate_array_with_parameters` promoted a `Float32` state to `Float64`, rejected a matrix of more than one column with one `μ`, and the loader and the batch checked their counts with `@assert`.
 
 - **`WideResNetLayer`, `ParametricResNetLayer` and `ParametricResNet` add a ResNet whose hidden layer has a width of its own, and a baseline whose blocks take the parameters of the system.** The three names are `public` and not exported, so a caller qualifies them or names them in a `using GeometricMachineLearning: …`. `WideResNetLayer(dim, width, activation)` computes `x + σ.(W_down σ.(W_up x .+ b_up) .+ b)` for a vector, a matrix, a tensor or a `(q, p)` pair, and `ResNet(dim, n_blocks, width, activation)` and `ResNet(dl, n_blocks, width)` build it wherever `width` differs from the system dimension; a width equal to the dimension keeps the plain `ResNetLayer`, so `ResNet(dim, n_blocks, activation)` and its recorded outputs are unchanged. `ParametricResNet(dim; width = dim, n_blocks = 1, activation = tanh, parameters)` appends the flattened `μ` below the state at the input of every upscaling weight, through `ParametricResNetLayer`s that pass `μ` on to the next layer, and the `NeuralNetwork` built around it gains the functor `nn(x, μ)`, which forwards to `apply_parametric` and so takes one sample with one `μ` or a batch with a vector of `μ`; an input shape that `apply_parametric` does not take raises its `MethodError`. `parameterlength` of both layers counts what `initialparameters` allocates: `width·(dim + 1) + dim·(width + 1)` without system parameters, and `width·(dim + |μ|) + width + dim·(width + 1)` with them. `scripts/reproduction/parametric_resnet.jl` trains the baseline on the analytic data of the harmonic oscillator, with the time as the parameter of the system, through `ParametricDataLoader`, `ParametricLoss` and `train!`. This is the port of pull request #207 by Benedikt Brantner (relates to #207), with `parameterlength(::ParametricResNetLayer)` corrected -- #207's `(width + |μ|)(dim + 1) + dim(width + 1)` counts the entries only at `width = dim + 1` -- and the two defects that #207's own test turned up kept fixed: `ParametricResNet(::DataLoader, …)` forwards `parameters` and `activation` instead of dropping them, and `nn(x, μ)` is the functor of the network and not AbstractNeuralNetworks' generic one.
 
@@ -127,7 +129,7 @@ breaking release).
 
 - **The test suite runs on an Apple GPU in its `metal` group.** `test/devices/metal.jl` checks the device paths that were measured by hand and not tested before: in the testset "_mul and tensor_tensor_mul on Metal match the host", `_mul` with a dense and a `SymmetricMatrix` weight from the left, a dense one from the right, and `tensor_tensor_mul` on `MetalBackend()`; and in "tensor_cayley and VolumePreservingAttention on Metal", the Cayley kernel and the layer that uses it, at every sequence length from 2 to 8. Each result must stay on the device and match the host in `Float32`, with scalar indexing off. `Pkg.test(test_args = ["metal"])` runs them without the subject tests, and a missing device records a skip, which happens inside a sandbox. Metal (1.10 or later) is a test dependency on every platform; it installs and precompiles on Linux and Windows, and nothing there loads it. On an Apple M4 Max with Metal 1.11.1 and Julia 1.13.1, the `_mul` and `tensor_tensor_mul` testset passes in a cold process, 8 of 8 assertions.
 - **A `Metal` workflow runs the Metal tests on GitHub's `macos-15` runner**, for the `min` and `1` Julia versions. The runner's GPU is Apple's paravirtual device, which Metal.jl supports from 1.10 and on macOS 15 or later. The workflow includes a step "Metal is functional" before the tests that fails the job where the device is unavailable, so it cannot pass without having run the tests. The job is not a required check. However, the Metal tests also run in two required checks—CI.yml's `Julia min - macOS-latest - default` and `Julia 1 - macOS-latest - default`—which execute `test/runtests.jl` with no test arguments on Apple silicon, selecting the `metal` group.
-- **`test/quality/jet.jl` checks the kernel launchers with JET.** It runs in `core`, directly after `quality/aqua.jl`, and takes about 6 s. Each function that launches a kernel of `src/` on the `CPU()` backend, the pullbacks among them, and the `Batch` functor that `test/data_loader/batch_index_set.jl` measures with `@allocated` get one `JET.report_opt` line per element type that the tests use, with `target_modules = (GeometricMachineLearning,)`. Each `@kernel` method also gets a line that analyses its generated `cpu_<kernel>` function at the launcher's context, because JET drops the reports of a kernel body when it analyses the launcher. Its 50 assertions all pass on Julia 1.13.2 (JET 0.12.2), and none is `@test_broken`. Where JET does not work, the file records one `@test_skip`. JET joins `test/Project.toml` with no `[compat]` bound. `KNOWN_ISSUES.md` K2 records the Revise `EMFILE` messages that JET brings into the test log, and K3 the dispatch that the lines do not see: an unstable argument of a kernel launch.
+- **`test/quality/jet.jl` checks the kernel launchers with JET.** It runs in `core`, directly after `quality/aqua.jl`, and takes about 6 s. Each function that launches a kernel of `src/` on the `CPU()` backend, the pullbacks among them, gets one `JET.report_opt` line per element type that the tests use, with `target_modules = (GeometricMachineLearning,)`. Each `@kernel` method also gets a line that analyses its generated `cpu_<kernel>` function at the launcher's context, because JET drops the reports of a kernel body when it analyses the launcher. Its 34 assertions all pass on Julia 1.13.3 (JET 0.13.1), and none is `@test_broken`. Where JET does not work, the file records one `@test_skip`. JET joins `test/Project.toml` with no `[compat]` bound. `KNOWN_ISSUES.md` K2 records the Revise `EMFILE` messages that JET brings into the test log, and K3 the dispatch that the lines do not see: an unstable argument of a kernel launch.
 - **A test calls `integrate` on `(q,)` initial conditions.** No test reached the `NeuralNetworkIntegrator` `integrate` method that takes `ics::NamedTuple{(:q,), Tuple{AT}}`. `test/integration/docstrings/layers_and_architectures.jl` runs the docstring's `ResNet` from `(q = [1, 1, 1],)` for four points and expects the same trajectory as the vector form, `[1 2 4 8; 1 3 9 27; 1 3 7 15]`.
 
 - **`NetworkLoss(arch)` is defined for a `NeuralNetworkIntegrator` `arch` and returns `FeedForwardLoss()`.** It is the default loss of that architecture: `train!` takes it for a `Batch{:FeedForward}` when no `loss` is given. A caller can name the default, or pass it explicitly, without relying on `train!`'s keyword default.
diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md
index f30fa8aec..5951ec65a 100644
--- a/KNOWN_ISSUES.md
+++ b/KNOWN_ISSUES.md
@@ -105,13 +105,13 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
 
 ### K20 · The `DataLoader((q = …, p = …))` constructors require one array type for `q` and `p`.
 
-- **location:** `src/data_loader/data_loader.jl:216`
+- **location:** `ReducedComplexityModeling` `src/data_loader/data_loader.jl:220`
 - **kind:** defect
 - **found:** #353 (2026-10-23); the constructors are older than that pull request
 - **evidence:**
-  The three methods at `src/data_loader/data_loader.jl:216-245` take
-  `NamedTuple{(:q, :p), Tuple{AT, AT}}` (and `Tuple{VT, VT}`), so `q` and `p` must have one
-  array type, while the `TransformerLoss` method at `src/loss/losses.jl:92` and the
+  The three methods at `ReducedComplexityModeling` `src/data_loader/data_loader.jl:220-249`
+  take `NamedTuple{(:q, :p), Tuple{AT, AT}}` (and `Tuple{VT, VT}`), so `q` and `p` must have one
+  array type, while the `TransformerLoss` method at `src/loss/losses.jl:98` and the
   `MultiHeadAttention` call take `NamedTupleOfArrays`, which admits two. A view for `q` and a
   `Matrix` for `p` raises:
 
@@ -122,10 +122,47 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
   ```
 
   The same call with `q = A` gives a `DataLoader` with `input_dim == 4`. The docstring at
-  `src/data_loader/data_loader.jl:205` states the one-type signature. Closing it means widening the
-  three signatures to `NamedTupleOfArrays{T, (:q, :p)}`, and the docstring with them.
-  `ParametricDataLoader`'s `(q, p)` constructor has the same one-type restriction and widens with
-  it.
+  `ReducedComplexityModeling` `src/data_loader/data_loader.jl:208` states the one-type signature.
+  Closing it means widening the three signatures to `NamedTupleOfArrays{T, (:q, :p)}`, and the
+  docstring with them. `ParametricDataLoader`'s `(q, p)` constructor has the same one-type
+  restriction and widens with it.
+
+### K37 · `metal_abort_runs.jl audit` throws when a run's crash report is gone, although its log records the signature.
+
+- **location:** `scripts/utilities/metal_abort_runs.jl:293`
+- **kind:** defect
+- **found:** #365
+- **evidence:**
+  `read_run` calls `report_signature(report_path(path, index, file, reports))` for every `M8-IPS` line.
+  `report_path` returns the copy beside the log when it exists and else `joinpath(reports, file)`, and
+  `report_signature` is `read(path, String)` on it (`:101-102`), so a copy that is missing with the
+  original deleted from `~/Library/Logs/DiagnosticReports` raises a `SystemError` and the audit of the
+  arm stops. Each `M8-IPS` line already holds `signature=` (`:426-427`), which the audit could use when
+  the file is absent. Read from the code; the failure is not reproduced.
+
+### K38 · The session arm's reader stamps `R0` after `yield()`, so `R0` is not the instant its first read finished.
+
+- **location:** `scripts/utilities/metal_abort_session.jl:134-137`
+- **kind:** defect
+- **found:** #365
+- **evidence:**
+  The docstring of `reader` (`:117-123`) says that `R0` is the instant the first read finished, and the
+  loop assigns it after `sum(v)`, `yield()` and `reads += 1`. When `yield()` hands the thread to another
+  task for a long stretch, `R0` lies later than the read, and `T0 < R0 < T1` can fail for a run whose
+  first read overlapped the training call, so the run counts as void and not as concurrent. The error
+  is on the side of a void run. Read from the code; not measured.
+
+### K39 · The session arm shares `M8_T0`, `M8_T1` and `M8_TRAINER` between two threads with no synchronization.
+
+- **location:** `scripts/utilities/metal_abort_session.jl:102-106`
+- **kind:** not verified
+- **found:** #365
+- **evidence:**
+  `training()` assigns the three globals and `reader()` polls and reads them (`:127`, `:133`,
+  `:141-142`) on another thread, with no lock and no atomic. Julia gives no guarantee for such reads of
+  a global that another thread writes. CodeRabbit marks the finding as one it could not verify, and no
+  wrong value has been observed: each global is a single `Int`, assigned around the training
+  call.
 
 ## Upstream
 
@@ -276,7 +313,7 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
 
 ### K33 · An epoch of `train!` on a Metal array advances the CPU random stream, so a two-epoch run draws its second epoch's minibatches differently from the host's.
 
-- **location:** `src/data_loader/batch.jl:195`
+- **location:** `ReducedComplexityModeling` `src/data_loader/batch.jl:186`
 - **kind:** not verified
 - **found:** 2026-10-26, part M7
 - **evidence:**
@@ -385,17 +422,38 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
 - **found:** 2026-10-01
 - **evidence:**
   Julia 1.13.2, `mutate.jl <repository> <file> <from> <to> quality/jet.jl`, one mutant per run.
-  Three mutants that put `Base.inferencebarrier` on an argument of a kernel launch SURVIVED, each
-  with 50 pass, 0 fail: `kernel!(c, Base.inferencebarrier(a), b, ndrange = size(c))` in
-  `src/kernels/tensor_tensor_mul.jl`; `tensor_mat_skew_sym_assign_k!(C,
+  Two mutants that put `Base.inferencebarrier` on an argument of a kernel launch SURVIVED, with no
+  assertion failing: `kernel!(c, Base.inferencebarrier(a), b, ndrange = size(c))` in
+  `src/kernels/tensor_tensor_mul.jl`, and `tensor_mat_skew_sym_assign_k!(C,
   Base.inferencebarrier(Z), A, ndrange = size(C))` in
-  `src/kernels/inverses/tensor_mat_skew_sym_assign.jl`; and `gather_kernel!(backend)(gathered,
-  Base.inferencebarrier(data), indices, offset, ndrange = size(gathered))` in
-  `src/data_loader/batch.jl`. At `target_modules = (GeometricMachineLearning,)` a barrier on a
-  launch argument gives 0 reports, and the same barrier in an `@assert` of the same launcher is
-  CAUGHT with 3 reports (47 pass, 3 fail): the kernel functor takes `args...`, so the dispatch sits
-  in KernelAbstractions' frames, which the filter drops. A barrier on a value inside the kernel
-  body of `tensor_tensor_mul_kernel!` is CAUGHT (47 pass, 3 fail).
+  `src/kernels/inverses/tensor_mat_skew_sym_assign.jl`. At
+  `target_modules = (GeometricMachineLearning,)` a barrier on a launch argument gives 0 reports, and
+  the same barrier in an `@assert` of the same launcher is CAUGHT: the kernel functor takes
+  `args...`, so the dispatch sits in KernelAbstractions' frames, which the filter drops. A barrier
+  on a value inside the kernel body of `tensor_tensor_mul_kernel!` is CAUGHT.
+
+  The loader's gather launch is outside this file's reach. Its launcher,
+  `gather_kernel!(backend)(gathered, data, indices, offset, ndrange = size(gathered))`, is
+  `ReducedComplexityModeling`'s `src/data_loader/batch.jl`, which `target_modules =
+  (GeometricMachineLearning,)` filters out, and gating it belongs to that package's
+  `test/quality/jet.jl`.
+
+### K36 · Three files of `test/data_loader/` test code that is no longer in `src/data_loader/`.
+
+- **location:** `test/data_loader/batch.jl`
+- **kind:** defect
+- **found:** #365
+- **evidence:**
+  `src/data_loader/` holds `train.jl` alone, since the loader is `ReducedComplexityModeling`'s. Three
+  test files of `test/data_loader/` have no source file of that name:
+  `batch.jl` tests the `Batch{:Transformer}` method of `integrate`, which is in
+  `src/architectures/transformer_integrator.jl`; `data_loader.jl` tests `solve!` of a `PSDArch`,
+  which is in `src/architectures/psd.jl`; and `parametric_data_loader.jl` tests `train!` on a
+  `ParametricDataLoader`, which is in `src/data_loader/train.jl`. The layout rule is that
+  `test/<path>.jl` tests `src/<path>.jl`; `test-layout.jl --check` prints nothing for the branch, so
+  the check does not see the three. Moving them to `test/architectures/transformer_integrator.jl`,
+  `test/architectures/psd.jl` and `test/data_loader/train.jl` also changes the lines that list them
+  in `test/runtests.jl`.
 
 ### K14 · `train!` infers its history as an abstract `Vector`, although the run-time history is `Vector{T}`.
 
diff --git a/Project.toml b/Project.toml
index f250397d4..47c3b04ca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,13 +11,13 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeometricBase = "9a0b12b7-583b-4f04-aa1f-d8551b6addc9"
 GeometricEquations = "c85262ba-a08a-430a-b926-d29770767bf2"
 GeometricOptimizers = "fc236c15-5557-4942-aa65-b650f329279e"
-GeometricSolutions = "7843afe4-64f4-4df4-9231-049495c56661"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 NeuralNetworkParameters = "67f4d93a-60e9-472b-8cdd-1ccf6005724a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReducedComplexityModeling = "abb4722c-1d8e-48eb-a3bd-45b3116bfaba"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
@@ -28,12 +28,12 @@ ForwardDiff = "1"
 GeometricBase = "0.15"
 GeometricEquations = "0.21.5"
 GeometricOptimizers = "0.9"
-GeometricSolutions = "0.6.6"
 KernelAbstractions = "0.9"
 LinearAlgebra = "1"
 NNlib = "0.8, 0.9"
 NeuralNetworkParameters = "0.5"
 ProgressMeter = "1"
 Random = "1"
+ReducedComplexityModeling = "0.4.1"
 Zygote = "0.7"
 julia = "1.12"
diff --git a/docs/check_references.jl b/docs/check_references.jl
index f49926161..a6683d1af 100644
--- a/docs/check_references.jl
+++ b/docs/check_references.jl
@@ -43,9 +43,8 @@ function resolve(text::AbstractString, mod::Module = Main; exact_signature::Bool
     DocSystem.iskeyword(binding) && return (true, "")
     DocSystem.defined(binding) || return (false, "undefined binding `$binding`")
     # A `@docs` entry has to match a specific method, but a `@ref` resolves to whatever anchor the
-    # manual created for the binding -- `[`accuracy(::Chain, ::Tuple, ::DataLoader)`](@ref)` links to
-    # `@docs accuracy` even though no method has that exact signature. Checking a `@ref` against the
-    # signature would reject those.
+    # manual created for the binding, so a `@ref` whose signature no method has exactly still links
+    # to the binding's `@docs` block. Checking a `@ref` against the signature would reject those.
     typesig = if exact_signature
         try
             Core.eval(mod, DocSystem.signature(ex, String(text)))
diff --git a/docs/src/architectures/transformer.md b/docs/src/architectures/transformer.md
index e3e298342..c76638aa7 100644
--- a/docs/src/architectures/transformer.md
+++ b/docs/src/architectures/transformer.md
@@ -58,7 +58,6 @@ StandardTransformerIntegrator
 Transformer
 ClassificationTransformer
 ClassificationLayer
-GeometricMachineLearning.assign_output_estimate
 ```
 
 ```@raw latex
diff --git a/docs/src/data_loader/data_loader.md b/docs/src/data_loader/data_loader.md
index 340f4fb6d..a34712c1d 100644
--- a/docs/src/data_loader/data_loader.md
+++ b/docs/src/data_loader/data_loader.md
@@ -1,6 +1,8 @@
 # The Data Loader 
 
-The `DataLoader` in `GeometricMachineLearning` is designed to make training convenient. 
+The `DataLoader` that `GeometricMachineLearning` re-exports from
+[ReducedComplexityModeling.jl](https://github.com/JuliaRCM/ReducedComplexityModeling.jl) is designed
+to make training convenient. 
 
 The data loader can be called with various types of arrays as input, for example a [snapshot matrix](@ref "Snapshot Matrix"):
 
@@ -49,7 +51,7 @@ dl.input_dim
 
 ## Drawing Batches with `GeometricMachineLearning` 
 
-If we want to draw mini batches from a data set, we need to allocate an instance of [`Batch`](@ref). If we call the corresponding functor on an instance of [`DataLoader`](@ref) we get the following result[^1]:
+If we want to draw mini batches from a data set, we need to allocate an instance of `Batch`. If we call the corresponding functor on an instance of `DataLoader` we get the following result[^1]:
 
 [^1]: We first demonstrate how to sample data on the example of a matrix. The case of sampling from a tensor is slightly more complicated and is explained below.
 
@@ -146,7 +148,7 @@ Sampling from a tensor is done the following way (``\mathcal{I}_i`` again denote
 3. ``\mathcal{I}_i \leftarrow \mathtt{complete\_indices[}(i - 1) \cdot \mathtt{batch\_size} + 1 : i \cdot \mathtt{batch\_size]}\text{ for }i=1, \ldots, (\mathrm{last} -1),``
 4. ``\mathcal{I}_\mathrm{last} \leftarrow \mathtt{complete\_indices[}(\mathrm{last} - 1) \cdot \mathtt{batch\_size} + 1\mathtt{:end]}.``
 
-Note that we supplied two additional arguments to the [`Batch`](@ref) constructor here: `seq_length` and `prediction_window`. These two arguments specify how many time are considered in one mini batch and how long the prediction runs into the future respectively. These two numbers are explained when we talk about [structure on product spaces](@ref "How is Structure Preserved?").
+Note that we supplied two additional arguments to the `Batch` constructor here: `seq_length` and `prediction_window`. These two arguments specify how many time are considered in one mini batch and how long the prediction runs into the future respectively. These two numbers are explained when we talk about [structure on product spaces](@ref "How is Structure Preserved?").
 
 This algorithm can be visualized the following way (here `batch_size = 4`):
 
@@ -157,13 +159,13 @@ Here the sampling is performed over the second axis (the *time step dimension*)
 
 ## Library Functions
 
-```@docs
-DataLoader
-Batch
-Batch(::Int)
-Batch(::Int, ::Int, ::Int)
-GeometricMachineLearning.number_of_batches
-GeometricMachineLearning.convert_input_and_batch_indices_to_array
-GeometricMachineLearning.accuracy
-GeometricMachineLearning.ParametricDataLoader
-```
+`DataLoader`, `Batch` and the loader's routines — `number_of_batches`,
+`convert_input_and_batch_indices_to_array`, `accuracy`, `assign_output_estimate` and
+`ParametricDataLoader` — are `ReducedComplexityModeling`'s, and
+[that package's manual](https://juliarcm.github.io/ReducedComplexityModeling.jl/stable/) documents
+them. This package re-exports the first two, so `using GeometricMachineLearning` alone gives both.
+
+[`train!`](@ref) is that package's generic as well, re-exported here. This package defines the
+method that trains a network on a `DataLoader` or a `ParametricDataLoader`, and this manual's
+`train!` entry is its docstring; the generic's own docstring, `train!(model, …)`, is in that
+package's manual.
diff --git a/docs/src/data_loader/snapshot_matrix.md b/docs/src/data_loader/snapshot_matrix.md
index 4d28f3707..8918ecfc6 100644
--- a/docs/src/data_loader/snapshot_matrix.md
+++ b/docs/src/data_loader/snapshot_matrix.md
@@ -31,9 +31,9 @@ The snapshot tensor fulfills the same role as the snapshot matrix but has a thir
 ![The three axes of the snapshot tensor are (i) dimension, (ii) time and (iii) parameters.](https://juliagni.github.io/GeometricFigures.jl/figures/data/tensor/tensor_light.svg)
 ![The three axes of the snapshot tensor are (i) dimension, (ii) time and (iii) parameters.](https://juliagni.github.io/GeometricFigures.jl/figures/data/tensor/tensor_dark.svg)
 
-When drawing samples from the snapshot tensor to train a neural networks we also need to specify a *sequence length* (as an argument to the [`Batch`](@ref) struct). When sampling a batch from the snapshot tensor we sample over the second axis (the *time dimension*) and the third axis of the tensor (the *parameter dimension*). The total number of batches[^2] is 
+When drawing samples from the snapshot tensor to train a neural networks we also need to specify a *sequence length* (as an argument to the `Batch` struct). When sampling a batch from the snapshot tensor we sample over the second axis (the *time dimension*) and the third axis of the tensor (the *parameter dimension*). The total number of batches[^2] is 
 
-[^2]: The number of batches shown here is for the case ``\mathtt{prediction\_window} = 0.`` If ``\mathtt{prediction\_window} \neq 0,`` this number may be smaller. The corresponding function is [`GeometricMachineLearning.number_of_batches`](@ref).
+[^2]: The number of batches shown here is for the case ``\mathtt{prediction\_window} = 0.`` If ``\mathtt{prediction\_window} \neq 0,`` this number may be smaller. The corresponding function is `ReducedComplexityModeling.number_of_batches`.
 
 ```math
 \lceil\mathtt{(dl.input\_time\_steps - batch.seq\_length) * dl.n\_params / batch.batch\_size}\rceil,
diff --git a/docs/src/introduction.md b/docs/src/introduction.md
index 4b58b0a47..0b9c7c9b9 100644
--- a/docs/src/introduction.md
+++ b/docs/src/introduction.md
@@ -22,7 +22,7 @@ Closely linked to the research presented here is the development of a software p
 train!
 ```
 
-So the docstring shows the name of the type or method, in most cases how to call it and then gives some information explaining what it does and potentially hyperlinks to other similar docstrings ([`DataLoader`](@ref) in this case); all of this information is indented by a tab. Docstrings may include other information under subheaders **Arguments** (showing the arguments the method can be supplied with), **Examples** (giving more detailed examples (including results) of how to use the method) and **Implementation** (giving details on how the method is implemented). When we reference a docstring it is always printed in blue (e.g. [`train!`](@ref)), indicating a hyperlink. In addition there is an *index of docstrings* showing all docstrings in chronological order with the associated page number.
+So the docstring shows the name of the type or method, in most cases how to call it and then gives some information explaining what it does and potentially hyperlinks to other similar docstrings ([`TransformerLoss`](@ref) in this case); all of this information is indented by a tab. Docstrings may include other information under subheaders **Arguments** (showing the arguments the method can be supplied with), **Examples** (giving more detailed examples (including results) of how to use the method) and **Implementation** (giving details on how the method is implemented). When we reference a docstring it is always printed in blue (e.g. [`train!`](@ref)), indicating a hyperlink. In addition there is an *index of docstrings* showing all docstrings in chronological order with the associated page number.
 
 Similar to **Library Functions**, which is included in most sections, almost every chapter concludes with a section **Chapter Summary** and an additional section **References** that shows further related reading material. The **Chapter Summary** recaps the important aspects of the corresponding chapter, states again what is new (this may be mathematical or software aspects) and gives information to what other parts of the dissertation the contents of the present chapter are relevant.
 
diff --git a/docs/src/tutorials/symplectic_autoencoder.md b/docs/src/tutorials/symplectic_autoencoder.md
index eb6ba06a0..d21f7aafb 100644
--- a/docs/src/tutorials/symplectic_autoencoder.md
+++ b/docs/src/tutorials/symplectic_autoencoder.md
@@ -15,7 +15,7 @@ Main.remark(raw"As with any neural network we have to make the following choices
 " * Main.indentation * raw"5. choose a *number of epochs*,
 " * Main.indentation * raw"where points 1 and 3 depend on a variable number of hyperparameters.")
 ```
-For the symplectic autoencoder point 1 is done by calling [`SymplecticAutoencoder`](@ref), point 2 is done by calling `NeuralNetwork`, point 3 is done by calling [`TrainingOptimizer`](@extref GeometricOptimizers GeometricOptimizers.TrainingOptimizer) and point 4 is done by calling [`Batch`](@ref).
+For the symplectic autoencoder point 1 is done by calling [`SymplecticAutoencoder`](@ref), point 2 is done by calling `NeuralNetwork`, point 3 is done by calling [`TrainingOptimizer`](@extref GeometricOptimizers GeometricOptimizers.TrainingOptimizer) and point 4 is done by calling `Batch`.
 
 ## The system
 
@@ -85,16 +85,16 @@ sol = integrate(pr, ImplicitMidpoint())
 nothing # hide
 ```
 
-We then put the data in the correct format by calling [`DataLoader`](@ref)[^1]:
+We then put the data in the correct format by calling `DataLoader`[^1]:
 
-[^1]: For more information on [`DataLoader`](@ref) see the [corresponding section](@ref "The Data Loader").
+[^1]: For more information on `DataLoader` see the [corresponding section](@ref "The Data Loader").
 
 ```@example toda_lattice
 dl_cpu = DataLoader(sol; autoencoder = true, suppress_info = true)
 nothing # hide
 ```
 
-Also note that the keyword `autoencoder` was set to true when calling [`DataLoader`](@ref). The keyword argument `supress_info` determines whether data loader provides some additional information on the data it is called on.
+Also note that the keyword `autoencoder` was set to true when calling `DataLoader`. The keyword argument `supress_info` determines whether data loader provides some additional information on the data it is called on.
 
 ## Train the network 
 
diff --git a/docs/src/tutorials/sympnet_tutorial.md b/docs/src/tutorials/sympnet_tutorial.md
index fa1a63954..2e9936539 100644
--- a/docs/src/tutorials/sympnet_tutorial.md
+++ b/docs/src/tutorials/sympnet_tutorial.md
@@ -35,7 +35,7 @@ solution = integrate(ho_problem, ImplicitMidpoint())
 nothing # hide
 ```
 
-We call [`DataLoader`](@ref) in order to conveniently handle the data:
+We call `DataLoader` in order to conveniently handle the data:
 
 ```@example sympnet
 using GeometricMachineLearning # hide
diff --git a/scripts/Project.toml b/scripts/Project.toml
index 26b087aab..bb6e27018 100644
--- a/scripts/Project.toml
+++ b/scripts/Project.toml
@@ -18,6 +18,7 @@ NeuralNetworkParameters = "67f4d93a-60e9-472b-8cdd-1ccf6005724a"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReducedComplexityModeling = "abb4722c-1d8e-48eb-a3bd-45b3116bfaba"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
diff --git a/scripts/reproduction/symplectic_autoencoders/online_sympnet.jl b/scripts/reproduction/symplectic_autoencoders/online_sympnet.jl
index b31fed1df..9e57b379d 100644
--- a/scripts/reproduction/symplectic_autoencoders/online_sympnet.jl
+++ b/scripts/reproduction/symplectic_autoencoders/online_sympnet.jl
@@ -4,6 +4,7 @@ using CairoMakie
 using GeometricProblems.TodaLattice: hodeensemble
 using GeometricMachineLearning
 using GeometricMachineLearning: solve!
+using ReducedComplexityModeling: number_of_batches
 using JLD2
 import Random
 
@@ -41,7 +42,7 @@ const batch_size = 4096
 sae_pairing = AdamOptimizerWithDecay(n_epochs)
 o = TrainingOptimizer(sae_nn.params; sae_pairing...)
 
-println("Number of batches: ", GeometricMachineLearning.number_of_batches(dl, Batch(batch_size)))
+println("Number of batches: ", number_of_batches(dl, Batch(batch_size)))
 
 psd_error = solve!(psd_nn, dl)
 sae_error = train!(sae_nn, o, dl, Batch(batch_size), n_epochs)
diff --git a/scripts/utilities/parametric_data_helpers.jl b/scripts/utilities/parametric_data_helpers.jl
index c9bf7b7d7..32641fc63 100644
--- a/scripts/utilities/parametric_data_helpers.jl
+++ b/scripts/utilities/parametric_data_helpers.jl
@@ -7,12 +7,12 @@
 # `GMLDatasets`' job, https://github.com/JuliaGNI/GMLDatasets.jl/issues/5 tracks moving them there --
 # so treat this file as a staging post.
 
-using GeometricMachineLearning: ParametricDataLoader
+using ReducedComplexityModeling: ParametricDataLoader
 
 """
     turn_parameters_into_correct_format(t, IC)
 
-The vector of `NamedTuple`s to be used by [`ParametricDataLoader`](@ref): one `(t = …,)` per time
+The vector of `NamedTuple`s to be used by `ParametricDataLoader`: one `(t = …,)` per time
 step of the trajectory, repeated once per initial condition of `IC`.
 
 There is one entry fewer than there are time steps in `t`, because a sample is a pair of consecutive
@@ -29,7 +29,7 @@ end
 Turn a `NamedTuple` of ``(q, p)`` data into the two tensors of the correct format.
 
 This is the tricky part, as the structure of the arrays has to conform to the structure of the
-parameters. A [`ParametricDataLoader`](@ref) reads a ``(\mathrm{dim}, nt, m)`` tensor as `m`
+parameters. A `ParametricDataLoader` reads a ``(\mathrm{dim}, nt, m)`` tensor as `m`
 trajectories of `nt` time steps and takes *one* parameter for each of the `m` trajectories, so a
 parameter that depends on the time -- as the time does here -- has to make every sample a trajectory
 of its own: a pair of consecutive states.
@@ -60,7 +60,7 @@ end
 """
     load_time_dependent_harmonic_oscillator_with_parametric_data_loader(qp, t, IC)
 
-The [`ParametricDataLoader`](@ref) of a harmonic oscillator whose parameters depend on the time,
+The `ParametricDataLoader` of a harmonic oscillator whose parameters depend on the time,
 with one `(t = …,)` per sample: `qp` is the `(q, p)` data of the trajectories in `IC`, and `t` their
 time grid.
 """
diff --git a/src/GeometricMachineLearning.jl b/src/GeometricMachineLearning.jl
index 2cf264525..d904a71b1 100644
--- a/src/GeometricMachineLearning.jl
+++ b/src/GeometricMachineLearning.jl
@@ -9,16 +9,14 @@ using AbstractNeuralNetworks
 import NeuralNetworkParameters: NetworkParameters
 using NeuralNetworkParameters: mapstorage
 # `flatten` puts a system's parameters into the input of a network that takes them
-# (`src/parametric.jl`); `parameter_eltype` and `parameterlayout` check their element type and
-# shape, and `map_cotangent` and `storage_gradient` turn a cotangent into a gradient.
-using NeuralNetworkParameters: flatten, parameter_eltype, parameterlayout, map_cotangent,
-                               storage_gradient
+# (`src/parametric.jl`), and `map_cotangent` and `storage_gradient` turn a cotangent into a
+# gradient.
+using NeuralNetworkParameters: flatten, map_cotangent, storage_gradient
 using ChainRulesCore
 using GeometricBase
 # `integrate` is GeometricBase's generic, and importing it is what makes the trajectory methods
 # below methods of that function rather than of a second one of the name.
 import GeometricBase: integrate
-using GeometricSolutions: GeometricSolution, EnsembleSolution, TimeSeries
 using GeometricEquations: HODEProblem, HODEEnsemble
 using KernelAbstractions
 using LinearAlgebra
@@ -78,6 +76,19 @@ import AbstractNeuralNetworks: networkbackend
 # The input and output types of a network: `NamedTupleOfArrays{T, (:q, :p)}` is data in `(q, p)`
 # form, and `ArrayOrNamedTuple{T}` is an array or a named tuple of arrays with any keys.
 using AbstractNeuralNetworks: NamedTupleOfArrays, ArrayOrNamedTuple
+# The data loader is `ReducedComplexityModeling`'s, and this package defines no second `DataLoader`,
+# `Batch` or `train!`: the three names are imported, so that `using GeometricMachineLearning`
+# reaches the objects the other package exports and the `train!` below is a method of that generic.
+# `ParametricDataLoader` is imported and not re-exported: the parametric family is this package's,
+# the loader is not, so a caller names it here or takes it from its owner.
+import ReducedComplexityModeling: DataLoader, Batch, train!, ParametricDataLoader
+# The loader's routines this package calls: `assign_output_estimate` crops a classifier's output,
+# `convert_input_and_batch_indices_to_array` gathers one minibatch for `train!`, and
+# `data_tensors_from_geometric_solution` puts a solution into the tensors that the reduced system's
+# errors are computed on.
+using ReducedComplexityModeling: assign_output_estimate,
+                                 convert_input_and_batch_indices_to_array,
+                                 data_tensors_from_geometric_solution
 # `save` and `load` are `NeuralNetworkParameters`' generics, which `AbstractNeuralNetworks`
 # re-binds; reach for them where they are defined.
 import NeuralNetworkParameters: save, load
@@ -95,16 +106,13 @@ export NetworkParameters
 export σ, sigmoid, softmax
 
 # `GeometricBase` defines `description` and does not export it, so re-exporting it needs the import.
-# `StateVariable` is taken from here, its owner: `GeometricSolutions` holds it only through a bare
-# `using GeometricBase`, and importing a name from a module that does not own it is what
-# `check_all_explicit_imports_via_owners` reports.
+# `StateVariable` is imported from here too, from the module that defines it, which is what
+# `check_all_explicit_imports_via_owners` asks of an explicit import.
 import GeometricBase: description, StateVariable
 export description
 
 include("utils.jl")
 
-include("data_loader/data_loader.jl")
-
 # The structured matrix types come from GeometricOptimizers, which also carries every method GML
 # needs on them; this package defines none of its own. They are re-exported here, so that `using
 # GeometricMachineLearning` on its own still gives a caller the matrix types its layers are
@@ -219,17 +227,15 @@ include("loss/losses.jl")
 export AbstractPullback, ZygotePullback
 include("pullbacks/zygote_pullback.jl")
 
-export DataLoader
-export Batch, train!
-include("data_loader/tensor_assign.jl")
-include("data_loader/batch.jl")
-include("data_loader/parametric_data_loader.jl")
+# The three names a training loop is written with, re-exported from their owner: `DataLoader`,
+# `Batch` and `train!` are `ReducedComplexityModeling`'s, and the method below on `train!` extends
+# that generic.
+export DataLoader, Batch, train!
 include("data_loader/train.jl")
-# after the loader, because a constructor of `ParametricResNet` takes a `ParametricDataLoader`
 include("architectures/parametric_resnet.jl")
 
 # The training of a network that takes the parameters of the system beside its input.
-public apply_parametric, ParametricLoss, ParametricDataLoader
+public apply_parametric, ParametricLoss
 
 include("reduced_system/reduced_system.jl")
 
diff --git a/src/architectures/forced_sympnet.jl b/src/architectures/forced_sympnet.jl
index ac7c9ed6b..b431dcf89 100644
--- a/src/architectures/forced_sympnet.jl
+++ b/src/architectures/forced_sympnet.jl
@@ -29,7 +29,7 @@ gradient layers and one forcing layer per step that changes ``p``.
 dimension raises an `ArgumentError` at construction.
 
 There exists an additional constructor that can be called by supplying an instance of
-[`DataLoader`](@ref).
+`DataLoader`.
 
 # Arguments
 
diff --git a/src/architectures/parametric_resnet.jl b/src/architectures/parametric_resnet.jl
index f8bc1c23a..b7f007ad3 100644
--- a/src/architectures/parametric_resnet.jl
+++ b/src/architectures/parametric_resnet.jl
@@ -30,8 +30,8 @@ width, and the linear output layer applies the identity activation.
 ParametricResNet(dl, n_blocks, width)
 ```
 
-where `dl` is an instance of [`DataLoader`](@ref) and `width` defaults to `dl.input_dim`, or an
-instance of [`ParametricDataLoader`](@ref). The latter gives the system dimension, the rows of its
+where `dl` is an instance of `DataLoader` and `width` defaults to `dl.input_dim`, or an
+instance of `ParametricDataLoader`. The latter gives the system dimension, the rows of its
 data or of `q` and `p` together, which is also the default of `width`, and the system parameters,
 whose first `μ` the architecture keeps for its shape; it therefore takes no `parameters` keyword.
 
diff --git a/src/architectures/sympnet.jl b/src/architectures/sympnet.jl
index 9ecc269c5..abbc6e935 100644
--- a/src/architectures/sympnet.jl
+++ b/src/architectures/sympnet.jl
@@ -22,7 +22,7 @@ const la_activation_default = tanh
 
 Make an ``LA``-SympNet with dimension ``d.``
 
-There exists an additional constructor that can be called by supplying an instance of [`DataLoader`](@ref).
+There exists an additional constructor that can be called by supplying an instance of `DataLoader`.
 
 # Examples
 
@@ -79,7 +79,7 @@ const g_init_upper_default = true
 
 Make a ``G``-SympNet with dimension ``d.``
 
-There exists an additional constructor that can be called by supplying an instance of [`DataLoader`](@ref) (see [`LASympNet`](@ref) for an example of using this constructor).
+There exists an additional constructor that can be called by supplying an instance of `DataLoader` (see [`LASympNet`](@ref) for an example of using this constructor).
 
 # Arguments
 
diff --git a/src/architectures/transformer_integrator.jl b/src/architectures/transformer_integrator.jl
index 1d539bc68..5d50a39b6 100644
--- a/src/architectures/transformer_integrator.jl
+++ b/src/architectures/transformer_integrator.jl
@@ -132,3 +132,13 @@ function integrate(nn::NeuralNetwork{<:TransformerIntegrator},
 
     valuation[:, 1:n_points]
 end
+
+# The three-argument form a training script writes: the prediction window comes from the `Batch`,
+# so that `integrate(nn, ics, batch)` gives the trajectory `integrate(nn, ics; prediction_window =
+# batch.prediction_window)` gives. `Batch` is `ReducedComplexityModeling`'s, so the annotation on
+# `nn` is what keeps this a method of this package's and not piracy on another package's type.
+function integrate(
+        nn::NeuralNetwork{<:TransformerIntegrator}, ics, batch::Batch{:Transformer};
+        n_points = 100)
+    integrate(nn, ics; n_points = n_points, prediction_window = batch.prediction_window)
+end
diff --git a/src/architectures/transformer_neural_network.jl b/src/architectures/transformer_neural_network.jl
index 3fe5de48e..e893b3697 100644
--- a/src/architectures/transformer_neural_network.jl
+++ b/src/architectures/transformer_neural_network.jl
@@ -1,9 +1,9 @@
 @doc raw"""
     ClassificationTransformer(dl)
 
-Make an instance of the `ClassificationTransformer` based on an instance of [`DataLoader`](@ref).
+Make an instance of the `ClassificationTransformer` based on an instance of `DataLoader`.
 
-This is a transformer neural network for classification purposes. `dl` has to store its input as a tensor and its target as a tensor of labels; [`GMLDatasets`](https://github.com/JuliaGNI/GMLDatasets.jl) builds such a [`DataLoader`](@ref) from an image data set, but any classification problem works.
+This is a transformer neural network for classification purposes. `dl` has to store its input as a tensor and its target as a tensor of labels; [`GMLDatasets`](https://github.com/JuliaGNI/GMLDatasets.jl) builds such a `DataLoader` from an image data set, but any classification problem works.
 
 # Arguments
 
diff --git a/src/data_loader/batch.jl b/src/data_loader/batch.jl
deleted file mode 100644
index 8a773d1e6..000000000
--- a/src/data_loader/batch.jl
+++ /dev/null
@@ -1,277 +0,0 @@
-@doc raw"""
-    Batch
-
-`Batch` is a struct whose functor acts on an instance of `DataLoader` to produce a sequence of training samples for training for one epoch. 
-
-See [`Batch(::Int)`](@ref) and [`Batch(::Int, ::Int, ::Int)`](@ref) for the different constructors.
-
-# The functor 
-
-An instance of `Batch` can be called on an instance of `DataLoader` to produce a sequence of samples that contain all the input data, i.e. for training for one epoch. 
-
-The output of applying `batch:Batch` to `dl::DataLoader` is a vector of vectors of integers. Each of these vectors contains two integers: the first is the *time index* and the second one is the *parameter index*.
-
-# Examples
-
-Consider the following example for drawing batches of size 2 for an instance of `DataLoader` constructed with a vector:
-
-```jldoctest
-using GeometricMachineLearning
-import Random
-
-rng = Random.TaskLocalRNG()
-Random.seed!(rng, 123)
-
-dl = DataLoader(rand(rng, 5); autoencoder = true)
-batch = Batch(2)
-
-batches = batch(dl)
-
-@assert length(batches) == 3
-@assert length.(batches) == [2, 2, 1]
-@assert sort(collect(union(batches...))) == [(1, i) for i in 1:5]
-
-true
-
-# output
-
-[ Info: You have provided a matrix as input. The axes will be interpreted as (i) system dimension and (ii) number of parameters.
-true
-```
-
-Here the first index is always 1 (the time dimension). We get a total number of 3 batches. 
-The last batch is only of size 1 because we *sample without replacement*.
-Also see the docstring for [`DataLoader(::AbstractVector)`](@ref).
-"""
-struct Batch{BatchType}
-    batch_size::Int
-    seq_length::Int
-    prediction_window::Int
-end
-
-# add an additional constructor for `TransformerLoss` taking batch as input.
-function TransformerLoss(batch::Batch{:Transformer})
-    TransformerLoss(batch.seq_length, batch.prediction_window)
-end
-
-function integrate(nn::NeuralNetwork, ics, batch::Batch{:Transformer}; n_points = 100)
-    integrate(nn, ics; n_points = n_points, prediction_window = batch.prediction_window)
-end
-
-@doc raw"""
-    Batch(batch_size)
-
-Make an instance of `Batch` for a specific batch size.
-
-This is, among others, used to train neural networks of [`NeuralNetworkIntegrator`](@ref) type (as opposed to [`TransformerIntegrator`](@ref)).
-"""
-function Batch(batch_size::Int)
-    Batch{:FeedForward}(batch_size, 1, 1)
-end
-
-@doc raw"""
-    Batch(batch_size, seq_length)
-
-Make an instance of `Batch` for a specific batch size and a sequence length.
-
-This is used to train neural networks of [`TransformerIntegrator`](@ref) type.
-
-Optionally the prediction window can also be specified by calling:
-
-```jldoctest
-using GeometricMachineLearning
-
-batch_size = 2
-seq_length = 3
-prediction_window = 2
-
-Batch(batch_size, seq_length, prediction_window)
-
-# output
-
-Batch{:Transformer}(2, 3, 2)
-```
-
-Note that here the batch is of type `:Transformer`.
-"""
-function Batch(batch_size::Int, seq_length::Int, prediction_window::Int = seq_length)
-    Batch{:Transformer}(batch_size, seq_length, prediction_window)
-end
-
-# The classifier's batch: one draw per sample, with the whole input sequence and the whole output.
-# The loader's fifth type parameter is `Int`, because both lengths come from a loader with an output.
-function Batch(batch_size::Int,
-        dl::DataLoader{T, AT, OT, K, Int}) where {T, AT, OT <: AbstractArray, K}
-    Batch(batch_size, dl.input_time_steps, dl.output_time_steps)
-end
-
-@doc raw"""
-    number_of_batches(dl, batch)
-
-Compute the number of batches.
-
-One epoch draws every start position of every parameter index once, in minibatches of the batch
-size, so the count is `cld(n · n_params, batch_size)` with `n` the number of start positions of one
-parameter index. No two indices are drawn twice, and the last minibatch holds what is left.
-
-# Examples
-
-```jldoctest
-using GeometricMachineLearning
-using GeometricMachineLearning: number_of_batches
-import Random
-
-Random.seed!(123)
-
-dat = [1, 2, 3, 4, 5]
-dl₁ = DataLoader(dat; autoencoder = false, suppress_info = true) # time series-like
-dl₂ = DataLoader(dat; autoencoder = true, suppress_info = true) # autoencoder-like
-batch = Batch(3)
-
-batches₁ = batch(dl₁)
-batches₂ = batch(dl₂)
-
-@assert number_of_batches(dl₁, batch) == 2
-@assert number_of_batches(dl₂, batch) == 2
-@assert length.(batches₁) == [3, 1]
-@assert length.(batches₂) == [3, 2]
-@assert sort(collect(union(batches₁...))) == [(i, 1) for i in 1:4]
-@assert sort(collect(union(batches₂...))) == [(1, i) for i in 1:5]
-
-true
-
-# output
-
-true
-```
-
-Here we see that in the *autoencoder case* that last minibatch has an additional element.
-"""
-number_of_batches(dl::DataLoader, batch::Batch) = cld(
-    number_of_start_positions(dl, batch) * dl.n_params, batch.batch_size)
-
-# The start positions of one parameter index: the input window must fit the input's time axis, and
-# the target must fit its own. The target is the output at offset `0` where the loader holds one,
-# the input at offset `seq_length` where it is a time series without an output, and there is none
-# where the loader is a set of samples without an output. A `Batch(batch_size)` pairs the input and
-# the output step by step, so it needs one length on the two time axes.
-function number_of_start_positions(dl::DataLoader{T, AT, OT, K},
-        batch::Batch) where {T, AT, OT, K}
-    input_window = dl.input_time_steps - batch.seq_length + 1
-    target_window = if OT <: AbstractArray
-        if batch isa Batch{:FeedForward} && dl.input_time_steps != dl.output_time_steps
-            throw(ArgumentError(
-                "`Batch(batch_size)` pairs the input and the output time step by time step, " *
-                "and the `DataLoader` has $(dl.input_time_steps) time steps and " *
-                "$(dl.output_time_steps) in its output; `Batch(batch_size, dl)` draws the " *
-                "whole sequence of each sample."))
-        end
-        dl.output_time_steps - batch.prediction_window + 1
-    elseif K === :TimeSeries
-        input_window - batch.prediction_window
-    else
-        input_window
-    end
-    check_window(min(input_window, target_window), dl, batch)
-end
-
-# A loader or a batch that leaves no minibatch is a caller's mistake: an epoch of none trains
-# nothing and says nothing, so the numbers that decide it are named here. The message is built only
-# where it is thrown, so the check itself allocates nothing.
-function check_window(n::Int, dl::DataLoader, batch::Batch)
-    n ≥ 1 && dl.n_params ≥ 1 || throw(ArgumentError(
-        "the `DataLoader` has $(dl.input_time_steps) time steps" *
-        (dl.output === nothing ? "" : " and $(dl.output_time_steps) in its output") *
-        ", with `n_params = $(dl.n_params)`, and the `Batch` asks for a sequence length of " *
-        "$(batch.seq_length) and a prediction window of $(batch.prediction_window), which " *
-        "leaves no minibatch."))
-    n
-end
-
-# The index set of one epoch: the product of the shuffled start positions and parameter indices, cut
-# into minibatches of the batch size.
-function batch_over_two_axes(batch::Batch, number_columns::Int, third_dim::Int)
-    complete_indices = vec(
-        collect(Iterators.product(shuffle(1:number_columns), shuffle(1:third_dim))))
-    map(collect, Iterators.partition(complete_indices, batch.batch_size))
-end
-
-function (batch::Batch)(dl::DataLoader)
-    batch_over_two_axes(batch, number_of_start_positions(dl, batch), dl.n_params)
-end
-
-@kernel function gather_kernel!(gathered, data, indices, offset)
-    i, j, k = @index(Global, NTuple)
-
-    gathered[i, j, k] = data[i, indices[1, k] + offset + j - 1, indices[2, k]]
-end
-
-# A kernel cannot index a `Vector{Tuple{Int, Int}}`, so the batch indices go in as a `2 × n` matrix:
-# the host builds it and the backend takes it in one copy.
-function convert_vector_of_tuples_to_matrix(backend::Backend,
-        batch_indices_tuple::Vector{Tuple{Int, Int}})
-    map_to_new_backend([index[i] for i in 1:2, index in batch_indices_tuple], backend)
-end
-
-# The gather kernel, applied to a tensor or to the `q` and `p` of a `(q, p)` pair. The gathered
-# tensor holds `width` positions from `offset` on, for every drawn start position and parameter
-# index, and has the element type of the data it reads — so one integer target stays integer.
-function gather(data::AbstractArray{T, 3}, indices::AbstractMatrix{Int}, offset::Int,
-        width::Int) where {T}
-    backend = networkbackend(data)
-    gathered = KernelAbstractions.allocate(
-        backend, T, size(data, 1), width, size(indices, 2))
-    gather_kernel!(backend)(gathered, data, indices, offset, ndrange = size(gathered))
-    gathered
-end
-
-function gather(data::NamedTuple{(:q, :p), Tuple{AT, AT}}, indices::AbstractMatrix{Int},
-        offset::Int, width::Int) where {AT <: AbstractArray{<:Any, 3}}
-    map(tensor -> gather(tensor, indices, offset, width), data)
-end
-
-@doc raw"""
-    convert_input_and_batch_indices_to_array(dl, batch, batch_indices)
-
-Assign batch data based on (i) input and (ii) batch indices.
-
-# Examples
-
-```jldoctest
-using GeometricMachineLearning
-
-dl = DataLoader([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9];
-    autoencoder = true, suppress_info = true)
-batch = Batch(3)
-batch_indices = [(1, 1), (1, 3), (1, 5)]
-
-GeometricMachineLearning.convert_input_and_batch_indices_to_array(dl, batch, batch_indices)
-
-# output
-
-1×1×3 Array{Float64, 3}:
-[:, :, 1] =
- 0.1
-
-[:, :, 2] =
- 0.3
-
-[:, :, 3] =
- 0.5
-```
-"""
-function convert_input_and_batch_indices_to_array(
-        dl::DataLoader{T, AT, OT, K}, batch::Batch,
-        batch_indices_tuple::Vector{Tuple{Int, Int}}) where {T, AT, OT, K}
-    indices = convert_vector_of_tuples_to_matrix(networkbackend(dl), batch_indices_tuple)
-    input = gather(dl.input, indices, 0, batch.seq_length)
-    # a loader that holds an output is a set of samples whatever its kind: the target is its output;
-    # without one, a time series predicts its own continuation and a set of samples has no target
-    if OT <: AbstractArray
-        (input, gather(dl.output, indices, 0, batch.prediction_window))
-    elseif K === :TimeSeries
-        (input, gather(dl.input, indices, batch.seq_length, batch.prediction_window))
-    else
-        input
-    end
-end
diff --git a/src/data_loader/data_loader.jl b/src/data_loader/data_loader.jl
deleted file mode 100644
index 01985b98e..000000000
--- a/src/data_loader/data_loader.jl
+++ /dev/null
@@ -1,447 +0,0 @@
-@doc raw"""
-    DataLoader(data)
-
-Make an instance based on a data set.
-
-This is designed to make training convenient.
-
-# Fields of `DataLoader`
-
-The fields of the `DataLoader` struct are the following:
-- `input`: The input data with axes (i) system dimension, (ii) number of time steps and (iii) number of parameters.
-- `output`: The tensor that contains the output (supervised learning) - this may be of type `Nothing` if the constructor is only called with one tensor (unsupervised learning).
-- `input_dim`: The *dimension* of the system, i.e. what is taken as input by a regular neural network.
-- `input_time_steps`: The length of the entire time series (length of the second axis).
-- `n_params`: The number of parameters that are present in the data set (length of third axis)
-- `output_dim`: The dimension of the output tensor (first axis). If `output` is of type `Nothing`, then this is also of type `Nothing`.
-- `output_time_steps`: The size of the second axis of the output tensor. If `output` is of type `Nothing`, then this is also of type `Nothing`.
-
-# Implementation
-
-Even though `DataLoader` can be called with inputs of various forms, internally it always stores tensors with three axes.
-
-```jldoctest
-using GeometricMachineLearning
-
-data = [1 2 3; 4 5 6]
-dl = DataLoader(data)
-dl.input
-
-# output
-
-[ Info: You have provided a matrix as input. The axes will be interpreted as (i) system dimension and (ii) number of time steps.
-2×3×1 Array{Int64, 3}:
-[:, :, 1] =
- 1  2  3
- 4  5  6
-```
-"""
-struct DataLoader{T, AT <: Union{NamedTuple, AbstractArray{T}},
-    OT <: Union{AbstractArray, Nothing}, K, S <: Union{Int, Nothing}}
-    input::AT
-    output::OT
-    input_dim::Int
-    input_time_steps::Int
-    n_params::Int
-    output_dim::S
-    output_time_steps::S
-end
-
-# The fifth parameter is the type of the two output fields, which is `Nothing` for a loader without
-# an output and `Int` with one. A caller that names the first four — `DataLoader{T, AT, OT,
-# :TimeSeries}(…)`, which is how `GMLDatasets` builds a loader — gets it filled in from the fields,
-# and the five counts are converted to `Int` whatever integer type the caller passes.
-function DataLoader{T, AT, OT, K}(input::AT, output::OT, input_dim, input_time_steps,
-        n_params, output_dim, output_time_steps) where {T, AT, OT, K}
-    S = output_dim === nothing ? Nothing : Int
-    DataLoader{T, AT, OT, K, S}(input, output, input_dim, input_time_steps, n_params,
-        output_dim, output_time_steps)
-end
-
-# The two decisions `autoencoder` makes at every constructor that takes a data set: the kind it
-# builds, and the two lengths of a matrix or a `(q, p)` pair read as a set of samples, whose second
-# axis is then the parameter index and whose time axis has length one.
-data_set_kind(autoencoder::Bool) = autoencoder ? :RegularData : :TimeSeries
-data_set_shape(autoencoder::Bool, n::Int) = autoencoder ? (1, n) : (n, 1)
-
-"""
-    DataLoader(data::AbstractArray{<:Number, 3})
-
-Make an instance of DataLoader for data that are in tensor format.
-
-# Arguments
-
-There are two optional keyword arguments:
-- `autoencoder = false` and
-- `suppress_info = false`.
-
-By default the data are stored as `TimeSeries` type. If you want to train an [`AutoEncoder`](@ref) with your data call:
-
-```julia
-    DataLoader(data; autoencoder = true)
-```
-
-The default is equivalent to `autoencoder = false`.
-
-By default we have:
-
-```jldoctest
-using GeometricMachineLearning
-
-data = [ 1;  2;  3;;
-         4;  5;  6;;;
-         7;  8;  9;;
-        10; 11; 12]
-
-DataLoader(data)
-
-# output
-
-┌ Info: You have provided a tensor with three axes as input. They will be interpreted as
-└  (i) system dimension, (ii) number of time steps and (iii) number of params.
-DataLoader{Int64, Array{Int64, 3}, Nothing, :TimeSeries, Nothing}([1 4; 2 5; 3 6;;; 7 10; 8 11; 9 12], nothing, 3, 2, 2, nothing, nothing)
-```
-
-But if we write
-
-```jldoctest
-using GeometricMachineLearning
-
-data = [ 1;  2;  3;;
-         4;  5;  6;;;
-         7;  8;  9;;
-        10; 11; 12]
-
-DataLoader(data; suppress_info = true)
-
-# output
-
-DataLoader{Int64, Array{Int64, 3}, Nothing, :TimeSeries, Nothing}([1 4; 2 5; 3 6;;; 7 10; 8 11; 9 12], nothing, 3, 2, 2, nothing, nothing)
-```
-the `@info` statement is not printed.
-"""
-Base.@constprop :aggressive function DataLoader(
-        data::AbstractArray{T, 3}; autoencoder::Bool = false, suppress_info = false) where {T}
-    if !suppress_info
-        @info "You have provided a tensor with three axes as input. They will be interpreted as \n (i) system dimension, (ii) number of time steps and (iii) number of params."
-    end
-    input_dim, input_time_steps, n_params = size(data)
-
-    DataLoader{T, typeof(data), Nothing, data_set_kind(autoencoder)}(
-        data, nothing, input_dim, input_time_steps, n_params, nothing, nothing)
-end
-
-"""
-    DataLoader(data::AbstractMatrix)
-
-Make an instance of [`DataLoader`](@ref) based on a matrix.
-
-# Arguments
-
-See [`DataLoader(::AbstractArray{<:Number, 3})`](@ref) for details.
-
-# Implementation
-
-The second axis of the matrix is the number of time steps, and internally the data are reshaped to
-a tensor of shape `(size(data)..., 1)`: one trajectory of one parameter. With `autoencoder = true`
-the matrix is a set of samples instead, so that axis is the parameter index and the time axis has
-length one.
-"""
-Base.@constprop :aggressive function DataLoader(
-        data::AbstractMatrix{T}; autoencoder::Bool = false, suppress_info = false) where {T}
-    if !suppress_info
-        @info "You have provided a matrix as input. The axes will be interpreted as (i) system dimension and (ii) $(autoencoder ? "number of parameters" : "number of time steps")."
-    end
-    input_dim, second_dim = size(data)
-    time_steps, n_params = data_set_shape(autoencoder, second_dim)
-    DataLoader(reshape(data, input_dim, time_steps, n_params);
-        autoencoder = autoencoder, suppress_info = true)
-end
-
-"""
-    DataLoader(data::AbstractVector)
-
-Make an instance of `DataLoader` based on a vector.
-
-# Extended help
-
-If the input to `DataLoader` is a vector, it is assumed that this vector represents one-dimensional time-series data and is therefore processed as:
-
-```julia
-    DataLoader(data::AbstractVector; autoencoder=false) = DataLoader(reshape(data, 1, length(data)); autoencoder = autoencoder)
-```
-"""
-function DataLoader(data::AbstractVector; autoencoder::Bool = false, suppress_info = false)
-    DataLoader(reshape(data, 1, length(data)); autoencoder = autoencoder, suppress_info = suppress_info)
-end
-
-# A loader that holds an output is a set of samples whatever its kind: its target is the output
-# itself and never the input shifted by the sequence length. The output's element type is the
-# input's or an integer, so a real input can carry an integer one-hot target — a classifier's data,
-# and `accuracy`'s argument — and any other pair of element types is a `MethodError`.
-function DataLoader(input::AbstractArray{T, 3}, output::AbstractArray{T1, 3};
-        suppress_info = false) where {T, T1 <: Union{T, Integer}}
-    @assert size(input, 3) == size(output, 3)
-    if !suppress_info
-        @info "You have provided an input and an output."
-    end
-
-    DataLoader{T, typeof(input), typeof(output), :RegularData}(
-        input, output, size(input, 1), size(input, 2),
-        size(input, 3), size(output, 1), size(output, 2))
-end
-
-function DataLoader(input::AbstractMatrix{T}, output::AbstractMatrix{T}; suppress_info = false) where {T}
-    DataLoader(reshape(input, size(input)..., 1),
-        reshape(output, size(output)..., 1); suppress_info = suppress_info)
-end
-
-function DataLoader(input::AbstractVector{T}, output::AbstractVector{T}; suppress_info = false) where {T}
-    DataLoader(reshape(input, 1, length(input)),
-        reshape(output, 1, length(output)); suppress_info = suppress_info)
-end
-
-@doc raw"""
-    DataLoader(data::NamedTuple{(:q, :p), Tuple{AT, AT}})
-
-Make an instance of `DataLoader` based on ``(q, p)`` data.
-
-# Implementation
-
-In this case the field `input_dim` of `DataLoader` is interpreted as the sum of the ``q``- and ``p``-dimensions, i.e. if ``q`` and ``p`` both evolve on ``\mathbb{R}^n``, then `input_dim` is ``2n``.
-
-Apart from this the input is treated similarly as if it were an `Array`, i.e. everything is converted to tensors internally.
-See e.g. [`DataLoader{::AbstractArray{<:Number, 3}}`](@ref).
-"""
-Base.@constprop :aggressive function DataLoader(data::NamedTuple{(:q, :p), Tuple{AT, AT}};
-        autoencoder::Bool = false,
-        suppress_info = false) where {T, AT <: AbstractMatrix{T}}
-    if !suppress_info
-        @info "You have provided a NamedTuple with keys q and p; the data are matrices. This is interpreted as *symplectic data*."
-    end
-    dim2, second_dim = size(data.q)
-    time_steps, n_params = data_set_shape(autoencoder, second_dim)
-    DataLoader(map(tensor -> reshape(tensor, dim2, time_steps, n_params), data);
-        autoencoder = autoencoder, suppress_info = true)
-end
-
-Base.@constprop :aggressive function DataLoader(data::NamedTuple{(:q, :p), Tuple{AT, AT}};
-        autoencoder::Bool = false,
-        suppress_info = false) where {T, AT <: AbstractArray{T, 3}}
-    if !suppress_info
-        @info "You have provided a NamedTuple with keys q and p; the data are tensors. This is interpreted as *symplectic data*."
-    end
-    dim2, time_steps, n_params = size(data.q)
-
-    DataLoader{T, typeof(data), Nothing, data_set_kind(autoencoder)}(
-        data, nothing, dim2 * 2, time_steps, n_params, nothing, nothing)
-end
-
-function DataLoader(data::NamedTuple{(:q, :p), Tuple{VT, VT}};
-        autoencoder::Bool = false, suppress_info = false) where {VT <: AbstractVector}
-    DataLoader(
-        (q = reshape(data.q, 1, length(data.q)), p = reshape(data.p, 1, length(data.p)));
-        autoencoder = autoencoder, suppress_info = suppress_info)
-end
-
-# The 1-based tensor of one coordinate over a vector of trajectories — one solution, or the members
-# of an ensemble — with the axes `(system dimension, time steps, trajectories)`. The `q` of a
-# solution carries the `DataSeries` time axis `0:n`, so the states are copied in one at a time.
-function data_tensor(trajectories::AbstractVector)
-    states = first(trajectories)
-    data = Array{eltype(first(states))}(
-        undef, length(first(states)), length(states), length(trajectories))
-    for (k, trajectory) in enumerate(trajectories)
-        for (t, state) in enumerate(trajectory)
-            data[:, t, k] = state
-        end
-    end
-    data
-end
-
-function data_tensors_from_geometric_solution(
-        solution::GeometricSolution{
-        <:Number, <:Number, <:Any, <:NamedTuple{(:t, :q, :p, :q̇, :ṗ)}})
-    (q = data_tensor([solution.q]), p = data_tensor([solution.p]))
-end
-
-function data_tensors_from_geometric_solution(solution::GeometricSolution{
-        <:Number, <:Number,
-        <:Any, <:Union{NamedTuple{(:t, :q, :v)}, NamedTuple{(:t, :q, :q̇)}}})
-    data_tensor([solution.q])
-end
-
-"""
-    DataLoader(solution)
-
-Make an instance of `DataLoader` for a `GeometricSolution`.
-
-`GeometricSolution`s are imported from the the package [`GeometricSolutions.jl`](https://github.com/JuliaGNI/GeometricSolutions.jl).
-
-# Arguments
-
-This functor for [`DataLoader`](@ref) also has the keyword arguments
-- `autoencoder = false` and
-- `suppress_info = false`.
-See the docstring for [`DataLoader(::AbstractArray{<:Number, 3})`](@ref).
-
-# Implementation
-
-Internally this stores the data as a tensor where the third axis has length 1.
-"""
-function DataLoader(solution::GeometricSolution{T, T2, TT, NT};
-        autoencoder::Bool = false, suppress_info = false) where {T <: Number, T2 <: Number,
-        TT <: TimeSeries{T2}, TuT,
-        NT <: Union{NamedTuple{(:t, :q, :p, :q̇, :ṗ), TuT},
-            NamedTuple{(:t, :q, :v), TuT}, NamedTuple{(:t, :q, :q̇), TuT}}}
-    DataLoader(data_tensors_from_geometric_solution(solution);
-        autoencoder = autoencoder, suppress_info = suppress_info)
-end
-
-function DataLoader(ensemble_solution::EnsembleSolution{T, T1, Vector{ST}};
-        autoencoder::Bool = false, suppress_info = false) where {T, T1, TuT,
-        TT <: TimeSeries{T1}, ST <: Union{
-            GeometricSolution{T, T1, TT, NamedTuple{(:t, :q, :v), TuT}},
-            GeometricSolution{T, T1, TT, NamedTuple{(:t, :q, :q̇), TuT}}}}
-    DataLoader(data_tensor([s.q for s in ensemble_solution.s]);
-        autoencoder = autoencoder, suppress_info = suppress_info)
-end
-
-"""
-    DataLoader(ensemble_solution)
-
-Make an instance of `DataLoader` for a `EnsembleSolution`.
-
-`EnsembleSolution`s are imported from the the package [`GeometricSolutions.jl`](https://github.com/JuliaGNI/GeometricSolutions.jl).
-
-# Arguments
-
-This functor for `DataLoader` also has the keyword arguments
-- `autoencoder = false` and
-- `suppress_info = false`.
-
-See the docstring for [`DataLoader(::AbstractArray{<:Number, 3})`](@ref).
-
-# Implementation
-
-Internally this stores the data as a tensor where the third axis has length equal to the number of solutions in the ensemble.
-"""
-function DataLoader(ensemble_solution::EnsembleSolution{T, T1, Vector{ST}};
-        autoencoder::Bool = false, suppress_info = false) where {T, T1, TuT,
-        TT <: TimeSeries{T1}, ST <: GeometricSolution{T, T1, TT,
-            NamedTuple{(:t, :q, :p, :q̇, :ṗ), TuT}}}
-    DataLoader(
-        (q = data_tensor([s.q for s in ensemble_solution.s]),
-            p = data_tensor([s.p for s in ensemble_solution.s]));
-        autoencoder = autoencoder, suppress_info = suppress_info)
-end
-
-# An array already on the backend is returned as it is. `Base.copyto!` carries a host array to any
-# backend; `KernelAbstractions.copyto!` has a method for the CPU alone, and a device backend such as
-# `JLArray`'s has none.
-function map_to_new_backend(input::AbstractArray{T},
-        backend::KernelAbstractions.Backend) where {T}
-    networkbackend(input) == backend && return input
-    Base.copyto!(KernelAbstractions.allocate(backend, T, size(input)...), input)
-end
-
-# `map` over a named tuple of arrays, such as `(q, p)` data, reaches each of its tensors; over a
-# plain array it would map the array's elements instead, which is why an array has a method of its
-# own.
-function map_to_new_backend(input::NamedTupleOfArrays, backend::KernelAbstractions.Backend)
-    map(tensor -> map_to_new_backend(tensor, backend), input)
-end
-
-map_to_type(input::AbstractArray, ::Type{T}) where {T} = T.(input)
-
-function map_to_type(input::NamedTupleOfArrays, ::Type{T}) where {T}
-    map(tensor -> map_to_type(tensor, T), input)
-end
-
-@doc raw"""
-    DataLoader(dl, backend)
-
-Make a new instance of `DataLoader` based on an existing instance of `DataLoader` for a new `backend`.
-
-This allocates new memory of the same size as is used for the original `dl` and then copies the data.
-
-By default the data type remains unchanged, i.e. `eltype(DataLoader(dl, backend)) == eltype(dl)` is `true`.
-
-If you want to change the data type write e.g.
-
-```julia
-    DataLoader(dl, backend, Float32)
-```
-
-# Arguments
-
-There is an optional keyword argument
-- `autoencoder = nothing`.
-
-By default this inherits the autoencoder property form `dl`.
-
-See the docstring for [`DataLoader(data::AbstractArray{<:Number, 3})`](@ref).
-"""
-function DataLoader(dl::DataLoader{T1, <:ArrayOrNamedTuple, Nothing, K},
-        backend::KernelAbstractions.Backend = networkbackend(dl),
-        ::Type{T} = T1;
-        autoencoder::Union{Nothing, Bool} = nothing
-) where {T1, K, T}
-    kind = isnothing(autoencoder) ? K : data_set_kind(autoencoder)
-    input = T == T1 ? dl.input : map_to_type(dl.input, T)
-    input = map_to_new_backend(input, backend)
-
-    DataLoader{T, typeof(input), Nothing, kind}(
-        input, nothing, dl.input_dim, dl.input_time_steps, dl.n_params, nothing, nothing)
-end
-
-function DataLoader(dl::DataLoader, ::Type{T}; kwargs...) where {T}
-    DataLoader(dl, networkbackend(dl), T; kwargs...)
-end
-
-@doc raw"""
-    accuracy(model, ps, dl)
-
-Compute the accuracy of a neural network classifier.
-
-This needs an instance of [`DataLoader`](@ref) that stores the *test data*. The result is one minus
-the fraction of errors, where a wrongly labelled sample is one error and a sample that carries no
-label is half an error. Its type is `float(T)`, for `T` the element type of the loader's input.
-"""
-function accuracy(model::Chain,
-        ps::NetworkParameters,
-        dl::DataLoader{T, AT, BT}) where {
-        T, T1 <: Integer, AT <: AbstractArray{T}, BT <: AbstractArray{T1}}
-    output_tensor = model(dl.input, ps)
-    output_estimate = assign_output_estimate(output_tensor, dl.output_time_steps)
-    backend = networkbackend(output_estimate)
-    tensor_of_maximum_elements = KernelAbstractions.zeros(backend, T1, size(output_estimate)...)
-    ind = argmax(output_estimate, dims = 1)
-    # get tensor of maximum elements
-    tensor_of_maximum_elements[ind] .= T1(1)
-    # The count is halved, not floored: a wrongly labelled sample differs from the target in two
-    # entries of its one-hot column and a sample that carries no label differs in one, so half the
-    # difference is the number of errors. The fraction is then read in `float(T)`, for `T` the
-    # element type of the loader's input; reading it in the output's would widen a `Float32`
-    # loader's result, and reading it in an integer `T` fails on a half error.
-    n = size(dl.output, 3)
-    errors = sum(abs.(dl.output - tensor_of_maximum_elements)) / 2
-    F = float(T)
-    F(n - errors) / F(n)
-end
-
-@doc raw"""
-    accuracy(nn, dl)
-
-Compute the accuracy of a neural network classifier.
-
-This is like [`accuracy(::Chain, ::Tuple, ::DataLoader)`](@ref), but for a `NeuralNetwork`.
-"""
-accuracy(nn::NeuralNetwork, dl::DataLoader) = accuracy(nn.model, params(nn), dl)
-
-Base.eltype(::DataLoader{T}) where {T} = T
-
-networkbackend(dl::DataLoader) = networkbackend(dl.input)
-function networkbackend(dl::DataLoader{T, <:NamedTupleOfArrays{T}}) where {T}
-    networkbackend(first(dl.input))
-end
diff --git a/src/data_loader/parametric_data_loader.jl b/src/data_loader/parametric_data_loader.jl
deleted file mode 100644
index 8353d640e..000000000
--- a/src/data_loader/parametric_data_loader.jl
+++ /dev/null
@@ -1,110 +0,0 @@
-@doc raw"""
-    ParametricDataLoader(data, parameters)
-    ParametricDataLoader(ensemble_solution)
-
-Make a loader of trajectories of a system together with the parameters ``\mu`` of the system that
-each trajectory belongs to.
-
-`data` is a `2d × nt × m` array, or a `(q, p)` named tuple of two `d × nt × m` arrays of one
-array type: `m` trajectories of `nt` time steps, with `q` above `p` where they are stacked.
-`parameters` is a vector of `m` named tuples, the ``\mu`` of each trajectory. Every ``\mu`` has the
-keys of the first, in the same order and at every level, and the same sizes of its leaves. The data
-and every ``\mu`` have one element type. An `EnsembleSolution` of a problem with `(q, p)` states
-gives its `q` and `p` and the parameters of its problem.
-
-# The batch
-
-A `Batch(batch_size)` draws, in one epoch, each pair of a start position and a trajectory once, in
-minibatches of `batch_size`; the last minibatch is shorter. A drawn minibatch is the triple
-`(input, output, μs)`: the state at the start position, the state one time step later, and the
-``\mu`` of each sample. [`train!`](@ref) on a `ParametricDataLoader` trains on
-[`ParametricLoss`](@ref)`()` when no `loss` is given, so the network is applied with
-[`apply_parametric`](@ref), one time step per sample. A `Batch` with a sequence length or a
-prediction window other than one is an `ArgumentError`.
-
-A problem without system parameters is trained with a [`DataLoader`](@ref).
-"""
-struct ParametricDataLoader{
-    T, AT <: ArrayOrNamedTuple{T}, PT <: AbstractVector{<:NamedTuple}}
-    input::AT
-    parameters::PT
-    input_time_steps::Int
-    n_params::Int
-
-    function ParametricDataLoader(input::AT,
-            parameters::PT) where {T,
-            AT <: Union{AbstractArray{T, 3},
-                NamedTuple{(:q, :p), Tuple{AT3, AT3}} where {AT3 <: AbstractArray{T, 3}}},
-            PT <: AbstractVector{<:NamedTuple}}
-        input_time_steps, n_params = data_size(input)
-        n_params ≥ 1 && length(parameters) == n_params || throw(ArgumentError(
-            "the data hold $(n_params) trajectories and there are $(length(parameters)) " *
-            "system parameters; a `ParametricDataLoader` takes one for each trajectory, " *
-            "and at least one trajectory."))
-        layout = parameterlayout(first(parameters))
-        for (j, μ) in enumerate(parameters)
-            S = parameter_eltype(μ)
-            S === T || S === Union{} ||
-                throw(ArgumentError(
-                    "the data have the element type $(T) and a system parameter has the " *
-                    "element type $(S); a `ParametricDataLoader` takes one element type."))
-            parameterlayout(μ) == layout || throw(ArgumentError(
-                "system parameter $(j) differs from system parameter 1 in its keys, their " *
-                "order or the sizes of its leaves; a `ParametricDataLoader` takes one shape."))
-        end
-        new{T, AT, PT}(input, parameters, input_time_steps, n_params)
-    end
-end
-
-# The number of time steps and the number of trajectories of the data.
-data_size(input::AbstractArray{<:Any, 3}) = (size(input, 2), size(input, 3))
-
-function data_size(input::NamedTuple{(:q, :p)})
-    size(input.q) == size(input.p) || throw(ArgumentError(
-        "`q` has the size $(size(input.q)) and `p` has the size $(size(input.p))."))
-    data_size(input.q)
-end
-
-function ParametricDataLoader(ensemble_solution::EnsembleSolution{T,
-        T1,
-        Vector{ST}}) where {
-        T, T1, TuT, TT <: TimeSeries{T1},
-        ST <: GeometricSolution{T, T1, TT, NamedTuple{(:t, :q, :p, :q̇, :ṗ), TuT}}}
-    ParametricDataLoader(
-        (q = data_tensor([s.q for s in ensemble_solution.s]),
-            p = data_tensor([s.p for s in ensemble_solution.s])),
-        ensemble_solution.problem.parameters)
-end
-
-function networkbackend(dl::ParametricDataLoader{T, <:AbstractArray}) where {T}
-    networkbackend(dl.input)
-end
-function networkbackend(dl::ParametricDataLoader{T, <:NamedTuple}) where {T}
-    networkbackend(dl.input.q)
-end
-
-# The start positions are those with a target `prediction_window` steps later; the network takes one
-# time step per sample and returns one, so the sequence length and the prediction window are one.
-function (batch::Batch)(dl::ParametricDataLoader)
-    batch.seq_length == 1 || throw(ArgumentError(
-        "a `ParametricDataLoader` draws one time step per sample, and the `Batch` asks for a " *
-        "sequence length of $(batch.seq_length)."))
-    batch.prediction_window == 1 || throw(ArgumentError(
-        "a `ParametricDataLoader` draws a target of one time step per sample, and the `Batch` " *
-        "asks for a prediction window of $(batch.prediction_window)."))
-    n = dl.input_time_steps - batch.prediction_window
-    n ≥ 1 || throw(ArgumentError(
-        "the `ParametricDataLoader` has $(dl.input_time_steps) time steps and the `Batch` asks " *
-        "for a prediction window of $(batch.prediction_window), which leaves no minibatch."))
-    batch_over_two_axes(batch, n, dl.n_params)
-end
-
-# The input at the start position, the target `seq_length` steps later, and the `μ` of the trajectory
-# of each sample, which row 2 of the batch indices names.
-function convert_input_and_batch_indices_to_array(dl::ParametricDataLoader, batch::Batch,
-        batch_indices_tuple::Vector{Tuple{Int, Int}})
-    indices = convert_vector_of_tuples_to_matrix(networkbackend(dl), batch_indices_tuple)
-    (gather(dl.input, indices, 0, batch.seq_length),
-        gather(dl.input, indices, batch.seq_length, batch.prediction_window),
-        [dl.parameters[trajectory] for (_, trajectory) in batch_indices_tuple])
-end
diff --git a/src/data_loader/tensor_assign.jl b/src/data_loader/tensor_assign.jl
deleted file mode 100644
index 435b471b2..000000000
--- a/src/data_loader/tensor_assign.jl
+++ /dev/null
@@ -1,28 +0,0 @@
-@doc raw"""
-    assign_output_estimate(full_output, prediction_window)
-
-Crop the output to get the correct number of output vectors.
-
-The function `assign_output_estimate` is closely related to the [`Transformer`](@ref). 
-It takes the last `prediction_window` columns of the output and uses them for the prediction.
-
-i.e.
-```math
-\mathbb{R}^{N\times{}T}\to\mathbb{R}^{N\times\mathtt{pw}}, 
-\begin{bmatrix} 
-    z^{(1)}_1               & \cdots & z^{(T)}_1 \\ 
-    \cdots                  & \cdots & \cdots    \\ 
-    z^{(1)}_n               & \cdots & z^{(T)}_n
-    \end{bmatrix} \mapsto 
-    \begin{bmatrix} 
-    z^{(T - \mathtt{pw} + 1)}_1 & \cdots      & z^{(T)}_1 \\ 
-    \cdots                  & \cdots      & \cdots \\ 
-    z^{(T - \mathtt{pw} + 1)}_n & \cdots      & z^{(T)}_n\end{bmatrix}     
-``` 
-
-If `prediction_window` is equal to `sequence_length`, then this is not needed.
-"""
-function assign_output_estimate(
-        full_output::AbstractArray{T, 3}, prediction_window::Int) where {T}
-    full_output[:, (end - prediction_window + 1):end, :]
-end
diff --git a/src/data_loader/train.jl b/src/data_loader/train.jl
index b3583eef7..87d6a8c1b 100644
--- a/src/data_loader/train.jl
+++ b/src/data_loader/train.jl
@@ -2,10 +2,10 @@
     train!(nn, opt, dl, batch, n_epochs; loss, show_progress = true)
 
 Train the neural network `nn` with the `GeometricOptimizers.TrainingOptimizer` `opt` on the data of
-the [`DataLoader`](@ref) or [`ParametricDataLoader`](@ref) `dl` for `n_epochs` epochs, and return
+the `DataLoader` or `ParametricDataLoader` `dl` for `n_epochs` epochs, and return
 the loss history.
 
-Each epoch draws the minibatches of `dl` with the [`Batch`](@ref) `batch` and takes one step per
+Each epoch draws the minibatches of `dl` with the `Batch` `batch` and takes one step per
 minibatch: `loss` computes the loss on the minibatch and its gradient, and
 `GeometricOptimizers.optimization_step!` writes the new parameters into `params(nn)`, with the whole
 gradient, a zero leaf included. Entry `i` of the history is the mean of the minibatch losses of epoch
@@ -70,13 +70,19 @@ _as_pullback(p::AbstractPullback) = p
 
 # A `DataLoader` takes the default loss of the architecture and the batch; a `ParametricDataLoader`
 # draws the system parameters of each sample, which `ParametricLoss` hands to the network.
+# An architecture and a batch with no `NetworkLoss` method below have no default, and `default_loss`
+# says so: a method of `NetworkLoss` that raises for every other pair would name only foreign types.
 function default_loss(nn::NeuralNetwork, ::DataLoader, batch::Batch)
-    NetworkLoss(nn.architecture, batch)
+    arch = nn.architecture
+    applicable(NetworkLoss, arch, batch) ||
+        throw(ArgumentError("`train!` has no default loss for a $(nameof(typeof(arch))) with a " *
+                            "$(typeof(batch)); pass the loss as the keyword `loss`."))
+    NetworkLoss(arch, batch)
 end
 default_loss(::NeuralNetwork, ::ParametricDataLoader, ::Batch) = ParametricLoss()
 
-# The default loss of an architecture and a batch. The two-argument methods live here and not beside
-# `NetworkLoss(::AutoEncoder)` in `src/loss/losses.jl`, because `Batch` is included after `src/loss/`.
+# The default loss of an architecture and a batch: what `default_loss` above builds for a network
+# trained on a `DataLoader`.
 function AbstractNeuralNetworks.NetworkLoss(
         ::TransformerIntegrator, batch::Batch{:Transformer})
     TransformerLoss(batch)
@@ -85,7 +91,3 @@ function AbstractNeuralNetworks.NetworkLoss(
         arch::Union{NeuralNetworkIntegrator, AutoEncoder}, ::Batch{:FeedForward})
     NetworkLoss(arch)
 end
-function AbstractNeuralNetworks.NetworkLoss(arch::Architecture, batch::Batch)
-    throw(ArgumentError("`train!` has no default loss for a $(nameof(typeof(arch))) with a " *
-                        "$(typeof(batch)); pass the loss as the keyword `loss`."))
-end
diff --git a/src/loss/losses.jl b/src/loss/losses.jl
index 0820e5d28..4ec69e7d9 100644
--- a/src/loss/losses.jl
+++ b/src/loss/losses.jl
@@ -58,6 +58,12 @@ end
 
 TransformerLoss(seq_length::Int) = TransformerLoss(seq_length, seq_length)
 
+# The loss a `Batch{:Transformer}` implies: `NetworkLoss(arch, batch)` reads the sequence length and
+# the prediction window off the batch instead of taking the two numbers.
+function TransformerLoss(batch::Batch{:Transformer})
+    TransformerLoss(batch.seq_length, batch.prediction_window)
+end
+
 # This crops the output array of the neural network so that it conforms with the output it should be compared to. This is needed for the transformer loss.
 function crop_array_for_transformer_loss(nn_output::AT,
         output::BT) where {T, T2, AT <: AbstractArray{T, 3}, BT <: AbstractArray{T2, 3}}
@@ -242,7 +248,7 @@ system beside its input.
 
 This is `FeedForwardLoss` with the system parameters passed to the network. `μ` is one
 `OptionalParameters` for one sample, or a vector of them, one per sample of a batch, as a
-[`ParametricDataLoader`](@ref) draws it. `train!` on a `ParametricDataLoader` trains on this loss
+`ParametricDataLoader` draws it. `train!` on a `ParametricDataLoader` trains on this loss
 when no `loss` is given.
 """
 struct ParametricLoss <: NetworkLoss end
diff --git a/src/parametric.jl b/src/parametric.jl
index 1866cd4b2..19cce65fb 100644
--- a/src/parametric.jl
+++ b/src/parametric.jl
@@ -46,7 +46,7 @@ of the system.
 `x` is one sample — a vector, or a `(q, p)` named tuple of vectors — with one `μ`, an
 `OptionalParameters`. Or `x` is a batch of `k` samples — a `2d × 1 × k` tensor, or a `(q, p)` named
 tuple of `d × 1 × k` tensors — with a vector of `k` values of `μ`, one per sample, as a
-[`ParametricDataLoader`](@ref) draws them. A batch holds at least one sample. The output has the
+`ParametricDataLoader` draws them. A batch holds at least one sample. The output has the
 shape of the input.
 
 The `Chain` receives the pair `(x, μ)`, so its first layer takes one sample and its `μ`, and a layer
diff --git a/test/Project.toml b/test/Project.toml
index 113884ead..6fba01fa5 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -20,6 +20,7 @@ Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
 NeuralNetworkParameters = "67f4d93a-60e9-472b-8cdd-1ccf6005724a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReducedComplexityModeling = "abb4722c-1d8e-48eb-a3bd-45b3116bfaba"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
diff --git a/test/architectures/parametric_resnet.jl b/test/architectures/parametric_resnet.jl
index f913bc626..507bdd157 100644
--- a/test/architectures/parametric_resnet.jl
+++ b/test/architectures/parametric_resnet.jl
@@ -5,8 +5,8 @@
 
 using GeometricMachineLearning
 using GeometricMachineLearning: ParametricResNet, ParametricResNetLayer, NullParameters,
-                                ParametricDataLoader, params, parameterlength,
-                                apply_parametric
+                                params, parameterlength, apply_parametric
+using ReducedComplexityModeling: ParametricDataLoader
 using Test
 import Random
 
diff --git a/test/data_loader/accuracy.jl b/test/data_loader/accuracy.jl
deleted file mode 100644
index cc4c2446f..000000000
--- a/test/data_loader/accuracy.jl
+++ /dev/null
@@ -1,171 +0,0 @@
-# `accuracy` is documented API -- `docs/src/data_loader/data_loader.md` renders both its methods --
-# and no other function of the package and no script calls it. A documented function with no
-# caller and no test is a function nobody notices breaking, so this file exercises it.
-#
-# The loader comes from `DataLoader(input, output)`. The output's element type is the input's
-# or an integer, so a real input and an integer one-hot target are one call; any other pair of
-# element types is a `MethodError`.
-
-using GeometricMachineLearning
-using Random
-using Test
-using Zygote
-
-include("../helpers/eltypes.jl")
-
-Random.seed!(1234)
-
-const input_dim = 4
-const seq_length = 3
-const n_classes = 3
-const n_params = 6
-
-const input = rand(Float32, input_dim, seq_length, n_params)
-
-# One `DataLoader` per label assignment and element type, over the same input. A label of `0` marks a
-# sample that carries none: its one-hot column stays zero, so the target differs from the prediction
-# in one entry where a wrong label differs in two.
-function classification_loader(classes::AbstractVector{<:Integer}, ::Type{T} = Float32) where {T}
-    output = zeros(Int64, n_classes, 1, n_params)
-    for (k, c) in pairs(classes)
-        iszero(c) && continue
-        output[c, 1, k] = 1
-    end
-    DataLoader(T.(input), output; suppress_info = true)
-end
-
-model = Chain(Dense(input_dim, n_classes, tanh))
-nn = NeuralNetwork(model, Float32)
-ps = GeometricMachineLearning.params(nn)
-
-# The labels are read off the network rather than drawn, so the two assertions below are exact and
-# do not depend on the seed: label every sample with what this network already predicts and the
-# accuracy is 1, rotate every label one class on and it is 0.
-output_estimate = GeometricMachineLearning.assign_output_estimate(model(input, ps), 1)
-predicted = [argmax(output_estimate[:, 1, k]) for k in 1:n_params]
-rotated = [mod1(c + 1, n_classes) for c in predicted]
-
-@testset "accuracy is 1 when every label is the predicted class, and 0 when none is" begin
-    @test GeometricMachineLearning.accuracy(model, ps, classification_loader(predicted)) ==
-          1
-    @test GeometricMachineLearning.accuracy(model, ps, classification_loader(rotated)) == 0
-end
-
-# The fraction is read in the element type of the loader's input, so a `Float32` run of the model
-# gives a `Float32` measure and does not silently widen at the last step. A loader of the model's own
-# element type cannot tell the loader's type from the model's output's; the mixed pair below does.
-@testset "the accuracy has the element type of the loader's input ($T)" for T in REAL_ELTYPES
-    dl = classification_loader(predicted, T)
-    @test eltype(dl) === T
-    @test GeometricMachineLearning.accuracy(model, ps, dl) isa T
-    @test GeometricMachineLearning.accuracy(model, ps, dl) == one(T)
-
-    wrong = classification_loader(rotated, T)
-    @test GeometricMachineLearning.accuracy(model, ps, wrong) isa T
-    @test iszero(GeometricMachineLearning.accuracy(model, ps, wrong))
-end
-
-# The loader's element type decides the result, not the element type of the model or of its output: a
-# `Float64` model on a `Float32` loader gives a `Float32`. The two networks do not agree on every
-# label, so the mixed case labels the loader from the `Float64` network's own predictions.
-@testset "a Float64 model on a Float32 loader still gives a Float32" begin
-    nn64 = NeuralNetwork(model, Float64)
-    ps64 = GeometricMachineLearning.params(nn64)
-    dl32 = classification_loader(predicted, Float32)
-
-    # the output this loader makes is `Float64` here, so the two element types differ
-    estimate64 = GeometricMachineLearning.assign_output_estimate(nn64.model(dl32.input, ps64), 1)
-    @test eltype(estimate64) === Float64
-
-    @test GeometricMachineLearning.accuracy(nn64, dl32) isa Float32
-    @test GeometricMachineLearning.accuracy(nn64.model, ps64, dl32) isa Float32
-
-    # labelled with what the `Float64` network predicts, every label is correct
-    predicted64 = [argmax(estimate64[:, 1, k]) for k in 1:n_params]
-    @test GeometricMachineLearning.accuracy(nn64,
-        classification_loader(predicted64, Float32)) == one(Float32)
-
-    @test GeometricMachineLearning.accuracy(nn64,
-        classification_loader(predicted, Float64)) isa Float64
-end
-
-# A sample that carries no label is half an error: the difference between the target and the
-# prediction is 1 in its one-hot column, where a wrong label gives 2, and the count is halved rather
-# than floored. A floored count drops that sample and reports every label as correct.
-@testset "a sample with no label halves the error count ($T)" for T in REAL_ELTYPES
-    classes = copy(predicted)
-    classes[1] = 0
-    a = GeometricMachineLearning.accuracy(model, ps, classification_loader(classes, T))
-    @test a isa T
-    @test a == (T(n_params) - T(1) / 2) / T(n_params)
-    @test a < one(T)
-end
-
-# An integer input reads the fraction in `float(Int64)`: a half error is not an integer, so the
-# integer type itself cannot hold the result.
-@testset "an integer input gives a Float64, half errors included" begin
-    int_input = rand(0:3, input_dim, 1, n_params)
-    int_model = Chain(Dense(input_dim, n_classes, tanh))
-    int_ps = GeometricMachineLearning.params(NeuralNetwork(int_model, Float64))
-    estimate = GeometricMachineLearning.assign_output_estimate(int_model(int_input, int_ps), 1)
-    labels = [argmax(estimate[:, 1, k]) for k in 1:n_params]
-    function int_loader(classes)
-        output = zeros(Int64, n_classes, 1, n_params)
-        for (k, c) in pairs(classes)
-            iszero(c) || (output[c, 1, k] = 1)
-        end
-        DataLoader(int_input, output; suppress_info = true)
-    end
-    @test eltype(int_loader(labels)) === Int64
-
-    a = GeometricMachineLearning.accuracy(int_model, int_ps, int_loader(labels))
-    @test a isa Float64
-    @test a == 1
-
-    unlabelled = copy(labels)
-    unlabelled[1] = 0
-    a = GeometricMachineLearning.accuracy(int_model, int_ps, int_loader(unlabelled))
-    @test a isa Float64
-    @test a == (n_params - 1 / 2) / n_params
-end
-
-# A `Bool` one-hot target is an integer target too: the count is halved in floating point, not in the
-# label type, so a half error does not have to fit in a `Bool`.
-@testset "a Bool target gives the accuracy of the Int64 one" begin
-    classes = copy(predicted)
-    classes[1] = 0
-    int_dl = classification_loader(classes)
-    bool_dl = DataLoader(int_dl.input, Bool.(int_dl.output); suppress_info = true)
-    a = GeometricMachineLearning.accuracy(model, ps, bool_dl)
-    @test a isa Float32
-    @test a == (Float32(n_params) - Float32(1) / 2) / Float32(n_params)
-    @test a == GeometricMachineLearning.accuracy(model, ps, int_dl)
-end
-
-@testset "the NeuralNetwork method agrees with the Chain method" begin
-    for classes in (predicted, rotated)
-        dl = classification_loader(classes)
-        @test GeometricMachineLearning.accuracy(nn, dl) ==
-              GeometricMachineLearning.accuracy(model, ps, dl)
-    end
-end
-
-# The crop of `assign_output_estimate` is the last `prediction_window` columns, and its gradient is
-# Zygote's `getindex` adjoint — zero on the cropped columns, and a cotangent of the primal's element
-# type whatever the seed's.
-@testset "assign_output_estimate crops the last columns ($T)" for T in REAL_ELTYPES
-    x = T.(reshape(1:(3 * 5 * 4), 3, 5, 4))
-    @test GeometricMachineLearning.assign_output_estimate(x, 2) == x[:, 4:5, :]
-    @test GeometricMachineLearning.assign_output_estimate(x, 5) == x
-    @test eltype(GeometricMachineLearning.assign_output_estimate(x, 2)) === T
-end
-
-@testset "the gradient of assign_output_estimate is Zygote's ($T)" for T in REAL_ELTYPES
-    x = T.(reshape(1:(3 * 5 * 4), 3, 5, 4))
-    crop(z) = GeometricMachineLearning.assign_output_estimate(z, 2)
-    gradient = only(Zygote.gradient(z -> sum(abs2, crop(z)), x))
-    @test all(iszero, gradient[:, 1:3, :])
-    @test gradient[:, 4:5, :] == 2 .* x[:, 4:5, :]
-    _, back = Zygote.pullback(crop, x)
-    @test eltype(only(back(ones(Float64, 3, 2, 4)))) === T
-end
diff --git a/test/data_loader/autoencoder_keyword.jl b/test/data_loader/autoencoder_keyword.jl
deleted file mode 100644
index 4e9b6dc51..000000000
--- a/test/data_loader/autoencoder_keyword.jl
+++ /dev/null
@@ -1,42 +0,0 @@
-# Data-set constructors accept `autoencoder::Bool`. A value outside `Bool`, such as `nothing`,
-# is rejected at the keyword with a `TypeError` rather than interpreted as a loader kind.
-
-using GeometricMachineLearning
-using Test
-import Random
-
-Random.seed!(1234)
-
-@testset "DataLoader(::AbstractArray{<:Number, 3})" begin
-    data = rand(Float32, 2, 3, 4)
-    @test_throws TypeError DataLoader(data; autoencoder = nothing, suppress_info = true)
-end
-
-@testset "DataLoader(::AbstractMatrix)" begin
-    data = rand(Float32, 2, 3)
-    @test_throws TypeError DataLoader(data; autoencoder = nothing, suppress_info = true)
-end
-
-@testset "DataLoader(::NamedTuple{(:q, :p)}) of matrices" begin
-    data = (q = rand(Float32, 2, 3), p = rand(Float32, 2, 3))
-    @test_throws TypeError DataLoader(data; autoencoder = nothing, suppress_info = true)
-    # the vector pair reaches the same branch through the matrices above
-    @test_throws TypeError DataLoader((q = rand(Float32, 3), p = rand(Float32, 3));
-        autoencoder = nothing, suppress_info = true)
-end
-
-@testset "DataLoader(::NamedTuple{(:q, :p)}) of tensors" begin
-    data = (q = rand(Float32, 2, 3, 4), p = rand(Float32, 2, 3, 4))
-    @test_throws TypeError DataLoader(data; autoencoder = nothing, suppress_info = true)
-end
-
-# `DataLoader(dl, backend)` accepts `autoencoder::Union{Nothing, Bool}`: `nothing` inherits the
-# loader's kind, a Boolean selects it, and a value outside the union raises a `TypeError`.
-@testset "DataLoader(dl, backend): nothing means inherit, and nothing falls through" begin
-    dl = DataLoader(rand(Float32, 2, 3, 4); autoencoder = true, suppress_info = true)
-
-    inherited = DataLoader(dl, GeometricMachineLearning.CPU())
-    @test typeof(inherited).parameters[4] === typeof(dl).parameters[4] === :RegularData
-
-    @test_throws TypeError DataLoader(dl, GeometricMachineLearning.CPU(); autoencoder = "x")
-end
diff --git a/test/data_loader/batch.jl b/test/data_loader/batch.jl
index 3d55d39a1..ddf2d2df1 100644
--- a/test/data_loader/batch.jl
+++ b/test/data_loader/batch.jl
@@ -1,9 +1,11 @@
-# The batching rules of `src/data_loader/batch.jl`. The data are ten samples
-# `x = T.(reshape(1:20, 2, 10))` with `y = 3 .* x`, so that `output == 3 .* input` identifies the
-# pair of one sample and not the pair of another, and a drawn index pair names the sample it draws.
+# The `Batch` method of `integrate` exists to carry one number: it is the three-argument form a
+# training script writes, and it has to give the trajectory the keyword form gives. The window is
+# chosen so that the two disagree when it is dropped: a standard transformer of system dimension
+# 5 on a four-column input takes two predictions of two columns for `n_points = 8`, where the
+# default window takes one of four columns. Both forms carry the input's element type, so a
+# trajectory copied to the host or plotted keeps its precision.
 
 using GeometricMachineLearning
-using GeometricMachineLearning: convert_input_and_batch_indices_to_array, number_of_batches
 using Test
 import Random
 
@@ -11,211 +13,10 @@ include("../helpers/eltypes.jl")
 
 Random.seed!(1234)
 
-# The kind of a loader is its fourth type parameter.
-kind(dl::DataLoader) = typeof(dl).parameters[4]
-
-samples(T) = T.(reshape(1:20, 2, 10))
-
-# What one minibatch holds, for every minibatch of an epoch.
-gathered(dl, batch) = map(batch(dl)) do indices
-    convert_input_and_batch_indices_to_array(dl, batch, indices)
-end
-
-@testset "a supervised loader samples every pair, the last included ($T)" for T in REAL_ELTYPES
-    x = samples(T)
-    y = 3 .* x
-    loaders = ("input and output" => DataLoader(x, y; suppress_info = true),
-        "a view output" => DataLoader(x, view(y, :, :); suppress_info = true))
-    for (name, dl) in loaders
-        @test kind(dl) === :RegularData
-        @test (dl.input_dim, dl.input_time_steps, dl.n_params) == (2, 10, 1)
-        @test (dl.output_dim, dl.output_time_steps) == (2, 10)
-
-        batch = Batch(3)
-        minibatches = batch(dl)
-        # every sample is drawn, once, the tenth included
-        @test sort(vcat(minibatches...)) == [(t, 1) for t in 1:10]
-        @test length(minibatches) == 4
-        @test number_of_batches(dl, batch) == 4
-
-        data = gathered(dl, batch)
-        inputs = first.(data)
-        @test sort(size.(inputs, 3)) == [1, 3, 3, 3]
-        @test all(size(input)[1:2] == (2, 1) for input in inputs)
-        for (input, output) in data
-            @test eltype(input) === T
-            @test eltype(output) === T
-            @test size(output) == size(input)
-            @test output == 3 .* input
-        end
-        # the gathered inputs are the columns of `x`, in some order
-        @test sort(eachcol(hcat((input[:, 1, k]
-        for input in inputs for k in axes(input, 3))...))) ==
-              sort(eachcol(Matrix(x)))
-    end
-end
-
-# The `(q, p)` loader of the same ten samples: `autoencoder = true` reads the second axis of a
-# matrix pair as the parameter index, so its one time step carries all ten samples.
-@testset "a (q, p) loader of samples draws every pair ($T)" for T in REAL_ELTYPES
-    x = samples(T)
-    dl = DataLoader((q = x, p = 3 .* x); autoencoder = true, suppress_info = true)
-    @test kind(dl) === :RegularData
-    @test (dl.input_dim, dl.input_time_steps, dl.n_params) == (4, 1, 10)
-
-    batch = Batch(3)
-    drawn = batch(dl)
-    @test sort(vcat(drawn...)) == [(1, p) for p in 1:10]
-    minibatches = map(
-        indices -> convert_input_and_batch_indices_to_array(dl, batch, indices),
-        drawn)
-    @test length(minibatches) == 4
-    @test sort([size(mb.q, 3) for mb in minibatches]) == [1, 3, 3, 3]
-    for (indices, mb) in zip(drawn, minibatches)
-        @test mb isa NamedTuple{(:q, :p)}
-        @test eltype(mb.q) === T
-        @test eltype(mb.p) === T
-        @test mb.p == 3 .* mb.q
-        # each slot holds the sample its own drawn index names, so a `q` and a `p` read from the
-        # same other column shows here
-        @test all(mb.q[:, 1, k] == x[:, indices[k][2]] for k in axes(mb.q, 3))
-    end
-    # the gathered columns are the columns of `x`, in some order
-    @test sort(eachcol(hcat((mb.q[:, 1, k]
-    for mb in minibatches for k in axes(mb.q, 3))...))) ==
-          sort(eachcol(Matrix(x)))
-end
-
-# One sample is a loader the training loop can draw from.
-@testset "one sample is one minibatch ($T)" for T in REAL_ELTYPES
-    x = T.(reshape(1:2, 2, 1))
-    dl = DataLoader(x, x; suppress_info = true)
-    batch = Batch(1)
-    minibatches = batch(dl)
-    @test length(minibatches) == 1
-    @test only(minibatches) == [(1, 1)]
-    input, output = only(gathered(dl, batch))
-    @test size(input) == (2, 1, 1)
-    @test eltype(input) === T
-    @test eltype(output) === T
-    @test output == input
-end
-
-# One count, `cld(n · n_params, batch_size)` with `n` the number of start positions, and an
-# `ArgumentError` where no window fits or the parameter axis is empty.
-@testset "one count, and a window that fits no start position" begin
-    trajectory = DataLoader(rand(Float64, 2, 10); suppress_info = true)
-    @test kind(trajectory) === :TimeSeries
-    @test length(Batch(3, 9, 1)(trajectory)) == 1
-    @test_throws ArgumentError Batch(3, 10, 1)(trajectory)
-
-    x = samples(Float64)
-    dl = DataLoader(x, 3 .* x; suppress_info = true)
-    @test length(Batch(3, 10, 10)(dl)) == 1
-    @test_throws ArgumentError Batch(3, 11, 1)(dl)
-
-    # the message names the time steps, the sequence length and the prediction window: three
-    # different numbers here, so that each one's presence is a separate reading
-    raised = try
-        Batch(3, 11, 5)(trajectory)
-    catch e
-        e
-    end
-    @test raised isa ArgumentError
-    for word in ("time steps", "sequence length", "prediction window")
-        @test occursin(word, raised.msg)
-    end
-    @test occursin("$(trajectory.input_time_steps)", raised.msg)
-    @test occursin("11", raised.msg)
-    @test occursin("5", raised.msg)
-
-    # a loader whose parameter axis is empty draws no minibatch either
-    @test_throws ArgumentError Batch(1)(DataLoader(zeros(Float64, 2, 0);
-        autoencoder = true, suppress_info = true))
-    @test_throws ArgumentError Batch(1)(DataLoader(zeros(Float64, 2, 3, 0); suppress_info = true))
-
-    # where the output's time steps bind the window, the message names them
-    with_output = DataLoader(rand(Float64, 2, 10, 3), rand(Float64, 2, 1, 3);
-        suppress_info = true)
-    raised = try
-        Batch(2, 1, 2)(with_output)
-    catch e
-        e
-    end
-    @test raised isa ArgumentError
-    @test occursin("10 time steps and 1 in its output", raised.msg)
-end
-
-# The classifier's loader, and `Batch(batch_size, dl)` — one draw per sample with the whole
-# sequence and the whole output, which is what GMLDatasets calls.
-@testset "Batch(batch_size, dl) draws every sample with the whole sequence ($T)" for T in REAL_ELTYPES
-    n_classes, n_patches, n_images = 10, 4, 10
-    input = rand(T, 9, n_patches, n_images)
-    output = zeros(Int, n_classes, 1, n_images)
-    for (k, c) in pairs(rand(1:n_classes, n_images))
-        output[c, 1, k] = 1
-    end
-    dl = DataLoader(input, output; suppress_info = true)
-    @test kind(dl) === :RegularData
-    @test (dl.input_dim, dl.input_time_steps, dl.n_params) == (9, 4, 10)
-    @test (dl.output_dim, dl.output_time_steps) == (10, 1)
-
-    batch = Batch(2, dl)
-    @test batch == Batch(2, 4, 1)
-    minibatches = batch(dl)
-    @test length(minibatches) == 5
-    # `Batch(2)` would pair the first patch with the one output column, so it raises
-    @test_throws ArgumentError Batch(2)(dl)
-    @test_throws "Batch(batch_size, dl)" number_of_batches(dl, Batch(2))
-    for indices in minibatches
-        input_batch, output_batch = convert_input_and_batch_indices_to_array(dl, batch, indices)
-        @test size(input_batch) == (9, 4, 2)
-        @test size(output_batch) == (10, 1, 2)
-        @test eltype(input_batch) === T
-        @test eltype(output_batch) === Int
-        @test output_batch == output[:, :, [i[2] for i in indices]]
-    end
-end
-
-# The target is the output when there is one, whatever kind the loader is typed as. A
-# `DataLoader{T, A3, Array{Int, 3}, :TimeSeries}` is what GMLDatasets and two tests build by naming
-# the type, and it must batch like the `:RegularData` loader of the same data.
-@testset "a loader with an output is samples whatever its kind ($T)" for T in REAL_ELTYPES
-    input = rand(T, 9, 4, 10)
-    output = zeros(Int, 10, 1, 10)
-    for (k, c) in pairs(rand(1:10, 10))
-        output[c, 1, k] = 1
-    end
-    dl = DataLoader(input, output; suppress_info = true)
-    dl_time_series = DataLoader{T, typeof(input), typeof(output), :TimeSeries}(
-        input, output, 9, 4, 10, 10, 1)
-
-    for batch in (Batch(2, dl), Batch(3, 1, 1))
-        Random.seed!(11)
-        asked = batch(dl_time_series)
-        Random.seed!(11)
-        reference = batch(dl)
-        @test asked == reference
-        for (i, j) in zip(asked, reference)
-            here = convert_input_and_batch_indices_to_array(dl_time_series, batch, i)
-            @test here == convert_input_and_batch_indices_to_array(dl, batch, j)
-            @test eltype(here[1]) === T
-            @test eltype(here[2]) === Int
-        end
-    end
-end
-
-# The `Batch` method of `integrate` exists to carry one number: it is the three-argument form a
-# training script writes, and it has to give the trajectory the keyword form gives. The window is
-# chosen so that the two disagree when it is dropped: a standard transformer of system dimension
-# 5 on a four-column input takes two predictions of two columns for `n_points = 8`, where the
-# default window takes one of four columns. Both forms carry the input's element type, so a
-# trajectory copied to the host or plotted keeps its precision.
 @testset "the Batch method of integrate forwards its prediction window ($T)" for T in REAL_ELTYPES
     nn = NeuralNetwork(StandardTransformerIntegrator(5), CPU(), T)
     ics = rand(T, 5, 4)
 
-    @test Batch(3, 5, 2) == Batch{:Transformer}(3, 5, 2)
     @test eltype(integrate(nn, ics, Batch(3, 5, 2); n_points = 5)) === T
     @test eltype(integrate(nn, ics; n_points = 5, prediction_window = 2)) === T
     @test integrate(nn, ics, Batch(3, 5, 2); n_points = 5) ==
@@ -225,3 +26,14 @@ end
     @test integrate(nn, ics, Batch(3, 5, 2); n_points = 8) ≠
           integrate(nn, ics; n_points = 8)
 end
+
+# The three-argument method is a method of `GeometricBase.integrate` on `ReducedComplexityModeling`'s
+# `Batch`, so its `nn` names this package's `TransformerIntegrator`; a signature without that
+# annotation is piracy. `Aqua`'s check does not see this one: it treats a symbol type parameter as a
+# type of the package's own, so no signature written on `Batch{:Transformer}` reads as foreign. That
+# is why the narrowing has a gate here and the piracies Aqua does see have one in
+# `test/quality/aqua.jl`.
+@testset "the Batch method of integrate is narrowed to a TransformerIntegrator" begin
+    @test only(methods(integrate, (NeuralNetwork, Any, Batch{:Transformer}))).sig <:
+          Tuple{typeof(integrate), NeuralNetwork{<:TransformerIntegrator}, Any, Any}
+end
diff --git a/test/data_loader/batch_functor.jl b/test/data_loader/batch_functor.jl
deleted file mode 100644
index 704e9301e..000000000
--- a/test/data_loader/batch_functor.jl
+++ /dev/null
@@ -1,59 +0,0 @@
-using Test
-using GeometricMachineLearning
-import Random
-
-Random.seed!(123)
-
-const data_set_raw = 1:10
-const data_set = reshape(data_set_raw, 1, 10)
-
-const dl = DataLoader(data_set; autoencoder = true)
-
-function test_batching_over_one_axis(batch_size::Int = 3)
-    batch = Batch(batch_size)
-    tuple_indices = batch(dl)
-    @test length(tuple_indices) == Int(ceil(10 / batch_size))
-    tuple_indices_sorted = sort(vcat(tuple_indices...))
-    indices_sorted = [tuple_index[2] for tuple_index in tuple_indices_sorted]
-    @test indices_sorted == data_set_raw
-end
-
-test_batching_over_one_axis()
-test_batching_over_one_axis(5)
-
-const data_set_raw₂ = hcat(1:10, 21:30)
-const data_set₂ = reshape(data_set_raw₂, 1, 10, 2)
-
-const dl₂ = DataLoader(data_set₂)
-
-function convert_vector_of_tuples_to_array(indices::Vector{<:Tuple})
-    elem_in_tuple = length(indices[1])
-    index_array = zeros(Int, length(indices), elem_in_tuple)
-    for (i, elem) in zip(axes(indices, 1), indices)
-        index_array[i, :] .= elem
-    end
-    index_array
-end
-
-sort_according_to_first_column(A::Matrix) = A[sortperm(A[:, 1]), :]
-
-function test_batching_over_two_axis_with_seq_length(
-        batch_size::Int = 3, seq_length::Int = 2, prediction_window::Int = 1)
-    batch = Batch(batch_size, seq_length, prediction_window)
-    tuple_indices₂ = batch(dl₂)
-    @test length(tuple_indices₂) ==
-          Int(ceil(2 * (10 - (seq_length - 1) - batch.prediction_window) / batch_size)) # have to multiply with two because we deal with tuples
-    num_elems = 10 - (batch.seq_length - 1) - batch.prediction_window
-    indices₁ = hcat(1:num_elems, 1 * ones(Int, num_elems))
-    indices₂ = hcat(1:num_elems, 2 * ones(Int, num_elems))
-    true_indices₁ = sort_according_to_first_column(vcat(indices₁, indices₂))
-    true_indices₂ = sort_according_to_first_column(vcat(indices₂, indices₁))
-    # convert indices to a format such that we can compare them to the true indices
-    indices_for_comparison = convert_vector_of_tuples_to_array(vcat(tuple_indices₂...))
-    # sort the indices
-    indices_for_comparison = sort_according_to_first_column(indices_for_comparison)
-    @test indices_for_comparison == true_indices₁ || indices_for_comparison == true_indices₂
-end
-
-test_batching_over_two_axis_with_seq_length()
-test_batching_over_two_axis_with_seq_length(3, 2, 2)
diff --git a/test/data_loader/batch_index_set.jl b/test/data_loader/batch_index_set.jl
deleted file mode 100644
index 24335f388..000000000
--- a/test/data_loader/batch_index_set.jl
+++ /dev/null
@@ -1,30 +0,0 @@
-# `batch_over_two_axes` builds the minibatch index set. These two assertions hold it to the two
-# properties that keep it cheap: its return type is concrete, and four times the minibatches cost
-# about four times the memory. A tuple grown by splatting has neither — it infers as
-# `Tuple{Vararg{Vector{Tuple{Int, Int}}}}` and costs the square of the minibatch count. See
-# `CHANGELOG.md`.
-
-using GeometricMachineLearning
-using GeometricMachineLearning: batch_over_two_axes
-using Test
-import Random
-
-Random.seed!(123)
-
-@test isconcretetype(Base.infer_return_type(
-    batch_over_two_axes, Tuple{Batch, Int, Int}))
-
-function index_set_allocations(n_columns::Int, batch_size::Int = 8)
-    dl = DataLoader(rand(Float32, 4, n_columns); autoencoder = false, suppress_info = true)
-    batch = Batch(batch_size)
-    batch(dl)
-    @allocated batch(dl)
-end
-
-const allocations_250_batches = index_set_allocations(2000)
-const allocations_1000_batches = index_set_allocations(8000)
-
-# Four times the minibatches. Linear growth puts the ratio at four; a quadratic index set puts it
-# at 12.25 here. The bound is wide enough for the fixed per-call cost and tight enough to fail on
-# quadratic growth.
-@test allocations_1000_batches / allocations_250_batches < 6
diff --git a/test/data_loader/data_loader.jl b/test/data_loader/data_loader.jl
index 7d0810a8b..c7bffbb0d 100644
--- a/test/data_loader/data_loader.jl
+++ b/test/data_loader/data_loader.jl
@@ -1,211 +1,18 @@
-# The constructors of `src/data_loader/data_loader.jl`: the `autoencoder` default and the keyword
-# set of every method, the fifth type parameter of the struct (#211), the copy loops out of a
-# solution and an ensemble, the element-type conversion of `DataLoader(dl, T)`, and the
-# `(input, output)` method that lets a classifier's data be built without naming the parametric type.
+# The reproduction script `scripts/reproduction/symplectic_autoencoders/training.jl` has a `solve!`
+# of a PSD architecture on a snapshot matrix, and `solve!` takes a `:RegularData` loader. That is
+# why its `DataLoader(data)` passes `autoencoder = true`, and why a matrix that does not is a
+# `MethodError` rather than a training run on input-shifted targets.
 
 using GeometricMachineLearning
 using GeometricMachineLearning: solve!
 using Test
 import Random
 
-using GeometricIntegrators: integrate, ImplicitMidpoint
-using GeometricProblems.HarmonicOscillator: hodeensemble, hodeproblem
-using GeometricProblems.RigidBody: odeensemble, odeproblem, default_parameters
-
-include("../helpers/eltypes.jl")
-
 Random.seed!(1234)
 
-# The kind of a loader is its fourth type parameter.
-kind(dl::DataLoader) = typeof(dl).parameters[4]
-
-@testset "every constructor defaults autoencoder = false ($T)" for T in REAL_ELTYPES
-    matrix = rand(T, 3, 10)
-    tensor = rand(T, 3, 10, 2)
-    @test kind(DataLoader(matrix; suppress_info = true)) === :TimeSeries
-    @test kind(DataLoader(rand(T, 10); suppress_info = true)) === :TimeSeries
-    @test kind(DataLoader(tensor; suppress_info = true)) === :TimeSeries
-    @test kind(DataLoader((q = matrix, p = matrix); suppress_info = true)) === :TimeSeries
-    # A literal keyword fixes the loader's kind, so the tensor constructors infer a concrete type.
-    @test (@inferred (() -> DataLoader(tensor; suppress_info = true))()) isa DataLoader
-    @test (@inferred (() -> DataLoader((q = tensor, p = tensor);
-        suppress_info = true))()) isa DataLoader
-    @test (@inferred (() -> DataLoader(tensor; autoencoder = true,
-        suppress_info = true))()) isa DataLoader
-    @test (@inferred (() -> DataLoader((q = tensor, p = tensor); autoencoder = true,
-        suppress_info = true))()) isa DataLoader
-    # ... and so do the matrix and vector constructors, which forward to them
-    vector = rand(T, 10)
-    @test (@inferred (() -> DataLoader(matrix; suppress_info = true))()) isa DataLoader
-    @test (@inferred (() -> DataLoader(vector; autoencoder = true,
-        suppress_info = true))()) isa DataLoader
-    @test (@inferred (() -> DataLoader((q = matrix, p = matrix);
-        suppress_info = true))()) isa DataLoader
-    @test (@inferred (() -> DataLoader((q = vector, p = vector); autoencoder = true,
-        suppress_info = true))()) isa DataLoader
-    @test kind(DataLoader((q = tensor, p = tensor); suppress_info = true)) === :TimeSeries
-    @test kind(DataLoader((q = rand(T, 10), p = rand(T, 10)); suppress_info = true)) ===
-          :TimeSeries
-    # ... and `autoencoder = true` is what makes a matrix a set of samples
-    @test kind(DataLoader(matrix; autoencoder = true, suppress_info = true)) ===
-          :RegularData
-    @test kind(DataLoader((q = rand(T, 10), p = rand(T, 10));
-        autoencoder = true, suppress_info = true)) === :RegularData
-    # a loader with an output is a set of samples, whatever the two element types are
-    @test kind(DataLoader(matrix, matrix; suppress_info = true)) === :RegularData
-end
-
-# The matrix message names the second axis by the kind actually built.
-@testset "the matrix message names the axis by the kind" begin
-    @test_logs (:info, r"time steps") DataLoader(rand(Float32, 3, 10))
-    @test_logs (:info, r"number of parameters") DataLoader(
-        rand(Float32, 3, 10); autoencoder = true)
-end
-
-# `suppress_info` is a keyword in every method of a data set, a solution and an ensemble.
-@testset "suppress_info is a keyword everywhere" begin
-    solution = integrate(
-        hodeproblem(; timespan = (0.0, 1.0), timestep = 0.25), ImplicitMidpoint())
-    @test DataLoader(solution; suppress_info = true) isa DataLoader
-    @test_throws MethodError DataLoader(solution, true)
-
-    ensemble = integrate(
-        odeensemble([(q = [0.0, 0.2, 1.0],), (q = [0.1, 0.3, 1.0],)];
-            timespan = (0.0, 0.4), timestep = 0.2, parameters = default_parameters()),
-        ImplicitMidpoint())
-    @test DataLoader(ensemble; suppress_info = true) isa DataLoader
-    @test_throws MethodError DataLoader(ensemble, true)
-end
-
-# The copy loops of `data_tensors_from_geometric_solution` build one 1-based tensor per coordinate.
-# The rigid-body `odeensemble` is the case whose `dataser` carries `(t, q, q̇)` and no `p`, which is
-# the case its own method is for.
-@testset "a solution and an ensemble become 1-based tensors" begin
-    ics = [(q = [sin(v), 0.0, cos(v)],) for v in 0.1:0.5:2.0]
-    ensemble = integrate(
-        odeensemble(ics; timespan = (0.0, 1.0), timestep = 0.2,
-            parameters = default_parameters()), ImplicitMidpoint())
-    dl = DataLoader(ensemble; suppress_info = true)
-    @test dl.input isa Array{Float64, 3}
-    @test kind(dl) === :TimeSeries
-    @test size(dl.input) == (3, length(ensemble.t), length(ensemble.s))
-    for (i, solution) in pairs(ensemble.s), t in 0:(length(solution.t) - 1)
-
-        @test dl.input[:, t + 1, i] == solution.q[t]
-    end
-
-    # one solution of the same problem
-    solution = integrate(
-        odeproblem(; timespan = (0.0, 1.0), timestep = 0.2), ImplicitMidpoint())
-    dl_one = DataLoader(solution; suppress_info = true)
-    @test dl_one.input isa Array{Float64, 3}
-    @test size(dl_one.input) == (3, length(solution.t), 1)
-    for t in 0:(length(solution.t) - 1)
-        @test dl_one.input[:, t + 1, 1] == solution.q[t]
-    end
-
-    # and a `(q, p)` ensemble: the two coordinates are not swapped
-    ensemble_qp = integrate(
-        hodeensemble(; timespan = (0.0, 0.5), timestep = 0.25), ImplicitMidpoint())
-    dl_qp = DataLoader(ensemble_qp; suppress_info = true)
-    @test dl_qp.input isa NamedTuple{(:q, :p)}
-    @test size(dl_qp.input.q) == (1, length(ensemble_qp.t), length(ensemble_qp.s))
-    for (i, sol) in pairs(ensemble_qp.s), t in 0:(length(sol.t) - 1)
-
-        @test dl_qp.input.q[:, t + 1, i] == sol.q[t]
-        @test dl_qp.input.p[:, t + 1, i] == sol.p[t]
-    end
-
-    # and one `(q, p)` solution: the two coordinates are not swapped
-    solution_qp = integrate(
-        hodeproblem(; timespan = (0.0, 1.0), timestep = 0.25), ImplicitMidpoint())
-    dl_qp_one = DataLoader(solution_qp; suppress_info = true)
-    @test dl_qp_one.input isa NamedTuple{(:q, :p)}
-    @test size(dl_qp_one.input.q) == (1, length(solution_qp.t), 1)
-    for t in 0:(length(solution_qp.t) - 1)
-        @test dl_qp_one.input.q[:, t + 1, 1] == solution_qp.q[t]
-        @test dl_qp_one.input.p[:, t + 1, 1] == solution_qp.p[t]
-    end
-end
-
-# The reproduction script `scripts/reproduction/symplectic_autoencoders/training.jl` has a `solve!`
-# of a PSD architecture on a snapshot matrix, and `solve!` takes a `:RegularData` loader. That is
-# why its `DataLoader(data)` passes `autoencoder = true`, and why a matrix that does not is a
-# `MethodError` rather than a training run on input-shifted targets.
 @testset "solve! of a PSD architecture needs a samples loader" begin
     data = rand(Float32, 6, 20)
     nn = NeuralNetwork(PSDArch(6, 2), CPU(), Float32)
     @test isfinite(solve!(nn, DataLoader(data; autoencoder = true, suppress_info = true)))
     @test_throws MethodError solve!(nn, DataLoader(data; suppress_info = true))
 end
-
-# #211: the two fields that are `nothing` for an unsupervised loader are concrete.
-@testset "#211: the output fields are concrete in every loader" begin
-    loaders = (DataLoader(rand(Float32, 2, 3, 4); suppress_info = true),
-        DataLoader(rand(Float32, 2, 3, 4), rand(Float32, 2, 3, 4); suppress_info = true),
-        DataLoader(rand(Float32, 2, 3, 4); autoencoder = true, suppress_info = true),
-        DataLoader(rand(Float32, 2, 3); autoencoder = true, suppress_info = true))
-    for dl in loaders
-        @test isconcretetype(fieldtype(typeof(dl), :output_dim))
-        @test isconcretetype(fieldtype(typeof(dl), :output_time_steps))
-        @test (fieldtype(typeof(dl), :output_dim) === Nothing) == (dl.output === nothing)
-    end
-end
-
-# The four-parameter constructor fills the fifth type parameter in, and it takes the counts as any
-# integer type: `GMLDatasets` builds a loader by naming the first four and handing it the `size`s.
-@testset "the four-parameter constructor converts its counts" begin
-    input = rand(Float32, 2, 3, 4)
-    output = zeros(Int, 3, 1, 4)
-    dl = DataLoader{Float32, typeof(input), typeof(output), :TimeSeries}(input, output,
-        Int32(2), Int32(3), Int32(4), Int32(3), Int32(1))
-    @test dl isa DataLoader{Float32, Array{Float32, 3}, Array{Int, 3}, :TimeSeries, Int}
-    @test (dl.input_dim, dl.input_time_steps, dl.n_params) == (2, 3, 4)
-    @test (dl.output_dim, dl.output_time_steps) == (3, 1)
-    @test dl.output === output
-
-    no_output = DataLoader{Float32, typeof(input), Nothing, :RegularData}(input, nothing,
-        Int8(2), Int8(3), Int8(4), nothing, nothing)
-    @test no_output isa
-          DataLoader{Float32, Array{Float32, 3}, Nothing, :RegularData, Nothing}
-    @test (no_output.output_dim, no_output.output_time_steps) == (nothing, nothing)
-end
-
-# The conversion to another element type keeps the kind and the counts, and the other backend
-# branch is the identity.
-@testset "DataLoader(dl, T) converts and keeps the kind ($T)" for T in REAL_ELTYPES
-    other = T === Float32 ? Float64 : Float32
-    dl = DataLoader(rand(other, 2, 3, 4); suppress_info = true)
-    converted = DataLoader(dl, T)
-    @test eltype(converted) === T
-    @test converted.input isa Array{T, 3}
-    @test converted.input == T.(dl.input)
-    @test kind(converted) === kind(dl)
-    @test (converted.input_dim, converted.input_time_steps, converted.n_params) ==
-          (dl.input_dim, dl.input_time_steps, dl.n_params)
-    @test converted.output === nothing
-    # the same backend and the same element type is the same tensor
-    @test DataLoader(dl, other).input === dl.input
-end
-
-# A classifier's real input and integer one-hot output, without naming the parametric type.
-@testset "a loader of an input and an integer output" begin
-    input = rand(Float32, 4, 3, 6)
-    output = zeros(Int64, 3, 1, 6)
-    dl = DataLoader(input, output; suppress_info = true)
-    @test kind(dl) === :RegularData
-    @test (dl.input_dim, dl.input_time_steps, dl.n_params) == (4, 3, 6)
-    @test (dl.output_dim, dl.output_time_steps) == (3, 1)
-    @test eltype(dl) === Float32
-    @test dl.output === output
-    # the same element type in both arguments still reaches a method, and so does an `Int`/`Int` pair
-    @test DataLoader(input, rand(Float32, 3, 1, 6); suppress_info = true) isa DataLoader
-    @test DataLoader(zeros(Int, 2, 3, 4), zeros(Int64, 2, 3, 4); suppress_info = true) isa
-          DataLoader{Int, Array{Int, 3}, Array{Int64, 3}, :RegularData}
-    # two element types of which neither is an integer stay a `MethodError`
-    @test_throws MethodError DataLoader(
-        rand(Float32, 2, 3, 4), rand(Float64, 2, 3, 4); suppress_info = true)
-    # the third axis is still the parameter index of both
-    @test_throws AssertionError DataLoader(
-        rand(Float32, 2, 3, 4), zeros(Int, 3, 1, 5); suppress_info = true)
-end
diff --git a/test/data_loader/draw_batch_for_tensor_test.jl b/test/data_loader/draw_batch_for_tensor_test.jl
deleted file mode 100644
index 46e04fb8b..000000000
--- a/test/data_loader/draw_batch_for_tensor_test.jl
+++ /dev/null
@@ -1,26 +0,0 @@
-using GeometricMachineLearning, Test
-using GeometricMachineLearning: convert_input_and_batch_indices_to_array
-import Random
-
-Random.seed!(1234)
-
-function test_data_loader_for_qp_tensor(
-        system_dim2::Int, input_time_steps::Int, n_params::Int, batch_size::Int,
-        seq_length::Int, prediction_window::Union{Int, Nothing} = nothing)
-    _prediction_window = isnothing(prediction_window) ? seq_length : prediction_window
-    dummy_data = (q = rand(system_dim2, input_time_steps, n_params),
-        p = rand(system_dim2, input_time_steps, n_params))
-
-    dl = DataLoader(dummy_data)
-    batch = isnothing(prediction_window) ? Batch(batch_size, seq_length) :
-            Batch(batch_size, seq_length, prediction_window)
-
-    batch_indices_all = batch(dl)
-    for batch_indices in batch_indices_all
-        batched_array = convert_input_and_batch_indices_to_array(dl, batch, batch_indices)
-        @test size(batched_array[1].q) == (system_dim2, seq_length, batch_size)
-        @test size(batched_array[2].q) == (system_dim2, _prediction_window, batch_size)
-    end
-end
-
-test_data_loader_for_qp_tensor(5, 50, 20, 10, 16)
diff --git a/test/data_loader/parametric_data_loader.jl b/test/data_loader/parametric_data_loader.jl
index 9557d0b30..ca651a219 100644
--- a/test/data_loader/parametric_data_loader.jl
+++ b/test/data_loader/parametric_data_loader.jl
@@ -1,8 +1,14 @@
+# The one testset of the loader's files that needs this package's own names: `train!` on a
+# `ParametricDataLoader` trains on `ParametricLoss` by default, and the network is the parametric
+# affine layer the rest of `test/parametric.jl` builds. The loader and its own tests are
+# `ReducedComplexityModeling`'s.
+#
+# The network is `(x, μ) ↦ W [x; flatten(μ)] + b`. `train!` without `loss` and with
+# `loss = ParametricLoss()` draw the same network and the same minibatches under one seed.
+
 using GeometricMachineLearning
-using GeometricMachineLearning: ParametricDataLoader, ParametricLoss, params,
-                                convert_input_and_batch_indices_to_array
-using GeometricProblems.CoupledHarmonicOscillator: hodeensemble, default_parameters
-using GeometricIntegrators: ImplicitMidpoint, integrate
+using GeometricMachineLearning: ParametricLoss, params
+using ReducedComplexityModeling: ParametricDataLoader
 using Test
 import Random
 
@@ -13,153 +19,9 @@ include("../helpers/parametric_affine.jl")
 const d, nt, m = 1, 5, 3
 
 # The data of trajectory `j` at time step `t`, in every entry, and its system parameters `μⱼ`.
-constant_data(::Type{T}) where {T} = [T(j) for _ in 1:(2d), t in 1:nt, j in 1:m]
 timed_data(::Type{T}) where {T} = [T(100j + t) for _ in 1:(2d), t in 1:nt, j in 1:m]
 system_parameters(::Type{T}) where {T} = [(a = T(j),) for j in 1:m]
 
-# The loader of a `2d × nt × m` array, and of the `(q, p)` pair of its two halves.
-function loaders(data, μs)
-    (ParametricDataLoader(data, μs),
-        ParametricDataLoader((q = data[1:d, :, :], p = data[(d + 1):end, :, :]), μs))
-end
-
-# Entry `k` of a drawn batch, as the vector of the state it holds.
-sample(x::AbstractArray, k) = x[:, 1, k]
-sample(x::NamedTuple, k) = vcat(x.q[:, 1, k], x.p[:, 1, k])
-
-# Every sample of one epoch: its start position and trajectory, its input, its output and its `μ`.
-function epoch(dl, batch)
-    samples = []
-    for indices in batch(dl)
-        input, output, μs = convert_input_and_batch_indices_to_array(dl, batch, indices)
-        @test length(μs) == length(indices)
-        for (k, index) in enumerate(indices)
-            push!(samples, (index, sample(input, k), sample(output, k), μs[k]))
-        end
-    end
-    samples
-end
-
-all_pairs = sort(vec([(t, j) for t in 1:(nt - 1), j in 1:m]))
-
-@testset "one epoch of Batch(b) on a ParametricDataLoader ($T)" for T in REAL_ELTYPES
-    Random.seed!(1)
-    for dl in loaders(constant_data(T), system_parameters(T))
-        samples = epoch(dl, Batch(4))
-        # each pair of a start position and a trajectory once
-        @test sort([index for (index, _, _, _) in samples]) == all_pairs
-        for (_, input, output, μ) in samples
-            @test μ.a == input[1]
-            @test all(==(input[1]), input)
-        end
-        # the last minibatch is shorter: 12 pairs in minibatches of 5
-        @test length.(Batch(5)(dl)) == [5, 5, 2]
-    end
-    for dl in loaders(timed_data(T), system_parameters(T))
-        samples = epoch(dl, Batch(4))
-        @test sort([index for (index, _, _, _) in samples]) == all_pairs
-        for ((t, j), input, output, μ) in samples
-            @test input == fill(T(100j + t), 2d)
-            # the output is the input one time step later
-            @test output == input .+ 1
-            @test μ.a == T(j)
-            @test eltype(input) === eltype(output) === T
-        end
-    end
-end
-
-@testset "a ParametricDataLoader of an EnsembleSolution" begin
-    parameters = [default_parameters(), map(v -> v .+ 1, default_parameters())]
-    solution = integrate(hodeensemble(; parameters = parameters), ImplicitMidpoint())
-    dl = ParametricDataLoader(solution)
-    @test dl.parameters == parameters
-    @test dl.n_params == 2
-    n_steps = length(solution.t)
-    @test dl.input_time_steps == n_steps
-    for j in 1:2, t in 1:n_steps
-
-        @test dl.input.q[:, t, j] == solution.s[j].q[t - 1]
-        @test dl.input.p[:, t, j] == solution.s[j].p[t - 1]
-    end
-    Random.seed!(2)
-    samples = epoch(dl, Batch(7))
-    @test sort([index for (index, _, _, _) in samples]) ==
-          sort(vec([(t, j) for t in 1:(n_steps - 1), j in 1:2]))
-    for ((t, j), input, output, μ) in samples
-        @test μ == parameters[j]
-        @test input == vcat(dl.input.q[:, t, j], dl.input.p[:, t, j])
-        @test output == vcat(dl.input.q[:, t + 1, j], dl.input.p[:, t + 1, j])
-    end
-end
-
-# The message of the `ArgumentError` that `f` throws.
-function argument_error_message(f)
-    try
-        f()
-    catch error
-        error isa ArgumentError && return error.msg
-        rethrow()
-    end
-    "no error"
-end
-
-@testset "the edges of a ParametricDataLoader" begin
-    data = timed_data(Float64)
-    μs = system_parameters(Float64)
-    dl = ParametricDataLoader(data, μs)
-
-    # one time step per sample: a sequence length other than one is an error before any step
-    message = argument_error_message(() -> Batch(2, 3, 1)(dl))
-    @test occursin("sequence length of 3", message)
-    # a prediction window other than one is an error before any step, with the time steps to spare
-    message = argument_error_message(() -> Batch{:FeedForward}(2, 1, 2)(dl))
-    @test occursin("prediction window of 2", message)
-
-    # a parameter vector of the wrong length, and no trajectory
-    message = argument_error_message(() -> ParametricDataLoader(data, μs[1:2]))
-    @test occursin("hold 3 trajectories", message) &&
-          occursin("there are 2 system parameters", message)
-    message = argument_error_message(
-        () -> ParametricDataLoader(zeros(2d, nt, 0), NamedTuple{(:a,), Tuple{Float64}}[]))
-    @test occursin("hold 0 trajectories", message)
-
-    # too few time steps for the batch
-    message = argument_error_message(() -> Batch(2)(ParametricDataLoader(data[:, 1:1, :], μs)))
-    @test occursin("1 time step", message)
-
-    # one element type for the data and every `μ`
-    message = argument_error_message(() -> ParametricDataLoader(Float32.(data), μs))
-    @test occursin("Float32", message) && occursin("Float64", message)
-    @test ParametricDataLoader(data, [(a = 1.0, b = [2.0, 3.0]) for _ in 1:m]).n_params == m
-    # an empty `μ` has no element type, and concatenates nothing
-    @test ParametricDataLoader(data, [(;) for _ in 1:m]).n_params == m
-
-    # `q` and `p` of one size
-    message = argument_error_message(
-        () -> ParametricDataLoader((q = data[1:1, :, :], p = data[:, :, :]), μs))
-    @test occursin("(1, $(nt), $(m))", message) && occursin("(2, $(nt), $(m))", message)
-
-    # `q` and `p` of one array type, as `gather` batches them
-    @test_throws MethodError ParametricDataLoader(
-        (q = data[1:d, :, :], p = view(data, (d + 1):(2d), :, :)), μs)
-    views = ParametricDataLoader(
-        (q = view(data, 1:d, :, :), p = view(data, (d + 1):(2d), :, :)), μs)
-    input, output, _ = convert_input_and_batch_indices_to_array(views, Batch(2), [(1, 2)])
-    @test sample(input, 1) == data[:, 1, 2] && sample(output, 1) == data[:, 2, 2]
-
-    # every `μ` with the keys of the first, in its order and at every level, and its leaf sizes
-    for parameters in ([(a = 1.0,), (b = 2.0,), (a = 3.0,)],
-        [(a = 1.0, b = 2.0), (b = 3.0, a = 4.0), (a = 5.0, b = 6.0)],
-        [(a = 1.0,), (a = 2.0, b = 3.0), (a = 3.0,)],
-        [(a = (c = 1.0,),), (a = (e = 2.0,),), (a = (c = 3.0,),)],
-        [(a = [1.0, 2.0],), (a = [1.0, 2.0, 3.0],), (a = [1.0, 2.0],)])
-        message = argument_error_message(() -> ParametricDataLoader(data, parameters))
-        @test occursin("system parameter 2", message)
-    end
-end
-
-# The network is `(x, μ) ↦ W [x; flatten(μ)] + b`. `train!` without `loss` and with
-# `loss = ParametricLoss()` draw the same network and the same minibatches under one seed.
 @testset "train! on a ParametricDataLoader trains on ParametricLoss ($T)" for T in REAL_ELTYPES
     dl = ParametricDataLoader(timed_data(T) ./ 100, system_parameters(T))
     function trained(; kwargs...)
diff --git a/test/data_loader/train.jl b/test/data_loader/train.jl
index 0486d25ae..13a54c815 100644
--- a/test/data_loader/train.jl
+++ b/test/data_loader/train.jl
@@ -2,6 +2,7 @@ using GeometricMachineLearning
 using GeometricMachineLearning: NetworkLoss, params, ZygotePullback
 using AbstractNeuralNetworks: ArrayOrNamedTuple
 using NeuralNetworkParameters: NetworkParameters
+import ReducedComplexityModeling
 using Test
 import Random
 
@@ -97,7 +98,7 @@ end
     p₀ = deepcopy(params(nn))
     Random.seed!(8)
     losses = map(batch(dl)) do indices
-        data = GeometricMachineLearning.convert_input_and_batch_indices_to_array(
+        data = ReducedComplexityModeling.convert_input_and_batch_indices_to_array(
             dl, batch, indices)
         loss(nn.model, params(nn), data...)
     end
@@ -125,10 +126,11 @@ end
     @test NetworkLoss(GSympNet(4; n_layers = 2), Batch(4)) isa FeedForwardLoss
     autoencoder = SymplecticAutoencoder(10, 4)
     @test NetworkLoss(autoencoder, Batch(4)) isa AutoEncoderLoss
-    # a transformer with a `Batch{:FeedForward}` and an unknown architecture both reach the fallback
+    # a transformer with a `Batch{:FeedForward}` and an unknown architecture have no `NetworkLoss`
+    # method; `train!` turns that into an `ArgumentError` that names the keyword
     for (arch, batch) in ((transformer, Batch(4)),
         (NeuralNetwork(Chain(Dense(2, 2, identity))).architecture, Batch(4)))
-        @test_throws ArgumentError NetworkLoss(arch, batch)
+        @test_throws MethodError NetworkLoss(arch, batch)
     end
 
     nn = NeuralNetwork(Chain(Dense(2, 2, identity)))
@@ -140,6 +142,14 @@ end
     end
     @test raised isa ArgumentError
     @test occursin("`loss`", raised.msg)
+    nn = NeuralNetwork(transformer)
+    raised = try
+        train!(nn, TrainingOptimizer(params(nn)), dl, Batch(5), 1; show_progress = false)
+    catch e
+        e
+    end
+    @test raised isa ArgumentError
+    @test occursin("`loss`", raised.msg)
 
     # the same default trains on, through `train!` itself
     nn = NeuralNetwork(SymplecticAutoencoder(10, 4))
diff --git a/test/integration/docstrings/data_loader.jl b/test/integration/docstrings/data_loader.jl
deleted file mode 100644
index ca846430f..000000000
--- a/test/integration/docstrings/data_loader.jl
+++ /dev/null
@@ -1,47 +0,0 @@
-using Test
-using GeometricMachineLearning
-using GeometricMachineLearning: convert_input_and_batch_indices_to_array, number_of_batches
-import Random
-
-Random.seed!(1234)
-
-@testset "Data loader docstring examples" begin
-    data = [1 2 3; 4 5 6]
-    @test DataLoader(data; autoencoder = true, suppress_info = true).input ==
-          reshape(data, 2, 1, 3)
-    @test DataLoader(data; suppress_info = true).input == reshape(data, 2, 3, 1)
-
-    tensor = [1; 2; 3;; 4; 5; 6;;; 7; 8; 9;; 10; 11; 12]
-    expected_tensor = [1 4; 2 5; 3 6;;; 7 10; 8 11; 9 12]
-    @test DataLoader(tensor; suppress_info = true).input == expected_tensor
-    @test DataLoader(tensor; suppress_info = true) isa DataLoader
-
-    dl = DataLoader(rand(5); autoencoder = true, suppress_info = true)
-    batches = Batch(2)(dl)
-    @test length(batches) == 3
-    @test length.(batches) == [2, 2, 1]
-    @test sort(collect(union(batches...))) == [(1, i) for i in 1:5]
-    @test Batch(2, 3, 2) isa Batch{:Transformer}
-
-    dat = [1, 2, 3, 4, 5]
-    dl₁ = DataLoader(dat; autoencoder = false, suppress_info = true)
-    dl₂ = DataLoader(dat; autoencoder = true, suppress_info = true)
-    batch = Batch(3)
-    batches₁ = batch(dl₁)
-    batches₂ = batch(dl₂)
-    @test number_of_batches(dl₁, batch) == number_of_batches(dl₂, batch) == 2
-    @test length.(batches₁) == [3, 1]
-    @test length.(batches₂) == [3, 2]
-    @test sort(collect(union(batches₁...))) == [(i, 1) for i in 1:4]
-    @test sort(collect(union(batches₂...))) == [(1, i) for i in 1:5]
-
-    dl = DataLoader(collect(0.1:0.1:0.9); autoencoder = true, suppress_info = true)
-    indices = [(1, 1), (1, 3), (1, 5)]
-    @test convert_input_and_batch_indices_to_array(dl, Batch(3), indices) ==
-          reshape([0.1, 0.3, 0.5], 1, 1, 3)
-    @test number_of_batches(DataLoader([1, 2, 3, 4, 5]; suppress_info = true), Batch(2)) ==
-          2
-    @test number_of_batches(
-        DataLoader([1, 2, 3, 4, 5]; autoencoder = true, suppress_info = true), Batch(2)) ==
-          3
-end
diff --git a/test/integration/qp_dispatch.jl b/test/integration/qp_dispatch.jl
index b8252920f..a962a6bb1 100644
--- a/test/integration/qp_dispatch.jl
+++ b/test/integration/qp_dispatch.jl
@@ -6,8 +6,8 @@
 # `networkbackend`, which only map over a named tuple or read its first array, take any keys.
 
 using GeometricMachineLearning
-using GeometricMachineLearning: ResNetLayer, map_to_type, map_to_new_backend, _diff,
-                                ZygotePullback
+using GeometricMachineLearning: ResNetLayer, _diff, ZygotePullback
+using ReducedComplexityModeling: map_to_type, map_to_new_backend
 using AbstractNeuralNetworks: networkbackend, params
 using Test
 import Random
diff --git a/test/optimizers/optimizer.jl b/test/optimizers/optimizer.jl
index d95b17ace..b3c758284 100644
--- a/test/optimizers/optimizer.jl
+++ b/test/optimizers/optimizer.jl
@@ -1,5 +1,6 @@
 using GeometricMachineLearning
-using GeometricMachineLearning: params, convert_input_and_batch_indices_to_array
+using GeometricMachineLearning: params
+using ReducedComplexityModeling: convert_input_and_batch_indices_to_array
 using Test
 import Random
 import Zygote
diff --git a/test/parametric.jl b/test/parametric.jl
index f848179ed..31e67dfa4 100644
--- a/test/parametric.jl
+++ b/test/parametric.jl
@@ -179,7 +179,7 @@ end
 end
 
 @testset "the parametric names are public and not exported" begin
-    for name in (:apply_parametric, :ParametricLoss, :ParametricDataLoader)
+    for name in (:apply_parametric, :ParametricLoss)
         @test Base.ispublic(GeometricMachineLearning, name)
         @test !Base.isexported(GeometricMachineLearning, name)
     end
diff --git a/test/quality/aqua.jl b/test/quality/aqua.jl
index 120f9a583..69fe18beb 100644
--- a/test/quality/aqua.jl
+++ b/test/quality/aqua.jl
@@ -13,8 +13,25 @@ using Test
 # own, and it is also what lets the count at the end be a further assertion rather than a
 # top-level `@test` that aborts the file the moment it fails.
 #
-# `piracies` runs. The 3-tensor functors of `AbstractNeuralNetworks`' `Dense` and `Linear`, and the
-# `save`/`load` methods on `NeuralNetwork`, are `AbstractNeuralNetworks`' own.
+# `piracies` runs through `test_all` with `treat_as_own = [Batch]`, the one type this package
+# re-exports from `ReducedComplexityModeling` that its deliberate methods on another package's
+# function need. Two methods are reported without the list, and each names a `Batch`: this
+# package's `train!` on a network (`src/data_loader/train.jl:34-35`) and the keyword wrapper Julia
+# derives from that method. Both are deliberate -- the training loop is this package's, and the
+# other package owns the loader and the generic -- so naming `Batch` as this package's own keeps
+# Aqua's check on and passing rather than switching it off. An architecture and a batch with no
+# default loss are reported by `default_loss`, this package's own function, and not by a
+# `NetworkLoss` method that names only foreign types.
+#
+# `Batch` alone is the whole list. `train!` also names `DataLoader`, but in a `Union` with
+# `ParametricDataLoader`, which the list does not hold, so `DataLoader` in the list exempts no
+# method and hides a real piracy: a `Base.zero(::DataLoader)` fails this file under `[Batch]` and
+# passes under `[DataLoader, Batch]`. The check stays blind to a signature that names only
+# `Batch`, and to a symbol type parameter, which Aqua reads as a type of the package's own:
+# `Batch{:Transformer}` is such a parameter, so the method of `integrate` on that batch is gated in
+# `test/data_loader/batch.jl` instead. A piracy of any other foreign function, and a further method
+# of `train!` or `NetworkLoss` that names neither `Batch` nor a type of this package, still fails
+# here.
 #
 # `ambiguities` stays off, and the assertion below replaces it. Aqua's check runs in a process of its
 # own and asks about whatever that process loaded, so its count moves with the versions of the
@@ -23,7 +40,8 @@ using Test
 # switch stays off rather than becoming a gate on an upstream count that no line of this package
 # controls.
 @testset "Aqua" begin
-    Aqua.test_all(GeometricMachineLearning; ambiguities = false)
+    Aqua.test_all(GeometricMachineLearning; ambiguities = false,
+        piracies = (treat_as_own = [Batch],))
 
     @test isempty(Test.detect_ambiguities(GeometricMachineLearning; recursive = true))
 end
@@ -77,7 +95,7 @@ const NON_PUBLIC_QUALIFIED = (
     :JacobianResult, :jacobian, :value)                          # DiffResults
 const NON_PUBLIC_IMPORTS = (
     :AbstractExplicitLayer, :Architecture, :_compute_loss, :dim, :model, # AbstractNeuralNetworks
-    :load, :save)                                                        # NeuralNetworkParameters
+    :load, :save)                                                       # NeuralNetworkParameters
 
 @testset "every name reached in GeometricOptimizers is public" begin
     @test check_all_qualified_accesses_are_public(GeometricMachineLearning;
diff --git a/test/quality/exports.jl b/test/quality/exports.jl
index e589629fe..d380cc4f2 100644
--- a/test/quality/exports.jl
+++ b/test/quality/exports.jl
@@ -1,5 +1,6 @@
 using GeometricMachineLearning
 using Test
+import ReducedComplexityModeling
 
 # Julia only errors on a dangling `export` when the name is actually *resolved*, so an exported name
 # that nothing defines is silent: the package loads, the docs build, the suite passes, and
@@ -56,3 +57,42 @@ end
         @test !isdefined(GeometricMachineLearning, alias)
     end
 end
+
+# The data loader is `ReducedComplexityModeling`'s, and this package re-exports the three names a
+# training loop is written with. Each of them is that package's object and not a second binding of
+# the name, so a caller who loads both resolves one `DataLoader`, one `Batch` and one `train!`, and
+# the method this package defines on `train!` lands on the generic that package owns.
+# `ParametricDataLoader` is the one loader name this package does not hand to `using`: the
+# parametric family is this package's, the loader is not, so a caller names it
+# `GeometricMachineLearning.ParametricDataLoader` or takes it from its owner.
+@testset "DataLoader, Batch and train! are ReducedComplexityModeling's" begin
+    for name in (:DataLoader, :Batch, :train!)
+        @test getfield(GeometricMachineLearning, name) ===
+              getfield(ReducedComplexityModeling, name)
+        @test parentmodule(getfield(GeometricMachineLearning, name)) ===
+              ReducedComplexityModeling
+        @test Base.isexported(GeometricMachineLearning, name)
+        @test name in names(GeometricMachineLearning)
+    end
+
+    # a caller who loads this package alone reaches all three
+    m = Module()
+    Core.eval(m, :(using GeometricMachineLearning))
+    for name in (:DataLoader, :Batch, :train!)
+        @test Core.eval(m, name) === getfield(ReducedComplexityModeling, name)
+    end
+
+    @test GeometricMachineLearning.ParametricDataLoader ===
+          ReducedComplexityModeling.ParametricDataLoader
+    @test !Base.isexported(GeometricMachineLearning, :ParametricDataLoader)
+    @test !Base.ispublic(GeometricMachineLearning, :ParametricDataLoader)
+    @test !(:ParametricDataLoader in names(GeometricMachineLearning))
+
+    # `networkbackend` is `AbstractNeuralNetworks`' generic, and the loader's method of it is the
+    # other package's too
+    @test nameof(parentmodule(GeometricMachineLearning.networkbackend)) ===
+          :AbstractNeuralNetworks
+    @test GeometricMachineLearning.networkbackend(
+        ReducedComplexityModeling.DataLoader(reshape(Float32.(1:20), 2, 10);
+        suppress_info = true)) == CPU()
+end
diff --git a/test/quality/jet.jl b/test/quality/jet.jl
index 8e321d301..72671dff4 100644
--- a/test/quality/jet.jl
+++ b/test/quality/jet.jl
@@ -10,13 +10,12 @@ const TARGET = (GML,)
 const rrule = GML.ChainRulesCore.rrule
 
 # The entry points are the functions that launch a kernel of `src/` on KernelAbstractions' `CPU()`
-# backend, the pullbacks among them, and the `Batch` functor that
-# `test/data_loader/batch_index_set.jl` asserts with `@allocated`. A launcher line analyses one
-# entry point at the argument types of one element type. Where a test outside `test/quality/`
-# calls the entry point directly, the element types are those of the direct calls; an element
-# type that reaches the entry point only through another function gets no line. Where no test
-# calls the entry point directly, the element types are those at which the tests reach it through
-# its callers. Another container type at the same element type gets no line.
+# backend, the pullbacks among them. A launcher line analyses one entry point at the argument types
+# of one element type. Where a test outside `test/quality/` calls the entry point directly, the
+# element types are those of the direct calls; an element type that reaches the entry point only
+# through another function gets no line. Where no test calls the entry point directly, the element
+# types are those at which the tests reach it through its callers. Another container type at the
+# same element type gets no line.
 #
 # JET drops the reports of a kernel body when it analyses the launcher, so each `@kernel` method
 # also gets a line per element type of its launchers, which analyses the generated
@@ -87,72 +86,6 @@ const A3{T} = Array{T, 3}
             @test isempty(kernel_body_reports(GML.tensor_mat_skew_sym_assign_kernel!,
                 (4, 4, 2), zeros(T, 4, 4, 2), zeros(T, 3, 4, 2), zeros(T, 3, 3)))
         end
-
-        # src/data_loader/batch.jl
-        indices = ones(Int, 2, 4)
-        TimeSeriesQP = DataLoader{
-            Float64, @NamedTuple{q::A3{Float64}, p::A3{Float64}}, Nothing, :TimeSeries, Nothing}
-        Indices = Vector{Tuple{Int, Int}}
-        # test/data_loader/draw_batch_for_tensor_test.jl
-        @test isempty(JET.get_reports(JET.report_opt(
-            GML.convert_input_and_batch_indices_to_array,
-            (TimeSeriesQP, Batch{:Transformer}, Indices); target_modules = TARGET)))
-        # the one gather kernel of a `(q, p)` input: `gather` launches it once per tensor
-        @test isempty(kernel_body_reports(GML.gather_kernel!,
-            (2, 3, 4), zeros(2, 3, 4), zeros(2, 8, 4), indices, 3))
-        # test/integration/docstrings/data_loader.jl
-        @test isempty(JET.get_reports(JET.report_opt(
-            GML.convert_input_and_batch_indices_to_array,
-            (DataLoader{Float64, A3{Float64}, Nothing, :RegularData, Nothing},
-                Batch{:FeedForward},
-                Indices); target_modules = TARGET)))
-        for (T, BatchType) in ((Float32, Batch{:Transformer}), (
-            Float64, Batch{:FeedForward}))
-            # the optimizer functor in test/data_loader/training_history_eltype.jl
-            @test isempty(JET.get_reports(JET.report_opt(
-                GML.convert_input_and_batch_indices_to_array,
-                (DataLoader{T, A3{T}, Nothing, :TimeSeries, Nothing},
-                    Batch{:FeedForward},
-                    Indices);
-                target_modules = TARGET)))
-            # input and output, of kind `:TimeSeries`: no loader of the suite has this type.
-            # `DataLoader(input, output)` builds `:RegularData`, and each `:TimeSeries` loader
-            # with an output that a test builds by naming the type has an `Array{Int, 3}` output
-            @test isempty(JET.get_reports(JET.report_opt(
-                GML.convert_input_and_batch_indices_to_array,
-                (DataLoader{T, A3{T}, A3{T}, :TimeSeries, Int}, BatchType, Indices);
-                target_modules = TARGET)))
-            # the gather kernel, once per element type
-            @test isempty(kernel_body_reports(GML.gather_kernel!,
-                (2, 3, 4), zeros(T, 2, 3, 4), zeros(T, 2, 8, 4), indices, 3))
-        end
-        # the `Batch` functor under `@allocated` in test/data_loader/batch_index_set.jl, and
-        # the other element types of its direct calls: `Int` in
-        # test/integration/docstrings/data_loader.jl, `Float64` in
-        # test/data_loader/draw_batch_for_tensor_test.jl
-        @test isempty(JET.get_reports(JET.report_opt(Batch(8),
-            (DataLoader{Float32, A3{Float32}, Nothing, :TimeSeries, Nothing},);
-            target_modules = TARGET)))
-        @test isempty(JET.get_reports(JET.report_opt(Batch(3),
-            (DataLoader{Int, A3{Int}, Nothing, :TimeSeries, Nothing},);
-            target_modules = TARGET)))
-        @test isempty(JET.get_reports(JET.report_opt(Batch(10, 16), (TimeSeriesQP,);
-            target_modules = TARGET)))
-        # `Batch(batch_size, dl)` reads the loader's own two lengths, and
-        # test/data_loader/batch.jl calls it at the element types below
-        @test isempty(JET.get_reports(JET.report_opt(Batch(2),
-            (DataLoader{Float32, A3{Float32}, A3{Int}, :RegularData, Int},);
-            target_modules = TARGET)))
-        @test isempty(JET.get_reports(JET.report_opt(Batch(2),
-            (DataLoader{Float32, A3{Float32}, A3{Float32}, :TimeSeries, Int},);
-            target_modules = TARGET)))
-
-        # src/data_loader/tensor_assign.jl: `assign_output_estimate` in
-        # test/data_loader/accuracy.jl, in both precisions
-        for T in REAL_ELTYPES
-            @test isempty(JET.get_reports(JET.report_opt(GML.assign_output_estimate,
-                (A3{T}, Int); target_modules = TARGET)))
-        end
     else
         @test_skip "JET does not work on Julia $(VERSION)"  # aviatesk/JET.jl#681
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index a342f08ba..fce3c1274 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -91,19 +91,13 @@ if "core" in GROUPS
     @safetestset "Test data loader for q and p data" include("data_loader/batch_data_loader_qp_test.jl")
     @safetestset "Test the data loader in combination with optimization_step!" include("data_loader/data_loader_optimization_step.jl")
     @safetestset "One epoch of train! with Adam on a data loader" include("data_loader/optimizer_functor_with_adam.jl")
-    @safetestset "Test data loader for a tensor (q and p data)" include("data_loader/draw_batch_for_tensor_test.jl")
-    @safetestset "Batch functor(s)" include("data_loader/batch_functor.jl")
     @safetestset "DataLoader for input and output" include("data_loader/data_loader_for_input_and_output.jl")
-    @safetestset "An unexpected `autoencoder` keyword throws rather than returning nothing" include("data_loader/autoencoder_keyword.jl")
     @safetestset "A Float32 network's training history stays Float32" include("data_loader/training_history_eltype.jl")
     @safetestset "train!: convergence, resume, history type, n_epochs" include("data_loader/train.jl")
     @safetestset "train! takes a TrainingOptimizer's algorithm and linesearch" include("optimizers/optimizer.jl")
-    @safetestset "The minibatch index set is concrete and grows linearly" include("data_loader/batch_index_set.jl")
-    @safetestset "Classifier accuracy" include("data_loader/accuracy.jl")
-    @safetestset "DataLoader batching" include("data_loader/batch.jl")
-    @safetestset "DataLoader constructors" include("data_loader/data_loader.jl")
-    @safetestset "ParametricDataLoader batching and training" include("data_loader/parametric_data_loader.jl")
-    @safetestset "Data loader docstring examples" include("integration/docstrings/data_loader.jl")
+    @safetestset "The Batch method of integrate forwards its prediction window" include("data_loader/batch.jl")
+    @safetestset "solve! of a PSD architecture needs a samples loader" include("data_loader/data_loader.jl")
+    @safetestset "train! on a ParametricDataLoader trains on ParametricLoss" include("data_loader/parametric_data_loader.jl")
     @safetestset "Layer and architecture docstring examples" include("integration/docstrings/layers_and_architectures.jl")
     @safetestset "Loss docstring examples" include("integration/docstrings/losses.jl")
     @safetestset "Manifold docstring examples" include("integration/docstrings/manifolds.jl")