sample on rectangles

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.9.3"
+version = "0.9.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,8 +29,8 @@ ArrayLayouts = "1.11"
 BandedMatrices = "1"
 BlockArrays = "1.0"
 BlockBandedMatrices = "0.13"
-ClassicalOrthogonalPolynomials = "0.15"
-ContinuumArrays = "0.19.5"
+ClassicalOrthogonalPolynomials = "0.15.8"
+ContinuumArrays = "0.19.6"
 DomainSets = "0.7"
 FastTransforms = "0.17"
 FillArrays = "1.0"
@@ -40,6 +40,7 @@ InfiniteLinearAlgebra = "0.9, 0.10"
 LazyArrays = "2.3.1"
 LazyBandedMatrices = "0.11.3"
 QuasiArrays = "0.12.2"
+Random = "1.11.0"
 RecurrenceRelationships = "0.2"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"
@@ -48,8 +49,9 @@ julia = "1.10"
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "ForwardDiff", "RecipesBase", "Test"]
+test = ["Base64", "ForwardDiff", "Random", "RecipesBase", "Test"]

src/MultivariateOrthogonalPolynomials.jl

@@ -21,7 +21,7 @@ import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout, PaddedLayo
 import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes, DiagTravLayout, invdiagtrav, ApplyBandedBlockBandedLayout, krontrav
 import InfiniteArrays: InfiniteCardinal, OneToInf
 
-import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw, OPLayout
+import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw, OPLayout, sample, sample_layout
 import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
                                           AngularMomentum, angularmomentum, BlockOneTo, BlockRange1, interlace,
                                           MultivariateOPLayout, AbstractMultivariateOPLayout, MAX_PLOT_BLOCKS

src/rect.jl

@@ -202,4 +202,14 @@ end
 function vec_layout(::KronExpansionLayout{OPLayout, OPLayout}, f)
     A,C,Bc = arguments(f)
     RectPolynomial(A, Bc') * DiagTrav(C)
-end
+end
+
+function sample_layout(::ExpansionLayout{MultivariateOrthogonalPolynomials.KronOPLayout{2}}, f::AbstractQuasiVector, n)
+    F = reshape(f)
+    x = sample(sum(F; dims=2), n)
+    # x = sample(F[:,y]) # TODO: this should work
+    y = [sample(F[x,:]) for x in x]
+    map(SVector, x, y)
+end
+
+sample_layout(lay::ExpansionLayout{MultivariateOrthogonalPolynomials.KronOPLayout{2}}, f::AbstractQuasiVector) = only(sample_layout(lay, f, 1))

test/test_rect.jl

@@ -1,9 +1,11 @@
-using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, ArrayLayouts, Test
-using ClassicalOrthogonalPolynomials: expand, coefficients, recurrencecoefficients
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, ArrayLayouts, Random, Test
+using ClassicalOrthogonalPolynomials: expand, coefficients, recurrencecoefficients, sample
 using MultivariateOrthogonalPolynomials: weaklaplacian, ClenshawKron
 using ContinuumArrays: plotgridvalues, ExpansionLayout
 using Base: oneto
 
+Random.seed!(3242)
+
 @testset "RectPolynomial" begin
     @testset "Evaluation" begin
         T = ChebyshevT()
@@ -161,7 +163,7 @@ using Base: oneto
         f = expand(P, splat((x,y) -> 1))
         @test diff(f,(1,0))[SVector(0.1,0.2)] == diff(f,(0,1))[SVector(0.1,0.2)] == 0.0
         f = expand(P, splat((x,y) -> x))
-        @test diff(f,(1,0))[SVector(0.1,0.2)] == 1.0
+        @test diff(f,(1,0))[SVector(0.1,0.2)] ≈ 1.0
         @test diff(f,(0,1))[SVector(0.1,0.2)] == 0.0
         f = expand(P, splat((x,y) -> cos(x*exp(y))))
         @test diff(f,(1,0))[SVector(0.1,0.2)] ≈ -sin(0.1*exp(0.2))*exp(0.2)
@@ -247,5 +249,8 @@ using Base: oneto
         F = reshape(f)
         @test sum(F; dims=1)[1,0.2] ≈ 2.346737615950585
         @test sum(F; dims=2)[0.1,1] ≈ 2.1748993079723618
+        
+        x,y = coordinates(P)
+        @test sum(sample(f, 100_000))/100_000 ≈ [sum(x .* f)/sum(f),sum(y .* f)/sum(f)] rtol=1E-1
     end
 end