Compile-time precision set, Pi->PI, inverse-var weighed average & red…

…ucer (#20)

* Allow compile-time setting of printing precision (number of places for
the error; `cligen/strUt.fmtUncertain` rounds value to same decimal).
810bdfb changed the default here from
2 to 3, but some users might even want 1.

* `Pi` -> `PI` because that is how `std/math.nim` spells it.  This lets
users do `nim c --hint:Name:on --styleCheck=error --styleCheck=usages`.

* Add inverse-var weighted averaging operator & varargs/openArray reducer
with a URI justifying the weights.

* Update README per new compile-time override of uncertain precision.

* Fix typo in spelling of `import measuremancer` in unchained example.

* Improve doc comment of `mean` & add one for `+/` adapting extant text.

README.org

@@ -90,7 +90,7 @@ logic to types other than =float=!
 Feel free to perform error propagation on =unchained= types for
 example:
 #+begin_src nim
-import unchained, meausuremancer
+import unchained, measuremancer
 
 let m = 1.0.kg ± 0.1.kg
 let a = 9.81.m•s⁻² ± 0.05.m•s⁻²
@@ -127,5 +127,5 @@ avoid having to set this setting manually every time.
 If compiled with `-d:useCligen`, uncertain number pair formatting is the
 default [`cligen/strUt`](https://github.com/c-blake/cligen/blob/master/cligen/strUt.nim).`fmtUncertain`
 That uses the uncertainty to limit precision of both value & uncertainty to the
-same decimal place (2 decimals of the uncertainty by default, as per one common
-convention).
+same decimal place (3 decimals of the uncertainty by default, as per one common
+convention, but you can compile with `-d:mmErrPrec=N` to use `N` places).

measuremancer.nim

@@ -283,7 +283,8 @@ else:
     when not (T is float):
       result.add " " & $T
 
-proc `$`*[T: FloatLike](m: Measurement[T]): string = pretty(m, 3)
+const mmErrPrec* {.intdefine.} = 3  # Digits of err; -d:mmErrPrec=N to edit
+proc `$`*[T: FloatLike](m: Measurement[T]): string = pretty(m, mmErrPrec)
 
 template print*(arg: untyped): untyped =
   echo astToStr(arg), ": ", $arg
@@ -519,24 +520,24 @@ defineSupportedFunctions:
   sec         ->  sec(x) * tan(x)
   csc         -> -csc(x) * cot(x)
   cot         -> -(1.0 + (cot(x)^2))
-  #sind        ->  Pi / 180.0 * cosd(x)
-  #cosd        -> -Pi / 180.0 * sind(x)
-  #tand        ->  Pi / 180.0 * (1.0 + (tand(x)^2))
-  #secd        ->  Pi / 180.0 * secd(x) * tand(x)
-  #cscd        -> -Pi / 180.0 * cscd(x) * cotd(x)
-  #cotd        -> -Pi / 180.0 * (1.0 + (cotd(x)^2))
+  #sind        ->  PI / 180.0 * cosd(x)
+  #cosd        -> -PI / 180.0 * sind(x)
+  #tand        ->  PI / 180.0 * (1.0 + (tand(x)^2))
+  #secd        ->  PI / 180.0 * secd(x) * tand(x)
+  #cscd        -> -PI / 180.0 * cscd(x) * cotd(x)
+  #cotd        -> -PI / 180.0 * (1.0 + (cotd(x)^2))
   arcsin      ->  1.0 / sqrt(1.0 - (x^2))
   arccos      -> -1.0 / sqrt(1.0 - (x^2))
   arctan      ->  1.0 / (1.0 + (x^2))
   arcsec      ->  1.0 / abs(x) / sqrt(x^2 - 1.0)
   arccsc      -> -1.0 / abs(x) / sqrt(x^2 - 1.0)
   arccot      -> -1.0 / (1.0 + (x^2))
-  #arcsind     ->  180.0 / Pi / sqrt(1.0 - (x^2))
-  #arccosd     -> -180.0 / Pi / sqrt(1.0 - (x^2))
-  #arctand     ->  180.0 / Pi / (1.0 + (x^2))
-  #arcsecd     ->  180.0 / Pi / abs(x) / sqrt(x^2 - 1.0)
-  #arccscd     -> -180.0 / Pi / abs(x) / sqrt(x^2 - 1.0)
-  #arccotd     -> -180.0 / Pi / (1.0 + (x^2))
+  #arcsind     ->  180.0 / PI / sqrt(1.0 - (x^2))
+  #arccosd     -> -180.0 / PI / sqrt(1.0 - (x^2))
+  #arctand     ->  180.0 / PI / (1.0 + (x^2))
+  #arcsecd     ->  180.0 / PI / abs(x) / sqrt(x^2 - 1.0)
+  #arccscd     -> -180.0 / PI / abs(x) / sqrt(x^2 - 1.0)
+  #arccotd     -> -180.0 / PI / (1.0 + (x^2))
   sinh        ->  cosh(x)
   cosh        ->  sinh(x)
   tanh        ->  sech(x)^2
@@ -549,13 +550,13 @@ defineSupportedFunctions:
   arcsech     -> -1.0 / x / sqrt(1.0 - (x^2))
   arccsch     -> -1.0 / abs(x) / sqrt(1.0 + (x^2))
   arccoth     ->  1.0 / (1.0 - (x^2))
-  deg2rad     ->  Pi / 180.0
-  rad2deg     ->  180.0 / Pi
-  erf         ->  2.0 * exp(-x*x) / sqrt(Pi)
-  erfinv      ->  0.5 * sqrt(Pi) * exp(erfinv(x) * erfinv(x))
-  erfc        -> -2.0 * exp(-x*x) / sqrt(Pi)
-  erfcinv     -> -0.5 * sqrt(Pi) * exp(erfcinv(x) * erfcinv(x))
-  erfi        ->  2.0 * exp(x*x) / sqrt(Pi)
+  deg2rad     ->  PI / 180.0
+  rad2deg     ->  180.0 / PI
+  erf         ->  2.0 * exp(-x*x) / sqrt(PI)
+  erfinv      ->  0.5 * sqrt(PI) * exp(erfinv(x) * erfinv(x))
+  erfc        -> -2.0 * exp(-x*x) / sqrt(PI)
+  erfcinv     -> -0.5 * sqrt(PI) * exp(erfcinv(x) * erfcinv(x))
+  erfi        ->  2.0 * exp(x*x) / sqrt(PI)
   #gamma       ->  digamma(x) * gamma(x)
   #lgamma      ->  digamma(x)
   #digamma     ->  trigamma(x)
@@ -569,7 +570,7 @@ defineSupportedFunctions:
   #besselj1    ->  (besselj0(x) - besselj(2.0, x)) / 2.0
   #bessely0    -> -bessely1(x)
   #bessely1    ->  (bessely0(x) - bessely(2.0, x)) / 2.0
-  #erfcx       ->  (2.0 * x * erfcx(x) - 2.0 / sqrt(Pi))
+  #erfcx       ->  (2.0 * x * erfcx(x) - 2.0 / sqrt(PI))
   #dawson      ->  (1.0 - 2.0 * x * dawson(x))
 
 func signbit*[T: FloatLike](m: Measurement[T]): bool = m.val.signbit
@@ -596,6 +597,21 @@ template comp(fn: untyped): untyped =
 comp(`<`)
 comp(`<=`)
 
+proc `+/`*[F](x, y: Measurement[F]): Measurement[F] =
+  ## en.wikipedia.org/wiki/Weighted_arithmetic_mean#Variance-defined_weights
+  ## suggests inverse-variance weights via max.likelihood over normal distros
+  ## theory which also makes sense from "variance linearity" perspective.
+  ## This binary operator does this for a pair of measurements.
+  let (wx, wy) = (x.error^(-2.0), y.error^(-2.0))
+  (wx*x + wy*y)/(wx + wy)   # ws are scalars but x,y are Measurements
+
+proc mean*[F](x: varargs[Measurement[F]]): Measurement[F] =
+  ## Average however many uncertain values in `x` into an overall estimate using
+  ## the binary `+/` operator.  `x` can also be any `openArray[Measurement[F]]`.
+  if x.len > 0:
+    result = x[0]
+    for i in 1 ..< x.len: result = result +/ x[i]
+
 when isMainModule:
   proc foo[T](x: T, μ, σ: float): T =
     let val = 2 * σ * σ