JuliaApproximation / ClassicalOrthogonalPolynomials.jl

A Julia package for classical orthogonal polynomials and expansions
MIT License
38 stars 6 forks source link

Shifted Legendre{BigFloat}() error #28

Open TSGut opened 3 years ago

TSGut commented 3 years ago

I'm seeing some strange behavior for Legendre polynomials with BigFloat precision shifted to 0..1. Quite possible that something about my syntax isn't how it is intended to be used but it might also be bug.

Here's a simple example that reproduces the error:

julia> # float64

julia> P = Legendre()[affine(Inclusion(0..1), axes(Legendre(),1)), :]
Legendre{Float64} affine mapped to 0..1

julia> x = axes(P,1)
Inclusion(0..1)

julia> P \ sin.(x)
ℵ₀-element LazyArrays.CachedArray{Float64,1,Array{Float64,1},Zeros{Float64,1,Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
  0.4596976941318605
  0.42791899124296007
 -0.03924363301542014
 -0.007212191228055754
  0.0002821450243386184
  2.8742703093836133e-5
 -7.146539529763204e-7
 -5.0363463994812964e-8
  9.178337836736593e-10
  4.944523575181289e-11
 -7.110265998998354e-13
 -3.106481507574827e-14
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  ⋮

julia> # bigfloat

julia> Pbig = Legendre{BigFloat}()[affine(Inclusion(BigInt(0)..BigInt(1)), axes(Legendre{BigFloat}(),1)), :]
Legendre{BigFloat} affine mapped to 0..1

julia> xbig = axes(Pbig,1)
Inclusion(0..1)

julia> Pbig \ sin.(xbig)
ERROR: cannot resize array with shared data
Stacktrace:
 [1] _deleteend! at ./array.jl:901 [inlined]
 [2] resize! at ./array.jl:1090 [inlined]
 [3] chop!(::Array{BigFloat,1}, ::BigFloat) at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:71
 [4] adaptivetransform_ldiv(::SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}, ::BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}) at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:104
 [5] transform_ldiv at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:64 [inlined]
 [6] transform_ldiv at /home/timon/.julia/packages/ContinuumArrays/D2vwQ/src/bases/bases.jl:178 [inlined]
 [7] copy at /home/timon/.julia/packages/ContinuumArrays/D2vwQ/src/bases/bases.jl:180 [inlined]
 [8] materialize(::Ldiv{ContinuumArrays.MappedBasisLayout,LazyArrays.BroadcastLayout{typeof(sin)},SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false},BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}}) at /home/timon/.julia/packages/ArrayLayouts/7LZ91/src/ldiv.jl:22
 [9] ldiv at /home/timon/.julia/packages/ArrayLayouts/7LZ91/src/ldiv.jl:86 [inlined]
 [10] \(::SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}, ::BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}) at /home/timon/.julia/packages/QuasiArrays/Kytgm/src/matmul.jl:34
 [11] top-level scope at REPL[77]:1

julia> # but this appears to work

julia> Pbig[:,1:20] \ sin.(xbig)
20-element Array{BigFloat,1}:
  0.459697694131860282599063392557024291715022853422108373904832662310088291651591
  0.4279189912429598877122041074528643054476364481970978558557517258602675648095019
 -0.03924363301542034337341694172731222575488871592208301467640058895035480424366479
 -0.007212191228055742704936278556706084716636823072871192246094339888763058783659867
  0.000282145024338535280236421222451671910237190823901368171354926761601708489205559
  2.874270309358435249957361731260692515947881520105291315102420408358148420958176e-05
 -7.146539530749001084413538229440980152040709571024099798175596892344455384798109e-07
 -5.036346369111653274248451118754811843127054009898515258792902259125311582439218e-08
  9.178336969399690650953052093298209802341083885853023198941221702769575165208465e-10
  4.944521184304420861236925934788833915496435733298095933603120844076503786764353e-11
 -7.112115015695451912101596871391910849417223830934967825666746546919715622639409e-13
 -3.101024774116237937183881067754829492183211029272491093718995995982405455400059e-14
  3.684185721261002718333096225757674226208193004456196754596387429774500100182845e-16
  1.349202245083837855737893815880041795386546710643032360063629130209568728057782e-17
 -1.365328931605316752701716780091452534277698239440733900315583084467897549345079e-19
 -4.310049800103575884451734511728530188560316941533188009045818391107703584309708e-21
  3.798549690418122116191713913664570774313410824956465599814127761088843678732203e-23
  1.053718404267767678750397050992597701028600476518485629180358929672362729284939e-24
 -8.224287969317306291875699136774388048268971744012539584301297942139600288084304e-27
 -2.035023861243751730600511331372270370404905435888710306620612749522021466923814e-28
dlfivefifty commented 3 years ago

This is a weird one I haven't seen before...the issue is clearly that resize! doesn't always work with BigFloat. Here's a MWE showing its probably a bug in StdLib:

julia> resize!(qr(BigFloat.(randn(5,5))) \ BigFloat.(randn(5)),3)
ERROR: cannot resize array with shared data
Stacktrace:
 [1] _deleteend!
   @ ./array.jl:893 [inlined]
 [2] resize!(a::Vector{BigFloat}, nl::Int64)
   @ Base ./array.jl:1109
 [3] top-level scope
   @ REPL[24]:1

See https://github.com/JuliaLang/julia/issues/40108

For now can you try changing adaptivetransform_ldiv to just call an allocating chop (which you might have to write yourself)?

dlfivefifty commented 3 years ago

PS at the moment the transform is falling back to the generic slow transform based on inverting the Vandermonde matrix. We should change it to call FastTransforms.jl