Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Shifted Legendre{BigFloat}() error #28

Open
TSGut opened this issue Mar 19, 2021 · 2 comments
Open

Shifted Legendre{BigFloat}() error #28

TSGut opened this issue Mar 19, 2021 · 2 comments

Comments

@TSGut
Copy link
Member

TSGut commented Mar 19, 2021

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
Copy link
Member

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
JuliaLang/LinearAlgebra.jl#823

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

@dlfivefifty
Copy link
Member

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants