diff --git a/CHANGELOG.md b/CHANGELOG.md
index f2fc4308e..6a56ff986 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -17,6 +17,8 @@ breaking release).
 
 ### Removed (breaking)
 
+- **The Hamiltonian and the Lagrangian neural networks move to two new packages, and SymbolicNeuralNetworks and Symbolics leave `[deps]` and `[compat]`.** [HamiltonianNeuralNetworks.jl](https://github.com/JuliaGML/HamiltonianNeuralNetworks.jl) takes `HamiltonianArchitecture`, `StandardHamiltonianArchitecture`, `GeneralizedHamiltonianArchitecture`, `HNNLoss`, `SymplecticEulerLoss` and `SymbolicPullback(::HamiltonianArchitecture)`, with their tests, the architecture page and the tutorial of the Hamiltonian neural network, and the scripts `hnn_pendulum.jl`, `hnn_pendulum_simulation.jl` and `symplectic_euler_loss_modified_hamiltonian.jl`. [LagrangianNeuralNetworks.jl](https://github.com/JuliaGML/LagrangianNeuralNetworks.jl) takes `LagrangianNeuralNetwork`, `LNNLoss` and `VariationalMidpointLoss`, with their tests, the architecture page of the Lagrangian neural network and the script `lnn_pendulum.jl`. GeometricMachineLearning no longer defines or exports these names, nor re-exports `SymbolicPullback`, and `scripts/utilities/plots.jl` and `scripts/utilities/data_problem.jl`, which only those scripts included, go with them, as do `∇H`, `dH` and `get_data_set` of `scripts/utilities/pendulum.jl`, which only `hnn_pendulum.jl` and `hnn_pendulum_simulation.jl` called. Migrate code that uses these names to `using GeometricMachineLearning, HamiltonianNeuralNetworks` or `using GeometricMachineLearning, LagrangianNeuralNetworks`; `DataLoader`, `Batch` and `train!` stay here. Known issues C7 and K15 are in HamiltonianNeuralNetworks.jl now; C20, a missing Knet.jl example of the Hamiltonian neural network, and K21, the field name `act` of `LagrangianNeuralNetwork`, leave this package's `KNOWN_ISSUES.md` because the code they describe leaves it; neither new package records them yet. Where another entry of this section names these types, losses or `SymbolicPullback`, that part describes code that ships in the two new packages, not in GeometricMachineLearning, and the `SymbolicNeuralNetworks` bound of the `[compat]` entry below does not apply.
+
 - **`NetworkLoss(::TransformerIntegrator, prediction_window)` is gone; the default loss is `NetworkLoss(nn.architecture, batch)`.** `train!`'s default is now the two-argument `NetworkLoss(arch, batch)`, so a transformer's prediction window comes from the `Batch` it trains on instead of from the architecture. The old method took `prediction_window = 1`, which against `Batch`'s default `seq_length` was a conflict; `NetworkLoss(transformer, Batch(seq, pred))` replaces it, and `NetworkLoss(transformer, n::Integer)` is a `MethodError`. Migrate `NetworkLoss(transformer)` and `NetworkLoss(transformer, n)` to `NetworkLoss(transformer, batch)` or `TransformerLoss(batch)`.
 
 - **`train!`'s `pullback` keyword is gone.** `loss` takes a `NetworkLoss` or an `AbstractPullback`, and a `NetworkLoss` is trained through a `ZygotePullback`. Migrate `train!(nn, o, dl, batch, n; loss = HNNLoss(arch), pullback = SymbolicPullback(arch))` to `train!(nn, o, dl, batch, n; loss = SymbolicPullback(arch))`. `train!` builds the pullback once per call, not once per step.
@@ -122,11 +124,11 @@ breaking release).
 
 - `Pkg.test()` no longer runs the doctests. `test/quality/doctests.jl` is now the group `doctests`, which an empty `ARGS` does not run; `Pkg.test(test_args = ["doctests"])` runs it. In CI the Doctests job stays their runner, so the test matrix no longer runs them a second time.
 - **CI uploads coverage from the `Julia 1 - ubuntu-latest` job instead of `Julia min`, and a test job saves the Julia cache only when it succeeds.**
-- **The package needs Julia 1.12 or later** (`[compat] julia = "1.12"`, was `"1.11"`). On Julia 1.11.9, the gradient of an `LNNLoss` with respect to the network parameters stops with `StackOverflowError` in type inference ("Internal error: during type inference of … Encountered stack overflow") on macOS and Windows; the `Julia min` job on Linux passes. The overflow comes from SymbolicUtils 4.49.0, whose `rrule` for `SymbolicUtils.Code.create_array` returns an exactly typed tuple of tangents, one per array element; Zygote's reverse pass through the generated Hessian of `lagrangian_acceleration` is then inferred far deeper. With the 4.48.0 rule restored, or in a task with a 64 MB stack, the same gradient passes on 1.11.9; Julia 1.12.7 and 1.13.1 pass it on the root task. This closes #327: the `Julia min` jobs on macOS and Windows had failed in the *Lagrangian Neural Network* testset since 2026-10-01, and every testset after it was never run there.
+- **The package needs Julia 1.12 or later** (`[compat] julia = "1.12"`, was `"1.11"`). The floor comes from `[compat] GeometricOptimizers = "0.9"`, whose v0.9.0 `Project.toml` has `julia = "1.12"`. It also closes #327, a `StackOverflowError` in type inference of the gradient of an `LNNLoss` on Julia 1.11.9, where the overflow came from SymbolicUtils 4.49.0's `rrule` for `SymbolicUtils.Code.create_array`; it failed the `Julia min` jobs on macOS and Windows in the *Lagrangian Neural Network* testset from 2026-10-01, and every testset after it never ran there. That loss now ships in LagrangianNeuralNetworks.jl, not here, so the floor does not rest on it.
 - **`test/Project.toml` has no `[compat]` entry for a dependency of the package.** The entries for AbstractNeuralNetworks, GeometricOptimizers, KernelAbstractions, LinearAlgebra, NeuralNetworkParameters, Random and Zygote are gone: the test environment contains the package, so the bounds of `Project.toml` apply there, and an entry in `test/Project.toml` can only repeat or narrow them. `NeuralNetworkParameters = "0.3"` was narrower than the package's `"0.3, 0.4"`.
-- **The test suite follows one flat layout.** `test/runtests.jl` is one list of `@safetestset "<label>" include("<path>")` lines in the groups `core`, `slow` and `metal`, selected by `ARGS`; with no arguments it runs `core` and `slow`, and `metal` too on an Apple-silicon Mac, so `Pkg.test()` and CI run every test they ran before. The 14 nested `test/<subject>/runtests.jl` drivers are gone. A test file sits where the source it tests sits: `test/attention/` and `test/transformers/` went into `test/layers/`, `test/architectures/` and `test/kernels/`; `test/reduced_order_modeling/` into `test/architectures/` and `test/reduced_system/`; `test/losses/training_method_losses.jl` into `test/loss/`; the tests that span the package, the docstring examples among them, into `test/integration/`; and `aqua.jl`, `exports.jl` and `reachability.jl` into `test/quality/`. `test/docstrings/layers_and_architectures.jl`, named in *Added* above, is now `test/integration/docstrings/layers_and_architectures.jl`. `test/layers/sympnet_upscaling.jl` is in `slow`. Thirteen files that drew random numbers without a seed now call `Random.seed!(1234)`. The shared helper of the two gradient-structure tests is `test/helpers/network_parameters_gradient_structure.jl`, and its `@test_broken` names issue #319.
+- **The test suite follows one flat layout.** `test/runtests.jl` is one list of `@safetestset "<label>" include("<path>")` lines in the groups `core`, `slow` and `metal`, selected by `ARGS`; with no arguments it runs `core` and `slow`, and `metal` too on an Apple-silicon Mac, so `Pkg.test()` and CI run every test they ran before. The 14 nested `test/<subject>/runtests.jl` drivers are gone. A test file sits where the source it tests sits: `test/attention/` and `test/transformers/` went into `test/layers/`, `test/architectures/` and `test/kernels/`; `test/reduced_order_modeling/` into `test/architectures/` and `test/reduced_system/`; the tests that span the package, the docstring examples among them, into `test/integration/`; and `aqua.jl`, `exports.jl` and `reachability.jl` into `test/quality/`. `test/docstrings/layers_and_architectures.jl`, named in *Added* above, is now `test/integration/docstrings/layers_and_architectures.jl`. `test/layers/sympnet_upscaling.jl` is in `slow`. Thirteen files that drew random numbers without a seed now call `Random.seed!(1234)`. The shared helper of the two gradient-structure tests is `test/helpers/network_parameters_gradient_structure.jl`, and its `@test_broken` names issue #319.
 - **The doctests run in the test suite.** `test/quality/doctests.jl`, in the `slow` group, runs `doctest(GeometricMachineLearning)` as the `Doctests` job of `CI.yml` does. `Documenter` joins `test/Project.toml`, with no `[compat]` bound, as a test-only dependency.
-- **The test suite compiles less.** `test/layers/sympnet_upscaling.jl` checks its three symplectic identities on 12 pairs `(N, N2)` instead of all 120 with `N` in `2:2:20` and `N2` in `(2N):2:(4N)`: `(2, 4)`, `(2, 6)`, `(2, 8)`, `(4, 10)`, `(6, 12)`, `(6, 20)`, `(8, 24)`, `(10, 30)`, `(12, 36)`, `(14, 42)`, `(20, 40)` and `(20, 80)`. They keep the smallest pair, `N2 = 2N`, `N2 = 4N`, `N2` not a multiple of `N`, and the largest `N`. The layers carry the sizes in their types, so each pair is a new compilation (about 1.3 s per pair, not once per file as the [0.8.0] entry assumed); the file has 36 assertions instead of 360, and takes 19.7 s cold; it stays in `slow`. The two training tests of `test/loss/training_method_losses.jl` use the network shape of the rest of the file, `StandardHamiltonianArchitecture(2n, 4, 1)` and `LagrangianNeuralNetwork(2; width = 4, nhidden = 1)`, in place of `(2n, 5, 2)` and `width = 5, nhidden = 2`; each still asserts that training reduces the loss. `test/architectures/lagrangian_neural_network_tests.jl` moves from `core` to `slow`: it takes 122.7 s cold, above the 60 s of a `core` file. It still runs in CI and in `Pkg.test()`, and leaves only the local `core` run. Every other `core` file takes at most 37.9 s cold. The times are cold, from one fresh Julia process that runs the `core` files in the order of `test/runtests.jl` (and from one fresh process for each of the two files named), on an Apple M4 Max with Julia 1.13.1, at a load average below 4 at the start.
+- **The test suite compiles less.** `test/layers/sympnet_upscaling.jl` checks its three symplectic identities on 12 pairs `(N, N2)` instead of all 120 with `N` in `2:2:20` and `N2` in `(2N):2:(4N)`: `(2, 4)`, `(2, 6)`, `(2, 8)`, `(4, 10)`, `(6, 12)`, `(6, 20)`, `(8, 24)`, `(10, 30)`, `(12, 36)`, `(14, 42)`, `(20, 40)` and `(20, 80)`. They keep the smallest pair, `N2 = 2N`, `N2 = 4N`, `N2` not a multiple of `N`, and the largest `N`. The layers carry the sizes in their types, so each pair is a new compilation (about 1.3 s per pair, not once per file as the [0.8.0] entry assumed); the file has 36 assertions instead of 360, and takes 19.7 s cold; it stays in `slow`. Every `core` file takes at most 37.9 s cold. The times are cold, from one fresh Julia process that runs the `core` files in the order of `test/runtests.jl` (and from one fresh process for the file named), on an Apple M4 Max with Julia 1.13.1, at a load average below 4 at the start.
 - **The optimizer reads `show_progress` as a `Bool`** instead of comparing it with `== true`. A non-`Bool` value now raises a `TypeError`: before, `nothing`, `0` or a string hid the progress bar with no error, and `1` or `1.0` showed it. Every call in the tests and the manual passes `false`, and every call in `scripts/` passes `true`.
 - **fatou reports no finding in `src/` or `test/`.** Only tests and one lint suppression changed; no behaviour of the package changed. `test/layers/multi_head_attention_stiefel_setup.jl` defined `check_setup(ps::NetworkParameters)` twice, and the later method silently replaced the `mapparameters` one, which never ran. The later method is deleted, and the test now asserts that the matrix method of `check_setup` runs once per leaf of the parameters, counted with `foldparameters`. The unused destructured values `Σ`, `Vt` and `Ũ₂` in `test/optimizers/psd_optim.jl` and `test/optimizers/svd_optim.jl` are now `_`; those files check `Ũ₁` and `Ũ₃` and never checked `Ũ₂`. The unreachable `new{…}(…)` after the `error` in the `GeneralizedHamiltonianArchitecture` constructor stays, because it records the field order for the future constructor, under `# fatou-ignore unreachable-code`.
 - **Three test files left the top level of `test/`**, because a test file sits there only where it mirrors `src/<name>.jl`. `test/norm_eltype.jl` tests `_norm` of `src/utils.jl` and is now `test/utils.jl`; `test/map_to_cpu_tests.jl` tests `src/map_to_cpu.jl` and is now `test/map_to_cpu.jl`; `test/documented_defaults.jl` tests the keyword defaults of `src/layers/symplectic_attention.jl`, `src/architectures/symplectic_transformer.jl` and `src/architectures/standard_transformer_integrator.jl`, whose deepest common directory is `src/` itself, and is now `test/integration/documented_defaults.jl`. Only the paths in `test/runtests.jl` changed; the labels, groups and order are the same.
diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md
index e43f11ad6..63e2b2c18 100644
--- a/KNOWN_ISSUES.md
+++ b/KNOWN_ISSUES.md
@@ -15,23 +15,6 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
   pattern, so the working tree and the net diff are clean, but the blobs are still reachable.
   Rewriting the branch would remove them.
 
-### C7 · `SymbolicPullback(::HamiltonianArchitecture)` reaches one private SymbolicNeuralNetworks name.
-
-- **location:** `src/pullbacks/symbolic_hnn_pullback.jl`
-- **kind:** upstream
-- **found:** 2026-08-16
-- **evidence:**
-  `SymbolicNeuralNetworks.SymbolicPullback(nn, loss)` cannot be used, because it sizes the symbolic
-  target of the loss with `output_dimension(nn.model)`, which for a Hamiltonian network is `1` — the
-  scalar Hamiltonian — while `HNNLoss` compares against the Hamiltonian *vector field*, which has the
-  dimension of the network's input. GML's construction therefore reaches
-  `SymbolicNeuralNetworks.symbolic_parameter_gradient`, which no exported name replaces:
-  `flat_parameter_gradient` is flat, so it wants an `unflatten` per step, and `ParameterGradient` is
-  not used. The pullback is the upstream type `SymbolicPullback(loss, fun)`, with a call method of
-  this package for `SymbolicPullback{<:HNNLoss, <:_HNNParameterGradient}` that is more specific
-  than the upstream call in the pullback and the parameters, and as specific in the data. That
-  private name is the one an upstream refactor breaks.
-
 ### C14 · `evaluate_vf_and_compute_∇Ψ` evaluates the decoder twice.
 
 - **location:** `src/reduced_system/reduced_system.jl`
@@ -58,24 +41,15 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
 - **kind:** upstream
 - **found:** 2026-09-18
 - **evidence:**
-  `no_implicit_imports` is off: **37** names arrive through a bare `using` across thirteen packages,
+  `no_implicit_imports` is off: **33** names arrive through a bare `using` across twelve packages,
   and naming them is a rewrite of the module header with a judgement per name.
   `all_explicit_imports_are_public` and `all_qualified_accesses_are_public` run with a list of the
-  names this package reaches in other packages that those packages do not declare `public`: **10**
-  explicit imports and **9** qualified accesses, none of them in `GeometricOptimizers`.
+  names this package reaches in other packages that those packages do not declare `public`: **9**
+  explicit imports and **4** qualified accesses, none of them in `GeometricOptimizers`.
   `Architecture`, `AbstractExplicitLayer`, `_compute_loss`, `dim` and the rest are used
   deliberately, and the resolution is upstream declaring them `public`. Each check's reason is
   written beside it in the test file.
 
-### C20 · There is no Knet.jl example for the Hamiltonian neural network.
-
-- **location:** `TODO.md`
-- **kind:** docs
-- **found:** 2026-09-19
-- **evidence:**
-  The deleted `TODO.md` asked
-  for one. Recorded rather than dropped; nothing depends on it.
-
 ### C21 · "An untrained autoencoder reduces worse than PSD" is not a property, and the suite no longer asserts it under an explicit integrator.
 
 - **location:** `test/reduced_system/reduced_system.jl`
@@ -141,28 +115,6 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
   on 1.13 and *rose* on 1.12". The fix states only the present constraint: a 32-epoch budget
   drives the learning rate to `η₂` almost at once, and 128 epochs train.
 
-### K15 · `HNNLoss` and `SymbolicPullback(arch)` accept `(q, p)` data in their signatures and fail on it inside the generated vector field.
-
-- **location:** `src/loss/hnn_loss.jl:33`
-- **kind:** found late
-- **found:** #344, review round 2 (2026-10-14); the `HNNLoss` signatures are older than that pull request
-- **evidence:**
-  `HNNLoss`'s two functors (`src/loss/hnn_loss.jl:31-42`) and the call method of
-  `SymbolicPullback(arch)` (`src/pullbacks/symbolic_hnn_pullback.jl:38-40`) type the data as
-  `ArrayOrNamedTuple`, which admits a `(q = …, p = …)` `NamedTuple`. The loss hands that data to
-  the generated `hvf`, which calls `ndims` on it. With `(q = randn(1, 5), p = randn(1, 5))` as
-  input and output, the symbolic pullback, `HNNLoss` and `ZygotePullback(HNNLoss(arch))` all
-  raise:
-
-  ```
-  MethodError: no method matching ndims(::@NamedTuple{q::Matrix{Float64}, p::Matrix{Float64}})
-  ```
-
-  Narrowing the pullback's method alone changes nothing: the data then reaches the upstream
-  `SymbolicPullback` call, which also accepts a `NamedTuple` and fails in the same place. The fix
-  belongs with the data signatures of the losses as a whole: either each loss accepts only the data
-  its vector field can take, or `(q, p)` data is converted before it reaches the field.
-
 ### K17 · `assign_q_and_p` and `PoissonTensor * v` drop the rows from `2N + 1` on without an error.
 
 - **location:** `src/kernels/assign_q_and_p.jl`
@@ -225,22 +177,6 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
   `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.
 
-### K21 · `LagrangianNeuralNetwork` stores its activation in the field `act`, where every other architecture with one activation names it `activation`.
-
-- **location:** `src/architectures/lagrangian_neural_network.jl:5`
-- **kind:** found late
-- **found:** 2026-10-20, by the critics of part M5; the field is the same on `origin/main` at `4c4a8320`
-- **evidence:**
-  ```
-  $ grep -rn " act::\|\.act\b" src/architectures/
-  src/architectures/lagrangian_neural_network.jl:5:    act::AT
-  src/architectures/lagrangian_neural_network.jl:16:        [Dense(nn.width, nn.width, nn.act) for _ in 1:nn.nhidden]
-  src/architectures/lagrangian_neural_network.jl:20:        Dense(nn.dimin, nn.width, nn.act),
-  ```
-
-  The constructor's keyword is `activation`; only the field keeps the short name. The fix renames
-  the field, which changes a public field name.
-
 ### K22 · A comment in `symplectic_attention.jl` names a test file under `test/parameters/`, which does not exist.
 
 - **location:** `src/layers/symplectic_attention.jl:131`
diff --git a/Project.toml b/Project.toml
index 07a5b1313..ef30d3a90 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,8 +18,6 @@ NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 NeuralNetworkParameters = "67f4d93a-60e9-472b-8cdd-1ccf6005724a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SymbolicNeuralNetworks = "aed23131-dcd0-47ca-8090-d21e605652e3"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
@@ -45,8 +43,6 @@ NNlib = "0.8, 0.9"
 NeuralNetworkParameters = "0.5"
 ProgressMeter = "1"
 Random = "1"
-SymbolicNeuralNetworks = "0.9"
-Symbolics = "7"
 TimerOutputs = "0.5, 1"
 Zygote = "0.7"
 julia = "1.12"
diff --git a/README.md b/README.md
index 77e6c9ad8..c29772fc8 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,7 @@
 [![Build Status](https://github.com/JuliaGNI/GeometricMachineLearning.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaGNI/GeometricMachineLearning.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/JuliaGNI/GeometricMachineLearning.jl/branch/main/graph/badge.svg?token=CFT76RROW2)](https://codecov.io/gh/JuliaGNI/GeometricMachineLearning.jl)
 
-`GeometricMachineLearning.jl` offers a flexible tool for designing neural networks for dynamical systems with geometric structure, such as Hamiltonian (symplectic) or Lagrangian (variational) systems.
+`GeometricMachineLearning.jl` offers a flexible tool for designing neural networks for dynamical systems with geometric structure, such as Hamiltonian (symplectic) or divergence-free (volume-preserving) systems. Hamiltonian and Lagrangian neural networks are in [HamiltonianNeuralNetworks.jl](https://juliagml.github.io/HamiltonianNeuralNetworks.jl/latest/) and [LagrangianNeuralNetworks.jl](https://juliagml.github.io/LagrangianNeuralNetworks.jl/dev/).
 
 At its core every neural network comprises three components: a neural network architecture, a loss function and an optimizer. 
 
diff --git a/docs/make.jl b/docs/make.jl
index 879235b12..9d8c1e01d 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -210,8 +210,6 @@ _architectures = "Architectures" => [
     "Using Architectures with `NeuralNetwork`" => "architectures/abstract_neural_networks.md",
     "Symplectic Autoencoders" => "architectures/symplectic_autoencoder.md",
     "Neural Network Integrators" => "architectures/neural_network_integrators.md",
-    "Hamiltonian Neural Network" => "architectures/hamiltonian_neural_network.md",
-    "Lagrangian Neural Network" => "architectures/lagrangian_neural_network.md",
     "SympNet" => "architectures/sympnet.md",
     "Volume-Preserving FeedForward" => "architectures/volume_preserving_feedforward.md",
     "Standard Transformer" => "architectures/transformer.md",
@@ -227,7 +225,6 @@ _data_loader = "Data Loader" => [
 
 _tutorials = "Tutorials" => [
     "SympNets" => "tutorials/sympnet_tutorial.md",
-    "Hamiltonian Neural Network" => "tutorials/hamiltonian_neural_network.md",
     "Symplectic Autoencoders" => "tutorials/symplectic_autoencoder.md",
     "Grassmann Manifold" => "tutorials/grassmann_layer.md",
     "Volume-Preserving Attention" => "tutorials/volume_preserving_attention.md",
diff --git a/docs/src/architectures/hamiltonian_neural_network.md b/docs/src/architectures/hamiltonian_neural_network.md
deleted file mode 100644
index 36feb569b..000000000
--- a/docs/src/architectures/hamiltonian_neural_network.md
+++ /dev/null
@@ -1,59 +0,0 @@
-# [Hamiltonian Neural Network](@id hnn_architecture)
-
-The Hamiltonian Neural Network (HNN) [greydanus2019hamiltonian](@cite) aims at building a [Hamiltonian vector field](@ref "Symplectic Systems")  with a neural network. We recall that a *canonical Hamiltonian vector field* on ``\mathbb{R}^{2d}`` is one that can be written as:
-
-```math
-    X_H(z) = \mathbb{J}_{2d}\nabla_zH,
-```
-where ``\mathbb{J}_{2d}`` is the [`PoissonTensor`](@ref). The idea behind a Hamiltonian neural network is to learn a vector field of this form, i.e. to learn:
-
-```math
-    X_{\mathcal{NN}}(z) = \mathbb{J}_{2d}\nabla_z\mathcal{NN},
-```
-where ``\mathcal{NN}:\mathbb{R}^{2d}\to\mathbb{R}`` is a neural network that approximates the Hamiltonian. There are then two different options to define a *HNN loss*, depending on the format in which the data are given.
-
-## HNN Loss for Vector Field Data
-
-For the first loss, we assume that the given data describe the vector field of the HNN:
-
-```math
-\mathcal{L}_\mathrm{HNN} = \sqrt{\sum_{i=1}^d\left(\Big|\Big|\frac{\partial\mathcal{NN}}{\partial{}q_i} + \dot{p}_i \Big|\Big|_2^2 + \Big|\Big| \frac{\partial\mathcal{NN}}{\partial{}p_i} - \dot{q}_i\Big|\Big|_2^2\right)}
-```
-
-## HNN Loss for Phase Space Data
-
-For the second loss, we assume that the given data describe points in phase space associated to a Hamiltonian system. For this approach we also need to specify a *symplectic integrator* [hairer2006geometric](@cite) in order to train the neural network. In the following we use [SymplecticEulerB](https://juliagni.github.io/GeometricIntegrators.jl/latest/modules/integrators/#GeometricIntegrators.Integrators.SymplecticEulerB) to define this loss. This integrator does the following:
-
-```math
-\mathrm{SymplecticEulerB}: (q^{(t)}, p^{(t)}) \mapsto (q^{(t+1)}, p^{(t+1)})
-```
-
-Note that this integrator is implicit in general.
-
-```math
-\mathcal{L}_\mathrm{HNN} = \sqrt{\sum_t\sum_{i=1}^d \Big|\Big| \frac{\partial\mathcal{NN}}{\partial{}q_i}(q^{(t)}, p^{(t+1)}) + \frac{p_i^{t+1} - p_i^{(t)}}{\Delta{}t} \Big|\Big|_2^2 + \Big|\Big| \frac{\partial\mathcal{NN}}{\partial{}p_i}(q^{(t)}, p^{(t+1)}) + \frac{q_i^{t+1} - q_i^{(t)}}{\Delta{}t} \Big|\Big|_2^2}
-```
-
-Here the derivatives (i.e. vector field data) ``\dot{q}_i^{(t)}`` and ``\dot{p}_i^{(t)}`` are approximated with finite differences: 
-
-A loss of this kind is [`GeometricMachineLearning.SymplecticEulerLoss`](@ref). Its `:B` variant evaluates the Hamiltonian vector field at ``(q^{(t)}, p^{(t+1)})``, as above, and its `:A` variant at ``(q^{(t+1)}, p^{(t)})``; it compares that field against the finite difference, with the sign convention ``X_H = (\partial{}H/\partial{}p, -\partial{}H/\partial{}q)``.
-
-!!! warning
-    The displayed formula is not what that loss computes. Under the convention above, the second term reads ``-(q_i^{t+1} - q_i^{(t)})/\Delta{}t``, and the formula has a plus. The sign here predates the implementation and is left as it stands; the implementation is the one to read.
-
-!!! info
-    Usually we use [`Zygote`](https://github.com/FluxML/Zygote.jl) for computing derivatives in `GeometricMachineLearning`, but as the [`Zygote` documentation](https://fluxml.ai/Zygote.jl/dev/limitations/#Second-derivatives-1) itself points out: "Often using a different AD system over Zygote is a better solution [for computing second-order derivatives]." For this reason we compute the loss of the HNN with [`SymbolicNeuralNetworks`](https://github.com/JuliaGNI/SymbolicNeuralNetworks.jl) and optionally also its gradient.
-
-## Library Functions
-
-```@docs
-GeometricMachineLearning.hamiltonian_vector_field(::HamiltonianArchitecture)
-GeometricMachineLearning.HamiltonianArchitecture
-GeometricMachineLearning.StandardHamiltonianArchitecture
-GeometricMachineLearning.HNNLoss
-GeometricMachineLearning.SymplecticEulerLoss
-GeometricMachineLearning._symplectic_euler_evaluation_point
-GeometricMachineLearning.symbolic_hamiltonian_vector_field(::GeometricMachineLearning.SymbolicNeuralNetwork)
-GeometricMachineLearning.SymbolicPullback(::HamiltonianArchitecture)
-GeometricMachineLearning.GeneralizedHamiltonianArchitecture
-```
\ No newline at end of file
diff --git a/docs/src/architectures/lagrangian_neural_network.md b/docs/src/architectures/lagrangian_neural_network.md
deleted file mode 100644
index 3bdc22d88..000000000
--- a/docs/src/architectures/lagrangian_neural_network.md
+++ /dev/null
@@ -1,75 +0,0 @@
-# Lagrangian Neural Network
-
-A *Lagrangian neural network* learns a scalar Lagrangian ``L(q, \dot{q})`` instead of the dynamics
-themselves. The dynamics are then whatever that Lagrangian implies, through the Euler–Lagrange
-equations
-
-```math
-    \frac{d}{dt}\frac{\partial{}L}{\partial\dot{q}} = \frac{\partial{}L}{\partial{}q}.
-```
-
-This is the Lagrangian counterpart of the [Hamiltonian neural network](@ref hnn_architecture): there
-the network is a Hamiltonian and the dynamics come from a Poisson tensor, here the network is a
-Lagrangian and the dynamics come from a variational principle.
-
-## The loss
-
-Expanding the Euler–Lagrange equations along a trajectory gives
-
-```math
-    \nabla_{\dot{q}}\nabla_{\dot{q}}L\,\ddot{q}
-        + (\nabla_q\nabla_{\dot{q}}L)^T\dot{q}
-        = \nabla_qL,
-```
-
-so the acceleration the learned Lagrangian predicts is
-
-```math
-    \ddot{q} = (\nabla_{\dot{q}}\nabla_{\dot{q}}L)^{-1}
-               \left( \nabla_qL - (\nabla_q\nabla_{\dot{q}}L)^T\dot{q} \right),
-```
-
-which is what [`GeometricMachineLearning.lagrangian_acceleration`](@ref) computes and what
-[`GeometricMachineLearning.LNNLoss`](@ref) compares against the data.
-
-The transpose is not cosmetic. ``\nabla_q\nabla_{\dot{q}}L`` is indexed
-``[i, j] = \partial^2L/\partial{}q_i\partial\dot{q}_j``, while the chain rule contracts the *first*
-index with ``\dot{q}``. It disappears only for a Lagrangian whose position–velocity coupling happens
-to be symmetric, which is why the test for this asserts a closed form whose coupling deliberately is
-not.
-
-## Why the derivatives are symbolic
-
-The loss needs the *second* derivative of the network with respect to its input, and then its own
-derivative with respect to the network parameters. Taking the inner one with `Zygote` inside the
-loss does not work: the parameter gradient then fails with
-`MethodError: no method matching getindex(::IdDict{Any, Any})`. The gradient and Hessian are
-therefore built as symbolic expressions with
-[`SymbolicNeuralNetworks`](https://github.com/JuliaGNI/SymbolicNeuralNetworks.jl) and compiled, so
-that only one differentiation is left for `Zygote` to do. This is the same route
-`hamiltonian_vector_field` takes, and for the same reason.
-
-## Training on positions alone
-
-[`GeometricMachineLearning.LNNLoss`](@ref) needs ``\ddot{q}`` in the data. When only positions are
-available, [`GeometricMachineLearning.VariationalMidpointLoss`](@ref) trains the same architecture
-through the discrete Euler–Lagrange equations of the *midpoint discrete Lagrangian*
-
-```math
-    L_d(q_n, q_{n+1}) = \Delta{}t\,
-        L\!\left( \frac{q_n + q_{n+1}}{2}, \frac{q_{n+1} - q_n}{\Delta{}t} \right),
-```
-
-which ask that ``D_2L_d(q_n, q_{n+1}) + D_1L_d(q_{n+1}, q_{n+2}) = 0`` along the data. Only the
-first derivative of the network is needed for this, because the chain rule through the midpoint is
-written out rather than taken with a nested `Zygote.gradient`.
-
-## Library Functions
-
-```@docs
-GeometricMachineLearning.LNNLoss
-GeometricMachineLearning.lagrangian_acceleration
-GeometricMachineLearning._euler_lagrange_acceleration
-GeometricMachineLearning.VariationalMidpointLoss
-GeometricMachineLearning._discrete_lagrangian_derivatives
-```
diff --git a/docs/src/architectures/sympnet.md b/docs/src/architectures/sympnet.md
index 9683d4462..b8c15bf56 100644
--- a/docs/src/architectures/sympnet.md
+++ b/docs/src/architectures/sympnet.md
@@ -24,7 +24,7 @@ SympNets can be viewed as a *symplectic integrator* or symplectic one-step metho
 ```
 The tilde in the above equation indicates *predicted data*. With standard SympNets[^2] the time step between predictions is not a parameter we can choose but is related to the *temporal frequency of the training data*. This means that if data is recorded in an interval of e.g. 0.1 seconds, then this will be the time step of our integrator.
 
-[^2]: Recently an approach [horn2025generalized](@cite) has been proposed that makes explicitly specifying the time step possible by viewing SympNets as a subclass of so-called "Generalized Hamiltonian Neural Networks".
+[^2]: Recently an approach [horn2025generalized](@cite) has been proposed that makes explicitly specifying the time step possible by viewing SympNets as a subclass of so-called "Generalized Hamiltonian Neural Networks". Hamiltonian neural networks are implemented in [HamiltonianNeuralNetworks.jl](https://juliagml.github.io/HamiltonianNeuralNetworks.jl/latest/).
 
 SympNets preserve symplecticity by exploiting the ``(q, p)`` structure of the system. This is visualized below:
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 95fe4ca0d..08c3aadfc 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -4,7 +4,7 @@ CurrentModule = GeometricMachineLearning
 
 # Geometric Machine Learning
 
-`GeometricMachineLearning` is a package for *structure-preserving scientific machine learning*. It contains models that can learn dynamical systems with geometric structure, such as Hamiltonian (symplectic) or Lagrangian (variational) systems.
+`GeometricMachineLearning` is a package for *structure-preserving scientific machine learning*. It contains models that can learn dynamical systems with geometric structure, such as Hamiltonian (symplectic) or divergence-free (volume-preserving) systems. Hamiltonian and Lagrangian neural networks are in [HamiltonianNeuralNetworks.jl](https://juliagml.github.io/HamiltonianNeuralNetworks.jl/latest/) and [LagrangianNeuralNetworks.jl](https://juliagml.github.io/LagrangianNeuralNetworks.jl/dev/).
 
 In that regard its aim is similar to traditional *geometric numerical integration* [hairer2006geometric, Kraus:2020:GeometricIntegrators](@cite) in that it models maps that share properties with the analytic solution of a differential equation:
 
@@ -20,7 +20,7 @@ In that regard its aim is similar to traditional *geometric numerical integratio
 
 ## Architectures
 
-Some of the neural network architectures in `GeometricMachineLearning` [brantner2023symplectic, brantner2025volume](@cite) have emerged in connection with developing this package[^1], other have existed before [jin2020sympnets, greydanus2019hamiltonian](@cite).
+Some of the neural network architectures in `GeometricMachineLearning` [brantner2023symplectic, brantner2025volume](@cite) have emerged in connection with developing this package[^1], others have existed before [jin2020sympnets](@cite).
 
 [^1]: The work on this software package was done in connection with a PhD thesis. You can read its [introduction](@ref "Introduction and Outline") and [conclusion](@ref "Conclusion") here.
 
diff --git a/docs/src/structure_preservation/structure_preserving_neural_networks.md b/docs/src/structure_preservation/structure_preserving_neural_networks.md
index bfa6daa84..012822eb4 100644
--- a/docs/src/structure_preservation/structure_preserving_neural_networks.md
+++ b/docs/src/structure_preservation/structure_preserving_neural_networks.md
@@ -52,7 +52,7 @@ Main.example(raw"We say that a neural network is **symplectic** if ``\mathcal{NN
 
 If we have ``m = n`` then we can use [SympNets](@ref "SympNet Architecture") to realize such architectures; SympNets furthermore are universal approximators for the set of canonical symplectic maps[^3] [jin2020sympnets](@cite). If ``m \neq n`` we can use [symplectic autoencoders](@ref "The Symplectic Autoencoder") to realize such an architecture. A different class of neural networks are *volume-preserving neural networks*:
 
-[^3]: Other neural network architectures that were developed with the same aim are *Hamiltonian neural networks* [greydanus2019hamiltonian](@cite), *Hénon nets* [burby2020fast](@cite) and *generalized Hamiltonian neural networks* [horn2025generalized](@cite). [celledoni2021structure](@cite) gives an overview over *structure preserving neural networks*.
+[^3]: Other neural network architectures that were developed with the same aim are *Hamiltonian neural networks* [greydanus2019hamiltonian](@cite), *Hénon nets* [burby2020fast](@cite) and *generalized Hamiltonian neural networks* [horn2025generalized](@cite); Hamiltonian neural networks are implemented in [HamiltonianNeuralNetworks.jl](https://juliagml.github.io/HamiltonianNeuralNetworks.jl/latest/). [celledoni2021structure](@cite) gives an overview over *structure preserving neural networks*.
 
 ```@eval
 Main.example(raw"We say that a neural network is **volume-preserving** if ``\mathcal{NN}_\Theta:\mathbb{R}^n\to\mathbb{R}^n`` is such that: 
diff --git a/docs/src/tutorials/hamiltonian_neural_network.md b/docs/src/tutorials/hamiltonian_neural_network.md
deleted file mode 100644
index 136af1c75..000000000
--- a/docs/src/tutorials/hamiltonian_neural_network.md
+++ /dev/null
@@ -1,90 +0,0 @@
-# [Hamiltonian Neural Network](@id hnn_tutorial)
-
-In this tutorial we build a [Hamiltonian neural network](@ref hnn_architecture). 
-
-## Training a HNN Based on VectorField Data
-
-We first train a HNN [based on vector field data](@ref "HNN Loss for Vector Field Data"):
-
-```@example hnn
-using GeometricMachineLearning # hide
-import Random # hide
-Random.seed!(1234) # hide
-
-𝕁 = PoissonTensor(2)
-vf(z) = 𝕁 * z
-domain = [[q, p] for q in -1:.1:1 for p in -1:.1:1]
-vf_data = vf.(domain)
-domain_matrix = hcat(domain...)
-vf_matrix = hcat(vf_data...)
-dl = DataLoader(domain_matrix, vf_matrix)
-nothing # hide
-```
-
-We then build the neural network:
-
-```@example hnn
-const intermediate_dim = 5
-hnn_arch = StandardHamiltonianArchitecture(2, intermediate_dim)
-hnn = NeuralNetwork(hnn_arch)
-nothing # hide
-```
-
-Next we define the loss function
-
-```@example hnn
-loss = HNNLoss(hnn_arch)
-nothing # hide
-```
-
-We can now train the network
-
-```@example hnn
-batch = Batch(10)
-n_epochs = 100
-o = TrainingOptimizer(hnn.params; algorithm = Adam())
-loss_array = train!(hnn, o, dl, batch, n_epochs; loss)
-using CairoMakie # hide
-lines(loss_array) # hide
-```
-
-!!! info
-    Usually we use [`Zygote`](https://github.com/FluxML/Zygote.jl) for computing derivatives in `GeometricMachineLearning`, but as the [`Zygote` documentation](https://fluxml.ai/Zygote.jl/dev/limitations/#Second-derivatives-1) itself points out: "Often using a different AD system over Zygote is a better solution [for computing second-order derivatives]." For this reason we compute the loss of the HNN with [`SymbolicNeuralNetworks`](https://github.com/JuliaGNI/SymbolicNeuralNetworks.jl) and optionally also its gradient.
-
-## Training a HNN Based on Phase Space Data
-
-We now train a HNN on the same system [based on phase space data](@ref "HNN Loss for Phase Space Data"). The data are not vector fields, but pairs of points a fixed timestep apart. We produce such pairs by applying the exact flow of `vf`, which is ``\exp(\Delta{}t\mathbb{J})``, to the points of the domain:
-
-```@example hnn
-const Δt = .1
-next_matrix = exp(Δt * Matrix(𝕁)) * domain_matrix
-dl_pairs = DataLoader(domain_matrix, next_matrix)
-nothing # hide
-```
-
-We start from a new network of the architecture we built above:
-
-```@example hnn
-hnn_pairs = NeuralNetwork(hnn_arch)
-nothing # hide
-```
-
-The loss needs the timestep, because the finite difference it compares the vector field against does:
-
-```@example hnn
-loss_pairs = SymplecticEulerLoss(hnn_arch, Δt)
-nothing # hide
-```
-
-[`GeometricMachineLearning.SymplecticEulerLoss`](@ref) defaults to the `:A` variant, which evaluates the vector field at ``(q^{(t+1)}, p^{(t)})``. Passing `variant = :B` evaluates it at ``(q^{(t)}, p^{(t+1)})`` instead, which is the variant the [loss formula](@ref "HNN Loss for Phase Space Data") is written for. That formula states the evaluation points of the `:B` variant; the warning beside it records that its sign differs from the implementation.
-
-!!! info
-    The data comes from the exact flow, but the loss holds the network to one step of the symplectic Euler method. The two agree only to first order in ``\Delta{}t``. The target of this training is therefore not the Hamiltonian of `vf`, which is ``H(q, p) = (q^2 + p^2) / 2``, but the Hamiltonian whose symplectic Euler step reproduces the exact flow; the two differ at order ``\Delta{}t``. So ``H`` itself does not make this loss vanish, and `scripts/verification/symplectic_euler_loss_modified_hamiltonian.jl` shows the residual it leaves halves when ``\Delta{}t`` does. The first section compares the network against the vector field directly and has no such offset.
-
-We can now train the network:
-
-```@example hnn
-o_pairs = TrainingOptimizer(hnn_pairs.params; algorithm = Adam())
-loss_array_pairs = train!(hnn_pairs, o_pairs, dl_pairs, batch, n_epochs; loss = loss_pairs)
-lines(loss_array_pairs) # hide
-```
diff --git a/ext/GPUArraysCoreExt.jl b/ext/GPUArraysCoreExt.jl
index cf18844d1..972a678af 100644
--- a/ext/GPUArraysCoreExt.jl
+++ b/ext/GPUArraysCoreExt.jl
@@ -2,8 +2,9 @@ module GPUArraysCoreExt
 
 # `PoissonTensor`'s three array methods of `Base.:*` take a `Strided…` right-hand side. That bound
 # is what keeps them from colliding with every `*(::AbstractMatrix, ::X)` that ArrayLayouts,
-# FillArrays, Symbolics and GeometricOptimizers define for their own array type `X` -- seventeen
-# ambiguities, and the reason for the bound is set out in `src/arrays/poisson_tensor.jl`.
+# FillArrays, Symbolics (which HamiltonianNeuralNetworks.jl loads) and GeometricOptimizers define
+# for their own array type `X` -- seventeen ambiguities, and the reason for the bound is set out in
+# `src/arrays/poisson_tensor.jl`.
 #
 # A GPU array is itself strided, because `CuArray` and `MtlArray` are `DenseArray`s, so the fast
 # path already covers it. A *wrapped* one is not: `view(A, [1, 2, 3, 4], :)` and `A'` are a
diff --git a/scripts/reproduction/hnn_pendulum.jl b/scripts/reproduction/hnn_pendulum.jl
deleted file mode 100644
index 7f8faee83..000000000
--- a/scripts/reproduction/hnn_pendulum.jl
+++ /dev/null
@@ -1,45 +0,0 @@
-# Train a Hamiltonian neural network on the pendulum and plot what it learned.
-using GeometricMachineLearning
-
-# this contains the Hamiltonian and the functions for generating the training data
-include("../utilities/pendulum.jl")
-
-# `smoke_size(full, smoke)` returns the second when GML_SMOKE is set. See scripts/README.md.
-include("../utilities/smoke.jl")
-
-# layer dimension/width
-const ld = 5
-
-# hidden layers
-const ln = 3
-
-# number of inputs/dimension of system
-const ninput = 2
-
-# number of epochs
-const nepochs = smoke_size(2000, 2)
-
-# `StandardHamiltonianArchitecture` is named explicitly, because calling the abstract
-# `HamiltonianArchitecture` warns and defaults to it. The constructor takes
-# `(dim, width, nhidden, activation)` positionally.
-const hnn = StandardHamiltonianArchitecture(ninput, ld, ln)
-
-nn = NeuralNetwork(hnn, CPU(), Float64)
-
-# A grid of phase-space points and the symplectic gradient at each.
-input, output = get_data_set()
-dl = DataLoader(input, output)
-
-# `HNNLoss` asks that the Hamiltonian vector field of the learned Hamiltonian match the vector field
-# in the data. It is the loss `NetworkLoss(hnn)` returns, so naming it is for the reader.
-loss = HNNLoss(hnn)
-
-# The step size is the `TrainingOptimizer`'s `linesearch`, not the method's: the same `Adam()` trains
-# at any learning rate. `1e-3` is already Adam's default here and is named for the reader.
-optimizer = TrainingOptimizer(nn.params; algorithm = Adam(), linesearch = 1e-3)
-total_loss = train!(
-    nn, optimizer, dl, Batch(10), nepochs; loss = loss, show_progress = true)
-
-# plot results
-include("../utilities/plots.jl")
-plot_hnn(H, x -> only(nn(x)), total_loss; filename = "hnn_pendulum.png")
diff --git a/scripts/reproduction/hnn_pendulum_simulation.jl b/scripts/reproduction/hnn_pendulum_simulation.jl
deleted file mode 100644
index 2fd637503..000000000
--- a/scripts/reproduction/hnn_pendulum_simulation.jl
+++ /dev/null
@@ -1,59 +0,0 @@
-# Train a Hamiltonian neural network on the pendulum, then *integrate* the vector field it learned
-# and compare the trajectory and the energy drift against the true Hamiltonian.
-#
-# This is the half of the HNN story that `hnn_pendulum.jl` does not tell. That script plots what the
-# learned Hamiltonian looks like on a grid; this one asks what happens when you step with it. The
-# two questions come apart: a network can match `H` closely on the training grid and still drift in
-# energy over a long integration, and no test in the suite checks the second.
-using GeometricMachineLearning
-using GeometricIntegrators: ODEProblem, ImplicitMidpoint, integrate
-
-# the Hamiltonian, its symplectic gradient, and the training data
-include("../utilities/pendulum.jl")
-
-# `smoke_size(full, smoke)` returns the second when GML_SMOKE is set. See scripts/README.md.
-include("../utilities/smoke.jl")
-
-const ld = 5
-const ln = 3
-const ninput = 2
-const nepochs = smoke_size(2000, 2)
-
-# The integration the comparison is about. `n_steps` stays the same in a smoke run -- stepping is
-# cheap next to training, and a handful of steps would show no drift to compare.
-const timestep = 0.1
-const n_steps = 100
-const x₀ = [0.0, 1.0]
-
-const hnn = StandardHamiltonianArchitecture(ninput, ld, ln)
-nn = NeuralNetwork(hnn, CPU(), Float64)
-
-input, output = get_data_set()
-dl = DataLoader(input, output)
-
-# `HNNLoss` holds the Hamiltonian vector field it compares against the data. Naming the loss keeps
-# that field for the integration below, where constructing a second `HNNLoss` would build the
-# symbolic derivative again.
-const loss = HNNLoss(hnn)
-
-optimizer = TrainingOptimizer(nn.params; algorithm = Adam(), linesearch = 1e-3)
-total_loss = train!(
-    nn, optimizer, dl, Batch(10), nepochs; loss = loss, show_progress = true)
-
-# The learned Hamiltonian and the vector field built from it. `loss.hvf` is what `HNNLoss` compares
-# against the data, so this integrates exactly the quantity that was trained -- the same object,
-# not a second derivative taken again here.
-H̃(x) = only(nn(x))
-dH̃(x) = loss.hvf(x, nn.params)
-
-# `ODEProblem` calls its vector field as `v(v, t, q, params)`.
-reference_field(v, t, q, params) = (v .= dH(q))
-learned_field(v, t, q, params) = (v .= dH̃(q))
-
-const timespan = (0.0, n_steps * timestep)
-sol_ref = integrate(ODEProblem(reference_field, timespan, timestep, x₀), ImplicitMidpoint())
-sol_hnn = integrate(ODEProblem(learned_field, timespan, timestep, x₀), ImplicitMidpoint())
-
-include("../utilities/plots.jl")
-plot_network_sim(
-    H, H̃, sol_ref, sol_hnn, total_loss; filename = "hnn_pendulum_simulation.png")
diff --git a/scripts/reproduction/lnn_pendulum.jl b/scripts/reproduction/lnn_pendulum.jl
deleted file mode 100644
index e0abd6080..000000000
--- a/scripts/reproduction/lnn_pendulum.jl
+++ /dev/null
@@ -1,57 +0,0 @@
-# Train a Lagrangian neural network on the pendulum.
-#
-# The data are `(q, q̇)` on a grid and the acceleration the *exact* Lagrangian implies at each point.
-# `LNNLoss` solves the Euler-Lagrange equations of the *learned* Lagrangian for its own acceleration
-# and compares the two, so a decreasing loss means the network is recovering a Lagrangian whose
-# dynamics match the pendulum's.
-using GeometricMachineLearning
-
-# `smoke_size(full, smoke)` returns the second when GML_SMOKE is set, which is how the CI
-# job runs this script end to end in seconds. See scripts/README.md.
-include("../utilities/smoke.jl")
-
-# this contains the Lagrangians and the functions for generating the training data
-include("../utilities/data_problem.jl")
-
-const problem = :pendulum
-
-const _, n_dim = dict_problem_L[problem]
-
-# number of inputs/dimension of system
-const ninput = 2 * n_dim
-
-# Layer dimension/width, and the number of hidden layers. These are the sizes that cost time here,
-# and `nepochs` is not: `LNNLoss` reaches the Euler-Lagrange acceleration through a *symbolic*
-# gradient and Hessian of the chain, and the run's wall clock is the one-time compilation of their
-# pullback, which grows with the chain and falls in the first epoch -- in a smoke run 106s at
-# `(5, 3)`, 49s at `(5, 1)` and 21s at `(2, 1)`, against 1.8s for a second epoch. The smoke sizes
-# below are `LagrangianNeuralNetwork`'s own defaults, `width = dimin` and `nhidden = 1`, so a
-# smoke run trains through the same symbolic path as the full size.
-const ld = smoke_size(5, ninput)
-const ln = smoke_size(3, 1)
-
-# number of epochs
-const nepochs = smoke_size(200, 2)
-
-const arch = LagrangianNeuralNetwork(ninput; nhidden = ln, width = ld)
-
-nn = NeuralNetwork(arch, CPU(), Float64)
-
-input, output = get_LNN_data(problem)
-dl = DataLoader(input, output)
-
-# `LNNLoss(arch)` is what `NetworkLoss(arch)` returns; naming it is for the reader.
-loss = LNNLoss(arch)
-
-# The step size is the `TrainingOptimizer`'s `linesearch`, not the method's. `LNNLoss` needs a larger
-# one than `HNNLoss` does: its gradient reaches the parameters through a solve against the velocity
-# Hessian, where `HNNLoss` is linear in the gradient of the learned Hamiltonian.
-optimizer = TrainingOptimizer(nn.params; algorithm = Adam(), linesearch = 1e-2)
-total_loss = train!(
-    nn, optimizer, dl, Batch(10), nepochs; loss = loss, show_progress = true)
-
-# plot results
-include("../utilities/plots.jl")
-L, _ = dict_problem_L[problem]
-plot_hnn(L, x -> only(nn(x)), total_loss;
-    filename = "lnn_pendulum.png", xmin = -1.2, xmax = +1.2, ymin = -1.2, ymax = +1.2)
diff --git a/scripts/runscripts.jl b/scripts/runscripts.jl
index 53242c7a2..4786ad7e5 100644
--- a/scripts/runscripts.jl
+++ b/scripts/runscripts.jl
@@ -39,7 +39,7 @@ const PROJECT = Base.active_project()
 const TIMEOUT_SECONDS = Dict("verification" => 900, "reproduction" => 300)
 
 # The budget for a whole mode, inside `Scripts.yml`'s `timeout-minutes: 90` less what that job
-# spends instantiating. Without it the per-script ceiling does not keep the promise above: 23
+# spends instantiating. Without it the per-script ceiling does not keep the promise above: 20
 # reproduction scripts each entitled to the ceiling outlast any runner, and the job then dies at the
 # runner's timeout with no verdict at all. The driver stops first, and names what it did not reach.
 const BUDGET_SECONDS = 3600
diff --git a/scripts/utilities/data_problem.jl b/scripts/utilities/data_problem.jl
deleted file mode 100644
index 43586c4d1..000000000
--- a/scripts/utilities/data_problem.jl
+++ /dev/null
@@ -1,290 +0,0 @@
-using Zygote
-#using  GeometricIntegrators
-using LinearAlgebra
-using GeometricMachineLearning
-
-#################################################################################################
-# Hamiltomian Dictionary
-
-dict_problem_H = Dict(
-    :pendulum => (x -> x[2]^2 / 2 + (1-cos(x[1])), 1),
-    :Hénon_Heiles =>
-        (x -> x[3]^2/(2) + x[4]^2/(2) + x[1]^2/2 + x[2]^2/2 + (x[1]^2*x[2]-x[2]^3/3), 2)
-    #:three_body_system =>(,)
-)
-
-# Lagrangian DictionaryHH
-
-dict_problem_L = Dict(
-    :pendulum => (x -> x[2]^2 / 2 - (1-cos(x[1])), 1),
-)
-
-###########################################################################
-# get data of a problem for HNN (results is (data = (q,p), target = (q̇,ṗ)))
-
-function get_HNN_data(nameproblem, num = 10, qpmin = -1.2, qpmax = +1.2)
-
-    #get the Hamiltonien corresponding to name_problem
-    H, n_dim = dict_problem_H[nameproblem]
-
-    # compute vector field
-    ∇H(x) = Zygote.gradient(H, x)[1]
-
-    I = Diagonal(ones(n_dim))
-    Z = zeros(n_dim, n_dim)
-    symplectic_matrix = [Z I; -I Z]
-
-    dH(x) = symplectic_matrix * ∇H(x)
-
-    #range in which the data should be in
-    rang = range(qpmin, stop = qpmax, length = num)
-
-    pre_data = collect(Iterators.product(fill(rang, 2*n_dim)...))
-    data = [[x...] for x in pre_data]
-
-    #compute the value of the vector field 
-    target = dH.(data)
-
-    return (data, target)
-end
-
-###############################################################################
-# get data of a problem for LNN (results is (data = (q,q̇), target = (qdotdot)))
-
-function get_LNN_data(nameproblem, num = 10, qq̇min = -1.2, qq̇max = +1.2)
-
-    #get the Hamiltonien corresponding to name_problem
-    L, n_dim = dict_problem_L[nameproblem]
-
-    #compute gradient and hessian
-    ∇qL(x) = Zygote.gradient(x->L(x), x)[1][1:(length(x) ÷ 2)]
-    ∇∇L(x) = Zygote.hessian(L, x)
-    ∇q̇∇q̇L(x) = ∇∇L(x)[(1 + length(x) ÷ 2):end, (1 + length(x) ÷ 2):end]
-    ∇q∇q̇L(x) = ∇∇L(x)[1:(length(x) ÷ 2), (1 + length(x) ÷ 2):end]
-
-    # Solve the Euler-Lagrange equations for the acceleration. The transpose is load-bearing:
-    # `∇q∇q̇L` is indexed `[i, j] = ∂²L/∂qᵢ∂q̇ⱼ` and the chain rule contracts the *first* index with
-    # `q̇`. It vanishes for the pendulum, whose position and velocity do not couple, and does not for
-    # a Lagrangian that couples them. `GeometricMachineLearning.LNNLoss` solves the same equation.
-    Qdotdot(q, q̇) = inv(∇q̇∇q̇L([q..., q̇...])) *
-                    (∇qL([q..., q̇...]) - ∇q∇q̇L([q..., q̇...])' * q̇)
-
-    #range in which the data should be in
-    rang = range(qq̇min, stop = qq̇max, length = num)
-
-    pre_data = collect(Iterators.product(fill(rang, 2*n_dim)...))
-    data = [[[x[1:n_dim]...], [x[(1 + n_dim):(2 * n_dim)]...]] for x in pre_data]
-
-    #compute target associate to data
-    target = [Qdotdot(qq̇[1], qq̇[2]) for qq̇ in data]
-
-    # `(input, output)` for `DataLoader`: `input` stacks q on q̇ and is `2 n_dim × N`, `output` is
-    # the acceleration and is `n_dim × N`.
-    (reduce(hcat, [vcat(qq̇[1], qq̇[2]) for qq̇ in vec(data)]), reduce(hcat, vec(target)))
-end
-
-########################################################################################
-# get data of a problem for Sypmnet (results is (data = target =(q,p))) from Hamiltonian
-
-function get_phase_space_data(nameproblem, q₀, p₀, timespan = (0.0, 100.0), timestep = 0.1)
-
-    # get the Hamiltonien corresponding to name_problem   
-    H_problem, n_dim = dict_problem_H[nameproblem]
-
-    q, p = compute_phase_space(H_problem, q₀, p₀, timespan, timestep)
-
-    return (q, p)
-end
-
-function get_phase_space_multiple_trajectoy(
-        nameproblem; singlematrix = true, n_trajectory = 1, n_points = 10,
-        timestep = 0.1, qmin = -0.2, pmin = -0.2, qmax = 0.2, pmax = 0.2)
-
-    # get the Hamiltonien corresponding to name_problem   
-    H_problem, n_dim = dict_problem_H[nameproblem]
-
-    #define timespan
-    timespan=(0.0, n_points*timestep)
-
-    #compute phase space for each trajectory staring from a random point
-    trajectory_q = [zeros(n_points+1, n_dim) for _ in 1:n_trajectory]
-    trajectory_p = [zeros(n_points+1, n_dim) for _ in 1:n_trajectory]
-    for i in 1:n_trajectory
-        q₀ = [rand()*(qmax-qmin)+qmin for _ in 1:n_dim]
-        p₀ = [rand()*(pmax-pmin)+pmin for _ in 1:n_dim]
-        trajectory_q[i], trajectory_p[i] = compute_phase_space(
-            H_problem, q₀, p₀, timespan, timestep)
-    end
-
-    if singlematrix
-        data_q = vcat(trajectory_q...)
-        data_p = vcat(trajectory_p...)
-        return data_q, data_p
-    else
-        return trajectory_q, trajectory_p
-    end
-end
-
-###############################################################################
-# compute phase space from the Hamiltonian
-
-function compute_phase_space(H_problem, q₀, p₀, timespan = (0.0, 100.0), timestep = 0.1)
-    n_dim = length(q₀)
-
-    H2(q, p) = H_problem([q..., p...])
-
-    gradH(q, p) = gradient(H2, q, p)
-
-    function v(v, t, q, p, params)
-        v .= gradH(q, p)[2]
-    end
-
-    function f(f, t, q, p, params)
-        f .= - gradH(q, p)[1]
-    end
-
-    h(t, q, p, params) = H2(q, p)
-
-    # simulate data with geometric Integrators
-    ode = HODEProblem(v, f, h, timespan, timestep, q₀, p₀)
-
-    #return sol = integrate(ode, SymplecticEulerA())
-    return sol = integrate(ode, SymplecticTableau(TableauExplicitEuler()))
-
-    # results are give into a matrix where q[i,j] is the ith step of the jth component of q 
-    q = hcat([sol.q[:, i] for i in 1:n_dim]...)
-    p = hcat([sol.p[:, i] for i in 1:n_dim]...)
-
-    return (q, p)
-end
-
-###############################################################################
-# get structure data for multiple trajectory
-
-struct storing_data{T}
-    Δt::Float64
-    nb_trajectory::Int
-    data::T
-end
-
-function get_multiple_trajectory_structure(
-        nameproblem; n_trajectory = 1, n_points = 10, timestep = 0.1,
-        qmin = -1.2, pmin = -1.2, qmax = 1.2, pmax = 1.2)
-
-    # get the Hamiltonien corresponding to name_problem   
-    H_problem, n_dim = dict_problem_H[nameproblem]
-
-    #define timespan
-    timespan=(0.0, n_points*timestep)
-
-    #compute phase space for each trajectory staring from a random point
-    pre_data = NamedTuple()
-
-    for i in 1:n_trajectory
-        q₀ = [rand()*(qmax-qmin)+qmin for _ in 1:n_dim]
-        p₀ = [rand()*(pmax-pmin)+pmin for _ in 1:n_dim]
-        q, p = compute_phase_space(H_problem, q₀, p₀, timespan, timestep)
-
-        Data = [(q[n], p[n]) for n in 1:size(q, 1)]
-
-        nt = NamedTuple{(Symbol("Trajectory_"*string(i)),)}(((
-            data = Data, len = n_points+1),))
-
-        pre_data = merge(pre_data, nt)
-    end
-
-    data = storing_data(timestep, n_trajectory, pre_data)
-
-    return data #data_trajectory(data, Get_nb_trajectory, Get_length_trajectory, Get_q, Get_p, Get_Δt)
-end
-
-function get_multiple_trajectory_structure_with_target(
-        nameproblem; n_trajectory = 1, n_points = 10, timestep = 0.1,
-        qmin = -0.2, pmin = -0.2, qmax = 0.2, pmax = 0.2)
-
-    # get the Hamiltonien corresponding to name_problem   
-    H_problem, n_dim = dict_problem_H[nameproblem]
-
-    # compute vector field
-    ∇H(x) = Zygote.gradient(H, x)[1]
-
-    I = Diagonal(ones(n_dim))
-    Z = zeros(n_dim, n_dim)
-    symplectic_matrix = [Z I; -I Z]
-
-    dH(x) = symplectic_matrix * ∇H(x)
-
-    #define timespan
-    timespan=(0.0, n_points*timestep)
-
-    #compute phase space for each trajectory staring from a random point
-    pre_data = NamedTuple()
-    Target = NamedTuple()
-
-    for i in 1:n_trajectory
-        q₀ = [rand()*(qmax-qmin)+qmin for _ in 1:n_dim]
-        p₀ = [rand()*(pmax-pmin)+pmin for _ in 1:n_dim]
-        q, p = compute_phase_space(H_problem, q₀, p₀, timespan, timestep)
-
-        Data = [(q[n], p[n]) for n in 1:size(q, 1)]
-        data_calc = [[q[n]..., p[n]...] for n in 1:size(q, 1)]
-
-        nt = NamedTuple{(Symbol("Trajectory_"*string(i)),)}(((
-            data = Data, len = n_points+1),))
-
-        pre_data = merge(pre_data, nt)
-
-        calc = dH.(data_calc)
-        pre_target = [(calc[n][1:n_dim], calc[n][(1 + n_dim):end]) for n in 1:size(q, 1)]
-        ntt = NamedTuple{(Symbol("Trajectory_"*string(i)),)}(((
-            target = pre_target, len = n_points+1),))
-        Target = merge(Target, ntt)
-    end
-
-    data = storing_data(timestep, n_trajectory, pre_data)
-
-    Get_Δt(Data) = Data.Δt
-    Get_nb_trajectory(Data) = Data.nb_trajectory
-    Get_length_trajectory(Data, i) = Data.data[Symbol("Trajectory_"*string(i))][:len]
-    Get_q(Data, i, n) = Data.data[Symbol("Trajectory_"*string(i))][:data][n][1]
-    Get_p(Data, i, n) = Data.data[Symbol("Trajectory_"*string(i))][:data][n][2]
-
-    Get_q̇(Target, i, n) = Target[Symbol("Trajectory_"*string(i))][:target][n][1][1]
-    Get_ṗ(Target, i, n) = Target[Symbol("Trajectory_"*string(i))][:target][n][2][1]
-
-    return dataTarget(
-        data_trajectory(
-            data, Get_nb_trajectory, Get_length_trajectory, Get_q, Get_p, Get_Δt),
-        Target,
-        Get_q̇,
-        Get_ṗ)
-end
-
-function get_multiple_trajectory_structure_Lagrangian(
-        nameproblem; n_trajectory = 1, n_points = 10, timestep = 0.1,
-        qmin = -0.2, pmin = -0.2, qmax = 0.2, pmax = 0.2)
-
-    # get the Hamiltonien corresponding to name_problem   
-    H_problem, n_dim = dict_problem_L[nameproblem]
-
-    #define timespan
-    timespan=(0.0, n_points*timestep)
-
-    #compute phase space for each trajectory staring from a random point
-    pre_data = []
-    push!(pre_data, [timestep])
-    push!(pre_data, [n_trajectory])
-    push!(pre_data, [n_points+1])
-
-    for i in 1:n_trajectory
-        q₀ = [rand()*(qmax-qmin)+qmin for _ in 1:n_dim]
-        p₀ = [rand()*(pmax-pmin)+pmin for _ in 1:n_dim]
-        q, p = compute_phase_space(H_problem, q₀, p₀, timespan, timestep)
-
-        Data = [q[n] for n in 1:size(q, 1)]
-
-        pre_data = push!(pre_data, Data)
-    end
-
-    return pre_data
-end
diff --git a/scripts/utilities/pendulum.jl b/scripts/utilities/pendulum.jl
index 862794ae6..883300d2e 100644
--- a/scripts/utilities/pendulum.jl
+++ b/scripts/utilities/pendulum.jl
@@ -5,11 +5,6 @@ H(x) = x[2]^2 / 2 + (1-cos(x[1]))
 H(q, p) = H([q[1], p[1]])
 H(t, q, p, params) = H(q, p)
 
-# The gradient of `H` and the symplectic gradient built from it. Both are written out rather than
-# taken with automatic differentiation, so that building the training data below needs no AD.
-∇H(x) = [sin(x[1]), x[2]]
-dH(x) = [0 1; -1 0] * ∇H(x)
-
 # vector field methods
 function v(v, t, q, p, params)
     v[1] = p[1]
@@ -18,28 +13,6 @@ function f(f, t, q, p, params)
     f[1] = -sin(q[1])
 end
 
-"""
-    get_data_set(num, xymin, xymax)
-
-A grid of `num`² points in phase space, together with the symplectic gradient at each — the
-`(q, p)` and `(q̇, ṗ)` a Hamiltonian neural network trains on.
-
-Returns `(input, output)`, both `2 × num²` matrices, which is what `DataLoader(input, output)` takes
-and what `HNNLoss` compares.
-"""
-function get_data_set(num = 10, xymin = -1.2, xymax = +1.2)
-    #range in which the data should be in
-    rang = range(xymin, stop = xymax, length = num)
-
-    # all combinations of (x,y) points
-    points = [[x, y] for x in rang for y in rang]
-
-    #compute the value of the vector field
-    derivatives = dH.(points)
-
-    (reduce(hcat, points), reduce(hcat, derivatives))
-end
-
 @doc raw"""
 Generates data for a pendulum in 2d with optional arguments:
 - `T`: the type of the data (`Float32`, `Float64`, `Float16`, etc.)
diff --git a/scripts/utilities/plots.jl b/scripts/utilities/plots.jl
deleted file mode 100644
index 150899d45..000000000
--- a/scripts/utilities/plots.jl
+++ /dev/null
@@ -1,113 +0,0 @@
-using CairoMakie
-using LaTeXStrings
-
-# Fill colours for the contour panels: the Hamiltonian panels use a reversed viridis and the error
-# panel the unreversed one, so that the two read differently at a glance.
-const _reversed_colormap = Reverse(:viridis)
-const _default_colormap = :viridis
-
-"""
-    _contour_panel!(position, X, Y, Z; kwargs...)
-
-A filled contour of `Z` over `X × Y`, on its own axis at `position` in the enclosing figure layout.
-Returns the axis so a caller can draw a second, unfilled contour on top of it.
-"""
-function _contour_panel!(
-        position, X, Y, Z; levels = 7, colormap = _reversed_colormap, kwargs...)
-    ax = Axis(position; kwargs...)
-    contourf!(ax, X, Y, Z; levels = levels, colormap = colormap)
-    ax
-end
-
-function plot_hnn(H, H̃, total_loss; xmin = -1.2, xmax = +1.2, ymin = -1.2,
-        ymax = +1.2, nsamples = 100, filename = nothing)
-    #get offset of learned Hamiltonian
-    H̃₀ = H̃([0, 0])
-
-    X = range(xmin, stop = xmax, length = nsamples)
-    Y = range(ymin, stop = ymax, length = nsamples)
-
-    fig = Figure(size = (1000, 800))
-
-    # contour lines of the Hamiltonian, with those of the learned one on top
-    ax_cnt = _contour_panel!(fig[1, 1], X, Y, [H([x, y]) for x in X, y in Y];
-        xlabel = L"q", ylabel = L"p")
-    contour!(ax_cnt, X, Y, [H̃([x, y]) - H̃₀ for x in X, y in Y]; color = :black)
-
-    # contours of the error of the Hamiltonian, in per cent of its maximum
-    m = maximum(H([x, y]) for x in X, y in Y)
-    _contour_panel!(
-        fig[1, 2], X, Y, 100*[(H̃([x, y]) - H̃₀ - H([x, y]))/m for x in X, y in Y];
-        colormap = _default_colormap, xlabel = L"q", ylabel = L"p")
-
-    # total loss, across the full width below the two contour panels
-    ax_loss = Axis(fig[2, 1:2]; xlabel = "n(training)", ylabel = "Total Loss")
-    lines!(ax_loss, total_loss)
-
-    # the loss occupies the bottom 40% of the figure
-    rowsize!(fig.layout, 2, Relative(0.4))
-
-    if filename !== nothing
-        CairoMakie.save(filename, fig)
-    end
-
-    return fig
-end
-
-"""
-    plot_network_sim(H, H̃, sol_ref, sol_hnn, total_loss; kwargs...)
-
-Four panels comparing a learned Hamiltonian against the true one: its contours, the two phase-space
-trajectories, the training loss, and the energy drift along each trajectory.
-
-`sol_ref` and `sol_hnn` are the solutions of the true and the learned Hamiltonian vector field over
-the same time span. The caller integrates them, so this function needs no solver: the integration
-belongs to the script that chooses the method and the step size.
-"""
-function plot_network_sim(H, H̃, sol_ref, sol_hnn, total_loss; xmin = -1.2, xmax = +1.2,
-        ymin = -1.2, ymax = +1.2, nsamples = 100, filename = nothing)
-    # get offset of learned Hamiltonian
-    H̃₀ = H̃([0, 0])
-
-    X = range(xmin, stop = xmax, length = nsamples)
-    Y = range(ymin, stop = ymax, length = nsamples)
-
-    fig = Figure(size = (1000, 800))
-
-    # contour lines of the Hamiltonian, with those of the learned one on top
-    ax_cnt = _contour_panel!(fig[1, 1], X, Y, [H([x, y]) for x in X, y in Y];
-        title = L"$H(q,p)$", xlabel = L"$q$", ylabel = L"$p$")
-    contour!(ax_cnt, X, Y, [H̃([x, y]) - H̃₀ for x in X, y in Y]; color = :black)
-
-    # solutions. A `DataSeries` is indexed from 0 and `sol.q[:, k]` is the series of component `k`,
-    # so both axes are collected into plain vectors before they reach Makie.
-    ax_sim = Axis(fig[1, 2]; xlabel = L"$q$", ylabel = L"$p$")
-    lines!(ax_sim, collect(sol_ref.q[:, 1]), collect(sol_ref.q[:, 2]); label = "Reference")
-    lines!(ax_sim, collect(sol_hnn.q[:, 1]), collect(sol_hnn.q[:, 2]); label = "HNN")
-    axislegend(ax_sim)
-
-    # total loss
-    ax_loss = Axis(fig[2, 1]; xlabel = L"$n_{\mathrm{training}}$", ylabel = "Total Loss",
-        yscale = log10)
-    lines!(ax_loss, total_loss)
-
-    # Energy drift along each trajectory, each against its own initial value. The learned
-    # Hamiltonian is only determined up to a constant, so its drift is what is comparable, not its
-    # value.
-    # `parent` is not decoration: a comprehension over a `DataSeries` keeps its 0-based axes, where
-    # `collect(sol.t)` returns a 1-based vector, and Makie broadcasts the two against each other.
-    drift(f, sol) = parent([f(q) for q in sol.q]) .- f(sol.q[begin])
-    ax_err = Axis(fig[2, 2]; xlabel = L"$t$", ylabel = L"$\Delta H(q(t))$")
-    lines!(ax_err, collect(sol_ref.t), drift(H, sol_ref); label = "Reference")
-    lines!(ax_err, collect(sol_hnn.t), drift(H̃, sol_hnn); label = "HNN")
-    axislegend(ax_err)
-
-    # the contour and trajectory panels occupy the top 70% of the figure
-    rowsize!(fig.layout, 1, Relative(0.7))
-
-    if filename !== nothing
-        CairoMakie.save(filename, fig)
-    end
-
-    return fig
-end
diff --git a/scripts/verification/symplectic_euler_loss_modified_hamiltonian.jl b/scripts/verification/symplectic_euler_loss_modified_hamiltonian.jl
deleted file mode 100644
index e45e58f9f..000000000
--- a/scripts/verification/symplectic_euler_loss_modified_hamiltonian.jl
+++ /dev/null
@@ -1,66 +0,0 @@
-# Check that the true Hamiltonian does not minimise SymplecticEulerLoss on exact-flow data, and
-# that the residual it leaves is first order in the timestep.
-#
-# Run from the repo root:
-#   julia --startup-file=no --project=. scripts/verification/symplectic_euler_loss_modified_hamiltonian.jl
-#
-# The HNN tutorial's phase space section trains on pairs carried by the exact flow of vf(z) = 𝕁z,
-# and scores the network with SymplecticEulerLoss, which holds it to one symplectic Euler step.
-# The two agree only to first order, so the minimiser of that loss is not H but the Hamiltonian
-# whose symplectic Euler step reproduces the exact flow. This is what stops a trained loss below
-# the true Hamiltonian's residual from reading as a fit better than exact. The tutorial's
-# `!!! info` at docs/src/tutorials/hamiltonian_neural_network.md records the rule.
-#
-# The loss is reimplemented here from src/loss/symplectic_euler_loss.jl rather than called: the
-# package form takes a network and its parameters, and the point of the check is the value at the
-# exact Hamiltonian, which is not a network. The two must not drift apart — the expression below
-# is that file's final three lines, with loss.hvf replaced by the field of H.
-
-using LinearAlgebra: norm
-
-const 𝕁 = [0.0 1.0; -1.0 0.0]
-
-# X_H = (∂H/∂p, -∂H/∂q) for H = (q² + p²)/2, the Hamiltonian of vf(z) = 𝕁z.
-hamiltonian_vector_field(z) = vcat(z[2:2, :], -z[1:1, :])
-
-"""
-The tutorial's grid, carried one step of `timestep` by the exact flow `exp(timestep * 𝕁)`.
-"""
-function exact_flow_data(timestep)
-    input = hcat([[q, p] for q in -1:0.1:1 for p in -1:0.1:1]...)
-    input, exp(timestep * 𝕁) * input
-end
-
-"""
-`SymplecticEulerLoss` evaluated at the true Hamiltonian instead of at a network.
-"""
-function residual(timestep, variant)
-    input, output = exact_flow_data(timestep)
-    n = size(input, 1) ÷ 2
-    evaluation_point = variant === :A ?
-                       vcat(output[1:n, :], input[(n + 1):(2n), :]) :
-                       vcat(input[1:n, :], output[(n + 1):(2n), :])
-    difference_quotient = (output - input) / timestep
-    norm(hamiltonian_vector_field(evaluation_point) - difference_quotient) /
-    norm(difference_quotient)
-end
-
-const TIMESTEPS = (0.1, 0.05, 0.025)
-
-for variant in (:A, :B)
-    residuals = [residual(timestep, variant) for timestep in TIMESTEPS]
-    ratios = residuals[1:(end - 1)] ./ residuals[2:end]
-
-    println("variant :$(variant)")
-    for (timestep, r) in zip(TIMESTEPS, residuals)
-        println("  Δt = $(timestep): residual = $(r)")
-    end
-    println("  ratios on halving Δt: $(ratios)")
-
-    # The true Hamiltonian leaves a residual, so it is not the minimiser.
-    @assert all(>(1e-3), residuals) "the true Hamiltonian minimises the loss for variant :$(variant)"
-    # First order: halving the timestep halves the residual.
-    @assert all(r -> isapprox(r, 2.0; atol = 1e-2), ratios) "residual is not first order in Δt for variant :$(variant)"
-end
-
-println("\nBoth variants: the true Hamiltonian leaves a first-order residual.")
diff --git a/src/GeometricMachineLearning.jl b/src/GeometricMachineLearning.jl
index 838601ba2..f596c49e3 100644
--- a/src/GeometricMachineLearning.jl
+++ b/src/GeometricMachineLearning.jl
@@ -6,13 +6,11 @@ using AbstractNeuralNetworks
 # package also exports `flatten`/`unflatten` and the leaf protocol, none of which this package
 # extends — `GeometricOptimizers` carries the protocol for the structured matrices.
 #
-# Two things come from there rather than being written again here: `mapstorage`, which reaches
-# the storage of a structured leaf and rebuilds the leaf around the result (`src/map_to_cpu.jl`), and
-# `parameter_eltype`, which gives the element type a parameter set holds — what a seeded gradient
-# is converted to. A plain `NamedTuple` of layers is walked with `Base.map`, which needs nothing
-# from anybody.
+# `mapstorage` comes from there rather than being written again here: it reaches the storage of a
+# structured leaf and rebuilds the leaf around the result (`src/map_to_cpu.jl`). A plain
+# `NamedTuple` of layers is walked with `Base.map`, which needs nothing from anybody.
 import NeuralNetworkParameters: NetworkParameters
-using NeuralNetworkParameters: mapstorage, parameter_eltype
+using NeuralNetworkParameters: mapstorage
 using ChainRulesCore
 using GeometricBase
 using GeometricSolutions: GeometricSolution, EnsembleSolution, TimeSeries
@@ -26,10 +24,6 @@ using Zygote
 using ForwardDiff
 using InteractiveUtils
 using TimerOutputs
-import SymbolicNeuralNetworks
-import SymbolicNeuralNetworks: SymbolicPullback
-using SymbolicNeuralNetworks: derivative, SymbolicNeuralNetwork
-import Symbolics
 
 # The manifolds, the structured matrix types, the global sections and the retractions are
 # `GeometricOptimizers`': the same objects, so upstream's generic machinery dispatches on them and
@@ -69,8 +63,6 @@ import GeometricOptimizers: solve!
 import AbstractNeuralNetworks: Architecture, AbstractExplicitLayer, NeuralNetwork,
                                FeedForwardLoss
 import AbstractNeuralNetworks: Chain
-# `input_dimension`/`output_dimension` are AbstractNeuralNetworks' since v0.6.4; the `Chain`
-# methods GML uses are added to them by SymbolicNeuralNetworks.
 import AbstractNeuralNetworks: input_dimension, output_dimension
 import AbstractNeuralNetworks: Dense, Linear
 # `update!` is deliberately not among these: `AbstractNeuralNetworks`' is a different generic
@@ -188,8 +180,7 @@ include("optimizers/optimizer.jl")
 
 export NeuralNetwork
 
-export NetworkLoss, TransformerLoss, FeedForwardLoss, AutoEncoderLoss, ReducedLoss, HNNLoss,
-       LNNLoss, SymplecticEulerLoss, VariationalMidpointLoss
+export NetworkLoss, TransformerLoss, FeedForwardLoss, AutoEncoderLoss, ReducedLoss
 
 #INCLUDE ARCHITECTURES
 include("architectures/neural_network_integrator.jl")
@@ -200,34 +191,23 @@ include("architectures/sympnet.jl")
 include("architectures/autoencoder.jl")
 include("architectures/symplectic_autoencoder.jl")
 include("architectures/psd.jl")
-include("architectures/hamiltonian_neural_network.jl")
-include("architectures/lagrangian_neural_network.jl")
 include("architectures/transformer_neural_network.jl")
 include("architectures/volume_preserving_feedforward.jl")
 include("architectures/volume_preserving_transformer.jl")
 
-export HamiltonianArchitecture
-export LagrangianNeuralNetwork
 export SympNet, LASympNet, GSympNet
 export ClassificationTransformer, ClassificationLayer
 export VolumePreservingFeedForward
 export SymplecticAutoencoder, PSDArch
-export HamiltonianArchitecture, StandardHamiltonianArchitecture,
-       GeneralizedHamiltonianArchitecture
 
 export solve!, encoder, decoder
 
 export iterate
 
 include("loss/losses.jl")
-include("loss/hnn_loss.jl")
-include("loss/lnn_loss.jl")
-include("loss/symplectic_euler_loss.jl")
-include("loss/variational_midpoint_loss.jl")
 
-export AbstractPullback, ZygotePullback, SymbolicPullback
+export AbstractPullback, ZygotePullback
 include("pullbacks/zygote_pullback.jl")
-include("pullbacks/symbolic_hnn_pullback.jl")
 
 export DataLoader
 export Batch, train!
diff --git a/src/architectures/hamiltonian_neural_network.jl b/src/architectures/hamiltonian_neural_network.jl
deleted file mode 100644
index b8e91c2a4..000000000
--- a/src/architectures/hamiltonian_neural_network.jl
+++ /dev/null
@@ -1,136 +0,0 @@
-"""
-    HamiltonianArchitecture <: Architecture
-
-See [`StandardHamiltonianArchitecture`](@ref) and [`GeneralizedHamiltonianArchitecture`](@ref).
-"""
-abstract type HamiltonianArchitecture{AT <: Activation} <: Architecture end
-
-const HNN_nhidden_default = 1
-const HNN_activation_default = AbstractNeuralNetworks.TanhActivation()
-
-function HamiltonianArchitecture(dim::Integer, width::Integer, nhidden::Integer, activation)
-    @warn "You called the abstract type `HamiltonianArchitecture` as a constructor. This is defaulting to `StandardHamiltonianArchitecture`."
-    StandardHamiltonianArchitecture(dim, width, nhidden, activation)
-end
-
-"""
-    StandardHamiltonianArchitecture <: HamiltonianArchitecture
-
-A realization of the standard Hamiltonian neural network (HNN) [greydanus2019hamiltonian](@cite).
-
-Also see [`GeneralizedHamiltonianArchitecture`](@ref).
-
-# Constructor
-
-The constructor takes the following input arguments:
-1. `dim`: system dimension,
-2. `width = dim`: width of the hidden layer. By default this is equal to `dim`,
-3. `nhidden = $(HNN_nhidden_default)`: the number of hidden layers,
-4. `activation = $(HNN_activation_default)`: the activation function used in the HNN.
-"""
-struct StandardHamiltonianArchitecture{AT} <: HamiltonianArchitecture{AT}
-    dim::Int
-    width::Int
-    nhidden::Int
-    activation::AT
-
-    function StandardHamiltonianArchitecture(
-            dim, width = dim, nhidden = HNN_nhidden_default,
-            activation = HNN_activation_default)
-        new{typeof(activation)}(dim, width, nhidden, activation)
-    end
-end
-
-GHNN_integrator_default = nothing
-
-"""
-    GeneralizedHamiltonianArchitecture <: HamiltonianArchitecture
-
-A realization of generalized Hamiltonian neural networks (GHNNs) as introduced in [horn2025generalized](@cite).
-
-Also see [`StandardHamiltonianArchitecture`](@ref).
-
-# Constructor
-
-The constructor takes the following input arguments:
-1. `dim`: system dimension,
-2. `width = dim`: width of the hidden layer. By default this is equal to `dim`,
-3. `nhidden = $(HNN_nhidden_default)`: the number of hidden layers,
-4. `activation = $(HNN_activation_default)`: the activation function used in the GHNN,
-5. `integrator = $(GHNN_integrator_default)`: the integrator that is used to design the GHNN.
-"""
-struct GeneralizedHamiltonianArchitecture{AT, IT} <: HamiltonianArchitecture{AT}
-    dim::Int
-    width::Int
-    nhidden::Int
-    activation::AT
-    integrator::IT
-
-    function GeneralizedHamiltonianArchitecture(
-            dim, width = dim, nhidden = HNN_nhidden_default,
-            activation = HNN_activation_default, integrator = GHNN_integrator_default)
-        error("GHNN still has to be implemented!")
-        # fatou-ignore unreachable-code
-        new{typeof(activation), typeof(integrator)}(
-            dim, width, nhidden, activation, integrator)
-    end
-end
-
-@inline AbstractNeuralNetworks.dim(arch::HamiltonianArchitecture) = arch.dim
-
-"""
-    symbolic_hamiltonian_vector_field(nn::SymbolicNeuralNetwork)
-
-Get the symbolic expression for the vector field belonging to the HNN `nn`.
-
-# Implementation 
-
-This is calling `SymbolicNeuralNetworks.Jacobian` and then multiplies the result with a Poisson tensor.
-"""
-function symbolic_hamiltonian_vector_field(nn::SymbolicNeuralNetwork)
-    □ = SymbolicNeuralNetworks.Jacobian(nn)
-    n = input_dimension(nn.model) ÷ 2
-    # The Poisson tensor is built from *integers* on purpose: a `Float64` literal in the symbolic
-    # expression would widen the result of a `Float32` network. `PoissonTensor`, which has the same
-    # convention, is an `AbstractMatrix{Float64}` and would do exactly that.
-    𝕆 = zeros(Int, n, n)
-    𝕀 = Matrix(1I, n, n)
-    𝕁 = [𝕆 𝕀; -𝕀 𝕆]
-    # `Jacobian` uses the convention `□[i, j] = ∂fᵢ/∂xⱼ` and the HNN output is scalar, so the one
-    # row of `derivative(□)` is the gradient of the Hamiltonian. The vector field is built as a
-    # *vector*, so that the generated function returns what `HNNLoss` compares against: a vector
-    # for a single sample, and one column per sample for a batch.
-    ∇H = vec(derivative(□))
-    𝕁 * ∇H
-end
-
-"""
-    hamiltonian_vector_field(arch::HamiltonianArchitecture)
-
-Compute an executable expression of the Hamiltonian vector field of a [`HamiltonianArchitecture`](@ref).
-
-# Implementation
-
-This first computes a symbolic expression of the vector field using [`symbolic_hamiltonian_vector_field`](@ref).
-
-The function is built with `inplace = false`: [`HNNLoss`](@ref) wraps it and is differentiated with
-`Zygote`, and the in-place kernel `SymbolicNeuralNetworks.build_nn_function` builds by default
-*mutates* its result, which `Zygote` does not support.
-"""
-function hamiltonian_vector_field(arch::HamiltonianArchitecture)
-    nn = SymbolicNeuralNetwork(arch)
-    hvf = symbolic_hamiltonian_vector_field(nn)
-    SymbolicNeuralNetworks.build_nn_function(hvf, nn.params, nn.input; inplace = false)
-end
-
-function Chain(arch::HamiltonianArchitecture)
-    inner_layers = Tuple(
-        [Dense(arch.width, arch.width, arch.activation) for _ in 1:arch.nhidden]
-    )
-
-    Chain(
-        Dense(arch.dim, arch.width, arch.activation),
-        inner_layers...,
-        Linear(arch.width, 1; use_bias = false)
-    )
-end
diff --git a/src/architectures/lagrangian_neural_network.jl b/src/architectures/lagrangian_neural_network.jl
deleted file mode 100644
index 96bc0c06b..000000000
--- a/src/architectures/lagrangian_neural_network.jl
+++ /dev/null
@@ -1,24 +0,0 @@
-struct LagrangianNeuralNetwork{AT} <: Architecture
-    dimin::Int
-    width::Int
-    nhidden::Int
-    act::AT
-
-    function LagrangianNeuralNetwork(dimin; width = dimin, nhidden = 1, activation = tanh)
-        new{typeof(activation)}(dimin, width, nhidden, activation)
-    end
-end
-
-@inline AbstractNeuralNetworks.dim(arch::LagrangianNeuralNetwork) = arch.dimin
-
-function Chain(nn::LagrangianNeuralNetwork)
-    inner_layers = Tuple(
-        [Dense(nn.width, nn.width, nn.act) for _ in 1:nn.nhidden]
-    )
-
-    Chain(
-        Dense(nn.dimin, nn.width, nn.act),
-        inner_layers...,
-        Linear(nn.width, 1; use_bias = false)
-    )
-end
diff --git a/src/arrays/poisson_tensor.jl b/src/arrays/poisson_tensor.jl
index aea25916d..0111559cb 100644
--- a/src/arrays/poisson_tensor.jl
+++ b/src/arrays/poisson_tensor.jl
@@ -102,12 +102,13 @@ end
 # The right-hand side of the three array methods is `Strided…`, not `Abstract…`, and the bound is
 # load-bearing rather than cosmetic. `PoissonTensor <: AbstractMatrix{T}`, so a method claiming the
 # whole of `AbstractVecOrMat` on the right collides with every `*(::AbstractMatrix, ::X)` that
-# ArrayLayouts, FillArrays, Symbolics and GeometricOptimizers define for their own special array
-# type `X` -- seventeen ambiguities, each of which throws a `MethodError` on a call a user can
-# write, `PoissonTensor(4, Float32) * StiefelManifold(...)` among them. `StridedArray` excludes
-# every such `X` while still covering what this package produces and consumes: `Array`, a strided
-# `SubArray`, and the GPU arrays, because `CuArray` and `MtlArray` are `DenseArray`s. An `Adjoint`
-# of a `Matrix` is not strided and takes the generic path below.
+# ArrayLayouts, FillArrays, Symbolics (which HamiltonianNeuralNetworks.jl loads) and
+# GeometricOptimizers define for their own special array type `X` -- seventeen ambiguities, each of
+# which throws a `MethodError` on a call a user can write, `PoissonTensor(4, Float32) *
+# StiefelManifold(...)` among them. `StridedArray` excludes every such `X` while still covering
+# what this package produces and consumes: `Array`, a strided `SubArray`, and the GPU arrays,
+# because `CuArray` and `MtlArray` are `DenseArray`s. An `Adjoint` of a `Matrix` is not strided and
+# takes the generic path below.
 #
 # A non-strided matrix or vector on the CPU reaches Julia's generic `AbstractMatrix` multiply, which
 # gives the same answer: this type carries `getindex` and `size` and *is* the matrix it claims to
diff --git a/src/data_loader/train.jl b/src/data_loader/train.jl
index 8ca5c2a22..2ced7f23a 100644
--- a/src/data_loader/train.jl
+++ b/src/data_loader/train.jl
@@ -18,7 +18,10 @@ trained through a [`ZygotePullback`](@ref). Without it, a [`TransformerIntegrato
 autoencoder with a `Batch{:FeedForward}` on [`AutoEncoderLoss`](@ref)`()`; any other combination
 needs `loss`, and `train!` raises an `ArgumentError` without it.
 
-A Hamiltonian network trains on its symbolic gradient with `loss = SymbolicPullback(arch)`.
+A network of HamiltonianNeuralNetworks.jl trains on vector-field data with
+`loss = SymbolicPullback(arch)` of that package, and on consecutive states with a loss such as its
+`SymplecticEulerLoss(arch, timestep)`; see its
+[manual](https://juliagml.github.io/HamiltonianNeuralNetworks.jl/latest/).
 
 `opt` is built for `params(nn)`, with `TrainingOptimizer(params(nn); algorithm, linesearch,
 retraction)`. It keeps the iteration number, the moments and the global sections, so a second call
diff --git a/src/loss/hnn_loss.jl b/src/loss/hnn_loss.jl
deleted file mode 100644
index f5a4fbea9..000000000
--- a/src/loss/hnn_loss.jl
+++ /dev/null
@@ -1,42 +0,0 @@
-"""
-    HNNLoss <: NetworkLoss
-
-The loss for a Hamiltonian neural network.
-
-# Constructor
-
-This can be called with a `NeuralNetwork`, built with a [`HamiltonianArchitecture`](@ref), as the only input arguemtn, i.e.:
-```julia
-HNNLoss(nn)
-```
-where `nn` is a `NeuralNetwork`, that is built with a [`HamiltonianArchitecture`](@ref), gives the corresponding Hamiltonian loss.
-
-# Functor
-
-```julia
-loss(c, ps, input, output)
-loss(ps, input, output) # equivalent to the above
-```
-"""
-struct HNNLoss{FT} <: NetworkLoss
-    hvf::FT
-end
-
-function HNNLoss(arch::HamiltonianArchitecture)
-    HNNLoss(hamiltonian_vector_field(arch))
-end
-
-AbstractNeuralNetworks.NetworkLoss(arch::HamiltonianArchitecture) = HNNLoss(arch)
-
-function (loss::HNNLoss)(::Union{Chain, AbstractExplicitLayer},
-        ps::NetworkParameters,
-        input::ArrayOrNamedTuple,
-        output::ArrayOrNamedTuple)
-    loss(ps, input, output)
-end
-
-function (loss::HNNLoss)(ps::NetworkParameters,
-        input::ArrayOrNamedTuple,
-        output::ArrayOrNamedTuple)
-    norm(loss.hvf(input, ps) - output) / norm(output)
-end
diff --git a/src/loss/lnn_loss.jl b/src/loss/lnn_loss.jl
deleted file mode 100644
index b27b87eae..000000000
--- a/src/loss/lnn_loss.jl
+++ /dev/null
@@ -1,135 +0,0 @@
-@doc raw"""
-    lagrangian_acceleration(arch::LagrangianNeuralNetwork)
-
-Compute an executable expression of the acceleration that the Lagrangian `arch` predicts.
-
-The Euler–Lagrange equations of a Lagrangian ``L(q, \dot{q})`` are
-
-```math
-    \frac{d}{dt}\frac{\partial{}L}{\partial\dot{q}} = \frac{\partial{}L}{\partial{}q},
-```
-
-which, expanded along a trajectory, reads
-
-```math
-    \nabla_{\dot{q}}\nabla_{\dot{q}}L\,\ddot{q}
-        + (\nabla_q\nabla_{\dot{q}}L)^T\dot{q}
-        = \nabla_qL,
-```
-
-so the acceleration is
-
-```math
-    \ddot{q} = (\nabla_{\dot{q}}\nabla_{\dot{q}}L)^{-1}
-               \left( \nabla_qL - (\nabla_q\nabla_{\dot{q}}L)^T\dot{q} \right).
-```
-
-The transpose is not cosmetic. ``\nabla_q\nabla_{\dot{q}}L`` is indexed
-``[i, j] = \partial^2L/\partial{}q_i\partial\dot{q}_j``, while the chain rule contracts the *first*
-index with ``\dot{q}``. It disappears only for a Lagrangian whose position–velocity coupling happens
-to be symmetric.
-
-# Implementation
-
-The gradient and the Hessian of the network with respect to its input are built as *symbolic*
-expressions with `SymbolicNeuralNetworks.Jacobian` — applied once, and then again to its own result —
-and compiled with `build_nn_function`. They are not taken with `Zygote` inside the loss:
-differentiating a loss that itself contains a `Zygote.gradient` with respect to the network
-parameters fails with `MethodError: no method matching getindex(::IdDict{Any, Any})`. This is the
-route [`hamiltonian_vector_field`](@ref) takes, and for the same reason.
-
-The compiled functions flatten a batch along the columns — ``(1, 2n\cdot{}B)`` for the gradient and
-``(2n, 2n\cdot{}B)`` for the Hessian — so both are reshaped before use. Trailing input dimensions
-beyond the first are flattened on the way in and restored on the way out, so an input of size
-``(2n, s, B)``, which is what [`Batch`](@ref) hands the loss, yields an acceleration of size
-``(n, s, B)``.
-"""
-function lagrangian_acceleration(arch::LagrangianNeuralNetwork)
-    snn = SymbolicNeuralNetwork(arch)
-    ∇L = derivative(SymbolicNeuralNetworks.Jacobian(snn))
-    ∇∇L = derivative(SymbolicNeuralNetworks.Jacobian(∇L, snn))
-    grad = SymbolicNeuralNetworks.build_nn_function(
-        ∇L, snn.params, snn.input; inplace = false)
-    hess = SymbolicNeuralNetworks.build_nn_function(
-        ∇∇L, snn.params, snn.input; inplace = false)
-    n = dim(arch) ÷ 2
-
-    function (input, ps)
-        input_matrix = reshape(input, 2n, :)
-        batch_size = size(input_matrix, 2)
-        g = reshape(grad(input_matrix, ps), 2n, batch_size)
-        H = reshape(hess(input_matrix, ps), 2n, 2n, batch_size)
-        acceleration = reduce(hcat,
-            [_euler_lagrange_acceleration(H[:, :, b], g[:, b],
-                 input_matrix[(n + 1):(2n), b])
-             for b in 1:batch_size])
-        reshape(acceleration, n, size(input)[2:end]...)
-    end
-end
-
-@doc raw"""
-    _euler_lagrange_acceleration(H, g, q̇)
-
-Solve the Euler–Lagrange equations of one sample for the acceleration.
-
-`H` is the full ``2n\times{}2n`` Hessian of the Lagrangian with respect to its input, `g` the full
-``2n`` gradient, and `q̇` the ``n`` velocity components of that same input.
-"""
-function _euler_lagrange_acceleration(H, g, q̇)
-    n = length(q̇)
-    H[(n + 1):(2n), (n + 1):(2n)] \ (g[1:n] - H[1:n, (n + 1):(2n)]' * q̇)
-end
-
-@doc raw"""
-    LNNLoss <: NetworkLoss
-
-The loss for a Lagrangian neural network.
-
-The network output is a scalar Lagrangian ``L(q, \dot{q})``. The loss compares the acceleration its
-Euler–Lagrange equations predict against the acceleration in the data:
-
-```math
-    \mathtt{loss}(\mathcal{NN}, \mathtt{input}, \mathtt{output})
-        = ||\ddot{q}(\mathtt{input}) - \mathtt{output}|| \,/\, ||\mathtt{output}||,
-```
-
-where ``\mathtt{input}`` stacks ``q`` on ``\dot{q}`` and is therefore ``2n`` rows tall, and
-``\mathtt{output}`` is ``\ddot{q}`` and is ``n`` rows tall. See [`lagrangian_acceleration`](@ref)
-for ``\ddot{q}``.
-
-# Constructor
-
-This can be called with a `LagrangianNeuralNetwork` as its only argument:
-```julia
-LNNLoss(arch)
-```
-
-# Functor
-
-```julia
-loss(model, ps, input, output)
-loss(ps, input, output) # equivalent to the above
-```
-"""
-struct LNNLoss{FT} <: NetworkLoss
-    acceleration::FT
-end
-
-function LNNLoss(arch::LagrangianNeuralNetwork)
-    LNNLoss(lagrangian_acceleration(arch))
-end
-
-AbstractNeuralNetworks.NetworkLoss(arch::LagrangianNeuralNetwork) = LNNLoss(arch)
-
-function (loss::LNNLoss)(::Union{Chain, AbstractExplicitLayer},
-        ps::NetworkParameters,
-        input::AbstractArray,
-        output::AbstractArray)
-    loss(ps, input, output)
-end
-
-function (loss::LNNLoss)(ps::NetworkParameters,
-        input::AbstractArray,
-        output::AbstractArray)
-    norm(loss.acceleration(input, ps) - output) / norm(output)
-end
diff --git a/src/loss/symplectic_euler_loss.jl b/src/loss/symplectic_euler_loss.jl
deleted file mode 100644
index b3bdf6dbc..000000000
--- a/src/loss/symplectic_euler_loss.jl
+++ /dev/null
@@ -1,92 +0,0 @@
-@doc raw"""
-    SymplecticEulerLoss <: NetworkLoss
-
-The loss that trains a Hamiltonian neural network on a *trajectory*, through one step of a
-symplectic Euler method.
-
-[`HNNLoss`](@ref) needs ``(\dot{q}, \dot{p})`` in the data. This one does not: it needs two
-consecutive states and the timestep between them, and it asks that one symplectic Euler step of the
-learned Hamiltonian carry the first state to the second.
-
-Variant `:A` evaluates the Hamiltonian vector field at ``(q_{n+1}, p_n)`` and variant `:B` at
-``(q_n, p_{n+1})``. Writing ``X_H`` for that vector field and ``\Delta{}t`` for the timestep, the
-residual of variant `:A` is
-
-```math
-    X_H(q_{n+1}, p_n) - \frac{1}{\Delta{}t}
-        \begin{pmatrix} q_{n+1} - q_n \\ p_{n+1} - p_n \end{pmatrix},
-```
-
-and the loss is its norm relative to the norm of the finite difference, so it is scale invariant in
-the same way [`HNNLoss`](@ref) is.
-
-# Constructor
-
-```julia
-SymplecticEulerLoss(arch, timestep)                 # variant :A
-SymplecticEulerLoss(arch, timestep; variant = :B)
-```
-
-where `arch` is a [`HamiltonianArchitecture`](@ref).
-
-# Functor
-
-```julia
-loss(model, ps, input, output)
-loss(ps, input, output) # equivalent to the above
-```
-
-`input` is ``(q_n, p_n)`` stacked and `output` is ``(q_{n+1}, p_{n+1})`` stacked, so both have
-``2n`` rows.
-"""
-struct SymplecticEulerLoss{V, FT, T} <: NetworkLoss
-    hvf::FT
-    timestep::T
-
-    function SymplecticEulerLoss(arch::HamiltonianArchitecture, timestep::Real;
-            variant::Symbol = :A)
-        variant in (:A, :B) ||
-            throw(ArgumentError("variant must be :A or :B, got :$(variant)"))
-        hvf = hamiltonian_vector_field(arch)
-        new{variant, typeof(hvf), typeof(timestep)}(hvf, timestep)
-    end
-end
-
-@doc raw"""
-    _symplectic_euler_evaluation_point(loss, input, output)
-
-The staggered state at which the Hamiltonian vector field is evaluated: ``(q_{n+1}, p_n)`` for
-variant `:A`, ``(q_n, p_{n+1})`` for variant `:B`.
-
-Taking one half of the state from `input` and the other from `output` is what makes the step
-implicit, and is the whole difference from an explicit Euler residual.
-"""
-function _symplectic_euler_evaluation_point(
-        ::SymplecticEulerLoss{:A}, input::AbstractArray, output::AbstractArray)
-    n = size(input, 1) ÷ 2
-    vcat(selectdim(output, 1, 1:n), selectdim(input, 1, (n + 1):(2n)))
-end
-
-function _symplectic_euler_evaluation_point(
-        ::SymplecticEulerLoss{:B}, input::AbstractArray, output::AbstractArray)
-    n = size(input, 1) ÷ 2
-    vcat(selectdim(input, 1, 1:n), selectdim(output, 1, (n + 1):(2n)))
-end
-
-function (loss::SymplecticEulerLoss)(::Union{Chain, AbstractExplicitLayer},
-        ps::NetworkParameters, input::AbstractArray, output::AbstractArray)
-    loss(ps, input, output)
-end
-
-function (loss::SymplecticEulerLoss)(
-        ps::NetworkParameters, input::AbstractArray, output::AbstractArray)
-    evaluation_point = _symplectic_euler_evaluation_point(loss, input, output)
-    # The timestep is converted to the data's element type rather than promoted against it: a
-    # `Float64` timestep must not widen the loss of a `Float32` network. `HNNLoss` has no constant
-    # to convert, and `symbolic_hamiltonian_vector_field` builds its Poisson tensor from `Int` for
-    # this same reason.
-    difference_quotient = (output - input) / eltype(input)(loss.timestep)
-    field = reshape(
-        loss.hvf(reshape(evaluation_point, size(input, 1), :), ps), size(input))
-    norm(field - difference_quotient) / norm(difference_quotient)
-end
diff --git a/src/loss/variational_midpoint_loss.jl b/src/loss/variational_midpoint_loss.jl
deleted file mode 100644
index 3e9b5da8c..000000000
--- a/src/loss/variational_midpoint_loss.jl
+++ /dev/null
@@ -1,107 +0,0 @@
-@doc raw"""
-    VariationalMidpointLoss <: NetworkLoss
-
-The loss that trains a Lagrangian neural network on positions alone, through the discrete
-Euler–Lagrange equations of the midpoint discrete Lagrangian.
-
-It needs three consecutive positions and the timestep between them — no velocities and no
-accelerations.
-
-The midpoint discrete Lagrangian is
-
-```math
-    L_d(q_n, q_{n+1}) = \Delta{}t\,
-        L\!\left( \frac{q_n + q_{n+1}}{2}, \frac{q_{n+1} - q_n}{\Delta{}t} \right),
-```
-
-and the discrete Euler–Lagrange equations ask that
-
-```math
-    D_2L_d(q_n, q_{n+1}) + D_1L_d(q_{n+1}, q_{n+2}) = 0 .
-```
-
-The loss is the norm of that residual, relative to the norm of the two derivatives it is built from,
-so a Lagrangian scaled by a constant gives the same number.
-
-# Implementation
-
-The two slot derivatives are written out with the chain rule rather than taken with a nested
-`Zygote.gradient`:
-
-```math
-    D_1L_d = \frac{\Delta{}t}{2}\nabla_qL - \nabla_{\dot{q}}L, \qquad
-    D_2L_d = \frac{\Delta{}t}{2}\nabla_qL + \nabla_{\dot{q}}L,
-```
-
-both evaluated at the midpoint. Only the *first* derivative of the network with respect to its input
-is then needed, and that comes from a compiled symbolic expression exactly as
-[`hamiltonian_vector_field`](@ref)'s does. Written the other way — a `Zygote.gradient` inside the
-loss — the parameter gradient fails with
-`MethodError: no method matching getindex(::IdDict{Any, Any})`.
-
-# Constructor
-
-```julia
-VariationalMidpointLoss(arch, timestep)
-```
-
-where `arch` is a `LagrangianNeuralNetwork`.
-
-# Functor
-
-```julia
-loss(model, ps, input, output)
-loss(ps, input, output) # equivalent to the above
-```
-
-`input` stacks ``q_n`` on ``q_{n+1}`` and so has ``2n`` rows; `output` is ``q_{n+2}`` and has ``n``.
-"""
-struct VariationalMidpointLoss{FT, T} <: NetworkLoss
-    gradient::FT
-    timestep::T
-
-    function VariationalMidpointLoss(arch::LagrangianNeuralNetwork, timestep::Real)
-        snn = SymbolicNeuralNetwork(arch)
-        ∇L = derivative(SymbolicNeuralNetworks.Jacobian(snn))
-        grad = SymbolicNeuralNetworks.build_nn_function(
-            ∇L, snn.params, snn.input; inplace = false)
-        new{typeof(grad), typeof(timestep)}(grad, timestep)
-    end
-end
-
-@doc raw"""
-    _discrete_lagrangian_derivatives(loss, ps, qa, qb)
-
-Return `(D₁Lᵈ, D₂Lᵈ)` of the midpoint discrete Lagrangian on the pair `(qa, qb)`, one column per
-sample.
-"""
-function _discrete_lagrangian_derivatives(loss::VariationalMidpointLoss, ps,
-        qa::AbstractMatrix, qb::AbstractMatrix)
-    n, batch_size = size(qa)
-    # The timestep is converted to the data's element type rather than promoted against it: a
-    # `Float64` timestep must not widen the loss of a `Float32` network.
-    Δt = eltype(qa)(loss.timestep)
-    midpoint = vcat((qa + qb) / 2, (qb - qa) / Δt)
-    g = reshape(loss.gradient(midpoint, ps), 2n, batch_size)
-    ∇qL = g[1:n, :]
-    ∇q̇L = g[(n + 1):(2n), :]
-    (Δt / 2) * ∇qL - ∇q̇L, (Δt / 2) * ∇qL + ∇q̇L
-end
-
-function (loss::VariationalMidpointLoss)(::Union{Chain, AbstractExplicitLayer},
-        ps::NetworkParameters, input::AbstractArray, output::AbstractArray)
-    loss(ps, input, output)
-end
-
-function (loss::VariationalMidpointLoss)(
-        ps::NetworkParameters, input::AbstractArray, output::AbstractArray)
-    n = size(input, 1) ÷ 2
-    qₙ = reshape(selectdim(input, 1, 1:n), n, :)
-    qₙ₊₁ = reshape(selectdim(input, 1, (n + 1):(2n)), n, :)
-    qₙ₊₂ = reshape(output, n, :)
-
-    _, D₂ = _discrete_lagrangian_derivatives(loss, ps, qₙ, qₙ₊₁)
-    D₁, _ = _discrete_lagrangian_derivatives(loss, ps, qₙ₊₁, qₙ₊₂)
-
-    norm(D₂ + D₁) / norm(vcat(D₂, D₁))
-end
diff --git a/src/pullbacks/symbolic_hnn_pullback.jl b/src/pullbacks/symbolic_hnn_pullback.jl
deleted file mode 100644
index 01cd9a023..000000000
--- a/src/pullbacks/symbolic_hnn_pullback.jl
+++ /dev/null
@@ -1,50 +0,0 @@
-"""
-    SymbolicPullback(arch::HamiltonianArchitecture)
-
-Make a `SymbolicPullback` of [`HNNLoss`](@ref) for a [`HamiltonianArchitecture`](@ref).
-
-Its value is `HNNLoss` on the batch, ``\\ell = \\|R\\|/\\|Y\\|`` with ``R = X_H(X) - Y``. Its
-gradient is ``\\Delta \\sum_i \\nabla_\\theta \\tfrac12\\|r_i\\|^2 / (\\|R\\|\\,\\|Y\\|)``: the generated
-per-sample gradient of the surrogate ``\\tfrac12\\|r_i\\|^2``, summed over the batch and scaled by the
-two norms of the value's residual. The gradient has the parameters' element type, and a zero
-residual gives a zero gradient.
-
-`SymbolicNeuralNetworks.SymbolicPullback(nn, loss)` cannot be used: it sizes the target of the loss
-with `output_dimension(nn.model)`, the scalar Hamiltonian, but `HNNLoss` compares against the vector
-field, which has the dimension of the input.
-"""
-function SymbolicPullback(arch::HamiltonianArchitecture)
-    nn = SymbolicNeuralNetwork(arch)
-    hvf = symbolic_hamiltonian_vector_field(nn)
-    # `inplace = false`, as in `hamiltonian_vector_field`: `HNNLoss` is also differentiated by `Zygote`.
-    loss = HNNLoss(SymbolicNeuralNetworks.build_nn_function(
-        hvf, nn.params, nn.input; inplace = false))
-    soutput = Symbolics.variables(:y, 1:input_dimension(nn.model))
-    gradient = SymbolicNeuralNetworks.symbolic_parameter_gradient(
-        sum(abs2, hvf - soutput) / 2, nn)
-    gradient_function = SymbolicNeuralNetworks.build_nn_function(
-        gradient, nn.params, nn.input, soutput; reduce = +)
-    SymbolicPullback(loss, _HNNParameterGradient(gradient_function))
-end
-
-# The `fun` of the HNN pullback; it holds the generated gradient and selects the call method below.
-struct _HNNParameterGradient{FT}
-    gradient_function::FT
-end
-
-# More specific than `SymbolicNeuralNetworks`' call, which computes the value and the gradient
-# separately and would evaluate `hvf` twice: narrower in the pullback and the parameters, and the
-# same in the data.
-function (pb::SymbolicPullback{<:HNNLoss, <:_HNNParameterGradient})(
-        ps::NetworkParameters, model,
-        (input, output)::Tuple{<:ArrayOrNamedTuple, <:ArrayOrNamedTuple})
-    norm_R = norm(pb.loss.hvf(input, ps) - output)
-    norm_Y = norm(output)
-    # `pinv`: a zero residual gives the zero gradient, not `NaN`. Each norm divides each entry on its
-    # own, because `‖R‖‖Y‖` overflows, and its inverse underflows, where neither norm does.
-    inv_R = pinv(norm_R)
-    TP = parameter_eltype(ps)
-    gradient_function = pb.fun.gradient_function
-    norm_R / norm_Y,
-    Δ -> mapstorage(x -> TP.(x .* Δ .* inv_R ./ norm_Y), gradient_function(input, output, ps))
-end
diff --git a/test/architectures/forward_inference.jl b/test/architectures/forward_inference.jl
index 4b94736a4..54ba6a113 100644
--- a/test/architectures/forward_inference.jl
+++ b/test/architectures/forward_inference.jl
@@ -1,6 +1,5 @@
 # The forward pass of every architecture that the package exports infers its result type, and the
-# result has the element type of the network. `GeneralizedHamiltonianArchitecture` has no chain:
-# its constructor raises "GHNN still has to be implemented!".
+# result has the element type of the network.
 
 using GeometricMachineLearning
 using Test
@@ -18,7 +17,6 @@ const ARCHITECTURES = (
     ("StandardTransformerIntegrator, default heads",
         StandardTransformerIntegrator(32), 32, 3),
     ("VolumePreservingTransformer", VolumePreservingTransformer(4, 3), 4, 3),
-    ("LagrangianNeuralNetwork", LagrangianNeuralNetwork(4), 4, 3),
     ("LASympNet", LASympNet(4), 4, 3),
     ("GSympNet", GSympNet(4), 4, 3),
     ("ClassificationTransformer",
@@ -26,7 +24,6 @@ const ARCHITECTURES = (
     ("VolumePreservingFeedForward", VolumePreservingFeedForward(4), 4, 3),
     ("SymplecticAutoencoder", SymplecticAutoencoder(8, 4), 8, 3),
     ("PSDArch", PSDArch(8, 4), 8, 3),
-    ("StandardHamiltonianArchitecture", StandardHamiltonianArchitecture(4), 4, 3),
     ("LinearSymplecticTransformer", LinearSymplecticTransformer(4, 3), 4, 3),
     ("SymplecticTransformer", SymplecticTransformer(4), 4, 3))
 
diff --git a/test/architectures/hamiltonian_neural_network_tests.jl b/test/architectures/hamiltonian_neural_network_tests.jl
deleted file mode 100644
index d4ecbdf64..000000000
--- a/test/architectures/hamiltonian_neural_network_tests.jl
+++ /dev/null
@@ -1,63 +0,0 @@
-using GeometricMachineLearning
-using Test
-using Zygote
-import Random
-Random.seed!(1234)
-
-const GMLA = GeometricMachineLearning.Activation
-function allocate_network_and_data_loader(dim::Integer, width::Integer, nhidden::Integer, activation::GMLA)
-    arch = StandardHamiltonianArchitecture(dim, width, nhidden, activation)
-    loss = HNNLoss(arch)
-    input = rand(dim, 10)
-    output = rand(dim, 10)
-    dl = DataLoader(input, output)
-    NeuralNetwork(arch), loss, dl
-end
-
-"""
-This tests if we can call the HNN loss.
-"""
-function test_hnn_loss(dim::Integer = 2,
-        width::Integer = dim,
-        nhidden::Integer = 1,
-        activation::GMLA = GeometricMachineLearning.SigmoidActivation())
-    nn, loss, dl = allocate_network_and_data_loader(dim, width, nhidden, activation)
-    @test typeof(loss(nn.params, dl.input, dl.output)) <: Real
-    @test typeof(loss) <: NetworkLoss
-end
-
-test_hnn_loss()
-
-"""
-This tests if we can differentiate the HNN loss.
-"""
-function test_hnn_loss_derivative(dim::Integer = 2,
-        width::Integer = dim,
-        nhidden::Integer = 1,
-        activation::GMLA = GeometricMachineLearning.SigmoidActivation())
-    nn, loss, dl = allocate_network_and_data_loader(dim, width, nhidden, activation)
-    dp = Zygote.gradient(ps -> loss(ps, dl.input, dl.output), nn.params)[1]
-    @test typeof(dp) <: NetworkParameters
-    @test keys(dp) == keys(nn.params)
-end
-
-test_hnn_loss_derivative()
-
-"""
-This tests if we can build the symbolic pullback of the HNN loss and evaluate it.
-"""
-function test_hnn_symbolic_pullback(dim::Integer = 2,
-        width::Integer = dim,
-        nhidden::Integer = 1,
-        activation::GMLA = GeometricMachineLearning.SigmoidActivation())
-    arch = StandardHamiltonianArchitecture(dim, width, nhidden, activation)
-    nn = NeuralNetwork(arch)
-    pb = SymbolicPullback(arch)
-    input = rand(dim, 10)
-    output = rand(dim, 10)
-    loss_value, back = pb(nn.params, nn.model, (input, output))
-    @test typeof(loss_value) <: Real
-    @test keys(back(1)) == keys(nn.params)
-end
-
-test_hnn_symbolic_pullback()
diff --git a/test/architectures/lagrangian_neural_network_tests.jl b/test/architectures/lagrangian_neural_network_tests.jl
deleted file mode 100644
index 660578e35..000000000
--- a/test/architectures/lagrangian_neural_network_tests.jl
+++ /dev/null
@@ -1,163 +0,0 @@
-using GeometricMachineLearning
-using GeometricMachineLearning: _euler_lagrange_acceleration, lagrangian_acceleration,
-                                params
-using LinearAlgebra
-using Test
-using Zygote
-import Random
-
-Random.seed!(1234)
-
-@doc raw"""
-The Euler–Lagrange solve must *be* the Euler–Lagrange equation, not merely run.
-
-For ``L = \tfrac12\dot{q}^TM\dot{q} + q^TC\dot{q} - \tfrac12q^TKq`` the equations
-``\tfrac{d}{dt}\partial{}L/\partial\dot{q} = \partial{}L/\partial{}q`` read
-``M\ddot{q} + C^T\dot{q} = C\dot{q} - Kq``, so
-
-```math
-    \ddot{q} = M^{-1}\left((C - C^T)\dot{q} - Kq\right),
-```
-
-in closed form and independently of any network. `C` is deliberately *not* symmetric: with a
-symmetric `C` the transpose in the implementation cancels and the test cannot see it.
-"""
-function test_euler_lagrange_acceleration_against_closed_form(n::Integer)
-    M = (A = rand(n, n); A' * A + n * I)
-    K = (B = rand(n, n); B' * B + n * I)
-    C = rand(n, n)
-    q = rand(n)
-    q̇ = rand(n)
-
-    H = [-K C; C' M]
-    g = vcat(C * q̇ - K * q, M * q̇ + C' * q)
-    reference = M \ ((C - C') * q̇ - K * q)
-
-    @test _euler_lagrange_acceleration(H, g, q̇) ≈ reference
-
-    # The same expression with the transpose dropped is a different one, so the assertion above
-    # is not passing for a reason unrelated to the index convention.
-    without_transpose = M \ (g[1:n] - C * q̇)
-    @test !isapprox(without_transpose, reference)
-end
-
-for n in (2, 3, 4)
-    test_euler_lagrange_acceleration_against_closed_form(n)
-end
-
-"""
-The compiled symbolic derivatives must agree with `Zygote`'s, per sample and in a batch.
-"""
-function test_lagrangian_acceleration_against_zygote(n::Integer, batch_size::Integer)
-    arch = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    acceleration = lagrangian_acceleration(arch)
-
-    input = randn(2n, batch_size)
-    predicted = acceleration(input, params(nn))
-    @test size(predicted) == (n, batch_size)
-
-    for b in 1:batch_size
-        x = input[:, b]
-        H = Zygote.hessian(y -> sum(nn(y, params(nn))), x)
-        g = Zygote.gradient(y -> sum(nn(y, params(nn))), x)[1]
-        @test predicted[:, b] ≈ _euler_lagrange_acceleration(H, g, x[(n + 1):(2n)])
-    end
-end
-
-test_lagrangian_acceleration_against_zygote(1, 4)
-test_lagrangian_acceleration_against_zygote(2, 3)
-
-"""
-Trailing input dimensions must be carried through. `Batch` hands the loss a `(2n, s, B)` array, not
-a matrix, so the documented reshape is on the path every training run takes.
-"""
-function test_lagrangian_acceleration_keeps_trailing_dimensions(n::Integer = 2)
-    arch = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    acceleration = lagrangian_acceleration(arch)
-
-    input = randn(2n, 3, 5)
-    predicted = acceleration(input, params(nn))
-    @test size(predicted) == (n, 3, 5)
-    @test reshape(predicted, n, :) ≈ acceleration(reshape(input, 2n, :), params(nn))
-end
-
-test_lagrangian_acceleration_keeps_trailing_dimensions()
-
-@doc raw"""
-The solve inverts ``\nabla_{\dot{q}}\nabla_{\dot{q}}L``, so the acceleration is only as trustworthy
-as that block's conditioning. This is what decided the solved form over an un-inverted residual, so
-it is asserted rather than left in a commit message.
-"""
-function test_velocity_hessian_is_well_conditioned(n::Integer; draws = 10, points = 10)
-    # The architecture matches the tests above, so `Zygote.hessian` is already specialised for it.
-    arch = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    worst = 0.0
-    for _ in 1:draws
-        nn = NeuralNetwork(arch)
-        for _ in 1:points
-            H = Zygote.hessian(y -> sum(nn(y, params(nn))), randn(2n))
-            worst = max(worst, cond(H[(n + 1):(2n), (n + 1):(2n)]))
-        end
-    end
-    @test worst < 1.0e8
-end
-
-for n in (1, 2)
-    test_velocity_hessian_is_well_conditioned(n)
-end
-
-"""
-`LNNLoss` must be the default loss of the architecture, be callable, and be differentiable with
-respect to the parameters.
-"""
-function test_lnn_loss(n::Integer = 1)
-    arch = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    loss = LNNLoss(arch)
-
-    @test typeof(loss) <: NetworkLoss
-    @test typeof(NetworkLoss(arch)) <: LNNLoss
-
-    input = randn(2n, 8)
-    output = randn(n, 8)
-    @test typeof(loss(params(nn), input, output)) <: Real
-    @test loss(nn.model, params(nn), input, output) == loss(params(nn), input, output)
-
-    dp = Zygote.gradient(ps -> loss(ps, input, output), params(nn))[1]
-    @test typeof(dp) <: NetworkParameters
-    @test keys(dp) == keys(params(nn))
-    @test !all(iszero, dp.L1.W)
-end
-
-test_lnn_loss(1)
-test_lnn_loss(2)
-
-"""
-The loss must reach zero on data a Lagrangian can actually reproduce, so train against the
-accelerations the *network's own* Euler–Lagrange equations give at its initial parameters. Those
-are reachable by construction, which makes "the loss decreases" a statement about training rather
-than about how hard the target was.
-"""
-function test_lnn_training(; n = 1, n_epochs = 30)
-    # The width and depth deliberately match the tests above. A different `Chain` length is a fresh
-    # `Zygote` specialisation of the whole loss, which costs 61 s of this file's wall time and buys
-    # nothing: the assertion is that training reduces the loss, and that does not depend on either.
-    arch = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    loss = LNNLoss(arch)
-
-    input = randn(2n, 60)
-    output = loss.acceleration(input, params(nn))
-
-    # Start the training from a different draw, so the target is not already met.
-    nn = NeuralNetwork(arch)
-    dl = DataLoader(input, output; suppress_info = true)
-    optimizer = TrainingOptimizer(nn.params; algorithm = Adam())
-    loss_array = train!(nn, optimizer, dl, Batch(10), n_epochs; loss, show_progress = false)
-
-    @test loss_array[end] < loss_array[1]
-end
-
-test_lnn_training()
diff --git a/test/data_loader/train.jl b/test/data_loader/train.jl
index dfa59311b..0486d25ae 100644
--- a/test/data_loader/train.jl
+++ b/test/data_loader/train.jl
@@ -6,7 +6,6 @@ using Test
 import Random
 
 include("../helpers/eltypes.jl")
-include("../helpers/leaves.jl")
 
 # The target is in the model class, so the loss goes to zero; the factor 100 leaves margin in
 # `Float32`.
@@ -111,36 +110,6 @@ end
     @test history[1] ≈ sum(losses) / 3
 end
 
-# `train!` steps with the gradient of the `loss` it is given. One epoch of one minibatch with
-# `GradientMethod` moves the parameters by `-η` times the symbolic gradient at `p₀`, and the history
-# is the pullback's loss value. `test/pullbacks/symbolic_hnn_pullback.jl` tests the gradient itself
-# against `ZygotePullback`'s (issue #343).
-@testset "train! trains a Hamiltonian network with SymbolicPullback ($T)" for T in REAL_ELTYPES
-    Random.seed!(9)
-    arch = StandardHamiltonianArchitecture(2, 4, 1, GeometricMachineLearning.SigmoidActivation())
-    dl = DataLoader(randn(T, 2, 8), randn(T, 2, 8); suppress_info = true)
-    nn = NeuralNetwork(arch, CPU(), T)
-    p₀ = deepcopy(params(nn))
-    pullback = SymbolicPullback(arch)
-    batch = Batch(8)
-    Random.seed!(2)
-    data = GeometricMachineLearning.convert_input_and_batch_indices_to_array(
-        dl, batch, only(batch(dl)))
-    value, back = pullback(p₀, nn.model, data)
-    g = back(one(value))
-    η = T(0.1)
-    opt = TrainingOptimizer(params(nn); algorithm = GradientMethod(), linesearch = η)
-    Random.seed!(2)
-    history = train!(nn, opt, dl, batch, 1; loss = pullback, show_progress = false)
-    @test history == [value]
-    @test value isa T
-    @test eltype(history) === T
-    foreach_leaf(params(nn), p₀, g) do layer, layer₀, gradient
-        @test layer ≈ layer₀ .- η .* gradient
-        @test layer != layer₀
-    end
-end
-
 # Without `loss`, `train!` trains on the loss that matches the architecture and the batch, and on
 # any other combination it says that it needs one.
 #
@@ -156,11 +125,9 @@ end
     @test NetworkLoss(GSympNet(4; n_layers = 2), Batch(4)) isa FeedForwardLoss
     autoencoder = SymplecticAutoencoder(10, 4)
     @test NetworkLoss(autoencoder, Batch(4)) isa AutoEncoderLoss
-    # a Hamiltonian network, a transformer with a `Batch{:FeedForward}` and an unknown architecture
-    # all reach the fallback
-    for (arch, batch) in ((StandardHamiltonianArchitecture(2, 4, 1), Batch(4)),
-        (transformer, Batch(4)), (
-        NeuralNetwork(Chain(Dense(2, 2, identity))).architecture, Batch(4)))
+    # a transformer with a `Batch{:FeedForward}` and an unknown architecture both reach the fallback
+    for (arch, batch) in ((transformer, Batch(4)),
+        (NeuralNetwork(Chain(Dense(2, 2, identity))).architecture, Batch(4)))
         @test_throws ArgumentError NetworkLoss(arch, batch)
     end
 
diff --git a/test/integration/train_manual_architectures.jl b/test/integration/train_manual_architectures.jl
index 63649db55..98206a051 100644
--- a/test/integration/train_manual_architectures.jl
+++ b/test/integration/train_manual_architectures.jl
@@ -13,7 +13,6 @@ include("../helpers/eltypes.jl")
 function manual_architectures(T)
     feedforward = DataLoader(rand(T, 4, 12, 3); suppress_info = true)
     transformer_batch = Batch(4, 3)
-    hnn = StandardHamiltonianArchitecture(2, 4)
     attention_data = DataLoader(rand(T, 1, 12, 2); suppress_info = true)
     attention_batch = Batch(4, 3, 1)
     attention_chain(skew_sym) = Chain(Dense(1, 2, tanh),
@@ -38,9 +37,6 @@ function manual_architectures(T)
             SymplecticAutoencoder(4, 2; n_encoder_blocks = 2, n_decoder_blocks = 2),
             DataLoader(rand(T, 4, 12); autoencoder = true, suppress_info = true), Batch(4),
             nothing),
-        ("StandardHamiltonianArchitecture", hnn,
-            DataLoader(rand(T, 2, 12), rand(T, 2, 12); suppress_info = true), Batch(4),
-            HNNLoss(hnn)),
         ("VolumePreservingAttention, skew_sym = true", attention_chain(true),
             attention_data, attention_batch, TransformerLoss(attention_batch)),
         ("VolumePreservingAttention, skew_sym = false", attention_chain(false),
diff --git a/test/loss/losses.jl b/test/loss/losses.jl
index 2d75aaad6..bd6853995 100644
--- a/test/loss/losses.jl
+++ b/test/loss/losses.jl
@@ -1,6 +1,5 @@
 using GeometricMachineLearning
-using GeometricMachineLearning: AutoEncoderLoss, FeedForwardLoss, HNNLoss, LNNLoss,
-                                SymplecticEulerLoss, VariationalMidpointLoss
+using GeometricMachineLearning: AutoEncoderLoss, FeedForwardLoss
 using AbstractNeuralNetworks: Chain, Dense, NeuralNetwork, params
 using NeuralNetworkParameters: NetworkParameters
 using Test
@@ -29,35 +28,20 @@ end
 
 # A loss's `ps` is a `NetworkParameters`, so a bare `NamedTuple` reaches no method of a
 # GML loss rather than a silently different gradient: which of the two AD hands the loss is what an
-# untyped argument hides. The four-argument methods of a training loss carry the same annotation,
-# and their model slot is reached by `AbstractNeuralNetworks`' own untyped fallback where the wrapper
-# is missing.
+# untyped argument hides. The four-argument methods of a loss carry the annotation, and their model
+# slot is reached by `AbstractNeuralNetworks`' own untyped fallback where the wrapper is missing.
 #
 # Not numeric: the element type of the data does not reach the method selection, so this is not run
 # over `REAL_ELTYPES`.
 @testset "a loss takes a NetworkParameters" begin
     Random.seed!(23)
-    hnn_arch = StandardHamiltonianArchitecture(2, 4, 1)
-    lnn_arch = LagrangianNeuralNetwork(2; width = 4, nhidden = 1)
 
-    hnn_input, hnn_output = randn(2, 3), randn(2, 3)
-    lnn_input, lnn_output = randn(2, 3), randn(1, 3)
-
-    losses = (
-        (HNNLoss(hnn_arch), hnn_arch, hnn_input, hnn_output),
-        (LNNLoss(lnn_arch), lnn_arch, lnn_input, lnn_output),
-        (SymplecticEulerLoss(hnn_arch, 0.1), hnn_arch, hnn_input, hnn_output),
-        (VariationalMidpointLoss(lnn_arch, 0.05), lnn_arch, lnn_input, lnn_output)
-    )
-    for (loss, arch, input, output) in losses
-        nn = NeuralNetwork(arch, CPU(), Float64)
-        ps = params(nn)
-        # the wrapper is accepted, and the same number comes back through the `NeuralNetwork`
-        # method, which forwards to the wrapper
-        @test loss(nn, input, output) == loss(ps, input, output)
-        # the bare `NamedTuple` under the wrapper is not
-        @test_throws MethodError loss(params(ps), input, output)
-    end
+    # the three-argument method: the wrapper is accepted, and the bare `NamedTuple` under it is not
+    model = Chain(Dense(4, 2, tanh), Dense(2, 4, tanh))
+    ae = NeuralNetwork(model, CPU(), Float64)
+    x = randn(4, 3)
+    @test AutoEncoderLoss()(ae, x) == AutoEncoderLoss()(model, ae.params, x)
+    @test_throws MethodError AutoEncoderLoss()(model, params(ae.params), x)
 
     # `TransformerLoss` has no three-argument method, so its four-argument ones are what the
     # annotation reaches: the wrapper is accepted, and a bare `NamedTuple` under it reaches no method
diff --git a/test/loss/training_method_losses.jl b/test/loss/training_method_losses.jl
deleted file mode 100644
index d94134b1a..000000000
--- a/test/loss/training_method_losses.jl
+++ /dev/null
@@ -1,349 +0,0 @@
-using GeometricMachineLearning
-using GeometricMachineLearning: _symplectic_euler_evaluation_point,
-                                _discrete_lagrangian_derivatives,
-                                hamiltonian_vector_field, params
-using LinearAlgebra
-using Test
-using Zygote
-import Random
-
-Random.seed!(1234)
-
-# ---------------------------------------------------------------- ambiguity with `NetworkLoss`
-
-@doc raw"""
-Each three-argument loss functor `(loss)(ps, input, output)` must not collide with
-`AbstractNeuralNetworks`' own `(::NetworkLoss)(::NeuralNetwork, input, output)`. With an untyped
-first argument the two are ambiguous -- same argument count, neither more specific -- and
-`loss(nn, input, output)` is a `MethodError: ... is ambiguous` rather than a call. Typing the first
-parameter `::NetworkParameters` excludes a `NeuralNetwork`, so the call resolves to upstream's
-method, which itself forwards to `loss(nn.params, input, output)` -- the two must therefore return
-the same number.
-"""
-function test_loss_of_a_network_matches_loss_of_its_parameters()
-    hnn_arch = StandardHamiltonianArchitecture(2, 4, 1)
-    nn = NeuralNetwork(hnn_arch)
-    loss = HNNLoss(hnn_arch)
-    input, output = randn(2, 5), randn(2, 5)
-    @test loss(nn, input, output) == loss(params(nn), input, output)
-
-    lnn_arch = LagrangianNeuralNetwork(2; width = 4, nhidden = 1)
-    nn = NeuralNetwork(lnn_arch)
-    loss = LNNLoss(lnn_arch)
-    input, output = randn(2, 5), randn(1, 5)
-    @test loss(nn, input, output) == loss(params(nn), input, output)
-
-    se_arch = StandardHamiltonianArchitecture(2, 4, 1)
-    nn = NeuralNetwork(se_arch)
-    loss = SymplecticEulerLoss(se_arch, 0.1)
-    input, output = randn(2, 5), randn(2, 5)
-    @test loss(nn, input, output) == loss(params(nn), input, output)
-
-    vm_arch = LagrangianNeuralNetwork(2; width = 4, nhidden = 1)
-    nn = NeuralNetwork(vm_arch)
-    loss = VariationalMidpointLoss(vm_arch, 0.05)
-    input, output = randn(2, 5), randn(1, 5)
-    @test loss(nn, input, output) == loss(params(nn), input, output)
-end
-
-test_loss_of_a_network_matches_loss_of_its_parameters()
-
-# ---------------------------------------------------------------- SymplecticEulerLoss
-
-@doc raw"""
-The loss must vanish on a trajectory the learned Hamiltonian actually generates.
-
-Take the network's *own* Hamiltonian vector field, step it with symplectic Euler A, and feed the
-resulting pair back in: the residual is then zero by construction, so a non-zero value would mean
-the loss is not the method it claims to be.
-
-Symplectic Euler A is implicit in ``q``, so the step is taken by fixed-point iteration.
-"""
-function symplectic_euler_a_step(hvf, ps, state, Δt; iterations = 60)
-    n = length(state) ÷ 2
-    q, p = state[1:n], state[(n + 1):(2n)]
-    q_next = copy(q)
-    for _ in 1:iterations
-        field = hvf(vcat(q_next, p), ps)
-        q_next = q + Δt * field[1:n]
-    end
-    field = hvf(vcat(q_next, p), ps)
-    vcat(q_next, p + Δt * field[(n + 1):(2n)])
-end
-
-function test_symplectic_euler_loss_vanishes_on_its_own_step(n::Integer = 1; Δt = 1.0e-3)
-    arch = StandardHamiltonianArchitecture(2n, 4, 1)
-    nn = NeuralNetwork(arch)
-    loss = SymplecticEulerLoss(arch, Δt)
-    hvf = hamiltonian_vector_field(arch)
-
-    input = reduce(hcat, [randn(2n) for _ in 1:6])
-    output = reduce(hcat,
-        [symplectic_euler_a_step(hvf, params(nn), input[:, b], Δt) for b in axes(input, 2)])
-
-    @test loss(params(nn), input, output) < 1.0e-8
-
-    # A trajectory that is *not* the method's own step must not pass, so the assertion above is
-    # not vacuous.
-    @test loss(params(nn), input, output + 0.1 * randn(size(output))) > 1.0e-3
-end
-
-test_symplectic_euler_loss_vanishes_on_its_own_step(1)
-test_symplectic_euler_loss_vanishes_on_its_own_step(2)
-
-"""
-The two variants must differ, and each must take its half of the state from the right array.
-"""
-function test_symplectic_euler_variants(n::Integer = 2)
-    arch = StandardHamiltonianArchitecture(2n, 4, 1)
-    nn = NeuralNetwork(arch)
-    input = randn(2n, 5)
-    output = randn(2n, 5)
-
-    loss_a = SymplecticEulerLoss(arch, 0.1)
-    loss_b = SymplecticEulerLoss(arch, 0.1; variant = :B)
-
-    point_a = _symplectic_euler_evaluation_point(loss_a, input, output)
-    point_b = _symplectic_euler_evaluation_point(loss_b, input, output)
-    @test point_a == vcat(output[1:n, :], input[(n + 1):(2n), :])
-    @test point_b == vcat(input[1:n, :], output[(n + 1):(2n), :])
-    @test loss_a(params(nn), input, output) != loss_b(params(nn), input, output)
-
-    @test_throws ArgumentError SymplecticEulerLoss(arch, 0.1; variant = :C)
-end
-
-test_symplectic_euler_variants()
-
-"""
-Both variants must be differentiable with respect to the parameters. Training is the whole point of
-a loss, and a residual built on a symbolic vector field is exactly the kind of expression an AD
-system can refuse.
-"""
-function test_symplectic_euler_loss_gradient(n::Integer = 1)
-    arch = StandardHamiltonianArchitecture(2n, 4, 1)
-    nn = NeuralNetwork(arch)
-    input = randn(2n, 5)
-    output = randn(2n, 5)
-    for variant in (:A, :B)
-        loss = SymplecticEulerLoss(arch, 0.1; variant = variant)
-        dp = Zygote.gradient(ps -> loss(ps, input, output), params(nn))[1]
-        @test typeof(dp) <: NetworkParameters
-        @test keys(dp) == keys(params(nn))
-        @test !all(iszero, dp.L1.W)
-    end
-end
-
-test_symplectic_euler_loss_gradient()
-
-"""
-A `Float64` timestep must not widen the loss of a `Float32` network. The timestep is a constant of
-the method, not data, so the data's element type is what decides the result's.
-"""
-function test_timestep_does_not_widen_the_element_type(n::Integer = 1)
-    hamiltonian = StandardHamiltonianArchitecture(2n, 4, 1)
-    ps = params(NeuralNetwork(hamiltonian, Float32))
-    input = randn(Float32, 2n, 5)
-    output = randn(Float32, 2n, 5)
-    @test SymplecticEulerLoss(hamiltonian, 0.1)(ps, input, output) isa Float32
-    @test SymplecticEulerLoss(hamiltonian, 0.1f0)(ps, input, output) isa Float32
-
-    lagrangian = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    psl = params(NeuralNetwork(lagrangian, Float32))
-    @test VariationalMidpointLoss(lagrangian, 0.05)(
-        psl, randn(Float32, 2n, 5), randn(Float32, n, 5)) isa Float32
-end
-
-test_timestep_does_not_widen_the_element_type()
-
-function test_symplectic_euler_training(; n = 1, Δt = 1.0e-2, n_epochs = 20)
-    # The shape of the tests above. The assertion is that training reduces the loss, which does
-    # not depend on the width and depth, and a new shape is a fresh `Zygote` specialisation of the
-    # whole loss.
-    arch = StandardHamiltonianArchitecture(2n, 4, 1)
-    nn = NeuralNetwork(arch)
-    loss = SymplecticEulerLoss(arch, Δt)
-    hvf = hamiltonian_vector_field(arch)
-
-    input = reduce(hcat, [randn(2n) for _ in 1:60])
-    output = reduce(hcat,
-        [symplectic_euler_a_step(hvf, params(nn), input[:, b], Δt) for b in axes(input, 2)])
-
-    nn = NeuralNetwork(arch)
-    dl = DataLoader(input, output; suppress_info = true)
-    loss_array = train!(
-        nn, TrainingOptimizer(nn.params; algorithm = Adam()), dl, Batch(10), n_epochs;
-        loss, show_progress = false)
-    @test loss_array[end] < loss_array[1]
-end
-
-test_symplectic_euler_training()
-
-# ---------------------------------------------------------------- VariationalMidpointLoss
-
-@doc raw"""
-The slot derivatives must be the derivatives of the midpoint discrete Lagrangian.
-
-`_discrete_lagrangian_derivatives` writes the chain rule out by hand, so it has to be checked
-against a `Zygote` gradient of ``L_d`` itself — which is the thing it exists in order not to call.
-"""
-function test_discrete_lagrangian_derivatives(n::Integer; Δt = 0.07)
-    arch = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    loss = VariationalMidpointLoss(arch, Δt)
-
-    qa = randn(n, 1)
-    qb = randn(n, 1)
-    D₁, D₂ = _discrete_lagrangian_derivatives(loss, params(nn), qa, qb)
-
-    Ld(a, b) = Δt * sum(nn(vcat((a + b) / 2, (b - a) / Δt), params(nn)))
-    @test vec(D₁) ≈ Zygote.gradient(a -> Ld(a, qb), qa)[1]
-    @test vec(D₂) ≈ Zygote.gradient(b -> Ld(qa, b), qb)[1]
-end
-
-for n in (1, 2, 3)
-    test_discrete_lagrangian_derivatives(n)
-end
-
-@doc raw"""
-The loss must *be* the discrete Euler–Lagrange residual, normalised.
-
-Compare it against the same quantity assembled from `Zygote` gradients of ``L_d`` itself, which is
-the call `VariationalMidpointLoss` exists in order not to make. This asserts what the loss computes
-without needing a discrete trajectory to be solved for first.
-"""
-function test_variational_loss_is_the_del_residual(n::Integer; Δt = 0.05, batch_size = 4)
-    arch = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    loss = VariationalMidpointLoss(arch, Δt)
-
-    qₙ = randn(n, batch_size)
-    qₙ₊₁ = randn(n, batch_size)
-    qₙ₊₂ = randn(n, batch_size)
-
-    Ld(a, b) = Δt * sum(nn(vcat((a + b) / 2, (b - a) / Δt), params(nn)))
-    D₂ = reduce(hcat,
-        [Zygote.gradient(b -> Ld(qₙ[:, i], b), qₙ₊₁[:, i])[1] for i in 1:batch_size])
-    D₁ = reduce(hcat,
-        [Zygote.gradient(a -> Ld(a, qₙ₊₂[:, i]), qₙ₊₁[:, i])[1] for i in 1:batch_size])
-
-    @test loss(params(nn), vcat(qₙ, qₙ₊₁), qₙ₊₂) ≈ norm(D₂ + D₁) / norm(vcat(D₂, D₁))
-end
-
-for n in (1, 2, 3)
-    test_variational_loss_is_the_del_residual(n)
-end
-
-@doc raw"""
-The loss must vanish on a trajectory that satisfies the discrete Euler–Lagrange equations.
-
-At one degree of freedom the residual is scalar, so ``q_{n+2}`` is found by bracketing a sign change
-and bisecting. Bisection cannot diverge, which Newton on this residual does: a single overshoot puts
-the midpoint velocity deep into `tanh` saturation, where the gradient is constant to machine
-precision and a numerical Jacobian is exactly zero.
-"""
-function variational_next_position(loss, ps, qₙ, qₙ₊₁)
-    _, D₂ = _discrete_lagrangian_derivatives(loss, ps, qₙ, qₙ₊₁)
-    residual(q) = only(first(_discrete_lagrangian_derivatives(
-        loss, ps, qₙ₊₁, reshape([q], 1, 1))) + D₂)
-
-    centre = only(qₙ₊₁)
-    candidates = range(centre - 5.0, centre + 5.0; length = 1001)
-    values = residual.(candidates)
-    crossing = findfirst(i -> values[i] * values[i + 1] < 0, 1:(length(values) - 1))
-    crossing === nothing && return nothing
-
-    lo, hi = candidates[crossing], candidates[crossing + 1]
-    for _ in 1:60
-        mid = (lo + hi) / 2
-        if residual(lo) * residual(mid) ≤ 0
-            hi = mid
-        else
-            lo = mid
-        end
-    end
-    reshape([(lo + hi) / 2], 1, 1)
-end
-
-"""
-    a_discrete_trajectory(loss, ps; Δt, attempts)
-
-A triple `(qₙ, qₙ₊₁, qₙ₊₂)` that solves the discrete Euler–Lagrange equations of `loss`, or
-`nothing` if no starting pair out of `attempts` admits one.
-
-Not every starting pair does. The discrete Euler–Lagrange equations of an arbitrary neural-network
-Lagrangian are a root-finding problem that need not have a solution near the starting pair, so this
-searches rather than assuming.
-"""
-function a_discrete_trajectory(loss, ps; Δt = 0.05, attempts = 20)
-    for _ in 1:attempts
-        qₙ = randn(1, 1)
-        qₙ₊₁ = qₙ + Δt * randn(1, 1)
-        qₙ₊₂ = variational_next_position(loss, ps, qₙ, qₙ₊₁)
-        qₙ₊₂ === nothing || return (qₙ, qₙ₊₁, qₙ₊₂)
-    end
-    nothing
-end
-
-function test_variational_loss_vanishes_on_a_discrete_trajectory(; Δt = 0.05)
-    arch = LagrangianNeuralNetwork(2; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    loss = VariationalMidpointLoss(arch, Δt)
-
-    trajectory = a_discrete_trajectory(loss, params(nn); Δt = Δt)
-    @test trajectory !== nothing
-    trajectory === nothing && return
-    qₙ, qₙ₊₁, qₙ₊₂ = trajectory
-
-    @test loss(params(nn), vcat(qₙ, qₙ₊₁), qₙ₊₂) < 1.0e-10
-    # A position that does not solve the discrete Euler–Lagrange equations must not pass.
-    @test loss(params(nn), vcat(qₙ, qₙ₊₁), qₙ₊₂ .+ 0.1) > 1.0e-3
-end
-
-test_variational_loss_vanishes_on_a_discrete_trajectory()
-
-function test_variational_loss_gradient(n::Integer = 1)
-    arch = LagrangianNeuralNetwork(2n; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    loss = VariationalMidpointLoss(arch, 0.05)
-    input = randn(2n, 6)
-    output = randn(n, 6)
-
-    @test typeof(loss) <: NetworkLoss
-    @test loss(nn.model, params(nn), input, output) == loss(params(nn), input, output)
-
-    dp = Zygote.gradient(ps -> loss(ps, input, output), params(nn))[1]
-    @test typeof(dp) <: NetworkParameters
-    @test keys(dp) == keys(params(nn))
-    @test !all(iszero, dp.L1.W)
-end
-
-test_variational_loss_gradient()
-
-@doc raw"""
-Training must reduce the loss on a real trajectory.
-
-The data are the harmonic oscillator, ``q(t) = \cos(t)``, sampled at a fixed step. Asking the
-network for a Lagrangian whose discrete Euler–Lagrange equations that trajectory satisfies is the
-intended use of this loss, and it does not depend on the network's own initial draw the way a
-self-generated trajectory would.
-"""
-function test_variational_training(; Δt = 0.05, samples = 60, n_epochs = 20)
-    # The shape of the tests above, for the reason `test_symplectic_euler_training` gives.
-    arch = LagrangianNeuralNetwork(2; width = 4, nhidden = 1)
-    nn = NeuralNetwork(arch)
-    loss = VariationalMidpointLoss(arch, Δt)
-
-    t = range(0.0, Δt * (samples + 1); step = Δt)
-    q = cos.(t)
-    qₙ = reshape(q[1:samples], 1, samples)
-    qₙ₊₁ = reshape(q[2:(samples + 1)], 1, samples)
-    qₙ₊₂ = reshape(q[3:(samples + 2)], 1, samples)
-
-    dl = DataLoader(vcat(qₙ, qₙ₊₁), qₙ₊₂; suppress_info = true)
-    loss_array = train!(
-        nn, TrainingOptimizer(nn.params; algorithm = Adam()), dl, Batch(10), n_epochs;
-        loss, show_progress = false)
-    @test loss_array[end] < loss_array[1]
-end
-
-test_variational_training()
diff --git a/test/pullbacks/symbolic_hnn_pullback.jl b/test/pullbacks/symbolic_hnn_pullback.jl
deleted file mode 100644
index fbc2260fe..000000000
--- a/test/pullbacks/symbolic_hnn_pullback.jl
+++ /dev/null
@@ -1,173 +0,0 @@
-using GeometricMachineLearning
-using GeometricMachineLearning: HNNLoss, ZygotePullback, SigmoidActivation
-using AbstractNeuralNetworks: NeuralNetwork, params
-using LinearAlgebra: norm
-using NeuralNetworkParameters: NetworkParameters
-using Test
-import Random
-
-include("../helpers/eltypes.jl")
-include("../helpers/leaves.jl")
-
-# About 1.8e19 in `Float32` and 1.8e154 in `Float64`: the scale at which the product of the two
-# norms overflows and each norm does not.
-extreme_scale(::Type{Float32}) = Float32(1e20)
-extreme_scale(::Type{Float64}) = Float64(1e155)
-
-# The scale at which the seed scalar `1/(‖R‖‖Y‖)` underflows to zero while a gradient entry does
-# not: `Zygote`'s own `norm` rule underflows there too, so the row asserts only that the gradient
-# survives.
-extreme_scale_underflow(::Type{Float32}) = Float32(1e23)
-extreme_scale_underflow(::Type{Float64}) = Float64(1e170)
-
-# None of these depends on the element type, so every testset shares one build. The reference has a
-# loss of its own, so it does not depend on the code under test.
-const arch = StandardHamiltonianArchitecture(2, 4, 1, SigmoidActivation())
-const pullback = SymbolicPullback(arch)
-const reference = ZygotePullback(HNNLoss(arch))
-const hvf = GeometricMachineLearning.hamiltonian_vector_field(arch)
-
-# The symbolic pullback is the gradient of `HNNLoss` on the batch: the generated per-sample
-# gradients of the surrogate `½‖rᵢ‖²`, scaled by `1/(‖R‖‖Y‖)`. Two batch sizes: at `N = 1` the
-# sum is a single per-sample gradient and the scale still divides, while `N = 5` is where the
-# gradient of the per-sample sum and the gradient of the batch loss differ.
-@testset "SymbolicPullback of a Hamiltonian network is HNNLoss's gradient ($T)" for T in REAL_ELTYPES
-    for N in (1, 5)
-        Random.seed!(21)
-        nn = NeuralNetwork(arch, CPU(), T)
-        ps = params(nn)
-        X = randn(Random.Xoshiro(11), T, 2, N)
-        Y = randn(Random.Xoshiro(12), T, 2, N)
-
-        value, back = pullback(ps, nn.model, (X, Y))
-        g = back(one(value))
-        zygote_value, zygote_back = reference(ps, nn.model, (X, Y))
-        zg = zygote_back(one(zygote_value))
-
-        @test value isa T
-        @test value ≈ zygote_value rtol = T(√eps(T))
-        @test g isa NetworkParameters
-        @test keys(g) == keys(ps)
-        foreach_leaf(g, zg) do leaf, ref
-            @test eltype(leaf) === T
-            @test leaf ≈ ref rtol = T(√eps(T))
-        end
-
-        # the value is `HNNLoss` itself, so a wrong scale cannot hide behind a different value
-        @test value ≈ norm(hvf(X, ps) - Y) / norm(Y) rtol = T(√eps(T))
-    end
-end
-
-# The two norms are divided by separately, and each of them divides a gradient entry rather than
-# the seed: their product overflows where neither norm does, and `inv` of the product is then a zero
-# gradient. The scalar `1/(‖R‖‖Y‖)` is zero at the second scale, where the per-entry divisions keep
-# every leaf finite and non-zero. Neither scale is a power of two, so round-off shows.
-@testset "the scale divides by each norm ($T)" for T in REAL_ELTYPES
-    Random.seed!(22)
-    nn = NeuralNetwork(arch, CPU(), T)
-    ps = params(nn)
-    X = randn(Random.Xoshiro(13), T, 2, 5)
-
-    # (the scale, whether the two pullbacks still agree there): `Zygote`'s `norm` rule forms the
-    # seed scalar too, so underflowing it takes both paths down together
-    rows = ((extreme_scale(T), true), (extreme_scale_underflow(T), false))
-    for (scale, agreement) in rows
-        Y = randn(Random.Xoshiro(14), T, 2, 5) .* scale
-
-        value, back = pullback(ps, nn.model, (X, Y))
-        g = back(one(value))
-
-        @test value isa T
-        @test isfinite(value)
-        foreach_leaf(g) do leaf
-            @test eltype(leaf) === T
-            @test all(isfinite, leaf)
-            @test any(!iszero, leaf)
-        end
-
-        agreement || continue
-        zygote_value, zygote_back = reference(ps, nn.model, (X, Y))
-        zg = zygote_back(one(zygote_value))
-        @test value ≈ zygote_value rtol = T(√eps(T))
-        foreach_leaf(g, zg) do leaf, ref
-            @test leaf ≈ ref rtol = T(√eps(T))
-        end
-    end
-end
-
-# Data the network's own vector field generates: the residual is zero, the value is zero, and the
-# gradient is the zero set rather than `NaN` from an unguarded `inv`.
-@testset "a zero residual gives a zero gradient ($T)" for T in REAL_ELTYPES
-    Random.seed!(23)
-    nn = NeuralNetwork(arch, CPU(), T)
-    ps = params(nn)
-    X = randn(Random.Xoshiro(17), T, 2, 5)
-    Y = hvf(X, ps)
-
-    value, back = pullback(ps, nn.model, (X, Y))
-    g = back(one(value))
-
-    @test value isa T
-    @test value == zero(T)
-    @test g isa NetworkParameters
-    foreach_leaf(g) do leaf
-        @test eltype(leaf) === T
-        @test all(iszero, leaf)
-    end
-end
-
-# The gradient is linear in the seed, as `ZygotePullback`'s is, and has the parameters' element type
-# whatever the seed's.
-@testset "the seed scales the gradient ($T)" for T in REAL_ELTYPES
-    Random.seed!(24)
-    nn = NeuralNetwork(arch, CPU(), T)
-    ps = params(nn)
-    X = randn(Random.Xoshiro(19), T, 2, 5)
-    Y = randn(Random.Xoshiro(20), T, 2, 5)
-
-    value, back = pullback(ps, nn.model, (X, Y))
-    _, zygote_back = reference(ps, nn.model, (X, Y))
-
-    @test value isa T
-    g₁ = back(one(value))
-    g₂ = back(T(2))
-    zg₂ = zygote_back(T(2))
-
-    foreach_leaf(g₂, g₁, zg₂) do leaf, leaf₁, ref
-        @test eltype(leaf) === T
-        @test eltype(leaf₁) === T
-        @test leaf ≈ 2 .* leaf₁ rtol = T(√eps(T))
-        @test leaf ≈ ref rtol = T(√eps(T))
-    end
-
-    # a seed of the other float type gives the same gradient, in the parameters' element type
-    g_other = back(T === Float32 ? 1.0 : 1.0f0)
-    foreach_leaf(g_other, g₁) do leaf, leaf₁
-        @test eltype(leaf) === T
-        @test leaf ≈ leaf₁ rtol = T(√eps(T))
-    end
-end
-
-# Data of the other float type: the gradient keeps the parameters' element type, as
-# `ZygotePullback`'s does, so `train!` can step with it. The two gradients agree to the precision of
-# the coarser type.
-@testset "data of the other float type gives a gradient of the parameters' type ($T)" for T in REAL_ELTYPES
-    S = T === Float32 ? Float64 : Float32
-    Random.seed!(26)
-    nn = NeuralNetwork(arch, CPU(), T)
-    ps = params(nn)
-    X = randn(Random.Xoshiro(27), S, 2, 5)
-    Y = randn(Random.Xoshiro(28), S, 2, 5)
-
-    value, back = pullback(ps, nn.model, (X, Y))
-    g = back(one(value))
-    zygote_value, zygote_back = reference(ps, nn.model, (X, Y))
-    zg = zygote_back(one(zygote_value))
-
-    @test g isa NetworkParameters
-    foreach_leaf(g, zg) do leaf, ref
-        @test eltype(leaf) === T
-        @test eltype(ref) === T
-        @test leaf ≈ ref rtol = √eps(Float32)
-    end
-end
diff --git a/test/quality/aqua.jl b/test/quality/aqua.jl
index ee504e437..2576d7c2e 100644
--- a/test/quality/aqua.jl
+++ b/test/quality/aqua.jl
@@ -18,16 +18,17 @@ using Test
 #
 # `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 ArrayLayouts',
-# FillArrays', Symbolics', GeometricOptimizers', BandedMatrices' and BlockArrays' versions, and with
-# *which* of them a given process had loaded. `detect_ambiguities` on this package, in this process,
+# FillArrays', GeometricOptimizers', BandedMatrices' and BlockArrays' versions, and with *which* of
+# them a given process had loaded. `detect_ambiguities` on this package, in this process,
 # asks only about the methods it defines; the count is zero, and no method of this package claims an
 # argument position wholesale on `PoissonTensor`, which is what keeps it there. That type is an
 # `AbstractMatrix{T}`, so such a method meets every special-array method another package writes
 # against `AbstractMatrix`:
 #
 #   `*(𝕁::PoissonTensor{T}, v::Strided…)` takes a `Strided…` right-hand side, which excludes each
-#       special array type `X` that ArrayLayouts, FillArrays, Symbolics and GeometricOptimizers
-#       define a `*(::AbstractMatrix, ::X)` for, and still covers everything the package builds.
+#       special array type `X` that ArrayLayouts, FillArrays, GeometricOptimizers and
+#       Symbolics (which HamiltonianNeuralNetworks.jl loads) define a `*(::AbstractMatrix, ::X)`
+#       for, and still covers everything the package builds.
 #   `getindex(𝕁::PoissonTensor, i::Int, j::Int)` types its indices, which excludes the `Block`,
 #       `BlockIndex` and `BandRangeType` methods BandedMatrices and BlockArrays add to
 #       `AbstractMatrix`. Neither of those two loads with this package alone, so an untyped index
@@ -47,7 +48,7 @@ end
 # `all_qualified_accesses_via_owners` and `no_self_qualified_accesses`. Two more run in the testset
 # below with a list of names, and `no_implicit_imports` is off:
 #
-#   no_implicit_imports              37 names arrive through a bare `using`, across thirteen packages
+#   no_implicit_imports              33 names arrive through a bare `using`, across twelve packages
 #                                    -- `norm`, `@kernel`, `rrule`, `pullback` and the rest. Naming
 #                                    every one is a change to the whole module header and a judgement
 #                                    per name, not a by-product of a dead-code pass.
@@ -73,14 +74,11 @@ test_explicit_imports(GeometricMachineLearning;
 # those packages do not declare public; each is deliberate, and the fix for each is upstream
 # declaring it. No `GeometricOptimizers` name is on either list.
 const NON_PUBLIC_QUALIFIED = (
-    :Initializer, :TanhActivation, :dim,                         # AbstractNeuralNetworks
-    :Jacobian, :symbolic_parameter_gradient,                     # SymbolicNeuralNetworks
-    :variables,                                                  # Symbolics
-    :copyto!, :zeros,                                            # KernelAbstractions
+    :Initializer, :dim,                                          # AbstractNeuralNetworks
+    :zeros,                                                      # KernelAbstractions
     :jacobian)                                                   # ForwardDiff
 const NON_PUBLIC_IMPORTS = (
     :AbstractExplicitLayer, :Architecture, :_compute_loss, :dim, :model, # AbstractNeuralNetworks
-    :derivative,                                                         # SymbolicNeuralNetworks
     :load, :save,                                                        # NeuralNetworkParameters
     :_vcat, :assign_q_and_p)                       # this package, imported by `GPUArraysCoreExt`
 
diff --git a/test/quality/jet.jl b/test/quality/jet.jl
index ad5934a0d..baccabe85 100644
--- a/test/quality/jet.jl
+++ b/test/quality/jet.jl
@@ -127,9 +127,9 @@ const A3{T} = Array{T, 3}
                     Batch{:FeedForward},
                     Indices);
                 target_modules = TARGET)))
-            # input and output: test/data_loader/data_loader_for_input_and_output.jl
-            # (`Float32`, `Batch(10, 1, 1)`) and
-            # test/architectures/lagrangian_neural_network_tests.jl (`Float64`, `Batch(10)`)
+            # 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);
diff --git a/test/runtests.jl b/test/runtests.jl
index ca1821161..8a25eef3d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -48,7 +48,6 @@ if "core" in GROUPS
     @safetestset "Linear Symplectic Transformer" include("architectures/linear_symplectic_transformer.jl")
     @safetestset "Symplectic Transformer chain construction" include("architectures/symplectic_transformer_chain.jl")
     @safetestset "Check parameterlength" include("architectures/check_parameterlengths.jl")
-    @safetestset "Hamiltonian Neural Network" include("architectures/hamiltonian_neural_network_tests.jl")
     @safetestset "SympNet integrator" include("architectures/sympnet_integrator.jl")
     @safetestset "PSD tests" include("architectures/psd_architecture_tests.jl")
     @safetestset "SymplecticAutoencoder tests" include("architectures/symplectic_autoencoder_tests.jl")
@@ -58,7 +57,6 @@ if "core" in GROUPS
     @safetestset "The forward pass of every exported architecture infers" include("architectures/forward_inference.jl")
     @safetestset "Check reduced model" include("reduced_system/reduced_system.jl")
     @safetestset "_norm keeps the element type of its argument" include("utils.jl")
-    @safetestset "Symplectic Euler and variational midpoint losses" include("loss/training_method_losses.jl")
     @safetestset "Test NetworkLoss + Optimizer" include("integration/losses_and_optimization.jl")
     @safetestset "ReducedLoss through train!" include("integration/reduced_loss_optimization.jl")
     @safetestset "Optimizer #2" include("optimizers/optimization_step.jl")
@@ -71,7 +69,6 @@ if "core" in GROUPS
     @safetestset "_mul: NetworkParameters gradient structure" include("layers/double_multiplication_network_parameters_gradient.jl")
     @safetestset "Symplectic attention: NetworkParameters gradient structure" include("layers/symplectic_attention_network_parameters_gradient.jl")
     @safetestset "ZygotePullback returns a NetworkParameters" include("pullbacks/zygote_pullback.jl")
-    @safetestset "SymbolicPullback of HNNLoss is the gradient of the batched loss" include("pullbacks/symbolic_hnn_pullback.jl")
     @safetestset "AutoEncoderLoss and a loss's NetworkParameters" include("loss/losses.jl")
     @safetestset "map_to_cpu" include("map_to_cpu.jl")
     @safetestset "changebackend" include("integration/changebackend_tests.jl")
@@ -101,7 +98,6 @@ if "core" in GROUPS
     @safetestset "docs/make.jl stops on a bad @ref before makedocs" include("quality/docs_make_checks_references.jl")
 end
 if "slow" in GROUPS
-    @safetestset "Lagrangian Neural Network" include("architectures/lagrangian_neural_network_tests.jl")
     @safetestset "Test symplecticity of the sympnet upscaling layer" include("layers/sympnet_upscaling.jl")
     @safetestset "train! trains every architecture of the manual" include("integration/train_manual_architectures.jl")
 end