From 0eaa1ae241a92f6d2742d514dc51a28292e985ab Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Thu, 29 Oct 2026 20:24:42 +0100 Subject: [PATCH 1/7] Resolve the Parameter ambiguity, document every name, drop export Array Aqua runs with no broken keyword: a Parameter method for a string name and a sample count resolves the ambiguity of issue #38. A test asserts that every exported or public name has a docstring with text and is owned by the package; the module and twelve names gain docstrings, and export Array goes. docs/make.jl sets checkdocs = :all. Tests pin the element type of Parameter and of each sampler in Float32 and Float64. KNOWN_ISSUES records K16 to K19, found on the way. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 23 +++++++++++++ KNOWN_ISSUES.md | 51 +++++++++++++++++++++++++++-- docs/make.jl | 1 + src/Models/ReducedBasis.jl | 28 ++++++++++++++++ src/ReducedComplexityModeling.jl | 10 +++++- src/TrainingData/GeometricData.jl | 29 ++++++++++++++++ src/parameters/parameter.jl | 24 ++++++++++++++ src/parameters/parametersampler.jl | 9 +++++ test/integration/exports.jl | 19 +++++++++++ test/parameters/parameter.jl | 20 +++++++++++ test/parameters/parametersampler.jl | 14 ++++++++ test/quality/aqua.jl | 2 +- 12 files changed, 226 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0f1d5b6..cb93728 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,29 @@ 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("μ", minimum, maximum, n)`, with a string name and a sample count, returns a + `Parameter` with `n` samples. It threw a `MethodError` before, because two methods matched it + equally well. This was the method ambiguity of issue #38, so Aqua's ambiguity check runs and + passes in the test suite. + +### Changed + +- 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 docstring of `Independent` states that its layout + throws an `ErrorException` (KNOWN_ISSUES K16). + ## [0.4.1] — 2026-10-29 ### Changed diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 16baac8..1a45f8d 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,50 @@ third element is a fallback that no build reads. - **Kind:** defect - **Found:** 2026-10-28 (critic round 1b) + +### K16 · `GeometricIntegratorEnsembleData(solutions, Independent())` throws + +- **Location:** `src/TrainingData/GeometricData.jl`, `Base.Array(solutions, ::Independent)`, which + ends in `Array(Z')` on the three-dimensional array `Z`. +- **Evidence:** with `sol` the Lorenz solution of `test/helpers/problems.jl`, + `GeometricIntegratorEnsembleData(GeometricSolution[sol, sol], Independent())` throws + `ErrorException: adjoint not defined for Array{Float64, 3}`; the same call with `Dependent()` + returns a `3 × 20002` snapshot. No test constructs the `Independent` layout. +- **Kind:** defect +- **Found:** 2026-10-29 + +### K17 · `reduce!` with `CotangentLiftEVD` throws on every call + +- **Location:** `src/Models/ReducedBasis.jl`, `reduce!(model, k, ::CotangentLiftEVD)`, the line + `U, _, _ = svd(cotangent_snap)[:, 1:k]`. +- **Evidence:** `reduce!(ReducedBasisModel(alg, randn(6, 10)), 2, alg)` with + `alg = ReducedComplexityModeling.CotangentLiftEVD()` throws `MethodError: no method matching + getindex(::LinearAlgebra.SVD{…}, ::Colon, ::UnitRange{Int64})`. No test calls the method, and + the mutant that makes its body `error("mutant")` passes `integration/pod_lorenz.jl` and + `TrainingProblem.jl`. +- **Kind:** defect +- **Found:** 2026-10-29 + +### K18 · `reduce!` with `ComplexSVD` throws on every call + +- **Location:** `src/Models/ReducedBasis.jl`, `reduce!(model, k, ::ComplexSVD)`, which assigns + `ϕ, Ψ = real.(U), imag.(U)` and then reads `Φ`. +- **Evidence:** `reduce!(ReducedBasisModel(alg, randn(6, 10)), 2, alg)` with + `alg = ReducedComplexityModeling.ComplexSVD()` throws `UndefVarError: Φ not defined in + ReducedComplexityModeling`. No test calls the method, and the mutant that makes its body + `error("mutant")` passes `integration/pod_lorenz.jl` and `TrainingProblem.jl`. +- **Kind:** defect +- **Found:** 2026-10-29 + +### K19 · The parameter tests write `temp.h5` in the working directory + +- **Location:** `test/parameters/parameter.jl`, `h5file = "temp.h5"`, and + `test/parameters/parametersampler.jl`, the testset `h5save and h5load of a randomly sampled + ParameterSpace`. +- **Evidence:** in a `mutate.jl --warm … --jobs 4` run of four mutants started at once, two logs + show `Parameter: Error During Test` with `unable to determine if temp.h5 is accessible in the + HDF5 format (file may not exist)` and then `read_parameters from a file path: Test Failed` at + `@test !isfile(h5file)`: the processes share one working directory, so one test removes or + finds the file of another. A test run in one process does not show it. +- **Kind:** defect +- **Found:** 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..569f467 100644 --- a/src/Models/ReducedBasis.jl +++ b/src/Models/ReducedBasis.jl @@ -1,3 +1,9 @@ +""" + 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 @@ -20,6 +26,13 @@ mutable struct ReducedBasisModel <: AbstractModel end end +""" + reduce!(model::ReducedBasisModel, k::Int, ::POD) + +Compute a reduced basis of dimension `k` 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. +""" function reduce!(model::ReducedBasisModel, k::Int, ::POD) U, _, _ = svd(model.snapshot) model.reduced_basis = U[:, 1:k] @@ -67,6 +80,14 @@ function reduce!(model::ReducedBasisModel, k::Int, ::ComplexSVD) 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 +105,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..f2ca264 100644 --- a/src/TrainingData/GeometricData.jl +++ b/src/TrainingData/GeometricData.jl @@ -11,6 +11,13 @@ function snapshot_matrix(sol::DataSeries) Array(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,8 +36,23 @@ 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. +""" 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. This layout is not implemented: +`GeometricIntegratorEnsembleData(solutions, Independent())` throws an `ErrorException`, because +the snapshot array has three dimensions and `adjoint` is not defined for it. +""" struct Independent <: AbstractDependence end function Base.Array(solutions::Vector{GeometricSolution}, ::Dependent) @@ -63,6 +85,13 @@ function Base.Array(solutions::Vector{GeometricSolution}, ::Independent) Array(Z') end +""" + GeometricIntegratorEnsembleData(solutions::Vector{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 `solutions` is +`GeometricSolution` itself; a vector with a concrete element type throws a `MethodError`. +""" struct GeometricIntegratorEnsembleData solutions::Vector{GeometricSolution} snapshot::Array diff --git a/src/parameters/parameter.jl b/src/parameters/parameter.jl index 648cc59..249ad4f 100644 --- a/src/parameters/parameter.jl +++ b/src/parameters/parameter.jl @@ -7,7 +7,21 @@ _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. `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 `n` points of `LinRange(minimum, maximum, n)`; their element type is that of the +bounds only for floating-point bounds, so integer bounds with a count throw a `MethodError`. 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 @@ -36,6 +50,11 @@ end function Parameter(name, minimum::DT, maximum::DT, n::Int) where {DT} Parameter(name, minimum, maximum, LinRange(minimum, maximum, n)) end +# A string name with a sample count matches the string-name method and the sample-count method +# above equally well; this method decides between them. +function Parameter(name::AbstractString, minimum::DT, maximum::DT, n::Int) where {DT} + Parameter(Symbol(name), minimum, maximum, n) +end function Parameter(name, samples::AbstractVector) Parameter(name, minimum(samples), maximum(samples), samples) end @@ -73,6 +92,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/test/integration/exports.jl b/test/integration/exports.jl index 2ffbc3b..2e7c99a 100644 --- a/test/integration/exports.jl +++ b/test/integration/exports.jl @@ -5,3 +5,22 @@ using Test @test isempty(filter( n -> !isdefined(ReducedComplexityModeling, n), names(ReducedComplexityModeling))) end + +# 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. +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. +@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 !isempty(docstring_text(ReducedComplexityModeling, name)) + @test parentmodule(getfield(ReducedComplexityModeling, name)) === + ReducedComplexityModeling + end +end diff --git a/test/parameters/parameter.jl b/test/parameters/parameter.jl index ab4d734..c085bfe 100644 --- a/test/parameters/parameter.jl +++ b/test/parameters/parameter.jl @@ -76,6 +76,26 @@ 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 "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..55db814 100644 --- a/test/parameters/parametersampler.jl +++ b/test/parameters/parametersampler.jl @@ -69,6 +69,20 @@ 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" ps = ParameterSpace(RandomParameterSampler(4, Xoshiro(1)), 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 From 126c6b03045c2f9b4cf6dfea10b2dadd080aecea Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Thu, 29 Oct 2026 21:04:36 +0100 Subject: [PATCH 2/7] Correct the Parameter count docstring; record K20 and K21 The count form keeps the bounds' type for every type that LinRange keeps, Rational included, not only for floats. K20 and K21 record two defects that the critics found and that exist on origin/main. Co-Authored-By: Claude Opus 5.5 --- KNOWN_ISSUES.md | 22 ++++++++++++++++++++++ src/parameters/parameter.jl | 7 ++++--- 2 files changed, 26 insertions(+), 3 deletions(-) diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 1a45f8d..210bd5c 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -158,3 +158,25 @@ finds the file of another. A test run in one process does not show it. - **Kind:** defect - **Found:** 2026-10-29 + +### K20 · The reduced-basis integration test asserts nothing + +- **Location:** `test/integration/pod_lorenz.jl`, the calls `reduce!(model, k, alg)`, + `compose(prob, model)` and `project(red_sol, model)`. +- **Evidence:** `grep -c '@test' test/integration/pod_lorenz.jl` gives `0`. The script fails only + when one of these calls throws, so a `reduce!`, `compose` or `project` that returns a wrong basis + or a wrong state passes it. A test that compares `project(red_sol, model)` with `sol.q` to a + stated tolerance would catch such a change. +- **Kind:** missing test +- **Found:** 2026-10-29 (critic round 1a) + +### K21 · The docstrings of `h5save` and `h5load` name only the `ParameterSpace` method + +- **Location:** `src/parameters/parameterspace.jl`, the docstrings `save ParameterSpace` and + `Load ParameterSpace`. +- **Evidence:** each name has three methods: the `ParameterSpace` method there, the `Parameter` + method in `src/parameters/parameter.jl` and the file-path method in + `src/parameters/h5routines.jl`. Only the `ParameterSpace` method has a docstring, and it gives + neither the arguments nor the other two forms. +- **Kind:** docs +- **Found:** 2026-10-29 (critic rounds 1a and 1b) diff --git a/src/parameters/parameter.jl b/src/parameters/parameter.jl index 249ad4f..3b03423 100644 --- a/src/parameters/parameter.jl +++ b/src/parameters/parameter.jl @@ -17,9 +17,10 @@ samples between them. `name` is a `Symbol` or a string, which is converted to a 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 `n` points of `LinRange(minimum, maximum, n)`; their element type is that of the -bounds only for floating-point bounds, so integer bounds with a count throw a `MethodError`. With -only a vector of samples, the bounds are its smallest and largest entry. +samples are the `n` points of `LinRange(minimum, maximum, n)`. This works for bounds of a type +that `LinRange` keeps, such as a float or a `Rational`; 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`. """ From a9d6a9026072935c5cbfdbfdc6eab6c2fdcfb7e1 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Thu, 29 Oct 2026 21:51:00 +0100 Subject: [PATCH 3/7] Fix the review findings of PR 64 Independent layout: permutedims gives an ny x nx x ns array whose slice i is the snapshot matrix of solution i; unequal lengths throw an ArgumentError. Dependent layout: allocate from the total number of entries. GeometricIntegratorEnsembleData and both Base.Array methods take AbstractVector{<:GeometricSolution}. reduce! with CotangentLiftEVD and ComplexSVD: take the U factor and cut it to k columns, fix the Phi spelling and the lower-right block index. Tests check orthonormality, the symplectic identity and the reconstruction of a rank-k snapshot. Parameter: the inner constructor takes a Symbol or a string name; the two string-name outer methods are removed. The docstring states that equal bounds give one sample. Docstrings of ReducedBasisModel, h5save, h5load and reduce! name their methods. The exports test requires each docstring to contain its name and drops the redundant isdefined testset. The POD Lorenz test asserts the projected solution on (0, 1). The parameter tests write into a temporary directory. KNOWN_ISSUES K16-K21 are removed. Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 49 +++++++++++++++----- KNOWN_ISSUES.md | 69 ----------------------------- src/Models/ReducedBasis.jl | 30 ++++++++----- src/TrainingData/GeometricData.jl | 32 +++++++------ src/parameters/parameter.jl | 25 ++++------- src/parameters/parameterspace.jl | 18 +++++++- test/Models/ReducedBasis.jl | 31 +++++++++++++ test/TrainingData/GeometricData.jl | 45 +++++++++++++++++++ test/helpers/problems.jl | 3 +- test/integration/exports.jl | 13 +++--- test/integration/pod_lorenz.jl | 17 ++++++- test/parameters/parameter.jl | 39 +++++++++++++++- test/parameters/parametersampler.jl | 2 +- test/parameters/parameterspace.jl | 2 +- test/runtests.jl | 2 + 15 files changed, 237 insertions(+), 140 deletions(-) create mode 100644 test/Models/ReducedBasis.jl create mode 100644 test/TrainingData/GeometricData.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index cb93728..f17a522 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,18 +22,43 @@ holds the changes since the last release until that release is cut. ### Bug Fixes -- `Parameter("μ", minimum, maximum, n)`, with a string name and a sample count, returns a - `Parameter` with `n` samples. It threw a `MethodError` before, because two methods matched it - equally well. This was the method ambiguity of issue #38, so Aqua's ambiguity check runs and - passes in the test suite. - -### Changed - -- 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 docstring of `Independent` states that its layout - throws an `ErrorException` (KNOWN_ISSUES K16). +- `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)`, without duplicates, so equal bounds give one 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: integer bounds with a sample count still throw a + `MethodError`, because the `LinRange` of integers 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`. + +- `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`. + +- `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), and the docstring of + `ReducedBasisModel` no longer names packages that are not dependencies of this package. ## [0.4.1] — 2026-10-29 diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 210bd5c..9f195e3 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -111,72 +111,3 @@ third element is a fallback that no build reads. - **Kind:** defect - **Found:** 2026-10-28 (critic round 1b) - -### K16 · `GeometricIntegratorEnsembleData(solutions, Independent())` throws - -- **Location:** `src/TrainingData/GeometricData.jl`, `Base.Array(solutions, ::Independent)`, which - ends in `Array(Z')` on the three-dimensional array `Z`. -- **Evidence:** with `sol` the Lorenz solution of `test/helpers/problems.jl`, - `GeometricIntegratorEnsembleData(GeometricSolution[sol, sol], Independent())` throws - `ErrorException: adjoint not defined for Array{Float64, 3}`; the same call with `Dependent()` - returns a `3 × 20002` snapshot. No test constructs the `Independent` layout. -- **Kind:** defect -- **Found:** 2026-10-29 - -### K17 · `reduce!` with `CotangentLiftEVD` throws on every call - -- **Location:** `src/Models/ReducedBasis.jl`, `reduce!(model, k, ::CotangentLiftEVD)`, the line - `U, _, _ = svd(cotangent_snap)[:, 1:k]`. -- **Evidence:** `reduce!(ReducedBasisModel(alg, randn(6, 10)), 2, alg)` with - `alg = ReducedComplexityModeling.CotangentLiftEVD()` throws `MethodError: no method matching - getindex(::LinearAlgebra.SVD{…}, ::Colon, ::UnitRange{Int64})`. No test calls the method, and - the mutant that makes its body `error("mutant")` passes `integration/pod_lorenz.jl` and - `TrainingProblem.jl`. -- **Kind:** defect -- **Found:** 2026-10-29 - -### K18 · `reduce!` with `ComplexSVD` throws on every call - -- **Location:** `src/Models/ReducedBasis.jl`, `reduce!(model, k, ::ComplexSVD)`, which assigns - `ϕ, Ψ = real.(U), imag.(U)` and then reads `Φ`. -- **Evidence:** `reduce!(ReducedBasisModel(alg, randn(6, 10)), 2, alg)` with - `alg = ReducedComplexityModeling.ComplexSVD()` throws `UndefVarError: Φ not defined in - ReducedComplexityModeling`. No test calls the method, and the mutant that makes its body - `error("mutant")` passes `integration/pod_lorenz.jl` and `TrainingProblem.jl`. -- **Kind:** defect -- **Found:** 2026-10-29 - -### K19 · The parameter tests write `temp.h5` in the working directory - -- **Location:** `test/parameters/parameter.jl`, `h5file = "temp.h5"`, and - `test/parameters/parametersampler.jl`, the testset `h5save and h5load of a randomly sampled - ParameterSpace`. -- **Evidence:** in a `mutate.jl --warm … --jobs 4` run of four mutants started at once, two logs - show `Parameter: Error During Test` with `unable to determine if temp.h5 is accessible in the - HDF5 format (file may not exist)` and then `read_parameters from a file path: Test Failed` at - `@test !isfile(h5file)`: the processes share one working directory, so one test removes or - finds the file of another. A test run in one process does not show it. -- **Kind:** defect -- **Found:** 2026-10-29 - -### K20 · The reduced-basis integration test asserts nothing - -- **Location:** `test/integration/pod_lorenz.jl`, the calls `reduce!(model, k, alg)`, - `compose(prob, model)` and `project(red_sol, model)`. -- **Evidence:** `grep -c '@test' test/integration/pod_lorenz.jl` gives `0`. The script fails only - when one of these calls throws, so a `reduce!`, `compose` or `project` that returns a wrong basis - or a wrong state passes it. A test that compares `project(red_sol, model)` with `sol.q` to a - stated tolerance would catch such a change. -- **Kind:** missing test -- **Found:** 2026-10-29 (critic round 1a) - -### K21 · The docstrings of `h5save` and `h5load` name only the `ParameterSpace` method - -- **Location:** `src/parameters/parameterspace.jl`, the docstrings `save ParameterSpace` and - `Load ParameterSpace`. -- **Evidence:** each name has three methods: the `ParameterSpace` method there, the `Parameter` - method in `src/parameters/parameter.jl` and the file-path method in - `src/parameters/h5routines.jl`. Only the `ParameterSpace` method has a docstring, and it gives - neither the arguments nor the other two forms. -- **Kind:** docs -- **Found:** 2026-10-29 (critic rounds 1a and 1b) diff --git a/src/Models/ReducedBasis.jl b/src/Models/ReducedBasis.jl index 569f467..d9e598a 100644 --- a/src/Models/ReducedBasis.jl +++ b/src/Models/ReducedBasis.jl @@ -10,11 +10,11 @@ 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 @@ -28,10 +28,16 @@ end """ reduce!(model::ReducedBasisModel, k::Int, ::POD) - -Compute a reduced basis of dimension `k` 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. + 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) @@ -44,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 @@ -70,13 +76,13 @@ 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 diff --git a/src/TrainingData/GeometricData.jl b/src/TrainingData/GeometricData.jl index f2ca264..8edf34f 100644 --- a/src/TrainingData/GeometricData.jl +++ b/src/TrainingData/GeometricData.jl @@ -41,7 +41,8 @@ abstract type AbstractDependence end 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. +first solution, then those of the second, and so on, as its columns. The solutions may have +different numbers of entries in `q`. """ struct Dependent <: AbstractDependence end @@ -49,17 +50,18 @@ 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. This layout is not implemented: -`GeometricIntegratorEnsembleData(solutions, Independent())` throws an `ErrorException`, because -the snapshot array has three dimensions and `adjoint` is not defined for it. +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`; 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) +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) + Z = zeros(eltype(solutions[1].q[1]), nx, ny) i = 1 for solution in solutions @@ -71,10 +73,12 @@ function Base.Array(solutions::Vector{GeometricSolution}, ::Dependent) Array(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]) + all(solution -> length(solution.q) == nx, solutions) || + throw(ArgumentError("the Independent layout needs solutions with equal numbers of entries")) Z = zeros(eltype(solutions[1].q[1]), nx, ny, ns) for (i, solution) in enumerate(solutions) @@ -82,21 +86,21 @@ function Base.Array(solutions::Vector{GeometricSolution}, ::Independent) Z[j, :, i] .= elem end end - Array(Z') + permutedims(Z, (2, 1, 3)) end """ - GeometricIntegratorEnsembleData(solutions::Vector{GeometricSolution}, deps) + 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 `solutions` is -`GeometricSolution` itself; a vector with a concrete element type throws a `MethodError`. +`deps` selects: [`Dependent`](@ref) or [`Independent`](@ref). """ 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 3b03423..fffbdb1 100644 --- a/src/parameters/parameter.jl +++ b/src/parameters/parameter.jl @@ -13,14 +13,15 @@ _sort(s::AbstractVector) = sort(unique(collect(s))) Parameter(name, samples::AbstractVector) A named parameter of a system: the bounds `minimum ≤ maximum` and, optionally, a vector of -samples between them. `name` is a `Symbol` or a string, which is converted to a `Symbol`. +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 `n` points of `LinRange(minimum, maximum, n)`. This works for bounds of a type -that `LinRange` keeps, such as a float or a `Rational`; 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: equal bounds give +one sample. This works for bounds of a type that `LinRange` keeps, such as a float or a +`Rational`; 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`. """ @@ -30,32 +31,24 @@ 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)) end -# A string name with a sample count matches the string-name method and the sample-count method -# above equally well; this method decides between them. -function Parameter(name::AbstractString, minimum::DT, maximum::DT, n::Int) where {DT} - Parameter(Symbol(name), minimum, maximum, n) -end function Parameter(name, samples::AbstractVector) Parameter(name, minimum(samples), maximum(samples), samples) end diff --git a/src/parameters/parameterspace.jl b/src/parameters/parameterspace.jl index 28c03ac..1208827 100644 --- a/src/parameters/parameterspace.jl +++ b/src/parameters/parameterspace.jl @@ -68,7 +68,15 @@ 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. """ function h5save(h5::H5DataStore, ps::ParameterSpace; path::AbstractString = "/") cols = columns(ps.samples) @@ -85,7 +93,13 @@ 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. """ 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/TrainingData/GeometricData.jl b/test/TrainingData/GeometricData.jl new file mode 100644 index 0000000..e0c1f13 --- /dev/null +++ b/test/TrainingData/GeometricData.jl @@ -0,0 +1,45 @@ +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()) +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 "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 2e7c99a..7a31f6e 100644 --- a/test/integration/exports.jl +++ b/test/integration/exports.jl @@ -1,13 +1,9 @@ using ReducedComplexityModeling using Test -@testset "Every exported name is defined" begin - @test isempty(filter( - n -> !isdefined(ReducedComplexityModeling, n), names(ReducedComplexityModeling))) -end - # 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. +# 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) @@ -15,11 +11,12 @@ function docstring_text(mod::Module, name::Symbol) 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. +# `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 !isempty(docstring_text(ReducedComplexityModeling, name)) + @test occursin(string(name), docstring_text(ReducedComplexityModeling, name)) @test parentmodule(getfield(ReducedComplexityModeling, name)) === ReducedComplexityModeling 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 c085bfe..ac46035 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) @@ -96,6 +96,43 @@ end 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 + @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 55db814..64c5450 100644 --- a/test/parameters/parametersampler.jl +++ b/test/parameters/parametersampler.jl @@ -84,7 +84,7 @@ 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/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") From 5ecd6046466a0a41066ced58eb925b159c14be37 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Thu, 29 Oct 2026 22:18:26 +0100 Subject: [PATCH 4/7] Fix test deps, equal bounds and snapshot eltype Declare LinearAlgebra in test/Project.toml: Pkg.test() failed to load it in test/Models/ReducedBasis.jl, which the TestEnv runs hid. Parameter with equal bounds and a positive count gives the one sample minimum, because LinRange(0.1, 0.1, 6) rounds above the bounds. The Independent layout rejects solutions with different state dimensions. Both layouts promote the element type across all solutions, so a Float64 value after a Float32 solution keeps its value. The h5save and h5load docstrings show a round trip as a doctest. Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 27 +++++++++++++++++++-------- src/TrainingData/GeometricData.jl | 18 +++++++++++++----- src/parameters/parameter.jl | 13 ++++++++----- src/parameters/parameterspace.jl | 25 +++++++++++++++++++++++++ test/Project.toml | 1 + test/TrainingData/GeometricData.jl | 27 +++++++++++++++++++++++++++ test/parameters/parameter.jl | 8 ++++++++ 7 files changed, 101 insertions(+), 18 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f17a522..bed5440 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,24 +24,31 @@ holds the changes since the last release until that release is cut. - `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)`, without duplicates, so equal bounds give one 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: integer bounds with a sample count still throw a - `MethodError`, because the `LinRange` of integers has a float element type, so - `Parameter(:a, 0, 10, 5)` throws. + `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`. + 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`. +- 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` @@ -57,7 +64,11 @@ holds the changes since the last release until that release is cut. `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), and the docstring 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 diff --git a/src/TrainingData/GeometricData.jl b/src/TrainingData/GeometricData.jl index 8edf34f..481f82a 100644 --- a/src/TrainingData/GeometricData.jl +++ b/src/TrainingData/GeometricData.jl @@ -53,15 +53,20 @@ The layout of the snapshot of a [`GeometricIntegratorEnsembleData`](@ref) in whi 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`; otherwise the constructor throws an -`ArgumentError`. +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 +# 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]), nx, ny) + Z = zeros(_snapshot_eltype(solutions), nx, ny) i = 1 for solution in solutions @@ -79,7 +84,9 @@ function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Independen ny = length(solutions[1].q[1]) all(solution -> length(solution.q) == nx, solutions) || throw(ArgumentError("the Independent layout needs solutions with equal numbers of entries")) - Z = zeros(eltype(solutions[1].q[1]), nx, ny, ns) + 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) @@ -93,7 +100,8 @@ 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). +`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} diff --git a/src/parameters/parameter.jl b/src/parameters/parameter.jl index fffbdb1..422b27b 100644 --- a/src/parameters/parameter.jl +++ b/src/parameters/parameter.jl @@ -18,10 +18,11 @@ 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 give -one sample. This works for bounds of a type that `LinRange` keeps, such as a float or a -`Rational`; 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. 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`. """ @@ -47,7 +48,9 @@ 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 = n > 0 && minimum == maximum ? [minimum] : LinRange(minimum, maximum, n) + Parameter(name, minimum, maximum, samples) end function Parameter(name, samples::AbstractVector) Parameter(name, minimum(samples), maximum(samples), samples) diff --git a/src/parameters/parameterspace.jl b/src/parameters/parameterspace.jl index 1208827..8e980ac 100644 --- a/src/parameters/parameterspace.jl +++ b/src/parameters/parameterspace.jl @@ -77,6 +77,17 @@ per parameter to the subgroup `samples`, and each of its parameters to the subgr 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) @@ -100,6 +111,20 @@ end 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/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 index e0c1f13..d213995 100644 --- a/test/TrainingData/GeometricData.jl +++ b/test/TrainingData/GeometricData.jl @@ -24,6 +24,33 @@ sol3 = oscillator_solution([3.0, 1.0], 0.05) @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 diff --git a/test/parameters/parameter.jl b/test/parameters/parameter.jl index ac46035..0bd5c90 100644 --- a/test/parameters/parameter.jl +++ b/test/parameters/parameter.jl @@ -129,6 +129,14 @@ end 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_throws MethodError Parameter(name, 0, 1, 3) end end From 4c32e1d939943e29bac669d23c9586f33d709aba Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Thu, 29 Oct 2026 22:27:53 +0100 Subject: [PATCH 5/7] Equal bounds with count 0; record K22 to K27 Equal bounds give fill(minimum, min(n, 1)) for every n >= 0, so Parameter(:a, 1, 1, 0) gives no sample instead of a MethodError from the Float64 LinRange. KNOWN_ISSUES records the findings that this PR does not fix: the four deferred review findings (K22 to K25), the Dependent layout without a state-dimension check (K26) and near-equal bounds (K27). Co-Authored-By: Claude Sonnet 5.5 --- KNOWN_ISSUES.md | 63 ++++++++++++++++++++++++++++++++++++ src/parameters/parameter.jl | 3 +- test/parameters/parameter.jl | 2 ++ 3 files changed, 67 insertions(+), 1 deletion(-) diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 9f195e3..15a6e4e 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -111,3 +111,66 @@ 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) + +### K26 · The `Dependent` layout accepts solutions with different state dimensions + +- **Location:** `src/TrainingData/GeometricData.jl`, `Base.Array(solutions, ::Dependent)`. +- **Evidence:** a solution with two components of `q` and one with one component, 11 entries each, + give a `2 × 22` snapshot without an error. The broadcast copies the one component into both + rows: column 12 is `[1.0, 1.0]`. The `Independent` layout throws an `ArgumentError` for the same + input. +- **Kind:** defect +- **Found:** PR #64 (2026-10-29) + +### 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/src/parameters/parameter.jl b/src/parameters/parameter.jl index 422b27b..39086a2 100644 --- a/src/parameters/parameter.jl +++ b/src/parameters/parameter.jl @@ -49,7 +49,8 @@ function Parameter(name, minimum::DT, maximum::DT) where {DT} 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 - samples = n > 0 && minimum == maximum ? [minimum] : LinRange(minimum, maximum, n) + 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) diff --git a/test/parameters/parameter.jl b/test/parameters/parameter.jl index 0bd5c90..f71074b 100644 --- a/test/parameters/parameter.jl +++ b/test/parameters/parameter.jl @@ -137,6 +137,8 @@ end 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 From a1e3cd44ca8a1b52a90e473f764f34a0b0613c92 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Fri, 30 Oct 2026 07:21:24 +0100 Subject: [PATCH 6/7] Throw an ArgumentError for Dependent solutions with different state dimensions (K26) Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 4 +++- KNOWN_ISSUES.md | 10 ---------- src/TrainingData/GeometricData.jl | 5 ++++- test/TrainingData/GeometricData.jl | 10 ++++++++++ 4 files changed, 17 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bed5440..ad353af 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,7 +42,9 @@ holds the changes since the last release until that release is cut. - `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`. + 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. - 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 diff --git a/KNOWN_ISSUES.md b/KNOWN_ISSUES.md index 15a6e4e..50aefd6 100644 --- a/KNOWN_ISSUES.md +++ b/KNOWN_ISSUES.md @@ -156,16 +156,6 @@ - **Kind:** not verified - **Found:** PR #64 (CodeRabbit review on a9d6a90) -### K26 · The `Dependent` layout accepts solutions with different state dimensions - -- **Location:** `src/TrainingData/GeometricData.jl`, `Base.Array(solutions, ::Dependent)`. -- **Evidence:** a solution with two components of `q` and one with one component, 11 entries each, - give a `2 × 22` snapshot without an error. The broadcast copies the one component into both - rows: column 12 is `[1.0, 1.0]`. The `Independent` layout throws an `ArgumentError` for the same - input. -- **Kind:** defect -- **Found:** PR #64 (2026-10-29) - ### 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)`. diff --git a/src/TrainingData/GeometricData.jl b/src/TrainingData/GeometricData.jl index 481f82a..5c00830 100644 --- a/src/TrainingData/GeometricData.jl +++ b/src/TrainingData/GeometricData.jl @@ -42,7 +42,8 @@ abstract type AbstractDependence end 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`. +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 @@ -66,6 +67,8 @@ 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) || + throw(ArgumentError("the Dependent layout needs solutions with equal state dimensions")) Z = zeros(_snapshot_eltype(solutions), nx, ny) i = 1 diff --git a/test/TrainingData/GeometricData.jl b/test/TrainingData/GeometricData.jl index d213995..67c081c 100644 --- a/test/TrainingData/GeometricData.jl +++ b/test/TrainingData/GeometricData.jl @@ -60,6 +60,16 @@ end 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}) From 7782b404a11315aa5a5e3b20bf0840eed6924393 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Fri, 30 Oct 2026 07:44:10 +0100 Subject: [PATCH 7/7] Change the snapshot layout with permutedims, not the conjugating adjoint Co-Authored-By: Claude Sonnet 5.5 --- CHANGELOG.md | 4 +++- src/TrainingData/GeometricData.jl | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ad353af..6506f3b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,7 +44,9 @@ holds the changes since the last release until that release is cut. 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. + 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 diff --git a/src/TrainingData/GeometricData.jl b/src/TrainingData/GeometricData.jl index 5c00830..bedff56 100644 --- a/src/TrainingData/GeometricData.jl +++ b/src/TrainingData/GeometricData.jl @@ -8,7 +8,7 @@ function snapshot_matrix(sol::DataSeries) for (i, elem) in enumerate(sol) Z[i, :] .= elem end - Array(Z') + permutedims(Z) end """ @@ -78,7 +78,7 @@ function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Dependent) i += 1 end end - Array(Z') + permutedims(Z) end function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Independent)