Fixes issues #6, #7 by extending power, misc fixes

changelog.org

@@ -1,3 +1,16 @@
+* v0.2.2
+- export ~to~ procedure (used to convert inner type of a
+  ~Measurement~)
+- fix ~**~, ~^~ for two Measurements
+- allow to disable printing of unicode symbols for string repr of
+  ~Measurement~ by compiling with ~-d:noUnicode~
+- fix issue #6 to allow ~measurement~ procedure to construct correct
+  type for given types other than ~float~
+- fix division of a literal with a measurement
+- fix power math operation for different argument types, including
+  fixing issue #7
+  - adds a new private submodule for a helper macro to generate
+    `power` calls for static integers (including negative)
 * v0.2.1
 - change import path in test file to relative import
 * v0.2.0

measuremancer.nim

@@ -18,17 +18,52 @@ type
   FloatLike* = concept x
     x.toFloat is float
 
+  ## A concept for any float like that is *not* an integer type
+  FloatLikeNotInt* = concept x, type T
+    x.toFloat is float
+    not (T is SomeInteger)
+
+  ## A concept for types that are float like and whose application of
+  ## `pow` compiles and yields the same type. E.g. `float32`, but not
+  ## an `unchained` unit (does not support `pow`)
+  FloatLikeSupportsPow* = concept x, type T
+    x.toFloat is float
+    pow(x, 1.0) is T
+
+  ## A concept for types that are float like and whose application of
+  ## `^` with a RT integer value compiles and yields the same type. E.g. `float32`,
+  ## but not an `unchained` unit (needs a `static int` to know the resulting
+  ## type at CT)
+  FloatLikeSupportsIntegerPowRT* = concept x, type T
+    x.toFloat is float
+    (x ^ 1) is T
+
   IdType = uint64 # "`uint64` should be enough for everyone... uhhh"
 
   DerivKey[T] = tuple[val, uncer: T, tag: IdType]
 
   Derivatives[T] = OrderedTable[DerivKey[T], T]
 
-  Measurement*[T: FloatLike] = object
-    val: T
-    uncer: T
-    id: IdType
-    der: Derivatives[T] # map of the derivatives
+  ## *NOTE*: When dealing with `unchained` units, the derivative returned has the
+  ## wrong units! Instead of `T` it should be the type of `∂a(x)/∂x`, but
+  ## the code gets too complicated for the time being, as it would require use to
+  ## store different types for the actual measurement and the gradient, turning
+  ## `Measurement` into a double generic.
+when (NimMajor, NimMinor, NimPatch) < (1, 7, 0):
+  type
+    Measurement*[T] = object
+      val: T
+      uncer: T
+      id: IdType
+      der: Derivatives[T] # map of the derivatives
+else:
+  type
+    Measurement*[T: FloatLike] = object
+      val: T
+      uncer: T
+      id: IdType
+      der: Derivatives[T] # map of the derivatives
+
 
 func value*[T: FloatLike](m: Measurement[T]): T {.inline.} = m.val
 func uncertainty*[T: FloatLike](m: Measurement[T]): T {.inline.} = m.uncer
@@ -70,7 +105,11 @@ proc `==`*[T: FloatLike](k1, k2: DerivKey[T]): bool =
 ## value and uncertainty.
 ## TODO: should we use `almostEqual` here? Or is that too imprecise for the purpose?
 proc `==`*[T: FloatLike](m1, m2: Measurement[T]): bool =
-  result = almostEqual(m1.val, m2.val) and almostEqual(m1.uncer, m2.uncer)
+  ## comparison of two measurements does not need to take into account the
+  ## type, as we require same type anyway. Hence use `toFloat` to compare as float
+  bind toFloat # bind the locally defined `toFloat`
+  result = almostEqual(m1.val.toFloat, m2.val.toFloat) and
+           almostEqual(m1.uncer.toFloat, m2.uncer.toFloat)
 
 proc `==`*[T: FloatLike](m: Measurement[T], x: T): bool =
   result = almostEqual(m.val, x) and m.uncer == 0.0
@@ -82,6 +121,7 @@ proc isInf*[T: FloatLike](m: Measurement[T]): bool = m.val == Inf
 proc isNegInf*[T: FloatLike](m: Measurement[T]): bool = m.val == -Inf
 
 proc toFloat*[T: SomeFloat](x: T): float = x.float # float is biggest type, so this should be fine
+proc toFloat*[T: SomeInteger](x: T): float = x.float # float is biggest type, so this should be fine
 
 proc initDerivatives[T](): Derivatives[T] = initOrderedTable[DerivKey[T], T]()
 
@@ -108,9 +148,10 @@ proc procRes[T](res: T, grad: T, arg: Measurement[T]): Measurement[T] =
       der[key] = (grad * val).T # force to `T` to avoid `unchained` errors
   # σ is 0 if input's σ is zero, else it's ∂(f(x))/∂x * Δx =^= grad * arg.uncer
   let σ = if arg.uncer == T(0.0): T(0.0) else: abs((grad * arg.uncer).T) # force to `T` to avoid `unchained` errors
-  result = initMeasurement[T](res, σ, der, 0'u32)
+  result = initMeasurement[T](res, σ, der, 0'u64)
 
 proc derivative[T](a: Measurement[T], key: DerivKey[T]): T =
+  ## See note about derivative in definition of `Measurement` type
   result = if key in a.der: a.der[key] else: T(0.0)
 
 ## A "type safe" solution required for `unchained` could be achieved with something like this.
@@ -138,7 +179,7 @@ when false:
               # add (σ_x•∂G/∂x)² to total uncertainty (squared), but only if deriv != 0
               resder[key] = ∂G_∂x
               err = err + ((σ_x * ∂G_∂x) * (σ_x * ∂G_∂x)) # convert product back to avoid `unchained` error
-    result = initMeasurement[T](res, sqrt(err), resder, 0'u32)
+    result = initMeasurement[T](res, sqrt(err), resder, 0'u64)
 
 proc procRes[T](res: T, grad: openArray[T], args: openArray[Measurement[T]]): Measurement[T] =
   ## Note: In the body of this we perform rather funny back and forth conversions between `T`
@@ -169,7 +210,7 @@ proc procRes[T](res: T, grad: openArray[T], args: openArray[Measurement[T]]): Me
             err = (err + ((σ_x * ∂G_∂x) * (σ_x * ∂G_∂x)).T).T # convert product back to avoid `unchained` error
   result = initMeasurement[T](res,
                               T(sqrt(err.float)), # convert to float and back to T to satisfy `unchained`
-                              resder, 0'u32)
+                              resder, 0'u64)
 
 proc `±`*[T: FloatLike](val, uncer: T): Measurement[T] =
   result = initMeasurement[T](val, uncer)
@@ -179,12 +220,18 @@ proc `±`*[T: FloatLike](val, uncer: T{lit}): Measurement[float] =
   result = initMeasurement[float](val.float, uncer.float)
 
 proc measurement*[T: FloatLike](value, uncertainty: T): Measurement[T] =
-  result = initMeasurement[float](value, uncertainty)
+  result = initMeasurement[T](value, uncertainty)
+
+proc measurement*[T: FloatLike](val, uncer: T{lit}): Measurement[float] =
+  result = initMeasurement[float](val.float, uncer.float)
 
 proc pretty*[T: FloatLike](m: Measurement[T], precision: int): string =
   let mval = m.val.float.formatBiggestFloat(precision = precision)
   let merr = m.uncer.float.formatBiggestFloat(precision = precision)
-  result = &"{mval} ± {merr}"
+  when not defined(noUnicode):
+    result = &"{mval} ± {merr}"
+  else:
+    result = &"{mval} +- {merr}"
   when not (T is float):
     result.add " " & $T
 
@@ -225,7 +272,7 @@ proc `-=`*[T: FloatLike](a: var Measurement[T], b: Measurement[T]) =
 
 
 ## Type conversion. TODO: make this more type safe, funnily enough
-proc to[T: FloatLike; U](m: Measurement[T], dtype: typedesc[U]): Measurement[U] =
+proc to*[T: FloatLike; U](m: Measurement[T], dtype: typedesc[U]): Measurement[U] =
   when T is U:
     result = m
   else:
@@ -290,29 +337,51 @@ proc `/`*[T: FloatLike; U: FloatLike](m: Measurement[T], x: U): auto =
   assign1(m.val / x, 1.0 / x, m)
 
 ## Overloads for literals that force same type as Measurement has
-proc `/`*[T: FloatLike; U: FloatLike](x: T{lit}, m: Measurement[U]): Measurement[U] =
-  result = procRes(U(U(x) / m.val), U(-x) / (m.val * m.val), m)
+proc `/`*[T: FloatLike; U: FloatLike](x: T{lit}, m: Measurement[U]): auto =
+  type A = typeof( x / m.val )
+  let arg: A = A(x / m.val)
+  let grad: A =  A(-x / (m.val * m.val))
+  result = procRes(arg, grad, m.to(A))
 proc `/`*[U: FloatLike; T: FloatLike](m: Measurement[U], x: T{lit}): Measurement[U] =
   result = procRes(U(m.val / U(x)), 1.0 / U(x), m)
 
 # Power `^`
 ## NOTE: Using any of the following exponentiation functions is dangerous. The Nim parser
 ## might include an additional prefix into the argument to `^` instead of keeping it as a,
 ## well, prefix. So `-(x - μ)^2` is parsed as `(-(x - μ))^2`!
-proc `**`*[T: FloatLike](m: Measurement[T], p: Natural): Measurement[T] =
-  result = procRes(m.val ^ p, p.float * m.val ^ (p - 1), m)
-proc `^`*[T: FloatLike](m: Measurement[T], p: Natural): Measurement[T] = m ** p
-
-proc `**`*[T: FloatLike](m: Measurement[T], p: FloatLike): Measurement[T] =
-  result = procRes(pow(m.val, p), p * pow(m.val, (p - 1.0)), m)
-proc `^`*[T: FloatLike](m: Measurement[T], p: FloatLike): Measurement[T] = m ** p
 
+from measuremancer / math_utils import power
+## -> for static exponents
+proc `**`*[T: FloatLike](m: Measurement[T], p: static SomeInteger): auto =
+  # Get the resulting type
+  type U = typeof(power(m.val, p))
+  # and convert what is necessary.
+  result = procRes(U(power(m.val, p)), U(p.float * power(m.val, (p - 1))), m.to(U))
+proc `^`*[T: FloatLike](m: Measurement[T], p: static SomeInteger): auto =
+  ## NOTE: If you import `unchained`, this version is not actually used, as `unchained`
+  ## defines a macro `^` for static integer exponents, which simply rewrites
+  ## the AST to an infix (or trivial) node!
+  m ** p
+
+## -> for RT natural exponents
+proc `**`*[T: FloatLikeSupportsIntegerPowRT](m: Measurement[T], p: SomeInteger): Measurement[T] =
+  result = procRes(power(m.val, p), p.float * power(m.val, (p - 1)), m)
+proc `^`*[T: FloatLikeSupportsIntegerPowRT](m: Measurement[T], p: SomeInteger): Measurement[T] =
+  m ** p
+
+## -> for explicit float exponents
+proc `**`*[T: FloatLikeSupportsPow; U: FloatLikeNotInt](
+  m: Measurement[T], p: U): Measurement[T] =
+    result = procRes(pow(m.val, p), p * pow(m.val, (p - 1.0)), m)
+proc `^`*[T: FloatLikeSupportsPow; U: FloatLikeNotInt](
+  m: Measurement[T], p: U): Measurement[T] =
+    m ** p
 
 proc `**`*[T: FloatLike](a, b: Measurement[T]): Measurement[T] =
   let powV = pow(a.val, b.val)
   result = procRes(powV,
-                   [pow(a.val, b.val - 1.0),
-                    powV * log(a.val)],
+                   [pow(a.val, b.val - 1.0) * b.val,
+                    powV * ln(a.val)],
                    [a, b])
 proc `^`*[T: FloatLike](a, b: Measurement[T]): Measurement[T] = a ** b
 
@@ -328,11 +397,12 @@ proc `^`*[T: FloatLike](a, b: Measurement[T]): Measurement[T] = a ** b
 
 import std / macros
 template genCall(fn, letStmt, x, m, deriv: untyped): untyped =
-  proc `fn`*[T](m: Measurement[T]): Measurement[T] =
+  proc `fn`*[T](m: Measurement[T]): auto =
     ## letStmt contains the definition of x, `let x = m.val`
     letStmt
     ## `deriv` is handed as an expression containing `x`
-    result = procRes(fn(x), deriv, m)
+    type U = typeof(fn(x))
+    result = procRes(fn(x), deriv, m.to(U))
 
 macro defineSupportedFunctions(body: untyped): untyped =
   result = newStmtList()

measuremancer/math_utils.nim

@@ -0,0 +1,43 @@
+import std / macros
+
+# Taken from `unchained`. This won't be exported from `measuremancer` though,
+# rather we will `bind` in the scope in which we use this.
+macro power*(x: typed, num: static int): untyped =
+  ## general purpose power using `^` for integers, which works for any
+  ## type by rewriting to a product of `*`.
+  ##
+  ## For the special cases of -1, 0, 1 we simply rewrite to the correct
+  ## result. Negative powers are written as `1 / x^p`
+  if num == 0:
+    result = quote do:
+      1.0
+  elif num == 1:
+    result = x
+  elif num == -1:
+    ## Assume that the type supports addition by a float!
+    result = quote do:
+      1.0 / `x`
+  else:
+    result = nnkInfix.newTree(ident"*")
+
+    proc addInfix(n, x: NimNode, num: int) =
+      var it = n
+      if num > 0:
+        it.add nnkInfix.newTree(ident"*")
+        it[1].addInfix(x, num - 1)
+      while it.len < 3:
+        it.add x
+
+    result.addInfix(x, abs(num) - 2)
+
+    # invert if num is negative
+    if num < -1:
+      ## Assume that the type supports addition by a float!
+      result = quote do:
+        1.0 / (`result`)
+
+import std/math
+proc power*[T](x: T, num: SomeInteger): T =
+  if num > 0: result = x ^ num
+  elif num == 0: result = T(1)
+  else: result = T(1) / (x ^ abs(num))

tests/tmeasuremancer.nim

@@ -66,6 +66,50 @@ suite "Measurement & measurement":
     # √( 0.5² / 1.0² + 5.0² · 0.1² / 1.0⁴ ) = √ 0.5 = 0.707106781187
     check y / x =~= 5.0 ± 0.7071067811865476
 
+  test "Power (integer literal)":
+    let x = 2.0 ± 0.2
+    # This should be:
+    # √( (2 · 2.0)² · 0.2² ) = 0.8
+    check x ^ 2 == 4.0 ± 0.8
+    check x ** 2 == 4.0 ± 0.8
+
+  test "Power (static integer)":
+    let x = 2.0 ± 0.2
+    # This should be:
+    # √( (2 · 2.0)² · 0.2² ) = 0.8
+    const y = 2
+    check x ^ y == 4.0 ± 0.8
+    check x ** y == 4.0 ± 0.8
+
+  test "Power (integer)":
+    let x = 2.0 ± 0.2
+    # This should be:
+    # √( (2 · 2.0)² · 0.2² ) = 0.8
+    let y = 2
+    check x ^ y == 4.0 ± 0.8
+    check x ** y == 4.0 ± 0.8
+
+  test "Power (float)":
+    let x = 2.0 ± 0.2
+    # This should be:
+    # √( [ 2.5 · 2.0^{1.5} ]² · 0.2² ) = √ ( 50 · 0.² ) = sqrt(2)
+    check x ^ 2.5 == pow(2.0, 2.5) ± sqrt(2.0)
+    check x ** 2.5 == pow(2.0, 2.5) ± sqrt(2.0)
+
+  test "Negative power (integer)":
+    let x = 2.0 ± 0.2
+    # This should be:
+    # √( (-2 / x³)² · Δx² ) = √( (-2 / 2³)² · 0.2² ) = 0.05
+    check x ^ (-2) == 1.0 / (x^2)
+    check x ** -2 == 1.0 / (x^2)
+    check x ^ -2 == 0.250 ± 0.0500
+
+  test "Power of two measurements":
+    let x = 2.0 ± 0.2
+    let y = 5.0 ± 0.5
+    check x ^ y == 32.0 ± 19.467818870203708
+    check x ** y == 32.0 ± 19.467818870203708
+
 suite "Measurement and same symbol":
   # Simple correlations due to same symbol being used are correctly handled.
   # i.e. subtraction and division results in a measurement without error.
@@ -138,6 +182,13 @@ suite "Measurements of other types (e.g. unchained)":
     check( (k1 + k2).value.type is keV )
     check( (k1 + k2).error.type is keV )
 
+  test "Construction of `measurement` via proc":
+    let k1 = measurement(5.0.keV, 1.0.keV)
+    let k2 = measurement(2.5.keV, 1.5.keV)
+    check k1 + k2 =~= 7.5.keV ± 1.802775637731995.keV
+    check( (k1 + k2).value.type is keV )
+    check( (k1 + k2).error.type is keV )
+
   test "Addition of incompatible units fails":
     let k1 = 5.0.keV ± 1.0.keV
     let m = 1.0.kg ± 0.1.kg
@@ -152,3 +203,55 @@ suite "Measurements of other types (e.g. unchained)":
     check m * a =~= 9.81.kg•m•s⁻² ± 1.101072658819571.kg•m•s⁻²
     check( (m * a).value.type is kg•m•s⁻² )
     check( (m * a).error.type is kg•m•s⁻² )
+
+  # integer powers can be supported for units. Float however not! (floats representing
+  # integer powers need to be converted to int manually!)
+  test "Power (integer literal)":
+    let x = 2.0.m ± 0.2.m
+    # This should be:
+    # √( (2 · 2.0)² · 0.2² ) = 0.8
+    check x ^ 2 == 4.0.m² ± 0.8.m²
+    check x ** 2 == 4.0.m² ± 0.8.m²
+
+  test "Power (static integer)":
+    let x = 2.0.m ± 0.2.m
+    # This should be:
+    # √( (2 · 2.0)² · 0.2² ) = 0.8
+    const y = 2
+    check x ^ y == 4.0.m² ± 0.8.m²
+    check x ** y == 4.0.m² ± 0.8.m²
+
+  test "Negative power with units":
+    let x = 2.0.m ± 0.2.m
+    defUnit(Meter⁻²)
+    # This should be:
+    # √( (-2 / x³)² · Δx² ) = √( (-2 / 2³)² · 0.2² ) = 0.05
+    check x ^ -2 == 1.0 / (x^2)
+    ## NOTE: As `**` does not use `unchained's` rewrite `^`, but rather
+    ## the regular `**` proc, the resulting unit is equivalent, but not
+    ## the *same Nim unit* (as they are `distinct` and `Meter⁻²` is not
+    ## declared in the measuremancer context. As such we cannot check the
+    ## sanity of the types at CT.
+    # TODO: investigate why this prints as `0.0500` instead of `0.050`
+    check $(x ** (-2)) == "0.250 ± 0.0500 Meter⁻²"
+    check (x ** -2).to(Meter⁻²) == to(1.0 / (x^2), Meter⁻²)
+
+  test "Trigonometric functions only work with UnitLess":
+    block UnitLessCheck:
+      ## the division is just a simple way to get a `UnitLess` value
+      ## as on `unchained` HEAD (2022-11-16) `UnitLess` is not exported.
+      let x = (1.0.kg / 1.0.kg) ± (0.1.kg / 1.0.kg)
+      # This should be:
+      # √ (cos²(1.0) · 0.1²)
+      check typeof(x.value) is UnitLess
+      # application of `sin` returns a `float`! Valid units for `sin` are those convertible
+      # to `float` via converter defined in unchained
+      check typeof(sin(x)) is Measurement[float]
+      let exp = sin(x.value) ± sqrt(cos(x.value)^2 * (x.error)^2)
+      check typeof(exp) is Measurement[float]
+      check sin(x) =~= exp
+    block RadianCheck: # converts to UnitLess
+      let x = 1.0.Radian ± 0.1.Radian
+      # This should be:
+      # √ (cos²(1.0) · 0.1²)
+      check sin(x) =~= sin(1.0.Radian) ± sqrt(cos(1.0.Radian)^2 * 0.1.Radian^2) # 0.05403023058681398