Symbolics.jl support

.github/workflows/ci.yml

@@ -21,7 +21,7 @@ jobs:
       fail-fast: false
       matrix:
         julia-version:
-          - "1.0"
+          - "1.6"
           - "1"
           - "nightly"
         os:

NEWS.md

@@ -1,5 +1,12 @@
 # History of Measurements.jl
 
+## v2.12.0 (????)
+
+### New features
+
+* Support for [symbolic variables](https://symbolics.juliasymbolics.org/stable/manual/variables/),
+   as provided by `Symbolics.jl` package.
+
 ## v2.11.0 (2023-11-03)
 
 ### New features

Project.toml

@@ -15,13 +15,15 @@ BaseType = "7fbed51b-1ef5-4d67-9085-a4a9b26f478c"
 Juno = "e5e0dc1b-0480-54bc-9374-aad01c23163d"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [extensions]
 MeasurementsBaseTypeExt = "BaseType"
 MeasurementsJunoExt = "Juno"
 MeasurementsRecipesBaseExt = "RecipesBase"
 MeasurementsSpecialFunctionsExt = "SpecialFunctions"
+MeasurementsSymbolicsExt = "Symbolics"
 MeasurementsUnitfulExt = "Unitful"
 
 [compat]
@@ -30,7 +32,7 @@ Calculus = "0.4.1, 0.5"
 LinearAlgebra = "<0.0.1, 1"
 RecipesBase = "0.6.0, 0.7, 0.8, 1.0"
 Requires = "0.5.0, 1"
-julia = "1"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
@@ -40,8 +42,9 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["Aqua", "BaseType", "QuadGK", "RecipesBase", "SpecialFunctions", "Statistics", "Test", "Unitful"]
+test = ["Aqua", "BaseType", "QuadGK", "RecipesBase", "SpecialFunctions", "Statistics", "Symbolics", "Test", "Unitful"]

README.md

@@ -52,6 +52,9 @@ easy-to-use calculator.
   [arbitrary precision](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic-1)
   (also called multiple precision) numbers with uncertainties.  This is useful
   for measurements with very low relative error
+* Support for numbers with uncertainties of
+  [symbolic variable](https://symbolics.juliasymbolics.org/stable/manual/variables/)
+  type, as provided by `Symbolics.jl` package.
 * Define arrays of measurements and perform calculations with them.  Some linear
   algebra functions work out-of-the-box
 * Propagate uncertainty for any function of real arguments (including functions

docs/src/appendix.md

@@ -11,7 +11,7 @@ The `Measurement` Type
 `Measurement` is a
 [composite](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1)
 [parametric](https://docs.julialang.org/en/v1/manual/types/#Parametric-Types-1)
-type, whose parameter is the `AbstractFloat` subtype of the nominal value and
+type, whose parameter is the `Real` subtype of the nominal value and
 the uncertainty of the measurement. `Measurement` type itself is subtype of
 `AbstractFloat`, thus `Measurement` objects can be used in any function taking
 `AbstractFloat` arguments without redefining it, and calculation of uncertainty
@@ -20,7 +20,7 @@ will be exact.
 In detail, this is the definition of the type:
 
 ```julia
-struct Measurement{T<:AbstractFloat} <: AbstractFloat
+struct Measurement{T<:Real} <: AbstractFloat
     val::T
     err::T
     tag::UInt64

docs/src/examples.md

@@ -272,6 +272,33 @@ hypot(a, b)
 log(2a) ^ b
 ```
 
+Symbolic Variables
+------------------
+
+You can perform any mathematical operation supported in `Measurements.jl` using
+[symbolic variables](https://symbolics.juliasymbolics.org/stable/manual/variables/),
+as provided by `Symbolics.jl` package:
+
+``` julia
+julia> using Measurements, Symbolics
+
+julia> @variables x_val, x_err, y_val, y_err
+4-element Vector{Num}:
+ x_val x_err y_val y_err
+
+julia> x = x_val ± x_err
+x_val ± x_err
+
+julia> y = y_val ± y_err
+y_val ± y_err
+
+julia> x + y
+x_val + y_val ± sqrt(abs2(x_err) + abs2(y_err))
+
+julia> x - x
+0 ± 0.0
+```
+
 Operations with Arrays and Linear Algebra
 -----------------------------------------
 

docs/src/index.md

@@ -48,6 +48,9 @@ The main features of the package are:
   precision](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic-1)
   (also called multiple precision) numbers with uncertainties. This is useful
   for measurements with very low relative error
+- Support for numbers with uncertainties of
+  [symbolic variable](https://symbolics.juliasymbolics.org/stable/manual/variables/)
+  type, as provided by `Symbolics.jl` package.
 - Define arrays of measurements and perform calculations with them.  Some linear
   algebra functions work out-of-the-box
 - Propagate uncertainty for any function of real arguments (including functions

docs/src/usage.md

@@ -20,11 +20,12 @@ where
 - `err` is its uncertainty, assumed to be a [standard
   deviation](https://en.wikipedia.org/wiki/Standard_deviation).
 
-They are both subtype of `AbstractFloat`. Some keyboard layouts provide an easy
+They are both subtype of `Real`. Some keyboard layouts provide an easy
 way to type the `±` sign, if your does not, remember you can insert it in Julia
 REPL with `\pm` followed by `TAB` key. You can provide `val` and `err` of any
-subtype of `Real` that can be converted to `AbstractFloat`. Thus,
-`measurement(42, 33//12)` and `pi ± 0.1` are valid.
+subtype of `Real` that can be converted to `AbstractFloat`, or are
+[symbolic variables](https://symbolics.juliasymbolics.org/stable/manual/variables/).
+Thus, `measurement(42, 33//12)` and `pi ± 0.1` are valid.
 
 `measurement(value)` creates a `Measurement` object with zero uncertainty, like
 mathematical constants. See below for further examples.

ext/MeasurementsSpecialFunctionsExt.jl

@@ -30,40 +30,40 @@ else
 end
 # Error function: erf, erfinv, erfc, erfcinv, erfcx, erfi, dawson
 
-function SpecialFunctions.erf(a::Measurement{T}) where {T<:AbstractFloat}
+function SpecialFunctions.erf(a::Measurement{T}) where {T<:Real}
     aval = a.val
     return result(erf(aval), 2*exp(-abs2(aval))/sqrt(T(pi)), a)
 end
 
-function SpecialFunctions.erfinv(a::Measurement{T}) where {T<:AbstractFloat}
+function SpecialFunctions.erfinv(a::Measurement{T}) where {T<:Real}
     res = erfinv(a.val)
     # For the derivative, see http://mathworld.wolfram.com/InverseErf.html
     return result(res, sqrt(T(pi)) * exp(abs2(res)) / 2, a)
 end
 
-function SpecialFunctions.erfc(a::Measurement{T}) where {T<:AbstractFloat}
+function SpecialFunctions.erfc(a::Measurement{T}) where {T<:Real}
     aval = a.val
     return result(erfc(aval), -2*exp(-abs2(aval))/sqrt(T(pi)), a)
 end
 
-function SpecialFunctions.erfcinv(a::Measurement{T}) where {T<:AbstractFloat}
+function SpecialFunctions.erfcinv(a::Measurement{T}) where {T<:Real}
     res = erfcinv(a.val)
     # For the derivative, see http://mathworld.wolfram.com/InverseErfc.html
     return result(res, -sqrt(T(pi)) * exp(abs2(res)) / 2, a)
 end
 
-function SpecialFunctions.erfcx(a::Measurement{T}) where {T<:AbstractFloat}
+function SpecialFunctions.erfcx(a::Measurement{T}) where {T<:Real}
     aval = a.val
     res = erfcx(aval)
     return result(res, 2 * (aval * res - inv(sqrt(T(pi)))), a)
 end
 
-function SpecialFunctions.erfi(a::Measurement{T}) where {T<:AbstractFloat}
+function SpecialFunctions.erfi(a::Measurement{T}) where {T<:Real}
     aval = a.val
     return result(erfi(aval), 2*exp(abs2(aval))/sqrt(T(pi)), a)
 end
 
-function SpecialFunctions.dawson(a::Measurement{T}) where {T<:AbstractFloat}
+function SpecialFunctions.dawson(a::Measurement{T}) where {T<:Real}
     aval = a.val
     res = dawson(aval)
     return result(res, one(T) - 2 * aval * res, a)

ext/MeasurementsSymbolicsExt.jl

@@ -0,0 +1,30 @@
+### MeasurementsSymbolicsExt.jl
+#
+# Copyright (C) 2023 Roman Lebedev.
+#
+# Maintainer: Mosè Giordano <mose AT gnu DOT org>
+# Keywords: uncertainty, error propagation, physics
+#
+# This file is a part of Measurements.jl.
+#
+# License is MIT "Expat".
+#
+### Commentary:
+#
+# This file contains integration with Symbolics.jl.
+#
+### Code:
+
+### Symbolics
+module MeasurementsSymbolicsExt
+
+import Measurements
+import Symbolics
+
+Measurements._is_symbolic(::Symbolics.Num) = true
+
+Measurements.measurement(val::Symbolics.Num) = Measurements._measurement(val)
+
+Measurements.measurement(val::Symbolics.Num, err::Symbolics.Num) = Measurements._measurement(val, err)
+
+end

src/Measurements.jl

@@ -27,6 +27,9 @@ using Calculus
 # Functions provided by this package and exposed to users
 export Measurement, measurement, ±
 
+# Query whether the value is of symbolic type.
+_is_symbolic(::Real) = false
+
 # Define the "Derivatives" type, used inside "Measurement" type.  This should be
 # a lightweight and immutable dictionary.
 include("derivatives-type.jl")
@@ -46,7 +49,7 @@ include("derivatives-type.jl")
 #     measurement and propagate the uncertainty in the case of functions with
 #     more than one argument (in order to deal with correlation between
 #     arguments).
-struct Measurement{T<:AbstractFloat} <: AbstractFloat
+struct Measurement{T<:Real} <: AbstractFloat
     val::T
     err::T
     tag::UInt64
@@ -72,16 +75,19 @@ function Measurement{T}(::S) where {T, S}
 end
 
 # Functions to quickly create an empty Derivatives object.
-@generated empty_der1(x::Measurement{T}) where {T<:AbstractFloat} = Derivatives{T}()
-@generated empty_der2(x::T) where {T<:AbstractFloat} = Derivatives{x}()
+@generated empty_der1(x::Measurement{T}) where {T<:Real} = Derivatives{T}()
+@generated empty_der2(x::T) where {T<:Real} = Derivatives{x}()
 
 # Start from 1, 0 is reserved to derived quantities
 const tag_counter = Threads.Atomic{UInt64}(1)
 
 measurement(x::Measurement) = x
-measurement(val::T) where {T<:AbstractFloat} = Measurement(val, zero(T), UInt64(0), empty_der2(val))
+
+_measurement(val::T) where {T<:Real} =  Measurement(val, zero(T), UInt64(0), empty_der2(val))
+measurement(val::AbstractFloat) = _measurement(val)
 measurement(val::Real) = measurement(float(val))
-function measurement(val::T, err::T) where {T<:AbstractFloat}
+
+function _measurement(val::T, err::T) where {T<:Real}
     newder = empty_der2(val)
     if iszero(err)
         Measurement{T}(val, err, UInt64(0), newder)
@@ -90,8 +96,12 @@ function measurement(val::T, err::T) where {T<:AbstractFloat}
         return Measurement{T}(val, err, tag, Derivatives(newder, (val, err, tag)=>one(T)))
     end
 end
-measurement(val::Real, err::Real) = measurement(promote(float(val), float(err))...)
+measurement(val::T, err::T) where {T<:AbstractFloat} = _measurement(val, err)
+measurement(val::T, err::T) where {T<:Real} = measurement(float(val), float(err))
+measurement(val::V, err::E) where {V<:Real, E<:Real} = measurement(promote(val, err)...)
+
 measurement(::Missing, ::Union{Real,Missing} = missing) = missing
+
 const ± = measurement
 
 """

src/conversions.jl

@@ -16,15 +16,15 @@
 #
 ### Code:
 
-Base.convert(::Type{Measurement{T}}, a::Irrational) where {T<:AbstractFloat} =
+Base.convert(::Type{Measurement{T}}, a::Irrational) where {T<:Real} =
     measurement(T(a))::Measurement{T}
-Base.convert(::Type{Measurement{T}}, a::Rational{<:Integer}) where {T<:AbstractFloat} =
+Base.convert(::Type{Measurement{T}}, a::Rational{<:Integer}) where {T<:Real} =
     measurement(T(a))::Measurement{T}
-Base.convert(::Type{Measurement{T}}, a::Real) where {T<:AbstractFloat} =
+Base.convert(::Type{Measurement{T}}, a::Real) where {T<:Real} =
     measurement(T(a))::Measurement{T}
-Base.convert(::Type{Measurement{T}}, a::Base.TwicePrecision) where {T<:AbstractFloat} =
+Base.convert(::Type{Measurement{T}}, a::Base.TwicePrecision) where {T<:Real} =
     measurement(T(a))::Measurement{T}
-Base.convert(::Type{Measurement{T}}, a::AbstractChar) where {T<:AbstractFloat} =
+Base.convert(::Type{Measurement{T}}, a::AbstractChar) where {T<:Real} =
     measurement(T(a))::Measurement{T}
 
 function Base.convert(::Type{Measurement{T}}, a::Complex) where {T}
@@ -35,9 +35,9 @@ function Base.convert(::Type{Measurement{T}}, a::Complex) where {T}
     end
 end
 
-Base.convert(::Type{Measurement{T}}, a::Measurement{T}) where {T<:AbstractFloat} = a
+Base.convert(::Type{Measurement{T}}, a::Measurement{T}) where {T<:Real} = a
 function Base.convert(::Type{Measurement{T}},
-                      a::Measurement{<:AbstractFloat}) where {T<:AbstractFloat}
+                      a::Measurement{<:Real}) where {T<:Real}
     newder = empty_der2(zero(T))
     for tag in keys(a.der)
         newder = Derivatives(newder, (T(tag[1]), T(tag[2]), tag[3])=>T(a.der[tag]))
@@ -55,10 +55,10 @@ function Base.convert(::Type{Int}, a::Measurement)
     return convert(Int, a.val)::Int
 end
 
-Base.promote_rule(::Type{Measurement{T}}, ::Type{S}) where {T<:AbstractFloat, S<:Real} =
+Base.promote_rule(::Type{Measurement{T}}, ::Type{S}) where {T<:Real, S<:Real} =
     Measurement{promote_type(T, S)}
 Base.promote_rule(::Type{Measurement{T}},
-                  ::Type{Measurement{S}}) where {T<:AbstractFloat, S<:AbstractFloat} =
+                  ::Type{Measurement{S}}) where {T<:Real, S<:Real} =
     Measurement{promote_type(T, S)}
 
 # adaptation of JuliaLang/julia#30952