Repository navigation
Aqua without broken, documented and owned exports, drop export Array #64
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Changes from all commits
0eaa1ae
126c6b0
a9d6a90
5ecd604
4c32e1d
a1e3cd4
7782b40
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,14 +1,20 @@ | ||
| """ | ||
| POD() | ||
|
|
||
| Proper orthogonal decomposition: `reduce!(model, k, POD())` takes the first `k` left singular | ||
| vectors of the snapshot matrix as the reduced basis. | ||
| """ | ||
| struct POD end | ||
| struct CotangentLiftSVD end | ||
| struct CotangentLiftEVD end | ||
| struct ComplexSVD end | ||
|
|
||
| """ | ||
| Reduced Basis | ||
|
|
||
| Uses ModelOrderReduction methods from ReducedBasisMethods.jl | ||
| to represent the data in latent space. | ||
| ReducedBasisModel(alg, snapshot) | ||
|
|
||
| A reduced-basis model of the snapshot matrix `snapshot`, whose columns are the snapshots, with | ||
| the reduction scheme `alg`, such as [`POD`](@ref). [`reduce!`](@ref) computes the reduced basis | ||
| and stores it in the field `reduced_basis`, which is undefined until then. | ||
| """ | ||
| mutable struct ReducedBasisModel <: AbstractModel | ||
| reduction_scheme::Any | ||
|
|
@@ -20,6 +26,19 @@ mutable struct ReducedBasisModel <: AbstractModel | |
| end | ||
| end | ||
|
|
||
| """ | ||
| reduce!(model::ReducedBasisModel, k::Int, ::POD) | ||
| reduce!(model::ReducedBasisModel, k::Int, ::CotangentLiftEVD) | ||
| reduce!(model::ReducedBasisModel, k::Int, ::ComplexSVD) | ||
|
|
||
| Compute a reduced basis from the snapshot matrix of `model`, whose columns are the snapshots, | ||
| store it in `model.reduced_basis` and return it. With [`POD`](@ref), the basis is the matrix of | ||
| the first `k` left singular vectors of the snapshot matrix. With `CotangentLiftEVD()` and | ||
| `ComplexSVD()`, which the package does not export, the snapshot has `2n` rows, `q` above `p`, and | ||
| the basis is an orthonormal and symplectic `2n × 2k` matrix: the cotangent lift of the first `k` | ||
| left singular vectors of `[q p]`, or the real form of the first `k` left singular vectors of | ||
| `q + im * p`. | ||
| """ | ||
| function reduce!(model::ReducedBasisModel, k::Int, ::POD) | ||
| U, _, _ = svd(model.snapshot) | ||
| model.reduced_basis = U[:, 1:k] | ||
|
|
@@ -31,7 +50,7 @@ function reduce!(model::ReducedBasisModel, k::Int, ::CotangentLiftEVD) | |
| n = floor(Int, size(snapshot, 1)/2) | ||
| q, p = snapshot[1:n, :], snapshot[(n + 1):end, :] | ||
| cotangent_snap = hcat(q, p) | ||
| U, _, _ = svd(cotangent_snap)[:, 1:k] | ||
| U = svd(cotangent_snap).U[:, 1:k] | ||
| reduced_basis = zeros(2*n, 2*k) | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Preserve the basis element type. Both allocations default to Also applies to: 81-81
|
||
| reduced_basis[1:n, 1:k] = U | ||
| reduced_basis[(n + 1):(2 * n), (k + 1):(2 * k)] = U | ||
|
|
@@ -57,16 +76,24 @@ function reduce!(model::ReducedBasisModel, k::Int, ::ComplexSVD) | |
| n = floor(Int, size(snapshot, 1)/2) | ||
| q, p = snapshot[1:n, :], snapshot[(n + 1):end, :] | ||
| csvd_snap = q + 1im * p | ||
| U, _, _ = svd(csvd_snap) | ||
| ϕ, Ψ = real.(U), imag.(U) | ||
| U = svd(csvd_snap).U[:, 1:k] | ||
| Φ, Ψ = real.(U), imag.(U) | ||
| reduced_basis = zeros(2*n, 2*k) | ||
| reduced_basis[1:n, 1:k] = Φ | ||
| reduced_basis[1:n, (k + 1):(2 * k)] = -Ψ | ||
| reduced_basis[(n + 1):(2 * n), 1:k] = Ψ | ||
| reduced_basis[((n + 1):(2 * n):(k + 1)):(2 * k)] = Φ | ||
| reduced_basis[(n + 1):(2 * n), (k + 1):(2 * k)] = Φ | ||
| model.reduced_basis = reduced_basis | ||
| end | ||
|
|
||
| """ | ||
| compose(prob::ODEProblem, model::ReducedBasisModel) | ||
|
|
||
| The reduced problem of `prob` in the reduced basis `R = model.reduced_basis`: an `ODEProblem` in | ||
| the reduced coordinates `u`, with the vector field `R' * v(R * u)`, where `v` is the vector field | ||
| of `prob`, and the initial condition `R' * q₀`. The time span, the time step and the parameters | ||
| are those of `prob`. Call [`reduce!`](@ref) on `model` first. | ||
| """ | ||
| function compose(prob::ODEProblem, model::ReducedBasisModel) | ||
| Rᵦ = model.reduced_basis | ||
|
|
||
|
|
@@ -84,6 +111,13 @@ function compose(prob::ODEProblem, model::ReducedBasisModel) | |
| red_prob | ||
| end | ||
|
|
||
| """ | ||
| project(rsol::GeometricSolution, model::ReducedBasisModel) | ||
|
|
||
| Map the solution `rsol` of a reduced problem (see [`compose`](@ref)) back to the full space. The | ||
| result is the matrix `model.reduced_basis * Q`, where the columns of `Q` are the entries of | ||
| `rsol.q`; its columns are the states in the full space. | ||
| """ | ||
| function project(rsol::GeometricSolution, model::ReducedBasisModel) | ||
| Rᵦ = model.reduced_basis | ||
| Rᵦ * snapshot_matrix(rsol.q) | ||
|
|
||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -8,9 +8,16 @@ function snapshot_matrix(sol::DataSeries) | |
| for (i, elem) in enumerate(sol) | ||
| Z[i, :] .= elem | ||
| end | ||
| Array(Z') | ||
| permutedims(Z) | ||
| end | ||
|
|
||
| """ | ||
| GeometricIntegratorData(solution::GeometricSolution) | ||
|
|
||
| The solution of a geometric integrator and its snapshot matrix: the field `snapshot` is a matrix | ||
| with one row per component of the state `q` and one column per entry of `solution.q`, the | ||
| initial state included. | ||
| """ | ||
| struct GeometricIntegratorData | ||
| solution::GeometricSolution | ||
| snapshot::Matrix | ||
|
|
@@ -29,15 +36,40 @@ end | |
|
|
||
| abstract type AbstractDependence end | ||
|
|
||
| """ | ||
| Dependent() | ||
|
|
||
| The layout of the snapshot of a [`GeometricIntegratorEnsembleData`](@ref) in which the solutions | ||
| form one snapshot matrix: one row per component of the state `q`, and the entries of `q` of the | ||
| first solution, then those of the second, and so on, as its columns. The solutions may have | ||
| different numbers of entries in `q`, but all must have the same number of components of `q`; | ||
| otherwise the constructor throws an `ArgumentError`. | ||
| """ | ||
| struct Dependent <: AbstractDependence end | ||
|
|
||
| """ | ||
| Independent() | ||
|
|
||
| The layout of the snapshot of a [`GeometricIntegratorEnsembleData`](@ref) in which each solution | ||
| keeps its own snapshot, along a third dimension: an array of size `ny × nx × ns`, for `ns` | ||
| solutions with `nx` entries of a state `q` with `ny` components. Its slice `[:, :, i]` is the | ||
| snapshot matrix of solution `i`, in the layout of [`GeometricIntegratorData`](@ref). All | ||
| solutions must have the same number of entries in `q` and the same number of components of `q`; | ||
| otherwise the constructor throws an `ArgumentError`. | ||
| """ | ||
| struct Independent <: AbstractDependence end | ||
|
|
||
| function Base.Array(solutions::Vector{GeometricSolution}, ::Dependent) | ||
| ns = length(solutions) | ||
| nx = length(solutions[1].q) | ||
| # the element type that holds the state `q` of every solution without loss | ||
| function _snapshot_eltype(solutions) | ||
| mapreduce(solution -> eltype(solution.q[1]), promote_type, solutions) | ||
| end | ||
|
|
||
| function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Dependent) | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Handle offset-indexed solution vectors. The new Use the first element without assuming its index. Also normalize the constructor input to a one-based vector before storing it. Also applies to: 76-76, 103-103
|
||
| nx = sum(solution -> length(solution.q), solutions) | ||
| ny = length(solutions[1].q[1]) | ||
|
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Read the initial entry when a solution has no stored steps. GeometricSolutions permits Read Also applies to: 79-79
|
||
| Z = zeros(eltype(solutions[1].q[1]), ns*nx, ny) | ||
| all(solution -> length(solution.q[1]) == ny, solutions) || | ||
| throw(ArgumentError("the Dependent layout needs solutions with equal state dimensions")) | ||
| Z = zeros(_snapshot_eltype(solutions), nx, ny) | ||
|
|
||
| i = 1 | ||
| for solution in solutions | ||
|
|
@@ -46,28 +78,40 @@ function Base.Array(solutions::Vector{GeometricSolution}, ::Dependent) | |
| i += 1 | ||
| end | ||
| end | ||
| Array(Z') | ||
| permutedims(Z) | ||
| end | ||
|
|
||
| function Base.Array(solutions::Vector{GeometricSolution}, ::Independent) | ||
| function Base.Array(solutions::AbstractVector{<:GeometricSolution}, ::Independent) | ||
| ns = length(solutions) | ||
| nx = length(solutions[1].q) | ||
| ny = length(solutions[1].q[1]) | ||
| Z = zeros(eltype(solutions[1].q[1]), nx, ny, ns) | ||
| all(solution -> length(solution.q) == nx, solutions) || | ||
| throw(ArgumentError("the Independent layout needs solutions with equal numbers of entries")) | ||
|
coderabbitai[bot] marked this conversation as resolved.
Resolved
|
||
| all(solution -> length(solution.q[1]) == ny, solutions) || | ||
| throw(ArgumentError("the Independent layout needs solutions with equal state dimensions")) | ||
| Z = zeros(_snapshot_eltype(solutions), nx, ny, ns) | ||
|
|
||
| for (i, solution) in enumerate(solutions) | ||
| for (j, elem) in enumerate(solution.q) | ||
| Z[j, :, i] .= elem | ||
| end | ||
| end | ||
| Array(Z') | ||
| permutedims(Z, (2, 1, 3)) | ||
| end | ||
|
|
||
| """ | ||
| GeometricIntegratorEnsembleData(solutions::AbstractVector{<:GeometricSolution}, deps) | ||
|
|
||
| The solutions of an ensemble of geometric integrations and their snapshot, in the layout that | ||
| `deps` selects: [`Dependent`](@ref) or [`Independent`](@ref). The element type of the snapshot | ||
| is the promotion of the element types of the states `q` of all solutions. | ||
| """ | ||
| struct GeometricIntegratorEnsembleData | ||
| solutions::Vector{GeometricSolution} | ||
| snapshot::Array | ||
|
|
||
| function GeometricIntegratorEnsembleData(solutions::Vector{GeometricSolution}, deps::AbstractDependence) | ||
| function GeometricIntegratorEnsembleData( | ||
| solutions::AbstractVector{<:GeometricSolution}, deps::AbstractDependence) | ||
| new(solutions, Array(solutions, deps)) | ||
| end | ||
|
|
||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Handle basis sizes larger than the snapshot count.
The default SVD returns only
min(rows, columns)left singular vectors. For example,ComplexSVDwith a6 × 1snapshot andk = 2throws aBoundsError, although a6 × 4basis is possible.CotangentLiftEVDhas the same issue whenkexceeds twice the snapshot count. Use a full SVD when the thin result has too few columns. (docs.julialang.org)Also applies to: 79-79