diff --git a/.gitignore b/.gitignore
index fe386cc..2799016 100644
--- a/.gitignore
+++ b/.gitignore
@@ -38,3 +38,6 @@ Presentations/*.png
 Presentations/*.pdf
 runs/*.h5
 runs/*.pdf
+
+# Output of the scripts
+/results
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6506f3b..bb9dabe 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,8 +20,36 @@ holds the changes since the last release until that release is cut.
   sees it without the export, so `Array` keeps its meaning in a caller's code; what changes is the
   export list that `names(ReducedComplexityModeling)` returns.
 
+- `ReducedBasisModel` has type parameters for the reduction scheme, the element type of the basis
+  and the type of the snapshot, and its three fields are concretely typed. They were `Any`. The
+  constructor takes an `AbstractMatrix` snapshot: a `Vector` throws a `MethodError`. The basis has
+  the floating-point type of the snapshot's element type, so a `Float32` snapshot gives a `Float32`
+  basis, and an integer snapshot gives a `Float64` basis.
+
+- `GeometricIntegratorData` and `GeometricIntegratorEnsembleData` have type parameters, for the
+  solution type, the element type of the snapshot and, for the ensemble, the dimension of the
+  snapshot, and their fields are concretely typed. The fields `snapshot` were `Matrix` and `Array`,
+  and the field `solutions` was a `Vector{GeometricSolution}`. A call that names the type without
+  parameters, such as `::GeometricIntegratorData`, still matches.
+
+- `reduce!(model, k, CotangentLiftEVD())` takes the eigenvectors of the Gram matrix `q qᵀ + p pᵀ`.
+  It took the left singular vectors of `[q p]` before. In exact arithmetic the two bases span the
+  same space where the singular values have a gap at `k`, and the columns can differ in sign and,
+  for equal singular values, in order, so a result that compares the basis entry by entry changes.
+  Where the `k`-th and the next singular value are equal, the two methods can pick different
+  subspaces, and either one is a valid choice. The Gram matrix squares the condition
+  number: the eigenvectors resolve only the singular vectors whose singular values exceed about
+  `sqrt(eps)` times the largest. Where those directions matter, use `CotangentLiftSVD()`, which
+  takes the singular vectors of `[q p]` directly. The eigendecomposition of the `n × n` Gram
+  matrix takes `O(n³)` time and `O(n²)` memory, for a snapshot with `2n` rows.
+
 ### Bug Fixes
 
+- `Parameter(name, minimum, maximum, n)` takes bounds that are a few ULPs apart (K27). A point of
+  `LinRange` that rounds outside the bounds threw an `AssertionError`. It is moved onto the
+  nearest bound now. `Parameter(:a, 0.1, nextfloat(0.1), 22)` threw, and 368 of the counts in
+  `2:500` did for these bounds.
+
 - `Parameter` takes a `Symbol` or a string name in every constructor form. The call
   `Parameter("μ", minimum, maximum, n)` returns a `Parameter` whose samples are those of
   `LinRange(minimum, maximum, n)` when the bounds differ, without duplicates. Equal bounds and a
@@ -60,8 +88,14 @@ holds the changes since the last release until that release is cut.
 
 - `reduce!(model, k, CotangentLiftEVD())` (K17) and `reduce!(model, k, ComplexSVD())` (K18) no
   longer throw on every call. Each stores a `2n × 2k` basis in `model.reduced_basis`, orthonormal
-  and symplectic (`AᵀJA = J`). Neither algorithm type is exported: a caller reaches them as
-  `ReducedComplexityModeling.CotangentLiftEVD` and `ReducedComplexityModeling.ComplexSVD`.
+  and symplectic (`AᵀJA = J`). `CotangentLiftSVD()`, `CotangentLiftEVD()` and `ComplexSVD()` take
+  a real snapshot: a complex one throws a `MethodError`, because a complex snapshot gives a basis
+  that does not satisfy `AᵀJA = J`.
+  None of the three algorithm types is exported: a caller reaches them as
+  `ReducedComplexityModeling.CotangentLiftSVD`, `ReducedComplexityModeling.CotangentLiftEVD` and
+  `ReducedComplexityModeling.ComplexSVD`. `CotangentLiftEVD()` scales the snapshot by its largest
+  entry before it forms the Gram matrix, so a snapshot of very large or very small entries gives a
+  finite basis.
 
 - Help text: the module and every exported name have a docstring with text. `Parameter`,
   `ParameterSampler` and `CartesianParameterSampler` replace their empty docstrings, and `POD`,
@@ -75,6 +109,34 @@ holds the changes since the last release until that release is cut.
   `path = "/"`, `h5load(Parameter, file)` throws a `KeyError`. The docstring of
   `ReducedBasisModel` no longer names packages that are not dependencies of this package.
 
+- `reduce!` takes a `k` above the columns of the thin SVD (K22), for `POD()` and for
+  `CotangentLiftSVD()`, `CotangentLiftEVD()` and `ComplexSVD()`, the last three with any `k ≤ n`.
+  It threw a `BoundsError`; the basis is completed with orthonormal vectors now, and stays
+  symplectic. The basis of the three symplectic algorithms has the floating-point type of the
+  snapshot's element type (K23). It was `Float64` for a `Float32` snapshot, and `compose` then
+  promoted the reduced problem to `Float64`. The snapshot of one solution, `GeometricIntegratorData`,
+  has the element type of the solution too. It was `Float64` always.
+
+- `GeometricIntegratorEnsembleData` takes an offset-indexed vector of solutions (K24), and the
+  snapshot functions read a solution that holds only its initial entry (K25). A `step` larger than
+  the number of time steps stores the initial state alone, and `snapshot_matrix` threw a
+  `BoundsError` for it.
+
+### New Features
+
+- `reduce!(model, k, CotangentLiftSVD())` stores the cotangent lift of the first `k` left singular
+  vectors of `[q p]`. The algorithm type is not exported.
+- `scripts/pod_lorenz.jl` draws the Lorenz solution, the reduced solution and its projection with
+  CairoMakie, into the directory of `--results-dir`. It has its own environment in `scripts/`. The
+  original solution uses `LorenzAttractor.plot_phase_portrait` of GeometricProblems 0.10.
+
+### Changed
+
+- The tests no longer plot. `test/integration/pod_lorenz.jl` builds its problem with
+  `LorenzAttractor.odeproblem` of GeometricProblems, in place of a local copy, and
+  `test/helpers/problems.jl` is deleted: its other two functions were empty. `Plots` leaves the
+  test environment, `OffsetArrays` joins it, and the test compat of GeometricProblems is `0.10`.
+
 ## [0.4.1] — 2026-10-29
 
 ### Changed
diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md
index 50aefd6..9f195e3 100644
--- a/KNOWN_ISSUES.md
+++ b/KNOWN_ISSUES.md
@@ -111,56 +111,3 @@
   third element is a fallback that no build reads.
 - **Kind:** defect
 - **Found:** 2026-10-28 (critic round 1b)
-
-### K22 · `reduce!` with `CotangentLiftEVD` or `ComplexSVD` throws when `k` exceeds the columns of the thin SVD
-
-- **Location:** `src/Models/ReducedBasis.jl`, `reduce!(model, k, ::CotangentLiftEVD)` and
-  `reduce!(model, k, ::ComplexSVD)`, at `svd(...).U[:, 1:k]`.
-- **Evidence:** on a `6 × 1` snapshot, `reduce!(ReducedBasisModel(ComplexSVD(), S), 2,
-  ComplexSVD())` throws `BoundsError: attempt to access 3×1 Matrix{ComplexF64} at index [1:3,
-  1:2]`, and `CotangentLiftEVD` with `k = 3` throws `BoundsError: attempt to access 3×2
-  Matrix{Float64} at index [1:3, 1:3]`. The thin SVD has `min(n, columns)` left singular vectors,
-  while a `2n × 2k` basis exists for every `k ≤ n`.
-- **Kind:** defect
-- **Found:** PR #64 (CodeRabbit review on a9d6a90)
-
-### K23 · `reduce!` with `CotangentLiftEVD` or `ComplexSVD` returns a `Float64` basis for a `Float32` snapshot
-
-- **Location:** `src/Models/ReducedBasis.jl`, `reduce!(model, k, ::CotangentLiftEVD)` and
-  `reduce!(model, k, ::ComplexSVD)`, at `zeros(2*n, 2*k)`.
-- **Evidence:** on a `6 × 4` `Float32` snapshot, `eltype(reduce!(model, 1, alg))` is `Float64` for
-  `CotangentLiftEVD()` and `ComplexSVD()`, and `Float32` for `POD()`. `compose` then promotes the
-  reduced problem to `Float64`.
-- **Kind:** defect
-- **Found:** PR #64 (CodeRabbit review on a9d6a90)
-
-### K24 · `GeometricIntegratorEnsembleData` rejects an offset-indexed vector of solutions
-
-- **Location:** `src/TrainingData/GeometricData.jl`, both `Base.Array(solutions, ::Dependent)` and
-  `Base.Array(solutions, ::Independent)`, which read `solutions[1]`, and the constructor, which
-  stores a `Vector{GeometricSolution}`.
-- **Evidence:** `GeometricIntegratorEnsembleData(OffsetVector([s, s], 0:1), deps)` throws
-  `DimensionMismatch: axes must agree, got (Base.OneTo(2),) and
-  (OffsetArrays.IdOffsetRange(values=0:1, indices=0:1),)` for `Dependent()` and for
-  `Independent()`.
-- **Kind:** defect
-- **Found:** PR #64 (CodeRabbit review on a9d6a90)
-
-### K25 · The snapshot functions read `q[1]`, which a solution with only its initial entry does not have
-
-- **Location:** `src/TrainingData/GeometricData.jl`, `snapshot_matrix`, `_snapshot_eltype` and both
-  `Base.Array` methods, at `q[1]`.
-- **Evidence:** CodeRabbit states that GeometricSolutions permits a `step` larger than the number of
-  time steps, and that `q` then holds only its entry at index `0`, so `q[1]` throws. No such
-  solution was built to check this.
-- **Kind:** not verified
-- **Found:** PR #64 (CodeRabbit review on a9d6a90)
-
-### K27 · `Parameter` with a count throws for distinct bounds one or a few ULPs apart
-
-- **Location:** `src/parameters/parameter.jl`, `Parameter(name, minimum, maximum, n::Int)`.
-- **Evidence:** `Parameter(:a, 0.1, nextfloat(0.1), 22)` throws an `AssertionError`, because
-  `LinRange(0.1, nextfloat(0.1), n)` has a point outside the bounds for 368 of the counts in
-  `2:500`. The special case for equal bounds does not cover distinct bounds.
-- **Kind:** defect
-- **Found:** PR #64 (pre-PR verification, 2026-10-29)
diff --git a/scripts/Project.toml b/scripts/Project.toml
new file mode 100644
index 0000000..67241cc
--- /dev/null
+++ b/scripts/Project.toml
@@ -0,0 +1,12 @@
+[deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+GeometricIntegrators = "dcce2d33-59f6-5b8d-9047-0defad88ae06"
+GeometricProblems = "18cb22b4-ad41-5c80-9c5f-710df63fbdc9"
+ReducedComplexityModeling = "abb4722c-1d8e-48eb-a3bd-45b3116bfaba"
+
+[compat]
+CairoMakie = "0.15"
+GeometricProblems = "0.10"
+
+[sources]
+ReducedComplexityModeling = {path = ".."}
diff --git a/scripts/pod_lorenz.jl b/scripts/pod_lorenz.jl
new file mode 100644
index 0000000..0c3a5b0
--- /dev/null
+++ b/scripts/pod_lorenz.jl
@@ -0,0 +1,50 @@
+# Proper orthogonal decomposition of the Lorenz system: the full solution, the solution of the
+# reduced problem and its projection back to the full space.
+#
+#     julia --project=scripts -e 'using Pkg; Pkg.instantiate()'
+#     julia --project=scripts scripts/pod_lorenz.jl [--results-dir DIR]
+#
+# `DIR` defaults to `results`. The script writes `pod_lorenz_original.png`,
+# `pod_lorenz_reduced.png` and `pod_lorenz_remade.png` there.
+
+using CairoMakie
+using GeometricIntegrators
+using GeometricProblems.LorenzAttractor
+using ReducedComplexityModeling
+using ReducedComplexityModeling: POD, reduce!, compose
+
+function results_dir(args)
+    i = findfirst(==("--results-dir"), args)
+    i === nothing || i < lastindex(args) ||
+        throw(ArgumentError("--results-dir needs a directory"))
+    dir = i === nothing ? "results" : args[i + 1]
+    mkpath(dir)
+end
+
+dir = results_dir(ARGS)
+
+# The same setup as test/integration/pod_lorenz.jl: a short time span, because on (0, 100) the
+# Lorenz system is chaotic and grows the round-off of the reduced integration to O(1).
+prob = LorenzAttractor.odeproblem([1.0, 0.0, 0.0]; timespan = (0.0, 1.0))
+sol = integrate(prob, ExplicitEuler())
+data = GeometricIntegratorData(sol)
+
+# two of the three modes: the reduced problem has two components, and the figure shows both
+k = 2
+alg = POD()
+model = ReducedBasisModel(alg, data.snapshot)
+reduce!(model, k, alg)
+red_sol = integrate(compose(prob, model), ExplicitEuler())
+psol = project(red_sol, model)
+
+save(joinpath(dir, "pod_lorenz_original.png"), LorenzAttractor.plot_phase_portrait(sol))
+
+fig = Figure(size = (600, 600))
+ax = Axis(fig[1, 1]; title = "Reduced basis", xlabel = "u₁", ylabel = "u₂")
+lines!(ax, [u[1] for u in red_sol.q], [u[2] for u in red_sol.q])
+save(joinpath(dir, "pod_lorenz_reduced.png"), fig)
+
+fig = Figure(size = (600, 600))
+ax = Axis3(fig[1, 1]; title = "Remade", xlabel = "x", ylabel = "y", zlabel = "z")
+lines!(ax, psol[1, :], psol[2, :], psol[3, :])
+save(joinpath(dir, "pod_lorenz_remade.png"), fig)
diff --git a/src/Models/ReducedBasis.jl b/src/Models/ReducedBasis.jl
index d9e598a..f7530ee 100644
--- a/src/Models/ReducedBasis.jl
+++ b/src/Models/ReducedBasis.jl
@@ -16,73 +16,100 @@ A reduced-basis model of the snapshot matrix `snapshot`, whose columns are the s
 the reduction scheme `alg`, such as [`POD`](@ref). [`reduce!`](@ref) computes the reduced basis
 and stores it in the field `reduced_basis`, which is undefined until then.
 """
-mutable struct ReducedBasisModel <: AbstractModel
-    reduction_scheme::Any
-    snapshot::Any
-    reduced_basis::Any
+mutable struct ReducedBasisModel{AT, T, ST <: AbstractMatrix} <: AbstractModel
+    reduction_scheme::AT
+    snapshot::ST
+    reduced_basis::Matrix{T}
 
-    function ReducedBasisModel(alg, snaps)
-        new(alg, snaps)
+    function ReducedBasisModel(alg::AT, snaps::ST) where {AT, ST <: AbstractMatrix}
+        new{AT, float(eltype(ST)), ST}(alg, snaps)
     end
 end
 
 """
     reduce!(model::ReducedBasisModel, k::Int, ::POD)
+    reduce!(model::ReducedBasisModel, k::Int, ::CotangentLiftSVD)
     reduce!(model::ReducedBasisModel, k::Int, ::CotangentLiftEVD)
     reduce!(model::ReducedBasisModel, k::Int, ::ComplexSVD)
 
 Compute a reduced basis from the snapshot matrix of `model`, whose columns are the snapshots,
 store it in `model.reduced_basis` and return it. With [`POD`](@ref), the basis is the matrix of
-the first `k` left singular vectors of the snapshot matrix. With `CotangentLiftEVD()` and
-`ComplexSVD()`, which the package does not export, the snapshot has `2n` rows, `q` above `p`, and
-the basis is an orthonormal and symplectic `2n × 2k` matrix: the cotangent lift of the first `k`
-left singular vectors of `[q p]`, or the real form of the first `k` left singular vectors of
-`q + im * p`.
+the first `k` left singular vectors of the snapshot matrix. With `CotangentLiftSVD()`,
+`CotangentLiftEVD()` and `ComplexSVD()`, which the package does not export, the snapshot has `2n`
+rows, `q` above `p`, and the basis is an orthonormal and symplectic `2n × 2k` matrix, for
+`k ≤ n`. `CotangentLiftSVD()` and `CotangentLiftEVD()` give the cotangent lift of the first `k`
+left singular vectors of `[q p]`: the first takes them from the singular value decomposition of
+`[q p]`, the second from the eigendecomposition of the Gram matrix `q qᵀ + p pᵀ`. The Gram matrix
+squares the condition number: in the floating-point type `T` of the basis, `CotangentLiftEVD()`
+resolves only the singular vectors whose singular values exceed about `sqrt(eps(T))` times the
+largest, and its eigendecomposition of an `n × n` matrix takes `O(n³)` time and `O(n²)` memory.
+`ComplexSVD()` gives the real form of the first `k` left singular vectors of `q + im * p`. The
+three algorithms take a real snapshot: a complex one throws a `MethodError`. Where `k` exceeds
+the number of singular vectors, the basis is completed with orthonormal vectors. The basis has
+the floating-point type of the snapshot's element type: `Float64` for an integer snapshot.
 """
 function reduce!(model::ReducedBasisModel, k::Int, ::POD)
-    U, _, _ = svd(model.snapshot)
-    model.reduced_basis = U[:, 1:k]
+    model.reduced_basis = _left_singular_vectors(model.snapshot, k)
 end
 
-function reduce!(model::ReducedBasisModel, k::Int, ::CotangentLiftEVD)
-    snapshot = model.snapshot
-    @assert mod(size(snapshot, 1), 2)==0 "The solutions in the snapshot are not sympletic"
-    n = floor(Int, size(snapshot, 1)/2)
-    q, p = snapshot[1:n, :], snapshot[(n + 1):end, :]
-    cotangent_snap = hcat(q, p)
-    U = svd(cotangent_snap).U[:, 1:k]
-    reduced_basis = zeros(2*n, 2*k)
-    reduced_basis[1:n, 1:k] = U
-    reduced_basis[(n + 1):(2 * n), (k + 1):(2 * k)] = U
-    model.reduced_basis = reduced_basis
+# the first `k` left singular vectors of `A`, completed to an orthonormal set where `A` has
+# fewer than `k` of them
+function _left_singular_vectors(A, k)
+    U = svd(A).U
+    r = size(U, 2)
+    k ≤ r && return U[:, 1:k]
+    k ≤ size(U, 1) || throw(BoundsError(U, (:, 1:k)))
+    # the trailing columns of the full `Q` of `U` are orthogonal to `U`; this builds `N × k`
+    # and not the `N × N` factor of a full SVD
+    B = qr(U).Q * Matrix{eltype(U)}(I, size(U, 1), k)
+    B[:, 1:r] = U
+    B
+end
+
+# `n`, `q` and `p` of a snapshot with `2n` rows
+function _split_snapshot(snapshot)
+    @assert iseven(size(snapshot, 1)) "The solutions in the snapshot are not symplectic"
+    n = size(snapshot, 1) ÷ 2
+    n, snapshot[1:n, :], snapshot[(n + 1):end, :]
+end
+
+# the `2n × 2k` matrix `blockdiag(U, U)`
+function _cotangent_lift(U)
+    n, k = size(U)
+    A = zeros(eltype(U), 2n, 2k)
+    A[1:n, 1:k] = U
+    A[(n + 1):(2n), (k + 1):(2k)] = U
+    A
+end
+
+function reduce!(model::ReducedBasisModel{AT, <:Real}, k::Int, ::CotangentLiftSVD) where {AT}
+    _, q, p = _split_snapshot(model.snapshot)
+    model.reduced_basis = _cotangent_lift(_left_singular_vectors(hcat(q, p), k))
+end
+
+function reduce!(
+        model::ReducedBasisModel{AT, T}, k::Int, ::CotangentLiftEVD) where {AT, T <: Real}
+    n, q, p = _split_snapshot(model.snapshot)
+    # scaling leaves the eigenvectors unchanged, and keeps the products finite and nonzero
+    scale = maximum(abs, model.snapshot)
+    scale = iszero(scale) ? one(scale) : scale
+    q, p = q / scale, p / scale
+    gram = Symmetric(q * q' + p * p')
+    # `eigen` sorts the eigenvalues in ascending order: take the `k` largest, largest first.
+    # It computes a `Float16` matrix in `Float32`, so the vectors are converted to `T`.
+    U = eigen(gram, (n - k + 1):n).vectors[:, k:-1:1]
+    model.reduced_basis = _cotangent_lift(convert(Matrix{T}, U))
 end
 
-# function reduce!(model::ReducedBasisModel,k::Int,::CotangentLiftEVD)
-#     snapshot = model.snapshot
-#     @assert mod(size(snapshot,1),2)==0 "The solutions in the snapshot are not sympletic"
-#     n = floor(Int,size(snapshot,1)/2)
-#     q,p = snapshot[1:n,:],snapshot[n+1:end,:]
-#     cotangent_snap = hcat(q,p)
-#     U = eigen(cotangent_snap).vectors[:,1:k]
-#     reduced_basis = zeros(2*n,2*k)
-#     reduced_basis[1:n,1:k] = U
-#     reduced_basis[n+1:2*n,k+1:2*k] = U
-#     model.reduced_basis = reduced_basis
-# end
-
-function reduce!(model::ReducedBasisModel, k::Int, ::ComplexSVD)
-    snapshot = model.snapshot
-    @assert mod(size(snapshot, 1), 2) == 0 "The solutions in the snapshot are not symplectic"
-    n = floor(Int, size(snapshot, 1)/2)
-    q, p = snapshot[1:n, :], snapshot[(n + 1):end, :]
-    csvd_snap = q + 1im * p
-    U = svd(csvd_snap).U[:, 1:k]
+function reduce!(model::ReducedBasisModel{AT, <:Real}, k::Int, ::ComplexSVD) where {AT}
+    n, q, p = _split_snapshot(model.snapshot)
+    U = _left_singular_vectors(q + im * p, k)
     Φ, Ψ = real.(U), imag.(U)
-    reduced_basis = zeros(2*n, 2*k)
+    reduced_basis = zeros(eltype(Φ), 2n, 2k)
     reduced_basis[1:n, 1:k] = Φ
-    reduced_basis[1:n, (k + 1):(2 * k)] = -Ψ
-    reduced_basis[(n + 1):(2 * n), 1:k] = Ψ
-    reduced_basis[(n + 1):(2 * n), (k + 1):(2 * k)] = Φ
+    reduced_basis[1:n, (k + 1):(2k)] = -Ψ
+    reduced_basis[(n + 1):(2n), 1:k] = Ψ
+    reduced_basis[(n + 1):(2n), (k + 1):(2k)] = Φ
     model.reduced_basis = reduced_basis
 end
 
diff --git a/src/TrainingData/GeometricData.jl b/src/TrainingData/GeometricData.jl
index bedff56..bcb3763 100644
--- a/src/TrainingData/GeometricData.jl
+++ b/src/TrainingData/GeometricData.jl
@@ -3,8 +3,8 @@
 
 function snapshot_matrix(sol::DataSeries)
     n1 = length(sol)
-    n2 = length(sol[1])
-    Z = zeros(n1, n2)
+    n2 = length(first(sol))
+    Z = zeros(eltype(first(sol)), n1, n2)
     for (i, elem) in enumerate(sol)
         Z[i, :] .= elem
     end
@@ -18,12 +18,13 @@ The solution of a geometric integrator and its snapshot matrix: the field `snaps
 with one row per component of the state `q` and one column per entry of `solution.q`, the
 initial state included.
 """
-struct GeometricIntegratorData
-    solution::GeometricSolution
-    snapshot::Matrix
+struct GeometricIntegratorData{ST <: GeometricSolution, T}
+    solution::ST
+    snapshot::Matrix{T}
 
-    function GeometricIntegratorData(solution::GeometricSolution)
-        new(solution, snapshot_matrix(solution.q))
+    function GeometricIntegratorData(solution::ST) where {ST <: GeometricSolution}
+        snapshot = snapshot_matrix(solution.q)
+        new{ST, eltype(snapshot)}(solution, snapshot)
     end
 
     # function GeometricIntegratorData(solution::ODESolution)
@@ -61,13 +62,13 @@ struct Independent <: AbstractDependence end
 
 # the element type that holds the state `q` of every solution without loss
 function _snapshot_eltype(solutions)
-    mapreduce(solution -> eltype(solution.q[1]), promote_type, solutions)
+    mapreduce(solution -> eltype(first(solution.q)), promote_type, solutions)
 end
 
 function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Dependent)
     nx = sum(solution -> length(solution.q), solutions)
-    ny = length(solutions[1].q[1])
-    all(solution -> length(solution.q[1]) == ny, solutions) ||
+    ny = length(first(first(solutions).q))
+    all(solution -> length(first(solution.q)) == ny, solutions) ||
         throw(ArgumentError("the Dependent layout needs solutions with equal state dimensions"))
     Z = zeros(_snapshot_eltype(solutions), nx, ny)
 
@@ -83,11 +84,11 @@ end
 
 function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Independent)
     ns = length(solutions)
-    nx = length(solutions[1].q)
-    ny = length(solutions[1].q[1])
+    nx = length(first(solutions).q)
+    ny = length(first(first(solutions).q))
     all(solution -> length(solution.q) == nx, solutions) ||
         throw(ArgumentError("the Independent layout needs solutions with equal numbers of entries"))
-    all(solution -> length(solution.q[1]) == ny, solutions) ||
+    all(solution -> length(first(solution.q)) == ny, solutions) ||
         throw(ArgumentError("the Independent layout needs solutions with equal state dimensions"))
     Z = zeros(_snapshot_eltype(solutions), nx, ny, ns)
 
@@ -106,13 +107,17 @@ The solutions of an ensemble of geometric integrations and their snapshot, in th
 `deps` selects: [`Dependent`](@ref) or [`Independent`](@ref). The element type of the snapshot
 is the promotion of the element types of the states `q` of all solutions.
 """
-struct GeometricIntegratorEnsembleData
-    solutions::Vector{GeometricSolution}
-    snapshot::Array
+struct GeometricIntegratorEnsembleData{ST <: GeometricSolution, T, N}
+    solutions::Vector{ST}
+    snapshot::Array{T, N}
 
     function GeometricIntegratorEnsembleData(
-            solutions::AbstractVector{<:GeometricSolution}, deps::AbstractDependence)
-        new(solutions, Array(solutions, deps))
+            solutions::AbstractVector{S},
+            deps::AbstractDependence
+    ) where {S <: GeometricSolution}
+        stored = collect(solutions)
+        snapshot = Array(solutions, deps)
+        new{S, eltype(snapshot), ndims(snapshot)}(stored, snapshot)
     end
 
     # function GeometricIntegratorEnsembleData(solutions::Vector{ODESolution},deps::AbstractDependence)
diff --git a/src/parameters/parameter.jl b/src/parameters/parameter.jl
index 39086a2..18072a7 100644
--- a/src/parameters/parameter.jl
+++ b/src/parameters/parameter.jl
@@ -18,11 +18,12 @@ samples between them. In every form, `name` is a `Symbol` or a string, which is
 
 The bounds have one type `DT`, and the samples have element type `DT`. The samples are stored
 sorted and without duplicates. Without samples, `samples` is `nothing`. With a count `n`, the
-samples are the points of `LinRange(minimum, maximum, n)`, without duplicates. Equal bounds and a
-positive count give the one sample `minimum`, because `LinRange` can round above equal bounds.
-This works for bounds of a type that `LinRange` keeps, such as a float or a `Rational`; unequal
-integer bounds with a count throw a `MethodError`, because their `LinRange` has a float element
-type. With only a vector of samples, the bounds are its smallest and largest entry.
+samples are the points of `LinRange(minimum, maximum, n)`, without duplicates. A point that
+rounds outside the bounds, by a few ULPs, is moved onto the nearest bound. Equal bounds and a
+positive count give the one sample `minimum`. This works for bounds of a type that `LinRange`
+keeps, such as a float or a `Rational`; unequal integer bounds with a count throw a `MethodError`,
+because their `LinRange` has a float element type. With only a vector of samples, the bounds are
+its smallest and largest entry.
 
 `minimum > maximum`, or a sample outside `[minimum, maximum]`, throws an `AssertionError`.
 """
@@ -48,9 +49,10 @@ function Parameter(name, minimum::DT, maximum::DT) where {DT}
     Parameter(name, minimum, maximum, nothing)
 end
 function Parameter(name, minimum::DT, maximum::DT, n::Int) where {DT}
-    # `LinRange(0.1, 0.1, 6)` holds `0.10000000000000002`, so equal bounds give their one sample
+    # `LinRange(0.1, 0.1, 6)` holds `0.10000000000000002`, so equal bounds give their one sample.
+    # Distinct bounds can also round a point a few ULPs outside: the clamp puts it back.
     samples = minimum == maximum && n ≥ 0 ? fill(minimum, min(n, 1)) :
-              LinRange(minimum, maximum, n)
+              clamp.(LinRange(minimum, maximum, n), minimum, maximum)
     Parameter(name, minimum, maximum, samples)
 end
 function Parameter(name, samples::AbstractVector)
diff --git a/test/Models/ReducedBasis.jl b/test/Models/ReducedBasis.jl
index 9c0f85a..1d75719 100644
--- a/test/Models/ReducedBasis.jl
+++ b/test/Models/ReducedBasis.jl
@@ -15,8 +15,9 @@ function lowrank_snapshot(rng, n, k, m)
     vcat(V * randn(rng, k, m), V * randn(rng, k, m))
 end
 
-@testset "reduce! with $(nameof(typeof(alg)))" for alg in (
-    RCM.CotangentLiftEVD(), RCM.ComplexSVD())
+const SYMPLECTIC_ALGS = (RCM.CotangentLiftSVD(), RCM.CotangentLiftEVD(), RCM.ComplexSVD())
+
+@testset "reduce! with $(nameof(typeof(alg)))" for alg in SYMPLECTIC_ALGS
     n, k, m = 4, 2, 10
     X = lowrank_snapshot(Xoshiro(1), n, k, m)
     model = ReducedBasisModel(alg, X)
@@ -29,3 +30,72 @@ end
     @test A' * canonical_J(n) * A ≈ canonical_J(k) atol=1e-12
     @test A * (A' * X) ≈ X rtol=1e-12
 end
+
+@testset "The cotangent lift by SVD and by EVD span the same space" begin
+    n, k, m = 6, 3, 12
+    X = lowrank_snapshot(Xoshiro(2), n, k, m)
+    A_svd = reduce!(ReducedBasisModel(RCM.CotangentLiftSVD(), X), k, RCM.CotangentLiftSVD())
+    A_evd = reduce!(ReducedBasisModel(RCM.CotangentLiftEVD(), X), k, RCM.CotangentLiftEVD())
+    # the orthogonal projectors on the two column spaces agree
+    @test A_svd * A_svd' ≈ A_evd * A_evd' atol=1e-12
+end
+
+@testset "CotangentLiftEVD with a snapshot of extreme scale, $magnitude" for magnitude in (1e160, 1e-160)
+    n, k, m = 4, 2, 10
+    X = lowrank_snapshot(Xoshiro(5), n, k, m)
+    A_ref = reduce!(ReducedBasisModel(RCM.CotangentLiftSVD(), X), k, RCM.CotangentLiftSVD())
+    alg = RCM.CotangentLiftEVD()
+    A = reduce!(ReducedBasisModel(alg, magnitude .* X), k, alg)
+    @test all(isfinite, A)
+    @test A * A' ≈ A_ref * A_ref' atol=1e-12
+end
+
+@testset "The symplectic bases reject a complex snapshot, $(nameof(typeof(alg)))" for alg in SYMPLECTIC_ALGS
+    X = ComplexF64[1 0; im 0; 0 0; 0 0]
+    @test_throws MethodError reduce!(ReducedBasisModel(alg, X), 1, alg)
+end
+
+@testset "reduce! beyond the thin SVD, $(nameof(typeof(alg)))" for alg in SYMPLECTIC_ALGS
+    # a 6 × 1 snapshot has n = 3 and one snapshot column, so k = 2 and k = 3 need completion
+    X = reshape(Float64[1, 2, 3, 4, 5, 6], 6, 1)
+    for k in 1:3
+        A = reduce!(ReducedBasisModel(alg, X), k, alg)
+        @test size(A) == (6, 2k)
+        @test A' * A ≈ I(2k) atol=1e-12
+        @test A' * canonical_J(3) * A ≈ canonical_J(k) atol=1e-12
+    end
+end
+
+@testset "POD beyond the thin SVD" begin
+    X = reshape(Float64[1, 2, 3, 4, 5, 6], 6, 1)
+    A = reduce!(ReducedBasisModel(RCM.POD(), X), 2, RCM.POD())
+    @test size(A) == (6, 2)
+    @test A' * A ≈ I(2) atol=1e-12
+end
+
+@testset "An integer snapshot gives a Float64 basis, $(nameof(typeof(alg)))" for alg in (
+    RCM.POD(), SYMPLECTIC_ALGS...)
+    X = [1 2; 3 4; 5 6; 7 8]
+    model = ReducedBasisModel(alg, X)
+    @test model isa ReducedBasisModel{typeof(alg), Float64, Matrix{Int}}
+    @test eltype(reduce!(model, 1, alg)) === Float64
+end
+
+@testset "The basis has the element type of the snapshot, $(nameof(typeof(alg)))" for alg in (
+    RCM.POD(), SYMPLECTIC_ALGS...)
+    for T in (Float16, Float32)
+        X = T.(lowrank_snapshot(Xoshiro(3), 4, 2, 10))
+        model = ReducedBasisModel(alg, X)
+        A = reduce!(model, 2, alg)
+        @test A === model.reduced_basis
+        @test eltype(A) === T
+    end
+end
+
+@testset "The fields of ReducedBasisModel are concretely typed" begin
+    X = lowrank_snapshot(Xoshiro(4), 4, 2, 10)
+    model = ReducedBasisModel(RCM.POD(), X)
+    @test model isa ReducedBasisModel{RCM.POD, Float64, Matrix{Float64}}
+    @test isconcretetype(typeof(model))
+    @test all(isconcretetype, fieldtypes(typeof(model)))
+end
diff --git a/test/Project.toml b/test/Project.toml
index 1b12581..88399dd 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -8,8 +8,8 @@ GeometricProblems = "18cb22b4-ad41-5c80-9c5f-710df63fbdc9"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReducedComplexityModeling = "abb4722c-1d8e-48eb-a3bd-45b3116bfaba"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
@@ -18,7 +18,7 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 Documenter = "1.19.0"
-GeometricProblems = "0.8, 0.9"
+GeometricProblems = "0.10"
 
 [sources]
 ReducedComplexityModeling = {path = ".."}
diff --git a/test/TrainingData/GeometricData.jl b/test/TrainingData/GeometricData.jl
index 67c081c..f9ef9fe 100644
--- a/test/TrainingData/GeometricData.jl
+++ b/test/TrainingData/GeometricData.jl
@@ -1,4 +1,5 @@
 using GeometricIntegrators
+using OffsetArrays: OffsetVector
 using ReducedComplexityModeling
 using Test
 
@@ -53,6 +54,12 @@ end
     @test Z == hcat(RCM.snapshot_matrix(sol32.q), RCM.snapshot_matrix(sol64.q))
 end
 
+@testset "The snapshot of one solution has its element type" begin
+    sol32 = integrate(ODEProblem(oscillator!, (0.0f0, 0.1f0), 0.01f0, Float32[1, 0]),
+        ExplicitEuler())
+    @test eltype(GeometricIntegratorData(sol32).snapshot) === Float32
+end
+
 @testset "Dependent: the snapshot matrices side by side" begin
     for sols in ((sol1, sol2), (sol1, sol3), (sol3, sol1))
         Z = GeometricIntegratorEnsembleData(RCM.GeometricSolution[sols...], Dependent()).snapshot
@@ -70,6 +77,44 @@ end
         RCM.GeometricSolution[sol1d, sol1], Dependent())
 end
 
+@testset "An offset-indexed vector of solutions" begin
+    sols = OffsetVector([sol1, sol2], 0:1)
+    for deps in (Dependent(), Independent())
+        data = GeometricIntegratorEnsembleData(sols, deps)
+        @test data.snapshot ==
+              GeometricIntegratorEnsembleData(RCM.GeometricSolution[sol1, sol2], deps).snapshot
+        @test data.solutions == [sol1, sol2]
+    end
+end
+
+@testset "A solution with only its initial entry" begin
+    # a step larger than the number of time steps stores the initial state alone
+    prob = ODEProblem(oscillator!, (0.0, 0.1), 0.01, [1.0, 0.0])
+    sol0 = RCM.GeometricSolution(prob, 1000)
+    @test length(sol0.q) == 1
+
+    @test RCM.snapshot_matrix(sol0.q) == reshape([1.0, 0.0], 2, 1)
+    @test GeometricIntegratorData(sol0).snapshot == reshape([1.0, 0.0], 2, 1)
+    sols = RCM.GeometricSolution[sol0, sol0]
+    @test GeometricIntegratorEnsembleData(sols, Dependent()).snapshot == [1.0 1.0; 0.0 0.0]
+    @test GeometricIntegratorEnsembleData(sols, Independent()).snapshot ==
+          reshape([1.0, 0.0, 1.0, 0.0], 2, 1, 2)
+end
+
+@testset "The fields of the data types are concretely typed" begin
+    data = GeometricIntegratorData(sol1)
+    @test isconcretetype(typeof(data))
+    @test all(isconcretetype, fieldtypes(typeof(data)))
+    @test typeof(data.snapshot) === Matrix{Float64}
+
+    for (deps, N) in ((Dependent(), 2), (Independent(), 3))
+        ens = GeometricIntegratorEnsembleData([sol1, sol2], deps)
+        @test isconcretetype(typeof(ens))
+        @test all(isconcretetype, fieldtypes(typeof(ens)))
+        @test typeof(ens.snapshot) === Array{Float64, N}
+    end
+end
+
 @testset "A vector with a concrete element type" begin
     sols = [sol1, sol2]
     @test !(sols isa Vector{RCM.GeometricSolution})
diff --git a/test/helpers/problems.jl b/test/helpers/problems.jl
deleted file mode 100644
index 784f4a4..0000000
--- a/test/helpers/problems.jl
+++ /dev/null
@@ -1,21 +0,0 @@
-using GeometricIntegrators
-
-function lorenz_problem(; tspan = (0.0, 100.0))
-    function lorenz(du, t, u, p)
-        du[1] = 10.0(u[2]-u[1])
-        du[2] = u[1]*(28.0-u[3]) - u[2]
-        du[3] = u[1]*u[2] - (8/3)*u[3]
-    end
-    u0 = [1.0, 0.0, 0.0]
-    tstep = 0.01
-    prob = ODEProblem(lorenz, tspan, tstep, u0)
-    sol = integrate(prob, ExplicitEuler())
-    prob, sol
-end
-
-# Implicitly assume that all PDEs are solved in 2 dimensions.
-function heat_equation()
-end
-
-function elastic_wave_equations()
-end
diff --git a/test/integration/pod_lorenz.jl b/test/integration/pod_lorenz.jl
index 1aef9ed..86eff4f 100644
--- a/test/integration/pod_lorenz.jl
+++ b/test/integration/pod_lorenz.jl
@@ -1,4 +1,5 @@
-include("../helpers/problems.jl")
+using GeometricIntegrators
+using GeometricProblems.LorenzAttractor
 using LinearAlgebra
 using ReducedComplexityModeling
 using ReducedComplexityModeling: POD, reduce!, compose
@@ -12,7 +13,8 @@ begin
 
     # 1. Get data. A short time span: on (0, 100) the Lorenz system is chaotic and grows the
     # round-off of the reduced integration to O(1).
-    prob, sol = lorenz_problem(; tspan = (0.0, 1.0))
+    prob = LorenzAttractor.odeproblem([1.0, 0.0, 0.0]; timespan = (0.0, 1.0))
+    sol = integrate(prob, ExplicitEuler())
     data = GeometricIntegratorData(sol)
 
     # 2. Define a model with an appropriate algorithm.
@@ -36,15 +38,4 @@ begin
         @test size(psol) == size(data.snapshot)
         @test psol ≈ data.snapshot rtol=1e-12
     end
-
-    #------------------------------------------------------------------------------
-
-    # Visualization and Shennanigans
-    using Plots
-    atemp = 1:k
-    display(plot(red_sol.q, vars = (1, 2), title = "Reduced Basis"))
-    display(plot(psol[1, :], psol[2, :], psol[3, :], title = "Remade"))
-    display(plot(sol.q, vars = (1, 2, 3), title = "Original"))
 end
-
-# Heat Equation. Solver first. Then use the Reduced basis.
diff --git a/test/parameters/parameter.jl b/test/parameters/parameter.jl
index f71074b..17943aa 100644
--- a/test/parameters/parameter.jl
+++ b/test/parameters/parameter.jl
@@ -143,6 +143,19 @@ end
     end
 end
 
+@testset "Bounds a few ULPs apart give samples inside the bounds" begin
+    for T in (Float32, Float64), n in 2:200
+
+        for (a, b) in ((T(0.1), nextfloat(T(0.1))), (T(0.1), nextfloat(T(0.1), 3)))
+            p = Parameter(:a, a, b, n)
+            @test all(s -> a ≤ s ≤ b, collect(p))
+            @test issorted(collect(p))
+            @test first(collect(p)) == a
+            @test last(collect(p)) == b
+        end
+    end
+end
+
 @testset "Parameter show" begin
     p = Parameter(:μ, 0.0, 1.0, 3)
     @test sprint(show, MIME"text/plain"(), p) ==