Keno / SIUnits.jl

Efficient unit-checked computation
Other
70 stars 26 forks source link

Linear least-squares fits for vectors with units #112

Open rawlik opened 7 years ago

rawlik commented 7 years ago

I am trying to make it possible to make linear least-squares fits:

x = linspace(0, 50)Volt
y = randn(length(x))Ampere
[x ones(x)] \ y

The first step was #111. Now I get an error:

ERROR: MethodError: Cannot `convert` an object of type SIUnits.SIQuantity{Float64,4,2,-6,-2,0,0,0,0,0} to an object of type Float64
This may have arisen from a call to the constructor Float64(...),
since type constructors fall back to convert methods.
 in generic_vecnorm2(::SubArray{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},1,Array{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},2},Tuple{UnitRange{Int64},Int64},true}) at ./linalg/generic.jl:130
 in vecnorm(::SubArray{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},1,Array{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},2},Tuple{UnitRange{Int64},Int64},true}, ::Int64) at ./linalg/generic.jl:185
 in indmaxcolumn(::SubArray{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},2,Array{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}) at ./linalg/qr.jl:41
 in qrfactPivotedUnblocked!(::Array{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},2}) at ./linalg/qr.jl:60
 in qrfact!(::Array{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},2}, ::Type{Val{true}}) at ./linalg/qr.jl:94
 in qrfact(::Array{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},2}, ::Type{T}) at ./linalg/qr.jl:164
 in \(::Array{SIUnits.SIQuantity{Float64,2,1,-3,-1,0,0,0,0,0},2}, ::Array{SIUnits.SIQuantity{Float64,0,0,0,1,0,0,0,0,0},1}) at ./linalg/generic.jl:363

So I tried defining:

import SIUnits.SIQuantity
import Base.vecnorm
function vecnorm{T,m,kg,s,A,K,mol,cd,rad,sr}(B::AbstractArray{SIQuantity{T,m,kg,s,A,K,mol,cd,rad,sr}}, p::Real = 2)
    valvector = [x.val for x in B]
    valnorm = vecnorm(valvector, p)
    SIQuantity{T,p*m,p*kg,p*s,p*A,p*K,p*mol,p*cd,p*rad,p*sr}(valnorm)
end

Now this works:

julia> vecnorm(x)
205.16295166207246 kg²m⁴s⁻⁶A⁻²
julia> vecnorm(@view x[2:end])
205.16295166207246 kg²m⁴s⁻⁶A⁻²

But still it crashes with the same error when I try [x ones(x)] \ y. I don't understand why is still vecnorm from ./linalg/generic.jl:185 called.

tomasaschan commented 7 years ago

Does the same error show if you restart the Julia session (or start a new one), then define the new vecnorm before trying to [x ones(x)] \ y?