diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e9eba34..86d6a33 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -12,7 +12,7 @@ jobs: - "1" - "1.6" steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} diff --git a/Project.toml b/Project.toml index 6686110..8b0f863 100644 --- a/Project.toml +++ b/Project.toml @@ -4,22 +4,26 @@ authors = ["Bowen S. Zhu and contributors"] version = "0.1.1" [deps] +Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +PolyChaos = "8d666b04-775d-5f6e-b778-5ac7c70f65a3" RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e" [compat] +Combinatorics = "1.0" +Distributions = "0.25" DocStringExtensions = "0.8, 0.9" ModelingToolkit = "8.21" +PolyChaos = "0.2" RandomizedLinAlg = "0.1" Setfield = "0.8, 1" -SymbolicUtils = "0.19" Symbolics = "4.10" TSVD = "0.4" julia = "1.6" diff --git a/docs/Project.toml b/docs/Project.toml index af1c16a..36ebf64 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -6,6 +6,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PolyChaos = "8d666b04-775d-5f6e-b778-5ac7c70f65a3" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] Documenter = "0.27" diff --git a/docs/make.jl b/docs/make.jl index 0e20c36..0f797da 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,4 +1,4 @@ -using Documenter, ModelOrderReduction +using Documenter, ModelOrderReduction, Symbolics, PolyChaos include("pages.jl") diff --git a/docs/pages.jl b/docs/pages.jl index 2a7c950..451d264 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,5 +1,12 @@ pages = [ "Home" => "index.md", - "Functions" => "functions.md", - "Tutorials" => Any["tutorials/deim_FitzHugh-Nagumo.md"], + "Polynomial Chaos Expansion (PCE)" => [ + "Introduction" => "pce/pce.md", + "pce/stochastic_galerkin.md", + ], + "Discrete Empirical Interpolation Method (DEIM)" => [ + "Introduction" => "deim/deim.md", + "deim/deim_FitzHugh-Nagumo.md", + ], + "internals.md", ] diff --git a/docs/src/deim/deim.md b/docs/src/deim/deim.md new file mode 100644 index 0000000..4c45eaa --- /dev/null +++ b/docs/src/deim/deim.md @@ -0,0 +1,5 @@ +# Discrete Empirical Interpolation Method (DEIM) + +```@docs +deim +``` diff --git a/docs/src/tutorials/deim_FitzHugh-Nagumo.md b/docs/src/deim/deim_FitzHugh-Nagumo.md similarity index 99% rename from docs/src/tutorials/deim_FitzHugh-Nagumo.md rename to docs/src/deim/deim_FitzHugh-Nagumo.md index 649a592..51487e6 100644 --- a/docs/src/tutorials/deim_FitzHugh-Nagumo.md +++ b/docs/src/deim/deim_FitzHugh-Nagumo.md @@ -1,4 +1,4 @@ -# Discrete Empirical Interpolation Method (DEIM) +# Example: FitzHugh-Nagumo System This section illustrates how ModelOrderReduction.jl can be used to build a reduced order model via Petrov-Galerkin projection using the Proper Orthogonal Decomposition (POD) and diff --git a/docs/src/functions.md b/docs/src/functions.md deleted file mode 100644 index 0aa1de4..0000000 --- a/docs/src/functions.md +++ /dev/null @@ -1,6 +0,0 @@ -```@index -``` - -```@docs -deim -``` diff --git a/docs/src/internals.md b/docs/src/internals.md new file mode 100644 index 0000000..d35ad2e --- /dev/null +++ b/docs/src/internals.md @@ -0,0 +1,6 @@ +# Internal Documentation + +```@autodocs +Modules = [ModelOrderReduction] +Public = false +``` \ No newline at end of file diff --git a/docs/src/pce/pce.md b/docs/src/pce/pce.md new file mode 100644 index 0000000..7cdcb5f --- /dev/null +++ b/docs/src/pce/pce.md @@ -0,0 +1,6 @@ +# Polynomial Chaos Expansion (PCE) + +```@docs +PCE +PCE(::AbstractVector{Num}, ::AbstractVector{Num}, ::AbstractVector{Num}, ::AbstractVector{<:Union{AbstractOrthoPoly, AbstractCanonicalOrthoPoly}}) +``` diff --git a/docs/src/pce/stochastic_galerkin.md b/docs/src/pce/stochastic_galerkin.md new file mode 100644 index 0000000..bf271f6 --- /dev/null +++ b/docs/src/pce/stochastic_galerkin.md @@ -0,0 +1,5 @@ +# Stochastic Galerkin Method (SGM) + +```@docs +stochastic_galerkin +``` diff --git a/examples/PCE/mean.png b/examples/PCE/mean.png new file mode 100644 index 0000000..b2337a3 Binary files /dev/null and b/examples/PCE/mean.png differ diff --git a/examples/PCE/traces.png b/examples/PCE/traces.png new file mode 100644 index 0000000..4bc5295 Binary files /dev/null and b/examples/PCE/traces.png differ diff --git a/examples/PCE/tutorial.jl b/examples/PCE/tutorial.jl new file mode 100644 index 0000000..93eff3a --- /dev/null +++ b/examples/PCE/tutorial.jl @@ -0,0 +1,123 @@ +using ModelOrderReduction, ModelingToolkit, DifferentialEquations, + PolyChaos, Plots, LinearAlgebra + +# Reaction system +# k[1](1+θ[1]/2)/k[2]: A + B ⇌ C +# k[3](1+θ[2]/2): C → D +# k[4]: B → D +# +# state: c[1:4] = [A, B, C, D] +# certain parameters: k[1:4] +# uncertain parameters: θ[1:2] ∼ U[-1,1] + +@parameters k[1:4], θ[1:2] +@variables t, c(t)[1:4] +D = Differential(t) +reactor = [D(c[1]) ~ -k[1] * (1 + 0.5 * θ[1]) * c[1] * c[2] + k[2] * c[3]; + D(c[2]) ~ -k[1] * (1 + 0.5 * θ[1]) * c[1] * c[2] - k[4] * c[2] + k[2] * c[3]; + D(c[3]) ~ k[1] * c[1] * c[2] - k[3] * c[3]; + D(c[4]) ~ k[3] * (1 + 0.5 * θ[2]) * c[3] + k[4] * c[2]]; + +@named reactor_model = ODESystem(reactor, t, c, vcat(k, θ)) + +# 1. choose/generate PCE of the system state +d_pce = 3 +pce = PCE(c, [θ[i] => Uniform_11OrthoPoly(d_pce) for i in eachindex(θ)]) + +# 2. generate moment equations +moment_eqs, pce_eval = moment_equations(reactor_model, pce) + +# 3. solve the moment equations +# 3a. find initial moments via Galerkin projection (note that c0 could depend on an uncertain parameter) +c0 = [1.0, 2.0, 0.0, 0.0] +z0 = reduce(vcat, pce_galerkin(c0, pce)) + +# 3b. solve the initial value problem for the moment equations +moment_prob = ODEProblem(moment_eqs, z0, (0.0, 10.0), + [k[1] => 1.0, k[2] => 0.2, k[3] => 12, k[4] => 0.1]) +moment_sol = solve(moment_prob, Tsit5()) + +# 3c. define a function approximating the parametric solution to the IVP. +pce_sol(t, θ) = pce_eval(moment_sol(t), θ) + +# 4. now let's compare to true solution +# take 100 uniform samples +n = 5 +θsamples = Iterators.product(range(-1, 1, length = n), range(-1, 1, length = n)) +ps = [k[1] => 1.0, k[2] => 0.2, k[3] => 12, k[4] => 0.1, θ[1] => 0.0, θ[2] => 0.0] +pvals = [p[2] for p in ps] +reactor_problem = ODEProblem(reactor_model, c0, (0.0, 10.0), ps) +t_range = range(0.0, 10.0, length = 100) +sols = [] +for θval in θsamples + pvals[end - 1] = θval[1] + pvals[end] = θval[2] + _prob = remake(reactor_problem, p = pvals) + push!(sols, Array(solve(_prob, Tsit5(), saveat = t_range))) +end + +species = ["A", "B", "C", "D"] +color = [:red, :green, :blue, :orange] +plots = [plot(title = "Species $(species[i])", xlabel = "time", ylabel = "concentration", + legend = false) for i in 1:4] +for sol in sols + for i in 1:4 + plot!(plots[i], t_range, sol[i, :], color = color[i]) + end +end +for θval in θsamples + pce_predictions = [pce_sol(t, [θval...]) for t in t_range] + for i in 1:4 + plot!(plots[i], t_range, [c[i] for c in pce_predictions], color = "black", + linestyle = :dash) + end +end + +fig = plot(plots...) +savefig(fig, string(@__DIR__, "/traces.png")) + +n = 100 +θsamples = Iterators.product(range(-1, 1, length = n), range(-1, 1, length = n)) +ps = [k[1] => 1.0, k[2] => 0.2, k[3] => 12, k[4] => 0.1, θ[1] => 0.0, θ[2] => 0.0] +pvals = [p[2] for p in ps] +reactor_problem = ODEProblem(reactor_model, c0, (0.0, 10.0), ps) +t_range = range(0.0, 10.0, length = 100) +sols = [] +for θval in θsamples + pvals[end - 1] = θval[1] + pvals[end] = θval[2] + _prob = remake(reactor_problem, p = pvals) + push!(sols, Array(solve(_prob, Tsit5(), saveat = t_range))) +end + +mean_solution = mean(sols) +var_solution = mean([sol .^ 2 for sol in sols]) - mean_solution .^ 2 + +L = length(pce.sym_basis) +mean_PCE = [[moment_sol(t)[((i - 1) * L + 1):(i * L)][1] for i in 1:4] for t in t_range] +var_weightings = computeSP2(pce.tensor_basis) +var_PCE = [[dot(moment_sol(t)[((i - 1) * L + 1):(i * L)] .^ 2, var_weightings) for i in 1:4] + for t in t_range] .- [m .^ 2 for m in mean_PCE] + +species = ["A", "B", "C", "D"] +color = [:red, :green, :blue, :orange] +plots = [plot(title = "Species $(species[i])", xlabel = "time", ylabel = "mean", + legend = false) for i in 1:4] + +for i in 1:4 + plot!(plots[i], t_range, [mean_solution[i, k] for k in axes(mean_solution, 2)], + color = color[i]) + plot!(plots[i], t_range, [s[i] for s in mean_PCE], color = "black", linestyle = :dash) +end +fig = plot(plots...) +savefig(fig, string(@__DIR__, "/mean.png")) + +plots = [plot(title = "Species $(species[i])", xlabel = "time", ylabel = "variance", + legend = false) for i in 1:4] +for i in 1:4 + plot!(plots[i], t_range, [var_solution[i, k] for k in axes(var_solution, 2)], + color = color[i]) + plot!(plots[i], t_range, [s[i] for s in var_PCE], color = "black", linestyle = :dash) +end +fig = plot(plots...) +savefig(fig, string(@__DIR__, "/var.png")) diff --git a/examples/PCE/tutorial.md b/examples/PCE/tutorial.md new file mode 100644 index 0000000..8c1dc36 --- /dev/null +++ b/examples/PCE/tutorial.md @@ -0,0 +1,109 @@ +# Polynomial Chaos Expansion for Random ODEs +## Sketch of the Theory +This section reviews the basic concept behind the use of Polynomial Chaos Expansion (PCE) for approximate solution to random ODEs as well as how `ModelOrderReduction` facilitates its use within the SciML ecosystem. + +In essence, PCE is nothing more than a function expansion technique for square integrable functions where the basis functions are chosen to be orthogonal polynomials; for a detailed treatment of the mathematical underpinnings the interested reader is referred to Gautschi (2004) [1], for example. For the sake of this tutorial we will remain much more informal and discuss only aspects that are critical for basic use of PCE in computational practice. To that end, we will consider the problem of finding a parametric solution to the following *random* ODE: +```math +\begin{aligned} + \frac{dx}{dt}(p) = f(x,p,t), \quad x(0) = x_0(p) +\end{aligned} +``` +where the parameters $p$ are *apriori* only known to follow a given distribution $\mu$. Specifically, we will seek to find an approximation to the solution $x(t,p)$ in the form of a PCE that matches the statistics of the random process $x(t,p)$ with $p \sim \mu$. To that end, we approximate +```math +\begin{aligned} + x(t,p) \approx \hat{x}(t,p) = \sum_{i=0}^L z_i(t) \zeta_i(p) +\end{aligned} +``` +where $\zeta_0(p), \dots, \zeta_L(p)$ are mutually orthogonal polynomial with respect to the distribution $\mu$, i.e., +```math +\begin{aligned} + \left\langle \zeta_i, \zeta_j \right\rangle := \int \zeta_i(p) \zeta_j(p) d\mu(p) = \delta_{ij} +\end{aligned} +``` +where $\delta_{ij}$ is the Kronecker delta. These polynomials can be efficiently constructed for many distributions. In particular, they are known analytically for many standard distributions including but not limited to Gaussians, uniform or beta distributions; moreover, they can be approximated numerically for distributions with known smooth densities and even empirically known distributions. For more details, please review the documentation of [PolyChaos.jl](github.com/SciML/PolyChaos.jl) or Gautschi (2004) [1] for instance. + +Inserting the PCE-Ansatz $\hat{x}$ into the random ODE of interest yields +```math +\begin{aligned} + \sum_{i=0}^L \frac{dz_i}{dt}(t) \zeta_i(p) = f\left(\sum_{i=0}^L z_i(t) \zeta_i(p), p, t \right), \sum_{i=0}^L z_i(0) \zeta_i(p) +\end{aligned} +``` +To derive evolution equations for the coefficients (or *moments*) $z_i$ we simply demand that the residual of the above equations to be orthogonal to the space spanned by $\zeta_0, \dots, \zeta_L$. Orthogonality here shall be understood with respect to the inner product as described above. This is equivalent to demanding minimal expected mismatch when evaluating the above equation for $p \sim \mu$. Mathematically, this condition translates into a set of deterministic ODEs for the moments $z_j$: +```math +\begin{aligned} + \left\langle\sum_{i=0}^L \frac{dz_i}{dt}(t) \zeta_i(p), \zeta_j(p) \right\rangle = \left \langle f\left(\sum_{i=0}^L z_i(t) \zeta_i(p), p, t \right), \sum_{i=0}^L z_i(0) \zeta_i(p), \zeta_j(p) \right\rangle +\end{aligned} +``` +and since $\left\langle \zeta_i, \zeta_j\right\rangle = \delta_{ij}$ it follows that +```math +\begin{aligned} + \frac{dz_0}{dt} &= \left \langle f\left(\sum_{i=0}^L z_i(t) \zeta_i(p), p, t \right), \sum_{i=0}^L z_i(0) \zeta_i(p), \zeta_0(p) \right\rangle\\ + & \vdots\\ + \frac{dz_L}{dt} &= \left \langle f\left(\sum_{i=0}^L z_i(t) \zeta_i(p), p, t \right), \sum_{i=0}^L z_i(0) \zeta_i(p), \zeta_L(p) \right\rangle +\end{aligned} +``` +with initial condition +```math +\begin{aligned} + z_i(0) = \left\langle x_0, \zeta_i \right\rangle. +\end{aligned} +``` +These equations are referred to as the moment equations and can in principle be solved with any suitable IVP solution technique. + +## Example +In this section, we will showcase how `ModelOrderReduction.jl` may be used to apply PCE for the approximate solution of random ODEs as described in the previous section. To that end, we will consider a simple nonlinear reactor in which the following reaction occurs +```math +\begin{aligned} + &A + B \xleftrightharpoons[k_1(1-\frac{\theta_1}{2})]{k_2} C\\ + &C \xrightarrow{k_3(1-\frac{\theta_2}{2})} D\\ + &B \xrightarrow{k_4} D +\end{aligned} +``` +where the rate parameters $\theta_1 \sim U[-1,1]$, $\theta_2 \sim U[-1,1]$ are uncertain. Using [ModelingToolkit.jl]{github.com/SciML/ModelingToolkit.jl} the following code creates the associated model. +``` +using ModelOrderReduction, ModelingToolkit, DifferentialEquations, PolyChaos, Plots + +@parameters k[1:4], θ[1:2] +@variables t, c(t)[1:4] +D = Differential(t) +reactor = [D(c[1]) ~ -k[1]*(1+0.5*θ[1])*c[1]*c[2] + k[2]*c[3]; + D(c[2]) ~ -k[1]*(1+0.5*θ[1])*c[1]*c[2] - k[4]*c[2] + k[2]*c[3]; + D(c[3]) ~ k[1]*c[1]*c[2] - k[3]*c[3]; + D(c[4]) ~ k[3]*(1+0.5*θ[2])*c[3] + k[4]*c[2]]; + +@named reactor_model = ODESystem(reactor, t, c, vcat(k, θ)) +``` + +Next, we are going to define the PCE Ansatz with the help of which we are going to approximate the parametric solution. To that end, we heavily rely on [PolyChaos.jl]{github.com/SciML/PolyChaos.jl}. In this example, we choose a PCE of degree 3. +``` +d_pce = 3 +pce = PCE(c, [θ[i] => Uniform_11OrthoPoly(d_pce) for i in eachindex(θ)]) +``` + +Next, we are constructing the corresponding moment equations via a simple function call of `moment_equations`. This function returns an `ModelingToolkit.ODESystem` describing the deterministic evolution equation of the PCE moments and a function that maps the state of this model via the PCE Ansatz back to the original model state (in this example, the concentration of the different species). +``` +moment_eqs, pce_eval = moment_equations(reactor_model, pce) +``` + +Next, we are solving the evolution equation for the PCE moments. The initial condition for the moments is also obtained via Galerkin projection onto the space spanned by the basis functions in the PCE Ansatz. +``` +c0 = [1.0, 2.0, 0.0, 0.0] +z0 = reduce(vcat, pce_galerkin(c0, pce)) +moment_prob = ODEProblem(moment_eqs, z0, (0.0, 10.0), [k[1] => 1.0, k[2] => 0.2, k[3] => 12, k[4] => 0.1]) +moment_sol = solve(moment_prob, Tsit5()) +``` + +Lastly, we wish to demonstrate that this approximation can be remarkably accurate in practice. The figure below shows the PCE approximation (dashed black lines) compared to the full model for several realization of the uncertain paramters. + +![traces](traces.png) + +While PCE can even provide a good parametric surrogate model of the true solution of the random ODE, it is really designed to capture the statistics of the solution of the random ODE well. Moreover, any statistics that can be computed in terms of the moments of the states of the random ODE are remarkably cheap to compute from the PCE coefficients/moments alone. The figures below shows an example for the mean and variances of the different species. + + +![mean](mean.png) +![var](var.png) + + + +## References +Gautschi, Walter. Orthogonal polynomials: computation and approximation. OUP Oxford, 2004. \ No newline at end of file diff --git a/examples/PCE/var.png b/examples/PCE/var.png new file mode 100644 index 0000000..903a460 Binary files /dev/null and b/examples/PCE/var.png differ diff --git a/src/DataReduction/POD.jl b/src/DataReduction/POD.jl index 21b5be3..2c5d6e8 100644 --- a/src/DataReduction/POD.jl +++ b/src/DataReduction/POD.jl @@ -1,6 +1,3 @@ -using TSVD: tsvd -using RandomizedLinAlg: rsvd - function matricize(VoV::Vector{Vector{T}}) where {T} reduce(hcat, VoV) end diff --git a/src/ModelOrderReduction.jl b/src/ModelOrderReduction.jl index fa2422d..ccde7d1 100644 --- a/src/ModelOrderReduction.jl +++ b/src/ModelOrderReduction.jl @@ -1,23 +1,31 @@ module ModelOrderReduction +using Combinatorics: powerset, combinations +using Distributions using DocStringExtensions - -using Symbolics -using ModelingToolkit using LinearAlgebra - +using ModelingToolkit +using PolyChaos +using RandomizedLinAlg: rsvd using Setfield +using SparseArrays +using Symbolics +using TSVD: tsvd include("utils.jl") include("Types.jl") include("ErrorHandle.jl") - include("DataReduction/POD.jl") -export SVD, TSVD, RSVD -export POD, reduce! +export POD, reduce!, matricize include("deim.jl") export deim +include("pce/pce.jl") +export PCE +include("pce/pce_metadata.jl") +include("pce/stochastic_galerkin.jl") +export stochastic_galerkin + end diff --git a/src/pce/pce.jl b/src/pce/pce.jl new file mode 100644 index 0000000..68c8031 --- /dev/null +++ b/src/pce/pce.jl @@ -0,0 +1,582 @@ +""" +$(TYPEDSIGNATURES) +Returns the number of monomials in the multinomial expansion +``(x_1 + x_2 + ⋯ + x_m)^n``. +""" +multi_indices_size(m::Integer, n::Integer)::Int = binomial(m + n - 1, n) +""" +$(TYPEDSIGNATURES) +Compute the count of solutions ``(a_1, a_2, …, a_m)`` with non-negative integers +``a_i`` such that +```math +\\begin{gather*} +\\operatorname{mn} ≤ a_1 + a_2 + ⋯ + a_m ≤ \\operatorname{mx} \\\\ +0 ≤ a_i ≤ r_i \\quad ∀ i +\\end{gather*} +``` +The count for equality +```math +a_1 + a_2 + ⋯ + a_m = n +``` +by using the inclusion–exclusion principle is +```math +∑_{S⊆\\{1,2,…,m\\}}(-1)^{|S|}\\binom{n+m-1-∑_{i∈S}(r_i+1)}{m-1} +``` +where ``\\binom{p}{q}=0`` when ``p < 0``. + +The final result is the sum of this formula with ``n`` ranging from `mn` to `mx`. +""" +function multi_indices_size(r::AbstractVector{<:Integer}, mn::Integer, mx::Integer)::Int + m = length(r) + m₁ = m - 1 + sum(powerset(1:m); init = 0) do S + l = length(S) + temp = m₁ - l - sum(@view r[S]; init = 0) + res = sum(binomial(a, m₁) for a in (mx + temp):-1:max(1, mn + temp); init = 0) + iseven(l) ? res : -res + end +end +""" +$(TYPEDSIGNATURES) +Compute the count of solutions ``(a_1, a_2, …, a_m)`` with non-negative integers +``a_i`` such that +```math +\\begin{gather*} +a_1 + a_2 + ⋯ + a_m ≤ \\max_i a_i \\\\ +0 ≤ a_i ≤ r_i \\quad ∀ i +\\end{gather*} +``` +""" +function multi_indices_size(r::AbstractVector{<:Integer})::Int + m = length(r) + mn, mx = extrema(r) + res = sum(multi_indices_size(m, n) for n in 0:mn) + if mn == mx + return res + end + res + multi_indices_size(r, mn + 1, mx) +end + +""" +$(TYPEDSIGNATURES) +Return a matrix where the columns are the degrees ``(d_1, d_2, …, d_m)`` of monomials +``x_1^{d_1}x_2^{d_2}⋯x_m^{d_m}`` in the graded lexicographic order such that +```math +\\begin{gather*} +d_1 + d_2 + ⋯ + d_m ≤ \\max_i r_i \\\\ +0 ≤ d_i ≤ r_i \\quad ∀ i +\\end{gather*} +``` + +As `$(FUNCTIONNAME)` has exponential space complexity, [`multi_indices_size`](@ref) is used +to compute the matrix size in order to reduce the number of allocations. +""" +function grlex(r::AbstractVector{<:Integer})::Matrix{Int} + mn, mx = extrema(r) + indices_size = multi_indices_size(r) + n_term = length(r) + res = zeros(Int, n_term, indices_size) + indices_i = 2 + @inbounds for total_degree in 1:mn + for stars in combinations(1:(n_term + total_degree - 1), total_degree) + degree = @view res[:, indices_i] + for (i, s) in enumerate(stars) + degree[s - i + 1] += 1 + end + indices_i += 1 + end + end + @inbounds for total_degree in (mn + 1):mx + for stars in combinations(1:(n_term + total_degree - 1), total_degree) + degree = @view res[:, indices_i] + for (i, s) in enumerate(stars) + degree[s - i + 1] += 1 + end + if any(degree .> r) + degree .= 0 + continue + end + indices_i += 1 + if indices_i > indices_size + return res + end + end + end + res +end + +""" +$(TYPEDEF) +`$(FUNCTIONNAME)` represents multivariate orthogonal polynomial bases formed as the tensor +product of univariate `PolyChaos.AbstractOrthoPoly`s. + +`$(FUNCTIONNAME)` is similar to `PolyChaos.MultiOrthoPoly`, but allows different degree +truncation for univariate bases. Multi-indices of `$(FUNCTIONNAME)` are stored in the +columns as Julia is column-major, while multi-indices are rows in the index matrix of +`PolyChaos.MultiOrthoPoly`. + +By default, the sum of multi-indices are restricted to the maximum degree among the +univariate bases. + +# Fields +$(TYPEDFIELDS) +""" +struct TensorProductOrthoPoly{OP <: Union{AbstractOrthoPoly, AbstractCanonicalOrthoPoly}} <: + AbstractOrthoPoly{ProductMeasure, AbstractQuad{Float64}} + "The degree truncation of each univariate orthogonal polynomials." + deg::Vector{Int} + "Multi-indices in the columns." + ind::Matrix{Int} + "Product measure." + measure::ProductMeasure + "Univariate orthogonal polynomials." + uni::Vector{OP} +end +""" +$(TYPEDSIGNATURES) +""" +function TensorProductOrthoPoly(ops::AbstractVector{ + <:Union{AbstractOrthoPoly, + AbstractCanonicalOrthoPoly}}) + if any(op.quad isa EmptyQuad for op in ops) + throw(InconsistencyError("at least one quadrature rule missing")) + end + degrees = deg.(ops) + ind = grlex(degrees) + measures = [op.measure for op in ops] + w(t) = prod(m.w(t) for m in measures) + measure = ProductMeasure(w, measures) + TensorProductOrthoPoly(degrees, ind, measure, Vector(ops)) +end + +PolyChaos.dim(tpop::TensorProductOrthoPoly) = size(tpop.ind, 2) + +function PolyChaos.computeTensorizedSP(dim::Integer, tpop::TensorProductOrthoPoly) + computeTensorizedSP(dim, tpop.uni, transpose(tpop.ind)) +end + +""" +$(TYPEDSIGNATURES) +Construct a `$(FUNCTIONNAME)` which is used to compute and store the results of inner +products +```math +⟨ϕ_{i_1}ϕ_{i_2}⋯ϕ_{i_{m-1}},ϕ_{i_m}⟩ +``` +where ``m`` is input argument `dim`. +""" +function PolyChaos.Tensor(dim::Int, tpop::TensorProductOrthoPoly) + tensor_entries = computeTensorizedSP(dim, tpop) + getfun(ind) = getentry(ind, tensor_entries, transpose(tpop.ind), dim) + Tensor(dim, tensor_entries, getfun, tpop) +end + +""" +$(TYPEDEF) +Represents product of orthogonal polynomial basis functions +``Ψ_{\\vec α_1} Ψ_{\\vec α_2} ⋯ Ψ_{\\vec α_m}``. + +This is used to combine product of basis functions into one single object and conveniently +take inner product tensor ``⟨⋅, Ψ_{\\vec β}⟩``. + +# Fields +$(TYPEDFIELDS) +""" +struct BasisProduct + "The indices of the multi-indices ``\\vec α_1, \\vec α_2, …, \\vec α_m``." + idx::Vector{Int} +end +function Base.:+(p::BasisProduct, s) + if SymbolicUtils.isadd(s) + d = SymbolicUtils._merge(+, s.dict, Base.ImmutableDict(p => 1); + filter = SymbolicUtils._iszero) + return SymbolicUtils.Add(Float64, s.coeff, d) + end + SymbolicUtils.Add(Float64, SymbolicUtils.makeadd(1, 0, s, p)...) +end +Base.:+(s, p::BasisProduct) = p + s +Base.:*(p::BasisProduct) = p +Base.:*(p::BasisProduct, q::BasisProduct) = BasisProduct(vcat(p.idx, q.idx)) +function Base.:*(p::BasisProduct, s) + if s isa Number + if iszero(s) + return s + elseif isone(s) + return p + end + return SymbolicUtils.Mul(Float64, SymbolicUtils.makemul(s, p)...) + elseif SymbolicUtils.ismul(s) + d = copy(s.dict) + for (k, v) in s.dict + if k isa BasisProduct + delete!(d, k) + d[p * k^v] = 1 + return SymbolicUtils.Mul(SymbolicUtils.symtype(s), s.coeff, d) + end + end + d[p] = 1 + return SymbolicUtils.Mul(SymbolicUtils.symtype(s), s.coeff, d) + end + SymbolicUtils.Mul(Float64, SymbolicUtils.makemul(1, s, p)...) +end +Base.:*(s, p::BasisProduct) = p * s +function Base.:^(p::BasisProduct, exp) + if exp isa Number + if iszero(exp) + return 1 + elseif isone(exp) + return p + elseif exp isa Integer && exp > 1 + return BasisProduct(repeat(p.idx, exp)) + end + end + SymbolicUtils.Pow(p, exp) +end +Base.:(==)(p::BasisProduct, q::BasisProduct) = p.idx == q.idx +SymbolicUtils.:<ₑ(p::BasisProduct, q::BasisProduct) = p.idx < q.idx +SymbolicUtils.issym(::BasisProduct) = true +SymbolicUtils.symtype(::BasisProduct) = Float64 +Base.nameof(p::BasisProduct) = Symbol(p) +Base.hash(p::BasisProduct, h::UInt) = hash(p.idx, h) + +""" +$(TYPEDEF) +Suppose a variable ``Y`` with finite variance is a function of ``n`` independent but not +identically distributed random variables ``X_1, …, X_n`` with joint density +``p(x_1, …, x_n) = p_1(x_1) p_2(x_2) ⋯ p_n(x_n)``. Then the Polynomial Chaos Expansion +(PCE) for ``Y = Y(X_i), X_i ∼ π_{x_i}, i = 1, …, n``, takes the form +```math +y(x_1, …, x_n) = ∑_{α_1=0}^∞ ∑_{α_2=0}^∞ ⋯ ∑_{α_n=0}^∞ +C_{(α_1, α_2, …, α_n)} Ψ_{(α_1, α_2, …, α_n)}(x_1, …, x_n) +``` +Here the summation runs across all possible combinations of the multi-index +``\\vec α = (α_1, …, α_n)``. + +The set of multivariate orthogonal polynomials is defined as +```math +Ψ_{\\vec α}(x_1, …, x_n) = ∏_{i=1}^n ψ_{α_i}^{(i)}(x_i) +``` +where ``\\{ψ_{α_i}^{(i)}(x_i)\\}_{α_i=0}^∞`` is the family of univariate orthogonal +polynomials with respect to ``p_i``. + +# Fields +$(TYPEDFIELDS) +""" +struct PCE + "The random variables ``\\mathbf Y`` that are represented by other random variables." + states::Vector{Num} + "The independent random variables ``\\mathbf X``." + parameters::Vector{Num} + """ + The tensor-product-based multivariate basis underpinning the PCE, which stores the + univariate orthogonal polynomial basis ``\\{ψ_{α_i}^{(i)}(x_i)\\}_{α_i=0}^{r_i}`` for + each ``X_i`` and multi-indices ``\\vec α``. + """ + tensor_basis::TensorProductOrthoPoly + "The coefficients ``C_{(α_1, α_2, …, α_n)}`` for each ``Y_i`` in the columns." + moments::Matrix + "The mapping ``Y => ∑_α C_α Ψ_α`` for each ``Y``." + ansatz::Dict{Num, Num} + "Results of tensor inner products ``⟨Ψ_{i_1}Ψ_{i_2}⋯Ψ_{i_{m-1}},Ψ_{i_m}⟩``." + tensors::Dict{Int, Tensor} +end +""" +$(TYPEDSIGNATURES) +Construct a `$(FUNCTIONNAME)` object. + +In practice, the infinite series of multi-indices must be truncated. By default, besides +the upper bound for the degree of each univariate orthogonal polynomial, the total degree +is restricted to the maximum degree among the univariate bases. That is +```math +\\begin{gather*} +α_1 + α_2 + ⋯ α_n ≤ \\max_i r_i \\\\ +0 ≤ α_i ≤ r_i \\quad ∀ i +\\end{gather*} +``` +where ``r_i`` is the degree upper bound for the univariate basis corresponding to ``X_i``. + +# Arguments +- `states`: Random vairables ``\\mathbf Y`` that are represented by other random variables. +- `ivs`: Independent vairables of ``\\mathbf Y``. Enter an empty vector if there is none. +- `parameters`: Independent random variables ``\\mathbf X``. +- `uni_basis`: Univariate orthogonal polynomial basis ``\\{ψ_{α_i}^{(i)}(x_i)\\}_{α_i=0}^{r_i}`` for each ``X_i``. +""" +function PCE(states::AbstractVector{Num}, ivs::AbstractVector{Num}, + parameters::AbstractVector{Num}, + uni_basis::AbstractVector{ + <:Union{AbstractOrthoPoly, AbstractCanonicalOrthoPoly + }}) + states = Symbolics.scalarize(states) + parameters = Symbolics.scalarize(parameters) + tensor_basis = TensorProductOrthoPoly(uni_basis) + C = Matrix{Num}(undef, dim(tensor_basis), length(states)) + Ψ = Vector{Num}(undef, dim(tensor_basis)) + snames = Symbolics.tosymbol.(states) + for i in axes(tensor_basis.ind, 2) + name = join(Symbolics.map_subscripts.(view(tensor_basis.ind, :, i)), "ˏ") + Ψname = Symbol(:Ψ, name) + # record index of multi-index α⃗ in symbolic metadata + Ψ[i] = first(@variables $Ψname [description = i]) + for (j, sname) in enumerate(snames) + Cname = Symbol(:C, sname, name) + C[i, j] = if isempty(ivs) + first(@variables $Cname [description = i]) + else + first(@variables $Cname(ivs...) [description = i]) + end + end + end + # X => x₀ + x₁ψ₁(X) + # where x₀ and x₁ are affince PCE coefficients of X + # ψ₁(x) = x - α₀ is the first-order monic basis polynomial of X + # TODO: Use mean and std of the distribution of X to compute x₀ and x₁ + x_dict = Dict(x => op.α[1] + Ψ[i + 1] + for (i, (x, op)) in enumerate(zip(parameters, uni_basis))) + # Y => ∑_α C_α Ψ_α + y_dict = Dict(y => dot(view(C, :, i), Ψ) for (i, y) in enumerate(states)) + t2 = Tensor(2, tensor_basis) + tensors = Dict(2 => t2) + PCE(states, parameters, tensor_basis, C, y_dict, tensors) +end +function PCE(sys::ModelingToolkit.AbstractSystem) + ps = parameters(sys) # all parameters + stochastic_ps = filter(hasdist, ps) # stochastic parameters + dists = getdist.(stochastic_ps) + # TODO: find a way for inputing maximum degrees + uni_basis = to_ortho_poly.(dists, 3) # orthogonal polynomial basis + tensor_basis = TensorProductOrthoPoly(uni_basis) + Ψ = map(BasisProduct, 1:dim(tensor_basis)) # indices of the multi-indices + x_dict = Dict{Num, BasisProduct}() + sizehint!(x_dict, length(stochastic_ps)) # reserve sufficient capacity to avoid growing + for (i, (p, d, op)) in enumerate(zip(stochastic_ps, dists, uni_basis)) + # TODO: deal with non-canonical ortho poly + x₀, x₁ = convert2affinePCE(mean(d), std(d), op) + x_dict[p] = x₀ + x₁ * Ψ[i + 1] + end + # TODO: transformation of states + # TDDO: tensor +end + +""" +$(TYPEDSIGNATURES) + +Construct univariate monic orthogonal polynomials accroding to the input continuous +univariate distribution. + +The input argument `deg` is the maximum degree. +""" +function to_ortho_poly(::Normal, deg::Integer) + GaussOrthoPoly(deg) +end +# TODO: add other distributions +function to_ortho_poly(dist::ContinuousUnivariateDistribution, deg::Integer) + error("unimplmented") # TODO + w(t) = pdf(dist, t) # weight function + supp = extrema(dist) # minimum and maximum + meas = Measure("measure", w, supp, false, Dict()) + OrthoPoly("op", deg, meas) +end + +# Compute the an ascending list of `n`-dimensional multi-indices with fixed `grade` (= sum of entries) +# in graded reverse lexicographic order. Constraints on the degrees considered can be incorporated. +function grevlex(n::Int, grade::Int) + if n == 1 + return reshape([grade], 1, 1) + end + if grade == 0 + return zeros(Int, 1, n) + end + sub_ind = grevlex(n - 1, grade) + ind = hcat(sub_ind, zeros(Int, size(sub_ind, 1))) + for k in 1:grade + sub_ind = grevlex(n - 1, grade - k) + ind = vcat(ind, hcat(sub_ind, k * ones(Int, size(sub_ind, 1)))) + end + return ind +end +function grevlex(n::Int, grades::AbstractVector{Int}) + return reduce(vcat, [grevlex(n, grade) for grade in grades]) +end +function grevlex(n::Int, grade::Int, max_degrees::Vector{Int}) + return grevlex(n, grade, [0:d for d in max_degrees]) +end +function grevlex(n::Int, grades::AbstractVector{Int}, max_degrees::Vector{Int}) + return reduce(vcat, [grevlex(n, grade, max_degrees) for grade in grades]) +end +function grevlex(n::Int, grade::Int, degree_constraints::Vector{<:AbstractVector}) + if n == 1 + return grade in degree_constraints[1] ? reshape([grade], 1, 1) : zeros(Int, 0, 1) + end + if grade == 0 + return all(0 in degs for degs in degree_constraints) ? zeros(Int, 1, n) : + zeros(Int, 0, n) + end + filtered_grades = filter(x -> x <= grade, degree_constraints[end]) + sub_ind = grevlex(n - 1, grade - filtered_grades[1], degree_constraints[1:(end - 1)]) + ind = hcat(sub_ind, filtered_grades[1] * ones(Int, size(sub_ind, 1))) + for k in filtered_grades[2:end] + sub_ind = grevlex(n - 1, grade - k, degree_constraints[1:(end - 1)]) + ind = vcat(ind, hcat(sub_ind, k * ones(Int, size(sub_ind, 1)))) + end + return ind +end +function grevlex(n::Int, grades::AbstractVector{Int}, + degree_constraints::Vector{<:AbstractVector}) + return reduce(vcat, [grevlex(n, grade, degree_constraints) for grade in grades]) +end + +# 1. apply PCE ansatz +""" +$(TYPEDSIGNATURES) + +Generate linear PCEs for the uncertain parameters. +""" +function generate_parameter_pce(pce::PCE) + par_dim = length(pce.parameters) + par_pce = Vector{Pair{eltype(pce.parameters), eltype(pce.sym_basis)}}(undef, par_dim) + for (i, basis) in enumerate(pce.uni_basis) + p, op = basis + par_pce[i] = p => pce.sym_basis[i + 1] + op.α[1] + end + return par_pce +end + +""" +$(TYPEDSIGNATURES) + +Substitute parameters shared between a set of symbolic `eqs` and the PCE `pce` for the corresponding linear PCEs. +""" +function substitute_parameters(eqs::AbstractVector, pce::PCE) + par_pce = generate_parameter_pce(pce) + subs_eqs = [substitute(eq, par_pce) for eq in eqs] + return subs_eqs +end + +""" +$(TYPEDSIGNATURES) + +Substitute PCE Ansatz defined in `pce` into a set of symbolic equations `eqs`. +""" +function substitute_pce_ansatz(eqs::AbstractVector, pce::PCE) + subs_eqs = [expand(expand(substitute(eq, pce.ansatz))) for eq in eqs] + return subs_eqs +end + +""" +$(TYPEDSIGNATURES) + +Apply PCE ansatz defined in `pce` to a given set of symbolic equations `eqs`. +""" +function apply_ansatz(eqs::AbstractVector, pce::PCE) + return substitute_pce_ansatz(substitute_parameters(eqs, pce), pce) +end + +# 2. extract PCE expansion coeffs +""" +$(TYPEDSIGNATURES) + +Given a set of symbolic equations `eqs` involving the basis functions of `pce`, +extract monomials of the basis functions and the corresponding coeffiecients. + +# Returns +`Vector` of `Dict`s mapping monomial of basis functions to its coefficient in the individual equations. +""" +function extract_basismonomial_coeffs(eqs::AbstractVector, pce::PCE) + basismonomial_coeffs = [extract_coeffs(eq, pce.sym_basis) for eq in eqs] + basismonomial_indices = [] + for coeffs in basismonomial_coeffs + union!(basismonomial_indices, + [mono => get_basis_indices(mono) for mono in keys(coeffs)]) + end + return basismonomial_coeffs, basismonomial_indices +end + +# 3. compute inner products +""" +$(TYPEDSIGNATURES) + +Evaluate scalar products between all basis functions in `pce` and +basis monomials as characterized by `mono_indices`. +""" +function eval_scalar_products(mono_indices, pce::PCE) + uni_degs = [deg(op) for op in pce.tensor_basis.uni] + max_degs = uni_degs + for (mono, id) in mono_indices + max_degs = max.(max_degs, vec(sum(pce.tensor_basis.ind[id .+ 1, :], dims = 1))) + end + quad_deg = max.(ceil.(Int, 0.5 * (max_degs + uni_degs .+ 1))) + + integrators = map((uni, deg) -> bump_degree(uni, deg), pce.tensor_basis.uni, quad_deg) + scalar_products = Dict() + for k in 1:dim(pce.tensor_basis) + scalar_products[k] = Dict(mono => computeSP(vcat(id, k - 1), pce.tensor_basis, + integrators) + for (mono, id) in mono_indices) + end + return scalar_products +end + +# 4. Galerkin projection +""" +$(TYPEDSIGNATURES) + +perform Galerkin projection of polynomial expressions characterized by `Dict`s mapping +basis monomials to coefficients. +""" +function galerkin_projection(bm_coeffs::Vector{<:Dict}, scalar_products::Dict, + pce::PCE) + projected_eqs = [] + scaling_factors = computeSP2(pce.tensor_basis) + for i in eachindex(bm_coeffs) + eqs = [] + for k in 1:dim(pce.tensor_basis) + push!(eqs, + sum(bm_coeffs[i][mono] * scalar_products[k][mono] + for mono in keys(bm_coeffs[i]))) + end + push!(projected_eqs, eqs ./ scaling_factors) + end + return projected_eqs +end + +# 5. combine everything +""" +$(TYPEDSIGNATURES) + +perform Galerkin projection onto the `pce`. +""" +function pce_galerkin(eqs::AbstractVector, pce::PCE) + expanded_eqs = apply_ansatz(eqs, pce) + basismono_coeffs, basismono_idcs = extract_basismonomial_coeffs(expanded_eqs, pce) + scalar_products = eval_scalar_products(basismono_idcs, pce) + projected_eqs = galerkin_projection(basismono_coeffs, scalar_products, pce) + return projected_eqs +end + +# 6. high-level interface +# 6a. apply pce to explicit ODE +""" +$(TYPEDSIGNATURES) + +Generate moment equations of an `ODESystem` from a given `PCE`-Ansatz via Galerkin projection. +""" +function moment_equations(sys::ODESystem, pce::PCE) + eqs = [eq.rhs for eq in equations(sys)] + projected_eqs = pce_galerkin(eqs, pce) + moment_eqs = reduce(vcat, projected_eqs) + iv = ModelingToolkit.get_iv(sys) + params = setdiff(parameters(sys), pce.parameters) + D = Differential(iv) + moments = reduce(vcat, pce.moments) + name = Symbol(String(nameof(sys)) * "_pce") + pce_system = ODESystem([D(moments[i]) ~ moment_eqs[i] for i in eachindex(moments)], + iv, moments, params, name = name) + + n_moments = dim(pce.tensor_basis) + n_states = length(states(sys)) + pce_eval = function (moment_vals, parameter_values) + shape_state = [moment_vals[(i * n_moments + 1):((i + 1) * n_moments)] + for i in 0:(n_states - 1)] + return pce(shape_state, parameter_values) + end + return pce_system, pce_eval +end diff --git a/src/pce/pce_metadata.jl b/src/pce/pce_metadata.jl new file mode 100644 index 0000000..3e24c41 --- /dev/null +++ b/src/pce/pce_metadata.jl @@ -0,0 +1,15 @@ +""" +$(TYPEDEF) + +An object used as the field `metadata` in `ModelingToolkit.AbstractSystem` to store the +information about the Polynomial Chaos Expansion of a stochastic system. + +# Fields +$(TYPEDFIELDS) +""" +struct PCEMetadata{S <: ModelingToolkit.AbstractSystem} + "A PCE object." + pce::PCE + "The original stochastic system on which PCE is applied." + sys::S +end diff --git a/src/pce/stochastic_galerkin.jl b/src/pce/stochastic_galerkin.jl new file mode 100644 index 0000000..38ae382 --- /dev/null +++ b/src/pce/stochastic_galerkin.jl @@ -0,0 +1,8 @@ +""" +$(TYPEDSIGNATURES) +""" +function stochastic_galerkin(sys::ODESystem, pce::PCE; + name::Symbol = Symbol(nameof(sys), :_pce), kwargs...) + iv = ModelingToolkit.get_iv(sys) + eqs = full_equations(sys; kwargs...) +end diff --git a/src/utils.jl b/src/utils.jl index 69b2b54..e19384d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,5 +1,3 @@ -using SparseArrays - """ $(TYPEDSIGNATURES) @@ -100,3 +98,61 @@ function linear_terms(exprs::AbstractVector, vars) linear = sparse(linear_I, linear_J, linear_V, length(exprs), length(vars)) return linear, const_terms, nonlinear_terms end + +# utiltiites to extracting coefficients of a polynomial in monomial basis in variables `vars` +function split_term(term::Symbolics.Mul, vars) + coeff = term.coeff + mono = Num(1.0) + for var in keys(term.dict) + if var in vars + mono *= var^term.dict[var] + else + coeff *= var^term.dict[var] + end + end + return coeff, mono +end +function split_term(term::Symbolics.Pow, vars) + if term.base in vars + return 1.0, term + else + return term, Val(1) + end +end +function split_term(term::T, vars) where {T <: Union{Symbolics.Term, Symbolics.Sym}} + if term in vars + return 1.0, term + else + return term, Val(1) + end +end + +function extract_coeffs(expr::Symbolics.Add, vars::Set) + coeffs = Dict() + if !iszero(expr.coeff) + coeffs[Val(1)] = expr.coeff + end + for term in keys(expr.dict) + num_coeff = expr.dict[term] + var_coeff, mono = split_term(term, vars) + try + coeffs[mono] += num_coeff * var_coeff + catch + coeffs[mono] = num_coeff * var_coeff + end + end + return coeffs +end +function extract_coeffs(expr::T, vars::Set) where {T <: Union{Symbolics.Sym, Symbolics.Mul}} + coeff, mono = split_term(expr, vars) + return Dict(mono => coeff) +end +extract_coeffs(expr::Num, vars::Set) = extract_coeffs(Symbolics.unwrap.(expr), vars) +function extract_coeffs(expr::Num, vars::AbstractArray) + extract_coeffs(Symbolics.unwrap.(expr), vars) +end +function extract_coeffs(expr, vars::AbstractArray{<:Num}) + extract_coeffs(expr, Symbolics.unwrap.(vars)) +end +extract_coeffs(expr, vars::AbstractArray) = extract_coeffs(expr, Set(vars)) +extract_coeffs(expr::Number, vars::Set) = Dict(Val(1) => expr) diff --git a/test/DataReduction.jl b/test/DataReduction.jl index da9ea80..811f550 100644 --- a/test/DataReduction.jl +++ b/test/DataReduction.jl @@ -1,5 +1,6 @@ using Test, ModelOrderReduction using DifferentialEquations +import ModelOrderReduction as MOR function lorenz_prob() function lorenz!(du, u, p, t) @@ -29,7 +30,7 @@ end solution = Array(sol) order = 2 - solver = SVD() + solver = MOR.SVD() matrix_reducer = POD(solution, order) snapshot_reducer = POD(sol.u, order) reduce!(matrix_reducer, solver) @@ -48,7 +49,7 @@ end @test reducer.nmodes == 1 order = 2 - solver = TSVD() + solver = MOR.TSVD() reducer = POD(solution, order) reduce!(reducer, solver) @@ -57,7 +58,7 @@ end @test reducer.renergy > 0.7 order = 2 - solver = RSVD() + solver = MOR.RSVD() reducer = POD(solution, order) reduce!(reducer, solver) diff --git a/test/Project.toml b/test/Project.toml index e8910b8..7e72c62 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,6 +2,8 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +PolyChaos = "8d666b04-775d-5f6e-b778-5ac7c70f65a3" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/pce.jl b/test/pce.jl new file mode 100644 index 0000000..c8ec689 --- /dev/null +++ b/test/pce.jl @@ -0,0 +1,205 @@ +using Test, ModelOrderReduction +using Symbolics, PolyChaos +using SpecialFunctions: beta +import ModelOrderReduction as MOR + +supp = (-1, 1) +w(t) = 1 + t +my_meas = Measure("my_meas", w, supp, false, Dict()) +deg = 4 +Nrec = 2 * deg +my_op = OrthoPoly("my_op", deg, my_meas; Nrec) + +@testset "PCE constructor" begin + @variables t x[1:4] y(t)[1:2] + degrees = [3, 4, 5, 4] + Nrec = 2 * minimum(degrees) + 1 + uni_basis = [ + HermiteOrthoPoly(degrees[1]; Nrec), + Uniform01OrthoPoly(degrees[2]; Nrec), + Uniform_11OrthoPoly(degrees[3]; Nrec), + OrthoPoly("my_op", degrees[4], my_meas; Nrec), + ] + pce = @test_nowarn PCE(y, [t], x, uni_basis) + @test isequal(pce.states, Symbolics.scalarize(y)) + @test isequal(pce.parameters, Symbolics.scalarize(x)) + multi_indices_dim = MOR.multi_indices_size(degrees) + @test size(pce.moments) == (multi_indices_dim, length(y)) + Nrec = 2 * minimum(degrees[1:3]) + 1 + uni_basis = [ + HermiteOrthoPoly(degrees[1]; Nrec), + Uniform01OrthoPoly(degrees[2]; Nrec), + Uniform_11OrthoPoly(degrees[3]; Nrec), + ] + @test_nowarn PCE(y, [t], x[1:3], uni_basis) + Nrec = 2 * minimum(degrees[1:2]) + 1 + uni_basis = [HermiteOrthoPoly(degrees[1]; Nrec), HermiteOrthoPoly(degrees[2]; Nrec)] + @test_nowarn PCE(y, [t], x[1:2], uni_basis) + uni_basis = [my_op, my_op] + @test_nowarn PCE(y, [t], x[1:2], uni_basis) +end + +@testset "TensorProductOrthoPoly" begin + ops = [GaussOrthoPoly(1)] + @test_nowarn MOR.TensorProductOrthoPoly(ops) + ops = [GaussOrthoPoly(4), GaussOrthoPoly(3)] + @test_nowarn MOR.TensorProductOrthoPoly(ops) + ops = [my_op] + @test_nowarn MOR.TensorProductOrthoPoly(ops) + ops = [my_op, my_op] + @test_nowarn MOR.TensorProductOrthoPoly(ops) + ops = [GaussOrthoPoly(4), Uniform01OrthoPoly(2), LaguerreOrthoPoly(3), my_op] + tpop = @test_nowarn MOR.TensorProductOrthoPoly(ops) + @test dim(tpop) == MOR.multi_indices_size([4, 2, 3, 4]) + ops = [GaussOrthoPoly(4; addQuadrature = false), Uniform01OrthoPoly(2)] + @test_throws InconsistencyError MOR.TensorProductOrthoPoly(ops) +end + +@testset "multi_indices_size" begin + @test MOR.multi_indices_size([0]) == 1 + @test MOR.multi_indices_size([1]) == 2 + @test MOR.multi_indices_size([5]) == 6 + @test MOR.multi_indices_size([2, 2]) == 6 + @test MOR.multi_indices_size([4, 2]) == 12 + @test MOR.multi_indices_size([2, 2, 2]) == 10 +end + +@testset "grlex" begin + @test MOR.grlex([0]) == zeros(Int, 1, 1) + @test MOR.grlex([1]) == [0 1] + @test MOR.grlex([5]) == [0 1 2 3 4 5] + @test MOR.grlex([2, 2]) == [0 1 0 2 1 0 + 0 0 1 0 1 2] + @test MOR.grlex([4, 2]) == [0 1 0 2 1 0 3 2 1 4 3 2 + 0 0 1 0 1 2 0 1 2 0 1 2] + @test MOR.grlex([2, 2, 2]) == + [0 1 0 0 2 1 1 0 0 0 + 0 0 1 0 0 1 0 2 1 0 + 0 0 0 1 0 0 1 0 1 2] + @test MOR.grlex([3, 4, 2]) == + [0 1 0 0 2 1 1 0 0 0 3 2 2 1 1 1 0 0 0 3 3 2 2 2 1 1 1 0 0 0 + 0 0 1 0 0 1 0 2 1 0 0 1 0 2 1 0 3 2 1 1 0 2 1 0 3 2 1 4 3 2 + 0 0 0 1 0 0 1 0 1 2 0 0 1 0 1 2 0 1 2 0 1 0 1 2 0 1 2 0 1 2] + + function check_grlex(r::AbstractVector{<:Integer}) + res = MOR.grlex(r) + mx = maximum(r) + @test allunique(view(res, :, i) for i in axes(res, 2)) + for i in axes(res, 2) + degree = @view res[:, i] + @test sum(degree) ≤ mx + @test all(degree .≤ r) + end + end + check_grlex([4, 1, 5]) + check_grlex([3, 2, 2]) + check_grlex([3, 2, 1, 2]) + check_grlex([2, 5, 3, 4, 2]) +end + +@testset "Tensor" begin + deg, n = 4, 20 + s_α, s_β = 2.1, 3.2 + beta_op = Beta01OrthoPoly(deg, s_α, s_β; Nrec = n, addQuadrature = true) + supp = (0, 1) + w(t) = (t^(s_α - 1) * (1 - t)^(s_β - 1) / beta(s_α, s_β)) + my_meas = Measure("my_meas", w, supp, false) + my_opq = OrthoPoly("my_op", deg, my_meas; Nrec = n, addQuadrature = true) + + ops = [beta_op, my_opq] + tpop = MOR.TensorProductOrthoPoly(ops) + mop = MultiOrthoPoly(ops, deg) + + for t in (2, 3) + tt = Tensor(t, tpop) + mt = Tensor(t, mop) + index = Vector{Int}(undef, t) + for i in 0:(dim(mop) - 1) + fill!(index, i) + @test tt.get(index) == mt.get(index) + end + end +end + +# testing GRevLex capabilities +@testset "PCE: GRevLeX ordering" begin + @test MOR.grevlex(1, 5) == reshape([5], 1, 1) + + ind = [2 0 0; + 1 1 0; + 0 2 0; + 1 0 1; + 0 1 1; + 0 0 2] + @test ind == MOR.grevlex(3, 2) + + ind = [1 1 0; + 1 0 1; + 0 1 1; + 0 0 2] + @test ind == MOR.grevlex(3, 2, [1, 1, 2]) + + ind = [1 1 0; + 0 1 1] + @test ind == MOR.grevlex(3, 2, [0:1, [1], 0:2]) + + @test zeros(Int, 0, 3) == MOR.grevlex(3, 0, [0:1, [1], 0:2]) + + ind0 = [0 0 0] + ind1 = [1 0 0; + 0 1 0; + 0 0 1] + ind2 = [2 0 0; + 1 1 0; + 0 2 0; + 1 0 1; + 0 1 1; + 0 0 2] + @test vcat(ind0, ind1, ind2) == MOR.grevlex(3, 0:2) + @test vcat(ind1, ind2) == MOR.grevlex(3, 1:2) + @test vcat(ind0, ind2) == MOR.grevlex(3, [0, 2]) + @test vcat(ind2, ind0) == MOR.grevlex(3, [2, 0]) + + degree_constraints = [0:1, [1], 1:2] + ind0 = zeros(Int, 0, 3) + ind1 = zeros(Int, 0, 3) + ind2 = [0 1 1] + @test vcat(ind0, ind1, ind2) == MOR.grevlex(3, 0:2, degree_constraints) + @test vcat(ind1, ind2) == MOR.grevlex(3, 1:2, degree_constraints) + @test vcat(ind0, ind2) == MOR.grevlex(3, [0, 2], degree_constraints) + @test vcat(ind2, ind0) == MOR.grevlex(3, [2, 0], degree_constraints) +end + +@testset "BasisProduct" begin + @test MOR.BasisProduct(Int[]) == MOR.BasisProduct(Int[]) + @test MOR.BasisProduct([1, 3]) == MOR.BasisProduct([1, 3]) + + @variables x y z + yv = Symbolics.value(y) + Ψ₁ = MOR.BasisProduct([4]) + Ψ₂ = MOR.BasisProduct([7]) + + expr1 = x * y * z + dict = Dict(x => Ψ₁, z => Ψ₂) + expr2 = substitute(expr1, dict) + expr3 = expand(expr2) + @test isequal(expr3, yv * MOR.BasisProduct([4, 7])) + + expr4 = (yv + Ψ₁)^2 + expr5 = expand(expr4) + @test isequal(expr5, yv^2 + 2yv * MOR.BasisProduct([4]) + MOR.BasisProduct([4, 4])) + + expr6 = Ψ₁ * yv + @test expr6 isa SymbolicUtils.Mul + @test expr6.coeff == 1 + @test length(expr6.dict) == 2 + @test (Ψ₁ => 1) in expr6.dict + @test (yv => 1) in expr6.dict + + expr7 = 2 * Ψ₁ * yv * Ψ₂ + @test expr7 isa SymbolicUtils.Mul + @test expr7.coeff == 2 + @test length(expr7.dict) == 2 + @test (MOR.BasisProduct([7, 4]) => 1) in expr7.dict + @test (yv => 1) in expr7.dict +end diff --git a/test/runtests.jl b/test/runtests.jl index 8d64d7f..c2ebc94 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,6 @@ using SafeTestsets -@safetestset "POD" begin include("DataReduction.jl") end -@safetestset "utils" begin include("utils.jl") end -@safetestset "DEIM" begin include("deim.jl") end +@time @safetestset "POD" begin include("DataReduction.jl") end +@time @safetestset "utils" begin include("utils.jl") end +@time @safetestset "DEIM" begin include("deim.jl") end +@time @safetestset "PCE" begin include("pce.jl") end diff --git a/test/utils.jl b/test/utils.jl index 2515cf8..9065e9a 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,5 +1,6 @@ using Test, ModelOrderReduction using Symbolics +import ModelOrderReduction as MOR @variables t w(t) x(t) y(t) z(t) @@ -10,7 +11,7 @@ using Symbolics 2.0z + 3.4w + 7.0 + sin(x) 9.8 + x * (1.0 - y) 5.6y + 1.3z^2] - A, c, n = ModelOrderReduction.linear_terms(exprs, vars) + A, c, n = MOR.linear_terms(exprs, vars) @test size(A) == (length(exprs), length(vars)) @test A == [3.0 4.5 0.0 0.0 0.0 2.0