Merge branch 'main' into dl/basissingularities

Author: dlfivefifty

.github/workflows/ci.yml

@@ -47,7 +47,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info

.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.13.7"
+version = "0.15.1"
+
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -22,40 +23,49 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
-RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
 RecurrenceRelationshipArrays = "b889d2dc-af3c-4820-88a8-238fa91d3518"
+RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
+[weakdeps]
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+
+[extensions]
+ClassicalOrthogonalPolynomialsMutableArithmeticsExt = "MutableArithmetics"
+
 [compat]
 ArrayLayouts = "1.3.1"
 BandedMatrices = "1"
 BlockArrays = "1"
 BlockBandedMatrices = "0.13"
-ContinuumArrays = "0.18.3"
+ContinuumArrays = "0.19"
 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"
-InfiniteLinearAlgebra = "0.8"
+InfiniteArrays = " 0.15"
+InfiniteLinearAlgebra = "0.10"
 IntervalSets = "0.7"
-LazyArrays = "2.2"
-LazyBandedMatrices = "0.10"
-QuasiArrays = "0.11"
-RecurrenceRelationships = "0.1"
-RecurrenceRelationshipArrays = "0.1"
+LazyArrays = "2.5.1"
+LazyBandedMatrices = "0.11"
+MutableArithmetics = "1"
+QuasiArrays = "0.12"
+RecurrenceRelationshipArrays = "0.1.2"
+RecurrenceRelationships = "0.2"
 SpecialFunctions = "1.0, 2"
 julia = "1.10"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "StaticArrays", "Random"]
+test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "StaticArrays", "Random", "DynamicPolynomials"]

docs/make.jl

@@ -1,5 +1,9 @@
 using Documenter, ClassicalOrthogonalPolynomials
 
+DocMeta.setdocmeta!(ClassicalOrthogonalPolynomials,
+    :DocTestSetup,
+    :(using ClassicalOrthogonalPolynomials))
+
 makedocs(
     modules = [ClassicalOrthogonalPolynomials],
     sitename="ClassicalOrthogonalPolynomials.jl",

docs/src/index.md

@@ -10,12 +10,23 @@ CurrentModule = ClassicalOrthogonalPolynomials
 We follow the [Digital Library of Mathematical Functions](https://dlmf.nist.gov/18.3),
 which defines the following classical orthogonal polynomials:
 
-1. Legendre: `P_n(x)`
-2. Chebyshev (1st kind, 2nd kind): `T_n(x)`, `U_n(x)`
-3. Ultraspherical: `C_n^{(λ)}(x)`
-4. Jacobi: `P_n^{(a,b)}(x)`
-5. Laguerre: `L_n^{(α)}(x)`
-6. Hermite: `H_n(x)`
+1. Legendre: $P_n(x)$, defined over $[-1, 1]$ with weight $w(x) = 1$.
+2. Chebyshev (1st kind, 2nd kind): $T_n(x)$ and $U_n(x)$, defined over $[-1, 1]$ with weights $w(x) = 1/\sqrt{1-x^2}$ and $w(x) = \sqrt{1-x^2}$, respectively.
+3. Ultraspherical: $C_n^{(\lambda)}(x)$, defined over $[-1, 1]$ with weight $w(x) = (1-x^2)^{\lambda-1/2}$.
+4. Jacobi: $P_n^{(a,b)}(x)$, defined over $[-1, 1]$ with weight $w(x) = (1-x)^a(1+x)^b$.
+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
 
@@ -126,17 +137,15 @@ U\T
 
 ```@docs
 ClassicalOrthogonalPolynomials.Chebyshev
-```
-```@docs
 ClassicalOrthogonalPolynomials.chebyshevt
-```
-```@docs
 ClassicalOrthogonalPolynomials.chebyshevu
 ```
 ```@docs
+ClassicalOrthogonalPolynomials.Legendre
 ClassicalOrthogonalPolynomials.legendrep
 ```
 ```@docs
+ClassicalOrthogonalPolynomials.Jacobi
 ClassicalOrthogonalPolynomials.jacobip
 ```
 ```@docs
@@ -147,8 +156,6 @@ ClassicalOrthogonalPolynomials.hermiteh
 ```
 
 
-
-
 ### Weights
 
 ```@docs
@@ -161,7 +168,9 @@ ClassicalOrthogonalPolynomials.HermiteWeight
 ClassicalOrthogonalPolynomials.Weighted
 ```
 ```@docs
+ClassicalOrthogonalPolynomials.LegendreWeight
 ClassicalOrthogonalPolynomials.ChebyshevWeight
+ClassicalOrthogonalPolynomials.JacobiWeight
 ```
 ```@docs
 ClassicalOrthogonalPolynomials.LaguerreWeight
@@ -170,6 +179,13 @@ ClassicalOrthogonalPolynomials.LaguerreWeight
 ClassicalOrthogonalPolynomials.HalfWeighted
 ```
 
+### Affine-mapped
+```@docs
+ClassicalOrthogonalPolynomials.legendre
+ClassicalOrthogonalPolynomials.jacobi
+ClassicalOrthogonalPolynomials.legendreweight
+ClassicalOrthogonalPolynomials.jacobiweight
+```
 
 ### Recurrences
 
@@ -190,32 +206,22 @@ ClassicalOrthogonalPolynomials.recurrencecoefficients
 ```
 
 
-
 ### Internal
 
 ```@docs
-ClassicalOrthogonalPolynomials.ShuffledIR2HC
-```
-```@docs
-ClassicalOrthogonalPolynomials.ShuffledR2HC
-```
-```@docs
+ClassicalOrthogonalPolynomials.ShuffledFFT
 ClassicalOrthogonalPolynomials.ShuffledIFFT
+ClassicalOrthogonalPolynomials.ShuffledR2HC
+ClassicalOrthogonalPolynomials.ShuffledIR2HC
 ```
 ```@docs
 ClassicalOrthogonalPolynomials.qr_jacobimatrix
-```
-```@docs
-ClassicalOrthogonalPolynomials.MappedOPLayout
-```
-```@docs
 ClassicalOrthogonalPolynomials.cholesky_jacobimatrix
 ```
 ```@docs
 ClassicalOrthogonalPolynomials.AbstractNormalizedOPLayout
-```
-```@docs
-ClassicalOrthogonalPolynomials.ShuffledFFT
+ClassicalOrthogonalPolynomials.MappedOPLayout
+ClassicalOrthogonalPolynomials.WeightedOPLayout
 ```
 ```@docs
 ClassicalOrthogonalPolynomials.legendre_grammatrix
@@ -224,9 +230,6 @@ ClassicalOrthogonalPolynomials.legendre_grammatrix
 ClassicalOrthogonalPolynomials.weightedgrammatrix
 ```
 ```@docs
-ClassicalOrthogonalPolynomials.WeightedOPLayout
-```
-```@docs
 ClassicalOrthogonalPolynomials.interlace!
 ```
 ```@docs

ext/ClassicalOrthogonalPolynomialsMutableArithmeticsExt.jl

@@ -0,0 +1,8 @@
+module ClassicalOrthogonalPolynomialsMutableArithmeticsExt
+using ClassicalOrthogonalPolynomials, MutableArithmetics
+import ClassicalOrthogonalPolynomials: initiateforwardrecurrence, recurrencecoefficients, _p0
+
+Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractMutable, n::Number) = initiateforwardrecurrence(n-1, recurrencecoefficients(P)..., x, _p0(P))[end]
+
+recurrencecoefficients(::Legendre{T}) where T<:AbstractMutable = recurrencecoefficients(Legendre())
+end

src/ClassicalOrthogonalPolynomials.jl

@@ -29,19 +29,19 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
-                    _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector,
+                    _getindex, layout_getindex, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector,
                     AbstractQuasiFill, equals_layout, QuasiArrayLayout, PolynomialLayout, diff_layout
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
-import InfiniteLinearAlgebra: chop!, chop, pad, choplength, compatible_resize!, partialcholesky!
-import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
+import InfiniteLinearAlgebra: chop!, chop, pad, choplength, compatible_resize!, partialcholesky!, SymTridiagonalConjugation, TridiagonalConjugation
+import ContinuumArrays: Basis, Weight, basis_axes, @simplify, AbstractAffineQuasiVector, ProjectionFactorization,
     grid, plotgrid, plotgrid_layout, plotvalues_layout, grid_layout, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
     AffineQuasiVector, AffineMap, AbstractWeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
     checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, broadcastbasis_layout,
     plan_grid_transform, plan_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout, plan_transform_layout, _cumsum
-import FastTransforms: Λ, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg
+import FastTransforms: Λ, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg, pochhammer
 import RecurrenceRelationships: forwardrecurrence, forwardrecurrence!, clenshaw, clenshaw!,
-                        check_clenshaw_recurrences
+                        check_clenshaw_recurrences, polynomialtype
 import RecurrenceRelationshipArrays: initiateforwardrecurrence, Clenshaw
 import FastGaussQuadrature: jacobimoment
 
@@ -90,6 +90,8 @@ represents an OP multiplied by its orthogonality weight.
 """
 struct WeightedOPLayout{Lay<:AbstractOPLayout} <: AbstractWeightedBasisLayout end
 
+grid_layout(::WeightedOPLayout, P, n) = grid(unweighted(P), n)
+
 isorthogonalityweighted(::WeightedOPLayout, _) = true
 function isorthogonalityweighted(::AbstractWeightedBasisLayout, wS)
     w,S = arguments(wS)
@@ -247,15 +249,14 @@ grammatrix_layout(::WeightedOPLayout{MappedOPLayout}, P) = grammatrix_layout(Map
 
 OrthogonalPolynomial(w::Weight) =error("Override for $(typeof(w))")
 
-@simplify *(B::Identity, C::OrthogonalPolynomial) = ApplyQuasiMatrix(*, C, jacobimatrix(C))
+@simplify *(B::QuasiDiagonal{<:Any,<:Inclusion}, C::OrthogonalPolynomial) = ApplyQuasiMatrix(*, C, jacobimatrix(C))
 
 function layout_broadcasted(::Tuple{PolynomialLayout,AbstractOPLayout}, ::typeof(*), x::Inclusion, C)
     x == axes(C,1) || throw(DimensionMismatch())
     C*jacobimatrix(C)
 end
 
 
-
 # function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::BroadcastQuasiVector, C::OrthogonalPolynomial)
 #     axes(a,1) == axes(C,1) || throw(DimensionMismatch())
 #     # re-expand in OP basis