sciml / polychaos.jl Goto Github PK

A Julia package to construct orthogonal polynomials, their quadrature rules, and use it with polynomial chaos expansions.

Home Page: https://docs.sciml.ai/PolyChaos/stable/

License: MIT License

Julia 92.99% MATLAB 7.01%

polynomial-chaos-expansions orthogonal-polynomials quadrature-rules quadrature uncertainty-quantification uncertainty-propagation julia-language quadrature-integration polynomials differential-equations

polychaos.jl's Introduction

PolyChaos -- Orthogonal Polynomials, Quadrature, and Polynomial Chaos

A Julia package to construct orthogonal polynomials, their quadrature rules, and use it with polynomial chaos expansions.

Tutorials and Documentation

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation, which contains the unreleased features.

The package requires Julia 1.3 or newer. In Julia switch to the package manager

using Pkg
Pkg.add("PolyChaos")

This will install PolyChaos and its dependencies. Once that is done, load the package:

using PolyChaos

That's it.

Let's take a look at a simple example. We would like to solve the integral

Exploiting the underlying uniform measure, the integration can be done exactly with a 3-point quadrature rule.

opq = Uniform01OrthoPoly(3)
integrate(x -> 6x^5, opq)

For more information please visit the documentation.

Citing

If you like PolyChaos.jl, consider citing our paper

@ARTICLE{2020arXiv200403970M,
       author = {{M{\"u}hlpfordt}, Tillmann and {Zahn}, Frederik and {Hagenmeyer}, Veit and {Faulwasser}, Timm},
        title = "{PolyChaos.jl -- A Julia Package for Polynomial Chaos in Systems and Control}",
      journal = {arXiv e-prints},
     keywords = {Electrical Engineering and Systems Science - Systems and Control, Mathematics - Numerical Analysis, Mathematics - Optimization and Control},
         year = 2020,
        month = apr,
          eid = {arXiv:2004.03970},
        pages = {arXiv:2004.03970},
archivePrefix = {arXiv},
       eprint = {2004.03970},
 primaryClass = {eess.SY},
}

polychaos.jl's People

Contributors

Stargazers

Watchers

polychaos.jl's Issues

On unifying Array into AbstractArray

Under the current function declarations of Array, it may bring some inconvenience when we try to include more usage scenarios with arrays affiliated to AbstractArray. See an example.

First let's define a variable in polynomial chaos expansion with type of OffsetArray.

using PolyChaos
L = 5; op = GaussOrthoPoly(L, addQuadrature=true)

using OffsetArrays
nx = 10; vChaos = OffsetArray{Float64}(undef, 0:nx+1, 1:L+1); vChaos .= 1.

We would like to evaluate the PCE onto quadrature points.

vQuad = OffsetArray{Float64}(undef, 0:nx+1, 1:op.quad.Nquad)
for i in axes(pce)[1]
    vQuad[i,:] = evaluatePCE(vChaos[i,:], op.quad.nodes, op)
end

Unfortunately, here it should turn out:

MethodError: no method matching evaluatePCE(::OffsetArray{Float64,1,Array{Float64,1}}, ::Array{Float64,1}, ::GaussOrthoPoly)
Closest candidates are:
  evaluatePCE(!Matched::Array{#s99,1} where #s99<:Real, ::Array{#s98,1} where #s98<:Real, ::AbstractOrthoPoly) at /home/tianbai/.julia/packages/PolyChaos/1jqDA/src/polynomial_chaos.jl:229
  evaluatePCE(!Matched::Array{#s99,1} where #s99<:Real, ::Array{#s98,1} where #s98<:Real, !Matched::Array{#s97,1} where #s97<:Real, !Matched::Array{#s96,1} where #s96<:Real) at /home/tianbai/.julia/packages/PolyChaos/1jqDA/src/polynomial_chaos.jl:216

Stacktrace:
 [1] top-level scope at ./In[17]:9

It is due to the incompatibility in the function declarations. A temporary solution is to reformulate an array from AbstractArray, then put it into functions.

vQuad1 = OffsetArray{Float64}(undef, 0:nx+1, 1:op.quad.Nquad)
for i in axes(pce)[1]
    vQuad1[i,1:op.quad.Nquad] = evaluatePCE(vChaos[i,1:L+1], op.quad.nodes, op)
end

vQuad1
12×5 OffsetArray(::Array{Float64,2}, 0:11, 1:5) with eltype Float64 with indices 0:11×1:5:
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168
 11.2059  -2.5914  3.0  -3.03138  46.4168

A more general solution is to modify the function declarations correspondingly.

Include jldoctests in documentation

jldoctests are tested when building the documentation with Documenter.jl. This should be favored to @example or its friends.

Underflow in stieltjes() when trying to compute recursion coefficients for very high degree

Hello,

I am trying to compute a high number of recursion coefficients (on the order of 500).

Using the OrthoPoly constructor as in the tutorial returns the error

ERROR: DomainError with 2.2250738585072014e-307:
Underflow in stieltjes() for k=255; try using `removeZeroWeights`

I do not see a way to use `removeZeroWeights'.
I would be grateful for any help or tips in general about obtaining such high number of coefficients.

Compilation error with Polychaos

Hello,
following the tutorial, when typing "using Polychaos" on version 0.2.2, i get the following error:

ERROR: LoadError: EOFError: read end of file
Stacktrace:
 [1] read(::IOStream, ::Type{Int64}) at ./iostream.jl:361
 [2] parse_cache_header(::IOStream) at ./loading.jl:1323
 [3] stale_cachefile(::String, ::String) at ./loading.jl:1402
 [4] _require_search_from_serialized(::Base.PkgId, ::String) at ./loading.jl:757
 [5] _require(::Base.PkgId) at ./loading.jl:1006
 [6] require(::Base.PkgId) at ./loading.jl:927
 [7] require(::Module, ::Symbol) at ./loading.jl:922
 [8] include(::Module, ::String) at ./Base.jl:377
 [9] top-level scope at none:2
 [10] eval at ./boot.jl:331 [inlined]
 [11] eval(::Expr) at ./client.jl:449
 [12] top-level scope at ./none:3
in expression starting at /Users/ludovicobessi/.julia/packages/PolyChaos/Dhyx4/src/PolyChaos.jl:8
ERROR: Failed to precompile PolyChaos [8d666b04-775d-5f6e-b778-5ac7c70f65a3] to /Users/ludovicobessi/.julia/compiled/v1.4/PolyChaos/bwvo3_lR8zj.ji.

There was probably a corrupted file in my .julia/compiled/ directory, I solved this by following: https://stackoverflow.com/questions/59570508/error-eoferror-read-end-of-file-when-using-package-after-installing-new-juli

computeTensorizedSP

The computeTensorizedSP seems to drop entries less than 1E-10,
This is not very reasonable when using Float64 and can cause weird behaviors when making use of tensor products.

MWE

using PolyChaos

op = Uniform_11OrthoPoly(20, Nrec = 40)
t2 = Tensor(2, op)

l1 = [t2.get([i, i]) for i=0:op.deg]
l2 = computeSP2(op)

l1 is not equal to l2

julia> l1
21-element Vector{Float64}:
 ⋮
 3.6023076661943104e-10
 0.0
 0.0
 0.0
 0.0

julia> l2
21-element Vector{Float64}:
 ⋮
 3.6023076661943104e-10
 9.013566368226456e-11
 2.2551316627840712e-11
 5.64173617647297e-12
 1.4113161166911746e-12

computeTensorizedSP is the default method used to construct Tensor, see

PolyChaos.jl/src/typesTensor.jl

Line 15 in e20be07

tensorEntries = computeTensorizedSP(dim, opq)

In the following example, Tensor.get() would not work at high orders as the denominators are zero.

function ODEgalerkin(du,u,p,t)
   du[:] = [ sum( p[j+1]*u[k+1]*t3.get([j,k,m])/t2.get([m,m]) for j=0:L for k=0:L) for m=0:L ] 
end

Distinguish between dependencies and extra packages in `project.toml`

This issue can be tackled as soon as there is a new version of Pkg.jl: JuliaLang/Pkg.jl#876

As we need certain packages only for testing, we should specify this in project.toml, see the documentation of Pkg.

The project.toml was created like this, see commit 97f3af6, but that broke the build, see this post.

Multivariate integration

Hey, I tried following the tutorial examples here but I can't seem to figure out the right syntax for multivariate integration. Here is what I am trying:

deg = [3, 5, 6, 4]
d = minimum(deg)
op1 = GaussOrthoPoly(deg[1], addQuadrature=true);
op2 = Uniform01OrthoPoly(deg[2], addQuadrature=true);
op3 = Beta01OrthoPoly(deg[3], 2, 1.2, addQuadrature=true);
ops = [op1, op2, op3];
mop = MultiOrthoPoly(ops, d);
integrate((a,b,c,d) -> a*b*c*d, mop)

however, I get the following error:

type MultiOrthoPoly has no field quad

Stacktrace:
 [1] getproperty(::MultiOrthoPoly{ProductMeasure,Quad{Float64,Array{Float64,1}},Array{AbstractCanonicalOrthoPoly{Array{Float64,1},M,Quad{Float64,Array{Float64,1}}} where M,1}}, ::Symbol) at .\Base.jl:33
 [2] integrate(::Function, ::MultiOrthoPoly{ProductMeasure,Quad{Float64,Array{Float64,1}},Array{AbstractCanonicalOrthoPoly{Array{Float64,1},M,Quad{Float64,Array{Float64,1}}} where M,1}}) at C:\Users\mrra\.julia\packages\PolyChaos\AlY80\src\auxfuns.jl:96
 [3] top-level scope at In[5]:8
 [4] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

Is there any way of doing this with the PolyChaos package?

Tag?

This is causing some downstream compat issues that seem fixed on master.

Show formulae of orthogonal polynomials

We represent orthogonal polynomials entirely in terms of their recurrence coefficients (a_i, b_i). The orthogonal polynomial of degree n can be evaluated at points x using evaluate(n,x,a,b). Sometimes, however, it may be interesting to look at the formula of the orthogonal polynomial of degree n, for example phi_n(x) = x^n + c_n-1 x^(n-1) + ... + c_1 x + c_0. So the task is to compute c_i from a_i, b_i and come up with a pretty output.

Inverse transform sampling

It would be helpful to transpile the inverse transform sampling method from here to our package. Yes, I am aware that there is a working implementation withing ApproxFun.jl, but we would like to avoid having to add this heavy dependency (load up time increases drastically.).

p.s.: the implementation withing ApproxFun.jl can be found here.

Bump compat for Distributions to 0.25

If an administrator could kick CompatHelper back to life, they can generate a PR (and see if there are any other out-of-data dependency compats)

This is needed downstream, for example for Surrogates.jl.

uniform distribution [-1,1]

Hi @timueh,

I tried to se a new default OrthoPoly to denote uniform distribution [-1,1].
I met an issue as follows.

using PolyChaos

Type_for_domain = Tuple{Float64, Float64}
struct Uniform_11Measure <: AbstractCanonicalMeasure
    w::Function
    dom::Type_for_domain
    symmetric::Bool

    function Uniform_11Measure()
        new(x->0.5, (-1,1), true)
    end
end

struct Uniform_11OrthoPoly{V, M, Q} <: AbstractCanonicalOrthoPoly#{V, M, Q}
    deg::Int          # maximum degree
    α::V  # recurrence coefficients
    β::V  # recurrence coefficients
    measure::M
    quad::Q

    # inner constructor
    function Uniform_11OrthoPoly(deg::Int;Nrec::Int=deg+1, addQuadrature::Bool = true)
        #_checkConsistency(deg, Nrec)
        α, β = r_scale(1., rm_legendre(Nrec)...)
        quadrature = addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad()
        new{promote_type(typeof(α), typeof(β)), Uniform_11Measure, typeof(quadrature)}(deg, α, β, Uniform_11Measure(), quadrature)
    end
end

a = Uniform_11OrthoPoly(5);
showpoly(0:5,a)

1
x
x^2 - 0.33
x^3 - 0.6x
x^4 - 0.86x^2 + 0.09
x^5 - 1.11x^3 + 0.24x

It's okay up to now. However, when I examined the recurrence coefficient, I found there might be something wrong. The variance

var([0,1],a)
0.6666666666666666

of a uniform distribution [-1,1] doesn't equal to its correct value 1/3.
I thins it's due to the coefficient,

a.β
2.0                
 0.3333333333333333 
 0.26666666666666666
 0.2571428571428571 
 0.25396825396825395
 0.25252525252525254

where the first element should be equal to 1.
Could we try to fix it?

Add macros for commonly used measures

For instance, op = OrthoPoly("gaussian",10) could be op = @gaussian 10. This really needs to be thought through -- what should be the syntax to distinguish between OrthoPoly and OrthoPolyQ?...

Recurrence coefficients for gamma distribution/laguerre polynomials

Construct bases for rates other than 1? This is currently a kind-of severe limitation.

Moments of order >= 2

probably needs IterTools or something comparable for an elegant solution.

Add tests for recurrence coefficients

We need to test the functions inside recurrence_coefficients_monic.jl.

We should take the functions from the OPQ suite as a reference.

Example in readme latex disagrees with code

The latex is a definite integration, but the example code doesn't mention the integration bounds.

Plotting using Plot.jl and recipes

non-exhaustive list of things that would be interesting to plot:

plot(w::Function,dom::Tuple(Float64,Float64))
plot(m::Measure)
histogram(x::Vector{Float64},op::Union{OrthoPolyQ, MultiOrthoPoly})
histogram(op::Union{OrthoPolyQ, MultiOrthoPoly)

This should be done with Plots.jl and recipes.

caveat: when plotting the w we need to consider cases of infinite support...

Code generation for random ODEs

Starting from a linear time-invariant ODE $\dot{x}(t) = A(\theta) x(t) + B(\theta) u(t)$ with initial condition $x(0) = x_0(\theta)$ and a random parameter $\theta$, auto-generate the Galerkin-projected ODE in terms of the PCE coefficients, similar to the example here where the Galerkin projection is done by hand.

Is evaluate() Automatic differentiable?

Hello everyone,

I am using PolyChaos.jl to build a new Surrogates at Surrogates.jl. Every surrogate is Automatic-differentiable at the moment, so I want to my PolynomialChaosSurrogate to be differentiable as well. However, I think evaluate() is not automatic differentiable, see this example:

PolyChaos.evaluate(collect(val), pcND.ortopolys)

Where val is a tuple like: val = (1.0,2.0,3.0)

I get the following mistake:

ERROR: LoadError: Mutating arrays is not supported
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] (::Zygote.var"#459#460")(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/lib/array.jl:67
 [3] (::Zygote.var"#1009#back#461"{Zygote.var"#459#460"})(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49
 [4] Array at ./array.jl:541 [inlined]
 [5] (::typeof(∂(Array{Float64,2})))(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface2.jl:0
 [6] Array at ./boot.jl:427 [inlined]
 [7] (::typeof(∂(Array{T,2} where T)))(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface2.jl:0
 [8] |> at ./operators.jl:823 [inlined]
 [9] (::typeof(∂(|>)))(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface2.jl:0
 [10] evaluate at /Users/ludovicobessi/.julia/packages/PolyChaos/Dhyx4/src/evaluate.jl:115 [inlined]
 [11] (::typeof(∂(evaluate)))(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface2.jl:0
 [12] (::Zygote.var"#175#176"{typeof(∂(evaluate)),Tuple{Tuple{Nothing,Nothing},Tuple{Nothing,Nothing}}})(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/lib/lib.jl:182
 [13] (::Zygote.var"#347#back#177"{Zygote.var"#175#176"{typeof(∂(evaluate)),Tuple{Tuple{Nothing,Nothing},Tuple{Nothing,Nothing}}}})(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49
 [14] evaluate at /Users/ludovicobessi/.julia/packages/PolyChaos/Dhyx4/src/evaluate.jl:118 [inlined]
 [15] (::typeof(∂(evaluate)))(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface2.jl:0
 [16] evaluate at /Users/ludovicobessi/.julia/packages/PolyChaos/Dhyx4/src/evaluate.jl:123 [inlined]
 [17] (::typeof(∂(evaluate)))(::Array{Float64,2}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface2.jl:0
 [18] PolynomialChaosSurrogate at /Users/ludovicobessi/.julia/dev/Surrogates/src/PolynomialChaos.jl:72 [inlined]
 [19] (::typeof(∂(λ)))(::Float64) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface2.jl:0
 [20] (::Zygote.var"#37#38"{typeof(∂(λ))})(::Float64) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface.jl:45
 [21] gradient(::Function, ::Tuple{Float64,Float64}) at /Users/ludovicobessi/.julia/packages/Zygote/1GXzF/src/compiler/interface.jl:54
 [22] (::var"#39#40")(::Tuple{Float64,Float64}) at /Users/ludovicobessi/.julia/dev/Surrogates/test/AD_compatibility.jl:244
 [23] top-level scope at /Users/ludovicobessi/.julia/dev/Surrogates/test/AD_compatibility.jl:245
 [24] include(::Module, ::String) at ./Base.jl:377
 [25] exec_options(::Base.JLOptions) at ./client.jl:288
 [26] _start() at ./client.jl:484
in expression starting at /Users/ludovicobessi/.julia/dev/Surrogates/test/AD_compatibility.jl:245

Register

@JuliaRegistrator register()

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Extend showbasis to multivariate bases

Error tagging new release

The REQUIRE file could not be found.
cc: @timueh

Add adaptive sparse polynomial chaos

see this paper and the section on it in the UQLab manual.

Sampling Uniform_11OrthoPoly

I am getting many warnings with the following message and Julia freezes indefinitely:

┌ Warning: ignoring keyword method; sampling from Distributions.jl instead
└ @ PolyChaos JULIA_DIR\PolyChaos\THVqe\src\polynomial_chaos.jl:208`

when executing the last command in the following code:

N = 100;
x = sampleMeasure(N, Uniform01OrthoPoly(5));
y = sampleMeasure(N, Uniform_11OrthoPoly(5))

Any ideas on what may be causing errors when sampling Uniform_11OrthoPoly but not when sampling Uniform01OrthoPoly?

Thanks

Error on Uniform01Measure

I tried to use the following command following the Numerical Integration tutorial:
Uniform01Measure(PolyChaos.w_uniform01, (0.0, 1.0), true)
However I get the following error message:

MethodError: no method matching Uniform01Measure(::typeof(w_uniform01), ::Tuple{Float64,Float64}, ::Bool)

Stacktrace:
 [1] top-level scope at In[7]:1
 [2] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

Error tagging new release

The REQUIRE file could not be found.
cc: @timueh

Cannot create orthogonal bases of degree 0

Problem description:
When creating a polynomial basis e.g. by op = GaussOrthoPoly(0) we get a BoundsError.

Expected Behavior:
Returns a polynomial basis type for the 0th-degree polynomial, which is defined as 1.

Fejer quadrature not working for n=1

Add tests for quadrature rules

We need to validate the quadature rules from quadrature_rules.jl.

nonlinear problem

Hi timueh,

Do you have any advice or comments in the current framework if I want to evaluate nonlinear moments, e.g. < exp(\sum f_i \Phi_i), \Phi_k > ?

404 on documentation link

Hello,

I'm interested in using this package for a project of mine; however, when I try and follow the documentation links either through the README or in the about section, they lead to 404s. This is the link I'm attemping to follow:

https://timueh.github.io/PolyChaos.jl/stable/

Did something break in a recent change?

Revise type hierarchy

The current type hierarchy is alright-ish, but not the cleanest-code version we can think of. It'd be cleaner to introduce abstract types and introduce specific types, e.g., for canonical uncertainties.

Something like this

abstract type AbstractMeasure end
abstract type AbstractCanonicalMeasure <: AbstractMeasure end

struct GaussMeasure <: AbstractCanonicalMeasure
    w::Function
    dom::Tuple{Real, Real}
    function GaussMeasure()
        new(gauss_density, (-Inf, Inf))
    end
end

# add other canonical measures below