|
1 | 1 | module MultivariateOrthogonalPolynomials |
2 | | -using Base, RecipesBase, ApproxFun, BandedMatrices, BlockArrays, BlockBandedMatrices, |
3 | | - FastTransforms, FastGaussQuadrature, StaticArrays, FillArrays, |
4 | | - LinearAlgebra, Libdl, SpecialFunctions, LazyArrays, InfiniteArrays, |
5 | | - DomainSets, ArrayLayouts |
| 2 | +using OrthogonalPolynomialsQuasi, FastTransforms, BlockBandedMatrices, BlockArrays, DomainSets, |
| 3 | + QuasiArrays, StaticArrays, ContinuumArrays, InfiniteArrays, InfiniteLinearAlgebra, |
| 4 | + LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts |
6 | 5 |
|
7 | | -# package code goes here |
8 | | -import Base: values,getindex,setindex!,*, +, -, ==,<,<=,>, |
9 | | - >=,/,^,\,∪,transpose, in, convert, issubset |
| 6 | +import Base: axes, in, ==, *, ^, \, copy, OneTo, getindex, size |
| 7 | +import DomainSets: boundary |
10 | 8 |
|
| 9 | +import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle |
| 10 | +import ContinuumArrays: @simplify, Weight, grid, TransformFactorization, Expansion |
11 | 11 |
|
12 | | -import BandedMatrices: inbands_getindex, inbands_setindex! |
| 12 | +import BlockArrays: block, blockindex, BlockSlice, viewblock |
| 13 | +import BlockBandedMatrices: _BandedBlockBandedMatrix |
| 14 | +import LinearAlgebra: factorize |
| 15 | +import LazyArrays: arguments, paddeddata |
13 | 16 |
|
14 | | -import BlockArrays: getblock, Block, blocksize |
| 17 | +import OrthogonalPolynomialsQuasi: jacobimatrix |
15 | 18 |
|
16 | | -import BlockBandedMatrices: blockbandwidths, subblockbandwidths |
| 19 | +export Triangle, JacobiTriangle, TriangleWeight, WeightedTriangle, PartialDerivative, Laplacian, MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial |
17 | 20 |
|
18 | | -# ApproxFun general import |
19 | | -import ApproxFunBase: BandedMatrix, |
20 | | - linesum,complexlength, BandedBlockBandedMatrix, |
21 | | - real, eps, isapproxinteger, FiniteRange, DFunction, |
22 | | - TransformPlan, ITransformPlan, plan_transform! |
| 21 | +######### |
| 22 | +# PartialDerivative{k} |
| 23 | +# takes a partial derivative in the k-th index. |
| 24 | +######### |
| 25 | + |
| 26 | + |
| 27 | +struct PartialDerivative{k,T,Ax<:Inclusion} <: LazyQuasiMatrix{T} |
| 28 | + axis::Ax |
| 29 | +end |
| 30 | + |
| 31 | +PartialDerivative{k,T}(axis::Inclusion) where {k,T} = PartialDerivative{k,T,typeof(axis)}(axis) |
| 32 | +PartialDerivative{k,T}(domain) where {k,T} = PartialDerivative{k,T}(Inclusion(domain)) |
| 33 | +PartialDerivative{k}(axis) where k = PartialDerivative{k,eltype(eltype(axis))}(axis) |
23 | 34 |
|
24 | | -# Domains import |
25 | | -import ApproxFunBase: fromcanonical, tocanonical, domainscompatible |
| 35 | +axes(D::PartialDerivative) = (D.axis, D.axis) |
| 36 | +==(a::PartialDerivative{k}, b::PartialDerivative{k}) where k = a.axis == b.axis |
| 37 | +copy(D::PartialDerivative{k}) where k = PartialDerivative{k}(copy(D.axis)) |
26 | 38 |
|
27 | | -# Operator import |
28 | | -import ApproxFunBase: bandwidths,SpaceOperator, ConversionWrapper, DerivativeWrapper, |
29 | | - rangespace, domainspace, InterlaceOperator, |
30 | | - promotedomainspace, CalculusOperator, interlace, Multiplication, |
31 | | - choosedomainspace, SubOperator, ZeroOperator, |
32 | | - Dirichlet, DirichletWrapper, Neumann, Laplacian, ConstantTimesOperator, Conversion, |
33 | | - Derivative, ConcreteMultiplication, ConcreteConversion, ConcreteLaplacian, |
34 | | - ConcreteDerivative, TimesOperator, MultiplicationWrapper, TridiagonalOperator |
| 39 | +struct Laplacian{T,D} <: LazyQuasiMatrix{T} |
| 40 | + axis::Inclusion{T,D} |
| 41 | +end |
35 | 42 |
|
| 43 | +Laplacian{T}(axis::Inclusion{<:Any,D}) where {T,D} = Laplacian{T,D}(axis) |
| 44 | +Laplacian{T}(domain) where T = Laplacian{T}(Inclusion(domain)) |
| 45 | +Laplacian(axis) = Laplacian{eltype(axis)}(axis) |
36 | 46 |
|
37 | | -# Spaces import |
38 | | -import ApproxFunBase: ConstantSpace, NoSpace, prectype, |
39 | | - SumSpace,PiecewiseSpace, ArraySpace, @containsconstants, |
40 | | - UnsetSpace, canonicalspace, canonicaldomain, domain, evaluate, |
41 | | - AnyDomain, plan_transform,plan_itransform, |
42 | | - transform,itransform,transform!,itransform!, |
43 | | - isambiguous, fromcanonical, tocanonical, checkpoints, ∂, spacescompatible, |
44 | | - mappoint, UnivariateSpace, setdomain, setcanonicaldomain, canonicaldomain, |
45 | | - Space, points, space, conversion_rule, maxspace_rule, |
46 | | - union_rule, coefficients, RealUnivariateSpace, PiecewiseSegment, rangetype, cfstype |
| 47 | +axes(D::Laplacian) = (D.axis, D.axis) |
| 48 | +==(a::Laplacian, b::Laplacian) = a.axis == b.axis |
| 49 | +copy(D::Laplacian) = Laplacian(copy(D.axis), D.k) |
47 | 50 |
|
48 | | -# Multivariate import |
49 | | -import ApproxFunBase: DirectSumSpace, AbstractProductSpace, factor, |
50 | | - BivariateFun, ProductFun, LowRankFun, lap, columnspace, |
51 | | - fromtensor, totensor, isbandedblockbanded, |
52 | | - Tensorizer, tensorizer, block, blockstart, blockstop, blocklengths, |
53 | | - domaintensorizer, rangetensorizer, blockrange, BlockRange1, |
54 | | - float |
| 51 | +^(D::PartialDerivative, k::Integer) = ApplyQuasiArray(^, D, k) |
| 52 | +^(D::Laplacian, k::Integer) = ApplyQuasiArray(^, D, k) |
55 | 53 |
|
56 | | -# Singularities |
57 | | -import ApproxFunBase: WeightSpace, weight |
58 | 54 |
|
59 | | -# Vec is for two points |
60 | | -import ApproxFunBase: Vec |
| 55 | +abstract type MultivariateOrthogonalPolynomial{d,T} <: Basis{T} end |
| 56 | +const BivariateOrthogonalPolynomial{T} = MultivariateOrthogonalPolynomial{2,T} |
| 57 | +const BlockOneTo = BlockRange{1,Tuple{OneTo{Int}}} |
61 | 58 |
|
62 | | -# Jacobi import |
63 | | -import ApproxFunOrthogonalPolynomials: jacobip, JacobiSD, PolynomialSpace, order |
64 | 59 |
|
65 | | -import ApproxFunFourier: polar, ipolar |
| 60 | +getindex(P::MultivariateOrthogonalPolynomial{<:Any,D}, xy::SVector{D}, JR::BlockOneTo) where D = |
| 61 | + error("Overload") |
| 62 | +getindex(P::MultivariateOrthogonalPolynomial{D}, xy::SVector{D}, J::Block{1}) where D = P[xy, Block.(OneTo(Int(J)))][J] |
| 63 | +getindex(P::MultivariateOrthogonalPolynomial{D}, xy::SVector{D}, JR::BlockRange{1}) where D = P[xy, Block.(OneTo(Int(maximum(JR))))][JR] |
| 64 | +getindex(P::MultivariateOrthogonalPolynomial{D}, xy::SVector{D}, Jj::BlockIndex{1}) where D = P[xy, block(Jj)][blockindex(Jj)] |
| 65 | +getindex(P::MultivariateOrthogonalPolynomial{D}, xy::SVector{D}, j::Integer) where D = P[xy, findblockindex(axes(P,2), j)] |
| 66 | +getindex(P::MultivariateOrthogonalPolynomial{D}, xy::SVector{D}, jr::AbstractVector) where D = P[xy, Block.(OneTo(Int(findblock(axes(P,2), maximum(jr)))))][jr] |
66 | 67 |
|
| 68 | +const FirstInclusion = BroadcastQuasiVector{<:Any, typeof(first), <:Tuple{Inclusion}} |
| 69 | +const LastInclusion = BroadcastQuasiVector{<:Any, typeof(last), <:Tuple{Inclusion}} |
67 | 70 |
|
68 | | -include("Triangle/Triangle.jl") |
69 | | -include("Disk/Disk.jl") |
70 | | -include("Cone/Cone.jl") |
| 71 | +function Base.broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::FirstInclusion, P::BivariateOrthogonalPolynomial) |
| 72 | + axes(x,1) == axes(P,1) || throw(DimensionMismatch()) |
| 73 | + P*jacobimatrix(Val(1), P) |
| 74 | +end |
| 75 | + |
| 76 | +function Base.broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::LastInclusion, P::BivariateOrthogonalPolynomial) |
| 77 | + axes(x,1) == axes(P,1) || throw(DimensionMismatch()) |
| 78 | + P*jacobimatrix(Val(2), P) |
| 79 | +end |
| 80 | + |
| 81 | +""" |
| 82 | + forwardrecurrence!(v, A, B, C, (x,y)) |
| 83 | +
|
| 84 | +evaluates the bivaraite orthogonal polynomials at points `(x,y)` , |
| 85 | +where `A`, `B`, and `C` are `AbstractVector`s containing the recurrence coefficients |
| 86 | +matrices. In particular, note that for any OPs we have |
| 87 | +
|
| 88 | + P[N+1,x] = (A[N]* [x*I; y*I] + B[N]) * P[N,x] - C[N] * P[N-1,x] |
| 89 | +
|
| 90 | +where A[N] is (N+1) x 2N, B[N] and C[N] are (N+1) x N. |
| 91 | +
|
| 92 | +""" |
| 93 | +# function forwardrecurrence!(v::AbstractBlockVector{T}, A::AbstractVector, B::AbstractVector, C::AbstractVector, xy) where T |
| 94 | +# N = blocklength(v) |
| 95 | +# N == 0 && return v |
| 96 | +# length(A)+1 ≥ N && length(B)+1 ≥ N && length(C)+1 ≥ N || throw(ArgumentError("A, B, C must contain at least $(N-1) entries")) |
| 97 | +# v[Block(1)] .= one(T) |
| 98 | +# N = 1 && return v |
| 99 | +# p_1 = view(v,Block(2)) |
| 100 | +# p_0 = view(v,Block(1)) |
| 101 | +# mul!(p_1, B[1], p_0) |
| 102 | +# xy_muladd!(xy, A[1], p_0, one(T), p_1) |
71 | 103 |
|
72 | | -include("clenshaw.jl") |
| 104 | +# @inbounds for n = 2:N-1 |
| 105 | +# p_2 = view(v,Block(n+1)) |
| 106 | +# mul!(p_2, B[n], p_1) |
| 107 | +# xy_muladd!(xy, A[n], p_1, one(T), p_2) |
| 108 | +# muladd!(-one(T), C[n], p_0, one(T), p_2) |
| 109 | +# p_1,p_0 = p_2,p_1 |
| 110 | +# end |
| 111 | +# v |
| 112 | +# end |
73 | 113 |
|
74 | | -# include("Sphere/SphericalHarmonics.jl") |
75 | 114 |
|
76 | | -include("show.jl") |
77 | | -include("plot.jl") |
| 115 | +# forwardrecurrence(N::Block{1}, A::AbstractVector, B::AbstractVector, C::AbstractVector, xy) = |
| 116 | +# forwardrecurrence!(PseudoBlockVector{promote_type(eltype(eltype(A)),eltype(eltype(B)),eltype(eltype(C)),eltype(xy))}(undef, 1:Int(N)), A, B, C, xy) |
| 117 | + |
| 118 | +# use block expansion |
| 119 | +ContinuumArrays.transform_ldiv(V::SubQuasiArray{<:Any,2,<:BivariateOrthogonalPolynomial,<:Tuple{Inclusion,BlockSlice{BlockOneTo}}}, B, _) = |
| 120 | + factorize(V) \ B |
| 121 | +function ContinuumArrays.transform_ldiv(V::SubQuasiArray{<:Any,2,<:BivariateOrthogonalPolynomial,<:Tuple{Inclusion,AbstractVector{Int}}}, B, _) |
| 122 | + P = parent(V) |
| 123 | + _,jr = parentindices(V) |
| 124 | + J = findblock(axes(P,2),maximum(jr)) |
| 125 | + (P[:,Block.(OneTo(Int(J)))] \ B)[jr] |
| 126 | +end |
| 127 | + |
| 128 | +# Make sure block structure matches. Probably should do this for all block mul |
| 129 | +LazyArrays.mul(A::MultivariateOrthogonalPolynomial, b::AbstractVector) = |
| 130 | + ApplyQuasiArray(*, A, PseudoBlockVector(b, (axes(A,2),))) |
| 131 | + |
| 132 | + |
| 133 | + |
| 134 | +include("Triangle/Triangle.jl") |
| 135 | + |
78 | 136 |
|
79 | 137 | end # module |
0 commit comments