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.

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

Conversation

j-fu
Copy link
Contributor

@j-fu j-fu commented Jun 9, 2023

Figure out what is needed to solve linear systems with dynamic quantities.

@j-fu
Copy link
Contributor Author

j-fu commented Jun 9, 2023

Not yet ready for merge.

We might need:

  • support adding/subtracting from zero(Quantity).
  • should abs be stripped of dimensions ? this would be needed for pivoting. However one could use RowNonZero pivoting. This needs to be decided.
  • How to route to the right implementation for dense linalg ?

@github-actions
Copy link
Contributor

github-actions bot commented Jun 9, 2023

Benchmark Results

main 182ed87... t[main]/t[182ed87...]
creation/Quantity(x) 3.7 ± 0.1 ns 3.8 ± 0.2 ns 0.974
creation/Quantity(x, length=y) 5.1 ± 0.1 ns 4.5 ± 0.1 ns 1.13
time_to_load 0.0822 ± 0.0056 s 0.0826 ± 0.00028 s 0.995
with_numbers/*real 4.3 ± 0.1 ns 4.1 ± 0.1 ns 1.05
with_numbers/^int 0.247 ± 0.012 μs 0.244 ± 0.012 μs 1.01
with_numbers/^int * real 0.244 ± 0.011 μs 0.241 ± 0.012 μs 1.01
with_quantity/+y 8.1 ± 0.001 ns 13.7 ± 0.1 ns 0.591
with_quantity//y 0.272 ± 0.0059 μs 0.273 ± 0.0043 μs 0.999
with_self/dimension 2 ± 0.2 ns 2.1 ± 0.1 ns 0.952
with_self/inv 10.7 ± 0.2 ns 10.7 ± 0.001 ns 1
with_self/ustrip 2 ± 0.1 ns 2.1 ± 0.1 ns 0.952

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

@j-fu
Copy link
Contributor Author

j-fu commented Jun 9, 2023

Sparspak works with PetrKryslUCSD/Sparspak.jl#25

Base.:+(l::Quantity, r::Quantity) = Quantity(l.value + r.value, l.dimensions, l.valid && r.valid && l.dimensions == r.dimensions)
Base.:-(l::Quantity, r::Quantity) = Quantity(l.value - r.value, l.dimensions, l.valid && r.valid && l.dimensions == r.dimensions)
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

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

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()).

@MilesCranmer
Copy link
Member

Not sure if this is relevant, but I've been experimenting with this "WildcardDimensionWrapper" type in SymbolicRegression.jl: https://github.com/MilesCranmer/SymbolicRegression.jl/blob/6ed49f166483ec130bc0350d69d6e2214666aeda/ext/SymbolicRegressionUnitfulExt.jl#L58-L69

It's basically a way of adding a regexp to a Quantity, so that you can solve for unknown constants (with their own physical dimensions).

The way it works is as follows. Any free constant (here: a or b) produces a "wildcard" (regexp Quantity) that can cancel out or create physical dimensions wherever necessary.

  • *//
    • Always valid, and propagates the wildcard (since a*x has arbitrary dimension, no matter the dimensions of x).
  • +/-
    • No wildcards => need same dimension
    • 1 wildcard => always valid, but loses the wildcard
    • 2 wildcards => always valid, and propagates a wildcard (since a*x + b*y is of arbitrary physical dimensions)
  • ^
    • No wildcard in power => need zero dimension in power
    • Wildcard in power => Always valid, but loses the wildcard
    • (Always propagates wildcard from base, since ((a*x)^(...) if of arbitrary physical dimension)

sqrt/cbrt/abs propagate wildcards.

Any other operator, if undefined for Quantity, would consume a wildcard. So you can have things like cos(a*x) be dimensionally valid for any x, because a free constant a produces a wildcard.

@j-fu
Copy link
Contributor Author

j-fu commented Jun 11, 2023

Just got the following idea: DynamicalQuantities provides
an abstract type AbstractDynamicalQuantity such that DynanmicalQuantities.Quantity <: AbstractDynamicalQuantity. All operations in math.jl would be defined for AbstractDynamicalQuantity,
and the modifications for linear algebra I am proposing could be implemented for another subtype with some room for experimentation.

@MilesCranmer
Copy link
Member

I think that’s a good idea — though perhaps just AbstractQuantity. And I suppose AbstractDimension would be nice as well in case someone wants to add or remove specific dimensions.

@j-fu
Copy link
Contributor Author

j-fu commented Jun 14, 2023

Possibly can be closed in favor of #24.

@MilesCranmer
Copy link
Member

Maybe we can put the XQuantity you are interested in inside a unittest for the abstract interface?

@MilesCranmer
Copy link
Member

I think most of these are superseded or already implemented now? Closing to clean the repo up

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

Successfully merging this pull request may close these issues.

2 participants