Add Unitful support #244
Add Unitful support #244
Conversation
... in the hopes of supporting Unitful numbers, since they are <: Number. This alone does not work though.
…ls instead of floats This makes working with units more exact, since the powers of Unitful units are always rational. Without this, we get powers such as 74981274124 / 1249141204412 for units.
Running a simple test
using Unitful: @u_str
import Interpolations as Ilib
import Interpolations
import FiniteDifferences as FDlib
tt = (0.0 : 0.01 : 1.1) .* u"s"
fn = 3u"J / s" .* tt
ii = Ilib.interpolate((tt,),fn, ComposedFunction(Ilib.Gridded, Ilib.Linear)())
cdm = FDlib.central_fdm(5,1)
@show dd = cdm(ii, 0.6u"s")
without this change in the order of operations causes the code to crash
in _limit_step, because of a dimension mismatch. This is a temporary "fix",
that actually not a fix. The given order of operations should really
work as-is...
There might still be a few places where units need to be enforced, since comparing quantities and pure numbers does not work, and my manual tests have just not hit those branches yet.
…strip and remove extra .* unit(x)
…it(::Type{Union{Missing,T}})
This is to resolve errors with uncluding units into FiniteDifferences.jl.
See JuliaDiff/FiniteDifferences.jl#244 for details.
Tests for Units.jl need to be updated, if this is to be accepted.
The totalness of unit introduced via unit(::Ay) might be a bit controversial,
but I don't think that should be a problem. In general, things don't have units.
Only Numbers do.
Closes JuliaPhysics#806.
For some reason, the input value x of the fn being differentiated should have the inverse units of the derivative, for the units to come out correctly. This should be lookied into...
…integers (although both work) Again, the input value x to fn needs to have non-sensible units for the desired units to come out of the finite difference method.
…tful tests These seem to at least produce correct units, but the numerical value is wrong for the adapt=0 case (not surprising, based on the documentation). The question is, why the heck do things go wrong, when you do not give an explicit adapt value, but leave it to the compiler-generated method that inserts the value adapt=1 automatically.
…bject instead of a container This means not relying on eltype.
…test Since the commit where the assumption about the function-likeness of differentiated objects was made in step estimation, the correct unit actually works. The step with adaptive=0 was also useless, as a method with no iterations does not converge, of course.
|
Note that this might rely on Unitful.jl merging JuliaPhysics/Unitful.jl#807, although that is not quaranteed, since the propagation of the Edit: nope, the current release of Unitful also works. |
This function itself should be called in a broadcasted manner, if broadcasting is desired.
|
The relaxation of types constraints is fine. I would at least like to see benchmarks showing that they all compile away fully for the case of normal |
|
The benchmarks are definitely still a concern of mine. My guestion is, what would constitute a good benchmark that would satisfy most people? I guess I could |
|
As a quick note, when called with numbers without units, the allocations seem to remain small, while unitful calls currently increase the number of allocations by some amount: I'll have to check whether the allocations could be avoided in the Unitful case, even if they are not increased for the unitless one. Also, note the incorrect units of the derivative of a unitless function Edit: nevermind, the units are correct. Working on multiple things at the same time again... |
|
The easiest thing might be to simply record what the units should be at the start based on the order of the derivative, and then just adding them back in at the end. |
Instead, attach units to step and acc with two separate calls to withUnit, and determine the number type N of method coefs intead of first adding units and then stripping them. A fair amount of allocations still remain for Unitful calls, though.
|
I would benchmark the following FDM methods: because those are the critical ones that people use, |
I think this is a good way to go about it. Once the units of the function value and inputs and the order of the derivative is know, the units of the output can be computed. |
|
Ah, but the problem with only assigning the units at the end is that unit information might be embedded into the function whose derivative is being estimated, as was done above: fn(t) = sin(t / u"t")For the purposes of computing the numerical value of a derivative first and then attaching the units, we can strip unit information from the number Now, I am still a bit baffled as to why there are more allocations in the current code than in the purely numerical one. Attaching and stripping units from scalars alone should not cause allocations, and there should be no dotted function calls in the WIP code anymore. Maybe it has something to do with type instability, since this is the output I get when adding a few debug prints to the code: julia> cdm = central_fdm(5,1)
FiniteDifferenceMethod:
order of method: 5
order of derivative: 1
grid: [-2, -1, 0, 1, 2]
coefficients: [0.08333333333333333, -0.6666666666666666, 0.0, 0.6666666666666666, -0.08333333333333333]
julia> cdm(t -> sin(t/u"s"),pi * u"s")
1
FiniteDifferences.AdaptedFiniteDifferenceMethod = FiniteDifferences.AdaptedFiniteDifferenceMethod
typeof(m) = FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}
x = π s
f = var"#85#86"()
estimate_step
x = 3.141592653589793 s
f(x) = 1.2246467991473532e-16
estimate_magnitudes
x = 3.141592653589793 s
f(x) = 1.2246467991473532e-16
_compute_step_acc_default
x = 3.141592653589793 s
_compute_step_acc
∇f_magnitude = 10.0
f_error = 2.220446049250313e-16
_limit_step
x = 3.141592653589793 s
step = 0.007231988488633704 s
_compute_step_acc_default
x = 3.141592653589793 s
_compute_step_acc
∇f_magnitude = 10.0
f_error = 2.220446049250313e-16
x = 3.141592653589793 s
f(x) = 1.2246467991473532e-16
_compute_estimate
x = 3.141592653589793 s
fs = [0.021694263403528785, 0.014463472655896089, 0.007231925447967751, 1.2246467991473532e-16, -0.007231925447967506, -0.014463472655895844, -0.02169426340352854]
_compute_estimate
x = 3.141592653589793 s
fs = [0.021694263403528785, 0.014463472655896089, 0.007231925447967751, 1.2246467991473532e-16, -0.007231925447967506, -0.014463472655895844, -0.02169426340352854]
_compute_estimate
x = 3.141592653589793 s
fs = [0.021694263403528785, 0.014463472655896089, 0.007231925447967751, 1.2246467991473532e-16, -0.007231925447967506, -0.014463472655895844, -0.02169426340352854]
_compute_step_acc
∇f_magnitude = 1.0000053814049705 s^-5
f_error = 3.469446951953614e-18
_limit_step
x = 3.141592653589793 s
step = 0.0004719677647195769 s
_compute_step_acc_default
x = 3.141592653589793 s
_compute_step_acc
∇f_magnitude = 10.0
f_error = 2.220446049250313e-16
1
FiniteDifferences.AdaptedFiniteDifferenceMethod = FiniteDifferences.AdaptedFiniteDifferenceMethod
step = 0.0004719677647195769 s
2
FiniteDifferences.AdaptedFiniteDifferenceMethod = FiniteDifferences.AdaptedFiniteDifferenceMethod
x = 3.141592653589793 s
f(x) = 1.2246467991473532e-16
_compute_estimate
x = 3.141592653589793 s
fs = [0.0009439353892626215, 0.00047196774719762154, 1.2246467991473532e-16, -0.00047196774719737657, -0.0009439353892623766]
-1.000000000000011 s^-1So initially |
For a given function |
|
I think the problem there is that if # The function f is passed to fdm and captured.
f_wrapper(x) = x -> ustrip(f(x))where |
|
Ok, it might not be an issue: julia> using Revise, Unitful, StaticArrays, BenchmarkTools, FiniteDifferences
julia> fn(t) = sin(t / u"s")
fn (generic function with 1 method)
julia> fnwrap(x) = ustrip(fn(x))
fnwrap (generic function with 1 method)
julia> @benchmark fnwrap(1.0u"s")
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
Range (min … max): 0.833 ns … 3.750 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 0.875 ns ┊ GC (median): 0.00%
Time (mean ± σ): 0.895 ns ± 0.062 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▇ ▁ ▁
█▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁▇▁▁▁▁█▁▁▁▁▁▇▁▁▁▁▁▆▁▁▁▁█▁▁▁▁▁▅▁▁▁▁▁▄▁▁▁▁█ █
0.833 ns Histogram: log(frequency) by time 1.25 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark fnwrap.((1.0u"s", 2u"s"))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
Range (min … max): 0.833 ns … 20.000 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 0.875 ns ┊ GC (median): 0.00%
Time (mean ± σ): 0.890 ns ± 0.194 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▆ ▁
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ █
0.833 ns Histogram: log(frequency) by time 0.959 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark fnwrap.(SVector((1.0u"s", 2u"s")))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
Range (min … max): 0.833 ns … 8.625 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 0.875 ns ┊ GC (median): 0.00%
Time (mean ± σ): 0.892 ns ± 0.097 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▆ ▁
█▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▆▁▁▁▁▁▁▁▆▁▁▁▁▁▁▁▆ █
0.833 ns Histogram: log(frequency) by time 1.12 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark fnwrap.([1.0u"s", 2.0u"s"])
BenchmarkTools.Trial: 10000 samples with 998 evaluations per sample.
Range (min … max): 18.996 ns … 1.018 μs ┊ GC (min … max): 0.00% … 97.05%
Time (median): 19.998 ns ┊ GC (median): 0.00%
Time (mean ± σ): 23.651 ns ± 36.816 ns ┊ GC (mean ± σ): 12.65% ± 7.91%
▃ █▂ ▄
▇█████▄▃▃▄▃▃▃▅██▅▇▇▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▂▂ ▃
19 ns Histogram: frequency by time 29 ns <
Memory estimate: 160 bytes, allocs estimate: 4. |
|
That looks promising! |
Makes the functions in
src/methodsacceptNumbersinstead of justReals ofAbstractFloats. This makes them almost compatible withQuantity <: Numberin Unitful.jl, however units also need to be handled in the function implementations. Towards this end, "dimensional analysis" is performed inside of functions likeestimate_step.Closes #243.