diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f86d22a..41ca884 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -18,13 +18,13 @@ jobs: fail-fast: false matrix: version: - - '1.6' + - '1.8' #- 'nightly' avoid for now os: - ubuntu-latest arch: - x64 - - x86 + #- x86 steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 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 + diff --git a/Project.toml b/Project.toml index f441bef..2139b9b 100644 --- a/Project.toml +++ b/Project.toml @@ -5,17 +5,18 @@ 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" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" 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" diff --git a/src/LowRankIntegrators.jl b/src/LowRankIntegrators.jl index d31a6ea..928b486 100644 --- a/src/LowRankIntegrators.jl +++ b/src/LowRankIntegrators.jl @@ -1,12 +1,20 @@ module LowRankIntegrators -using Reexport, UnPack +using Reexport, UnPack, DocStringExtensions, LinearAlgebra, ProgressMeter 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 +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..2db1c0d 100644 --- a/src/integrators.jl +++ b/src/integrators.jl @@ -1,5 +1,8 @@ +# basic catch all cases: +include("integrators/on_the_fly_interpolation.jl") 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/greedy_integrator.jl b/src/integrators/greedy_integrator.jl index 8b0dae1..c78b3ee 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 @@ -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) @@ -46,18 +46,6 @@ 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::MatrixHybridProblem, alg::GreedyIntegrator, dt) t0, tf = prob.tspan @assert tf > t0 "Integration in reverse time direction is not supported" diff --git a/src/integrators/on_the_fly_interpolation.jl b/src/integrators/on_the_fly_interpolation.jl new file mode 100644 index 0000000..7476af1 --- /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_svd.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/projector_splitting.jl b/src/integrators/projector_splitting.jl index 7459519..6290750 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,17 +37,19 @@ 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 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) @@ -57,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) @@ -67,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) @@ -77,11 +87,95 @@ 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) +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...) + 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) + 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...) + 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 + Π_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...) + step!(LIntegrator,0.01*dt,true) + set_t!(LIntegrator, t0) + + 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) + nothing, nothing, nothing, nothing, interpolation_cache) end function alg_cache(prob::MatrixDataProblem, alg::ProjectorSplitting, u, dt; t0 = prob.tspan[1]) @@ -101,17 +195,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) -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) + prob.y, ycurr, yprev, Δy, nothing) end function primal_LT_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem}) @@ -188,6 +272,90 @@ 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 + QRK = qr!(US) + u.U .= Matrix(QRK.Q) + + # S step + 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) + 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) + + # 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 + 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) @@ -208,5 +376,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/src/integrators/rank_adaptive_unconventional.jl b/src/integrators/rank_adaptive_unconventional.jl index 13ecfe0..6539954 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]) @@ -65,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) @@ -72,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) @@ -84,23 +93,107 @@ 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, KIntegrator, KProblem, r, alg.alg_params.tol, alg.alg_params.r_max, - nothing, nothing, nothing, nothing) + nothing, 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) +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...) + 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 + Π_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...) + 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 + Π_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...) + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) + + 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]) @@ -124,15 +217,15 @@ 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 - r = rank(u) + @unpack SProblem, KProblem, LProblem, tol, r_max, y, yprev, ycurr, Δy, interpolation_cache = cache + r = LowRankArithmetic.rank(u) n, m = size(u) US = zeros(n,r) @@ -147,25 +240,47 @@ function alg_recache(cache::RankAdaptiveUnconventionalAlgorithm_Cache, alg::Rank if isnothing(KProblem) KIntegrator = MatrixDataIntegrator(Δy, US, I, u.V, 1) else - 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...) + 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) else - 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...) + 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) else - 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...) + 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, - 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) @@ -220,7 +335,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])) diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index b1283d7..cd9c687 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,25 @@ 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 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,12 +45,20 @@ 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 +rank_DEIM(cache::UnconventionalAlgorithm_Cache) = rank_DEIM(cache.interpolation_cache) +rank_DEIM(::Nothing) = 0 +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::OnTheFlyInterpolation_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) @@ -57,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) @@ -67,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) @@ -77,13 +104,100 @@ 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} + 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 + 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...) + 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 + Π_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...) + 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 + Π_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, 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) return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, - nothing, nothing, nothing, nothing) + nothing, nothing, nothing, nothing, interpolation_cache) end -function alg_cache(prob::MatrixDataProblem, alg::UnconventionalAlgorithm,u,dt; t0 = prob.tspan[1]) + +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,19 +217,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) -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) + prob.y, ycurr, yprev, Δy, nothing) end function unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem}) @@ -130,9 +232,49 @@ 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 + @unpack US, VS, M, N, QRK, QRL, + KIntegrator, LIntegrator, SIntegrator = 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) + set_u!(SIntegrator, M*u.S*N') + step!(SIntegrator, dt, true) + 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 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) @@ -150,10 +292,17 @@ function unconventional_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) + # integration set_u!(SIntegrator, M*u.S*N') step!(SIntegrator, dt, true) u.S .= SIntegrator.u + + update_interpolation!(interpolation_cache, u, t, dt) end function step!(integrator::DLRIntegrator, ::UnconventionalAlgorithm, dt) diff --git a/src/primitives.jl b/src/primitives.jl index 3d5ae2d..7024e3c 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,13 +68,22 @@ end """ solves the given problem with the specified algorithm and step size """ -function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt) - integrator = init(prob, alg, dt) - 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, dt, t_int, save_idcs = init(prob, alg, dt, save_increment) + println("... initialization complete. Start Integration ...") + k = 2 + progressmeter = Progress(length(t_int)) + while integrator.iter < length(t_int) - 1 step!(integrator, alg, dt) - update_sol!(integrator) + next!(progressmeter) + flush(stdout) + if integrator.iter + 1 == save_idcs[k] + update_sol!(integrator, k) + k += 1 + end end + println("... integration complete.") return integrator.sol end function solve(prob::MatrixDataProblem, alg::AbstractDLRAlgorithm) @@ -79,26 +91,60 @@ 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, 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), dt, 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(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 -function init_sol(dt::Int, t0, tf, u0) +function init_sol(dt::Int, t0::Int, tf::Int, u0, save_increment::Int) # initialize solution object - steps = t0:dt:tf - return DLRSolution(Vector{typeof(u0)}(undef, length(steps)), collect(steps)) + t_int = t0:dt:tf + 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 + 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) + 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 + +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) end \ No newline at end of file diff --git a/src/sparse_interpolation.jl b/src/sparse_interpolation.jl new file mode 100644 index 0000000..2fc04bf --- /dev/null +++ b/src/sparse_interpolation.jl @@ -0,0 +1,459 @@ +import LinearAlgebra: adjoint, rank +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 + +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) + + 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) + + 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) + + 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 = InterpolatorCache(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 = 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.weights, 1) + dim_corange = size(interpolator.corange.weights, 1) + T = eltype(Π.cache.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 + +""" + $(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 + +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 + +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 <: 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 <: 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,::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,::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 = findlast(x -> x > alg.tol*alg.elasticity, S) + if isnothing(cutoff) + 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 \ No newline at end of file diff --git a/test/data_agnostic_DEIM_approximation.jl b/test/data_agnostic_DEIM_approximation.jl new file mode 100644 index 0000000..a1ae91a --- /dev/null +++ b/test/data_agnostic_DEIM_approximation.jl @@ -0,0 +1,256 @@ +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(init_range, alg) + col_idcs = index_selection(init_corange, 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 = 10 + 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, 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 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) diff --git a/test/runtests.jl b/test/runtests.jl index 0475d8f..8fda4f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,8 @@ -using LowRankIntegrators +using LinearAlgebra, 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") +include("data_agnostic_DEIM_approximation.jl") +include("sparse_interpolation.jl") diff --git a/test/sparse_interpolation.jl b/test/sparse_interpolation.jl new file mode 100644 index 0000000..da0607f --- /dev/null +++ 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