From 85275e83eb7694cb2aefb2dbaa392d78836399ce Mon Sep 17 00:00:00 2001 From: timueh Date: Wed, 11 Sep 2019 16:13:51 +0200 Subject: [PATCH 1/8] revised tests for auxfuns.jl --- test/auxfuns.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/auxfuns.jl b/test/auxfuns.jl index 0d967245..f02eddf0 100644 --- a/test/auxfuns.jl +++ b/test/auxfuns.jl @@ -1,5 +1,17 @@ using Test, PolyChaos +@test nw(EmptyQuad()) |> typeof == Array{Float64,2} +@test nw(EmptyQuad()) |> length == 0 + +nodes, weights = gauss(rm_legendre01(5)...) +myQuad = Quad("myQuad", length(nodes), nodes, weights) +foo(x) = 6*x^5 +res = 1 +@test isapprox(integrate(foo, nodes, weights), res) +@test isapprox(integrate(foo, myQuad), res) +@test_throws DomainError integrate(foo, EmptyQuad()) +@test isapprox(integrate(foo, Uniform01OrthoPoly(4)), res) + d, nunc = 5, 10 op = genLaguerreOrthoPoly(d,1.23,addQuadrature = false) opq = genLaguerreOrthoPoly(d,1.23) From 9fb2da9cdd555a1d62ce6699fce271de30fec81a Mon Sep 17 00:00:00 2001 From: timueh Date: Wed, 11 Sep 2019 16:14:40 +0200 Subject: [PATCH 2/8] revised tests for densities.jl --- src/densities.jl | 43 ++++++++++++++----------------------------- test/densities.jl | 18 ++++++++++++++++++ 2 files changed, 32 insertions(+), 29 deletions(-) create mode 100644 test/densities.jl diff --git a/src/densities.jl b/src/densities.jl index 12ce2195..70b7706c 100644 --- a/src/densities.jl +++ b/src/densities.jl @@ -1,8 +1,6 @@ export w_legendre, build_w_jacobi, - build_w_jacobi01, w_jacobi, - w_jacobi01, w_laguerre, w_hermite, build_w_genhermite, @@ -17,24 +15,18 @@ export w_legendre, build_w_beta, build_w_gamma +_throwError(t) = throw(DomainError(t, "not in support")) + function w_legendre(t) - -1. <= t <= 1. ? 1. : throw(error("$t not in support")) + -1. <= t <= 1. ? 1. : _throwError(t) end function build_w_jacobi(a,b) return t->w_jacobi(t,a,b) end -function w_jacobi(t,a,b) - -1. <=t<= 1. ? ((1-t)^a*(1+t)^b) : throw(error("$t not in support")) -end - -function build_w_jacobi01(a,b) - @assert a > -1. && b > -1. "Invalid shape parameters" - return t->w_jacobi01(t,a,b) -end -function w_jacobi01(t,a,b) - 0. <=t<= 1. ? ((1-t)^a*t^b) : throw(error("$t not in support")) +function w_jacobi(t,a,b) + -1. <= t <= 1. ? (1-t)^a*(1+t)^b : _throwError(t) end function w_hermite(t) @@ -45,7 +37,7 @@ function build_w_genhermite(mu) return t->w_genhermite(t,mu) end -function w_genhermite(t,μ::Float64) +function w_genhermite(t,μ) abs(t)^(2*μ)*exp(-t^2) end @@ -53,51 +45,44 @@ function build_w_genlaguerre(a) return t->w_genlaguerre(t,a) end -function w_genlaguerre(t,a) - t>=0. ? Float64(t)^a*exp(-Float64(t)) : throw(error("$t not in support")) -end - function w_laguerre(t) - t>=0. ? exp(-t) : throw(error("$t not in support")) + t >= 0. ? exp(-t) : _throwError(t) end function w_meixner_pollaczek(t,lambda,phi) - 1/(2*pi)*exp((2*phi-pi)*t)*abs(gamma(lambda+im*t))^2 + 1 / (2pi) * exp((2*phi - pi)*t) * abs(gamma(lambda+im*t))^2 end function build_w_meixner_pollaczek(lambda,phi) t->w_meixner_pollaczek(t,lambda,phi) end - ################################################## # probability density functions function w_gaussian(t) - 1/(sqrt(2*pi))*exp(-0.5*t^2) + 1 / (sqrt(2*pi))*exp(-0.5*t^2) end function build_w_beta(α,β) return t->w_beta(t,α,β) end + function w_beta(t,α,β) - t^(α-1)*(1-t)^(β-1)/beta(α,β) + -1 <= t <= 1 ? t^(α-1)*(1-t)^(β-1)/beta(α,β) : _throwError(t) end - - function build_w_gamma(α) return t->w_gamma(t,α) end + function w_gamma(t,α) - t>=0. ? (1/gamma(α)*Float64(t)^(α-1)*exp(-t)) : (error("$t not in support")) + t >= 0. ? (1/gamma(α)*Float64(t)^(α-1)*exp(-t)) : _throwError(t) end - function w_uniform01(t) - 0. <=t<= 1. ? 1. : (error("$t not in support")) + 0. <= t <= 1. ? 1. : _throwError(t) end - function w_logistic(t) 0.25*sech(0.5t)^2 end diff --git a/test/densities.jl b/test/densities.jl new file mode 100644 index 00000000..bac76391 --- /dev/null +++ b/test/densities.jl @@ -0,0 +1,18 @@ +using PolyChaos, Test + +outliers = [-1.3, 1.3] + +densities = [w_legendre, + w_uniform01, + build_w_jacobi(1.2,3.4), + build_w_beta(1.2,3.4) +] + +@testset "Supports of densities" begin + for ρ in densities, x in outliers + @test_throws DomainError ρ(x) + end +end + +@test_throws DomainError w_laguerre(-1) +@test_throws DomainError build_w_gamma(1)(-1) \ No newline at end of file From ab2ce6acccd17bdb5a515a007df87113fe8ea23e Mon Sep 17 00:00:00 2001 From: timueh Date: Wed, 11 Sep 2019 17:29:09 +0200 Subject: [PATCH 3/8] revised quadrature tests --- src/multiIndices.jl | 25 +++++++++++++------------ src/quadrature_rules.jl | 29 ++++++++++++++++------------- test/densities.jl | 2 ++ test/quadrature_rules.jl | 15 +++++++++++++++ test/runtests.jl | 1 + 5 files changed, 47 insertions(+), 25 deletions(-) diff --git a/src/multiIndices.jl b/src/multiIndices.jl index a50aa397..2fe1c701 100644 --- a/src/multiIndices.jl +++ b/src/multiIndices.jl @@ -1,4 +1,5 @@ -export calculateMultiIndices, findUnivariateIndices, calculateMultiIndices_interaction +export calculateMultiIndices, + findUnivariateIndices function calculateMultiIndices(d::Int, n::Int) # d denotes dimension of random variables/number of sources of uncertainty, @@ -133,14 +134,14 @@ function findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64} end ################################################################# -function calculateMultiIndices_interaction(nξ::Int,deg::Int,j::Int,p::Int) - inds = calculateMultiIndices(nξ,deg) - get_interaction(inds,j,p) -end - -function get_interaction(inds::Matrix{Int},j::Int,p::Int) - j < 0 && throw(error("interaction order must be non-negative")) - j > size(inds,2) && throw(error("interaction order cannot be greater than number of uncertainties")) - Iterators.filter(x -> count(!iszero,x) == j && sum(x) == p, eachrow(inds)) - # alternative -> collect -end +# function calculateMultiIndices_interaction(nξ::Int,deg::Int,j::Int,p::Int) +# inds = calculateMultiIndices(nξ,deg) +# get_interaction(inds,j,p) +# end + +# function get_interaction(inds::Matrix{Int},j::Int,p::Int) +# j < 0 && throw(error("interaction order must be non-negative")) +# j > size(inds,2) && throw(error("interaction order cannot be greater than number of uncertainties")) +# Iterators.filter(x -> count(!iszero,x) == j && sum(x) == p, eachrow(inds)) +# # alternative -> collect +# end diff --git a/src/quadrature_rules.jl b/src/quadrature_rules.jl index 3a366bd5..b1db43af 100644 --- a/src/quadrature_rules.jl +++ b/src/quadrature_rules.jl @@ -96,13 +96,13 @@ function golubwelsch(α::Vector{<:Real}, β::Vector{<:Real}, maxiter::Int=30) a[idx], (β0*w.^2)[idx] end -golubwelsch(op::OrthoPoly) = golubwelsch(op.α,op.β) +golubwelsch(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = golubwelsch(op.α,op.β) """ gauss(N::Int,α::Vector{<:Real},β::Vector{<:Real}) gauss(α::Vector{<:Real},β::Vector{<:Real}) - gauss(N::Int,op::OrthoPoly) - gauss(op::OrthoPoly) + gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) + gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) Gauss quadrature rule, also known as Golub-Welsch algorithm `gauss()` generates the `N` Gauss quadrature nodes and weights for a given weight function. @@ -124,15 +124,18 @@ function gauss(N::Int,α::Vector{<:Real},β::Vector{<:Real}) @assert N0 >= N "not enough recurrence coefficients" @inbounds golubwelsch(α[1:N],β[1:N]) end + gauss(α::Vector{<:Real},β::Vector{<:Real}) = gauss(length(α)-1,α,β) -gauss(N::Int,op::OrthoPoly) = gauss(N::Int,op.α,op.β) -gauss(op::OrthoPoly) = gauss(op.α,op.β) + +gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = gauss(N::Int,op.α,op.β) + +gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly}) = gauss(op.α,op.β) """ radau(N::Int,α::Vector{<:Real},β::Vector{<:Real},end0::Real) radau(α::Vector{<:Real},β::Vector{<:Real},end0::Real) - radau(N::Int,op::OrthoPoly,end0::Real) - radau(op::OrthoPoly,end0::Real) + radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) + radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) Gauss-Radau quadrature rule. Given a weight function encoded by the recurrence coefficients `(α,β)`for the associated orthogonal polynomials, the function generates the @@ -164,14 +167,14 @@ function radau(N::Int,α_::Vector{<:Real},β::Vector{<:Real},end0::Real) gauss(N+1,α,β) end radau(α::Vector{<:Real},β::Vector{<:Real},end0::Real) = radau(length(α)-2,α,β,end0) -radau(N::Int,op::OrthoPoly,end0::Real) = radau(N,op.α,op.β,end0) -radau(op::OrthoPoly,end0::Real) = radau(op.α,op.β,end0::Real) +radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) = radau(N,op.α,op.β,end0) +radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real) = radau(op.α,op.β,end0::Real) """ lobatto(N::Int,α::Vector{<:Real},β::Vector{<:Real},endl::Real,endr::Real) lobatto(α::Vector{<:Real},β::Vector{<:Real},endl::Real,endr::Real) - lobatto(N::Int,op::OrthoPoly,endl::Real,endr::Real) - lobatto(op::OrthoPoly,endl::Real,endr::Real) + lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) + lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) Gauss-Lobatto quadrature rule. Given a weight function encoded by the recurrence coefficients for the associated orthogonal polynomials, the function generates @@ -209,5 +212,5 @@ function lobatto(N::Int,α_::Vector{<:Real},β_::Vector{<:Real},endl::Real,endr: gauss(N+2,α,β) end lobatto(α::Vector{<:Real},β::Vector{<:Real},endl::Real,endr::Real) = lobatto(length(α)-3,α,β,endl,endr) -lobatto(N::Int,op::OrthoPoly,endl::Real,endr::Real) = lobatto(N,op.α,op.β,endl,endr) -lobatto(op::OrthoPoly,endl::Real,endr::Real) = lobatto(op.α,op.β,endl,endr) +lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) = lobatto(N,op.α,op.β,endl,endr) +lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real) = lobatto(op.α,op.β,endl,endr) diff --git a/test/densities.jl b/test/densities.jl index bac76391..65443464 100644 --- a/test/densities.jl +++ b/test/densities.jl @@ -2,6 +2,8 @@ using PolyChaos, Test outliers = [-1.3, 1.3] +@test_throws DomainError PolyChaos._throwError(3) + densities = [w_legendre, w_uniform01, build_w_jacobi(1.2,3.4), diff --git a/test/quadrature_rules.jl b/test/quadrature_rules.jl index 4af1e706..43aa1299 100644 --- a/test/quadrature_rules.jl +++ b/test/quadrature_rules.jl @@ -123,3 +123,18 @@ intervals = [(-5.,5.),(-Inf,Inf),(-Inf,-1.),(1.,Inf),(-3.,5.),(-10.,1.)] @test isapprox(exp.(-abs.(nod2))'*wei2-analytic,0;atol=1e-4) end end + +degree = 9 +N = degree + 1 +α, β = rm_legendre01(degree + 1) +op = Uniform01OrthoPoly(degree) +endl, endr = 0, 1 + +myMeas = Measure("myMeasure", t->1, (0, 1), true) +op_num = OrthoPoly("myOp", degree, myMeas) + +@testset "OrthoPoly overloading" begin + @test gauss(N-1, α, β) == gauss(α, β) == gauss(N-1, op) == gauss(op) + @test radau(N-2, α, β, endr) == radau(α, β, endr) == radau(N-2, op, endr) == radau(op, endr) + @test lobatto(N-3, α, β, endl, endr) == lobatto(α, β, endl, endr) == lobatto(N-3, op, endl, endr) == lobatto(op, endl, endr) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 0d2d2a76..014192ce 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,4 +7,5 @@ include("quadrature_rules.jl") include("show.jl") include("types.jl") include("auxfuns.jl") +include("densities.jl") include("polynomial_chaos.jl") From 475f8446964231c0e25ee0ff9e1a9609f0394467 Mon Sep 17 00:00:00 2001 From: timueh Date: Thu, 12 Sep 2019 16:34:45 +0200 Subject: [PATCH 4/8] revised tests for evaluate.jl --- test/evaluate.jl | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/test/evaluate.jl b/test/evaluate.jl index 506753c6..d14e3308 100644 --- a/test/evaluate.jl +++ b/test/evaluate.jl @@ -5,11 +5,23 @@ op = genLaguerreOrthoPoly(deg,1.23) n = 10 X = [ exp(3 + rand()), exp.(3 .+ rand(n)) ] + @testset "evaluate univariate basis at point(s)" begin for x in X @test evaluate(collect(0:deg),x,op.α,op.β) == evaluate(collect(0:deg),x,op) for d in 0:deg - @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op) + if typeof(x) <: Real + @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op) == evaluate(d,[x],op.α,op.β) + else + @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op) + end + end + + res = hcat(map(i->evaluate(i,x,op.α,op.β),collect(0:op.deg))...) + if typeof(x) <: Real + @test evaluate(x,op) == evaluate([x], op) == evaluate(collect(0:op.deg),[x],op) == evaluate(collect(0:op.deg),[x],op.α,op.β) == res + else + @test evaluate(x, op) == evaluate(collect(0:op.deg),x,op) == evaluate(collect(0:op.deg),x,op.α,op.β) == res end end end @@ -29,6 +41,12 @@ X = [ exp.(3 .+ rand(nop)) , exp.(3 .+ rand(n,nop)) ] @testset "evaluate multivariate basis at point(s)" begin for x in X @test evaluate(mop.ind,x,coeffs(mop)...) == evaluate(mop.ind,x,mop) == evaluate(x,mop) + if typeof(x) <: Vector{<:Real} + @test evaluate(mop.ind,x,coeffs(mop)...) == evaluate(mop.ind,x,mop) == evaluate(x,mop) == evaluate(mop.ind,reshape(x,1,length(x)),coeffs(mop)...) + @test evaluate(x,mop) == evaluate(mop.ind,reshape(x,1,length(x)),mop) == evaluate(mop.ind,reshape(x,1,length(x)),coeffs(mop)...) + else + @test evaluate(x,mop) == evaluate(mop.ind,x,mop) == evaluate(mop.ind,x,coeffs(mop)...) + end for d in eachrow(mop.ind[1:3,:]) @test evaluate(collect(d),x,coeffs(mop)...) == evaluate(collect(d),x,mop) end From 16fe8cd4deac9fb90cdab4e4aa9e0e21fe3eba93 Mon Sep 17 00:00:00 2001 From: timueh Date: Thu, 12 Sep 2019 16:49:57 +0200 Subject: [PATCH 5/8] revised test for recurrence coefficients, cleaned up code --- src/recurrence_coefficients_monic.jl | 75 +++++++++++++++------------- src/sparse_basis.jl | 60 ++++++++++++++++++++++ src/stieltjes.jl | 25 ---------- test/recurrence_coefficients.jl | 5 ++ 4 files changed, 105 insertions(+), 60 deletions(-) create mode 100644 src/sparse_basis.jl diff --git a/src/recurrence_coefficients_monic.jl b/src/recurrence_coefficients_monic.jl index aa0090d6..755b6419 100644 --- a/src/recurrence_coefficients_monic.jl +++ b/src/recurrence_coefficients_monic.jl @@ -194,17 +194,49 @@ that are orthogonal on ``(0,1)`` relative to ``w(t) = 1``. rm_legendre01(N::Int) = rm_jacobi01(N) -function rm_chebyshev1(N::Int) - @assert N>=0 "N has to be non-negative" - α = zeros(Float64,N) - if N == 1 - return α, [pi] - elseif N == 2 - return α, [pi; 0.5] +""" + rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real) + rm_meixner_pollaczek(N::Int,lambda::Real) + + Creates `N` recurrence coefficients for monic + Meixner-Pollaczek polynomials with parameters λ and ϕ. These are orthogonal on + ``[-\\infty,\\infty]`` relative to the weight function ``w(t)=(2 \\pi)^{-1} \\exp{(2 \\phi-\\pi)t} |\\Gamma(\\lambda+ i t)|^2``. + + The call `rm_meixner_pollaczek(n,lambda)` is the same as `rm_meixner_pollaczek(n,lambda,pi/2)`. +""" +function rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real) + @assert N>=0 && lambda>0. && phi>0. "parameter(s) out of range" + N == 0 && return Array{Float64,1}(undef,0), Array{Float64,1}(undef,0) + n=1:N; + sinphi=sin(phi); lam2=2*lambda; + ab=zeros(Float64,N,2) + if sinphi==1 + ab[:,1]=zeros(N); else - return α, pushfirst!(pushfirst!(0.25*ones(N-2),0.5),pi) + ab[:,1]= -(n .+ (lambda-1))/tan(phi); end + ab[:,2]= (n .- 1).*(n .+ (lam2-2))/(4*sinphi^2); + ab[1,2]=gamma(lam2)/(2*sinphi)^lam2; + return ab[:,1], ab[:,2] end +rm_meixner_pollaczek(N::Int,lambda::Real) = rm_meixner_pollaczek(N,lambda,pi/2) + + +################################################################### +################################################################### +################################################################### + +# function rm_chebyshev1(N::Int) +# @assert N>=0 "N has to be non-negative" +# α = zeros(Float64,N) +# if N == 1 +# return α, [pi] +# elseif N == 2 +# return α, [pi; 0.5] +# else +# return α, pushfirst!(pushfirst!(0.25*ones(N-2),0.5),pi) +# end +# end # rm_hahn(N::Int,a::Real,b::Real) # rm_hahn(N::Int,a::Real) # rm_hahn(N::Int) @@ -248,30 +280,3 @@ end # end # rm_hahn(N::Int,a::Real) = rm_hahn(N,a,a) # rm_hahn(N::Int) = rm_hahn(N,0.,0.) - -""" - rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real) - rm_meixner_pollaczek(N::Int,lambda::Real) - - Creates `N` recurrence coefficients for monic - Meixner-Pollaczek polynomials with parameters λ and ϕ. These are orthogonal on - ``[-\\infty,\\infty]`` relative to the weight function ``w(t)=(2 \\pi)^{-1} \\exp{(2 \\phi-\\pi)t} |\\Gamma(\\lambda+ i t)|^2``. - - The call `rm_meixner_pollaczek(n,lambda)` is the same as `rm_meixner_pollaczek(n,lambda,pi/2)`. -""" -function rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real) - @assert N>=0 && lambda>0. && phi>0. "parameter(s) out of range" - N == 0 && return Array{Float64,1}(undef,0), Array{Float64,1}(undef,0) - n=1:N; - sinphi=sin(phi); lam2=2*lambda; - ab=zeros(Float64,N,2) - if sinphi==1 - ab[:,1]=zeros(N); - else - ab[:,1]= -(n .+ (lambda-1))/tan(phi); - end - ab[:,2]= (n .- 1).*(n .+ (lam2-2))/(4*sinphi^2); - ab[1,2]=gamma(lam2)/(2*sinphi)^lam2; - return ab[:,1], ab[:,2] -end -rm_meixner_pollaczek(N::Int,lambda::Real) = rm_meixner_pollaczek(N,lambda,pi/2) diff --git a/src/sparse_basis.jl b/src/sparse_basis.jl new file mode 100644 index 00000000..597bc183 --- /dev/null +++ b/src/sparse_basis.jl @@ -0,0 +1,60 @@ +# In later versions of `PolyChaos.jl` we include basis-adaptive sparse PCE +# This file is a first take on this + +export orthosparse + +function regress(M,y) + c = pinv(M)*y + c, M*c +end + +function lack_of_fit(y::Vector{Float64},x::Vector{Float64},ymod::Vector{Float64}) + sum(δy^2 for δy in y - ymod) / sum(δy^2 for δy in y .- mean(y)) +end + +function orthosparse(y::Vector{Float64},x::Vector{Float64},name::String,p_max::Int64;ε_forward::Float64,ε_backward::Float64,accuracy::Float64) + p_index = 1:p_max + final_index = Int64[] + op = OrthoPoly(name,p_max) + y_bar = mean(y) + A, A_plus = [0], [0] + for i in p_index + i >= p_max && break + ########### FORWARD STEP + push!(A,i) + Φ = evaluate(A,x,op) + a_hat, ymod = regress(Φ,y) + R2 = 1 - lack_of_fit(y,x,ymod) + R2 <= ε_forward ? A_plus = filter(x -> x != i, A) : nothing + + ########### Backward STEP + R2_array = Float64[] + for b in A_plus[1:end-1] + Φ = evaluate(filter(x -> x ≠ b, A_plus),x,op) + a_tmp, ymod = regress(Φ,y) + r2 = 1 - lack_of_fit(y,x,ymod) + push!(R2_array,r2) + end + s::Int64 = 1 + for r2 in R2_array + if abs(R2 - r2) < ε_backward + filter!(x -> x ≠ A_plus[s], A_plus) + s -= 1 + end + s += 1 + end + A = A_plus + # CHECK TERMINATION CRITERION + if length(A) != 0 + Φ = evaluate(A,x,op) + a_tmp, ymod = regress(Φ,y) + r2 = 1 - lack_of_fit(y,x,ymod) + if r2 >= accuracy + println("Accuracy achieved. Breaking.\n") + return A + end + end + ##################################### + end + error("Algorithm terminated early; perhaps a pathological problem was provided.") +end diff --git a/src/stieltjes.jl b/src/stieltjes.jl index 9ee1989c..9039d8fe 100644 --- a/src/stieltjes.jl +++ b/src/stieltjes.jl @@ -88,31 +88,6 @@ function lanczos(N::Int64,nodes::Vector{Float64},weights::Vector{Float64};remove @inbounds p0[Base.OneTo(N)], p1[Base.OneTo(N)] end - -function mcdis_twologistics(n::Int,p1::Vector{Float64},p2::Vector{Float64};Mmax::Int=100,eps0::Real=1e-9) - M0 = n - Mcap = 0 - Mi = M0 - α, β = zeros(n), zeros(n) - for i=1:Mmax - nai, wai = quadrature_logistic(Mi,p1[1],p1[2]) - nbi, wbi = quadrature_logistic(Mi,p2[1],p2[2]) - wai, wbi = 0.5*wai, 0.5*wbi - ni, wi = [nai; nbi], [wai; wbi] - α_, β_ = stieltjes(n,ni,wi) - err = maximum( abs(β_[k]-β[k])/β_[k] for k=1:n ) - display("err = $err") - if err<=eps0 - Mcap = i-1 - return α_, β_ - else - α, β, = α_, β_ - i==1 ? Mi = M0+1 : Mi = Int(2^(floor((i-1)/5))*n) - end - end - warn("Algorithm did not terminate after $Mmax iterations.") -end - """ mcdiscretization(N::Int64,quads::Vector{},discretemeasure::Matrix{Float64}=zeros(0,2);discretization::Function=stieltjes,Nmax::Integer=300,ε::Float64=1e-8,gaussquad::Bool=false) This routine returns ``N`` recurrence coefficients of the polynomials that are diff --git a/test/recurrence_coefficients.jl b/test/recurrence_coefficients.jl index 451e5892..001328f0 100644 --- a/test/recurrence_coefficients.jl +++ b/test/recurrence_coefficients.jl @@ -47,6 +47,7 @@ end @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol) αβcom = coeffs(JacobiOrthoPoly(n-1, al, be, addQuadrature=false)) @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol) + @test rm_jacobi(n,al) == rm_jacobi(n,al,al) close(myfile) end end @@ -58,6 +59,7 @@ end myfile = open("dataRecCoeffs/jac01$(n)al$(al)be$(be).txt") αβref = parse.(Float64,readlines(myfile)) αβcom = rm_jacobi01(n,al,be) + @test rm_jacobi01(n,al) == rm_jacobi01(n,al,al) @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol) close(myfile) end @@ -89,6 +91,9 @@ meipol = 0.1:0.2:2 αβcom = coeffs(MeixnerPollaczekOrthoPoly(n-1, lambda, phi, addQuadrature=false)) @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol) close(myfile) + @test rm_meixner_pollaczek(n,lambda) == rm_meixner_pollaczek(n,lambda,pi/2) + @test rm_meixner_pollaczek(n,lambda,pi/2)[1] == zeros(n) end + end end From d6a010bf285e52498b077579f9b11db6180020c1 Mon Sep 17 00:00:00 2001 From: timueh Date: Thu, 12 Sep 2019 17:57:53 +0200 Subject: [PATCH 6/8] tests for multi_indices.jl + renaming --- src/PolyChaos.jl | 4 +- src/multi_indices.jl | 79 ++++++++++++++++++++++ src/scalar_product.jl | 149 ++++++++++++++++++++++++++++++++++++++++++ src/sparse_basis.jl | 60 ----------------- test/multi_indices.jl | 26 ++++++++ test/runtests.jl | 1 + 6 files changed, 257 insertions(+), 62 deletions(-) create mode 100644 src/multi_indices.jl create mode 100644 src/scalar_product.jl delete mode 100644 src/sparse_basis.jl create mode 100644 test/multi_indices.jl diff --git a/src/PolyChaos.jl b/src/PolyChaos.jl index 1e3ece82..4a967900 100644 --- a/src/PolyChaos.jl +++ b/src/PolyChaos.jl @@ -16,14 +16,14 @@ include("typesMeasures.jl") include("typesQuad.jl") include("typesOrthoPolys.jl") include("typesTensor.jl") -include("multiIndices.jl") +include("multi_indices.jl") include("evaluate.jl") include("quadrature_rules.jl") include("recurrence_coefficients_monic.jl") include("stieltjes.jl") include("densities.jl") include("auxfuns.jl") -include("scalarproduct.jl") +include("scalar_product.jl") include("tensor.jl") include("show.jl") include("polynomial_chaos.jl") diff --git a/src/multi_indices.jl b/src/multi_indices.jl new file mode 100644 index 00000000..1540b583 --- /dev/null +++ b/src/multi_indices.jl @@ -0,0 +1,79 @@ +export calculateMultiIndices, + findUnivariateIndices + +function calculateMultiIndices(d::Int, n::Int) + # d denotes dimension of random variables/number of sources of uncertainty, + # n the maximum degree of multivariate basis + # function to calculate indices of multivariate basis following the algorithm + # from "Spectral Methods for Uncertainty Quantification; Le Maitre, Knio;2014" p.516-517 + # No::BigInt = factorial(BigInt(d+n))/(factorial(Bigint(d))*factorial(BigInt(n))); #number of polynomials of multivariate basis + n < 0 && throw(DomainError(n, "maximum degree must be non-negative")) + d <= 0 && throw(DomainError(d, "number of uncertainties must be positive")) + # catch case n == 0 --> No-d==0 + n == 0 && return zeros(Int64,1,d) + # non-pathological cases begin here + No = numberPolynomials(d,n) + inds = vcat(zeros(Int64,1,d),Matrix(1I,d,d),zeros(Int64, No-d-1, d)); #initiate index matrix for basis + pi=ones(Int64,No,d); + + for k=2:No + g = 0; + for l=1:d + pi[k,l]=sum(pi[k-1,:])-g; + g=g+pi[k-1,l]; + end + end + + P=d+1; + for k=2:n + L = P; + for j=1:d, m=L-pi[k,j]+1:L + P=P+1; + inds[P,:]=inds[m,:]; + inds[P,j]=inds[P,j]+1; + end + end + + return inds +end + +""" + computes the number of polynomials with a multivariate basis + `(d+n)!/(d!+n!)` +""" + +function numberPolynomials(d::Int64, n::Int64) + x, y = max(d,n), min(d,n) + return UInt128(prod(UInt128(x+1):UInt128(d+n)) ÷ factorial(UInt128(y))) +end + +""" + findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64} +Given the multi-index `ind` this function returns all entries of the multivariate basis +that correspond to the `i`th univariate basis. +""" +function findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64} + l,p = size(ind) + @assert i<=p "basis is $p-variate, you requested $i-variate" + deg = ind[end,end] + @assert deg>=0 "invalid degree" + col = ind[:,i] + myind = zeros(Int64,deg) + for deg_ = 1:deg + myind[deg_] = findfirst(x->x==deg_,col) + end + pushfirst!(myind,1) +end + +################################################################# +# function calculateMultiIndices_interaction(nξ::Int,deg::Int,j::Int,p::Int) +# inds = calculateMultiIndices(nξ,deg) +# get_interaction(inds,j,p) +# end + +# function get_interaction(inds::Matrix{Int},j::Int,p::Int) +# j < 0 && throw(error("interaction order must be non-negative")) +# j > size(inds,2) && throw(error("interaction order cannot be greater than number of uncertainties")) +# Iterators.filter(x -> count(!iszero,x) == j && sum(x) == p, eachrow(inds)) +# # alternative -> collect +# end diff --git a/src/scalar_product.jl b/src/scalar_product.jl new file mode 100644 index 00000000..235d832a --- /dev/null +++ b/src/scalar_product.jl @@ -0,0 +1,149 @@ +export computeSP, + computeSP2 + +function computeSP(a_::Vector{<:Int}, + α::Vector{<:Vector{<:Real}},β::Vector{<:Vector{<:Real}}, + nodes::Vector{<:Vector{<:Real}},weights::Vector{<:Vector{<:Real}}, + ind::Matrix{<:Int}; + issymmetric::BitArray=falses(length(α)), + zerotol::Float64=1e-10) + minimum(a_) < 0 && throw(DomainError(minimum(a_),"no negative degrees allowed")) + l, p = size(ind) # p-variate basis + p == 1 && computeSP(a_,α[1],β[1],nodes[1],weights[1];issymmetric=issymmetric[1]) + l -= 1 + maximum(a_) > l && throw(DomainError(maximum(a_), "not enough elements in multi-index (requested: $(maximum(a)), max: $l)")) + !(length(α)==length(β)==length(nodes)==length(weights)==length(issymmetric)==p) && throm(InconsistencyError("inconsistent number of rec. coefficients and/or nodes/weights")) + + a = Vector{Int64}() + + for aa in a_ + !iszero(aa) && push!(a,aa) + end + + if iszero(length(a)) + return prod(β[i][1] for i in 1:p) + elseif isone(length(a)) + return 0. + else + inds_uni = multi2uni(a,ind) + val = 1. + for i in 1:p + # compute univariate scalar product + @inbounds v = computeSP(inds_uni[i,:],α[i],β[i],nodes[i],weights[i];issymmetric=issymmetric[i]) + isapprox(v,0,atol=zerotol) ? (return 0.) : (val*=v) + end + end + return val +end + +function computeSP(a::Vector{<:Int}, op::Vector{<:AbstractOrthoPoly}, ind::Matrix{<:Int}) + α, β = coeffs(op) + nodes, weights = nw(op) + computeSP(a, α, β, nodes, weights, ind; issymmetric=issymmetric.(op)) +end + +computeSP(a::Vector{<:Int},mop::MultiOrthoPoly) = computeSP(a, mop.uni, mop.ind) + +""" +__Univariate__ +``` +computeSP(a::Vector{<:Int},α::Vector{<:Real},β::Vector{<:Real},nodes::Vector{<:Real},weights::Vector{<:Real};issymmetric::Bool=false) +computeSP(a::Vector{<:Int},op::AbstractOrthoPoly;issymmetric=issymmetric(op)) +``` +__Multivariate__ +``` +computeSP( a::Vector{<:Int}, + α::Vector{<:Vector{<:Real}},β::Vector{<:Vector{<:Real}}, + nodes::Vector{<:Vector{<:Real}},weights::Vector{<:Vector{<:Real}}, + ind::Matrix{<:Int}; + issymmetric::BitArray=falses(length(α))) +computeSP(a::Vector{<:Int},op::Vector{<:AbstractOrthoPoly},ind::Matrix{<:Int}) +computeSP(a::Vector{<:Int},mOP::MultiOrthoPoly) +``` + +Computes the scalar product +```math +\\langle \\phi_{a_1},\\phi_{a_2},\\cdots,\\phi_{a_n} \\rangle, +``` +where `n = length(a)`. +This requires to provide the recurrence coefficients `(α,β)` and the quadrature rule +`(nodes,weights)`, as well as the multi-index `ind`. +If provided via the keyword `issymmetric`, symmetry of the weight function is exploited. +All computations of the multivariate scalar products resort back to efficient computations +of the univariate scalar products. +Mathematically, this follows from Fubini's theorem. + +The function is dispatched to facilitate its use with `AbstractOrthoPoly` and its quadrature rule `Quad`. + +!!! note + - Zero entries of ``a`` are removed automatically to simplify computations, which follows from + ```math + \\langle \\phi_i, \\phi_j, \\phi_0,\\cdots,\\phi_0 \\rangle = \\langle \\phi_i, \\phi_j \\rangle, + ``` + because `\\phi_0 = 1`. + - It is checked whether enough quadrature points are supplied to solve the integral exactly. +""" +function computeSP(a_::Vector{<:Int},α::Vector{<:Real},β::Vector{<:Real},nodes::Vector{<:Real},weights::Vector{<:Real};issymmetric::Bool=false) + minimum(a_) < 0 && throw(DomainError(minimum(a_), "no negative degrees allowed")) + # this works because ``<Φ_i,Φ_j,Φ_0,...,Φ_0> = <Φ_i,Φ_j>`` + # hence, avoiding unnecessary operations that introduce numerical gibberish + # a = Vector{Int64}() + # for aa in a_ + # !iszero(aa) && push!(a,aa) + # end + a = filter(x -> x > 0, a_) + # simplification in case the measure is symmetric w.r.t t=0 + # exploit symmetry of the density, see Theorem 1.17 in W. Gautschi's book + # and exploit the fact that + (issymmetric && isodd(sum(a))) && return 0. + # Gauss quadrature rules have exactness 2N-1 + !(Int(ceil(0.5*(sum(a)+1)))<=length(nodes)) && throw(InconsistencyError("not enough nodes to integrate exactly ($(length(nodes)) provided, where $(Int(ceil(0.5*(sum(a)+1)))) are needed)")) + + res = + if iszero(length(a)) + first(β) + elseif isone(length(a)) + 0. + elseif length(a) == 2 + @inbounds a[1] == a[2] ? computeSP2(first(a),β)[end] : 0. + else + f = ones(Float64,length(nodes)) + @simd for i = 1:length(a) + @inbounds f = f.*evaluate(a[i],nodes,α,β) + end + dot(weights,f) + end +end + +function computeSP(a::Vector{<:Int},op::AbstractOrthoPoly) + computeSP(a,op.α,op.β,op.quad.nodes,op.quad.weights;issymmetric=issymmetric(op)) +end + +""" + computeSP2(n::Int,β::Vector{<:Real}) + computeSP2(n::Int,op::AbstractOrthoPoly) = computeSP2(n,op.β) + computeSP2(op::AbstractOrthoPoly) = computeSP2(op.deg,op.β) + +Computes the `n` *regular* scalar products aka 2-norms of the orthogonal polynomials, namely +```math +\\|ϕ_i\\|^2 = \\langle \\phi_i,\\phi_i\\rangle \\quad \\forall i \\in \\{ 0,\\dots,n \\}. +``` +Notice that only the values of `β` of the recurrence coefficients `(α,β)` are required. +The computation is based on equation (1.3.7) from Gautschi, W. "Orthogonal Polynomials: Computation and Approximation". +Whenever there exists an analytic expressions for `β`, this function should be used. + +The function is multiply dispatched to facilitate its use with `AbstractOrthoPoly`. +""" +function computeSP2(n::Int,β::Vector{<:Real}) + n < 0 && throw(DomainError(n, "can only compute scalar products for non-negative degrees")) + length(β) < n + 1 && throw(InconsistencyError("inconsistent length of β")) + n == 0 && return β[1] + s = ones(Float64,n+1) + @inbounds s[1] = β[1] # order 0computeSP2 + for i in 2:n+1 + @inbounds s[i] = s[i-1]*β[i] + end + return s +end +computeSP2(n::Int,op::AbstractOrthoPoly) = computeSP2(n,op.β) +computeSP2(op::AbstractOrthoPoly) = computeSP2(op.deg,op.β) \ No newline at end of file diff --git a/src/sparse_basis.jl b/src/sparse_basis.jl deleted file mode 100644 index 597bc183..00000000 --- a/src/sparse_basis.jl +++ /dev/null @@ -1,60 +0,0 @@ -# In later versions of `PolyChaos.jl` we include basis-adaptive sparse PCE -# This file is a first take on this - -export orthosparse - -function regress(M,y) - c = pinv(M)*y - c, M*c -end - -function lack_of_fit(y::Vector{Float64},x::Vector{Float64},ymod::Vector{Float64}) - sum(δy^2 for δy in y - ymod) / sum(δy^2 for δy in y .- mean(y)) -end - -function orthosparse(y::Vector{Float64},x::Vector{Float64},name::String,p_max::Int64;ε_forward::Float64,ε_backward::Float64,accuracy::Float64) - p_index = 1:p_max - final_index = Int64[] - op = OrthoPoly(name,p_max) - y_bar = mean(y) - A, A_plus = [0], [0] - for i in p_index - i >= p_max && break - ########### FORWARD STEP - push!(A,i) - Φ = evaluate(A,x,op) - a_hat, ymod = regress(Φ,y) - R2 = 1 - lack_of_fit(y,x,ymod) - R2 <= ε_forward ? A_plus = filter(x -> x != i, A) : nothing - - ########### Backward STEP - R2_array = Float64[] - for b in A_plus[1:end-1] - Φ = evaluate(filter(x -> x ≠ b, A_plus),x,op) - a_tmp, ymod = regress(Φ,y) - r2 = 1 - lack_of_fit(y,x,ymod) - push!(R2_array,r2) - end - s::Int64 = 1 - for r2 in R2_array - if abs(R2 - r2) < ε_backward - filter!(x -> x ≠ A_plus[s], A_plus) - s -= 1 - end - s += 1 - end - A = A_plus - # CHECK TERMINATION CRITERION - if length(A) != 0 - Φ = evaluate(A,x,op) - a_tmp, ymod = regress(Φ,y) - r2 = 1 - lack_of_fit(y,x,ymod) - if r2 >= accuracy - println("Accuracy achieved. Breaking.\n") - return A - end - end - ##################################### - end - error("Algorithm terminated early; perhaps a pathological problem was provided.") -end diff --git a/test/multi_indices.jl b/test/multi_indices.jl new file mode 100644 index 00000000..f9b36799 --- /dev/null +++ b/test/multi_indices.jl @@ -0,0 +1,26 @@ +using Test, PolyChaos + +A = [ 0 0 0 0; + 1 0 0 0; + 0 1 0 0; + 0 0 1 0; + 0 0 0 1; + 2 0 0 0; + 1 1 0 0; + 1 0 1 0; + 1 0 0 1; + 0 2 0 0; + 0 1 1 0; + 0 1 0 1; + 0 0 2 0; + 0 0 1 1; + 0 0 0 2] + +@test calculateMultiIndices(4, 2) == A +@test calculateMultiIndices(3, 0) == zeros(Int64, 1, 3) +@test_throws DomainError calculateMultiIndices( 3, -1) +@test_throws DomainError calculateMultiIndices( 0, 1) +@test_throws DomainError calculateMultiIndices(-1, 1) + +@test PolyChaos.numberPolynomials(3,4) == factorial(3+4) ÷ ( factorial(3)*factorial(4) ) +@test findUnivariateIndices(2,A) == [1,3,10] \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 014192ce..cb45805b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,3 +9,4 @@ include("types.jl") include("auxfuns.jl") include("densities.jl") include("polynomial_chaos.jl") +include("multi_indices.jl") \ No newline at end of file From eb671416e24c35bf8f48bfb88521c00d96df8ff0 Mon Sep 17 00:00:00 2001 From: timueh Date: Thu, 12 Sep 2019 17:59:35 +0200 Subject: [PATCH 7/8] deleted superfluous files --- src/multiIndices.jl | 147 ------------------------------------------ src/scalarproduct.jl | 149 ------------------------------------------- 2 files changed, 296 deletions(-) delete mode 100644 src/multiIndices.jl delete mode 100644 src/scalarproduct.jl diff --git a/src/multiIndices.jl b/src/multiIndices.jl deleted file mode 100644 index 2fe1c701..00000000 --- a/src/multiIndices.jl +++ /dev/null @@ -1,147 +0,0 @@ -export calculateMultiIndices, - findUnivariateIndices - -function calculateMultiIndices(d::Int, n::Int) - # d denotes dimension of random variables/number of sources of uncertainty, - # n the maximum degree of multivariate basis - # function to calculate indices of multivariate basis following the algorithm - # from "Spectral Methods for Uncertainty Quantification; Le Maitre, Knio;2014" p.516-517 - # No::BigInt = factorial(BigInt(d+n))/(factorial(Bigint(d))*factorial(BigInt(n))); #number of polynomials of multivariate basis - n < 0 && throw(error("maximum degree must be non-negative")) - d <= 0 && throw(error("number of uncertainties must be positive")) - # catch case n == 0 --> No-d==0 - n == 0 && return zeros(Int64,1,d) - # non-pathological cases begin here - No = numberPolynomials(d,n) - inds = vcat(zeros(Int64,1,d),Matrix(1I,d,d),zeros(Int64, No-d-1, d)); #initiate index matrix for basis - pi=ones(Int64,No,d); - - for k=2:No - g = 0; - for l=1:d - pi[k,l]=sum(pi[k-1,:])-g; - g=g+pi[k-1,l]; - end - end - - P=d+1; - for k=2:n - L = P; - for j=1:d, m=L-pi[k,j]+1:L - P=P+1; - inds[P,:]=inds[m,:]; - inds[P,j]=inds[P,j]+1; - end - end - - return inds -end - -function calculateMultiIndices(a::Array{T,1}) where T<:Integer - #a contains maximum degrees of univariate bases - #function to calculate indices of multivariate basis following the algorithm - #from "Spectral Methods for Uncertainty Quantification; Le Maitre, Knio;2014" p.516-517 - d = length(a) - n = maximum(a) - No = numberPolynomials(d,n) #number of polynomials of multivariate basis - if d==0 - print("Input zero in indicesMulti!") - return zeros(Int64, 1,d) - end - - # catch case n == 0 --> No-d==0 - if n==0 - inds=zeros(Int64,1,d) - return inds - end - - inds = vcat(zeros(Int64,1,d),Matrix(1I,d,d),zeros(Int64, No-d-1, d)); #initiate index matrix for basis - pi=ones(Int64,No,d); - - for k=2:No - g = 0; - for l=1:d - pi[k,l]=sum(pi[k-1,:])-g; - g=g+pi[k-1,l]; - end - end - - P=d+1; - for k=2:n - L = P; - for j=1:d, m=L-pi[k,j]+1:L - P=P+1; - inds[P,:]=inds[m,:]; - inds[P,j]=inds[P,j]+1; - end - end - - #remove elements of high order - #TODO more efficient way to do that? - for j=1:length(a) - k=size(inds,1) - l=1 - while l<=k - if inds[l,j]>a[j] - inds=vcat(inds[1:l-1,:],inds[l+1:end,:]) - k-=1 - l-=1 - end - l+=1 - end - end - - return inds -end - - -""" - computes the number of polynomials with a multivariate basis - `(d+n)!/(d!+n!)` -""" -function numberPolynomials(a::Array{T,1}) where T<:Integer - d = length(a) - n = maximum(a) - - x = max(d,n) - y = min(d,n) - - return (prod(UInt128(x+1):UInt128(d+n))/factorial(UInt128(y))) -end - -function numberPolynomials(d::Int64, n::Int64) - x = max(d,n) - y = min(d,n) - return UInt128(prod(UInt128(x+1):UInt128(d+n))/factorial(UInt128(y))) -end - -""" - findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64} -Given the multi-index `ind` this function returns all entries of the multivariate basis -that correspond to the `i`th univariate basis. -""" -function findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64} - l,p = size(ind) - @assert i<=p "basis is $p-variate, you requested $i-variate" - deg = ind[end,end] - @assert deg>=0 "invalid degree" - col = ind[:,i] - myind = zeros(Int64,deg) - for deg_ = 1:deg - myind[deg_] = findfirst(x->x==deg_,col) - end - pushfirst!(myind,1) -end - -################################################################# -# function calculateMultiIndices_interaction(nξ::Int,deg::Int,j::Int,p::Int) -# inds = calculateMultiIndices(nξ,deg) -# get_interaction(inds,j,p) -# end - -# function get_interaction(inds::Matrix{Int},j::Int,p::Int) -# j < 0 && throw(error("interaction order must be non-negative")) -# j > size(inds,2) && throw(error("interaction order cannot be greater than number of uncertainties")) -# Iterators.filter(x -> count(!iszero,x) == j && sum(x) == p, eachrow(inds)) -# # alternative -> collect -# end diff --git a/src/scalarproduct.jl b/src/scalarproduct.jl deleted file mode 100644 index 235d832a..00000000 --- a/src/scalarproduct.jl +++ /dev/null @@ -1,149 +0,0 @@ -export computeSP, - computeSP2 - -function computeSP(a_::Vector{<:Int}, - α::Vector{<:Vector{<:Real}},β::Vector{<:Vector{<:Real}}, - nodes::Vector{<:Vector{<:Real}},weights::Vector{<:Vector{<:Real}}, - ind::Matrix{<:Int}; - issymmetric::BitArray=falses(length(α)), - zerotol::Float64=1e-10) - minimum(a_) < 0 && throw(DomainError(minimum(a_),"no negative degrees allowed")) - l, p = size(ind) # p-variate basis - p == 1 && computeSP(a_,α[1],β[1],nodes[1],weights[1];issymmetric=issymmetric[1]) - l -= 1 - maximum(a_) > l && throw(DomainError(maximum(a_), "not enough elements in multi-index (requested: $(maximum(a)), max: $l)")) - !(length(α)==length(β)==length(nodes)==length(weights)==length(issymmetric)==p) && throm(InconsistencyError("inconsistent number of rec. coefficients and/or nodes/weights")) - - a = Vector{Int64}() - - for aa in a_ - !iszero(aa) && push!(a,aa) - end - - if iszero(length(a)) - return prod(β[i][1] for i in 1:p) - elseif isone(length(a)) - return 0. - else - inds_uni = multi2uni(a,ind) - val = 1. - for i in 1:p - # compute univariate scalar product - @inbounds v = computeSP(inds_uni[i,:],α[i],β[i],nodes[i],weights[i];issymmetric=issymmetric[i]) - isapprox(v,0,atol=zerotol) ? (return 0.) : (val*=v) - end - end - return val -end - -function computeSP(a::Vector{<:Int}, op::Vector{<:AbstractOrthoPoly}, ind::Matrix{<:Int}) - α, β = coeffs(op) - nodes, weights = nw(op) - computeSP(a, α, β, nodes, weights, ind; issymmetric=issymmetric.(op)) -end - -computeSP(a::Vector{<:Int},mop::MultiOrthoPoly) = computeSP(a, mop.uni, mop.ind) - -""" -__Univariate__ -``` -computeSP(a::Vector{<:Int},α::Vector{<:Real},β::Vector{<:Real},nodes::Vector{<:Real},weights::Vector{<:Real};issymmetric::Bool=false) -computeSP(a::Vector{<:Int},op::AbstractOrthoPoly;issymmetric=issymmetric(op)) -``` -__Multivariate__ -``` -computeSP( a::Vector{<:Int}, - α::Vector{<:Vector{<:Real}},β::Vector{<:Vector{<:Real}}, - nodes::Vector{<:Vector{<:Real}},weights::Vector{<:Vector{<:Real}}, - ind::Matrix{<:Int}; - issymmetric::BitArray=falses(length(α))) -computeSP(a::Vector{<:Int},op::Vector{<:AbstractOrthoPoly},ind::Matrix{<:Int}) -computeSP(a::Vector{<:Int},mOP::MultiOrthoPoly) -``` - -Computes the scalar product -```math -\\langle \\phi_{a_1},\\phi_{a_2},\\cdots,\\phi_{a_n} \\rangle, -``` -where `n = length(a)`. -This requires to provide the recurrence coefficients `(α,β)` and the quadrature rule -`(nodes,weights)`, as well as the multi-index `ind`. -If provided via the keyword `issymmetric`, symmetry of the weight function is exploited. -All computations of the multivariate scalar products resort back to efficient computations -of the univariate scalar products. -Mathematically, this follows from Fubini's theorem. - -The function is dispatched to facilitate its use with `AbstractOrthoPoly` and its quadrature rule `Quad`. - -!!! note - - Zero entries of ``a`` are removed automatically to simplify computations, which follows from - ```math - \\langle \\phi_i, \\phi_j, \\phi_0,\\cdots,\\phi_0 \\rangle = \\langle \\phi_i, \\phi_j \\rangle, - ``` - because `\\phi_0 = 1`. - - It is checked whether enough quadrature points are supplied to solve the integral exactly. -""" -function computeSP(a_::Vector{<:Int},α::Vector{<:Real},β::Vector{<:Real},nodes::Vector{<:Real},weights::Vector{<:Real};issymmetric::Bool=false) - minimum(a_) < 0 && throw(DomainError(minimum(a_), "no negative degrees allowed")) - # this works because ``<Φ_i,Φ_j,Φ_0,...,Φ_0> = <Φ_i,Φ_j>`` - # hence, avoiding unnecessary operations that introduce numerical gibberish - # a = Vector{Int64}() - # for aa in a_ - # !iszero(aa) && push!(a,aa) - # end - a = filter(x -> x > 0, a_) - # simplification in case the measure is symmetric w.r.t t=0 - # exploit symmetry of the density, see Theorem 1.17 in W. Gautschi's book - # and exploit the fact that - (issymmetric && isodd(sum(a))) && return 0. - # Gauss quadrature rules have exactness 2N-1 - !(Int(ceil(0.5*(sum(a)+1)))<=length(nodes)) && throw(InconsistencyError("not enough nodes to integrate exactly ($(length(nodes)) provided, where $(Int(ceil(0.5*(sum(a)+1)))) are needed)")) - - res = - if iszero(length(a)) - first(β) - elseif isone(length(a)) - 0. - elseif length(a) == 2 - @inbounds a[1] == a[2] ? computeSP2(first(a),β)[end] : 0. - else - f = ones(Float64,length(nodes)) - @simd for i = 1:length(a) - @inbounds f = f.*evaluate(a[i],nodes,α,β) - end - dot(weights,f) - end -end - -function computeSP(a::Vector{<:Int},op::AbstractOrthoPoly) - computeSP(a,op.α,op.β,op.quad.nodes,op.quad.weights;issymmetric=issymmetric(op)) -end - -""" - computeSP2(n::Int,β::Vector{<:Real}) - computeSP2(n::Int,op::AbstractOrthoPoly) = computeSP2(n,op.β) - computeSP2(op::AbstractOrthoPoly) = computeSP2(op.deg,op.β) - -Computes the `n` *regular* scalar products aka 2-norms of the orthogonal polynomials, namely -```math -\\|ϕ_i\\|^2 = \\langle \\phi_i,\\phi_i\\rangle \\quad \\forall i \\in \\{ 0,\\dots,n \\}. -``` -Notice that only the values of `β` of the recurrence coefficients `(α,β)` are required. -The computation is based on equation (1.3.7) from Gautschi, W. "Orthogonal Polynomials: Computation and Approximation". -Whenever there exists an analytic expressions for `β`, this function should be used. - -The function is multiply dispatched to facilitate its use with `AbstractOrthoPoly`. -""" -function computeSP2(n::Int,β::Vector{<:Real}) - n < 0 && throw(DomainError(n, "can only compute scalar products for non-negative degrees")) - length(β) < n + 1 && throw(InconsistencyError("inconsistent length of β")) - n == 0 && return β[1] - s = ones(Float64,n+1) - @inbounds s[1] = β[1] # order 0computeSP2 - for i in 2:n+1 - @inbounds s[i] = s[i-1]*β[i] - end - return s -end -computeSP2(n::Int,op::AbstractOrthoPoly) = computeSP2(n,op.β) -computeSP2(op::AbstractOrthoPoly) = computeSP2(op.deg,op.β) \ No newline at end of file From 779573ce7c17b6c1ddf5563673a74a9206e178a3 Mon Sep 17 00:00:00 2001 From: timueh Date: Thu, 12 Sep 2019 19:34:19 +0200 Subject: [PATCH 8/8] updated dependencies + NEWS --- Manifest.toml | 34 +++++++++++++++++++--------------- NEWS.md | 5 +++++ Project.toml | 2 +- todos.md | 16 ---------------- 4 files changed, 25 insertions(+), 32 deletions(-) delete mode 100644 todos.md diff --git a/Manifest.toml b/Manifest.toml index 189abedd..be002f01 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -57,6 +57,11 @@ git-tree-sha1 = "9a11d428dcdc425072af4aea19ab1e8c3e01c032" uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" version = "1.3.0" +[[DataAPI]] +git-tree-sha1 = "8903f0219d3472543fc4b2f5ebaf675a07f817c0" +uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +version = "1.0.1" + [[DataStructures]] deps = ["InteractiveUtils", "OrderedCollections", "Random", "Serialization", "Test"] git-tree-sha1 = "ca971f03e146cf144a9e2f2ce59674f5bf0e8038" @@ -89,15 +94,15 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[Distributions]] deps = ["LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"] -git-tree-sha1 = "4ae0ff1a3a556383f7b1b004f5fb964b26ac96c9" +git-tree-sha1 = "baaf9e165ba8a2d11fb4fb3511782ee070ee3694" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.20.0" +version = "0.21.1" [[FFTW]] -deps = ["AbstractFFTs", "BinaryProvider", "Compat", "Conda", "Libdl", "LinearAlgebra", "Reexport", "Test"] -git-tree-sha1 = "29cda58afbf62f35b1a094882ad6c745a47b2eaa" +deps = ["AbstractFFTs", "BinaryProvider", "Conda", "Libdl", "LinearAlgebra", "Reexport", "Test"] +git-tree-sha1 = "03f8776fbdae28c20c0d1d2ae4e090cd1dfcd247" uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "0.2.4" +version = "1.0.0" [[ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"] @@ -139,10 +144,9 @@ deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" [[Missings]] -deps = ["SparseArrays", "Test"] -git-tree-sha1 = "f0719736664b4358aa9ec173077d4285775f8007" +git-tree-sha1 = "29858ce6c8ae629cf2d733bffa329619a1c843d0" uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "0.4.1" +version = "0.4.2" [[Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" @@ -167,9 +171,9 @@ version = "0.9.9" [[Parsers]] deps = ["Dates", "Test"] -git-tree-sha1 = "db2b35dedab3c0e46dc15996d170af07a5ab91c9" +git-tree-sha1 = "ef0af6c8601db18c282d092ccbd2f01f3f0cd70b" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "0.3.6" +version = "0.3.7" [[Pkg]] deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] @@ -187,9 +191,9 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" [[QuadGK]] deps = ["DataStructures", "LinearAlgebra"] -git-tree-sha1 = "389fd27b958a33df5234772cc464b663b208273d" +git-tree-sha1 = "438630b843c210b375b2a246329200c113acc61b" uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -version = "2.0.4" +version = "2.1.0" [[REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets"] @@ -251,10 +255,10 @@ deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [[StatsBase]] -deps = ["DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"] -git-tree-sha1 = "8a0f4b09c7426478ab677245ab2b0b68552143c7" +deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"] +git-tree-sha1 = "c53e809e63fe5cf5de13632090bc3520649c9950" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.30.0" +version = "0.32.0" [[StatsFuns]] deps = ["Rmath", "SpecialFunctions", "Test"] diff --git a/NEWS.md b/NEWS.md index 966f6f90..ba71c9a8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,10 @@ # PolyChaos Release Notes +## Version 0.2.1 +- increased code coverage +- examples work with new type hierarchy +- updated dependencies + ## Version 0.2.0 - new type hierarchy that involves four abstract types: - `AbstractMeasure` diff --git a/Project.toml b/Project.toml index 160a9c63..69653458 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PolyChaos" uuid = "8d666b04-775d-5f6e-b778-5ac7c70f65a3" authors = ["timueh "] -version = "0.2.0" +version = "0.2.1" [deps] AdaptiveRejectionSampling = "c75e803d-635f-53bd-ab7d-544e482d8c75" diff --git a/todos.md b/todos.md deleted file mode 100644 index 597f7a69..00000000 --- a/todos.md +++ /dev/null @@ -1,16 +0,0 @@ -# Todos for PolyChaos.jl - -## High priority - - documentation: - - make conversion to Latex work - - examples: - - quadrature rules (fejer, fejer2, clenshaw-curtis, gauss, radau, lobatto) - - ~~discretization procedures (stieltjes, lanczos), see `test/discretization.jl`~~ - - `quadgp()` - - optimizing Rosenbrock/quadratic function - - DC-OPF with stochastic uncertainties - - -## Low priority - - orthonormal polynomials - - arbitrary polynomials