Skip to content

Modifications for linear algebra #18

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.

Sign up for GitHub

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

Closed
wants to merge 3 commits into from
Closed

Modifications for linear algebra #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2b8db3f
Add some methods to support linear algebra
j-fu Jun 9, 2023
1f50c4f
start test code for linear algebra
j-fu Jun 9, 2023
182ed87
Merge branch 'main' into linalg
j-fu Jun 9, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 35 additions & 3 deletions src/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,28 @@ Base.:/(l::Number, r::Quantity) = l * inv(r)
Base.:/(l::Dimensions, r::Number) = Quantity(inv(r), l)
Base.:/(l::Number, r::Dimensions) = Quantity(l, inv(r))

Base.:+(l::Quantity, r::Quantity) = dimension(l) == dimension(r) ? Quantity(l.value + r.value, l.dimensions) : throw(DimensionError(l, r))
Base.:-(l::Quantity, r::Quantity) = dimension(l) == dimension(r) ? Quantity(l.value - r.value, l.dimensions) : throw(DimensionError(l, r))
function Base.:+(l::Quantity, r::Quantity)
if iszero(l)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't you want a DimensionError to still be raised even if iszero is true? If adding a 2D vector to a 3D vector, you would always want to raise an error, regardless of if one of them was zero, no?

It might not even be the case that iszero is defined for all possible Quantity, since the value of a Quantity might be an abstract object, rather than just a scalar real

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The use case is that there are places where a value is initialized with zero, and than += or -= are used to add other quantities to it.

This can be solved also e.g. in Sparspak. But IMHO it this belongs to the algebraic properties of Quantity. What is zero and what is one ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like other unit packages also raise DimensionErrors even if the value is zero. e.g., astropy:

from astropy import units as u

q = 10 * u.m # 10 meters
q2 = 10 * u.s # 10 seconds

q * 0 + q2
# ^ Raises UnitConversionsError

likewise with Unitful.jl

julia> using Unitful

julia> 0.0u"m" + 10.0u"s"
ERROR: DimensionError: 0.0 m and 10.0 s are not dimensionally compatible.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think "zero meters" + "one second" raising an error makes the most sense as a default? Thoughts?

Copy link
Member

@MilesCranmer MilesCranmer Jun 9, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In Sparspak maybe you can do:

x = iszero(x) ? q : x + q

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had it running like this, so for sparspak I certainly could switch to this. But sparspak may not be the only place where this functionality may be needed.

Fundamentally this leads to the question how units are working algebraically. I'll try to think more about it to find a less ad-hoc approach.

IMHO it would be good if DynamicQuantities could be used like Dual numbers, MultiFloats, Measurements and the like as a drop-in replacement for Float64.
But I am aware that this may be not everyone's version for this package.

Copy link
Member

@MilesCranmer MilesCranmer Jun 12, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think even if there are a few additional packages that need this, it is better left out of the main package, because it will slow down every call to + with 3x as many comparisons on the dimensions (i.e., rather than just dimension(l) == dimension(r), now it has to also compute dimension(l) == Dimensions() and dimension(r) == Dimensions()).

return r
elseif iszero(r)
return l
else
l.dimensions == r.dimensions ? Quantity(l.value + r.value, l.dimensions) : throw(DimensionError(l, r))
end
end


function Base.:-(l::Quantity, r::Quantity)
if iszero(l)
return -r
elseif iszero(r)
return l
else
l.dimensions == r.dimensions ? Quantity(l.value - r.value, l.dimensions) : throw(DimensionError(l, r))
end
end

Base.:-(l::Quantity) = Quantity(- l.value, l.dimensions)

_pow(l::Dimensions{R}, r::R) where {R} = @map_dimensions(Base.Fix1(*, r), l)
_pow(l::Quantity{T,R}, r::R) where {T,R} = Quantity(l.value^convert(T, r), _pow(l.dimensions, r))
Expand All @@ -32,4 +52,16 @@ Base.sqrt(q::Quantity) = Quantity(sqrt(q.value), sqrt(q.dimensions))
Base.cbrt(d::Dimensions{R}) where {R} = d^inv(convert(R, 3))
Base.cbrt(q::Quantity) = Quantity(cbrt(q.value), cbrt(q.dimensions))

Base.abs(q::Quantity) = Quantity(abs(q.value), q.dimensions)

#
# We need this for pivoting: we could introduce a pivoting type RowNonZero instead.
#
#Base.abs(q::Quantity) = Quantity(abs(q.value), q.dimensions, q.valid)
Base.abs(q::Quantity) = Quantity(abs(q.value))


Base.iszero(d::Dimensions) = d==Dimensions()
Base.isless(q::Quantity,r::Quantity) = q.value<r.value && q.dimensions==r.dimensions
Base.isone(q::Quantity{T}) where T = isone(q.value) && iszero(q.dimensions)
Base.iszero(q::Quantity{T}) where T = iszero(q.value)

27 changes: 27 additions & 0 deletions test/linsolve.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
module linsolve


using SparseArrays,Sparspak,DynamicQuantities, LinearAlgebra
using GenericLinearAlgebra

function makeproblem(n;dimA=DynamicQuantities.Dimensions(length=1),dimb=DynamicQuantities.Dimensions(tim=1))
A=-sprand(n,n,0.5)+100I
SparseMatrixCSC(size(A)...,A.colptr,A.rowval,Quantity.(A.nzval,dimA)), Quantity.(rand(n),dimb)

# [ Quantity(b0[i],length=1) for i=1:n]
#DynamicQuantities.Quantity.(rand(10),time=1)
end


function densetest(n)
As,b=makeproblem(n)
A=Matrix(As)
lu(A)\b
end

function sparsetest(n)
A,b=makeproblem(n)
sparspaklu(A)\b
end

end