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

Embrace view for subspaces #21

Closed
3 tasks
dlfivefifty opened this issue Sep 27, 2019 · 2 comments
Closed
3 tasks

Embrace view for subspaces #21

dlfivefifty opened this issue Sep 27, 2019 · 2 comments

Comments

@dlfivefifty
Copy link
Member

At the moment the following makes a lazy *, but I think it should be left as a view:

julia> S = LinearSpline{Float64}(0.0:3);

julia> S[:,2:end-1] # Lazy S * Project
QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{Spline{1,Float64},BandedMatrices.BandedMatrix{Int64,FillArrays.Ones{Int64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Base.OneTo{Int64}}}}(*, (Spline{1,Float64}([0.0, 1.0, 2.0, 3.0]), [0 0; 1 0; 0 1; 0 0]))

julia> view(S,:,2:size(S,2)-1) # proposal: returns this
QuasiArrays.SubQuasiArray{Float64,2,Spline{1,Float64},Tuple{Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}(Spline{1,Float64}([0.0, 1.0, 2.0, 3.0]), (Inclusion(0.0..3.0), 2:3), 0, 0)

I believe this is what you were asking for @jagot.

The steps to make it work are

  • Support lazy multiplication with views
  • Fix getindex for lazy multiplication of views
  • Lower D*view(S,:,kr) to view(D*S,:,kr)
@jagot
Copy link
Member

jagot commented Jan 7, 2020

I guess this is related to this issue; I'm adding restrictions to FiniteDifferencesQuasi.jl, and encounter the same problem as mentioned in #37:

julia> N = 10
10

julia> ρ = 1.0
1.0

julia> R = RadialDifferences(N, ρ)
Radial finite differences basis {Float64} on 0.0..10.5 with 10 points spaced by ρ = 1.0

julia> r = axes(R, 1)
Inclusion(0.0..10.5)

julia> D = Derivative(r)
Derivative{Float64,IntervalSets.Interval{:closed,:closed,Float64}}(Inclusion(0.0..10.5))

julia>= R[:,3:8]
Radial finite differences basis {Float64} on 0.0..10.5 with 10 points spaced by ρ = 1.0, restricted to basis functions 3..8  1..10

julia> materialize(applied(*, R̃', D, R))
6×10 BandedMatrices.BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
     -0.533333   0.0        0.514286                                                   
               -0.514286   0.0        0.507937                                         
                         -0.507937   0.0        0.505051                               
                                   -0.505051   0.0        0.503497                     
                                             -0.503497   0.0       0.502564            
                                                       -0.502564  0.0       0.501961   

julia> apply(*, R̃', D, R)
6×10 BandedMatrices.BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
     -0.533333   0.0        0.514286                                                   
               -0.514286   0.0        0.507937                                         
                         -0.507937   0.0        0.505051                               
                                   -0.505051   0.0        0.503497                     
                                             -0.503497   0.0       0.502564            
                                                       -0.502564  0.0       0.501961   

julia>'D*R
6×10 Array{Float64,2}:
 0.0  -0.533333   0.0        0.514286   0.0        0.0        0.0       0.0       0.0       0.0
 0.0   0.0       -0.514286   0.0        0.507937   0.0        0.0       0.0       0.0       0.0
 0.0   0.0        0.0       -0.507937   0.0        0.505051   0.0       0.0       0.0       0.0
 0.0   0.0        0.0        0.0       -0.505051   0.0        0.503497  0.0       0.0       0.0
 0.0   0.0        0.0        0.0        0.0       -0.503497   0.0       0.502564  0.0       0.0
 0.0   0.0        0.0        0.0        0.0        0.0       -0.502564  0.0       0.501961  0.0

My interpretation so far is that this occurs due to fullmaterialize:
https://github.com/JuliaApproximation/QuasiArrays.jl/blob/21d81858e496435176b811148c9b90e0d142219a/src/matmul.jl#L57-L59

Could an alternative approach to the steps you outlined above be to add a trait that decides whether to fullmaterialize or not? In conjunction with JuliaArrays/LazyArrays.jl#91, I would prefer the views to remain attached to the Basis and leave Derivative untouched, i.e. a string of multiplications a*b*c*... should be resolved to a Mul(a,b,c,...) which then gets materialized. Is that a good approach?

@dlfivefifty
Copy link
Member Author

I think this is done.

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