diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0f1d5b6..6506f3b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,69 @@ tags. It is named as a gap rather than reconstructed, because a changelog assemb
 fact loses exactly the reasoning that makes it worth keeping. The `[Unreleased]` section below
 holds the changes since the last release until that release is cut.
 
+## [Unreleased]
+
+### Breaking Changes
+
+- `ReducedComplexityModeling` no longer exports `Array`. The binding is `Core`'s, and every module
+  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.
+
+### Bug Fixes
+
+- `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
+  positive count give the one sample `minimum`. A count of 0 gives no sample. It threw a
+  `MethodError` before, because two methods matched it equally well. This relates to issue #38:
+  the two string-name methods that clashed are deleted, and
+  `detect_ambiguities(ReducedComplexityModeling)` reports no ambiguity. A gap remains: unequal
+  integer bounds with a sample count still throw a `MethodError`. Their `LinRange` has a float
+  element type, so `Parameter(:a, 0, 10, 5)` throws.
+
+- `GeometricIntegratorEnsembleData(solutions, Independent())` builds its snapshot (K16). It threw
+  an `ErrorException` before. The snapshot is an array of size `ny × nx × ns`, for `ns` solutions
+  with `nx` entries of a state `q` with `ny` components. Its slice `[:, :, i]` is the snapshot
+  matrix of solution `i`, in the layout of `GeometricIntegratorData`. Solutions with different
+  numbers of entries in `q` throw an `ArgumentError`. Solutions with different numbers of
+  components of `q` also throw an `ArgumentError`.
+
+- `GeometricIntegratorEnsembleData(solutions, Dependent())` sizes its snapshot from the total number
+  of entries of all solutions. It was sized from the first solution alone: a later solution with
+  fewer entries left zero-filled trailing columns, and one with more entries threw a `BoundsError`.
+  The solutions may now have different numbers of entries in `q`. Solutions with different numbers
+  of components of `q` throw an `ArgumentError` (K26). A solution with one component after one with
+  two was copied into both rows without an error before. The snapshot of one solution and the
+  `Dependent` snapshot change their layout with `permutedims`. `Array(Z')` conjugated complex
+  values. No complex state reaches them today, because `integrate` does not take one.
+
+- In both layouts, the snapshot of `GeometricIntegratorEnsembleData` has the `promote_type` of
+  the element types of the states `q` of all solutions. A `Float32` solution followed by a
+  `Float64` solution gives a `Float64` snapshot, so a value of `1e40`, which is above
+  `floatmax(Float32)`, is kept. The order of the columns is unchanged.
+
+- `GeometricIntegratorEnsembleData` and both `Base.Array(solutions, ::Dependent)` and
+  `Base.Array(solutions, ::Independent)` take a vector with a concrete element type, such as
+  `[sol, sol]`. `GeometricIntegratorEnsembleData([sol, sol], Dependent())` threw a `MethodError`
+  before.
+
+- `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`.
+
+- Help text: the module and every exported name have a docstring with text. `Parameter`,
+  `ParameterSampler` and `CartesianParameterSampler` replace their empty docstrings, and `POD`,
+  `reduce!`, `compose`, `project`, `hassamples`, `GeometricIntegratorData`,
+  `GeometricIntegratorEnsembleData`, `Dependent` and `Independent` gain one. The exports test
+  requires the docstring of every exported or public name to contain that name. The docstrings of
+  `h5save` and `h5load` name their three methods and their arguments (K21). Each shows a round
+  trip, as a doctest. `h5save` writes a `ParameterSpace`, and `h5load` reads it back.
+  `h5save(file, p)` writes a `Parameter` `p` to the group of its name.
+  `h5load(Parameter, file; path = string(p.name))` reads it back. With the default
+  `path = "/"`, `h5load(Parameter, file)` throws a `KeyError`. The docstring of
+  `ReducedBasisModel` no longer names packages that are not dependencies of this package.
+
 ## [0.4.1] — 2026-10-29
 
 ### Changed
diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md
index 16baac8..50aefd6 100644
--- a/KNOWN_ISSUES.md
+++ b/KNOWN_ISSUES.md
@@ -67,7 +67,7 @@
 
 ### K13 · Explicit import checks report implicit and non-public dependency accesses
 
-- **Location:** `src/ReducedComplexityModeling.jl:4`, the dependency imports, and their uses in
+- **Location:** `src/ReducedComplexityModeling.jl:13`, the dependency imports, and their uses in
   included source files.
 - **Evidence:** on Julia 1.13.3, ExplicitImports 1.15.0,
   `print_explicit_imports(ReducedComplexityModeling; report_non_public = true)` reports 25
@@ -75,7 +75,7 @@
   `ArrayOrNamedTuple`, `H5DataStore`, `NamedTupleOfArrays`), and 8 non-public qualified accesses.
   `test_explicit_imports(ReducedComplexityModeling)` fails the implicit-import, public-import and
   public-qualified-access checks. `ArrayOrNamedTuple` and `NamedTupleOfArrays` are imported at
-  `src/ReducedComplexityModeling.jl:26`, for the data loader's dispatch; the other import
+  `src/ReducedComplexityModeling.jl:35`, for the data loader's dispatch; the other import
   statements are present at merge base `3b0554e`. The `Chain` and `NeuralNetwork` imports have no
   diagnostic.
 - **Kind:** defect
@@ -111,3 +111,56 @@
   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/docs/make.jl b/docs/make.jl
index 7bb7e24..02a13c2 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -33,6 +33,7 @@ makedocs(;
     authors = "Michael Kraus",
     repo = "https://github.com/JuliaRCM/ReducedComplexityModeling.jl/blob/{commit}{path}#{line}",
     sitename = "ReducedComplexityModeling.jl",
+    checkdocs = :all,
     format = Documenter.HTML(;
         prettyurls = get(ENV, "CI", "false") == "true",
         canonical = "https://JuliaRCM.github.io/ReducedComplexityModeling.jl",
diff --git a/src/Models/ReducedBasis.jl b/src/Models/ReducedBasis.jl
index bb91f66..d9e598a 100644
--- a/src/Models/ReducedBasis.jl
+++ b/src/Models/ReducedBasis.jl
@@ -1,14 +1,20 @@
+"""
+    POD()
+
+Proper orthogonal decomposition: `reduce!(model, k, POD())` takes the first `k` left singular
+vectors of the snapshot matrix as the reduced basis.
+"""
 struct POD end
 struct CotangentLiftSVD end
 struct CotangentLiftEVD end
 struct ComplexSVD end
 
 """
-Reduced Basis
-
-Uses ModelOrderReduction methods from ReducedBasisMethods.jl
-to represent the data in latent space.
+    ReducedBasisModel(alg, snapshot)
 
+A reduced-basis model of the snapshot matrix `snapshot`, whose columns are the snapshots, with
+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
@@ -20,6 +26,19 @@ mutable struct ReducedBasisModel <: AbstractModel
     end
 end
 
+"""
+    reduce!(model::ReducedBasisModel, k::Int, ::POD)
+    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`.
+"""
 function reduce!(model::ReducedBasisModel, k::Int, ::POD)
     U, _, _ = svd(model.snapshot)
     model.reduced_basis = U[:, 1:k]
@@ -31,7 +50,7 @@ function reduce!(model::ReducedBasisModel, k::Int, ::CotangentLiftEVD)
     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)[:, 1:k]
+    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
@@ -57,16 +76,24 @@ function reduce!(model::ReducedBasisModel, k::Int, ::ComplexSVD)
     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)
-    ϕ, Ψ = real.(U), imag.(U)
+    U = svd(csvd_snap).U[:, 1:k]
+    Φ, Ψ = real.(U), imag.(U)
     reduced_basis = zeros(2*n, 2*k)
     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[(n + 1):(2 * n), (k + 1):(2 * k)] = Φ
     model.reduced_basis = reduced_basis
 end
 
+"""
+    compose(prob::ODEProblem, model::ReducedBasisModel)
+
+The reduced problem of `prob` in the reduced basis `R = model.reduced_basis`: an `ODEProblem` in
+the reduced coordinates `u`, with the vector field `R' * v(R * u)`, where `v` is the vector field
+of `prob`, and the initial condition `R' * q₀`. The time span, the time step and the parameters
+are those of `prob`. Call [`reduce!`](@ref) on `model` first.
+"""
 function compose(prob::ODEProblem, model::ReducedBasisModel)
     Rᵦ = model.reduced_basis
 
@@ -84,6 +111,13 @@ function compose(prob::ODEProblem, model::ReducedBasisModel)
     red_prob
 end
 
+"""
+    project(rsol::GeometricSolution, model::ReducedBasisModel)
+
+Map the solution `rsol` of a reduced problem (see [`compose`](@ref)) back to the full space. The
+result is the matrix `model.reduced_basis * Q`, where the columns of `Q` are the entries of
+`rsol.q`; its columns are the states in the full space.
+"""
 function project(rsol::GeometricSolution, model::ReducedBasisModel)
     Rᵦ = model.reduced_basis
     Rᵦ * snapshot_matrix(rsol.q)
diff --git a/src/ReducedComplexityModeling.jl b/src/ReducedComplexityModeling.jl
index 2452517..c844dcd 100644
--- a/src/ReducedComplexityModeling.jl
+++ b/src/ReducedComplexityModeling.jl
@@ -1,3 +1,12 @@
+"""
+    ReducedComplexityModeling
+
+Reduced models of the solutions of geometric integrators: the snapshot data of a solution or an
+ensemble of solutions, reduced bases ([`ReducedBasisModel`](@ref), [`POD`](@ref)), parameters
+and parameter spaces with their samplers, the fit of a model's parameters to data
+([`TrainingProblem`](@ref), [`train`](@ref)), and the loaders that draw batches of data for a
+neural network ([`DataLoader`](@ref), [`Batch`](@ref)).
+"""
 module ReducedComplexityModeling
 
 using ForwardDiff
@@ -62,7 +71,6 @@ export ObservableSpace, IntrinsicSpace, TrainingData
 
 include("Models/Models.jl")
 
-export Array
 export GeometricIntegratorData, GeometricIntegratorEnsembleData
 export Dependent, Independent
 export ReducedBasisModel
diff --git a/src/TrainingData/GeometricData.jl b/src/TrainingData/GeometricData.jl
index 10651de..bedff56 100644
--- a/src/TrainingData/GeometricData.jl
+++ b/src/TrainingData/GeometricData.jl
@@ -8,9 +8,16 @@ function snapshot_matrix(sol::DataSeries)
     for (i, elem) in enumerate(sol)
         Z[i, :] .= elem
     end
-    Array(Z')
+    permutedims(Z)
 end
 
+"""
+    GeometricIntegratorData(solution::GeometricSolution)
+
+The solution of a geometric integrator and its snapshot matrix: the field `snapshot` is a matrix
+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
@@ -29,15 +36,40 @@ end
 
 abstract type AbstractDependence end
 
+"""
+    Dependent()
+
+The layout of the snapshot of a [`GeometricIntegratorEnsembleData`](@ref) in which the solutions
+form one snapshot matrix: one row per component of the state `q`, and the entries of `q` of the
+first solution, then those of the second, and so on, as its columns. The solutions may have
+different numbers of entries in `q`, but all must have the same number of components of `q`;
+otherwise the constructor throws an `ArgumentError`.
+"""
 struct Dependent <: AbstractDependence end
 
+"""
+    Independent()
+
+The layout of the snapshot of a [`GeometricIntegratorEnsembleData`](@ref) in which each solution
+keeps its own snapshot, along a third dimension: an array of size `ny × nx × ns`, for `ns`
+solutions with `nx` entries of a state `q` with `ny` components. Its slice `[:, :, i]` is the
+snapshot matrix of solution `i`, in the layout of [`GeometricIntegratorData`](@ref). All
+solutions must have the same number of entries in `q` and the same number of components of `q`;
+otherwise the constructor throws an `ArgumentError`.
+"""
 struct Independent <: AbstractDependence end
 
-function Base.Array(solutions::Vector{GeometricSolution}, ::Dependent)
-    ns = length(solutions)
-    nx = length(solutions[1].q)
+# 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)
+end
+
+function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Dependent)
+    nx = sum(solution -> length(solution.q), solutions)
     ny = length(solutions[1].q[1])
-    Z = zeros(eltype(solutions[1].q[1]), ns*nx, ny)
+    all(solution -> length(solution.q[1]) == ny, solutions) ||
+        throw(ArgumentError("the Dependent layout needs solutions with equal state dimensions"))
+    Z = zeros(_snapshot_eltype(solutions), nx, ny)
 
     i = 1
     for solution in solutions
@@ -46,28 +78,40 @@ function Base.Array(solutions::Vector{GeometricSolution}, ::Dependent)
             i += 1
         end
     end
-    Array(Z')
+    permutedims(Z)
 end
 
-function Base.Array(solutions::Vector{GeometricSolution}, ::Independent)
+function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Independent)
     ns = length(solutions)
     nx = length(solutions[1].q)
     ny = length(solutions[1].q[1])
-    Z = zeros(eltype(solutions[1].q[1]), nx, ny, ns)
+    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) ||
+        throw(ArgumentError("the Independent layout needs solutions with equal state dimensions"))
+    Z = zeros(_snapshot_eltype(solutions), nx, ny, ns)
 
     for (i, solution) in enumerate(solutions)
         for (j, elem) in enumerate(solution.q)
             Z[j, :, i] .= elem
         end
     end
-    Array(Z')
+    permutedims(Z, (2, 1, 3))
 end
 
+"""
+    GeometricIntegratorEnsembleData(solutions::AbstractVector{<:GeometricSolution}, deps)
+
+The solutions of an ensemble of geometric integrations and their snapshot, in the layout that
+`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
 
-    function GeometricIntegratorEnsembleData(solutions::Vector{GeometricSolution}, deps::AbstractDependence)
+    function GeometricIntegratorEnsembleData(
+            solutions::AbstractVector{<:GeometricSolution}, deps::AbstractDependence)
         new(solutions, Array(solutions, deps))
     end
 
diff --git a/src/parameters/parameter.jl b/src/parameters/parameter.jl
index 648cc59..39086a2 100644
--- a/src/parameters/parameter.jl
+++ b/src/parameters/parameter.jl
@@ -7,7 +7,24 @@ _sort(s::Nothing) = s
 _sort(s::AbstractVector) = sort(unique(collect(s)))
 
 """
-
+    Parameter(name, minimum, maximum)
+    Parameter(name, minimum, maximum, samples)
+    Parameter(name, minimum, maximum, n::Int)
+    Parameter(name, samples::AbstractVector)
+
+A named parameter of a system: the bounds `minimum ≤ maximum` and, optionally, a vector of
+samples between them. In every form, `name` is a `Symbol` or a string, which is converted to a
+`Symbol`.
+
+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.
+
+`minimum > maximum`, or a sample outside `[minimum, maximum]`, throws an `AssertionError`.
 """
 struct Parameter{DT <: Number, ST <: AbstractSample{DT}}
     name::Symbol
@@ -15,26 +32,26 @@ struct Parameter{DT <: Number, ST <: AbstractSample{DT}}
     maximum::DT
     samples::ST
 
-    function Parameter(name::Symbol, minimum::DT, maximum::DT, samples::ST) where {
-            DT, ST <: AbstractSample{DT}}
+    function Parameter(name::Union{Symbol, AbstractString}, minimum::DT, maximum::DT,
+            samples::ST) where {DT, ST <: AbstractSample{DT}}
         @assert minimum ≤ maximum
         if typeof(samples) <: AbstractVector
             @assert all([minimum ≤ s for s in samples])
             @assert all([maximum ≥ s for s in samples])
         end
         _samples = _sort(samples)
-        new{DT, typeof(_samples)}(name, minimum, maximum, _samples)
+        new{DT, typeof(_samples)}(Symbol(name), minimum, maximum, _samples)
     end
 end
 
-function Parameter(name::AbstractString, minimum, maximum, samples)
-    Parameter(Symbol(name), minimum, maximum, samples)
-end
 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}
-    Parameter(name, minimum, maximum, LinRange(minimum, maximum, n))
+    # `LinRange(0.1, 0.1, 6)` holds `0.10000000000000002`, so equal bounds give their one sample
+    samples = minimum == maximum && n ≥ 0 ? fill(minimum, min(n, 1)) :
+              LinRange(minimum, maximum, n)
+    Parameter(name, minimum, maximum, samples)
 end
 function Parameter(name, samples::AbstractVector)
     Parameter(name, minimum(samples), maximum(samples), samples)
@@ -73,6 +90,11 @@ Base.length(p::ParameterWithoutSamples) = 0
 Base.size(p::ParameterWithSamples) = size(p.samples)
 Base.size(p::ParameterWithoutSamples) = (0,)
 
+"""
+    hassamples(p::Parameter)
+
+Whether the parameter `p` has a non-empty vector of samples.
+"""
 hassamples(p::ParameterWithSamples) = length(p.samples) > 0
 hassamples(p::ParameterWithoutSamples) = false
 
diff --git a/src/parameters/parametersampler.jl b/src/parameters/parametersampler.jl
index e04593b..437837c 100644
--- a/src/parameters/parametersampler.jl
+++ b/src/parameters/parametersampler.jl
@@ -1,5 +1,9 @@
 """
+    ParameterSampler
 
+The supertype of the samplers that [`sample`](@ref) draws the points of a parameter space with:
+[`CartesianParameterSampler`](@ref), [`RandomParameterSampler`](@ref) and
+[`QuasiRandomParameterSampler`](@ref).
 """
 abstract type ParameterSampler end
 
@@ -25,7 +29,12 @@ function sample(ps::ParameterSampler, parameters::Vararg{Parameter})
 end
 
 """
+    CartesianParameterSampler()
 
+Draws every combination of the samples stored in the parameters: the Cartesian product of their
+sample vectors, with the first parameter varying fastest. Each parameter must have samples
+(see [`hassamples`](@ref)); a parameter without them throws an `AssertionError`. The samples keep
+the element type of each parameter.
 """
 struct CartesianParameterSampler <: ParameterSampler end
 
diff --git a/src/parameters/parameterspace.jl b/src/parameters/parameterspace.jl
index 28c03ac..8e980ac 100644
--- a/src/parameters/parameterspace.jl
+++ b/src/parameters/parameterspace.jl
@@ -68,7 +68,26 @@ function ParameterSpace(fpath::AbstractString, path::AbstractString = "/")
 end
 
 """
-save ParameterSpace
+    h5save(h5::H5DataStore, ps::ParameterSpace; path = "/")
+    h5save(h5::H5DataStore, param::Parameter; path = string(param.name))
+    h5save(fpath::AbstractString, data, args...; mode = "w", kwargs...)
+
+Write `ps` or `param` to the group `path` of `h5`. A [`ParameterSpace`](@ref) writes one dataset
+per parameter to the subgroup `samples`, and each of its parameters to the subgroup `parameters`.
+A [`Parameter`](@ref) writes the datasets `minimum`, `maximum` and, if it has samples, `samples`.
+With a file path, open the file `fpath` in `mode` and write `data` to it. [`h5load`](@ref) reads
+the data back.
+
+```jldoctest
+julia> ps = ParameterSpace(Parameter(:μ, 0.0, 1.0, 3), Parameter(:ν, 1.0, 2.0, 2));
+
+julia> file = joinpath(mktempdir(), "parameterspace.h5");
+
+julia> h5save(file, ps);
+
+julia> h5load(ParameterSpace, file) == ps
+true
+```
 """
 function h5save(h5::H5DataStore, ps::ParameterSpace; path::AbstractString = "/")
     cols = columns(ps.samples)
@@ -85,7 +104,27 @@ function h5save(h5::H5DataStore, ps::ParameterSpace; path::AbstractString = "/")
 end
 
 """
-Load ParameterSpace
+    h5load(ParameterSpace, h5::H5DataStore; path = "/")
+    h5load(Parameter, h5::H5DataStore; path = "/")
+    h5load(T::Type, fpath::AbstractString, args...; mode = "r", kwargs...)
+
+Read a [`ParameterSpace`](@ref) or a [`Parameter`](@ref) from the group `path` of `h5`, in the
+layout that [`h5save`](@ref) writes. A parameter takes the name of its group. With a file path,
+open the file `fpath` in `mode` and read a `T` from it.
+
+`h5save` writes a parameter `p` to the group of its name, so read it back with
+`path = string(p.name)`:
+
+```jldoctest
+julia> p = Parameter(:μ, 0.0, 1.0, 3);
+
+julia> file = joinpath(mktempdir(), "parameter.h5");
+
+julia> h5save(file, p);
+
+julia> h5load(Parameter, file; path = string(p.name)) == p
+true
+```
 """
 function h5load(::Type{ParameterSpace}, h5::H5DataStore; path::AbstractString = "/")
     ParameterSpace(h5, path)
diff --git a/test/Models/ReducedBasis.jl b/test/Models/ReducedBasis.jl
new file mode 100644
index 0000000..9c0f85a
--- /dev/null
+++ b/test/Models/ReducedBasis.jl
@@ -0,0 +1,31 @@
+using LinearAlgebra
+using Random: Xoshiro
+using ReducedComplexityModeling
+using Test
+
+const RCM = ReducedComplexityModeling
+
+# The canonical symplectic matrix of dimension 2m.
+canonical_J(m) = [zeros(m, m) I(m); -I(m) zeros(m, m)]
+
+# A snapshot whose q and p columns lie in the span of the orthonormal columns of `V`, so a basis
+# of `size(V, 2)` vectors in q and in p reconstructs it.
+function lowrank_snapshot(rng, n, k, m)
+    V = Matrix(qr(randn(rng, n, k)).Q)
+    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())
+    n, k, m = 4, 2, 10
+    X = lowrank_snapshot(Xoshiro(1), n, k, m)
+    model = ReducedBasisModel(alg, X)
+    A = reduce!(model, k, alg)
+
+    @test A === model.reduced_basis
+    @test size(A) == (2n, 2k)
+    # the residuals are below 1e-15; a wrong basis gives O(1)
+    @test A' * A ≈ I(2k) atol=1e-12
+    @test A' * canonical_J(n) * A ≈ canonical_J(k) atol=1e-12
+    @test A * (A' * X) ≈ X rtol=1e-12
+end
diff --git a/test/Project.toml b/test/Project.toml
index c799742..1b12581 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -7,6 +7,7 @@ GeometricIntegrators = "dcce2d33-59f6-5b8d-9047-0defad88ae06"
 GeometricProblems = "18cb22b4-ad41-5c80-9c5f-710df63fbdc9"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/test/TrainingData/GeometricData.jl b/test/TrainingData/GeometricData.jl
new file mode 100644
index 0000000..67c081c
--- /dev/null
+++ b/test/TrainingData/GeometricData.jl
@@ -0,0 +1,82 @@
+using GeometricIntegrators
+using ReducedComplexityModeling
+using Test
+
+const RCM = ReducedComplexityModeling
+
+# The harmonic oscillator q̇₁ = q₂, q̇₂ = −q₁; each initial condition gives a different solution.
+oscillator!(v, t, q, params) = (v[1] = q[2]; v[2] = -q[1])
+
+function oscillator_solution(q0, tend)
+    integrate(ODEProblem(oscillator!, (0.0, tend), 0.01, q0), ExplicitEuler())
+end
+
+sol1 = oscillator_solution([1.0, 0.0], 0.1)
+sol2 = oscillator_solution([0.0, 2.0], 0.1)
+sol3 = oscillator_solution([3.0, 1.0], 0.05)
+
+@testset "Independent: slice i is the snapshot matrix of solution i" begin
+    Z = GeometricIntegratorEnsembleData(RCM.GeometricSolution[sol1, sol2], Independent()).snapshot
+    @test size(Z) == (2, length(sol1.q), 2)
+    @test Z[:, :, 1] == RCM.snapshot_matrix(sol1.q)
+    @test Z[:, :, 2] == RCM.snapshot_matrix(sol2.q)
+    @test Z[:, 1, 2] == [0.0, 2.0]
+
+    @test_throws ArgumentError GeometricIntegratorEnsembleData(
+        RCM.GeometricSolution[sol1, sol3], Independent())
+
+    # a solution with one component after one with two
+    decay!(v, t, q, params) = (v[1] = -q[1])
+    sol1d = integrate(ODEProblem(decay!, (0.0, 0.1), 0.01, [1.0]), ExplicitEuler())
+    @test length(sol1d.q) == length(sol1.q)
+    @test_throws ArgumentError GeometricIntegratorEnsembleData(
+        RCM.GeometricSolution[sol1, sol1d], Independent())
+end
+
+@testset "The snapshot promotes the element types of all solutions" begin
+    # 1e40 overflows to Inf in Float32
+    sol32 = integrate(ODEProblem(oscillator!, (0.0f0, 0.1f0), 0.01f0, Float32[1, 0]),
+        ExplicitEuler())
+    sol64 = integrate(ODEProblem(oscillator!, (0.0, 0.1), 0.01, [1e40, 0.0]), ExplicitEuler())
+    @test eltype(sol32.q[1]) === Float32
+    sols = RCM.GeometricSolution[sol32, sol64]
+
+    Z = GeometricIntegratorEnsembleData(sols, Independent()).snapshot
+    @test eltype(Z) === Float64
+    @test Z[1, 1, 2] == 1e40
+    @test Z[:, :, 1] == RCM.snapshot_matrix(sol32.q)
+    @test Z[:, :, 2] == RCM.snapshot_matrix(sol64.q)
+
+    Z = GeometricIntegratorEnsembleData(sols, Dependent()).snapshot
+    @test eltype(Z) === Float64
+    @test Z[1, length(sol32.q) + 1] == 1e40
+    @test Z == hcat(RCM.snapshot_matrix(sol32.q), RCM.snapshot_matrix(sol64.q))
+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
+        @test Z == hcat(RCM.snapshot_matrix(sols[1].q), RCM.snapshot_matrix(sols[2].q))
+    end
+end
+
+@testset "Dependent: solutions with different state dimensions throw" begin
+    decay!(v, t, q, params) = (v[1] = -q[1])
+    sol1d = integrate(ODEProblem(decay!, (0.0, 0.1), 0.01, [1.0]), ExplicitEuler())
+    @test length(sol1d.q) == length(sol1.q)
+    @test_throws ArgumentError GeometricIntegratorEnsembleData(
+        RCM.GeometricSolution[sol1, sol1d], Dependent())
+    @test_throws ArgumentError GeometricIntegratorEnsembleData(
+        RCM.GeometricSolution[sol1d, sol1], Dependent())
+end
+
+@testset "A vector with a concrete element type" begin
+    sols = [sol1, sol2]
+    @test !(sols isa Vector{RCM.GeometricSolution})
+    for deps in (Dependent(), Independent())
+        data = GeometricIntegratorEnsembleData(sols, deps)
+        @test data.snapshot ==
+              GeometricIntegratorEnsembleData(RCM.GeometricSolution[sol1, sol2], deps).snapshot
+        @test data.solutions == sols
+    end
+end
diff --git a/test/helpers/problems.jl b/test/helpers/problems.jl
index 00bb227..784f4a4 100644
--- a/test/helpers/problems.jl
+++ b/test/helpers/problems.jl
@@ -1,6 +1,6 @@
 using GeometricIntegrators
 
-function lorenz_problem()
+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]
@@ -8,7 +8,6 @@ function lorenz_problem()
     end
     u0 = [1.0, 0.0, 0.0]
     tstep = 0.01
-    tspan = (0.0, 100.0)
     prob = ODEProblem(lorenz, tspan, tstep, u0)
     sol = integrate(prob, ExplicitEuler())
     prob, sol
diff --git a/test/integration/exports.jl b/test/integration/exports.jl
index 2ffbc3b..7a31f6e 100644
--- a/test/integration/exports.jl
+++ b/test/integration/exports.jl
@@ -1,7 +1,23 @@
 using ReducedComplexityModeling
 using Test
 
-@testset "Every exported name is defined" begin
-    @test isempty(filter(
-        n -> !isdefined(ReducedComplexityModeling, n), names(ReducedComplexityModeling)))
+# The raw text of the docstrings of the binding `name` of `mod`, stripped. `Docs.hasdoc` is also
+# true for an empty docstring, so the text is what shows that a docstring says something: it names
+# the binding, as a signature line does.
+function docstring_text(mod::Module, name::Symbol)
+    binding = Docs.Binding(mod, name)
+    multidoc = get(Docs.meta(binding.mod), binding, nothing)
+    isnothing(multidoc) && return ""
+    return strip(join((join(string.(doc.text)) for doc in values(multidoc.docs)), "\n"))
+end
+
+# `names` returns the exported and the public names, and the name of the module itself. An
+# undefined name fails the loop, because `getfield` throws for it.
+@testset "Every exported or public name is documented and owned by the package" begin
+    @testset "$name" for name in names(ReducedComplexityModeling)
+        @test Docs.hasdoc(ReducedComplexityModeling, name)
+        @test occursin(string(name), docstring_text(ReducedComplexityModeling, name))
+        @test parentmodule(getfield(ReducedComplexityModeling, name)) ===
+              ReducedComplexityModeling
+    end
 end
diff --git a/test/integration/pod_lorenz.jl b/test/integration/pod_lorenz.jl
index ebe8775..1aef9ed 100644
--- a/test/integration/pod_lorenz.jl
+++ b/test/integration/pod_lorenz.jl
@@ -1,6 +1,8 @@
 include("../helpers/problems.jl")
+using LinearAlgebra
 using ReducedComplexityModeling
 using ReducedComplexityModeling: POD, reduce!, compose
+using Test
 
 begin
     # Lorenz problem with POD
@@ -8,8 +10,9 @@ begin
     # Reduced Basis Model Pipeline
     #------------------------------------------------------------------------------
 
-    # 1. Get data 
-    prob, sol = lorenz_problem()
+    # 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))
     data = GeometricIntegratorData(sol)
 
     # 2. Define a model with an appropriate algorithm.
@@ -24,6 +27,16 @@ begin
     red_sol = integrate(red_prob, ExplicitEuler())
     psol = project(red_sol, model)
 
+    # With k = 3 the basis is an orthogonal 3 × 3 matrix, and the explicit Euler steps of the
+    # reduced problem are those of the full problem in rotated coordinates, so the projected
+    # solution is the snapshot of the full integration. The relative error is 1.5e-15.
+    @testset "POD with k = n reproduces the full solution" begin
+        R = model.reduced_basis
+        @test R' * R ≈ I(k) atol=1e-12
+        @test size(psol) == size(data.snapshot)
+        @test psol ≈ data.snapshot rtol=1e-12
+    end
+
     #------------------------------------------------------------------------------
 
     # Visualization and Shennanigans
diff --git a/test/parameters/parameter.jl b/test/parameters/parameter.jl
index ab4d734..f71074b 100644
--- a/test/parameters/parameter.jl
+++ b/test/parameters/parameter.jl
@@ -2,7 +2,7 @@ using ReducedComplexityModeling
 using Test
 
 samples = [0.0, 1.0, 2.0, 3.0, 4.0]
-h5file = "temp.h5"
+h5file = joinpath(mktempdir(), "temp.h5")
 
 @testset "Parameter" begin
     @test_throws AssertionError Parameter(:μ, 1.0, 4.0, samples)
@@ -76,6 +76,73 @@ h5file = "temp.h5"
     @test p1 == p2
 end
 
+@testset "Parameter keeps the element type of its bounds ($T)" for T in (Float32, Float64)
+    for name in (:mu, "mu")
+        p = Parameter(name, zero(T), one(T), 5)
+        @test p.name === :mu
+        @test minimum(p) isa T
+        @test maximum(p) isa T
+        @test eltype(collect(p)) === T
+        @test collect(p) == T[0, 1 // 4, 1 // 2, 3 // 4, 1]
+    end
+
+    # bounds across the range of `T`; each gives `T` samples between the bounds
+    for (a, b) in ((-zero(T), nextfloat(zero(T))), (floatmin(T), eps(T)),
+        (-floatmax(T) / 2, floatmax(T) / 2), (-one(T), floatmax(T)),
+        ((T(10)^k .* (-one(T), one(T)) for k in -8:4:8)...))
+        p = Parameter(:mu, a, b, 3)
+        @test eltype(collect(p)) === T
+        @test all(a .≤ collect(p) .≤ b)
+    end
+end
+
+@testset "Every form takes a Symbol or a string name" begin
+    for name in (:mu, "mu")
+        p = Parameter(name, 0.0f0, 1.0f0, 5)
+        @test p isa Parameter{Float32, Vector{Float32}}
+        @test p.name === :mu
+        @test collect(p) == Float32[0, 1 // 4, 1 // 2, 3 // 4, 1]
+
+        p = Parameter(name, 0.0, 1.0, 3)
+        @test p isa Parameter{Float64, Vector{Float64}}
+        @test p.name === :mu
+        @test collect(p) == [0.0, 0.5, 1.0]
+
+        p = Parameter(name, 0.0, 1.0)
+        @test p isa Parameter{Float64, Nothing}
+        @test p.name === :mu
+
+        p = Parameter(name, [0.0, 1.0])
+        @test p isa Parameter{Float64, Vector{Float64}}
+        @test p.name === :mu
+        @test (minimum(p), maximum(p), collect(p)) == (0.0, 1.0, [0.0, 1.0])
+
+        p = Parameter(name, 0.0, 1.0, [0.5])
+        @test p isa Parameter{Float64, Vector{Float64}}
+        @test p.name === :mu
+        @test collect(p) == [0.5]
+    end
+    @test_throws MethodError Parameter(:a, 0, 10, 5)
+end
+
+@testset "A count gives the samples of LinRange without duplicates" begin
+    for name in (:μ, "μ")
+        @test collect(Parameter(name, 1.0, 1.0, 5)) == [1.0]
+        @test length(Parameter(name, 1.0, 1.0, 5)) == 1
+        # `LinRange(0.1, 0.1, n)` rounds above the bounds for these counts
+        for n in (6, 8, 11, 15)
+            p = Parameter(name, 0.1, 0.1, n)
+            @test p isa Parameter{Float64, Vector{Float64}}
+            @test collect(p) == [0.1]
+        end
+        @test collect(Parameter(name, 0.1f0, 0.1f0, 6)) == Float32[0.1]
+        @test collect(Parameter(name, 0.1, 0.1, 0)) == Float64[]
+        @test collect(Parameter(name, 1, 1, 3)) == [1]
+        @test collect(Parameter(name, 1, 1, 0)) == Int[]
+        @test_throws MethodError Parameter(name, 0, 1, 3)
+    end
+end
+
 @testset "Parameter show" begin
     p = Parameter(:μ, 0.0, 1.0, 3)
     @test sprint(show, MIME"text/plain"(), p) ==
diff --git a/test/parameters/parametersampler.jl b/test/parameters/parametersampler.jl
index 08ec098..64c5450 100644
--- a/test/parameters/parametersampler.jl
+++ b/test/parameters/parametersampler.jl
@@ -69,8 +69,22 @@ end
     @test_throws AssertionError sample(CartesianParameterSampler(), (a = params.μ,))
 end
 
+@testset "Each sampler keeps the element type of the parameters ($T)" for T in (Float32, Float64)
+    μ = Parameter(:μ, zero(T), one(T), 3)
+    ν = Parameter(:ν, -one(T), T(2), 2)
+
+    @testset "$(nameof(typeof(sampler)))" for sampler in (CartesianParameterSampler(),
+        RandomParameterSampler(5, Xoshiro(1)), QuasiRandomParameterSampler(5))
+        s = sample(sampler, μ, ν)
+        @test eltype(s.μ) === T
+        @test eltype(s.ν) === T
+        @test all(minimum(μ) .≤ s.μ .≤ maximum(μ))
+        @test all(minimum(ν) .≤ s.ν .≤ maximum(ν))
+    end
+end
+
 @testset "h5save and h5load of a randomly sampled ParameterSpace" begin
-    h5file = "temp.h5"
+    h5file = joinpath(mktempdir(), "temp.h5")
     ps = ParameterSpace(RandomParameterSampler(4, Xoshiro(1)),
         Parameter(:μ, 0.0, 1.0), Parameter(:ν, 0.0, 2.0))
 
diff --git a/test/parameters/parameterspace.jl b/test/parameters/parameterspace.jl
index b98efd2..cb91487 100644
--- a/test/parameters/parameterspace.jl
+++ b/test/parameters/parameterspace.jl
@@ -8,7 +8,7 @@ grid = [0.0 1.0 0.0
         0.0 1.0 4.0
         0.5 1.0 4.0
         1.0 1.0 4.0]
-h5file = "temp.h5"
+h5file = joinpath(mktempdir(), "temp.h5")
 
 @testset "ParameterSpace" begin
     μ = Parameter(:μ, 0.0, 1.0, 3)
diff --git a/test/quality/aqua.jl b/test/quality/aqua.jl
index d4a198e..9cfc5ac 100644
--- a/test/quality/aqua.jl
+++ b/test/quality/aqua.jl
@@ -3,7 +3,7 @@ using Aqua
 using Test
 
 @testset "Aqua" begin
-    Aqua.test_all(ReducedComplexityModeling; ambiguities = (; broken = true))   # issue #38
+    Aqua.test_all(ReducedComplexityModeling)
 end
 
 # The names that GeometricMachineLearning and ReducedBasisMethods reach: public, so that a check of
diff --git a/test/runtests.jl b/test/runtests.jl
index 85b0f27..13341f5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,6 +10,8 @@ if "core" in GROUPS
     @safetestset "Parameter samplers" include("parameters/parametersampler.jl")
     @safetestset "Parameter spaces" include("parameters/parameterspace.jl")
     @safetestset "TrainingData" include("TrainingData/TrainingData.jl")
+    @safetestset "Snapshots of geometric solutions" include("TrainingData/GeometricData.jl")
+    @safetestset "Reduced bases" include("Models/ReducedBasis.jl")
     @safetestset "TrainingProblem" include("TrainingProblem.jl")
     @safetestset "The loader's names and constructors" include("data_loader/data_loader.jl")
     @safetestset "Batching" include("data_loader/batch.jl")