Merge branch 'main' into dv/docdef

Author: dlfivefifty

.github/workflows/downstream.yml

@@ -0,0 +1,71 @@
+name: IntegrationTest
+on:
+  push:
+    branches: [main]
+    tags: [v*]
+  pull_request:
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
+
+jobs:
+  pre_job:
+    # continue-on-error: true # Uncomment once integration is finished
+    runs-on: ubuntu-latest
+    # Map a step output to a job output
+    outputs:
+      should_skip: ${{ steps.skip_check.outputs.should_skip }}
+    steps:
+      - id: skip_check
+        uses: fkirc/skip-duplicate-actions@v5
+  test:
+    needs: pre_job
+    if: needs.pre_job.outputs.should_skip != 'true'
+    name: ${{ matrix.package.group }}/${{ matrix.package.repo }}/${{ matrix.julia-version }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        julia-version: ['1']
+        os: [ubuntu-latest]
+        package:
+          - {repo: HarmonicOrthogonalPolynomials.jl, group: JuliaApproximation}
+          - {repo: MultivariateOrthogonalPolynomials.jl, group: JuliaApproximation}
+          - {repo: SemiclassicalOrthogonalPolynomials.jl, group: JuliaApproximation}
+          - {repo: PiecewiseOrthogonalPolynomials.jl, group: JuliaApproximation}
+          - {repo: SingularIntegrals.jl, group: JuliaApproximation}
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.julia-version }}
+          arch: x64
+      - uses: julia-actions/julia-buildpkg@latest
+      - name: Clone Downstream
+        uses: actions/checkout@v4
+        with:
+          repository: ${{ matrix.package.group }}/${{ matrix.package.repo }}
+          path: downstream
+      - name: Load this and run the downstream tests
+        shell: julia --color=yes --project=downstream {0}
+        run: |
+          using Pkg
+          try
+            # force it to use this PR's version of the package
+            Pkg.develop(PackageSpec(path="."))  # resolver may fail with main deps
+            Pkg.update()
+            Pkg.test(; coverage = true)  # resolver may fail with test time deps
+          catch err
+            err isa Pkg.Resolve.ResolverError || rethrow()
+            # If we can't resolve that means this is incompatible by SemVer and this is fine
+            # It means we marked this as a breaking change, so we don't need to worry about
+            # Mistakenly introducing a breaking change, as we have intentionally made one
+            @info "Not compatible with this release. No problem." exception=err
+            exit(0)  # Exit immediately, as a success
+          end
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v5
+        with:
+          token: ${{ secrets.CODECOV_TOKEN }}
+          files: lcov.info

Project.toml

@@ -1,7 +1,8 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.14.1"
+version = "0.14.4"
+
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -42,7 +43,7 @@ DomainSets = "0.6, 0.7"
 DynamicPolynomials = "0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "1"
-FastTransforms = "0.16.6"
+FastTransforms = "0.16.6, 0.17"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = " 0.14, 0.15"
@@ -53,7 +54,7 @@ LazyBandedMatrices = "0.10, 0.11"
 MutableArithmetics = "1"
 QuasiArrays = "0.11"
 RecurrenceRelationshipArrays = "0.1.2"
-RecurrenceRelationships = "0.1.1"
+RecurrenceRelationships = "0.1.1, 0.2"
 SpecialFunctions = "1.0, 2"
 julia = "1.10"
 

docs/src/index.md

@@ -17,6 +17,17 @@ which defines the following classical orthogonal polynomials:
 5. Laguerre: $L_n^{(\alpha)}(x)$, defined over $[0, ∞)$ with weight $w(x) = x^\alpha \mathrm{e}^{-x}$.
 6. Hermite: $H_n(x)$, defined over $(-∞, ∞)$ with weight $w(x) = \mathrm{e}^{-x^2}$.
 
+These special polynomials have many applications and can be used as a basis for any function given their domain conditions are met, however these polynomials have some advantages due to their formulation:
+
+- Because of their relation to Laplace’s equation, **Legendre polynomials** can be useful as a basis for functions with spherical symmetry.
+- **Chebyshev polynomials** are generally effective in reducing errors from numerical methods such as quadrature, interpolation, and approximation.
+- Due to the flexibility of its parameters, **Jacobi polynomials** are capable of tailoring the behavior of an approximation around its endpoints, making these polynomials particularly useful in boundary value problems.
+- **Ultraspherical polynomials** are advantageous in spectral methods for solving differential equations.
+- **Laguerre polynomials** have a semi-infinite domain, therefore they are beneficial for problems involving exponential decay.
+- Because of its weight function, **Hermite polynomials** can be useful in situations where functions display a Gaussian-like distribution.
+
+These are just a few applications of these polynomials. They have many more uses across mathematics, physics, and engineering.
+
 ## Evaluation
 
 The simplest usage of this package is to evaluate classical

src/choleskyQR.jl

@@ -1,5 +1,5 @@
 """
-Represent an Orthogonal polynomial which has a conversion operator from P, that is, Q = P * inv(U).
+Represents orthonormal polynomials defined via a conversion operator from P, that is, Q = P * inv(U).
 """
 struct ConvertedOrthogonalPolynomial{T, WW<:AbstractQuasiVector{T}, XX, UU, PP} <: OrthonormalPolynomial{T}
     weight::WW
@@ -8,15 +8,22 @@ struct ConvertedOrthogonalPolynomial{T, WW<:AbstractQuasiVector{T}, XX, UU, PP}
     P::PP
 end
 
-_p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P)
+_p0(Q::ConvertedOrthogonalPolynomial) = _p0(Q.P)/Q.U[1,1]
 
 axes(Q::ConvertedOrthogonalPolynomial) = axes(Q.P)
+
+
+struct ConvertedOPLayout <: AbstractNormalizedOPLayout end
 MemoryLayout(::Type{<:ConvertedOrthogonalPolynomial}) = ConvertedOPLayout()
+
+
+
+
 jacobimatrix(Q::ConvertedOrthogonalPolynomial) = Q.X
 orthogonalityweight(Q::ConvertedOrthogonalPolynomial) = Q.weight
 
 
-# transform to P * U if needed for differentiation, etc.
+# transform to P * inv(U) if needed for differentiation, etc.
 arguments(::ApplyLayout{typeof(*)}, Q::ConvertedOrthogonalPolynomial) = Q.P, ApplyArray(inv, Q.U)
 
 OrthogonalPolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w, orthogonalpolynomial(singularities(w)))

src/classical/jacobi.jl

@@ -217,8 +217,8 @@ basis_singularities(w::JacobiWeight) = Weighted(Jacobi(w.a, w.b))
 computes the `n`-th Jacobi polynomial, orthogonal with
 respec to `(1-x)^a*(1+x)^b`, at `z`.
 """
-jacobip(n::Integer, a, b, z) = Base.unsafe_getindex(Jacobi{polynomialtype(typeof(a), typeof(b), typeof(z))}(a,b), z, n+1)
-normalizedjacobip(n::Integer, a, b, z) = Base.unsafe_getindex(Normalized(Jacobi{polynomialtype(typeof(a), typeof(b), typeof(z))}(a,b)), z, n+1)
+jacobip(n::Integer, a, b, z) = Base.unsafe_getindex(Jacobi{polynomialtype(promote_type(typeof(a), typeof(b)), typeof(z))}(a,b), z, n+1)
+normalizedjacobip(n::Integer, a, b, z) = Base.unsafe_getindex(Normalized(Jacobi{polynomialtype(promote_type(typeof(a), typeof(b)), typeof(z))}(a,b)), z, n+1)
 
 OrthogonalPolynomial(w::JacobiWeight) = Jacobi(w.a, w.b)
 orthogonalityweight(P::Jacobi) = JacobiWeight(P.a, P.b)
@@ -294,7 +294,7 @@ grid(P::AbstractJacobi{T}, n::Integer) where T = ChebyshevGrid{1,T}(n)
 plotgrid(P::AbstractJacobi{T}, n::Integer) where T = ChebyshevGrid{2,T}(min(40n, MAX_PLOT_POINTS))
 
 plan_transform(::AbstractJacobi{T}, szs::NTuple{N,Int}, dims...) where {T,N} = error("Override")
-plan_transform(P::Jacobi{T}, szs::NTuple{N,Int}, dims...) where {T,N} = JacobiTransformPlan(FastTransforms.plan_th_cheb2jac!(T, szs, P.a, P.b, dims), plan_chebyshevtransform(T, szs, dims...))
+plan_transform(P::Jacobi{T}, szs::NTuple{N,Int}, dims...) where {T,N} = JacobiTransformPlan(FastTransforms.plan_th_cheb2jac!(T, szs, P.a, P.b, dims...), plan_chebyshevtransform(T, szs, dims...))
 
 ldiv(P::Jacobi{V}, f::Inclusion{T}) where {T,V} = _op_ldiv(P, f)
 ldiv(P::Jacobi{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)

src/classical/ultraspherical.jl

@@ -73,7 +73,7 @@ ultraspherical(λ, d::AbstractInterval{T}) where T = Ultraspherical{float(promot
 # transforms
 ###
 
-plan_transform(P::Ultraspherical{T}, szs::NTuple{N,Int}, dims...) where {T,N} = JacobiTransformPlan(FastTransforms.plan_th_ultra2ultra!(T, szs, one(P.λ), P.λ, dims), plan_chebyshevutransform(T, szs, dims...))
+plan_transform(P::Ultraspherical{T}, szs::NTuple{N,Int}, dims...) where {T,N} = JacobiTransformPlan(FastTransforms.plan_th_ultra2ultra!(T, szs, one(P.λ), P.λ, dims...), plan_chebyshevutransform(T, szs, dims...))
 
 ###
 # interrelationships

src/normalized.jl

@@ -35,7 +35,6 @@ isnormalized(_) = false
 represents OPs that are of the form P * R where P is another family of OPs and R is upper-triangular.
 """
 abstract type AbstractNormalizedOPLayout <: AbstractOPLayout end
-struct ConvertedOPLayout <: AbstractNormalizedOPLayout end
 struct NormalizedOPLayout{LAY<:AbstractBasisLayout} <: AbstractNormalizedOPLayout end
 
 MemoryLayout(::Type{<:Normalized{<:Any, OPs}}) where OPs = NormalizedOPLayout{typeof(MemoryLayout(OPs))}()
@@ -302,4 +301,4 @@ basismap(W::Weighted) = basismap(W.P)
 const MappedOPLayouts = Union{MappedOPLayout,WeightedOPLayout{MappedOPLayout}}
 diff_layout(::MappedOPLayouts, A, dims...) = diff_layout(MappedBasisLayout(), A, dims...)
 
-diff_layout(::NormalizedOPLayout, A, dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, dims...)
+diff_layout(::AbstractNormalizedOPLayout, A, dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, dims...)

test/test_choleskyQR.jl

@@ -1,5 +1,5 @@
 using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices, InfiniteLinearAlgebra
-import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial
+import ClassicalOrthogonalPolynomials: cholesky_jacobimatrix, qr_jacobimatrix, orthogonalpolynomial, _p0
 import LazyArrays: AbstractCachedMatrix, resizedata!
 
 @testset "CholeskyQR" begin
@@ -76,7 +76,7 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
                 # Comparison with Lanczos
                 @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
             end
-            
+
             @testset "w(x) = (1-x)*exp(x)" begin
                 P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
                 x = axes(P,1)
@@ -89,7 +89,7 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
                 # Comparison with Lanczos
                 @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
             end
-            
+
             @testset "w(x) = (1-x^2)*exp(x^2)" begin
                 P = Normalized(legendre(0..1))
                 x = axes(P,1)
@@ -102,7 +102,7 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
                 # Comparison with Lanczos
                 @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
             end
-            
+
             @testset "w(x) = x*(1-x^2)*exp(-x^2)" begin
                 P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
                 x = axes(P,1)
@@ -227,10 +227,9 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
         @test Q ≠ P
         @test Q == Q̃
         @test Q̃ == Q
-        
-        @test Q[0.1,1] ≈ 1/sqrt(2)
+
+        @test Q[0.1,1] ≈ _p0(Q) ≈ 1/sqrt(2)
         @test Q[0.1,1:10] ≈ Q̃[0.1,1:10]
-        # AWESOME, thanks TSGUT!!
         @test Q[0.1,10_000] ≈ Q̃[0.1,10_000]
 
         R = P \ Q
@@ -247,9 +246,21 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
     @testset "Chebyshev" begin
          U = ChebyshevU()
          Q = orthogonalpolynomial(x -> (1+x^2)*sqrt(1-x^2), U)
+         x = axes(U,1)
+
+         @test Q[0.1,1] ≈ _p0(Q) ≈ 1/sqrt(sum(expand(Weighted(U),x -> (1+x^2)*sqrt(1-x^2))))
          @test bandwidths(Q\U) == (0,2)
 
          Q̃ = OrthogonalPolynomial(x -> (1+x^2)*sqrt(1-x^2), U)
          @test jacobimatrix(Q)[1:10,1:10] == jacobimatrix(Q̃)[1:10,1:10]
     end
+
+    @testset "diff" begin
+        P = Legendre()
+        x = axes(P,1)
+        Q = OrthogonalPolynomial(1 .- x)
+        Q̃ = Normalized(Jacobi(1,0))
+        @test_skip diff(Q)[0.1,1:5] ≈ diff(Q̃)[0.1,1:5] # broken due to lazy array issues
+        @test [diff(Q)[0.1,k] for k=1:5] ≈ diff(Q̃)[0.1,1:5]
+    end
 end

test/test_gram.jl

@@ -0,0 +1,10 @@
+using ClassicalOrthogonalPolynomials, FastTransforms
+
+P = Legendre()
+x = axes(P,1)
+w = @.(1-x^2)
+
+μ = P'w
+X = jacobimatrix(P)
+n = 20
+@test GramMatrix(μ[1:2n], X[1:2n,1:2n]) ≈ (P' * (w .* P))[1:n,1:n]