From 5369e808c781b7fb2b92ca3918c1558d61c1b3b7 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Sun, 5 Feb 2023 22:45:31 -0500 Subject: [PATCH 01/35] print time steps --- src/primitives.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/primitives.jl b/src/primitives.jl index 3d5ae2d..2a76e83 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -71,6 +71,7 @@ function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt) while (prob.tspan[2]-integrator.t)/T > 1e-8 step!(integrator, alg, dt) update_sol!(integrator) + println(integrator.t) end return integrator.sol end From 918e6f752931f89d08eb2cdc39bd402d1f229b1b Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Sun, 5 Feb 2023 22:46:56 -0500 Subject: [PATCH 02/35] print out more informative log --- src/primitives.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/primitives.jl b/src/primitives.jl index 2a76e83..ad54beb 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -66,12 +66,14 @@ end solves the given problem with the specified algorithm and step size """ function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt) + println("Initialize integrator:") integrator = init(prob, alg, dt) + println("Initialization complete. Start Integration.") T = prob.tspan[2] - prob.tspan[1] while (prob.tspan[2]-integrator.t)/T > 1e-8 step!(integrator, alg, dt) update_sol!(integrator) - println(integrator.t) + println("t = $(integrator.t)") end return integrator.sol end From df0f1b3d7acb537adbaabb5cfa8f16eadb0a9536 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Wed, 22 Feb 2023 10:42:54 -0500 Subject: [PATCH 03/35] define sparse interpolation primitives --- src/sparse_interpolation.jl | 389 ++++++++++++++++++++++++++++++++++++ 1 file changed, 389 insertions(+) create mode 100644 src/sparse_interpolation.jl diff --git a/src/sparse_interpolation.jl b/src/sparse_interpolation.jl new file mode 100644 index 0000000..bf91332 --- /dev/null +++ b/src/sparse_interpolation.jl @@ -0,0 +1,389 @@ +using DocStringExtensions, UnPack, LinearAlgebra + +import LinearAlgebra.adjoint +import Base.size +""" + $(TYPEDEF) + + Abstract type identifying sparse interpolation methods. +""" +abstract type SparseInterpolator end + +""" + $(TYPEDEF) + + Abstract type identifying caches for sparse interpolation methods. +""" +abstract type AbstractInterpolatorCache end + +""" + $(TYPEDEF) + + Auxiliary cache for potentially fast in place updates. +""" +struct InterpolatorCache{T} <: AbstractInterpolatorCache + rows::Matrix{T} + cols::Matrix{T} + elements::Matrix{T} +end + +""" + $(TYPEDEF) + + type identifying (vector- or matrix-valued) functions which may be evaluated in components + (individual rows, columns or even elements). + + #Fields + * rows! function that evaluates (dx, x, p, t, row_idcs) -> F(x,p,t)[row_idcs, :] + and stores the result in dx, i.e., evaluation of isolated rows + * cols! function that evaluates (dx, x, p, t, col_idcs) -> F(x,p,t)[:, col_idcs] + and stores the result in dx, i.e., evaluation of isolated columns + * elements! function that evaluates (dx, x, p, t, row_idcs, col_idcs) -> F(x,p,t)[row_idcs, col_idcs] + and stores the result in dx, i.e., evaluation of isolated entries of a given index grid + * n_rows number of rows of the function output + * n_cols number of columns of the function output +""" +struct ComponentFunction{FR,FC,FE} + rows!::FR + cols!::FC + elements!::FE + n_rows::Int + n_cols::Int +end + +function ComponentFunction(n_rows, n_cols = 1; rows = nothing, cols=nothing, elements=nothing) + ComponentFunction(rows, cols, elements, n_rows, n_cols) +end + +""" + $(TYPEDEF) + + Discrete Empirical Interpolation Method (DEIM) Interpolators. +""" +struct DEIMInterpolator{T} <: SparseInterpolator + interpolation_idcs::Vector{Int} + weights::Matrix{T} +end + +function (Π::DEIMInterpolator)(A::AbstractMatrix) + return Π.weights * A[Π.interpolation_idcs,:] +end + +function (Π::DEIMInterpolator)(A::AbstractVector) + return Π.weights * A[Π.interpolation_idcs,:] +end + +""" + $(TYPEDEF) + + Lazy adjoint for `DEIMInterpolator`s. +""" +struct AdjointDEIMInterpolator{P <: DEIMInterpolator} <: SparseInterpolator + parent::P +end + +function (Π::AdjointDEIMInterpolator)(A::AbstractMatrix) + return A[:,Π.parent.interpolation_idcs] * Π.parent.weights' +end + +""" + $(TYPEDSIGNATURES) + + lazy adjoint for `DEIMInterpolator`s. We mean adjoint in the sense that the interpolation is + applied to the corange. +""" +adjoint(Π::DEIMInterpolator) = AdjointDEIMInterpolator(Π) +adjoint(Π::AdjointDEIMInterpolator) = Π.parent + +""" + $(TYPEDEF) + + Sparse interpolator for matrices encoded by interpolators for range and corange. +""" +struct SparseMatrixInterpolator{R <: SparseInterpolator, C <: SparseInterpolator} <: SparseInterpolator + range::R + corange::C +end + +function SparseMatrixInterpolator(row_idcs::Vector{Int}, col_idcs::Vector{Int}, + range_weights::Matrix{T}, corange_weights::Matrix{T}) where T + range = DEIMInterpolator(row_idcs, range_weights) + corange = DEIMInterpolator(col_idcs, corange_weights) + return SparseMatrixInterpolator(range, corange) +end + + +function (Π::SparseMatrixInterpolator)(A::AbstractMatrix) + @unpack range, corange = Π + return range.weights * A[range.interpolation_idcs,corange.interpolation_idcs] * corange.weights' +end + +""" + $(TYPEDEF) + + sparse interpolation for vector-valued function of the form `F!(dx,x,p,t)`. + All computation are assumed to be in-place. + + #Fields + * `F` `ComponentFunction` to be sparsely interpolated + * `interpolator` `SparseInterpolator` that carries information for how to interpolate + rows and columns (specific type determines how and which interpolation is performed) + * `cache` cache for auxiliary memory used for evaluation of rows! + + If interpolator is a `DEIMInterpolator` then it is assumed that this is intended to interpolate only for the range + If interpolator is a `AdjointDEIMInterpolator` then it is assumed that this intended to interpolate only for the corange + If interpolator is a `SparseMatrixInterpolator` then both range and corange are interpolated +""" +mutable struct SparseFunctionInterpolator{IP <: SparseInterpolator, fType <: ComponentFunction, T} <: SparseInterpolator + F::fType + interpolator::IP + cache::InterpolatorCache{T} +end + +function SparseFunctionInterpolator(F, interpolator::DEIMInterpolator; + output::Type = Float64) + r_rows = rank(interpolator) + cache = InterpolatorCache(Matrix{output}(undef, r_rows, F.n_cols), + Matrix{output}(undef, 0, 0), + Matrix{output}(undef, 0, 0)) + return SparseFunctionInterpolator(F, interpolator, cache) +end + +function SparseFunctionInterpolator(F, interpolator::AdjointDEIMInterpolator; + output::Type = Float64) + r_cols = rank(interpolator) + cache = InterpolatorCache(Matrix{output}(undef, 0, 0), + Matrix{output}(undef, F.n_rows, r_cols), + Matrix{output}(undef, 0, 0)) + return SparseFunctionInterpolator(F, interpolator, cache) +end + +function SparseFunctionInterpolator(F, interpolator::SparseMatrixInterpolator; + output::Type = Float64) + r_rows, r_cols = rank(interpolator.range), rank(interpolator.corange) + dim_range = size(interpolator.range.weights, 1) + dim_corange = size(interpolator.corange.weights, 1) + cache = InterpolatorCache(Matrix{output}(undef, r_rows, dim_corange), + Matrix{output}(undef, dim_range, r_cols), + Matrix{output}(undef, r_rows, r_cols)) + return SparseFunctionInterpolator(F, interpolator, cache) +end + +""" + $(TYPEDSIGNATURES) + + makes updates to interpolator information in `SparseFunctionInterpolator` and reinitializes any cached memory. +""" +function update_interpolator!(Π::SparseFunctionInterpolator, interpolator::DEIMInterpolator) + Π.interpolator = interpolator + r_rows = rank(interpolator) + T = eltype(Π.cache.cols) + Π.cache = InterpolationCache(Matrix{T}(undef, r_rows, F.n_cols), + Matrix{T}(undef, 0, 0), + Matrix{T}(undef, 0, 0)) +end + +function update_interpolator!(Π::SparseFunctionInterpolator, interpolator::AdjointDEIMInterpolator) + Π.interpolator = interpolator + r_cols = rank(interpolator) + T = eltype(Π.cache.cols) + Π.cache = InterpolationCache(Matrix{T}(undef, 0, 0), + Matrix{T}(undef, F.n_rows, r_cols), + Matrix{T}(undef, 0, 0)) +end + +function update_interpolator!(Π::SparseFunctionInterpolator, interpolator::SparseMatrixInterpolator) + Π.interpolator = interpolator + r_rows, r_cols = rank(interpolator.range), rank(interpolator.corange) + dim_range = size(interpolator.range.interpolator, 1) + dim_corange = size(interpolator.corange.weights, 1) + T = eltype(Π.cache.cols) + Π.cache = InterpolationCache(Matrix{T}(undef, r_rows, dim_corange), + Matrix{T}(undef, dim_range, r_cols), + Matrix{T}(undef, r_rows, r_cols)) +end + +""" + $(TYPEDSIGNATURES) + + in-place evaluation of interpolation indices. +""" +function eval!(Π::SparseFunctionInterpolator{IP, fType, T}, x, p, t) where {IP <: DEIMInterpolator, fType, T} + @unpack cache = Π + range = Π.interpolator + Π.F.rows!(cache.rows, x, p, t, range.interpolation_idcs) +end + +function eval!(Π::SparseFunctionInterpolator{IP, fType, T}, x, p, t) where {IP <: AdjointDEIMInterpolator, fType, T} + @unpack cache = Π + corange = Π.interpolator.parent + Π.F.cols!(cache.cols, x, p, t, corange.interpolation_idcs) +end + +function eval!(Π::SparseFunctionInterpolator{IP, fType, T}, x, p, t) where {IP <: SparseMatrixInterpolator, fType, T} + @unpack range, corange = Π.interpolator + @unpack cache = Π + Π.F.elements!(cache.elements, x, p, t, range.interpolation_idcs, corange.interpolation_idcs) +end + +""" + $(TYPEDSIGNATURES) + + functor for inplace evaluation of `SparseFunctionInterpolator`s. +""" +function (Π::SparseFunctionInterpolator{IP, fType, T})(dx, x, p, t) where {IP <: DEIMInterpolator, fType, T} + eval!(Π, x, p, t) + mul!(dx, Π.interpolator.weights, Π.cache.rows) +end + +function (Π::SparseFunctionInterpolator{IP, fType, T})(dx, x, p, t) where {IP <: AdjointDEIMInterpolator, fType, T} + eval!(Π, x, p, t) + mul!(dx, Π.cache.cols, Π.interpolator.parent.weights') +end + +function (Π::SparseFunctionInterpolator{IP, fType, T})(dx, x, p, t) where {IP <: SparseMatrixInterpolator, fType, T} + eval!(Π, x, p, t) + mul!(Π.cache.rows, Π.cache.elements, Π.interpolator.corange.weights') + mul!(dx, Π.interpolator.range.weights, Π.cache.rows) +end + + +""" + $(TYPEDSIGNATURES) + + returns rank of sparse interpolators. +""" +rank(Π::DEIMInterpolator) = length(Π.interpolation_idcs) +rank(Π::AdjointDEIMInterpolator) = rank(Π.parent) +rank(Π::SparseMatrixInterpolator) = min(rank(Π.range), rank(Π.corange)) +rank(Π::SparseFunctionInterpolator) = rank(Π.interpolator) + +""" + $(TYPEDSIGNATURES) + + returns dimension of sparse interpolators. +""" +dim(Π::DEIMInterpolator) = rank(Π) +dim(Π::AdjointDEIMInterpolator) = rank(Π.parent) +dim(Π::SparseMatrixInterpolator) = (rank(Π.range), rank(Π.corange)) +dim(Π::SparseFunctionInterpolator) = dim(Π.interpolator) + +""" + $(TYPEDSIGNATURES) + + returns size of interpolations. +""" +size(Π::DEIMInterpolator) = size(Π.weights,1) +size(Π::AdjointDEIMInterpolator) = size(Π.parent) +size(Π::SparseMatrixInterpolator) = (size(Π.range.weights,1), size(Π.corange.weights,1)) +size(Π::SparseFunctionInterpolator) = size(Π.interpolator) + +""" + $(TYPEDEF) + + abstract type characterizing index selection algorithms for sparse interpolation. +""" +abstract type IndexSelectionAlgorithm end + +""" + $(TYPEDEF) + + DEIM index selection algorithm. +""" +struct DEIM end + +""" + $(TYPEDEF) + + QDEIM index selection algorithm. +""" +struct QDEIM end + +""" + $(TYPEDSIGNATURES) + + returns interpolation indices for sparse interpolation. Supports DEIM and QDEIM index selection. +""" +function index_selection(U::Matrix,::DEIM) + m = size(U, 2) + indices = Vector{Int}(undef, m) + @inbounds @views begin + r = abs.(U[:, 1]) + indices[1] = argmax(r) + for l in 2:m + U = U[:, 1:(l - 1)] + P = indices[1:(l - 1)] + PᵀU = U[P, :] + uₗ = U[:, l] + Pᵀuₗ = uₗ[P, :] + c = vec(PᵀU \ Pᵀuₗ) + mul!(r, U, c) + @. r = abs(uₗ - r) + indices[l] = argmax(r) + end + end + return indices +end + +function index_selection(U::Matrix,::QDEIM) + QR = qr(U', ColumnNorm()) + return QR.p[1:size(U,2)] +end + +# a few lil tests +#= +function rows!(dx,x,p,t,idcs) + @views dx .= x[idcs, :] +end +function cols!(dx,x,p,t,idcs) + @views dx .= x[:, idcs] +end +function elements!(dx,x,p,t,idcsA, idcsB) + @views dx .= x[idcsA, idcsB] +end + +F = ComponentFunction(rows!, cols!, elements!, 5, 8) +A = Matrix(reshape(1.0:5*8, (5, 8))) + +row_idcs = [1, 5, 2] +range_weights = zeros(eltype(A), 5, 3) +range_weights[row_idcs,1:3] .= Matrix{eltype(A)}(I,3,3) + +col_idcs = [1, 8, 3, 6] +corange_weights = zeros(eltype(A),8, 4) +corange_weights[col_idcs,1:4] .= Matrix{eltype(A)}(I,4,4) + +Π_range = DEIMInterpolator(row_idcs,range_weights) +Π_corange = DEIMInterpolator(col_idcs,corange_weights) +Π_mat = SparseMatrixInterpolator(Π_range, Π_corange) + +test = all(Π_range(A)[row_idcs,:] .== A[row_idcs,:]) +test = all(Π_corange'(A)[:,col_idcs] .== A[:,col_idcs]) + +F_range_int = SparseFunctionInterpolator(F, Π_range, output = eltype(A)) +F_corange_int = SparseFunctionInterpolator(F, Π_corange', output = eltype(A)) +F_int = SparseFunctionInterpolator(F, Π_mat, output = eltype(A)) + +eval!(F_range_int, A, (), 0.0) +test = all(F_range_int.cache.rows .== A[row_idcs,:]) +eval!(F_corange_int, A, (), 0.0) +test = all(F_corange_int.cache.cols .== A[:,col_idcs]) +eval!(F_int, A, (), 0.0) +test = all(F_int.cache.elements .== A[row_idcs,col_idcs]) + +dA = similar(A) +F_range_int(dA, A, (), 0.0) +test = all(dA[row_idcs, :] .== A[row_idcs, :]) +F_corange_int(dA, A, (), 0.0) +test = all(dA[:, col_idcs] .== A[:, col_idcs]) +F_int(dA, A, (), 0.0) +test = all(dA[row_idcs, col_idcs] .== A[row_idcs, col_idcs]) + +using BenchmarkTools +# non-allocating! +parms = () +t = 0.0 +@btime $F_range_int($dA, $A, $params, $t) +@btime $F_corange_int($dA, $A, $params, $t) +@btime $F_int($dA, $A, $params, $t) +=# \ No newline at end of file From af0b6e293580db277e49ef9dc208fb2deac84b5f Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Wed, 22 Feb 2023 10:45:43 -0500 Subject: [PATCH 04/35] use rel. tolerance for rank adaptation --- src/integrators/rank_adaptive_unconventional.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrators/rank_adaptive_unconventional.jl b/src/integrators/rank_adaptive_unconventional.jl index 13ecfe0..890f72c 100644 --- a/src/integrators/rank_adaptive_unconventional.jl +++ b/src/integrators/rank_adaptive_unconventional.jl @@ -220,7 +220,7 @@ function rankadaptive_unconventional_step!(u, cache, t, dt) step!(SIntegrator, dt, true) U, S, V = svd(SIntegrator.u) - r_new = min(r_max, LowRankArithmetic.truncate_to_tolerance(S, tol)) + r_new = min(r_max, LowRankArithmetic.truncate_to_tolerance(S, tol, rel=true)) if r_new == r mul!(u.U,Uhat,U[:,1:r_new]) u.S .= Matrix(Diagonal(S[1:r_new])) From 8878e2055deaf890e32b4a3a124bdeff110bd993 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Wed, 22 Feb 2023 10:46:21 -0500 Subject: [PATCH 05/35] first steps toward allowing for continuous DEIM with the unconventional integrator --- src/integrators/unconventional.jl | 203 ++++++++++++++++++++++++++++-- 1 file changed, 194 insertions(+), 9 deletions(-) diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index b1283d7..f3ff9a4 100644 --- a/src/integrators/unconventional.jl +++ b/src/integrators/unconventional.jl @@ -1,4 +1,4 @@ -struct UnconventionalAlgorithm_Params{sType, lType, kType} +struct UnconventionalAlgorithm_Params{sType, lType, kType, iType} S_rhs # rhs of S step (core projected rhs) L_rhs # rhs of L step (range projected rhs) K_rhs # rhs of K step (corange projected rhs) @@ -8,6 +8,7 @@ struct UnconventionalAlgorithm_Params{sType, lType, kType} S_alg::sType L_alg::lType K_alg::kType + interpolation::iType end struct UnconventionalAlgorithm{sType, lType, kType} <: AbstractDLRAlgorithm @@ -16,11 +17,39 @@ end function UnconventionalAlgorithm(; S_rhs = nothing, L_rhs = nothing, K_rhs = nothing, S_kwargs = Dict(), L_kwargs = Dict(), K_kwargs = Dict(), S_alg = Tsit5(), L_alg = Tsit5(), K_alg = Tsit5()) - params = UnconventionalAlgorithm_Params(S_rhs, L_rhs, K_rhs, S_kwargs, L_kwargs, K_kwargs, S_alg, L_alg, K_alg) + params = UnconventionalAlgorithm_Params(S_rhs, L_rhs, K_rhs, + S_kwargs, L_kwargs, K_kwargs, + S_alg, L_alg, K_alg, + nothing) return UnconventionalAlgorithm(params) end +function UnconventionalAlgorithm(deim; S_rhs = nothing, L_rhs = nothing, K_rhs = nothing, + S_kwargs = Dict(), L_kwargs = Dict(), K_kwargs = Dict(), + S_alg = Tsit5(), L_alg = Tsit5(), K_alg = Tsit5()) + params = UnconventionalAlgorithm_Params(S_rhs, L_rhs, K_rhs, + S_kwargs, L_kwargs, K_kwargs, + S_alg, L_alg, K_alg, + deim) + return UnconventionalAlgorithm(params) +end + +struct SparseInterpolation + selection_alg + tol + rmin + rmax + init_range + init_corange +end + +function SparseInterpolation(selection_alg, init_range, init_corange; + rmin = 1, rmax = min(size(init_range,1), size(init_corange,1)), tol = eps(Float64)) + return SparseInterpolation(selection_alg, tol, rmin, rmax, init_range, init_corange) +end -struct UnconventionalAlgorithm_Cache{uType,SIntegratorType,LIntegratorType,KIntegratorType,yType} <: AbstractDLRAlgorithm_Cache +struct UnconventionalAlgorithm_Cache{uType, + SIntegratorType,LIntegratorType,KIntegratorType, + yFunc, yType, dType} <: AbstractDLRAlgorithm_Cache US::Matrix{uType} VS::Matrix{uType} M::Matrix{uType} @@ -30,10 +59,19 @@ struct UnconventionalAlgorithm_Cache{uType,SIntegratorType,LIntegratorType,KInte SIntegrator::SIntegratorType LIntegrator::LIntegratorType KIntegrator::KIntegratorType - y + y::yFunc ycurr::yType yprev::yType Δy::yType + interpolation_cache::dType +end + +struct UnconventionalDEIM_Cache + params::SparseInterpolation + Π::SparseFunctionInterpolator + Π_K::SparseFunctionInterpolator + Π_L::SparseFunctionInterpolator + Π_S::SparseFunctionInterpolator end function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) @@ -80,10 +118,90 @@ function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, - nothing, nothing, nothing, nothing) + nothing, nothing, nothing, nothing, nothing) end -function alg_cache(prob::MatrixDataProblem, alg::UnconventionalAlgorithm,u,dt; t0 = prob.tspan[1]) + +function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, + alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) where + {fType <: ComponentFunction} + # allocate memory for frequently accessed arrays + tspan = (t0,t0+dt) + + @unpack tol, rmin, rmax, + init_range, init_corange, + selection_alg = alg.alg_params.interpolation + + row_idcs = index_selection(init_range, selection_alg) + col_idcs = index_selection(init_corange, selection_alg) + + Π_corange = DEIMInterpolator(col_idcs, init_range/init_range[col_idcs,:]) + Π_range = DEIMInterpolator(row_idcs, init_corange/init_corange[row_idcs,:]) + Π = SparseFunctionInterpolator(F, SparseMatrixInterpolator(Π_range, Π_corange)) + + US = u.U*u.S + QRK = qr(US) + M = Matrix(QRK.Q)'*u.U + + VS = u.V*u.S' + QRL = qr(VS) + N = Matrix(QRL.Q)'*u.V + + if isnothing(alg.alg_params.K_rhs) + K_rhs = function (dK, K, p, t) + Π_K, V0, params = p + Π_K(dK, TwoFactorRepresentation(K, V0), params, t) + end + else + K_rhs = alg.alg_params.K_rhs + end + Π_K_corange = DEIMInterpolator(col_idcs, + V0'*Π.interpolator.corange.weights)' + Π_K = SparseFunctionInterpolator(prob.F, Π_K_corange) + p_K = (Π_K, u.V, ()) + KProblem = ODEProblem(K_rhs, US, tspan, p_K) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + + if isnothing(alg.alg_params.L_rhs) + L_rhs = function (dL, L, p, t) + Π_L, U0, params = p + Π_L(dL', TwoFactorRepresentation(U0,L), params, t) + end + else + L_rhs = alg.alg_params.L_rhs + end + Π_L_range = DEIMInterpolator(row_idcs, + U0'*Π.interpolator.range.weights) + Π_K = SparseFunctionInterpolator(prob.F, Π_L_range) + p_K = (Π_L, u.U, ()) + LProblem = ODEProblem(L_rhs, VS, tspan, u.U) + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + + if isnothing(alg.alg_params.S_rhs) + S_rhs = function (dS, S, p, t) + Π_S, U1, V1, params = p + Π_S(dS, SVDLikeRepresentation(U1,S,V1), params, t) + end + else + S_rhs = alg.alg_params.S_rhs + end + Π_S_mat = SparseMatrixInterpolator(row_dics, col_idcs, + u.U'*Π.interpolator.range.weights, + u.V'*Π.interpolator.corange.weights) + Π_S = SparseFunctionInterpolator(prob.F, Π_S_mat) + p_S = (Π_S, u.U, u.V, ()) + + SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, p_S) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + + interpolation_cache = UnconventionalDEIM_Cache(alg.alg_params.interpolation, Π, Π_K, Π_L, Π_S) + return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, + SIntegrator, LIntegrator, KIntegrator, + nothing, nothing, nothing, nothing, interpolation_cache) +end + + +function alg_cache(prob::MatrixDataProblem, alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) n, r = size(u.U) m = size(u.V, 1) @unpack y = prob @@ -103,7 +221,7 @@ function alg_cache(prob::MatrixDataProblem, alg::UnconventionalAlgorithm,u,dt; t return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, - prob.y, ycurr, yprev, Δy) + prob.y, ycurr, yprev, Δy, nothing) end function init(prob::AbstractDLRProblem, alg::UnconventionalAlgorithm, dt) @@ -131,7 +249,8 @@ function unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem}) end function unconventional_step!(u, cache, t, dt) - @unpack US, VS, M, N, QRK, QRL, KIntegrator, LIntegrator, SIntegrator = cache + @unpack US, VS, M, N, QRK, QRL, + KIntegrator, LIntegrator, SIntegrator = cache # K step mul!(US, u.U, u.S) @@ -156,9 +275,75 @@ function unconventional_step!(u, cache, t, dt) u.S .= SIntegrator.u end +function unconventional_deim_step!(u, cache, t, dt) + @unpack US, VS, M, N, QRK, QRL, + KIntegrator, LIntegrator, SIntegrator, + interpolation_cache = cache + @unpack Π, Π_L, Π_K, Π_S = interpolation_cache + + # K step + mul!(US, u.U, u.S) + set_u!(KIntegrator, US) + step!(KIntegrator, dt, true) + US .= KIntegrator.u + QRK = qr!(US) + mul!(M, Matrix(QRK.Q)', u.U) + + # L step + mul!(VS, u.V, u.S') + set_u!(LIntegrator, VS) + step!(LIntegrator, dt, true) + VS .= LIntegrator.u + QRL = qr!(VS) + mul!(N,Matrix(QRL.Q)',u.V) + u.V .= Matrix(QRL.Q) + u.U .= Matrix(QRK.Q) + + # update core interpolator + mul!(Π_S.range.weights, u.U', Π.interpolator.range.weights) + mul!(Π_S.corange.weights, u.V', Π.interpolator.corange.weights) + + # integration + set_u!(SIntegrator, M*u.S*N') + step!(SIntegrator, dt, true) + u.S .= SIntegrator.u + + update_interpolation!(interpolation_cache, u) +end + +function update_interpolation!(interpolation_cache, u) + @unpack params, Π, Π_K, Π_L, Π_S = interpolation_cache + @unpack selection_method, tol, rmin, rmax = params + + eval!(Π_L, u, (), SIntegrator.t) # rows + eval!(Π_K, u, (), SIntegrator.t) # columns + VF = truncated_svd(Π_L.cache.rows, tol=tol, rmin=rmin, rmax=rmax).V # corange from rows + UF = truncated_svd(Π_K.cache.cols, tol=tol, rmin=rmin, rmax=rmax).U # range from cols + + # find interpolation indices + row_idcs = index_selection(UF, selection_method) + col_idcs = index_selection(VF, selection_method) + + # find interpolation weights + @views range_weights = UF/UF[row_idcs,:] + @views corange_weights = VF/VF[col_idcs,:] + + # new interpolators + range = DEIMInterpolator(row_idcs, range_weights) + projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) + corange = DEIMInterpolator(col_idcs, corange_weights) + projected_corange = DEIMInterpolator(row_idcs, u.V'*corange_weights) + + # update function interpolators + update_interpolator!(Π_L, projected_range) + update_interpolator!(Π_K, projected_corange') + update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) + update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) +end + function step!(integrator::DLRIntegrator, ::UnconventionalAlgorithm, dt) @unpack u, t, iter, cache, probType = integrator unconventional_step!(u, cache, t, dt, probType) integrator.t += dt integrator.iter += 1 -end \ No newline at end of file +end From d6b0f4b868907482571476813eac168310a95901 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Sat, 11 Mar 2023 23:42:17 -0500 Subject: [PATCH 06/35] working status --- Project.toml | 1 + src/LowRankIntegrators.jl | 10 +- src/integrators.jl | 7 +- src/integrators/do_integrators.jl | 250 ++++++++++++++++++ src/integrators/greedy_integrator.jl | 22 +- src/integrators/projector_splitting.jl | 18 +- .../rank_adaptive_unconventional.jl | 28 +- src/integrators/unconventional.jl | 241 +++++++++++++---- src/primitives.jl | 91 +++++-- src/sparse_interpolation.jl | 168 +++++++++++- 10 files changed, 714 insertions(+), 122 deletions(-) create mode 100644 src/integrators/do_integrators.jl diff --git a/Project.toml b/Project.toml index f441bef..3879061 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.1.0" [deps] DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LowRankArithmetic = "8016d1d1-7a3b-4ad2-b79e-5984174ed314" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" diff --git a/src/LowRankIntegrators.jl b/src/LowRankIntegrators.jl index d31a6ea..40bb65d 100644 --- a/src/LowRankIntegrators.jl +++ b/src/LowRankIntegrators.jl @@ -1,5 +1,5 @@ module LowRankIntegrators -using Reexport, UnPack +using Reexport, UnPack, DocStringExtensions import DifferentialEquations: step!, set_u!, init @reexport using LinearAlgebra, DifferentialEquations, LowRankArithmetic @@ -7,6 +7,14 @@ import DifferentialEquations: step!, set_u!, init include("utils.jl") export orthonormalize!, GradientDescent, QR, SVD, SecondMomentMatching, normal_component +include("sparse_interpolation.jl") +export SparseFunctionInterpolator, SparseMatrixInterpolator, + DEIMInterpolator, AdjointDEIMInterpolator, + DEIM, QDEIM, LDEIM, + ComponentFunction, SparseInterpolation, + update_interpolator!, eval!, + index_selection + include("primitives.jl") export MatrixDEProblem, MatrixDataProblem, MatrixHybridProblem DLRIntegrator, DLRSolution, diff --git a/src/integrators.jl b/src/integrators.jl index c161d37..840cacd 100644 --- a/src/integrators.jl +++ b/src/integrators.jl @@ -1,5 +1,10 @@ +# basic catch all cases: +rank_DEIM(::AbstractDLRAlgorithm_Cache) = -1 +interpolation_indices(::AbstractDLRAlgorithm_Cache) = (Int[], Int[]) + include("integrators/projector_splitting.jl") include("integrators/unconventional.jl") include("integrators/data_integrator.jl") include("integrators/rank_adaptive_unconventional.jl") -include("integrators/greedy_integrator.jl") \ No newline at end of file +include("integrators/greedy_integrator.jl") + diff --git a/src/integrators/do_integrators.jl b/src/integrators/do_integrators.jl new file mode 100644 index 0000000..3505e6a --- /dev/null +++ b/src/integrators/do_integrators.jl @@ -0,0 +1,250 @@ +# implements basic capabilities for DO integrators +abstract type AbstractRetraction end +struct DirectTimeMarching <: AbstractRetraction end +struct TruncatedSVD <: AbstractRetraction end +struct AlternatingProjection{oType} <: AbstractRetraction + order::oType +end + +struct Primal end +struct Dual end +struct PrimalDual end + +mutable struct DOAlgorithm{ρType} <: AbstractDLRAlgorithm + ρ::ρType # retraction + alg_params # parameters +end + +#= +Direct time marching + U_{n+1} = U_{n} + dt * L_Z*(L_U'*U_{n}) + Z_{n+1} = Z_{n} + dt * ((I-U_{n}*U_{n}')*L_U)(L_Z'*Z*(Z'*Z)^-1) + orthogonalize(U_{n+1}, Z_{n+1}) + repeat +=# + +struct DirectTimeMarching_Params + F + orth_alg +end + +struct DirectTimeMarching_Cache + F + Y + LU + LZ + C +end + +function DirectTimeMarching(F; orthonormalization = QRFact()) + return DOAlgorithm(DirectTimeMarching(), DirectTimeMarching_Params(F, orthonormalization)) +end + +function init(prob::MatrixDEProblem, alg::DOAlgorithm{<:AbstractRetraction}, dt) + t0, tf = prob.tspan + u = deepcopy(prob.u0) + @assert tf > t0 "Integration in reverse time direction is not supported" + # number of steps + n = floor(Int,(tf-t0)/dt) + 1 + # compute more sensible dt # rest will be handled via interpolation/save_at + dt = (tf-t0)/(n-1) + # initialize solution object + sol = DLRSolution(Vector{typeof(prob.u0)}(undef, n), collect(range(t0, tf, length=n))) + sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object + # initialize cache + cache = alg_cache(prob, alg, u, dt) + return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) +end + +function alg_cache(prob::MatrixDEProblem, alg::DOAlgorithm{DirectTimeMarching}, u, dt) + @unpack alg_params = alg + F = (u,t,dt) -> alg_params.F(prob.f, u, t, dt) + L = F(u, prob.tspan[1] ,dt) + rL = rank(L) + r = rank(u) + LU = zeros(rL,r) + LZ = zeros(rL,r) + Y = deepcopy(u) + C = zeros(r, r) + return DirectTimeMarching_Cache(F, Y, LU, LZ, C) +end + +function step!(integrator::DLRIntegrator, alg::DOAlgorithm{DirectTimeMarching}, dt) + @unpack u, t, iter, cache = integrator + direct_time_marching_do_step!(u, cache, t, dt) + orthonormalize!(u, alg.alg_params.orth_alg) + cache.Y.Z .= u.Z + cache.Y.U .= u.U + integrator.t += dt + integrator.iter += 1 +end + +function direct_time_marching_do_step!(u, cache, t, dt) + @unpack Y, LU, LZ, C, F = cache + L = F(Y,t,dt) + + # Z update + mul!(LU, L.U', Y.U) + mul!(u.Z, L.Z, LU, dt, 1) + + # U update + #((I-U*U')*L_U)(L_Z'*Z*(Z'*Z)^-1) + mul!(L.U, Y.U, LU', -1, 1) #(I-U*U')*L_U + mul!(LZ, L.Z', Y.Z) # L_Z'*Z + mul!(C, Y.Z', Y.Z) # covariance computation + LU .= LZ / Symmetric(C) # revisit + mul!(u.U, L.U, LU, dt, 1) +end + +#= +Alternating projection + Primal step: + LRA_{n+1} = LRA_{n} + dt*F(LRA_{n},t) + solve min \|LRA{n+1} - U_{n} Z_{n+1}'\| => Z_{n+1}' = U_{n}'*LRA_{n+1} => Z_{n+1} = LRA_{n+1}'*U{n} + solve min \|LRA{n+1} - U_{n+1} Z_{n+1}'\| + => procrustes problem: Q, _, P = svd(LRA_{n+1}*Z_{n+1}), U_{n+1}= Q*P' + Dual step: + LRA_{n+1} = LRA_{n} + dt*F(LRA_{n},t) + solve min \|LRA{n+1} - U_{n+1} Z_{n}'\| + => procrustes problem: Q, _, P = svd(LRA_{n+1}*Z_{n}), U_{n+1}= Q*P' + solve min \|LRA{n+1} - U_{n+1} Z_{n+1}'\| => Z_{n+1} = U_{n+1}'*LRA_{n+1} +=# + +struct AlternatingProjection_Params + F + orth_alg +end + +struct AlternatingProjection_Cache + F + UtU + UQ + C + orth_alg +end + +function AlternatingProjection(order, F; orthonormalization = QRFact()) + return DOAlgorithm(AlternatingProjection(order), AlternatingProjection_Params(F, orthonormalization)) +end + +function alg_cache(prob::MatrixDEProblem, alg::DOAlgorithm{AlternatingProjection{T}}, u, dt) where T + @unpack alg_params = alg + F = (u,t,dt) -> alg_params.F(prob.f, u, t, dt) + Y = u + F(u, prob.tspan[1] ,dt) + rY = rank(Y) + r = rank(u) + UtU = zeros(rY, r) + C = zeros(rY, r) + UQ = zeros(size(Y,1), r) + return AlternatingProjection_Cache(F, UtU, UQ, C, alg_params.orth_alg) +end + +function step!(integrator::DLRIntegrator, alg::DOAlgorithm{AlternatingProjection{Primal}}, dt) + @unpack u, t, iter, cache = integrator + primal_alternating_projection_do_step!(u, cache, t, dt) + integrator.t += dt + integrator.iter += 1 +end + +function step!(integrator::DLRIntegrator, alg::DOAlgorithm{AlternatingProjection{Dual}}, dt) + @unpack u, t, iter, cache = integrator + dual_alternating_projection_do_step!(u, cache, t, dt) + integrator.t += dt + integrator.iter += 1 +end + +function step!(integrator::DLRIntegrator, alg::DOAlgorithm{AlternatingProjection{PrimalDual}}, dt) + @unpack u, t, iter, cache = integrator + primal_alternating_projection_do_step!(u, cache, t, dt/2) + dual_alternating_projection_do_step!(u, cache, t + dt/2, dt/2) + integrator.t += dt + integrator.iter += 1 +end + +function primal_alternating_projection_do_step!(u, cache, t, dt) + @unpack F, UtU, UQ, C, orth_alg = cache + Y = u + dt*F(u,t,dt) # investigate using the integrator interface of DifferentialEquations.jl + orthonormalize!(Y, orth_alg) + + # coefficient update + mul!(UtU, Y.U', u.U) + mul!(u.Z, Y.Z, UtU) + + # basis update (procrustes problem) + mul!(C, Y.Z', u.Z) + Q, _, P = svd!(C) + mul!(UQ, Y.U, Q) + mul!(u.U, UQ, P') +end + +function dual_alternating_projection_do_step!(u, cache, t, dt) + @unpack F, UtU, UQ, C, orth_alg = cache + Y = u + dt*F(u,t,dt) # investigate using the integrator interface of DifferentialEquations.jl + orthonormalize!(Y, orth_alg) + + # basis update (procrustes problem) + mul!(C, Y.Z', u.Z) + Q, _, P = svd!(C) + mul!(UQ, Y.U, Q) + mul!(u.U, UQ, P') + + # coefficient update + mul!(UtU, Y.U', u.U) + mul!(u.Z, Y.Z, UtU) +end + +#= +SVD => (svd, tsvd, gradient flow) + LRA_{n+1} = LRA_{n} + dt*F(LRA_{n},t) + LRA_{n+1} = truncated_svd(LRA_{n+1}, rmax = rank(LRA_{n+1})) + or + LRA_{n+1} = truncated_svd(LRA_{n+1}, tol = #tolerance, rmax = #max rank) +=# + +struct TruncatedSVD_Params + F + SVD_alg + orth_alg + tol + rmin + rmax +end + +struct TruncatedSVD_Cache + F +end + +function TruncatedSVD(F; SVD_alg = SVDFact(), orthonormalization = QRFact(), + rmin::Int=1, rmax::Int=Int(1e8), tol = sqrt(eps(Float64))) + return DOAlgorithm(TruncatedSVD(), TruncatedSVD_Params(F, SVD_alg, orthonormalization, tol, rmin, rmax)) +end + +function alg_cache(prob::MatrixDEProblem, alg::DOAlgorithm{TruncatedSVD}, u, dt) + @unpack alg_params = alg + F = (u,t,dt) -> alg_params.F(prob.f, u, t, dt) + return TruncatedSVD_Cache(F) +end + +function step!(integrator::DLRIntegrator, alg::DOAlgorithm{TruncatedSVD}, dt) + @unpack u, t, iter, cache = integrator + truncated_svd_do_step!(u, cache, t, dt, alg.alg_params) + integrator.t += dt + integrator.iter += 1 +end + +function truncated_svd_do_step!(u, cache, t, dt, params) + @unpack F = cache + @unpack tol, rmin, rmax, SVD_alg, orth_alg = params + Y = u + dt*F(u, t, dt) + retract!(Y, u, SVD_alg, tol, rmin, rmax, orth_alg) +end + +function retract!(Y, u, alg::T, tol, rmin, rmax, orth_alg) where T <: Union{SVDFact, TSVD} + Y = round(Y, alg, tol=tol, rmin = rmin, rmax = rmax, alg_orthonormalize = orth_alg) + u.Z = Y.Z + u.U = Y.U +end + +function retract!(Y, u, alg::GradientDescent, tol, rmin, rmax, orth_alg) + gd_truncated_svd!(u, Y, alg, orth_alg) +end \ No newline at end of file diff --git a/src/integrators/greedy_integrator.jl b/src/integrators/greedy_integrator.jl index 8b0dae1..4524260 100644 --- a/src/integrators/greedy_integrator.jl +++ b/src/integrators/greedy_integrator.jl @@ -46,17 +46,17 @@ function alg_cache(prob::MatrixHybridProblem, alg::GreedyIntegrator, u, dt) return GreedyIntegrator_Cache(prob.y, X, XZ, nothing, nothing, ZIntegrator) end -function init(prob::MatrixDataProblem, alg::GreedyIntegrator, dt) - t0, tf = prob.tspan - @assert tf > t0 "Integration in reverse time direction is not supported" - u = deepcopy(prob.u0) - # initialize solution - sol = init_sol(dt, t0, tf, prob.u0) - # initialize cache - cache = alg_cache(prob, alg, u, dt) - sol.Y[1] = deepcopy(prob.u0) - return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -end +# function init(prob::MatrixDataProblem, alg::GreedyIntegrator, dt) +# t0, tf = prob.tspan +# @assert tf > t0 "Integration in reverse time direction is not supported" +# u = deepcopy(prob.u0) +# # initialize solution +# sol = init_sol(dt, t0, tf, prob.u0) +# # initialize cache +# cache = alg_cache(prob, alg, u, dt) +# sol.Y[1] = deepcopy(prob.u0) +# return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) +# end function init(prob::MatrixHybridProblem, alg::GreedyIntegrator, dt) t0, tf = prob.tspan diff --git a/src/integrators/projector_splitting.jl b/src/integrators/projector_splitting.jl index 7459519..adc8cba 100644 --- a/src/integrators/projector_splitting.jl +++ b/src/integrators/projector_splitting.jl @@ -104,15 +104,15 @@ function alg_cache(prob::MatrixDataProblem, alg::ProjectorSplitting, u, dt; t0 = prob.y, ycurr, yprev, Δy) end -function init(prob::AbstractDLRProblem, alg::ProjectorSplitting, dt) - t0, tf = prob.tspan - @assert tf > t0 "Integration in reverse time direction is not supported" - u = deepcopy(prob.u0) - sol = init_sol(dt, t0, tf, prob.u0) - cache = alg_cache(prob, alg, u, dt, t0 = t0) - sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object - return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -end +# function init(prob::AbstractDLRProblem, alg::ProjectorSplitting, dt) +# t0, tf = prob.tspan +# @assert tf > t0 "Integration in reverse time direction is not supported" +# u = deepcopy(prob.u0) +# sol = init_sol(dt, t0, tf, prob.u0) +# cache = alg_cache(prob, alg, u, dt, t0 = t0) +# sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object +# return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) +# end function primal_LT_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem}) @unpack y, ycurr, yprev, Δy = cache diff --git a/src/integrators/rank_adaptive_unconventional.jl b/src/integrators/rank_adaptive_unconventional.jl index 890f72c..ecaf6d9 100644 --- a/src/integrators/rank_adaptive_unconventional.jl +++ b/src/integrators/rank_adaptive_unconventional.jl @@ -47,7 +47,7 @@ end function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) # allocate memory for frequently accessed arrays tspan = (t0, t0+dt) - r = rank(u) + r = LowRankArithmetic.rank(u) n, m = size(u) US = u.U*u.S VS = u.V*u.S' @@ -91,17 +91,17 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit nothing, nothing, nothing, nothing) end -function init(prob::AbstractDLRProblem, alg::RankAdaptiveUnconventionalAlgorithm, dt) - t0, tf = prob.tspan - @assert tf > t0 "Integration in reverse time direction is not supported" - u = deepcopy(prob.u0) - # initialize solution - sol = init_sol(dt, t0, tf, prob.u0) - # initialize cache - cache = alg_cache(prob, alg, u, dt; t0 = t0) - sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object - return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -end +# function init(prob::AbstractDLRProblem, alg::RankAdaptiveUnconventionalAlgorithm, dt) +# t0, tf = prob.tspan +# @assert tf > t0 "Integration in reverse time direction is not supported" +# u = deepcopy(prob.u0) +# # initialize solution +# sol = init_sol(dt, t0, tf, prob.u0) +# # initialize cache +# cache = alg_cache(prob, alg, u, dt; t0 = t0) +# sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object +# return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) +# end function alg_cache(prob::MatrixDataProblem, alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) # creates caches for frequently used arrays by performing the first time step @@ -132,7 +132,7 @@ end function alg_recache(cache::RankAdaptiveUnconventionalAlgorithm_Cache, alg::RankAdaptiveUnconventionalAlgorithm, u, t) @unpack SProblem, KProblem, LProblem, tol, r_max, y, yprev, ycurr, Δy = cache - r = rank(u) + r = LowRankArithmetic.rank(u) n, m = size(u) US = zeros(n,r) @@ -220,7 +220,7 @@ function rankadaptive_unconventional_step!(u, cache, t, dt) step!(SIntegrator, dt, true) U, S, V = svd(SIntegrator.u) - r_new = min(r_max, LowRankArithmetic.truncate_to_tolerance(S, tol, rel=true)) + r_new = min(r_max, LowRankArithmetic.truncate_to_tolerance(S, tol))#, rel=true)) if r_new == r mul!(u.U,Uhat,U[:,1:r_new]) u.S .= Matrix(Diagonal(S[1:r_new])) diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index f3ff9a4..0f974a3 100644 --- a/src/integrators/unconventional.jl +++ b/src/integrators/unconventional.jl @@ -26,27 +26,13 @@ end function UnconventionalAlgorithm(deim; S_rhs = nothing, L_rhs = nothing, K_rhs = nothing, S_kwargs = Dict(), L_kwargs = Dict(), K_kwargs = Dict(), S_alg = Tsit5(), L_alg = Tsit5(), K_alg = Tsit5()) - params = UnconventionalAlgorithm_Params(S_rhs, L_rhs, K_rhs, + params = UnconventionalAlgorithm_Params(S_rhs, L_rhs, K_rhs, S_kwargs, L_kwargs, K_kwargs, S_alg, L_alg, K_alg, deim) return UnconventionalAlgorithm(params) end -struct SparseInterpolation - selection_alg - tol - rmin - rmax - init_range - init_corange -end - -function SparseInterpolation(selection_alg, init_range, init_corange; - rmin = 1, rmax = min(size(init_range,1), size(init_corange,1)), tol = eps(Float64)) - return SparseInterpolation(selection_alg, tol, rmin, rmax, init_range, init_corange) -end - struct UnconventionalAlgorithm_Cache{uType, SIntegratorType,LIntegratorType,KIntegratorType, yFunc, yType, dType} <: AbstractDLRAlgorithm_Cache @@ -68,12 +54,20 @@ end struct UnconventionalDEIM_Cache params::SparseInterpolation + u_prev::SVDLikeRepresentation Π::SparseFunctionInterpolator Π_K::SparseFunctionInterpolator Π_L::SparseFunctionInterpolator Π_S::SparseFunctionInterpolator end +rank_DEIM(cache::UnconventionalAlgorithm_Cache) = rank_DEIM(cache.interpolation_cache) +rank_DEIM(::Nothing) = 0 +rank_DEIM(cache::UnconventionalDEIM_Cache) = rank(cache.Π)[1] +interpolation_indices(cache::UnconventionalAlgorithm_Cache) = interpolation_indices(cache.interpolation_cache) +interpolation_indices(::Nothing) = (Int[],Int[]) +interpolation_indices(cache::UnconventionalDEIM_Cache) = interpolation_indices(cache.Π) + function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) # allocate memory for frequently accessed arrays tspan = (t0,t0+dt) @@ -124,7 +118,7 @@ end function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) where - {fType <: ComponentFunction} + {fType <: ComponentFunction, uType, tType} # allocate memory for frequently accessed arrays tspan = (t0,t0+dt) @@ -135,9 +129,9 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, row_idcs = index_selection(init_range, selection_alg) col_idcs = index_selection(init_corange, selection_alg) - Π_corange = DEIMInterpolator(col_idcs, init_range/init_range[col_idcs,:]) - Π_range = DEIMInterpolator(row_idcs, init_corange/init_corange[row_idcs,:]) - Π = SparseFunctionInterpolator(F, SparseMatrixInterpolator(Π_range, Π_corange)) + Π_corange = DEIMInterpolator(col_idcs, init_corange/init_corange[col_idcs,:]) + Π_range = DEIMInterpolator(row_idcs, init_range/init_range[row_idcs,:]) + Π = SparseFunctionInterpolator(prob.f, SparseMatrixInterpolator(Π_range, Π_corange)) US = u.U*u.S QRK = qr(US) @@ -156,8 +150,8 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, K_rhs = alg.alg_params.K_rhs end Π_K_corange = DEIMInterpolator(col_idcs, - V0'*Π.interpolator.corange.weights)' - Π_K = SparseFunctionInterpolator(prob.F, Π_K_corange) + u.V'*Π.interpolator.corange.weights)' + Π_K = SparseFunctionInterpolator(prob.f, Π_K_corange) p_K = (Π_K, u.V, ()) KProblem = ODEProblem(K_rhs, US, tspan, p_K) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) @@ -171,10 +165,10 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, L_rhs = alg.alg_params.L_rhs end Π_L_range = DEIMInterpolator(row_idcs, - U0'*Π.interpolator.range.weights) - Π_K = SparseFunctionInterpolator(prob.F, Π_L_range) - p_K = (Π_L, u.U, ()) - LProblem = ODEProblem(L_rhs, VS, tspan, u.U) + u.U'*Π.interpolator.range.weights) + Π_L = SparseFunctionInterpolator(prob.f, Π_L_range) + p_L = (Π_L, u.U, ()) + LProblem = ODEProblem(L_rhs, VS, tspan, p_L) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) if isnothing(alg.alg_params.S_rhs) @@ -185,16 +179,17 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, else S_rhs = alg.alg_params.S_rhs end - Π_S_mat = SparseMatrixInterpolator(row_dics, col_idcs, + Π_S_mat = SparseMatrixInterpolator(row_idcs, col_idcs, u.U'*Π.interpolator.range.weights, u.V'*Π.interpolator.corange.weights) - Π_S = SparseFunctionInterpolator(prob.F, Π_S_mat) + Π_S = SparseFunctionInterpolator(prob.f, Π_S_mat) p_S = (Π_S, u.U, u.V, ()) SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, p_S) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) - interpolation_cache = UnconventionalDEIM_Cache(alg.alg_params.interpolation, Π, Π_K, Π_L, Π_S) + interpolation_cache = UnconventionalDEIM_Cache(alg.alg_params.interpolation, deepcopy(u), + Π, Π_K, Π_L, Π_S) return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, interpolation_cache) @@ -224,17 +219,17 @@ function alg_cache(prob::MatrixDataProblem, alg::UnconventionalAlgorithm, u, dt; prob.y, ycurr, yprev, Δy, nothing) end -function init(prob::AbstractDLRProblem, alg::UnconventionalAlgorithm, dt) - t0, tf = prob.tspan - @assert tf > t0 "Integration in reverse time direction is not supported" - u = deepcopy(prob.u0) - # initialize solution - sol = init_sol(dt, t0, tf, prob.u0) - # initialize cache - cache = alg_cache(prob, alg, u, dt, t0 = t0) - sol.Y[1] = deepcopy(prob.u0) - return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -end +# function init(prob::AbstractDLRProblem, alg::UnconventionalAlgorithm, dt) +# t0, tf = prob.tspan +# @assert tf > t0 "Integration in reverse time direction is not supported" +# u = deepcopy(prob.u0) +# # initialize solution +# sol = init_sol(dt, t0, tf, prob.u0) +# # initialize cache +# cache = alg_cache(prob, alg, u, dt, t0 = t0) +# sol.Y[1] = deepcopy(prob.u0) +# return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) +# end function unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem}) @unpack y, ycurr, yprev, Δy = cache @@ -248,6 +243,10 @@ function unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem}) unconventional_step!(u, cache, t, dt) end +function unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem{<:ComponentFunction}}) + unconventional_deim_step!(u, cache, t, dt) +end + function unconventional_step!(u, cache, t, dt) @unpack US, VS, M, N, QRK, QRL, KIntegrator, LIntegrator, SIntegrator = cache @@ -279,7 +278,13 @@ function unconventional_deim_step!(u, cache, t, dt) @unpack US, VS, M, N, QRK, QRL, KIntegrator, LIntegrator, SIntegrator, interpolation_cache = cache - @unpack Π, Π_L, Π_K, Π_S = interpolation_cache + @unpack params, u_prev, Π, Π_L, Π_K, Π_S = interpolation_cache + + if params.update_scheme == :avg_flow + u_prev.U = copy(u.U) + u_prev.S = copy(u.S) + u_prev.V = copy(u.V) + end # K step mul!(US, u.U, u.S) @@ -300,29 +305,32 @@ function unconventional_deim_step!(u, cache, t, dt) u.U .= Matrix(QRK.Q) # update core interpolator - mul!(Π_S.range.weights, u.U', Π.interpolator.range.weights) - mul!(Π_S.corange.weights, u.V', Π.interpolator.corange.weights) + mul!(Π_S.interpolator.range.weights, u.U', Π.interpolator.range.weights) + mul!(Π_S.interpolator.corange.weights, u.V', Π.interpolator.corange.weights) # integration set_u!(SIntegrator, M*u.S*N') step!(SIntegrator, dt, true) u.S .= SIntegrator.u - update_interpolation!(interpolation_cache, u) + update_interpolation!(interpolation_cache, u, t, dt) + println(rank(interpolation_cache.Π)) end -function update_interpolation!(interpolation_cache, u) +# without rank adaptation +#= +function update_interpolation!(interpolation_cache, u, t) @unpack params, Π, Π_K, Π_L, Π_S = interpolation_cache - @unpack selection_method, tol, rmin, rmax = params + @unpack selection_alg, tol, rmin, rmax = params - eval!(Π_L, u, (), SIntegrator.t) # rows - eval!(Π_K, u, (), SIntegrator.t) # columns + eval!(Π_L, u, (), t) # rows + eval!(Π_K, u, (), t) # columns VF = truncated_svd(Π_L.cache.rows, tol=tol, rmin=rmin, rmax=rmax).V # corange from rows UF = truncated_svd(Π_K.cache.cols, tol=tol, rmin=rmin, rmax=rmax).U # range from cols - + # find interpolation indices - row_idcs = index_selection(UF, selection_method) - col_idcs = index_selection(VF, selection_method) + row_idcs = index_selection(range_svd.U, selection_alg) + col_idcs = index_selection(corange_svd.V, selection_alg) # find interpolation weights @views range_weights = UF/UF[row_idcs,:] @@ -332,7 +340,136 @@ function update_interpolation!(interpolation_cache, u) range = DEIMInterpolator(row_idcs, range_weights) projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) corange = DEIMInterpolator(col_idcs, corange_weights) - projected_corange = DEIMInterpolator(row_idcs, u.V'*corange_weights) + projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) + + # update function interpolators + update_interpolator!(Π_L, projected_range) + update_interpolator!(Π_K, projected_corange') + update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) + update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) +end + +# alternative but max rank limited by 2*rank(u) +# ... not great +function update_interpolation!(interpolation_cache, u_old, u, dt) + @unpack params, Π, Π_K, Π_L, Π_S = interpolation_cache + @unpack selection_alg, tol, rmin, rmax = params + + dF = u - u_old + dF_svd = svd(dF) + #println(typeof(dF_svd.S)) + cutoff = findfirst(s -> s < tol/dt, diag(dF_svd.S) ) + if isnothing(cutoff) + # compute the full function and take derivative + cutoff= size(dF_svd.S,1) + else + println(dF_svd.S[cutoff]) + cutoff = max(rmin, min(cutoff, rmax)) + end + + if cutoff != rank(Π_K) + println("new DEIM rank: $cutoff") + end + + # find interpolation indices + @views row_idcs = index_selection(dF_svd.U[:,1:cutoff], selection_alg) + @views col_idcs = index_selection(dF_svd.V[:,1:cutoff], selection_alg) + + # find interpolation weights + @views range_weights = dF_svd.U/dF_svd.U[row_idcs,:] + @views corange_weights = dF_svd.V/dF_svd.V[col_idcs,:] + + # new interpolators + range = DEIMInterpolator(row_idcs, range_weights) + projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) + corange = DEIMInterpolator(col_idcs, corange_weights) + projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) + + # update function interpolators + update_interpolator!(Π_L, projected_range) + update_interpolator!(Π_K, projected_corange') + update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) + update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) +end +=# + +function select_idcs(F_rows::AbstractMatrix, F_cols::AbstractMatrix, selection_alg::IndexSelectionAlgorithm) + range = svd(F_cols) + corange = svd(F_rows) + S = max.(range.S, corange.S) + row_idcs = index_selection(range.U, S, selection_alg) + col_idcs = index_selection(Matrix(corange.V), S, selection_alg) + return range.U, corange.V, row_idcs, col_idcs +end + +function select_idcs(F::SVDLikeRepresentation, selection_alg::IndexSelectionAlgorithm) + F_svd = svd(F) + row_idcs = index_selection(F_svd.U, diag(F_svd.S), selection_alg) + col_idcs = index_selection(F_svd.V, diag(F_svd.S), selection_alg) + UF = F_svd.U[:, 1:length(row_idcs)] + VF = F_trunc.V[:, 1:length(col_idcs)] + return UF, VF, row_idcs, col_idcs +end + +function sparse_selection(interpolation_cache, u, t, dt) + @unpack params = interpolation_cache + if params.update_scheme == :avg_flow + @unpack u_prev = interpolation_cache + @unpack selection_alg = params + return select_idcs(u - u_prev, selection_alg) + elseif params.update_scheme == :last_iterate + @unpack Π_K, Π_L = interpolation_cache + @unpack selection_alg = params + eval!(Π_L, u, (), t+dt) + eval!(Π_K, u, (), t+dt) + UF, VF, row_idcs, col_idcs = select_idcs(Π_L.cache.rows, Π_K.cache.cols, selection_alg) + return UF, VF, row_idcs, col_idcs + else + error("Update scheme $(params.update_scheme) not supported") + end + + # eval index selection + +end + +function compute_weights(range, corange, row_idcs, col_idcs, u, t, dt, interpolation_cache) + @unpack params = interpolation_cache + if size(range, 2) != length(row_idcs) + @unpack Π = interpolation_cache + cols = Matrix{eltype(Π.cache.cols)}(undef, Π.F.n_rows, length(col_idcs)) + eval_cols!(cols, Π, u, (), t+dt, col_idcs) + UF = svd(cols).U + @views range_weights = UF/UF[row_idcs,:] + else + @views range_weights = range/range[row_idcs,:] + end + + if size(corange, 2) != length(col_idcs) + @unpack Π = interpolation_cache + rows = Matrix{eltype(Π.cache.cols)}(undef, length(row_idcs), Π.F.n_cols) + eval_rows!(rows, Π, u, (), t+dt, row_idcs) + VF = svd(rows).V + @views corange_weights = VF/VF[col_idcs,:] + else + @views corange_weights = corange/corange[col_idcs,:] + end + return range_weights, corange_weights +end + +function update_interpolation!(interpolation_cache, u, t, dt) + @unpack Π, Π_L, Π_K, Π_S, params = interpolation_cache + @unpack selection_alg = params + + # eval range/corange spans + UF, VF, row_idcs, col_idcs = sparse_selection(interpolation_cache, u, t, dt) + range_weights, corange_weights = compute_weights(UF, VF, row_idcs, col_idcs, u, t, dt, interpolation_cache) + + + # new interpolators + range = DEIMInterpolator(row_idcs, range_weights) + projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) + corange = DEIMInterpolator(col_idcs, corange_weights) + projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) # update function interpolators update_interpolator!(Π_L, projected_range) @@ -346,4 +483,4 @@ function step!(integrator::DLRIntegrator, ::UnconventionalAlgorithm, dt) unconventional_step!(u, cache, t, dt, probType) integrator.t += dt integrator.iter += 1 -end +end \ No newline at end of file diff --git a/src/primitives.jl b/src/primitives.jl index ad54beb..66ada46 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -46,6 +46,9 @@ end mutable struct DLRSolution{solType, tType} <: AbstractDLRSolution Y::Vector{solType} t::Vector{tType} + r::Vector{Int} + r_DEIM::Vector{Int} + idcs_DEIM::Vector{Tuple{Vector{Int}, Vector{Int}}} end """ @@ -65,16 +68,23 @@ end """ solves the given problem with the specified algorithm and step size """ -function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt) - println("Initialize integrator:") - integrator = init(prob, alg, dt) - println("Initialization complete. Start Integration.") - T = prob.tspan[2] - prob.tspan[1] - while (prob.tspan[2]-integrator.t)/T > 1e-8 +function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt; save_increment::Int=1) + println("Initialize integrator ...") + integrator, t_int, save_idcs = init(prob, alg, dt, save_increment) + println("... initialization complete. Start Integration ...") + dt = step(t_int) + disp_digits = abs(round(Int, log10(save_increment*dt))) + k = 2 + #while (prob.tspan[2]-integrator.t)/T > 1e-8 + while integrator.iter < length(t_int) - 1 step!(integrator, alg, dt) - update_sol!(integrator) - println("t = $(integrator.t)") + if integrator.iter + 1 == save_idcs[k] + update_sol!(integrator, k) + k += 1 + println("t = $(round(integrator.t,sigdigits=disp_digits))") + end end + println("... integration complete.") return integrator.sol end function solve(prob::MatrixDataProblem, alg::AbstractDLRAlgorithm) @@ -82,26 +92,65 @@ function solve(prob::MatrixDataProblem, alg::AbstractDLRAlgorithm) return solve(prob, alg, 1) end -function update_sol!(integrator::AbstractDLRIntegrator) - if integrator.iter <= length(integrator.sol.Y) - 1 - integrator.sol.Y[integrator.iter + 1] = deepcopy(integrator.u) - integrator.sol.t[integrator.iter + 1] = integrator.t - else - push!(integrator.sol.Y, deepcopy(integrator.u)) - push!(integrator.sol.t, integrator.t) - end +function init(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt, save_increment::Int) + t0, tf = prob.tspan + @assert tf > t0 "Integration in reverse time direction is not supported" + u = deepcopy(prob.u0) + # initialize solution + sol, t_int, save_idcs = init_sol(dt, t0, tf, prob.u0, save_increment) + # initialize cache + cache = alg_cache(prob, alg, u, dt, t0 = t0) + sol.Y[1] = deepcopy(prob.u0) + return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0), t_int, save_idcs end -function init_sol(dt, t0, tf, u0) +function init_sol(dt, t0, tf, u0, save_increment) n = floor(Int,(tf-t0)/dt) + 1 - # compute more sensible dt # rest will be handled via interpolation/save_at dt = (tf-t0)/(n-1) + t_int = t0:dt:tf + save_idcs = collect(1:save_increment:length(t_int)) + if save_idcs[end] != length(t_int) + push!(save_idcs, length(t_int)) + end + n_eff = length(save_idcs) # initialize & return solution object - return DLRSolution(Vector{typeof(u0)}(undef, n), collect(range(t0, tf, length=n))) + Y = Vector{typeof(u0)}(undef, n_eff) + r = Vector{Int}(undef, n_eff) + r_deim = Vector{Int}(undef, n_eff) + t = collect(range(t0, tf, length=n_eff)) + interpolation_idcs = Vector{Tuple{Vector{Int},Vector{Int}}}(undef, n_eff) + return DLRSolution(Y, t, r, r_deim, interpolation_idcs), t_int, save_idcs end function init_sol(dt::Int, t0, tf, u0) # initialize solution object - steps = t0:dt:tf - return DLRSolution(Vector{typeof(u0)}(undef, length(steps)), collect(steps)) + n = floor(Int,(tf-t0)/dt) + 1 + dt = (tf-t0)/(n-1) + t_int = t0:dt:tf + save_idcs = 1:save_increment:length(t_int) + if save_idcs[end] != length(t_int) + push!(save_idcs, length(t_int)) + end + n_eff = length(save_idcs) + Y = Vector{typeof(u0)}(undef, n_eff) + r = Vector{Int}(undef, n_eff) + r_deim = Vector{Int}(undef, n_eff) + interpolation_idcs = Vector{Tuple{Vector{Int},Vector{Int}}}(undef, n_eff) + return DLRSolution(Y, steps, r, r_deim, interpolation_idcs), t_int, save_idcs +end + +function update_sol!(integrator::AbstractDLRIntegrator, idx) + #if integrator.iter <= length(integrator.sol.Y) - 1 + integrator.sol.Y[idx] = deepcopy(integrator.u) + integrator.sol.t[idx] = integrator.t + integrator.sol.r[idx] = LowRankArithmetic.rank(integrator.u) + integrator.sol.r_DEIM[idx] = rank_DEIM(integrator.cache) + integrator.sol.idcs_DEIM[idx] = interpolation_indices(integrator.cache) + #else + # push!(integrator.sol.Y, deepcopy(integrator.u)) + # push!(integrator.sol.t, integrator.t) + # push!(integrator.sol.r[integrator.iter + 1], LowRankArithmetic.rank(integrator.u)) + # push!(integrator.sol.r_DEIM[integrator.iter + 1], rank_DEIM(integrator.cache)) + # push!(integrator.sol.r_DEIM[integrator.iter + 1], indices(integrator.cache)) + # end end \ No newline at end of file diff --git a/src/sparse_interpolation.jl b/src/sparse_interpolation.jl index bf91332..beee73c 100644 --- a/src/sparse_interpolation.jl +++ b/src/sparse_interpolation.jl @@ -1,7 +1,6 @@ -using DocStringExtensions, UnPack, LinearAlgebra - import LinearAlgebra.adjoint import Base.size + """ $(TYPEDEF) @@ -27,6 +26,14 @@ struct InterpolatorCache{T} <: AbstractInterpolatorCache elements::Matrix{T} end +function resize_rows(IC::InterpolatorCache{T}, n_rows::Int) where T + return InterpolatorCache(Matrix{T}(undef, n_rows, size(IC.rows,2)), IC.cols, IC.rows) +end + +function resize_cols(IC::InterpolatorCache{T}, n_cols::Int) where T + return InterpolatorCache(IC.rows, Matrix{T}(undef, size(IC.cols,2), n_cols), IC.rows) +end + """ $(TYPEDEF) @@ -168,6 +175,16 @@ function SparseFunctionInterpolator(F, interpolator::SparseMatrixInterpolator; Matrix{output}(undef, r_rows, r_cols)) return SparseFunctionInterpolator(F, interpolator, cache) end +""" + $(TYPEDSIGNATURES) + + helper function querying the interpolation indices. +""" +interpolation_indices(Π::DEIMInterpolator) = Π.interpolation_idcs +interpolation_indices(Π::AdjointDEIMInterpolator) = Π.parent.interpolation_idcs +interpolation_indices(Π::SparseMatrixInterpolator) = (interpolation_indices(Π.range), interpolation_indices(Π.corange)) +interpolation_indices(Π::SparseFunctionInterpolator) = interpolation_indices(Π.interpolator) + """ $(TYPEDSIGNATURES) @@ -178,7 +195,7 @@ function update_interpolator!(Π::SparseFunctionInterpolator, interpolator::DEIM Π.interpolator = interpolator r_rows = rank(interpolator) T = eltype(Π.cache.cols) - Π.cache = InterpolationCache(Matrix{T}(undef, r_rows, F.n_cols), + Π.cache = InterpolatorCache(Matrix{T}(undef, r_rows, Π.F.n_cols), Matrix{T}(undef, 0, 0), Matrix{T}(undef, 0, 0)) end @@ -187,20 +204,20 @@ function update_interpolator!(Π::SparseFunctionInterpolator, interpolator::Adjo Π.interpolator = interpolator r_cols = rank(interpolator) T = eltype(Π.cache.cols) - Π.cache = InterpolationCache(Matrix{T}(undef, 0, 0), - Matrix{T}(undef, F.n_rows, r_cols), + Π.cache = InterpolatorCache(Matrix{T}(undef, 0, 0), + Matrix{T}(undef, Π.F.n_rows, r_cols), Matrix{T}(undef, 0, 0)) end function update_interpolator!(Π::SparseFunctionInterpolator, interpolator::SparseMatrixInterpolator) Π.interpolator = interpolator r_rows, r_cols = rank(interpolator.range), rank(interpolator.corange) - dim_range = size(interpolator.range.interpolator, 1) + dim_range = size(interpolator.range.weights, 1) dim_corange = size(interpolator.corange.weights, 1) T = eltype(Π.cache.cols) - Π.cache = InterpolationCache(Matrix{T}(undef, r_rows, dim_corange), - Matrix{T}(undef, dim_range, r_cols), - Matrix{T}(undef, r_rows, r_cols)) + Π.cache = InterpolatorCache(Matrix{T}(undef, r_rows, dim_corange), + Matrix{T}(undef, dim_range, r_cols), + Matrix{T}(undef, r_rows, r_cols)) end """ @@ -226,6 +243,14 @@ function eval!(Π::SparseFunctionInterpolator{IP, fType, T}, x, p, t) where {IP Π.F.elements!(cache.elements, x, p, t, range.interpolation_idcs, corange.interpolation_idcs) end +function eval_rows!(dx, Π::SparseFunctionInterpolator, x, p, t, idcs) + Π.F.rows!(dx, x, p, t, idcs) +end + +function eval_cols!(dx, Π::SparseFunctionInterpolator, x, p, t, idcs) + Π.F.cols!(dx, x, p, t, idcs) +end + """ $(TYPEDSIGNATURES) @@ -290,21 +315,55 @@ abstract type IndexSelectionAlgorithm end DEIM index selection algorithm. """ -struct DEIM end +struct DEIM <: IndexSelectionAlgorithm + tol::Float64 + elasticity::Float64 + rmin::Int + rmax::Int +end + +function DEIM( ;tol = eps(Float64), rmin = 1, rmax = 2^62, elasticity = 0.1) + return DEIM(tol, elasticity, rmin, rmax) +end """ $(TYPEDEF) QDEIM index selection algorithm. """ -struct QDEIM end +struct QDEIM <: IndexSelectionAlgorithm + tol::Float64 + elasticity::Float64 + rmin::Int + rmax::Int +end + +function QDEIM( ;tol = eps(Float64), rmin = 1, rmax = 2^62, elasticity = 0.1) + return QDEIM(tol, elasticity, rmin, rmax) +end + +""" + $(TYPEDEF) + + LDEIM index selection algorithm. Selects `n` indices. +""" +struct LDEIM <: IndexSelectionAlgorithm + tol::Float64 + elasticity::Float64 + rmin::Int + rmax::Int +end + +function LDEIM( ;tol = eps(Float64), rmin = 1, rmax = 2^62, elasticity = 0.1) + return LDEIM(tol, elasticity, rmin, rmax) +end """ $(TYPEDSIGNATURES) returns interpolation indices for sparse interpolation. Supports DEIM and QDEIM index selection. """ -function index_selection(U::Matrix,::DEIM) +function index_selection(U,::DEIM) m = size(U, 2) indices = Vector{Int}(undef, m) @inbounds @views begin @@ -325,11 +384,94 @@ function index_selection(U::Matrix,::DEIM) return indices end -function index_selection(U::Matrix,::QDEIM) +function index_selection(U,::QDEIM) QR = qr(U', ColumnNorm()) return QR.p[1:size(U,2)] end +function index_selection(U::Matrix, r::Int, ::LDEIM) + n, k = size(U) + if k > r + return index_selection(view(U, :, 1:r), DEIM()) + else + Ψ = copy(U) + idcs = zeros(Int, r) + @views for i in 1:k-1 + idcs[i] = argmax(Iterators.map(abs, Ψ[:,i])) + corr = Ψ[idcs[1:i],1:i] \ Ψ[idcs[1:i], i+1] + Ψ[:,i+1] -= Ψ[:,1:i] * corr + end + @views idcs[k] = argmax(Iterators.map(abs, Ψ[:,k])) + @views l = [i in idcs ? -1.0 : norm(Ψ[i,:]) for i in 1:n] + idcs[k+1:end] = partialsortperm(l, 1:r-k, by=-) + return idcs + end +end + +function index_selection(U::Matrix, ::LDEIM) + return index_selection(U, DEIM()) +end + +function index_selection(U::Matrix, S::Vector, alg::IndexSelectionAlgorithm) + @assert size(U,2) == length(S) "Number of columns provided must coincide with number of singular values." + + cutoff = length(S) + for σ in reverse(S) + if σ > alg.tol + break + end + cutoff -= 1 + end + r = max(alg.rmin, min(alg.rmax, cutoff)) + return index_selection(U[:,1:r], alg) +end + +function index_selection(U::Matrix, S::Vector, alg::LDEIM) + @assert size(U,2) == length(S) "Number of columns provided must coincide with number of singular values." + + if S[end] > alg.tol + r = min(alg.rmax, length(S)+1) + return index_selection(U, r, alg) + else + #cutoff = findfirst(x -> x < alg.tol*alg.elasticity, S) + cutoff = findlast(x -> x > alg.tol*alg.elasticity, S) + if isnothing(cutoff) + cutoff = 1 + end + #cutoff = length(S) + # for σ in reverse(S) + # if σ > alg.tol*alg.elasticity + # break + # end + # cutoff -= 1 + # end + return index_selection(U, cutoff, alg) + end +end + + +""" + $(TYPEDEF) + + type for carrying all necessary information for continuous sparse interpolation + to be applied for dynamical low rank approximation. +""" +struct SparseInterpolation + selection_alg + update_scheme + tol + rmin + rmax + init_range + init_corange +end + +function SparseInterpolation(selection_alg, init_range, init_corange; update_scheme = :last_iterate, + rmin = 1, rmax = min(size(init_range,1), size(init_corange,1)), tol = eps(Float64)) + return SparseInterpolation(selection_alg, update_scheme, tol, rmin, rmax, init_range, init_corange) +end + + # a few lil tests #= function rows!(dx,x,p,t,idcs) From 5a924eb445f103619a9f235b3767e8515779f192 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Sun, 12 Mar 2023 21:11:25 -0400 Subject: [PATCH 07/35] ignore data file formats and figures --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index fa297c4..1d556a5 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,6 @@ examples/.ipynb_checkpoints *.mp4 products/ .DS_STORE +*.jld2 +*.pdf + From a83853fffc521569d5883d0dc4d5b6ef72e4ff85 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Sun, 12 Mar 2023 22:25:35 -0400 Subject: [PATCH 08/35] update testing environment and add data-informed test --- test/runtests.jl | 3 ++- test/sparse_interpolation.jl | 0 2 files changed, 2 insertions(+), 1 deletion(-) create mode 100644 test/sparse_interpolation.jl diff --git a/test/runtests.jl b/test/runtests.jl index 0475d8f..eeb9b1c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,4 +2,5 @@ using LowRankIntegrators using Test include("data_agnostic_approximation.jl") -include("data_driven_approximation.jl") \ No newline at end of file +include("data_driven_approximation.jl") +include("data_informed_approximation.jl") \ No newline at end of file diff --git a/test/sparse_interpolation.jl b/test/sparse_interpolation.jl new file mode 100644 index 0000000..e69de29 From 34f05c8e2e2380511be5c4f2a95b5ebb22439bf2 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Sun, 12 Mar 2023 23:15:57 -0400 Subject: [PATCH 09/35] update test env --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index f63ddf1..b36202c 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LowRankIntegrators = "a4d2327a-b40c-407e-b18c-dbb78b346b8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 56b59cbd11d3144a1f35433cc3c5f548c138d877 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Sun, 12 Mar 2023 23:16:29 -0400 Subject: [PATCH 10/35] fix bug in greedy integrator for data driven problem --- src/integrators/greedy_integrator.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/integrators/greedy_integrator.jl b/src/integrators/greedy_integrator.jl index 4524260..5dbd2ad 100644 --- a/src/integrators/greedy_integrator.jl +++ b/src/integrators/greedy_integrator.jl @@ -21,7 +21,7 @@ function GreedyIntegrator(; Z_alg = Tsit5(), Z_kwargs = Dict{Symbol,Any}()) return GreedyIntegrator(params) end -function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::TwoFactorRepresentation, dt) +function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::TwoFactorRepresentation, dt; t0 = prob.tspan[1]) X = zeros(size(u)) r = rank(u) n = size(X,1) @@ -29,7 +29,7 @@ function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::TwoFactorR return GreedyIntegrator_Cache(prob.y, X, XZ, nothing, nothing, nothing) end -function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::SVDLikeRepresentation, dt) +function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::SVDLikeRepresentation, dt; t0 = prob.tspan[1]) X = zeros(size(u)) r = rank(u) n, m = size(X) @@ -38,7 +38,7 @@ function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::SVDLikeRep return GreedyIntegrator_Cache(prob.y, X, nothing, XV, XU, nothing) end -function alg_cache(prob::MatrixHybridProblem, alg::GreedyIntegrator, u, dt) +function alg_cache(prob::MatrixHybridProblem, alg::GreedyIntegrator, u, dt; t0 = prob.tspan[1]) XZ = zeros(size(u,1), rank(u)) X = zeros(size(u)...) ZProblem = ODEProblem(prob.f, u.Z, prob.tspan, u.U) From ac6afca50ca143a35366ebfc8f2f2e75f38ffa8f Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 00:27:40 -0400 Subject: [PATCH 11/35] fix initialization change for data driven problems --- src/integrators/greedy_integrator.jl | 4 +- src/primitives.jl | 20 ++++----- src/sparse_interpolation.jl | 61 +--------------------------- 3 files changed, 13 insertions(+), 72 deletions(-) diff --git a/src/integrators/greedy_integrator.jl b/src/integrators/greedy_integrator.jl index 5dbd2ad..8ff2a5c 100644 --- a/src/integrators/greedy_integrator.jl +++ b/src/integrators/greedy_integrator.jl @@ -1,5 +1,5 @@ # greedy approach to fit -struct GreedyIntegrator_Cache +struct GreedyIntegrator_Cache <: AbstractDLRAlgorithm_Cache Y X XZ @@ -31,7 +31,7 @@ end function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::SVDLikeRepresentation, dt; t0 = prob.tspan[1]) X = zeros(size(u)) - r = rank(u) + r = LowRankArithmetic.rank(u) n, m = size(X) XV = zeros(n,r) XU = zeros(m,r) diff --git a/src/primitives.jl b/src/primitives.jl index 66ada46..264a043 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -70,9 +70,8 @@ end """ function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt; save_increment::Int=1) println("Initialize integrator ...") - integrator, t_int, save_idcs = init(prob, alg, dt, save_increment) + integrator, dt, t_int, save_idcs = init(prob, alg, dt, save_increment) println("... initialization complete. Start Integration ...") - dt = step(t_int) disp_digits = abs(round(Int, log10(save_increment*dt))) k = 2 #while (prob.tspan[2]-integrator.t)/T > 1e-8 @@ -97,11 +96,11 @@ function init(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt, save_incr @assert tf > t0 "Integration in reverse time direction is not supported" u = deepcopy(prob.u0) # initialize solution - sol, t_int, save_idcs = init_sol(dt, t0, tf, prob.u0, save_increment) + sol, dt, t_int, save_idcs = init_sol(dt, t0, tf, prob.u0, save_increment) # initialize cache cache = alg_cache(prob, alg, u, dt, t0 = t0) sol.Y[1] = deepcopy(prob.u0) - return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0), t_int, save_idcs + return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0), dt, t_int, save_idcs end function init_sol(dt, t0, tf, u0, save_increment) @@ -119,15 +118,16 @@ function init_sol(dt, t0, tf, u0, save_increment) r_deim = Vector{Int}(undef, n_eff) t = collect(range(t0, tf, length=n_eff)) interpolation_idcs = Vector{Tuple{Vector{Int},Vector{Int}}}(undef, n_eff) - return DLRSolution(Y, t, r, r_deim, interpolation_idcs), t_int, save_idcs + return DLRSolution(Y, t, r, r_deim, interpolation_idcs), dt, t_int, save_idcs end -function init_sol(dt::Int, t0, tf, u0) +function init_sol(dt::Int, t0::Int, tf::Int, u0, save_increment::Int) # initialize solution object - n = floor(Int,(tf-t0)/dt) + 1 - dt = (tf-t0)/(n-1) t_int = t0:dt:tf - save_idcs = 1:save_increment:length(t_int) + if !(tf in t_int) + t_int = vcat(collect(t_int), tf) + end + save_idcs = collect(1:save_increment:length(t_int)) if save_idcs[end] != length(t_int) push!(save_idcs, length(t_int)) end @@ -136,7 +136,7 @@ function init_sol(dt::Int, t0, tf, u0) r = Vector{Int}(undef, n_eff) r_deim = Vector{Int}(undef, n_eff) interpolation_idcs = Vector{Tuple{Vector{Int},Vector{Int}}}(undef, n_eff) - return DLRSolution(Y, steps, r, r_deim, interpolation_idcs), t_int, save_idcs + return DLRSolution(Y, steps, r, r_deim, interpolation_idcs), dt, t_int, save_idcs end function update_sol!(integrator::AbstractDLRIntegrator, idx) diff --git a/src/sparse_interpolation.jl b/src/sparse_interpolation.jl index beee73c..417e587 100644 --- a/src/sparse_interpolation.jl +++ b/src/sparse_interpolation.jl @@ -469,63 +469,4 @@ end function SparseInterpolation(selection_alg, init_range, init_corange; update_scheme = :last_iterate, rmin = 1, rmax = min(size(init_range,1), size(init_corange,1)), tol = eps(Float64)) return SparseInterpolation(selection_alg, update_scheme, tol, rmin, rmax, init_range, init_corange) -end - - -# a few lil tests -#= -function rows!(dx,x,p,t,idcs) - @views dx .= x[idcs, :] -end -function cols!(dx,x,p,t,idcs) - @views dx .= x[:, idcs] -end -function elements!(dx,x,p,t,idcsA, idcsB) - @views dx .= x[idcsA, idcsB] -end - -F = ComponentFunction(rows!, cols!, elements!, 5, 8) -A = Matrix(reshape(1.0:5*8, (5, 8))) - -row_idcs = [1, 5, 2] -range_weights = zeros(eltype(A), 5, 3) -range_weights[row_idcs,1:3] .= Matrix{eltype(A)}(I,3,3) - -col_idcs = [1, 8, 3, 6] -corange_weights = zeros(eltype(A),8, 4) -corange_weights[col_idcs,1:4] .= Matrix{eltype(A)}(I,4,4) - -Π_range = DEIMInterpolator(row_idcs,range_weights) -Π_corange = DEIMInterpolator(col_idcs,corange_weights) -Π_mat = SparseMatrixInterpolator(Π_range, Π_corange) - -test = all(Π_range(A)[row_idcs,:] .== A[row_idcs,:]) -test = all(Π_corange'(A)[:,col_idcs] .== A[:,col_idcs]) - -F_range_int = SparseFunctionInterpolator(F, Π_range, output = eltype(A)) -F_corange_int = SparseFunctionInterpolator(F, Π_corange', output = eltype(A)) -F_int = SparseFunctionInterpolator(F, Π_mat, output = eltype(A)) - -eval!(F_range_int, A, (), 0.0) -test = all(F_range_int.cache.rows .== A[row_idcs,:]) -eval!(F_corange_int, A, (), 0.0) -test = all(F_corange_int.cache.cols .== A[:,col_idcs]) -eval!(F_int, A, (), 0.0) -test = all(F_int.cache.elements .== A[row_idcs,col_idcs]) - -dA = similar(A) -F_range_int(dA, A, (), 0.0) -test = all(dA[row_idcs, :] .== A[row_idcs, :]) -F_corange_int(dA, A, (), 0.0) -test = all(dA[:, col_idcs] .== A[:, col_idcs]) -F_int(dA, A, (), 0.0) -test = all(dA[row_idcs, col_idcs] .== A[row_idcs, col_idcs]) - -using BenchmarkTools -# non-allocating! -parms = () -t = 0.0 -@btime $F_range_int($dA, $A, $params, $t) -@btime $F_corange_int($dA, $A, $params, $t) -@btime $F_int($dA, $A, $params, $t) -=# \ No newline at end of file +end \ No newline at end of file From c8028138a0e4cead43ff7d0060ee734eb09d37cc Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 00:33:11 -0400 Subject: [PATCH 12/35] fix bug in solution time range --> endpoint now correctly stored --- src/primitives.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/primitives.jl b/src/primitives.jl index 264a043..33fa8e3 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -116,7 +116,7 @@ function init_sol(dt, t0, tf, u0, save_increment) Y = Vector{typeof(u0)}(undef, n_eff) r = Vector{Int}(undef, n_eff) r_deim = Vector{Int}(undef, n_eff) - t = collect(range(t0, tf, length=n_eff)) + t = collect(t_int[save_idcs]) interpolation_idcs = Vector{Tuple{Vector{Int},Vector{Int}}}(undef, n_eff) return DLRSolution(Y, t, r, r_deim, interpolation_idcs), dt, t_int, save_idcs end @@ -135,8 +135,9 @@ function init_sol(dt::Int, t0::Int, tf::Int, u0, save_increment::Int) Y = Vector{typeof(u0)}(undef, n_eff) r = Vector{Int}(undef, n_eff) r_deim = Vector{Int}(undef, n_eff) + t = collect(t_int[save_idcs]) interpolation_idcs = Vector{Tuple{Vector{Int},Vector{Int}}}(undef, n_eff) - return DLRSolution(Y, steps, r, r_deim, interpolation_idcs), dt, t_int, save_idcs + return DLRSolution(Y, t, r, r_deim, interpolation_idcs), dt, t_int, save_idcs end function update_sol!(integrator::AbstractDLRIntegrator, idx) From bdf3c395e0f2dffa693da4b3dcabad941aa4df70 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 00:37:00 -0400 Subject: [PATCH 13/35] rm print statement --- src/integrators/unconventional.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index 0f974a3..44cfc1d 100644 --- a/src/integrators/unconventional.jl +++ b/src/integrators/unconventional.jl @@ -314,7 +314,6 @@ function unconventional_deim_step!(u, cache, t, dt) u.S .= SIntegrator.u update_interpolation!(interpolation_cache, u, t, dt) - println(rank(interpolation_cache.Π)) end # without rank adaptation From d6818d6e34a913039fb50f5a95fdf7796c16f6ff Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 00:40:10 -0400 Subject: [PATCH 14/35] fix tests for data driven reduction and include tests for data informed reduction --- test/data_driven_approximation.jl | 4 +--- test/data_informed_approximation.jl | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/test/data_driven_approximation.jl b/test/data_driven_approximation.jl index 6412291..070004c 100644 --- a/test/data_driven_approximation.jl +++ b/test/data_driven_approximation.jl @@ -29,6 +29,4 @@ discrete_sol = LowRankIntegrators.solve(discrete_prob, solver) @test all(Matrix(sol.Y[end]) .≈ Matrix(discrete_sol.Y[end])) end -end - -#sol = LowRankIntegrators.solve(prob, solver, 1e-2) \ No newline at end of file +end \ No newline at end of file diff --git a/test/data_informed_approximation.jl b/test/data_informed_approximation.jl index 671bde2..52fbb5a 100644 --- a/test/data_informed_approximation.jl +++ b/test/data_informed_approximation.jl @@ -55,7 +55,7 @@ true_sols = [zeros(n,m^2) for i in 1:length(t_range)] ode_prob = ODEProblem(burgers!, ones(n), (0.0, 1.0), (∇,Δ,C,D)) for (i,ξ) in enumerate(ξs) - println(i) + println("Uncertainty realization $i") _prob = remake(ode_prob, u0 = ub(x_range) + uprime(x_range, ξ, σ)) _sol = Array(DifferentialEquations.solve(_prob, saveat=t_range)) for k in 1:length(t_range) From d74b9a41244749229efba94945df71d88c9810cb Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 00:49:53 -0400 Subject: [PATCH 15/35] add sparse interpolation tests and clean up --- src/sparse_interpolation.jl | 2 +- test/runtests.jl | 3 ++- test/sparse_interpolation.jl | 49 ++++++++++++++++++++++++++++++++++++ 3 files changed, 52 insertions(+), 2 deletions(-) diff --git a/src/sparse_interpolation.jl b/src/sparse_interpolation.jl index 417e587..88ab517 100644 --- a/src/sparse_interpolation.jl +++ b/src/sparse_interpolation.jl @@ -1,4 +1,4 @@ -import LinearAlgebra.adjoint +import LinearAlgebra: adjoint, rank import Base.size """ diff --git a/test/runtests.jl b/test/runtests.jl index eeb9b1c..06ecca2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,4 +3,5 @@ using Test include("data_agnostic_approximation.jl") include("data_driven_approximation.jl") -include("data_informed_approximation.jl") \ No newline at end of file +include("data_informed_approximation.jl") +include("sparse_interpolation.jl") \ No newline at end of file diff --git a/test/sparse_interpolation.jl b/test/sparse_interpolation.jl index e69de29..da0607f 100644 --- a/test/sparse_interpolation.jl +++ b/test/sparse_interpolation.jl @@ -0,0 +1,49 @@ +# a few lil tests +@testset "Sparse Interpolation Tests" begin + function rows!(dx,x,p,t,idcs) + @views dx .= x[idcs, :] + end + function cols!(dx,x,p,t,idcs) + @views dx .= x[:, idcs] + end + function elements!(dx,x,p,t,idcsA, idcsB) + @views dx .= x[idcsA, idcsB] + end + + F = ComponentFunction(rows!, cols!, elements!, 5, 8) + A = Matrix(reshape(1.0:5*8, (5, 8))) + + row_idcs = [1, 5, 2] + range_weights = zeros(eltype(A), 5, 3) + range_weights[row_idcs,1:3] .= Matrix{eltype(A)}(I,3,3) + + col_idcs = [1, 8, 3, 6] + corange_weights = zeros(eltype(A),8, 4) + corange_weights[col_idcs,1:4] .= Matrix{eltype(A)}(I,4,4) + + Π_range = DEIMInterpolator(row_idcs,range_weights) + Π_corange = DEIMInterpolator(col_idcs,corange_weights) + Π_mat = SparseMatrixInterpolator(Π_range, Π_corange) + + @test all(Π_range(A)[row_idcs,:] .== A[row_idcs,:]) + @test all(Π_corange'(A)[:,col_idcs] .== A[:,col_idcs]) + + F_range_int = SparseFunctionInterpolator(F, Π_range, output = eltype(A)) + F_corange_int = SparseFunctionInterpolator(F, Π_corange', output = eltype(A)) + F_int = SparseFunctionInterpolator(F, Π_mat, output = eltype(A)) + + eval!(F_range_int, A, (), 0.0) + @test all(F_range_int.cache.rows .== A[row_idcs,:]) + eval!(F_corange_int, A, (), 0.0) + @test all(F_corange_int.cache.cols .== A[:,col_idcs]) + eval!(F_int, A, (), 0.0) + @test all(F_int.cache.elements .== A[row_idcs,col_idcs]) + + dA = similar(A) + F_range_int(dA, A, (), 0.0) + @test all(dA[row_idcs, :] .== A[row_idcs, :]) + F_corange_int(dA, A, (), 0.0) + @test all(dA[:, col_idcs] .== A[:, col_idcs]) + F_int(dA, A, (), 0.0) + @test all(dA[row_idcs, col_idcs] .== A[row_idcs, col_idcs]) +end From 8597b2fcd37889b7cd3d3dd42bb4cce2eec9fb80 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 01:13:21 -0400 Subject: [PATCH 16/35] rm reexport of LinearAlgebra => extension of rank works as intended now --- src/LowRankIntegrators.jl | 4 ++-- src/integrators/greedy_integrator.jl | 2 +- src/integrators/rank_adaptive_unconventional.jl | 2 +- test/runtests.jl | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/LowRankIntegrators.jl b/src/LowRankIntegrators.jl index 40bb65d..04558a1 100644 --- a/src/LowRankIntegrators.jl +++ b/src/LowRankIntegrators.jl @@ -1,8 +1,8 @@ module LowRankIntegrators -using Reexport, UnPack, DocStringExtensions +using Reexport, UnPack, DocStringExtensions, LinearAlgebra import DifferentialEquations: step!, set_u!, init -@reexport using LinearAlgebra, DifferentialEquations, LowRankArithmetic +@reexport using DifferentialEquations, LowRankArithmetic include("utils.jl") export orthonormalize!, GradientDescent, QR, SVD, SecondMomentMatching, normal_component diff --git a/src/integrators/greedy_integrator.jl b/src/integrators/greedy_integrator.jl index 8ff2a5c..ef2e79f 100644 --- a/src/integrators/greedy_integrator.jl +++ b/src/integrators/greedy_integrator.jl @@ -31,7 +31,7 @@ end function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::SVDLikeRepresentation, dt; t0 = prob.tspan[1]) X = zeros(size(u)) - r = LowRankArithmetic.rank(u) + r = rank(u) n, m = size(X) XV = zeros(n,r) XU = zeros(m,r) diff --git a/src/integrators/rank_adaptive_unconventional.jl b/src/integrators/rank_adaptive_unconventional.jl index ecaf6d9..19a5eef 100644 --- a/src/integrators/rank_adaptive_unconventional.jl +++ b/src/integrators/rank_adaptive_unconventional.jl @@ -47,7 +47,7 @@ end function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) # allocate memory for frequently accessed arrays tspan = (t0, t0+dt) - r = LowRankArithmetic.rank(u) + r = rank(u) n, m = size(u) US = u.U*u.S VS = u.V*u.S' diff --git a/test/runtests.jl b/test/runtests.jl index 06ecca2..81434bf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -using LowRankIntegrators +using LinearAlgebra, LowRankIntegrators using Test include("data_agnostic_approximation.jl") From d72312571e8ea618fe429f654ff2e7e6d109f9aa Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 01:21:24 -0400 Subject: [PATCH 17/35] clean up --- src/integrators/do_integrators.jl | 250 ------------------ src/integrators/greedy_integrator.jl | 12 - src/integrators/projector_splitting.jl | 10 - .../rank_adaptive_unconventional.jl | 12 - src/integrators/unconventional.jl | 12 - src/primitives.jl | 7 - src/sparse_interpolation.jl | 8 - 7 files changed, 311 deletions(-) delete mode 100644 src/integrators/do_integrators.jl diff --git a/src/integrators/do_integrators.jl b/src/integrators/do_integrators.jl deleted file mode 100644 index 3505e6a..0000000 --- a/src/integrators/do_integrators.jl +++ /dev/null @@ -1,250 +0,0 @@ -# implements basic capabilities for DO integrators -abstract type AbstractRetraction end -struct DirectTimeMarching <: AbstractRetraction end -struct TruncatedSVD <: AbstractRetraction end -struct AlternatingProjection{oType} <: AbstractRetraction - order::oType -end - -struct Primal end -struct Dual end -struct PrimalDual end - -mutable struct DOAlgorithm{ρType} <: AbstractDLRAlgorithm - ρ::ρType # retraction - alg_params # parameters -end - -#= -Direct time marching - U_{n+1} = U_{n} + dt * L_Z*(L_U'*U_{n}) - Z_{n+1} = Z_{n} + dt * ((I-U_{n}*U_{n}')*L_U)(L_Z'*Z*(Z'*Z)^-1) - orthogonalize(U_{n+1}, Z_{n+1}) - repeat -=# - -struct DirectTimeMarching_Params - F - orth_alg -end - -struct DirectTimeMarching_Cache - F - Y - LU - LZ - C -end - -function DirectTimeMarching(F; orthonormalization = QRFact()) - return DOAlgorithm(DirectTimeMarching(), DirectTimeMarching_Params(F, orthonormalization)) -end - -function init(prob::MatrixDEProblem, alg::DOAlgorithm{<:AbstractRetraction}, dt) - t0, tf = prob.tspan - u = deepcopy(prob.u0) - @assert tf > t0 "Integration in reverse time direction is not supported" - # number of steps - n = floor(Int,(tf-t0)/dt) + 1 - # compute more sensible dt # rest will be handled via interpolation/save_at - dt = (tf-t0)/(n-1) - # initialize solution object - sol = DLRSolution(Vector{typeof(prob.u0)}(undef, n), collect(range(t0, tf, length=n))) - sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object - # initialize cache - cache = alg_cache(prob, alg, u, dt) - return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -end - -function alg_cache(prob::MatrixDEProblem, alg::DOAlgorithm{DirectTimeMarching}, u, dt) - @unpack alg_params = alg - F = (u,t,dt) -> alg_params.F(prob.f, u, t, dt) - L = F(u, prob.tspan[1] ,dt) - rL = rank(L) - r = rank(u) - LU = zeros(rL,r) - LZ = zeros(rL,r) - Y = deepcopy(u) - C = zeros(r, r) - return DirectTimeMarching_Cache(F, Y, LU, LZ, C) -end - -function step!(integrator::DLRIntegrator, alg::DOAlgorithm{DirectTimeMarching}, dt) - @unpack u, t, iter, cache = integrator - direct_time_marching_do_step!(u, cache, t, dt) - orthonormalize!(u, alg.alg_params.orth_alg) - cache.Y.Z .= u.Z - cache.Y.U .= u.U - integrator.t += dt - integrator.iter += 1 -end - -function direct_time_marching_do_step!(u, cache, t, dt) - @unpack Y, LU, LZ, C, F = cache - L = F(Y,t,dt) - - # Z update - mul!(LU, L.U', Y.U) - mul!(u.Z, L.Z, LU, dt, 1) - - # U update - #((I-U*U')*L_U)(L_Z'*Z*(Z'*Z)^-1) - mul!(L.U, Y.U, LU', -1, 1) #(I-U*U')*L_U - mul!(LZ, L.Z', Y.Z) # L_Z'*Z - mul!(C, Y.Z', Y.Z) # covariance computation - LU .= LZ / Symmetric(C) # revisit - mul!(u.U, L.U, LU, dt, 1) -end - -#= -Alternating projection - Primal step: - LRA_{n+1} = LRA_{n} + dt*F(LRA_{n},t) - solve min \|LRA{n+1} - U_{n} Z_{n+1}'\| => Z_{n+1}' = U_{n}'*LRA_{n+1} => Z_{n+1} = LRA_{n+1}'*U{n} - solve min \|LRA{n+1} - U_{n+1} Z_{n+1}'\| - => procrustes problem: Q, _, P = svd(LRA_{n+1}*Z_{n+1}), U_{n+1}= Q*P' - Dual step: - LRA_{n+1} = LRA_{n} + dt*F(LRA_{n},t) - solve min \|LRA{n+1} - U_{n+1} Z_{n}'\| - => procrustes problem: Q, _, P = svd(LRA_{n+1}*Z_{n}), U_{n+1}= Q*P' - solve min \|LRA{n+1} - U_{n+1} Z_{n+1}'\| => Z_{n+1} = U_{n+1}'*LRA_{n+1} -=# - -struct AlternatingProjection_Params - F - orth_alg -end - -struct AlternatingProjection_Cache - F - UtU - UQ - C - orth_alg -end - -function AlternatingProjection(order, F; orthonormalization = QRFact()) - return DOAlgorithm(AlternatingProjection(order), AlternatingProjection_Params(F, orthonormalization)) -end - -function alg_cache(prob::MatrixDEProblem, alg::DOAlgorithm{AlternatingProjection{T}}, u, dt) where T - @unpack alg_params = alg - F = (u,t,dt) -> alg_params.F(prob.f, u, t, dt) - Y = u + F(u, prob.tspan[1] ,dt) - rY = rank(Y) - r = rank(u) - UtU = zeros(rY, r) - C = zeros(rY, r) - UQ = zeros(size(Y,1), r) - return AlternatingProjection_Cache(F, UtU, UQ, C, alg_params.orth_alg) -end - -function step!(integrator::DLRIntegrator, alg::DOAlgorithm{AlternatingProjection{Primal}}, dt) - @unpack u, t, iter, cache = integrator - primal_alternating_projection_do_step!(u, cache, t, dt) - integrator.t += dt - integrator.iter += 1 -end - -function step!(integrator::DLRIntegrator, alg::DOAlgorithm{AlternatingProjection{Dual}}, dt) - @unpack u, t, iter, cache = integrator - dual_alternating_projection_do_step!(u, cache, t, dt) - integrator.t += dt - integrator.iter += 1 -end - -function step!(integrator::DLRIntegrator, alg::DOAlgorithm{AlternatingProjection{PrimalDual}}, dt) - @unpack u, t, iter, cache = integrator - primal_alternating_projection_do_step!(u, cache, t, dt/2) - dual_alternating_projection_do_step!(u, cache, t + dt/2, dt/2) - integrator.t += dt - integrator.iter += 1 -end - -function primal_alternating_projection_do_step!(u, cache, t, dt) - @unpack F, UtU, UQ, C, orth_alg = cache - Y = u + dt*F(u,t,dt) # investigate using the integrator interface of DifferentialEquations.jl - orthonormalize!(Y, orth_alg) - - # coefficient update - mul!(UtU, Y.U', u.U) - mul!(u.Z, Y.Z, UtU) - - # basis update (procrustes problem) - mul!(C, Y.Z', u.Z) - Q, _, P = svd!(C) - mul!(UQ, Y.U, Q) - mul!(u.U, UQ, P') -end - -function dual_alternating_projection_do_step!(u, cache, t, dt) - @unpack F, UtU, UQ, C, orth_alg = cache - Y = u + dt*F(u,t,dt) # investigate using the integrator interface of DifferentialEquations.jl - orthonormalize!(Y, orth_alg) - - # basis update (procrustes problem) - mul!(C, Y.Z', u.Z) - Q, _, P = svd!(C) - mul!(UQ, Y.U, Q) - mul!(u.U, UQ, P') - - # coefficient update - mul!(UtU, Y.U', u.U) - mul!(u.Z, Y.Z, UtU) -end - -#= -SVD => (svd, tsvd, gradient flow) - LRA_{n+1} = LRA_{n} + dt*F(LRA_{n},t) - LRA_{n+1} = truncated_svd(LRA_{n+1}, rmax = rank(LRA_{n+1})) - or - LRA_{n+1} = truncated_svd(LRA_{n+1}, tol = #tolerance, rmax = #max rank) -=# - -struct TruncatedSVD_Params - F - SVD_alg - orth_alg - tol - rmin - rmax -end - -struct TruncatedSVD_Cache - F -end - -function TruncatedSVD(F; SVD_alg = SVDFact(), orthonormalization = QRFact(), - rmin::Int=1, rmax::Int=Int(1e8), tol = sqrt(eps(Float64))) - return DOAlgorithm(TruncatedSVD(), TruncatedSVD_Params(F, SVD_alg, orthonormalization, tol, rmin, rmax)) -end - -function alg_cache(prob::MatrixDEProblem, alg::DOAlgorithm{TruncatedSVD}, u, dt) - @unpack alg_params = alg - F = (u,t,dt) -> alg_params.F(prob.f, u, t, dt) - return TruncatedSVD_Cache(F) -end - -function step!(integrator::DLRIntegrator, alg::DOAlgorithm{TruncatedSVD}, dt) - @unpack u, t, iter, cache = integrator - truncated_svd_do_step!(u, cache, t, dt, alg.alg_params) - integrator.t += dt - integrator.iter += 1 -end - -function truncated_svd_do_step!(u, cache, t, dt, params) - @unpack F = cache - @unpack tol, rmin, rmax, SVD_alg, orth_alg = params - Y = u + dt*F(u, t, dt) - retract!(Y, u, SVD_alg, tol, rmin, rmax, orth_alg) -end - -function retract!(Y, u, alg::T, tol, rmin, rmax, orth_alg) where T <: Union{SVDFact, TSVD} - Y = round(Y, alg, tol=tol, rmin = rmin, rmax = rmax, alg_orthonormalize = orth_alg) - u.Z = Y.Z - u.U = Y.U -end - -function retract!(Y, u, alg::GradientDescent, tol, rmin, rmax, orth_alg) - gd_truncated_svd!(u, Y, alg, orth_alg) -end \ No newline at end of file diff --git a/src/integrators/greedy_integrator.jl b/src/integrators/greedy_integrator.jl index ef2e79f..c78b3ee 100644 --- a/src/integrators/greedy_integrator.jl +++ b/src/integrators/greedy_integrator.jl @@ -46,18 +46,6 @@ function alg_cache(prob::MatrixHybridProblem, alg::GreedyIntegrator, u, dt; t0 = return GreedyIntegrator_Cache(prob.y, X, XZ, nothing, nothing, ZIntegrator) end -# function init(prob::MatrixDataProblem, alg::GreedyIntegrator, dt) -# t0, tf = prob.tspan -# @assert tf > t0 "Integration in reverse time direction is not supported" -# u = deepcopy(prob.u0) -# # initialize solution -# sol = init_sol(dt, t0, tf, prob.u0) -# # initialize cache -# cache = alg_cache(prob, alg, u, dt) -# sol.Y[1] = deepcopy(prob.u0) -# return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -# end - function init(prob::MatrixHybridProblem, alg::GreedyIntegrator, dt) t0, tf = prob.tspan @assert tf > t0 "Integration in reverse time direction is not supported" diff --git a/src/integrators/projector_splitting.jl b/src/integrators/projector_splitting.jl index adc8cba..0d8b483 100644 --- a/src/integrators/projector_splitting.jl +++ b/src/integrators/projector_splitting.jl @@ -104,16 +104,6 @@ function alg_cache(prob::MatrixDataProblem, alg::ProjectorSplitting, u, dt; t0 = prob.y, ycurr, yprev, Δy) end -# function init(prob::AbstractDLRProblem, alg::ProjectorSplitting, dt) -# t0, tf = prob.tspan -# @assert tf > t0 "Integration in reverse time direction is not supported" -# u = deepcopy(prob.u0) -# sol = init_sol(dt, t0, tf, prob.u0) -# cache = alg_cache(prob, alg, u, dt, t0 = t0) -# sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object -# return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -# end - function primal_LT_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem}) @unpack y, ycurr, yprev, Δy = cache update_data!(ycurr, y, t, dt) diff --git a/src/integrators/rank_adaptive_unconventional.jl b/src/integrators/rank_adaptive_unconventional.jl index 19a5eef..905e884 100644 --- a/src/integrators/rank_adaptive_unconventional.jl +++ b/src/integrators/rank_adaptive_unconventional.jl @@ -91,18 +91,6 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit nothing, nothing, nothing, nothing) end -# function init(prob::AbstractDLRProblem, alg::RankAdaptiveUnconventionalAlgorithm, dt) -# t0, tf = prob.tspan -# @assert tf > t0 "Integration in reverse time direction is not supported" -# u = deepcopy(prob.u0) -# # initialize solution -# sol = init_sol(dt, t0, tf, prob.u0) -# # initialize cache -# cache = alg_cache(prob, alg, u, dt; t0 = t0) -# sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object -# return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -# end - function alg_cache(prob::MatrixDataProblem, alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) # creates caches for frequently used arrays by performing the first time step @unpack y = prob diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index 44cfc1d..333b5a5 100644 --- a/src/integrators/unconventional.jl +++ b/src/integrators/unconventional.jl @@ -219,18 +219,6 @@ function alg_cache(prob::MatrixDataProblem, alg::UnconventionalAlgorithm, u, dt; prob.y, ycurr, yprev, Δy, nothing) end -# function init(prob::AbstractDLRProblem, alg::UnconventionalAlgorithm, dt) -# t0, tf = prob.tspan -# @assert tf > t0 "Integration in reverse time direction is not supported" -# u = deepcopy(prob.u0) -# # initialize solution -# sol = init_sol(dt, t0, tf, prob.u0) -# # initialize cache -# cache = alg_cache(prob, alg, u, dt, t0 = t0) -# sol.Y[1] = deepcopy(prob.u0) -# return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0) -# end - function unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem}) @unpack y, ycurr, yprev, Δy = cache update_data!(ycurr, y, t, dt) diff --git a/src/primitives.jl b/src/primitives.jl index 33fa8e3..806135d 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -147,11 +147,4 @@ function update_sol!(integrator::AbstractDLRIntegrator, idx) integrator.sol.r[idx] = LowRankArithmetic.rank(integrator.u) integrator.sol.r_DEIM[idx] = rank_DEIM(integrator.cache) integrator.sol.idcs_DEIM[idx] = interpolation_indices(integrator.cache) - #else - # push!(integrator.sol.Y, deepcopy(integrator.u)) - # push!(integrator.sol.t, integrator.t) - # push!(integrator.sol.r[integrator.iter + 1], LowRankArithmetic.rank(integrator.u)) - # push!(integrator.sol.r_DEIM[integrator.iter + 1], rank_DEIM(integrator.cache)) - # push!(integrator.sol.r_DEIM[integrator.iter + 1], indices(integrator.cache)) - # end end \ No newline at end of file diff --git a/src/sparse_interpolation.jl b/src/sparse_interpolation.jl index 88ab517..7d6448d 100644 --- a/src/sparse_interpolation.jl +++ b/src/sparse_interpolation.jl @@ -433,18 +433,10 @@ function index_selection(U::Matrix, S::Vector, alg::LDEIM) r = min(alg.rmax, length(S)+1) return index_selection(U, r, alg) else - #cutoff = findfirst(x -> x < alg.tol*alg.elasticity, S) cutoff = findlast(x -> x > alg.tol*alg.elasticity, S) if isnothing(cutoff) cutoff = 1 end - #cutoff = length(S) - # for σ in reverse(S) - # if σ > alg.tol*alg.elasticity - # break - # end - # cutoff -= 1 - # end return index_selection(U, cutoff, alg) end end From a30bc7fe739df14ef58608e944df8a82a894cc71 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 14:10:11 -0400 Subject: [PATCH 18/35] update to LowRankArithmetic 0.1.3 and Julia 1.8 compats --- Project.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 3879061..d09dd9f 100644 --- a/Project.toml +++ b/Project.toml @@ -13,10 +13,9 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] DifferentialEquations = "7" -LowRankArithmetic = "0.1.2" Reexport = "1" UnPack = "1" -julia = "1.6" +julia = "1.8" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 75e56e7e3779f316a2267eaf04c6d2fc921702fe Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 15:31:13 -0400 Subject: [PATCH 19/35] add test for sparse interpolation integrator --- test/data_agnostic_DEIM_approximation.jl | 236 +++++++++++++++++++++++ test/runtests.jl | 3 +- 2 files changed, 238 insertions(+), 1 deletion(-) create mode 100644 test/data_agnostic_DEIM_approximation.jl diff --git a/test/data_agnostic_DEIM_approximation.jl b/test/data_agnostic_DEIM_approximation.jl new file mode 100644 index 0000000..20f2bce --- /dev/null +++ b/test/data_agnostic_DEIM_approximation.jl @@ -0,0 +1,236 @@ +using LowRankIntegrators, LinearAlgebra + +@testset "DLRA + on-the-fly DEIM for Burgers' equation" begin + # Model parameters + N = 1024 # spatial discretization + S = 256 # number of samples + σ_t = 0.01 # 0.01 # noise level in forcing + σ_x = 0.005 # noise level in initial condition + ν = 0.0025 #0.0025 # viscosity + d = 4 # dimension of uncertainty + + # the rest is computed + x_range = range(0+1/(N+1), 1-1/(N+1), length=N) # spatial grid + Δx = step(x_range) # step size + ξ = [randn(d) for i in 1:S-1] # random samples + pushfirst!(ξ, zeros(d)) # add mean + ξ_mat = reduce(hcat, ξ) + + # for initial condition uncertainty: + gauss_kernel(x,z,σ) = exp(-(x-z)^2/(2*σ^2)) + M = [gauss_kernel(x,z,σ_x) for x in x_range, z in x_range] + Λ, Ψ = eigen(M, sortby=-) + x0_mean = @. 0.5*sin(2π*x_range) * (exp(cos(2π*x_range)) - 1.5) + X0_noise = SVDLikeRepresentation(Ψ[:,1:d], diagm(sqrt.(Λ[1:d])), ξ_mat') + X0 = TwoFactorRepresentation(x0_mean, ones(S)) + X0_noise + + @inline function D²(xleft, xmid, xright, Δx) + return (xleft + xright - 2*xmid)/Δx^2 + end + + @inline function D(xleft, xmid, Δx) + return (xmid - xleft)/Δx + end + + function burgers_col!(dx, x, p, t, col) + Δx, ν, ξ, σ_t, d, n = p + n = size(x,1) + x_bc = x_left(t, ξ[col], σ_t, d) + dx[1] = ν*D²(x_bc, x[1], x[2], Δx) - + x[1]*D(x_bc, x[1], Δx) + for i in 2:n-1 + dx[i] = ν*D²(x[i-1], x[i], x[i+1], Δx) - + x[i]*D(x[i-1], x[i], Δx) + end + dx[n] = ν*D²(x[n-1], x[n], 0, Δx) - + x[n]*D(x[n-1], x[n], Δx) + end + + function burgers_cols!(dx, x::TwoFactorRepresentation, p, t, cols) + core_cache, col_cache, params = p + # cannot use threads ... or otherwise too much allocation + @views for (i,col) in enumerate(cols) + mul!(col_cache, x.U, x.Z[col,:]) + burgers_col!(dx[:,i], col_cache, params, t, col) #col_cache, params, t, col) + end + end + + function burgers_cols!(dx, x::SVDLikeRepresentation, p, t, cols) + core_cache, col_cache, params = p + mul!(core_cache, x.U, x.S) + @views for (i,col) in enumerate(cols) + mul!(col_cache, core_cache, x.V[col,:]) + @views burgers_col!(dx[:,i], col_cache, params, t, col) + end + end + + function burgers_row!(dx, x, p, t, row) + Δx, ν, ξ, σ_t, d, n = p + if row == 1 + for j in eachindex(dx) + x_bc = x_left(t, ξ[j], σ_t, d) + dx[j] = ν*D²(x_bc, x[1,j], x[2,j], Δx) - x[1,j]*D(x_bc, x[1,j], Δx) + end + elseif row == n + for j in eachindex(dx) + dx[j] = ν*D²(x[1,j], x[2,j], 0, Δx) - x[2,j]*D(x[1,j], x[2,j], Δx) + end + else + for j in eachindex(dx) + dx[j] = ν*D²(x[1,j], x[2,j], x[3,j], Δx) - x[2,j]*D(x[1,j], x[2,j], Δx) + end + end + end + + function burgers_rows!(dx, x::TwoFactorRepresentation, p, t, rows) + row_cache, col_cache, params = p + n = params[end] + @views for (i,row) in enumerate(rows) + if row == 1 + mul!(row_cache[1:2,:], x.U[1:2,:], x.Z') + elseif row == n + mul!(row_cache[1:2,:], x.U[n-1:n,:], x.Z') + else + mul!(row_cache, x.U[row-1:row+1,:], x.Z') + end + burgers_row!(dx[i,:], row_cache, params, t, row) + end + end + + function burgers_rows!(dx, x::SVDLikeRepresentation, p, t, rows) + row_cache, col_cache, params = p + n = params[end] + mul!(col_cache, x.S, x.V') + @views for (i,row) in enumerate(rows) + if row == 1 + mul!(row_cache[1:2,:], x.U[1:2,:], col_cache) + elseif row == n + mul!(row_cache[1:2,:], x.U[n-1:n,:], col_cache) + else + mul!(row_cache, x.U[row-1:row+1,:], col_cache) + end + burgers_row!(dx[i,:], row_cache, params, t, row) + end + end + + function burgers_elements!(dx, x::SVDLikeRepresentation, p, t, rows, cols) + col_cache, row_cache, params = p + n = params[end] + @views for (j, col) in enumerate(cols) + mul!(col_cache, x.S, x.V[col,:]) + for (i, row) in enumerate(rows) + if row == 1 + mul!(row_cache[1:2,:], x.U[1:2,:], col_cache) + elseif row == n + mul!(row_cache[1:2,:], x.U[n-1:n,:], col_cache) + else + mul!(row_cache, x.U[row-1:row+1,:], col_cache) + end + dx[i,j] = burgers_rhs(row_cache, params, t, row, col) + end + end + end + + function burgers_elements!(dx, x::TwoFactorRepresentation, p, t, rows, cols) + col_cache, row_cache, neighbor_cache, params = p + @views for (i, row) in enumerate(rows), (j, col) in enumerate(cols) + if row == 1 + mul!(row_cache[1:2,:], x.U[1:2,:], x.Z[col,:]) + elseif row == n + mul!(row_cache[1:2,:], x.U[n-1:n,:], x.Z[col,:]) + else + mul!(row_cache, x.U[row-1:row+1,:], x.Z[col,:]) + end + dx[i,j] = burgers_rhs(row_cache, params, t, row, col) + end + end + + function burgers_rhs(x, p, t, row, col) + Δx, ν, ξ, σ_t, d, n = p + if row == 1 + x_bc = x_left(t, ξ[col], σ_t, d) + return ν*D²(x_bc, x[1], x[2], Δx) - x[1]*D(x_bc, x[1], Δx) + elseif row == n + return ν*D²(x[1], x[2], 0, Δx) - x[2]*D(x[1], x[2], Δx) + else + return ν*D²(x[1], x[2], x[3], Δx) - x[2]*D(x[1], x[2], Δx) + end + end + + x_left(t, ξ, σ_t, d) = -0.4*sin(2π*t) + σ_t * sum(1/i^2*ξ[i]*sin(i*π*t) for i in 1:d) + + # solve parameters + r = 5 + r_deim = 15 + params = (Δx, ν, ξ, σ_t, d, N) + rows_cache = (zeros(3,S),zeros(r,S), params) + cols_cache = (zeros(N,r), zeros(N), params) + elements_cache = (zeros(r), zeros(3), params) + + rows! = (dx, x, p, t, rows) -> burgers_rows!(dx, x, rows_cache, t, rows) + cols! = (dx, x, p, t, cols) -> burgers_cols!(dx, x, cols_cache, t, cols) + elements! = (dx, x, p, t, rows, cols) -> burgers_elements!(dx, x, elements_cache, t, rows, cols) + burgers_components = ComponentFunction(rows!,cols!,elements!, N, S) + + X0_lr = truncated_svd(Matrix(X0),r) + dX0 = zeros(N, S) + @views for k in axes(dX0, 2) + burgers_components.cols!(dX0[:,k], X0_lr, (), 0.0, [k]) + end + + dX0_lr = truncated_svd(dX0, r_deim) + # first test projections + idx_selections = [DEIM(), QDEIM(), LDEIM()] + init_range, init_corange = dX0_lr.U, dX0_lr.V + U, V = X0_lr.U, X0_lr.V + dX_test = similar(dX0) + for alg in idx_selections + row_idcs = index_selection(dX0_lr.U, alg) + col_idcs = index_selection(dX0_lr.V, alg) + + # interpolators + Π_corange = DEIMInterpolator(col_idcs, init_corange/init_corange[col_idcs,:]) + Π_range = DEIMInterpolator(row_idcs, init_range/init_range[row_idcs,:]) + Π = SparseFunctionInterpolator(burgers_components, SparseMatrixInterpolator(Π_range, Π_corange)) + + # core approximation + dX_test = zeros(N,S) + Π(dX_test, X0_lr, (), 0.0) + @test norm(dX_test - dX0_lr)/norm(Matrix(dX0_lr)) < 1e-8 + + # low rank integration subproblem rhs + Π_K_corange = DEIMInterpolator(col_idcs,V'*Π.interpolator.corange.weights)' + Π_K = SparseFunctionInterpolator(burgers_components, Π_K_corange) + + # K step + dK_test = zeros(N,r) + K0 = X0_lr.U*X0_lr.S + Π_K(dK_test, TwoFactorRepresentation(K0, X0_lr.V), (), 0.0) + @test norm(dK_test - dX0_lr*V)/norm(Matrix(dX0_lr*V)) < 1e-8 + + # L step + Π_L_range = DEIMInterpolator(row_idcs, U'*Π.interpolator.range.weights) + Π_L = SparseFunctionInterpolator(burgers_components, Π_L_range) + dL_test = zeros(r,S) + L0 = X0_lr.V*X0_lr.S' + Π_L(dL_test, TwoFactorRepresentation(X0_lr.U,L0), (), 0.0) + @test norm(dL_test - U'*dX0_lr)/norm(Matrix(U'*dX0_lr)) < 1e-8 + + # S step + Π_S_mat = SparseMatrixInterpolator(row_idcs, col_idcs, + U'*Π.interpolator.range.weights, + V'*Π.interpolator.corange.weights) + Π_S = SparseFunctionInterpolator(burgers_components, Π_S_mat) + dS_test = zeros(r,r) + S0 = X0_lr.S + Π_S(dS_test, X0_lr, (), 0.0) + @test norm(dS_test - U'*dX0_lr*V)/norm(Matrix(U'*X0_lr*V)) < 1e-8 + end + + r_deim = 5 + interpolation = SparseInterpolation(DEIM(rmax=r_deim, rmin=r_deim, tol = 1.0, elasticity=1e-1), dX0_lr.U[:,1:r_deim], dX0_lr.V[:,1:r_deim]) + alg_deim = UnconventionalAlgorithm(interpolation) + + deim_prob = MatrixDEProblem(burgers_components, X0_lr, (0.0,1.0)) + deim_sol = LowRankIntegrators.solve(deim_prob, alg_deim, 1e-3, save_increment=10) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 81434bf..8fda4f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,4 +4,5 @@ using Test include("data_agnostic_approximation.jl") include("data_driven_approximation.jl") include("data_informed_approximation.jl") -include("sparse_interpolation.jl") \ No newline at end of file +include("data_agnostic_DEIM_approximation.jl") +include("sparse_interpolation.jl") From 11a29ac7bcd9c5d97c72de4888720dd62630d7a6 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 15:48:40 -0400 Subject: [PATCH 20/35] change CI to Julia 1.8 --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f86d22a..e421b31 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -18,7 +18,7 @@ jobs: fail-fast: false matrix: version: - - '1.6' + - '1.8' #- 'nightly' avoid for now os: - ubuntu-latest From f4debd76831ac99a1692978843acd5e6b7bdd568 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 16:29:44 -0400 Subject: [PATCH 21/35] fix sparse interpolation lr integration test --- test/data_agnostic_DEIM_approximation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/data_agnostic_DEIM_approximation.jl b/test/data_agnostic_DEIM_approximation.jl index 20f2bce..68bba09 100644 --- a/test/data_agnostic_DEIM_approximation.jl +++ b/test/data_agnostic_DEIM_approximation.jl @@ -185,8 +185,8 @@ using LowRankIntegrators, LinearAlgebra U, V = X0_lr.U, X0_lr.V dX_test = similar(dX0) for alg in idx_selections - row_idcs = index_selection(dX0_lr.U, alg) - col_idcs = index_selection(dX0_lr.V, alg) + row_idcs = index_selection(init_range, alg) + col_idcs = index_selection(init_corange, alg) # interpolators Π_corange = DEIMInterpolator(col_idcs, init_corange/init_corange[col_idcs,:]) From 0a0cb3a27903317ab3eaa6b04cbdcea37fb54360 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 17:04:59 -0400 Subject: [PATCH 22/35] debugging index selection errors --- test/runtests.jl | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 8fda4f6..908eea4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,8 +1,15 @@ using LinearAlgebra, LowRankIntegrators using Test -include("data_agnostic_approximation.jl") -include("data_driven_approximation.jl") -include("data_informed_approximation.jl") -include("data_agnostic_DEIM_approximation.jl") -include("sparse_interpolation.jl") +X = randn(10,10) +test = svd(X) +algs = [QDEIM(), DEIM(), LDEIM()] +for alg in algs + idcs = index_selection(test.U[:,1:5], alg) + @test length(idcs) == 5 +end +#include("data_agnostic_approximation.jl") +#include("data_driven_approximation.jl") +#include("data_informed_approximation.jl") +#include("data_agnostic_DEIM_approximation.jl") +#include("sparse_interpolation.jl") From 3026ba08f9cfec7cc1cf9058604a581973ee9a55 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 18:14:15 -0400 Subject: [PATCH 23/35] revert to normal tests --- test/runtests.jl | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 908eea4..8fda4f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,15 +1,8 @@ using LinearAlgebra, LowRankIntegrators using Test -X = randn(10,10) -test = svd(X) -algs = [QDEIM(), DEIM(), LDEIM()] -for alg in algs - idcs = index_selection(test.U[:,1:5], alg) - @test length(idcs) == 5 -end -#include("data_agnostic_approximation.jl") -#include("data_driven_approximation.jl") -#include("data_informed_approximation.jl") -#include("data_agnostic_DEIM_approximation.jl") -#include("sparse_interpolation.jl") +include("data_agnostic_approximation.jl") +include("data_driven_approximation.jl") +include("data_informed_approximation.jl") +include("data_agnostic_DEIM_approximation.jl") +include("sparse_interpolation.jl") From b03d5260fb278bd4dc1c89189ba2f9db76df69a0 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 18:41:42 -0400 Subject: [PATCH 24/35] fix progress printouts --- src/primitives.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/primitives.jl b/src/primitives.jl index 806135d..a198f5e 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -72,7 +72,7 @@ function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt; save_inc println("Initialize integrator ...") integrator, dt, t_int, save_idcs = init(prob, alg, dt, save_increment) println("... initialization complete. Start Integration ...") - disp_digits = abs(round(Int, log10(save_increment*dt))) + disp_digits = -floor(Int, log10(save_increment*dt)) k = 2 #while (prob.tspan[2]-integrator.t)/T > 1e-8 while integrator.iter < length(t_int) - 1 @@ -80,7 +80,7 @@ function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt; save_inc if integrator.iter + 1 == save_idcs[k] update_sol!(integrator, k) k += 1 - println("t = $(round(integrator.t,sigdigits=disp_digits))") + println("t = $(round(integrator.t,digits=disp_digits))") end end println("... integration complete.") From fd1ab53b9a9ffb409a2acc25c00792a036b0ee23 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 18:50:33 -0400 Subject: [PATCH 25/35] rm test env dependence for LRI --- test/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/test/Project.toml b/test/Project.toml index b36202c..f63ddf1 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,5 @@ [deps] DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -LowRankIntegrators = "a4d2327a-b40c-407e-b18c-dbb78b346b8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 01b7a9c4fceedc84c8ba17b23a39e3e0d4d05841 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 20:16:49 -0400 Subject: [PATCH 26/35] rename UnconventionalDEIM_Cache to OnTheFlyInterpolation_Cache --- src/integrators.jl | 9 +++++++++ src/integrators/unconventional.jl | 15 +++------------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/integrators.jl b/src/integrators.jl index 840cacd..50f2f01 100644 --- a/src/integrators.jl +++ b/src/integrators.jl @@ -2,6 +2,15 @@ rank_DEIM(::AbstractDLRAlgorithm_Cache) = -1 interpolation_indices(::AbstractDLRAlgorithm_Cache) = (Int[], Int[]) +struct OnTheFlyInterpolation_Cache + params::SparseInterpolation + u_prev::SVDLikeRepresentation + Π::SparseFunctionInterpolator + Π_K::SparseFunctionInterpolator + Π_L::SparseFunctionInterpolator + Π_S::SparseFunctionInterpolator +end + include("integrators/projector_splitting.jl") include("integrators/unconventional.jl") include("integrators/data_integrator.jl") diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index 333b5a5..47ed2dd 100644 --- a/src/integrators/unconventional.jl +++ b/src/integrators/unconventional.jl @@ -52,21 +52,12 @@ struct UnconventionalAlgorithm_Cache{uType, interpolation_cache::dType end -struct UnconventionalDEIM_Cache - params::SparseInterpolation - u_prev::SVDLikeRepresentation - Π::SparseFunctionInterpolator - Π_K::SparseFunctionInterpolator - Π_L::SparseFunctionInterpolator - Π_S::SparseFunctionInterpolator -end - rank_DEIM(cache::UnconventionalAlgorithm_Cache) = rank_DEIM(cache.interpolation_cache) rank_DEIM(::Nothing) = 0 -rank_DEIM(cache::UnconventionalDEIM_Cache) = rank(cache.Π)[1] +rank_DEIM(cache::OnTheFlyInterpolation_Cache) = rank(cache.Π)[1] interpolation_indices(cache::UnconventionalAlgorithm_Cache) = interpolation_indices(cache.interpolation_cache) interpolation_indices(::Nothing) = (Int[],Int[]) -interpolation_indices(cache::UnconventionalDEIM_Cache) = interpolation_indices(cache.Π) +interpolation_indices(cache::OnTheFlyInterpolation_Cache) = interpolation_indices(cache.Π) function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) # allocate memory for frequently accessed arrays @@ -188,7 +179,7 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, p_S) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) - interpolation_cache = UnconventionalDEIM_Cache(alg.alg_params.interpolation, deepcopy(u), + interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), Π, Π_K, Π_L, Π_S) return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, From c70a966efa73da0bed9e100da466800d76777bdd Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 21:26:43 -0400 Subject: [PATCH 27/35] restructure filesystem => on the fly interpolation into a separate file/chunk --- src/integrators.jl | 13 +- src/integrators/on_the_fly_interpolation.jl | 171 ++++++++++++++++++++ src/integrators/unconventional.jl | 161 ------------------ 3 files changed, 172 insertions(+), 173 deletions(-) create mode 100644 src/integrators/on_the_fly_interpolation.jl diff --git a/src/integrators.jl b/src/integrators.jl index 50f2f01..2db1c0d 100644 --- a/src/integrators.jl +++ b/src/integrators.jl @@ -1,16 +1,5 @@ # basic catch all cases: -rank_DEIM(::AbstractDLRAlgorithm_Cache) = -1 -interpolation_indices(::AbstractDLRAlgorithm_Cache) = (Int[], Int[]) - -struct OnTheFlyInterpolation_Cache - params::SparseInterpolation - u_prev::SVDLikeRepresentation - Π::SparseFunctionInterpolator - Π_K::SparseFunctionInterpolator - Π_L::SparseFunctionInterpolator - Π_S::SparseFunctionInterpolator -end - +include("integrators/on_the_fly_interpolation.jl") include("integrators/projector_splitting.jl") include("integrators/unconventional.jl") include("integrators/data_integrator.jl") diff --git a/src/integrators/on_the_fly_interpolation.jl b/src/integrators/on_the_fly_interpolation.jl new file mode 100644 index 0000000..df4cad2 --- /dev/null +++ b/src/integrators/on_the_fly_interpolation.jl @@ -0,0 +1,171 @@ +rank_DEIM(::AbstractDLRAlgorithm_Cache) = -1 +interpolation_indices(::AbstractDLRAlgorithm_Cache) = (Int[], Int[]) + +struct OnTheFlyInterpolation_Cache + params::SparseInterpolation + u_prev::SVDLikeRepresentation + Π::SparseFunctionInterpolator + Π_K::SparseFunctionInterpolator + Π_L::SparseFunctionInterpolator + Π_S::SparseFunctionInterpolator +end + +function select_idcs(F_rows::AbstractMatrix, F_cols::AbstractMatrix, selection_alg::IndexSelectionAlgorithm) + range = svd(F_cols) + corange = svd(F_rows) + S = max.(range.S, corange.S) + row_idcs = index_selection(range.U, S, selection_alg) + col_idcs = index_selection(Matrix(corange.V), S, selection_alg) + return range.U, corange.V, row_idcs, col_idcs +end + +function select_idcs(F::SVDLikeRepresentation, selection_alg::IndexSelectionAlgorithm) + F_svd = svd(F) + row_idcs = index_selection(F_svd.U, diag(F_svd.S), selection_alg) + col_idcs = index_selection(F_svd.V, diag(F_svd.S), selection_alg) + UF = F_svd.U[:, 1:length(row_idcs)] + VF = F_trunc.V[:, 1:length(col_idcs)] + return UF, VF, row_idcs, col_idcs +end + +function sparse_selection(interpolation_cache, u, t, dt) + @unpack params = interpolation_cache + if params.update_scheme == :avg_flow + @unpack u_prev = interpolation_cache + @unpack selection_alg = params + return select_idcs(u - u_prev, selection_alg) + elseif params.update_scheme == :last_iterate + @unpack Π_K, Π_L = interpolation_cache + @unpack selection_alg = params + eval!(Π_L, u, (), t+dt) + eval!(Π_K, u, (), t+dt) + UF, VF, row_idcs, col_idcs = select_idcs(Π_L.cache.rows, Π_K.cache.cols, selection_alg) + return UF, VF, row_idcs, col_idcs + else + error("Update scheme $(params.update_scheme) not supported") + end + + # eval index selection + +end + +function compute_weights(range, corange, row_idcs, col_idcs, u, t, dt, interpolation_cache) + @unpack params = interpolation_cache + if size(range, 2) != length(row_idcs) + @unpack Π = interpolation_cache + cols = Matrix{eltype(Π.cache.cols)}(undef, Π.F.n_rows, length(col_idcs)) + eval_cols!(cols, Π, u, (), t+dt, col_idcs) + UF = svd(cols).U + @views range_weights = UF/UF[row_idcs,:] + else + @views range_weights = range/range[row_idcs,:] + end + + if size(corange, 2) != length(col_idcs) + @unpack Π = interpolation_cache + rows = Matrix{eltype(Π.cache.cols)}(undef, length(row_idcs), Π.F.n_cols) + eval_rows!(rows, Π, u, (), t+dt, row_idcs) + VF = svd(rows).V + @views corange_weights = VF/VF[col_idcs,:] + else + @views corange_weights = corange/corange[col_idcs,:] + end + return range_weights, corange_weights +end + +function update_interpolation!(interpolation_cache, u, t, dt) + @unpack Π, Π_L, Π_K, Π_S, params = interpolation_cache + @unpack selection_alg = params + + # eval range/corange spans + UF, VF, row_idcs, col_idcs = sparse_selection(interpolation_cache, u, t, dt) + range_weights, corange_weights = compute_weights(UF, VF, row_idcs, col_idcs, u, t, dt, interpolation_cache) + + # new interpolators + range = DEIMInterpolator(row_idcs, range_weights) + projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) + corange = DEIMInterpolator(col_idcs, corange_weights) + projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) + + # update function interpolators + update_interpolator!(Π_L, projected_range) + update_interpolator!(Π_K, projected_corange') + update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) + update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) +end + +# without rank adaptation +#= +function update_interpolation!(interpolation_cache, u, t) + @unpack params, Π, Π_K, Π_L, Π_S = interpolation_cache + @unpack selection_alg, tol, rmin, rmax = params + + eval!(Π_L, u, (), t) # rows + eval!(Π_K, u, (), t) # columns + VF = truncated_svd(Π_L.cache.rows, tol=tol, rmin=rmin, rmax=rmax).V # corange from rows + UF = truncated_svd(Π_K.cache.cols, tol=tol, rmin=rmin, rmax=rmax).U # range from cols + + # find interpolation indices + row_idcs = index_selection(range_svd.U, selection_alg) + col_idcs = index_selection(corange_svd.V, selection_alg) + + # find interpolation weights + @views range_weights = UF/UF[row_idcs,:] + @views corange_weights = VF/VF[col_idcs,:] + + # new interpolators + range = DEIMInterpolator(row_idcs, range_weights) + projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) + corange = DEIMInterpolator(col_idcs, corange_weights) + projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) + + # update function interpolators + update_interpolator!(Π_L, projected_range) + update_interpolator!(Π_K, projected_corange') + update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) + update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) +end + +# alternative but max rank limited by 2*rank(u) +# ... not great +function update_interpolation!(interpolation_cache, u_old, u, dt) + @unpack params, Π, Π_K, Π_L, Π_S = interpolation_cache + @unpack selection_alg, tol, rmin, rmax = params + + dF = u - u_old + dF_svd = svd(dF) + #println(typeof(dF_svd.S)) + cutoff = findfirst(s -> s < tol/dt, diag(dF_svd.S) ) + if isnothing(cutoff) + # compute the full function and take derivative + cutoff= size(dF_svd.S,1) + else + println(dF_svd.S[cutoff]) + cutoff = max(rmin, min(cutoff, rmax)) + end + + if cutoff != rank(Π_K) + println("new DEIM rank: $cutoff") + end + + # find interpolation indices + @views row_idcs = index_selection(dF_svd.U[:,1:cutoff], selection_alg) + @views col_idcs = index_selection(dF_svd.V[:,1:cutoff], selection_alg) + + # find interpolation weights + @views range_weights = dF_svd.U/dF_svd.U[row_idcs,:] + @views corange_weights = dF_svd.V/dF_svd.V[col_idcs,:] + + # new interpolators + range = DEIMInterpolator(row_idcs, range_weights) + projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) + corange = DEIMInterpolator(col_idcs, corange_weights) + projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) + + # update function interpolators + update_interpolator!(Π_L, projected_range) + update_interpolator!(Π_K, projected_corange') + update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) + update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) +end +=# \ No newline at end of file diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index 47ed2dd..764be2e 100644 --- a/src/integrators/unconventional.jl +++ b/src/integrators/unconventional.jl @@ -295,167 +295,6 @@ function unconventional_deim_step!(u, cache, t, dt) update_interpolation!(interpolation_cache, u, t, dt) end -# without rank adaptation -#= -function update_interpolation!(interpolation_cache, u, t) - @unpack params, Π, Π_K, Π_L, Π_S = interpolation_cache - @unpack selection_alg, tol, rmin, rmax = params - - eval!(Π_L, u, (), t) # rows - eval!(Π_K, u, (), t) # columns - VF = truncated_svd(Π_L.cache.rows, tol=tol, rmin=rmin, rmax=rmax).V # corange from rows - UF = truncated_svd(Π_K.cache.cols, tol=tol, rmin=rmin, rmax=rmax).U # range from cols - - # find interpolation indices - row_idcs = index_selection(range_svd.U, selection_alg) - col_idcs = index_selection(corange_svd.V, selection_alg) - - # find interpolation weights - @views range_weights = UF/UF[row_idcs,:] - @views corange_weights = VF/VF[col_idcs,:] - - # new interpolators - range = DEIMInterpolator(row_idcs, range_weights) - projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) - corange = DEIMInterpolator(col_idcs, corange_weights) - projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) - - # update function interpolators - update_interpolator!(Π_L, projected_range) - update_interpolator!(Π_K, projected_corange') - update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) - update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) -end - -# alternative but max rank limited by 2*rank(u) -# ... not great -function update_interpolation!(interpolation_cache, u_old, u, dt) - @unpack params, Π, Π_K, Π_L, Π_S = interpolation_cache - @unpack selection_alg, tol, rmin, rmax = params - - dF = u - u_old - dF_svd = svd(dF) - #println(typeof(dF_svd.S)) - cutoff = findfirst(s -> s < tol/dt, diag(dF_svd.S) ) - if isnothing(cutoff) - # compute the full function and take derivative - cutoff= size(dF_svd.S,1) - else - println(dF_svd.S[cutoff]) - cutoff = max(rmin, min(cutoff, rmax)) - end - - if cutoff != rank(Π_K) - println("new DEIM rank: $cutoff") - end - - # find interpolation indices - @views row_idcs = index_selection(dF_svd.U[:,1:cutoff], selection_alg) - @views col_idcs = index_selection(dF_svd.V[:,1:cutoff], selection_alg) - - # find interpolation weights - @views range_weights = dF_svd.U/dF_svd.U[row_idcs,:] - @views corange_weights = dF_svd.V/dF_svd.V[col_idcs,:] - - # new interpolators - range = DEIMInterpolator(row_idcs, range_weights) - projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) - corange = DEIMInterpolator(col_idcs, corange_weights) - projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) - - # update function interpolators - update_interpolator!(Π_L, projected_range) - update_interpolator!(Π_K, projected_corange') - update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) - update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) -end -=# - -function select_idcs(F_rows::AbstractMatrix, F_cols::AbstractMatrix, selection_alg::IndexSelectionAlgorithm) - range = svd(F_cols) - corange = svd(F_rows) - S = max.(range.S, corange.S) - row_idcs = index_selection(range.U, S, selection_alg) - col_idcs = index_selection(Matrix(corange.V), S, selection_alg) - return range.U, corange.V, row_idcs, col_idcs -end - -function select_idcs(F::SVDLikeRepresentation, selection_alg::IndexSelectionAlgorithm) - F_svd = svd(F) - row_idcs = index_selection(F_svd.U, diag(F_svd.S), selection_alg) - col_idcs = index_selection(F_svd.V, diag(F_svd.S), selection_alg) - UF = F_svd.U[:, 1:length(row_idcs)] - VF = F_trunc.V[:, 1:length(col_idcs)] - return UF, VF, row_idcs, col_idcs -end - -function sparse_selection(interpolation_cache, u, t, dt) - @unpack params = interpolation_cache - if params.update_scheme == :avg_flow - @unpack u_prev = interpolation_cache - @unpack selection_alg = params - return select_idcs(u - u_prev, selection_alg) - elseif params.update_scheme == :last_iterate - @unpack Π_K, Π_L = interpolation_cache - @unpack selection_alg = params - eval!(Π_L, u, (), t+dt) - eval!(Π_K, u, (), t+dt) - UF, VF, row_idcs, col_idcs = select_idcs(Π_L.cache.rows, Π_K.cache.cols, selection_alg) - return UF, VF, row_idcs, col_idcs - else - error("Update scheme $(params.update_scheme) not supported") - end - - # eval index selection - -end - -function compute_weights(range, corange, row_idcs, col_idcs, u, t, dt, interpolation_cache) - @unpack params = interpolation_cache - if size(range, 2) != length(row_idcs) - @unpack Π = interpolation_cache - cols = Matrix{eltype(Π.cache.cols)}(undef, Π.F.n_rows, length(col_idcs)) - eval_cols!(cols, Π, u, (), t+dt, col_idcs) - UF = svd(cols).U - @views range_weights = UF/UF[row_idcs,:] - else - @views range_weights = range/range[row_idcs,:] - end - - if size(corange, 2) != length(col_idcs) - @unpack Π = interpolation_cache - rows = Matrix{eltype(Π.cache.cols)}(undef, length(row_idcs), Π.F.n_cols) - eval_rows!(rows, Π, u, (), t+dt, row_idcs) - VF = svd(rows).V - @views corange_weights = VF/VF[col_idcs,:] - else - @views corange_weights = corange/corange[col_idcs,:] - end - return range_weights, corange_weights -end - -function update_interpolation!(interpolation_cache, u, t, dt) - @unpack Π, Π_L, Π_K, Π_S, params = interpolation_cache - @unpack selection_alg = params - - # eval range/corange spans - UF, VF, row_idcs, col_idcs = sparse_selection(interpolation_cache, u, t, dt) - range_weights, corange_weights = compute_weights(UF, VF, row_idcs, col_idcs, u, t, dt, interpolation_cache) - - - # new interpolators - range = DEIMInterpolator(row_idcs, range_weights) - projected_range = DEIMInterpolator(row_idcs, u.U'*range_weights) - corange = DEIMInterpolator(col_idcs, corange_weights) - projected_corange = DEIMInterpolator(col_idcs, u.V'*corange_weights) - - # update function interpolators - update_interpolator!(Π_L, projected_range) - update_interpolator!(Π_K, projected_corange') - update_interpolator!(Π_S, SparseMatrixInterpolator(projected_range, projected_corange)) - update_interpolator!(Π, SparseMatrixInterpolator(range, corange)) -end - function step!(integrator::DLRIntegrator, ::UnconventionalAlgorithm, dt) @unpack u, t, iter, cache, probType = integrator unconventional_step!(u, cache, t, dt, probType) From 45b105ec4f435ab8fd16759a3afcf2498cafa426 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Mon, 13 Mar 2023 22:50:04 -0400 Subject: [PATCH 28/35] on the fly sparse interpolation for projector splitting algorithms --- src/integrators/projector_splitting.jl | 186 +++++++++++++++++++++-- test/data_agnostic_DEIM_approximation.jl | 12 +- 2 files changed, 185 insertions(+), 13 deletions(-) diff --git a/src/integrators/projector_splitting.jl b/src/integrators/projector_splitting.jl index 0d8b483..71c3a76 100644 --- a/src/integrators/projector_splitting.jl +++ b/src/integrators/projector_splitting.jl @@ -2,7 +2,7 @@ struct PrimalLieTrotter end struct DualLieTrotter end struct Strang end -struct ProjectorSplitting_Params{sType, lType, kType} +struct ProjectorSplitting_Params{sType, lType, kType, iType} S_rhs # rhs of S step (core projected rhs) L_rhs # rhs of L step (range projected rhs) K_rhs # rhs of K step (corange projected rhs) @@ -12,9 +12,12 @@ struct ProjectorSplitting_Params{sType, lType, kType} S_alg::sType L_alg::lType K_alg::kType + interpolation::iType end -struct ProjectorSplitting_Cache{uType,SIntegratorType,LIntegratorType,KIntegratorType,yType} <: AbstractDLRAlgorithm_Cache +struct ProjectorSplitting_Cache{uType, + SIntegratorType,LIntegratorType,KIntegratorType, + yFunc,yType,dType} <: AbstractDLRAlgorithm_Cache US::Matrix{uType} VS::Matrix{uType} QRK::LinearAlgebra.QRCompactWY{uType, Matrix{uType}} @@ -22,10 +25,11 @@ struct ProjectorSplitting_Cache{uType,SIntegratorType,LIntegratorType,KIntegrato SIntegrator::SIntegratorType LIntegrator::LIntegratorType KIntegrator::KIntegratorType - y + y::yFunc ycurr::yType yprev::yType Δy::yType + interpolation_cache::dType end struct ProjectorSplitting{oType, sType, lType, kType} <: AbstractDLRAlgorithm @@ -33,10 +37,13 @@ struct ProjectorSplitting{oType, sType, lType, kType} <: AbstractDLRAlgorithm alg_params::ProjectorSplitting_Params{sType, lType, kType} end -function ProjectorSplitting(order = PrimalLieTrotter();S_rhs = nothing, L_rhs = nothing, K_rhs = nothing, - S_kwargs = Dict(), L_kwargs = Dict(), K_kwargs = Dict(), - S_alg=Tsit5(), L_alg = Tsit5(), K_alg = Tsit5()) - params = ProjectorSplitting_Params(S_rhs,L_rhs,K_rhs,S_kwargs,L_kwargs,K_kwargs,S_alg,L_alg,K_alg) +function ProjectorSplitting(order = PrimalLieTrotter(), deim = nothing;S_rhs = nothing, L_rhs = nothing, K_rhs = nothing, + S_kwargs = Dict(), L_kwargs = Dict(), K_kwargs = Dict(), + S_alg=Tsit5(), L_alg = Tsit5(), K_alg = Tsit5()) + params = ProjectorSplitting_Params(S_rhs,L_rhs,K_rhs, + S_kwargs,L_kwargs,K_kwargs, + S_alg,L_alg,K_alg, + deim) return ProjectorSplitting(order, params) end @@ -81,7 +88,83 @@ function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = p return ProjectorSplitting_Cache(US, VS, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, - nothing, nothing, nothing, nothing) + nothing, nothing, nothing, nothing, nothing) +end + +function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, + alg::ProjectorSplitting, u, dt; t0 = prob.tspan[1]) where + {fType <: ComponentFunction, uType, tType} + tspan = (t0,t0+dt) + + @unpack tol, rmin, rmax, + init_range, init_corange, + selection_alg = alg.alg_params.interpolation + + row_idcs = index_selection(init_range, selection_alg) + col_idcs = index_selection(init_corange, selection_alg) + + Π_corange = DEIMInterpolator(col_idcs, init_corange/init_corange[col_idcs,:]) + Π_range = DEIMInterpolator(row_idcs, init_range/init_range[row_idcs,:]) + Π = SparseFunctionInterpolator(prob.f, SparseMatrixInterpolator(Π_range, Π_corange)) + + US = u.U*u.S + VS = u.V*u.S' + QRK = qr(US) + QRL = qr(VS) + + if isnothing(alg.alg_params.K_rhs) + K_rhs = function (dK, K, p, t) + Π_K, V0, params = p + Π_K(dK, TwoFactorRepresentation(K, V0), params, t) + end + else + K_rhs = alg.alg_params.K_rhs + end + Π_K_corange = DEIMInterpolator(col_idcs, + u.V'*Π.interpolator.corange.weights)' + Π_K = SparseFunctionInterpolator(prob.f, Π_K_corange) + p_K = (Π_K, u.V, ()) + KProblem = ODEProblem(K_rhs, US, tspan, p_K) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + + if isnothing(alg.alg_params.S_rhs) + S_rhs = function (dS, S, p, t) + Π_S, U1, V1, params = p + -Π_S(dS, SVDLikeRepresentation(U1,S,V1), params, t) + end + else + S_rhs = alg.alg_params.S_rhs + end + Π_S_mat = SparseMatrixInterpolator(row_idcs, col_idcs, + u.U'*Π.interpolator.range.weights, + u.V'*Π.interpolator.corange.weights) + Π_S = SparseFunctionInterpolator(prob.f, Π_S_mat) + p_S = (Π_S, u.U, u.V, ()) + + SProblem = ODEProblem(S_rhs, QRK.R, tspan, p_S) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + + if isnothing(alg.alg_params.L_rhs) + L_rhs = function (dL, L, p, t) + Π_L, U0, params = p + Π_L(dL', TwoFactorRepresentation(U0,L), params, t) + end + else + L_rhs = alg.alg_params.L_rhs + end + Π_L_range = DEIMInterpolator(row_idcs, + u.U'*Π.interpolator.range.weights) + Π_L = SparseFunctionInterpolator(prob.f, Π_L_range) + p_L = (Π_L, u.U, ()) + LProblem = ODEProblem(L_rhs, VS, tspan, p_L) + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + + + interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), + Π, Π_K, Π_L, Π_S) + return ProjectorSplitting_Cache(US, VS, QRK, QRL, + SIntegrator, LIntegrator, KIntegrator, + nothing, nothing, nothing, nothing, interpolation_cache) end function alg_cache(prob::MatrixDataProblem, alg::ProjectorSplitting, u, dt; t0 = prob.tspan[1]) @@ -101,7 +184,7 @@ function alg_cache(prob::MatrixDataProblem, alg::ProjectorSplitting, u, dt; t0 = return ProjectorSplitting_Cache(US, VS, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, - prob.y, ycurr, yprev, Δy) + prob.y, ycurr, yprev, Δy, nothing) end function primal_LT_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem}) @@ -178,6 +261,88 @@ function dual_LT_step!(u, cache, t, dt) u.S .= QRK.R end +function primal_LT_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem{<:ComponentFunction}}) + primal_LT_deim_step!(u, cache, t, dt) +end + +function primal_LT_deim_step!(u, cache, t, dt) + @unpack US, VS, QRK, QRL, + KIntegrator, LIntegrator, SIntegrator, + interpolation_cache = cache + @unpack params, u_prev, Π, Π_L, Π_K, Π_S = interpolation_cache + + if params.update_scheme == :avg_flow + u_prev.U = copy(u.U) + u_prev.S = copy(u.S) + u_prev.V = copy(u.V) + end + + # K step + mul!(US,u.U,u.S) + set_u!(KIntegrator, US) + step!(KIntegrator, dt, true) + US .= KIntegrator.u + QRL = qr!(US) + u.U .= Matrix(QRL.Q) + mul!(Π_S.interpolator.range.weights, u.U', Π.interpolator.range.weights) + + # S step + set_u!(SIntegrator, QRL.R) + step!(SIntegrator, dt, true) + + # L step + mul!(VS,u.V,SIntegrator.u') + set_u!(LIntegrator, VS) + step!(LIntegrator, dt, true) + VS .= LIntegrator.u + QRL = qr!(VS) + u.V .= Matrix(QRL.Q) + u.S .= QRL.R' + + update_interpolation!(interpolation_cache, u, t, dt) +end + +function dual_LT_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem{<:ComponentFunction}}) + dual_LT_deim_step!(u, cache, t, dt) +end + +function dual_LT_deim_step!(u, cache, t, dt) + @unpack US, VS, QRK, QRL, + KIntegrator, LIntegrator, SIntegrator, + interpolation_cache = cache + @unpack params, u_prev, Π, Π_L, Π_K, Π_S = interpolation_cache + + if params.update_scheme == :avg_flow + u_prev.U = copy(u.U) + u_prev.S = copy(u.S) + u_prev.V = copy(u.V) + end + + # L step + mul!(VS,u.V,u.S') + set_u!(LIntegrator, VS) + step!(LIntegrator, dt, true) + VS .= LIntegrator.u + QRL = qr!(VS) + u.V .= Matrix(QRL.Q) + mul!(Π_S.interpolator.corange.weights, u.V', Π.interpolator.corange.weights) + + # S step + set_u!(SIntegrator, Matrix(QRL.R')) + step!(SIntegrator, dt, true) + + # K step + mul!(US,u.U,SIntegrator.u) + set_u!(KIntegrator, US) + step!(KIntegrator, dt, true) + US .= KIntegrator.u + QRK = qr!(US) + u.U .= Matrix(QRK.Q) + u.S .= QRK.R + + update_interpolation!(interpolation_cache, u, t, dt) +end + function step!(integrator::DLRIntegrator, ::ProjectorSplitting{PrimalLieTrotter,S,L,K}, dt) where {S,L,K} @unpack u, t, iter, cache, probType = integrator primal_LT_step!(u, cache, t, dt, probType) @@ -198,5 +363,4 @@ function step!(integrator::DLRIntegrator, ::ProjectorSplitting{Strang,S,L,K}, dt dual_LT_step!(u, cache, t + dt/2, dt/2, probType) integrator.t += dt integrator.iter += 1 -end - +end \ No newline at end of file diff --git a/test/data_agnostic_DEIM_approximation.jl b/test/data_agnostic_DEIM_approximation.jl index 68bba09..6df071d 100644 --- a/test/data_agnostic_DEIM_approximation.jl +++ b/test/data_agnostic_DEIM_approximation.jl @@ -227,10 +227,18 @@ using LowRankIntegrators, LinearAlgebra @test norm(dS_test - U'*dX0_lr*V)/norm(Matrix(U'*X0_lr*V)) < 1e-8 end - r_deim = 5 + r_deim = 10 interpolation = SparseInterpolation(DEIM(rmax=r_deim, rmin=r_deim, tol = 1.0, elasticity=1e-1), dX0_lr.U[:,1:r_deim], dX0_lr.V[:,1:r_deim]) alg_deim = UnconventionalAlgorithm(interpolation) deim_prob = MatrixDEProblem(burgers_components, X0_lr, (0.0,1.0)) - deim_sol = LowRankIntegrators.solve(deim_prob, alg_deim, 1e-3, save_increment=10) + deim_sol_ref = LowRankIntegrators.solve(deim_prob, alg_deim, 1e-3, save_increment=10) + + algs = [ProjectorSplitting(PrimalLieTrotter(), interpolation), + ProjectorSplitting(DualLieTrotter(), interpolation), + ProjectorSplitting(Strang(), interpolation)] + for alg in algs + deim_sol = LowRankIntegrators.solve(deim_prob, alg_deim, 1e-3, save_increment=10) + @test norm(Matrix(deim_sol.Y[end] - deim_sol_ref.Y[end]))/norm(Matrix(deim_sol_ref.Y[end])) < 1e-3 + end end \ No newline at end of file From eed9d3877b41d2f8e43b86b41c8b0583cd49c2e0 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Tue, 14 Mar 2023 00:30:00 -0400 Subject: [PATCH 29/35] fix bug in index selection and strengthen test --- src/sparse_interpolation.jl | 8 +++--- test/data_agnostic_DEIM_approximation.jl | 34 ++++++++++++++++-------- 2 files changed, 27 insertions(+), 15 deletions(-) diff --git a/src/sparse_interpolation.jl b/src/sparse_interpolation.jl index 7d6448d..d569dd7 100644 --- a/src/sparse_interpolation.jl +++ b/src/sparse_interpolation.jl @@ -370,13 +370,13 @@ function index_selection(U,::DEIM) r = abs.(U[:, 1]) indices[1] = argmax(r) for l in 2:m - U = U[:, 1:(l - 1)] + Uₗ = U[:, 1:(l - 1)] P = indices[1:(l - 1)] - PᵀU = U[P, :] + PᵀUₗ = Uₗ[P, :] uₗ = U[:, l] Pᵀuₗ = uₗ[P, :] - c = vec(PᵀU \ Pᵀuₗ) - mul!(r, U, c) + c = vec(PᵀUₗ \ Pᵀuₗ) + mul!(r, Uₗ, c) @. r = abs(uₗ - r) indices[l] = argmax(r) end diff --git a/test/data_agnostic_DEIM_approximation.jl b/test/data_agnostic_DEIM_approximation.jl index 6df071d..a1ae91a 100644 --- a/test/data_agnostic_DEIM_approximation.jl +++ b/test/data_agnostic_DEIM_approximation.jl @@ -228,17 +228,29 @@ using LowRankIntegrators, LinearAlgebra end r_deim = 10 - interpolation = SparseInterpolation(DEIM(rmax=r_deim, rmin=r_deim, tol = 1.0, elasticity=1e-1), dX0_lr.U[:,1:r_deim], dX0_lr.V[:,1:r_deim]) - alg_deim = UnconventionalAlgorithm(interpolation) + dt = 1e-4 + sub_alg = SSPRK22() + sub_kwargs = Dict(:dt => dt) + kwargs = Dict(:K_alg => sub_alg, :S_alg => sub_alg, :L_alg => sub_alg, + :K_kwargs => sub_kwargs, :S_kwargs => sub_kwargs, :L_kwargs => sub_kwargs) + idx_selection_algs = [QDEIM(), + DEIM(), + LDEIM(rmax = 2*r_deim, rmin = 5, tol = 1.0, elasticity=1e-1)] + for idx_alg in idx_selection_algs + interpolation = SparseInterpolation(idx_alg, + dX0_lr.U[:,1:r_deim], dX0_lr.V[:,1:r_deim]) + alg_deim = UnconventionalAlgorithm(interpolation; kwargs...) - deim_prob = MatrixDEProblem(burgers_components, X0_lr, (0.0,1.0)) - deim_sol_ref = LowRankIntegrators.solve(deim_prob, alg_deim, 1e-3, save_increment=10) - - algs = [ProjectorSplitting(PrimalLieTrotter(), interpolation), - ProjectorSplitting(DualLieTrotter(), interpolation), - ProjectorSplitting(Strang(), interpolation)] - for alg in algs - deim_sol = LowRankIntegrators.solve(deim_prob, alg_deim, 1e-3, save_increment=10) - @test norm(Matrix(deim_sol.Y[end] - deim_sol_ref.Y[end]))/norm(Matrix(deim_sol_ref.Y[end])) < 1e-3 + deim_prob = MatrixDEProblem(burgers_components, X0_lr, (0.0,1.0)) + deim_sol_ref = LowRankIntegrators.solve(deim_prob, alg_deim, dt, save_increment=100) + + algs = [ProjectorSplitting(PrimalLieTrotter(), interpolation; kwargs...), + ProjectorSplitting(DualLieTrotter(), interpolation; kwargs...), + ProjectorSplitting(Strang(), interpolation; kwargs...)] + for alg in algs + deim_sol = LowRankIntegrators.solve(deim_prob, alg_deim, dt, save_increment=100) + mismatch = norm(Matrix(deim_sol.Y[end] - deim_sol_ref.Y[end]))/norm(Matrix(deim_sol_ref.Y[end])) + @test mismatch < 1e-8 + end end end \ No newline at end of file From bf2005fbb9498bc7025665a82b0263be28ce22d9 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Tue, 14 Mar 2023 01:00:40 -0400 Subject: [PATCH 30/35] disable x86 CI --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e421b31..41ca884 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -24,7 +24,7 @@ jobs: - ubuntu-latest arch: - x64 - - x86 + #- x86 steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 From 15aad94ef55e49917b235a1d9f335c1dd4212439 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Tue, 14 Mar 2023 01:01:01 -0400 Subject: [PATCH 31/35] rm docstring for functor --- src/sparse_interpolation.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/sparse_interpolation.jl b/src/sparse_interpolation.jl index d569dd7..2fc04bf 100644 --- a/src/sparse_interpolation.jl +++ b/src/sparse_interpolation.jl @@ -251,11 +251,6 @@ function eval_cols!(dx, Π::SparseFunctionInterpolator, x, p, t, idcs) Π.F.cols!(dx, x, p, t, idcs) end -""" - $(TYPEDSIGNATURES) - - functor for inplace evaluation of `SparseFunctionInterpolator`s. -""" function (Π::SparseFunctionInterpolator{IP, fType, T})(dx, x, p, t) where {IP <: DEIMInterpolator, fType, T} eval!(Π, x, p, t) mul!(dx, Π.interpolator.weights, Π.cache.rows) From e0eba3de0f7e20356280f08fb6e345dd6a44e5b7 Mon Sep 17 00:00:00 2001 From: FHoltorf <32248677+FHoltorf@users.noreply.github.com> Date: Wed, 26 Apr 2023 23:49:40 -0400 Subject: [PATCH 32/35] type clean up --- .../rank_adaptive_unconventional.jl | 135 +++++++++++++++--- 1 file changed, 118 insertions(+), 17 deletions(-) diff --git a/src/integrators/rank_adaptive_unconventional.jl b/src/integrators/rank_adaptive_unconventional.jl index 905e884..9cbf25a 100644 --- a/src/integrators/rank_adaptive_unconventional.jl +++ b/src/integrators/rank_adaptive_unconventional.jl @@ -22,7 +22,10 @@ function RankAdaptiveUnconventionalAlgorithm(tol; S_rhs = nothing, L_rhs = nothi return RankAdaptiveUnconventionalAlgorithm(params) end -struct RankAdaptiveUnconventionalAlgorithm_Cache{uType,SIntegratorType,LIntegratorType,KIntegratorType,yType} <: AbstractDLRAlgorithm_Cache +struct RankAdaptiveUnconventionalAlgorithm_Cache{uType, + SIntegratorType,LIntegratorType,KIntegratorType, + SProbType, LProbType, KProbType, + yFunc,yType,dType} <: AbstractDLRAlgorithm_Cache US::Matrix{uType} Uhat::Matrix{uType} VS::Matrix{uType} @@ -30,18 +33,19 @@ struct RankAdaptiveUnconventionalAlgorithm_Cache{uType,SIntegratorType,LIntegrat M::Matrix{uType} N::Matrix{uType} SIntegrator::SIntegratorType - SProblem + SProblem::SProbType LIntegrator::LIntegratorType - LProblem + LProblem::LProbType KIntegrator::KIntegratorType - KProblem + KProblem::KProbType r::Int # current rank - tol - r_max - y + tol::Float64 + r_max::Int + y::yFunc ycurr::yType yprev::yType Δy::yType + interpolation_cache::dType end function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) @@ -88,7 +92,95 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, SIntegrator, SProblem, LIntegrator, LProblem, KIntegrator, KProblem, r, alg.alg_params.tol, alg.alg_params.r_max, - nothing, nothing, nothing, nothing) + nothing, nothing, nothing, nothing, nothing) +end + +function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, + alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) where + {fType <: ComponentFunction, uType, tType} + # allocate memory for frequently accessed arrays + tspan = (t0, t0+dt) + r = rank(u) + n, m = size(u) + US = u.U*u.S + VS = u.V*u.S' + Uhat = zeros(n,2*r) + Vhat = zeros(m,2*r) + M = zeros(2*r,r) + N = zeros(2*r,r) + + @unpack tol, rmin, rmax, + init_range, init_corange, + selection_alg = alg.alg_params.interpolation + + row_idcs = index_selection(init_range, selection_alg) + col_idcs = index_selection(init_corange, selection_alg) + + Π_corange = DEIMInterpolator(col_idcs, init_corange/init_corange[col_idcs,:]) + Π_range = DEIMInterpolator(row_idcs, init_range/init_range[row_idcs,:]) + Π = SparseFunctionInterpolator(prob.f, SparseMatrixInterpolator(Π_range, Π_corange)) + + US = u.U*u.S + QRK = qr(US) + M = Matrix(QRK.Q)'*u.U + + VS = u.V*u.S' + QRL = qr(VS) + N = Matrix(QRL.Q)'*u.V + + if isnothing(alg.alg_params.K_rhs) + K_rhs = function (dK, K, p, t) + Π_K, V0, params = p + Π_K(dK, TwoFactorRepresentation(K, V0), params, t) + end + else + K_rhs = alg.alg_params.K_rhs + end + Π_K_corange = DEIMInterpolator(col_idcs, + u.V'*Π.interpolator.corange.weights)' + Π_K = SparseFunctionInterpolator(prob.f, Π_K_corange) + p_K = (Π_K, u.V, ()) + KProblem = ODEProblem(K_rhs, US, tspan, p_K) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + + if isnothing(alg.alg_params.L_rhs) + L_rhs = function (dL, L, p, t) + Π_L, U0, params = p + Π_L(dL', TwoFactorRepresentation(U0,L), params, t) + end + else + L_rhs = alg.alg_params.L_rhs + end + Π_L_range = DEIMInterpolator(row_idcs, + u.U'*Π.interpolator.range.weights) + Π_L = SparseFunctionInterpolator(prob.f, Π_L_range) + p_L = (Π_L, u.U, ()) + LProblem = ODEProblem(L_rhs, VS, tspan, p_L) + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + + if isnothing(alg.alg_params.S_rhs) + S_rhs = function (dS, S, p, t) + Π_S, U1, V1, params = p + Π_S(dS, SVDLikeRepresentation(U1,S,V1), params, t) + end + else + S_rhs = alg.alg_params.S_rhs + end + Π_S_mat = SparseMatrixInterpolator(row_idcs, col_idcs, + U_hat'*Π.interpolator.range.weights, + V_hat'*Π.interpolator.corange.weights) + Π_S = SparseFunctionInterpolator(prob.f, Π_S_mat) + p_S = (Π_S, U_hat, V_hat, ()) + + SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, p_S) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + + interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), + Π, Π_K, Π_L, Π_S) + return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, + SIntegrator, SProblem, LIntegrator, LProblem, + KIntegrator, KProblem, r, alg.alg_params.tol, alg.alg_params.r_max, + nothing, nothing, nothing, nothing, interpolation_cache) end function alg_cache(prob::MatrixDataProblem, alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) @@ -112,14 +204,14 @@ function alg_cache(prob::MatrixDataProblem, alg::RankAdaptiveUnconventionalAlgor LIntegrator = MatrixDataIntegrator(Δy', VS, I, u.U, 1) SIntegrator = MatrixDataIntegrator(Δy, zeros(2r,2r), Uhat, Vhat, 1) - return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, - SIntegrator, nothing, LIntegrator, nothing, - KIntegrator, nothing, r, alg.alg_params.tol, alg.alg_params.r_max, - y, ycurr, yprev, Δy) + return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, + SIntegrator, nothing, LIntegrator, nothing, + KIntegrator, nothing, r, alg.alg_params.tol, alg.alg_params.r_max, + y, ycurr, yprev, Δy, nothing) end function alg_recache(cache::RankAdaptiveUnconventionalAlgorithm_Cache, alg::RankAdaptiveUnconventionalAlgorithm, u, t) - @unpack SProblem, KProblem, LProblem, tol, r_max, y, yprev, ycurr, Δy = cache + @unpack SProblem, KProblem, LProblem, tol, r_max, y, yprev, ycurr, Δy, interpolation_cache = cache r = LowRankArithmetic.rank(u) n, m = size(u) @@ -134,26 +226,35 @@ function alg_recache(cache::RankAdaptiveUnconventionalAlgorithm_Cache, alg::Rank # the following sequence can probably somehow be handled via dispatch in a simpler way if isnothing(KProblem) KIntegrator = MatrixDataIntegrator(Δy, US, I, u.V, 1) - else + elseif isnothing(interpolation_cache) KProblem = remake(KProblem, u0 = US, p = u.V, tspan = (t, t+1.0)) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + else + KProblem = remake(KProblem, u0 = US, p = (interpolation_cache.Π_K, u.V, ()), tspan = (t, t+1.0)) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) end if isnothing(LProblem) LIntegrator = MatrixDataIntegrator(Δy', VS, I, u.U, 1) - else + elseif isnothing(interpolation_cache) LProblem = remake(LProblem, u0 = VS, p = u.U, tspan = (t, t+1.0)) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + else + LProblem = remake(LProblem, u0 = VS, p = (interpolation_cache.Π_L, u.U, ()), tspan = (t, t+1.0)) + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) end if isnothing(SProblem) SIntegrator = MatrixDataIntegrator(Δy, zeros(2r,2r), Uhat, Vhat, 1) - else + elseif isnothing(interpolation_cache) SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = [Uhat, Vhat], tspan = (t, t+1.0)) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + else + SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = (interpolation_cache.Π_S, U_hat, V_hat, ()), tspan = (t, t+1.0)) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) end return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, SIntegrator, SProblem, LIntegrator, LProblem, - KIntegrator, KProblem, r, tol, r_max, y, ycurr, yprev, Δy) + KIntegrator, KProblem, r, tol, r_max, y, ycurr, yprev, Δy, interpolation_cache) end function step!(integrator::DLRIntegrator, alg::RankAdaptiveUnconventionalAlgorithm, dt) From 5410dc627ea4690ad0734d7ec1eabb252e0d0e58 Mon Sep 17 00:00:00 2001 From: FHoltorf Date: Mon, 29 Jan 2024 23:53:23 -0500 Subject: [PATCH 33/35] adjust initialization for correct fsal updates --- src/integrators/on_the_fly_interpolation.jl | 2 +- src/integrators/projector_splitting.jl | 47 +++++++++------ .../rank_adaptive_unconventional.jl | 58 ++++++++++++++----- src/integrators/unconventional.jl | 30 ++++++---- 4 files changed, 93 insertions(+), 44 deletions(-) diff --git a/src/integrators/on_the_fly_interpolation.jl b/src/integrators/on_the_fly_interpolation.jl index df4cad2..7476af1 100644 --- a/src/integrators/on_the_fly_interpolation.jl +++ b/src/integrators/on_the_fly_interpolation.jl @@ -24,7 +24,7 @@ function select_idcs(F::SVDLikeRepresentation, selection_alg::IndexSelectionAlgo row_idcs = index_selection(F_svd.U, diag(F_svd.S), selection_alg) col_idcs = index_selection(F_svd.V, diag(F_svd.S), selection_alg) UF = F_svd.U[:, 1:length(row_idcs)] - VF = F_trunc.V[:, 1:length(col_idcs)] + VF = F_svd.V[:, 1:length(col_idcs)] return UF, VF, row_idcs, col_idcs end diff --git a/src/integrators/projector_splitting.jl b/src/integrators/projector_splitting.jl index 71c3a76..6290750 100644 --- a/src/integrators/projector_splitting.jl +++ b/src/integrators/projector_splitting.jl @@ -50,7 +50,6 @@ end function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = prob.tspan[1]) # allocate memory for frequently accessed arrays tspan = (t0,t0+dt) - US = u.U*u.S VS = u.V*u.S' QRK = qr(US) @@ -64,7 +63,9 @@ function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = p K_rhs = alg.alg_params.K_rhs end KProblem = ODEProblem(K_rhs, US, tspan, u.V) - KIntegrator = init(KProblem, alg.alg_params.K_alg, save_everystep=false, alg.alg_params.K_kwargs...) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) if isnothing(alg.alg_params.S_rhs) S_rhs = function (S, (U,V), t) @@ -74,7 +75,9 @@ function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = p S_rhs = alg.alg_params.S_rhs end SProblem = ODEProblem(S_rhs, QRK.R, tspan, (u.U, u.V)) - SIntegrator = init(SProblem, alg.alg_params.S_alg, save_everystep=false, alg.alg_params.S_kwargs...) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) if isnothing(alg.alg_params.L_rhs) L_rhs = function (VS, U, t) @@ -84,8 +87,10 @@ function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = p L_rhs = alg.alg_params.L_rhs end LProblem = ODEProblem(L_rhs, VS, tspan, u.U) - LIntegrator = init(LProblem, alg.alg_params.L_alg, save_everystep=false, alg.alg_params.L_kwargs...) - + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) + return ProjectorSplitting_Cache(US, VS, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, nothing) @@ -126,11 +131,13 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_K = (Π_K, u.V, ()) KProblem = ODEProblem(K_rhs, US, tspan, p_K) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) - + step!(KIntegrator,0.01*dt,true) + set_t!(KIntegrator, t0) + if isnothing(alg.alg_params.S_rhs) S_rhs = function (dS, S, p, t) Π_S, U1, V1, params = p - -Π_S(dS, SVDLikeRepresentation(U1,S,V1), params, t) + Π_S(dS, SVDLikeRepresentation(U1,S,V1), params, t) end else S_rhs = alg.alg_params.S_rhs @@ -143,7 +150,9 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, SProblem = ODEProblem(S_rhs, QRK.R, tspan, p_S) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) - + step!(SIntegrator,0.01*dt,true) + set_t!(SIntegrator, t0) + if isnothing(alg.alg_params.L_rhs) L_rhs = function (dL, L, p, t) Π_L, U0, params = p @@ -158,10 +167,12 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_L = (Π_L, u.U, ()) LProblem = ODEProblem(L_rhs, VS, tspan, p_L) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) - + step!(LIntegrator,0.01*dt,true) + set_t!(LIntegrator, t0) + interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), - Π, Π_K, Π_L, Π_S) + Π, Π_K, Π_L, Π_S) return ProjectorSplitting_Cache(US, VS, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, interpolation_cache) @@ -282,15 +293,16 @@ function primal_LT_deim_step!(u, cache, t, dt) set_u!(KIntegrator, US) step!(KIntegrator, dt, true) US .= KIntegrator.u - QRL = qr!(US) - u.U .= Matrix(QRL.Q) - mul!(Π_S.interpolator.range.weights, u.U', Π.interpolator.range.weights) - + QRK = qr!(US) + u.U .= Matrix(QRK.Q) + # S step - set_u!(SIntegrator, QRL.R) + mul!(Π_S.interpolator.range.weights, u.U', Π.interpolator.range.weights, -1.0, 0) + set_u!(SIntegrator, QRK.R) step!(SIntegrator, dt, true) # L step + mul!(Π_L.interpolator.weights, u.U', Π.interpolator.range.weights) mul!(VS,u.V,SIntegrator.u') set_u!(LIntegrator, VS) step!(LIntegrator, dt, true) @@ -325,14 +337,15 @@ function dual_LT_deim_step!(u, cache, t, dt) VS .= LIntegrator.u QRL = qr!(VS) u.V .= Matrix(QRL.Q) - mul!(Π_S.interpolator.corange.weights, u.V', Π.interpolator.corange.weights) - + # S step + mul!(Π_S.interpolator.corange.weights, u.V', Π.interpolator.corange.weights, -1, 0) set_u!(SIntegrator, Matrix(QRL.R')) step!(SIntegrator, dt, true) # K step mul!(US,u.U,SIntegrator.u) + mul!(Π_K.interpolator.parent.weights, u.V', Π.interpolator.corange.weights) set_u!(KIntegrator, US) step!(KIntegrator, dt, true) US .= KIntegrator.u diff --git a/src/integrators/rank_adaptive_unconventional.jl b/src/integrators/rank_adaptive_unconventional.jl index 9cbf25a..6539954 100644 --- a/src/integrators/rank_adaptive_unconventional.jl +++ b/src/integrators/rank_adaptive_unconventional.jl @@ -69,6 +69,9 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit end KProblem = ODEProblem(K_rhs, US, tspan, u.V) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) + if isnothing(alg.alg_params.L_rhs) L_rhs = function (VS, U, t) return Matrix(prob.f(TwoFactorRepresentation(U,VS),t)'*U) @@ -76,9 +79,11 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit else L_rhs = alg.alg_params.L_rhs end - LProblem = ODEProblem(L_rhs, VS, tspan, u.U) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) + if isnothing(alg.alg_params.S_rhs) S_rhs = function (S, (U,V), t) return Matrix(U'*prob.f(SVDLikeRepresentation(U,S,V),t)*V) @@ -88,6 +93,8 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit end SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, [Uhat, Vhat]) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, SIntegrator, SProblem, LIntegrator, LProblem, @@ -142,6 +149,8 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_K = (Π_K, u.V, ()) KProblem = ODEProblem(K_rhs, US, tspan, p_K) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) if isnothing(alg.alg_params.L_rhs) L_rhs = function (dL, L, p, t) @@ -157,6 +166,8 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_L = (Π_L, u.U, ()) LProblem = ODEProblem(L_rhs, VS, tspan, p_L) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) if isnothing(alg.alg_params.S_rhs) S_rhs = function (dS, S, p, t) @@ -174,6 +185,8 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, p_S) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), Π, Π_K, Π_L, Π_S) @@ -226,31 +239,44 @@ function alg_recache(cache::RankAdaptiveUnconventionalAlgorithm_Cache, alg::Rank # the following sequence can probably somehow be handled via dispatch in a simpler way if isnothing(KProblem) KIntegrator = MatrixDataIntegrator(Δy, US, I, u.V, 1) - elseif isnothing(interpolation_cache) - KProblem = remake(KProblem, u0 = US, p = u.V, tspan = (t, t+1.0)) - KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) else - KProblem = remake(KProblem, u0 = US, p = (interpolation_cache.Π_K, u.V, ()), tspan = (t, t+1.0)) - KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + if isnothing(interpolation_cache) + KProblem = remake(KProblem, u0 = US, p = u.V, tspan = (t, t+1.0)) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + else + KProblem = remake(KProblem, u0 = US, p = (interpolation_cache.Π_K, u.V, ()), tspan = (t, t+1.0)) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + end + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) end + if isnothing(LProblem) LIntegrator = MatrixDataIntegrator(Δy', VS, I, u.U, 1) - elseif isnothing(interpolation_cache) - LProblem = remake(LProblem, u0 = VS, p = u.U, tspan = (t, t+1.0)) - LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) else - LProblem = remake(LProblem, u0 = VS, p = (interpolation_cache.Π_L, u.U, ()), tspan = (t, t+1.0)) - LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + if isnothing(interpolation_cache) + LProblem = remake(LProblem, u0 = VS, p = u.U, tspan = (t, t+1.0)) + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + else + LProblem = remake(LProblem, u0 = VS, p = (interpolation_cache.Π_L, u.U, ()), tspan = (t, t+1.0)) + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + end + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) end if isnothing(SProblem) SIntegrator = MatrixDataIntegrator(Δy, zeros(2r,2r), Uhat, Vhat, 1) - elseif isnothing(interpolation_cache) - SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = [Uhat, Vhat], tspan = (t, t+1.0)) - SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) else - SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = (interpolation_cache.Π_S, U_hat, V_hat, ()), tspan = (t, t+1.0)) - SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + if isnothing(interpolation_cache) + SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = [Uhat, Vhat], tspan = (t, t+1.0)) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + else + SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = (interpolation_cache.Π_S, U_hat, V_hat, ()), tspan = (t, t+1.0)) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + end + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) end return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, SIntegrator, SProblem, LIntegrator, LProblem, diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index 764be2e..cd9c687 100644 --- a/src/integrators/unconventional.jl +++ b/src/integrators/unconventional.jl @@ -80,7 +80,9 @@ function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t end KProblem = ODEProblem(K_rhs, US, tspan, u.V) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) - + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) + if isnothing(alg.alg_params.L_rhs) L_rhs = function (VS, U, t) return Matrix(prob.f(TwoFactorRepresentation(U,VS),t)'*U) @@ -90,7 +92,9 @@ function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t end LProblem = ODEProblem(L_rhs, VS, tspan, u.U) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) - + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) + if isnothing(alg.alg_params.S_rhs) S_rhs = function (S, (U,V), t) return Matrix(U'*prob.f(SVDLikeRepresentation(U,S,V),t)*V) @@ -100,17 +104,17 @@ function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t end SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, (u.U, u.V)) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) - + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) + return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, nothing) end - function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) where {fType <: ComponentFunction, uType, tType} - # allocate memory for frequently accessed arrays tspan = (t0,t0+dt) @unpack tol, rmin, rmax, @@ -146,7 +150,9 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_K = (Π_K, u.V, ()) KProblem = ODEProblem(K_rhs, US, tspan, p_K) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) - + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) + if isnothing(alg.alg_params.L_rhs) L_rhs = function (dL, L, p, t) Π_L, U0, params = p @@ -161,7 +167,9 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_L = (Π_L, u.U, ()) LProblem = ODEProblem(L_rhs, VS, tspan, p_L) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) - + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) + if isnothing(alg.alg_params.S_rhs) S_rhs = function (dS, S, p, t) Π_S, U1, V1, params = p @@ -178,9 +186,11 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, p_S) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) - + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) + interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), - Π, Π_K, Π_L, Π_S) + Π, Π_K, Π_L, Π_S) return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, interpolation_cache) @@ -282,7 +292,7 @@ function unconventional_deim_step!(u, cache, t, dt) mul!(N,Matrix(QRL.Q)',u.V) u.V .= Matrix(QRL.Q) u.U .= Matrix(QRK.Q) - + # update core interpolator mul!(Π_S.interpolator.range.weights, u.U', Π.interpolator.range.weights) mul!(Π_S.interpolator.corange.weights, u.V', Π.interpolator.corange.weights) From 242db61faf8d07b77520f8e2e28e1c525c2a01c3 Mon Sep 17 00:00:00 2001 From: FHoltorf Date: Tue, 30 Jan 2024 23:00:56 -0500 Subject: [PATCH 34/35] add progressmeter --- Project.toml | 1 + src/LowRankIntegrators.jl | 2 +- src/primitives.jl | 5 ++--- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index d09dd9f..2139b9b 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LowRankArithmetic = "8016d1d1-7a3b-4ad2-b79e-5984174ed314" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" diff --git a/src/LowRankIntegrators.jl b/src/LowRankIntegrators.jl index 04558a1..928b486 100644 --- a/src/LowRankIntegrators.jl +++ b/src/LowRankIntegrators.jl @@ -1,5 +1,5 @@ module LowRankIntegrators -using Reexport, UnPack, DocStringExtensions, LinearAlgebra +using Reexport, UnPack, DocStringExtensions, LinearAlgebra, ProgressMeter import DifferentialEquations: step!, set_u!, init @reexport using DifferentialEquations, LowRankArithmetic diff --git a/src/primitives.jl b/src/primitives.jl index a198f5e..09b61a2 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -72,15 +72,14 @@ function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt; save_inc println("Initialize integrator ...") integrator, dt, t_int, save_idcs = init(prob, alg, dt, save_increment) println("... initialization complete. Start Integration ...") - disp_digits = -floor(Int, log10(save_increment*dt)) k = 2 - #while (prob.tspan[2]-integrator.t)/T > 1e-8 + progressmeter = Progress(length(t_int)) while integrator.iter < length(t_int) - 1 step!(integrator, alg, dt) + next!(progressmeter) if integrator.iter + 1 == save_idcs[k] update_sol!(integrator, k) k += 1 - println("t = $(round(integrator.t,digits=disp_digits))") end end println("... integration complete.") From 1cb2da019e28c6ede529a7e90eec4c5e78774f27 Mon Sep 17 00:00:00 2001 From: FHoltorf Date: Tue, 30 Jan 2024 23:14:43 -0500 Subject: [PATCH 35/35] flush progressmeter to stream --- src/primitives.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/primitives.jl b/src/primitives.jl index 09b61a2..7024e3c 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -77,6 +77,7 @@ function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt; save_inc while integrator.iter < length(t_int) - 1 step!(integrator, alg, dt) next!(progressmeter) + flush(stdout) if integrator.iter + 1 == save_idcs[k] update_sol!(integrator, k) k += 1