Fix the deferred findings K22-K25, add CotangentLiftSVD, move the Lorenz plots to a script

.gitignore

@@ -38,3 +38,6 @@ Presentations/*.png
 Presentations/*.pdf
 runs/*.h5
 runs/*.pdf
+
+# Output of the scripts
+/results

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

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)

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 = ".."}

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)

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