Cannot normalize MulQuasiArrays about continuumarrays.jl HOT 7 CLOSED

Does using broadcast notation help?

By the way, the choice of norm should be dictated by the axis. That is, if the axis is <:AbstractUnitRange its the standard norm, but if the axis is a float range it's as you described. This means that finite differences are accessed by floats not integers, and so to get the kth entry ones should do v[axes(v,1)[k]].

from continuumarrays.jl.

If you mean simply putting a dot in front of the /=, no it does not help, i.e. neither
v.mul.factors[2] ./= norm(v,p)
nor
v.mul.factors[2][:] ./= norm(v,p)
works.

Can I dispatch on axis? Or would it simple be an if axis isa AbstractUnitRange in the norm function?

Lastly, I am unsure what you meant with that I should access as v[axes(v,1)[k]]; which setting are you referring to?

I do have some troubles in Atoms.jl where I represent the different orbitals of the electrons as columns in a MulQuasiArray{T,2,<:Mul{<:Tuple,<:Tuple{<:B,<:AbstractMatrix}}}. I want to normalize each column separately, but if I index as atom.radial_orbitals[:,j], I get ContinuumArrays.QuasiArrays.SubQuasiArray back, which seems nasty to dispatch on. Instead, I index as atom.radial_orbitals[:,j:j]. However, it could be that my problem stems from here. What do you think?

from continuumarrays.jl.

Can I dispatch on axis?

In QuasiArrays module add something like

_norm(::OneTo, x) = norm(collect(x))
norm(x::AbstractQuasiArray) = _norm(axes(x), x)

and then you can override _norm based on the axes type. Though perhaps it doesn't help here since you are using the structure.

Lastly, I am unsure what you meant with that I should access as v[axes(v,1)[k]]; which setting are you referring to?

I was thinking that v lives naturally at the grid point. so v[0.1] would return the value at 0.1 if that's a grid point. Alternatively, the axis could still be OneTo but "weighted" to indicate that the inner product is different, this would support classical indexing but also change the norm.

if I index as atom.radial_orbitals[:,j], I get ContinuumArrays.QuasiArrays.SubQuasiArray

Indexing works by materializing a view:

from continuumarrays.jl.

I realized I should be passing a view of the MulQuasiMatrix to normalize!, so it works now: https://github.com/JuliaAtoms/Atoms.jl/blob/master/src/hydrogenic.jl#L78. Previously, I had been passing a slice which then of course was copied in the function call. Maybe we can close this issue now. However regarding your points above, I do think it would be worthwhile to have a general implementation in ContinuumArrays.QuasiArrays, and let different packages just provide a weight of the axis. This would not be enough for finite differences with varying discretization along the axis, but in that case you'd not be using a general implementation of an inner product anyway.

I tried your suggestion

_norm(::OneTo, x) = norm(collect(x))
norm(x::AbstractQuasiArray) = _norm(axes(x), x)

but collect fails with the following error:

┌ Warning: `length(i::ClosedInterval)` is deprecated, use `IntervalSets.duration(i)` instead.
│   caller = map(::typeof(length), ::Tuple{IntervalSets.Interval{:closed,:closed,Float64},Base.OneTo{Int64}}) at tuple.jl:163
└ @ Base ./tuple.jl:163
ERROR: MethodError: no method matching duration(::IntervalSets.Interval{:closed,:closed,Float64})
Closest candidates are:
  duration(::IntervalSets.TypedEndpointsInterval{:closed,:closed,Dates.Date}) at /Users/jagot/.julia/packages/IntervalSets/xr34V/src/IntervalSets.jl:201
  duration(::IntervalSets.TypedEndpointsInterval{:closed,:closed,T<:Integer}) where T<:Integer at /Users/jagot/.julia/packages/IntervalSets/xr34V/src/IntervalSets.jl:200
Stacktrace:
 [1] length(::IntervalSets.Interval{:closed,:closed,Float64}) at ./deprecated.jl:55
 [2] map(::typeof(length), ::Tuple{IntervalSets.Interval{:closed,:closed,Float64},Base.OneTo{Int64}}) at ./tuple.jl:163
 [3] size(::MulQuasiArray{Float64,2,LazyArrays.Mul{Tuple{LazyArrays.UnknownLayout,LazyArrays.DenseColumnMajor},Tuple{RadialDifferences{Float64,Int64},Array{Float64,2}}}}) at /Users/jagot/.julia/dev/ContinuumArrays/src/QuasiArrays/matmul.jl:78
 [4] length at /Users/jagot/.julia/dev/ContinuumArrays/src/QuasiArrays/abstractquasiarray.jl:127 [inlined]
 [5] _similar_for(::UnitRange{Int64}, ::Type, ::MulQuasiArray{Float64,2,LazyArrays.Mul{Tuple{LazyArrays.UnknownLayout,LazyArrays.DenseColumnMajor},Tuple{RadialDifferences{Float64,Int64},Array{Float64,2}}}}, ::Base.HasLength) at ./array.jl:532
 [6] _collect(::UnitRange{Int64}, ::MulQuasiArray{Float64,2,LazyArrays.Mul{Tuple{LazyArrays.UnknownLayout,LazyArrays.DenseColumnMajor},Tuple{RadialDifferences{Float64,Int64},Array{Float64,2}}}}, ::Base.HasEltype, ::Base.HasLength) at ./array.jl:563
 [7] collect(::MulQuasiArray{Float64,2,LazyArrays.Mul{Tuple{LazyArrays.UnknownLayout,LazyArrays.DenseColumnMajor},Tuple{RadialDifferences{Float64,Int64},Array{Float64,2}}}}) at ./array.jl:557
 [8] top-level scope at none:0

I guess collect(::MulQuasiArray) should return the coefficients of the basis functions in the AbstractQuasiMatrix, however, how many implementation of collect would we need? I can think of these as useful (B isa AbstractQuasiMatrix, v isa VecOrMat):

• collect(B'v) -> v
• collect(v'B) -> v'

from continuumarrays.jl.

I guess collect(::MulQuasiArray) should return the coefficients of the basis functions

Definitely not. collect should only be used for discrete axes with countable entries. For the continuous case, norm would need to be overriden.

from continuumarrays.jl.

Can you try again? There's a function that needs to be overriden for special types so that v[:] .= 4 actually modifies v, that is, becomes equivalent to view(v,:) .= 4. I forgot the name of this function, but I have the feeling I added it....

EDIT: Accidentally sent this! The function is dotview and I believe I added it:

https://discourse.julialang.org/t/overload-x-0-0-5-0-5-1-when-x-is-not-abstractvector/26980

from continuumarrays.jl.

I am not sure if this is exactly what you meant, but I believe it has been working for some time now,
without broadcast notation:

_norm(B::AbstractFiniteDifferences, c::AbstractArray, p::Real=2) =
    norm(c, p)*(step(B)^(inv(p)))

LinearAlgebra.norm(v::FDVecOrMat, p::Real=2) = _norm(v.args..., p)
LinearAlgebra.norm(v::Mul{<:Any, <:Tuple{<:AbstractFiniteDifferences, <:AbstractArray}},
                   p::Real=2) = _norm(v.args..., p)

function LinearAlgebra.normalize!(v::FDVecOrMat, p::Real=2)
    v.args[2][:] /= norm(v,p)
    v
end

function LinearAlgebra.normalize!(v::Mul{<:Any, <:Tuple{<:AbstractFiniteDifferences, <:AbstractArray}},
                                  p::Real=2)
    v.args[2][:] /= norm(v,p)
    v
end

In Atoms.jl, I switched to passing views to normalize!, but now it seems to work with normal arrays as well:

julia> using FiniteDifferencesQuasi, LinearAlgebra
[ Info: Recompiling stale cache file /Users/jagot/.julia/compiled/v1.2/FiniteDifferencesQuasi/QNWLo.ji for FiniteDifferencesQuasi [2e420920-f185-11e8-29af-4b87f0732ba2]

julia> R = FiniteDifferences(10, 0.1)
Finite differences basis {Float64} on 0.0..1.1 with 10 points spaced by Δx = 0.1

julia> M = R*rand(size(R,2),10)
QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{FiniteDifferences{Float64,Int64},Array{Float64,2}}}(*, (Finite differences basis {Float64} on 0.0..1.1 with 10 points spaced by Δx = 0.1, [0.8471444653563867 0.12774328892029807 … 0.07710046830927775 0.758295126787923; 0.42545721313797813 0.11192096320993916 … 0.8304955161873768 0.7894869955582433; … ; 0.43796248119348435 0.06855217265419622 … 0.1494057235848587 0.4627325180866715; 0.24401972006415673 0.1884716617424358 … 0.17941235320373683 0.7429602042961412]))

julia> norm(M)
1.8956676228306937

julia> normalize!(M)
QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{FiniteDifferences{Float64,Int64},Array{Float64,2}}}(*, (Finite differences basis {Float64} on 0.0..1.1 with 10 points spaced by Δx = 0.1, [0.4468844934384613 0.067386965616655 … 0.04067193393014117 0.4000148114866274; 0.22443660904155066 0.05904039392877001 … 0.4381018624706185 0.4164691035759459; … ; 0.23103337099754842 0.03616254866020823 … 0.07881430361814143 0.24410002709003337; 0.12872495005204332 0.09942231405577402 … 0.09464335996614769 0.39192535408012114]))

julia> norm(M)
1.0

from continuumarrays.jl.