michakraus · Oct 26, 2026 · Oct 26, 2026 · Oct 26, 2026 · Oct 26, 2026 · Oct 26, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md
@@ -68,18 +68,17 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
 
 ## Found late
 
-### K17 · `assign_q_and_p` and `PoissonTensor * v` drop the rows from `2N + 1` on without an error.
+### K17 · `assign_q_and_p` drops the rows from `2N + 1` on without an error.
 
 - **location:** `src/kernels/assign_q_and_p.jl`
 - **kind:** defect
 - **found:** #345, review of 2026-10-15; the behaviour is the same on `origin/main` at `2a63f26a`
 - **evidence:**
   `assign_q_and_p(x, N)` takes `x[1:N, …]` for `q` and the next `N` rows for `p`, and checks nothing
-  about `size(x, 1)`. An input with `size(x, 1) > 2N` loses its trailing rows silently, and
-  `PoissonTensor * v` reaches the same method.
+  about `size(x, 1)`. An input with `size(x, 1) > 2N` loses its trailing rows silently.
 
   Closing it means an error for `size(x, 1) != 2N`. That changes the public behaviour of
-  `assign_q_and_p` and of `PoissonTensor * v`, so it is a change of its own and not a cleanup.
+  `assign_q_and_p`, so it is a change of its own and not a cleanup.
 
 ### K19 · In `Float32` the Cayley kernel returns non-finite slices for a skew `C` of norm 1e8, where no method is accurate.
 
@@ -213,25 +212,24 @@ What is known to be wrong in GeometricMachineLearning and is not fixed.
 
 ## The test suite
 
-### K3 · `test/quality/jet.jl` does not see an unstable kernel-launch argument, nor two barriers in the Poisson-tensor kernel
+### K3 · `test/quality/jet.jl` does not see an unstable kernel-launch argument
 
 - **location:** `test/quality/jet.jl`
 - **kind:** missing test
 - **found:** 2026-10-01
 - **evidence:**
-  Julia 1.13.1, JET 0.12.2, `mutate.jl <repository> … quality/jet.jl`. Three mutants that put
-  `Base.inferencebarrier` on an argument of a kernel launch SURVIVED (142 pass, 3 broken):
-  `kernel!(c, Base.inferencebarrier(a), b, ndrange = size(c))` in
-  `src/kernels/tensor_tensor_mul.jl`; `augment_zeros!(zero_tensor,
-  Base.inferencebarrier(output_diff), seq_length, …)` in `src/data_loader/tensor_assign.jl`; and
-  `skew_mat_mul_k!(C, S, Base.inferencebarrier(B), n, …)` in `src/kernels/mat_tensor_mul.jl`. At
-  `target_modules = (GeometricMachineLearning,)` a barrier on a launch argument gives 0 reports,
-  and the same barrier in an `@assert` of the same launcher gives 3: the kernel functor takes
-  `args...`, so the dispatch sits in KernelAbstractions' frames, which the filter drops. In
-  `assign_ones_for_poisson_tensor_kernel!` (`src/arrays/poisson_tensor.jl:95`) two mutants
-  SURVIVED: a barrier on `n` in `map_index_for_poisson_tensor(i, n)`, and a barrier on the written
-  array `J`; a barrier on the stored value is CAUGHT. JET never reports a dynamic dispatch in the
-  kernel's own frame. A barrier on a value inside each kernel body is CAUGHT.
+  Julia 1.13.2, `mutate.jl <repository> <file> <from> <to> quality/jet.jl`, one mutant per run.
+  Three mutants that put `Base.inferencebarrier` on an argument of a kernel launch SURVIVED, each
+  with 50 pass, 0 fail: `kernel!(c, Base.inferencebarrier(a), b, ndrange = size(c))` in
+  `src/kernels/tensor_tensor_mul.jl`; `tensor_mat_skew_sym_assign_k!(C,
+  Base.inferencebarrier(Z), A, ndrange = size(C))` in
+  `src/kernels/inverses/tensor_mat_skew_sym_assign.jl`; and `gather_kernel!(backend)(gathered,
+  Base.inferencebarrier(data), indices, offset, ndrange = size(gathered))` in
+  `src/data_loader/batch.jl`. At `target_modules = (GeometricMachineLearning,)` a barrier on a
+  launch argument gives 0 reports, and the same barrier in an `@assert` of the same launcher is
+  CAUGHT with 3 reports (47 pass, 3 fail): the kernel functor takes `args...`, so the dispatch sits
+  in KernelAbstractions' frames, which the filter drops. A barrier on a value inside the kernel
+  body of `tensor_tensor_mul_kernel!` is CAUGHT (47 pass, 3 fail).
 
 ### K10 · The two training tests of `adam_with_learning_rate_decay.jl` run in `Float64` only.
 

diff --git a/Project.toml b/Project.toml
@@ -12,39 +12,28 @@ GeometricBase = "9a0b12b7-583b-4f04-aa1f-d8551b6addc9"
 GeometricEquations = "c85262ba-a08a-430a-b926-d29770767bf2"
 GeometricOptimizers = "fc236c15-5557-4942-aa65-b650f329279e"
 GeometricSolutions = "7843afe4-64f4-4df4-9231-049495c56661"
-InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 NeuralNetworkParameters = "67f4d93a-60e9-472b-8cdd-1ccf6005724a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
-[weakdeps]
-GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
-
-[extensions]
-GPUArraysCoreExt = "GPUArraysCore"
-
 [compat]
 AbstractNeuralNetworks = "0.9.2"
 ChainRulesCore = "1"
 DiffResults = "1.1"
 ForwardDiff = "1"
-GPUArraysCore = "0.2"
 GeometricBase = "0.15"
 GeometricEquations = "0.21.5"
 GeometricOptimizers = "0.9"
 GeometricSolutions = "0.6.6"
-InteractiveUtils = "1"
 KernelAbstractions = "0.9"
 LinearAlgebra = "1"
 NNlib = "0.8, 0.9"
 NeuralNetworkParameters = "0.5"
 ProgressMeter = "1"
 Random = "1"
-TimerOutputs = "0.5, 1"
 Zygote = "0.7"
 julia = "1.12"
diff --git a/README.md b/README.md
@@ -55,7 +55,7 @@ g_loss_array = train!(g_nn, g_opt, dl, Batch(batch_size), nepochs)
 # plot the result
 ics = (q=qp_data.q[:,1], p=qp_data.p[:,1])
 const steps_to_plot = 200
-g_trajectory = iterate(g_nn, ics; n_points = steps_to_plot)
+g_trajectory = integrate(g_nn, ics; n_points = steps_to_plot)
 fig = Figure()
 ax = Axis(fig[1, 1]; xlabel = "q", ylabel = "p")
 lines!(ax, vec(qp_data.q')[1:steps_to_plot], vec(qp_data.p')[1:steps_to_plot]; label = "training data")

diff --git a/docs/check_references.jl b/docs/check_references.jl
@@ -12,6 +12,7 @@
 #     julia --project=docs docs/check_references.jl
 
 using GeometricMachineLearning
+using GeometricOptimizers: solve! # the manual's `@docs solve!` entry resolves through `Main`, and this package does not export it
 using AbstractNeuralNetworks # signatures in `@docs` blocks name types from here
 using Documenter
 using Documenter: DocSystem

diff --git a/docs/make.jl b/docs/make.jl
@@ -4,6 +4,7 @@
 include(joinpath(@__DIR__, "check_references.jl"))
 
 using GeometricMachineLearning
+using GeometricOptimizers: solve! # the manual's `@docs solve!` entry names it, and this package does not export it
 using HDF5
 using AbstractNeuralNetworks
 using Documenter

diff --git a/docs/src/architectures/neural_network_integrators.md b/docs/src/architectures/neural_network_integrators.md
@@ -58,9 +58,9 @@ Main.remark(raw"For standard multi-step methods (that are not neural network-bas
 NeuralNetworkIntegrator
 ResNet
 GeometricMachineLearning.ResNetLayer
-iterate(::NeuralNetwork{<:NeuralNetworkIntegrator}, ::NamedTuple{(:q, :p), Tuple{AT, AT}}) where {T, AT<:AbstractVector{T}}
+integrate(::NeuralNetwork{<:NeuralNetworkIntegrator}, ::NamedTuple{(:q, :p), Tuple{AT, AT}}) where {T, AT<:AbstractVector{T}}
 TransformerIntegrator
-iterate(::NeuralNetwork{<:TransformerIntegrator}, ::NamedTuple{(:q, :p), Tuple{AT, AT}}) where {T, AT<:AbstractMatrix{T}}
+integrate(::NeuralNetwork{<:TransformerIntegrator}, ::NamedTuple{(:q, :p), Tuple{AT, AT}}) where {T, AT<:AbstractMatrix{T}}
 ```
 
 ```@raw latex

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -52,7 +52,7 @@ Many layers have been adapted in order to be used for problems in scientific mac
 
 `GeometricMachineLearning` allocates and computes through `KernelAbstractions.jl` [churavy2020kernel](@cite), so its layers and architectures are written against any backend that package supports: `CUDA.jl` [besard2018juliagpu](@cite), `AMDGPU.jl`, `Metal.jl` [besard2022metal](@cite) and `oneAPI.jl` [besard2022one](@cite).
 
-**Little of that is tested.** `test/devices/metal.jl` runs two testsets on an Apple GPU, in the `metal` test group, which the `Metal` job and the `macOS-latest` job of CI run: the products of a `PoissonTensor` with a wrapped device array; and the products of `_mul` — a dense and a `SymmetricMatrix` weight from the left, a dense one from the right — together with the `tensor_tensor_mul` kernel. A green test matrix says nothing about the rest of the GPU path. What is written below was measured by hand on an Apple M4 Max through `Metal.jl` 1.11.1, in `Float32`, because Apple GPUs have no `Float64` at all, and against `GeometricOptimizers` 0.9.0.
+**Little of that is tested.** `test/devices/metal.jl` runs two testsets on an Apple GPU, in the `metal` test group, which the `Metal` job and the `macOS-latest` job of CI run: the products of `_mul` — a dense and a `SymmetricMatrix` weight from the left, a dense one from the right — together with the `tensor_tensor_mul` kernel; and `tensor_cayley` with the layer that uses it. A green test matrix says nothing about the rest of the GPU path. What is written below was measured by hand on an Apple M4 Max through `Metal.jl` 1.11.1, in `Float32`, because Apple GPUs have no `Float64` at all, and against `GeometricOptimizers` 0.9.0.
 
 Everything tried on that device ran. Constructed, applied to both a matrix and a 3-tensor, and copied back with `map_to_cpu`: `GSympNet`, `LASympNet`, `StandardTransformerIntegrator`, `LinearSymplecticTransformer`, `SymplecticTransformer` at both its default `transformer_dim` and an upscaling one, `VolumePreservingFeedForward`, `VolumePreservingTransformer`, `SymplecticAutoencoder`, `PSDArch` and `Transformer(…; Stiefel = true)`. The product of a `StrictlyLowerTriangular` matrix with a 3-tensor runs. `train!` trains a `GSympNet` for two epochs with its data in `MtlArray`s. `ClassificationTransformer` was not tried.
 

diff --git a/docs/src/structure_preservation/symplecticity.md b/docs/src/structure_preservation/symplecticity.md
@@ -111,11 +111,6 @@ It is important to note that symplecticity is a very strong property[^2] that ma
 
 [^2]: Symplecticity imposes in general greater restrictions on the flow map than e.g. conservation of energy. The famous *Ge-Marsden theorem* [ge1988lie](@cite) says that one cannot achieve preservation of energy and preservation of symplecticity at the same time, unless the numerical method is the exact integral of the Hamiltonian ODE. In practice we hence almost always choose a symplectic over an energy-preserving scheme.
 
-## Library Functions
-```@docs
-PoissonTensor
-```
-
 ```@raw latex
 \begin{comment}
 ```

diff --git a/docs/src/tutorials/linear_symplectic_transformer.md b/docs/src/tutorials/linear_symplectic_transformer.md
@@ -156,9 +156,9 @@ function make_validation_plot(n_steps = n_steps; theme = :dark)
         xlabelcolor = textcolor,
         ylabelcolor = textcolor,
     )
-    prediction_standard = iterate(nn_standard, init_con; n_points = n_steps, prediction_window = seq_length)
-    prediction_symplectic = iterate(nn_symplectic, init_con; n_points = n_steps, prediction_window = seq_length)
-    prediction_sympnet = iterate(nn_sympnet, init_con_ff; n_points = n_steps)
+    prediction_standard = integrate(nn_standard, init_con; n_points = n_steps, prediction_window = seq_length)
+    prediction_symplectic = integrate(nn_symplectic, init_con; n_points = n_steps, prediction_window = seq_length)
+    prediction_sympnet = integrate(nn_sympnet, init_con_ff; n_points = n_steps)
 
     # we use linewidth  = 2
     lines!(ax, dl.input.q[1, 1:n_steps, index]; color = mblue, label = "Implicit midpoint", linewidth = 2)

diff --git a/docs/src/tutorials/matrix_softmax.md b/docs/src/tutorials/matrix_softmax.md
@@ -161,8 +161,8 @@ function make_validation_plot(n_steps = n_steps; theme = :dark)
         xlabelcolor = textcolor,
         ylabelcolor = textcolor,
     )
-    prediction_vector = iterate(nn1, init_con; n_points = n_steps, prediction_window = seq_length)
-    prediction_matrix = iterate(nn2, init_con; n_points = n_steps, prediction_window = seq_length)
+    prediction_vector = integrate(nn1, init_con; n_points = n_steps, prediction_window = seq_length)
+    prediction_matrix = integrate(nn2, init_con; n_points = n_steps, prediction_window = seq_length)
 
     # we use linewidth  = 2
     lines!(ax, dl.input.q[1, 1:n_steps, index]; color = mblue, label = "Implicit midpoint", linewidth = 2)

diff --git a/docs/src/tutorials/softmax_comparison.md b/docs/src/tutorials/softmax_comparison.md
@@ -157,8 +157,8 @@ function make_validation_plot(n_steps = n_steps; theme = :dark)
         xlabelcolor = textcolor,
         ylabelcolor = textcolor,
     )
-    prediction_vector = iterate(nn1, init_con; n_points = n_steps, prediction_window = seq_length)
-    prediction_matrix = iterate(nn2, init_con; n_points = n_steps, prediction_window = seq_length)
+    prediction_vector = integrate(nn1, init_con; n_points = n_steps, prediction_window = seq_length)
+    prediction_matrix = integrate(nn2, init_con; n_points = n_steps, prediction_window = seq_length)
 
     # we use linewidth  = 2
     lines!(ax, dl.input.q[1, 1:n_steps, index]; color = mblue, label = "Implicit midpoint", linewidth = 2)

diff --git a/docs/src/tutorials/symplectic_autoencoder.md b/docs/src/tutorials/symplectic_autoencoder.md
@@ -77,6 +77,7 @@ The training data can very easily be obtained by using the packages [`GeometricP
 ```@example toda_lattice
 using GeometricIntegrators: integrate, ImplicitMidpoint
 using GeometricMachineLearning 
+using GeometricMachineLearning: solve! # the PSD architecture is solved, not trained
 import Random # hide
 
 pr = tl.hodeproblem(; timespan = (0.0, 800.))
@@ -307,8 +308,8 @@ We can now evaluate the solution:
 ```@example toda_lattice
 integrator_nn = load(NeuralNetwork, "integrator_parameters.h5", integrator_architecture) # hide
 ics = encoder(sae_nn_cpu)((q = dl.input.q[:, 1:seq_length, 1], p = dl.input.p[:, 1:seq_length, 1])) # hide
-iterate(mtc(integrator_nn), ics; n_points = length(sol.t), prediction_window = seq_length) # hide
-@time "time stepping with transformer" time_series = iterate(mtc(integrator_nn), ics; n_points = length(sol.t), prediction_window = seq_length)
+integrate(mtc(integrator_nn), ics; n_points = length(sol.t), prediction_window = seq_length) # hide
+@time "time stepping with transformer" time_series = integrate(mtc(integrator_nn), ics; n_points = length(sol.t), prediction_window = seq_length)
 nothing # hide
 ```
 
@@ -335,7 +336,7 @@ function plot_transformer_validation!(fig, coordinates, t_steps::Integer=100; th
     lines!(ax_val, Ω, sae_rs.decoder((q = sol_sae_reduced.dataser.q[t_steps], p = sol_sae_reduced.dataser.p[t_steps])).q, 
         label = rich("SAE + Implicit Midpoint"; color = textcolor), color = mgreen)
 
-    time_series = iterate(mtc(integrator_nn), ics; n_points = t_steps, prediction_window = seq_length)
+    time_series = integrate(mtc(integrator_nn), ics; n_points = t_steps, prediction_window = seq_length)
     # prediction = (q = time_series.q[:, end], p = time_series.p[:, end])
     prediction = (q = time_series.q[:, end], p = time_series.p[:, end])
     sol = decoder(sae_nn_cpu)(prediction)
@@ -497,8 +498,8 @@ nothing # hide
 We again integrate the system and then plot the result:
 
 ```@example toda_lattice
-iterate(mtc(integrator_nn2), ics; n_points = length(sol.t), prediction_window = seq_length) # hide
-@time "time stepping with transformer" time_series2 = iterate(mtc(integrator_nn2), ics; n_points = length(sol.t), prediction_window = seq_length)
+integrate(mtc(integrator_nn2), ics; n_points = length(sol.t), prediction_window = seq_length) # hide
+@time "time stepping with transformer" time_series2 = integrate(mtc(integrator_nn2), ics; n_points = length(sol.t), prediction_window = seq_length)
 nothing # hide
 ```
 
@@ -536,7 +537,7 @@ function plot_validation!(fig, coordinates::Tuple, t_steps::Integer=100; theme =
     lines!(ax_val, Ω, psd_rs2.decoder((q = sol_psd_reduced2.dataser.q[t_steps], p = sol_psd_reduced2.dataser.p[t_steps])).q, 
         label = rich("PSD + Implicit Midpoint"; color = textcolor), color = morange)
 
-    time_series2 = iterate(mtc(integrator_nn2), ics; n_points = t_steps, prediction_window = seq_length)
+    time_series2 = integrate(mtc(integrator_nn2), ics; n_points = t_steps, prediction_window = seq_length)
     # prediction = (q = time_series.q[:, end], p = time_series.p[:, end])
     prediction2 = (q = time_series2.q[:, end], p = time_series2.p[:, end])
     sol = decoder(psd_nn2)(prediction2)

diff --git a/docs/src/tutorials/symplectic_transformer.md b/docs/src/tutorials/symplectic_transformer.md
@@ -156,9 +156,9 @@ function make_validation_plot(n_steps = n_steps; theme = :dark)
         xlabelcolor = textcolor,
         ylabelcolor = textcolor,
     )
-    prediction_standard = iterate(nn_standard, init_con; n_points = n_steps, prediction_window = seq_length)
-    prediction_symplectic = iterate(nn_symplectic, init_con; n_points = n_steps, prediction_window = seq_length)
-    prediction_sympnet = iterate(nn_sympnet, init_con_ff; n_points = n_steps)
+    prediction_standard = integrate(nn_standard, init_con; n_points = n_steps, prediction_window = seq_length)
+    prediction_symplectic = integrate(nn_symplectic, init_con; n_points = n_steps, prediction_window = seq_length)
+    prediction_sympnet = integrate(nn_sympnet, init_con_ff; n_points = n_steps)
 
     # we use linewidth  = 2
     lines!(ax, dl.input.q[1, 1:n_steps, index]; color = mblue, label = "Implicit midpoint", linewidth = 2)

diff --git a/docs/src/tutorials/sympnet_tutorial.md b/docs/src/tutorials/sympnet_tutorial.md
@@ -169,16 +169,16 @@ nothing # hide
 ![Training loss for SympNet.](sympnet_training_loss_light.png)
 ![Training loss for SympNet.](sympnet_training_loss_dark.png)
 
-Now we can make a prediction. We compare the initial data with a prediction starting from the same phase space point using the function [`GeometricMachineLearning.iterate`](@ref):
+Now we can make a prediction. We compare the initial data with a prediction starting from the same phase space point using the function [`GeometricMachineLearning.integrate`](@ref):
 
 ```@example sympnet
 ics = (q=dl.input.q[:, 1, 1], p=dl.input.p[:, 1, 1])
 
 steps_to_plot = 200
 
 #predictions
-la_trajectory = iterate(la_nn, ics; n_points = steps_to_plot)
-g_trajectory =  iterate(g_nn, ics; n_points = steps_to_plot)
+la_trajectory = integrate(la_nn, ics; n_points = steps_to_plot)
+g_trajectory =  integrate(g_nn, ics; n_points = steps_to_plot)
 nothing # hide
 ```
 
@@ -297,7 +297,7 @@ nothing # hide
 And we see that the loss is significantly lower than for both the ``LA``-SympNet and the ``G``-SympNet. We can also plot the prediction:
 
 ```@example sympnet
-rn_trajectory = iterate(rn_nn, ics; n_points = steps_to_plot)
+rn_trajectory = integrate(rn_nn, ics; n_points = steps_to_plot)
 nothing # hide
 ```
 
@@ -344,9 +344,9 @@ We see that the ResNet is slowly gaining energy which consitutes unphysical beha
 steps_to_plot = 400
 
 #predictions
-la_trajectory = iterate(la_nn, ics; n_points = steps_to_plot)
-g_trajectory =  iterate(g_nn, ics; n_points = steps_to_plot)
-rn_trajectory = iterate(rn_nn, ics; n_points = steps_to_plot)
+la_trajectory = integrate(la_nn, ics; n_points = steps_to_plot)
+g_trajectory =  integrate(g_nn, ics; n_points = steps_to_plot)
+rn_trajectory = integrate(rn_nn, ics; n_points = steps_to_plot)
 nothing # hide
 ```
 

diff --git a/docs/src/tutorials/volume_preserving_attention.md b/docs/src/tutorials/volume_preserving_attention.md
@@ -231,9 +231,9 @@ textcolor = theme == :dark ? :white : :black # hide
         xlabelcolor = textcolor,
         ylabelcolor = textcolor,
         )
-    validation_skew = iterate(nn_skew, initial_condition; n_points = n_points, prediction_window = 1)
-    validation_arb  = iterate(nn_arb,  initial_condition; n_points = n_points, prediction_window = 1)
-    validation_comp = iterate(nn_comp, initial_condition[:, 1]; n_points = n_points)
+    validation_skew = integrate(nn_skew, initial_condition; n_points = n_points, prediction_window = 1)
+    validation_arb  = integrate(nn_arb,  initial_condition; n_points = n_points, prediction_window = 1)
+    validation_comp = integrate(nn_comp, initial_condition[:, 1]; n_points = n_points)
 
     p2 = type == :cos ? lines!(dl.input[1, 1:n_points, 2], color = mpurple, label = "reference") : plot(dl.input[1, 1:n_points, 1], color = mpurple, label = "reference")