PoC for a generalized conversion of unit system A to system B

generalized_unit_system_conversion.nim

#[
This code solves the unit conversion from system A to system B generically.
If natural units are involed, the remaining units (e.g. Energy in particle physics
natural units) are combined with the constants set to 1 that replace real units.

It solves the following system of linear equations which arises from:

`Π_i,α A_i^α = Π_i^β B_i^β`

where `A_i` is the i-th unit of unit system `A` and `α` the power needed to raise
each unit (and `B_i, β` same for system `B`). The product expresses the product of
the quantities in each system, `A` the source and `B` the target system.

The system of linear equations arises by inserting the definitions of the
quantities in system `A` by their counterpart in `B`.

For example for standard natural units in particle physics (reduced to 3 components):

A:
  1. energy via eV    = J   = kg•m²•s⁻²
  2. action via h_bar = J•s = kg•m²•s⁻¹
  3. velocity via c   = m•s⁻¹

yields the system:
```
       kg  m  s
eV    ( 1  2 -2 )   ( α_1 )   ( β_1 ) kg
h_bar ( 1  2 -1 ) · ( α_2 ) = ( β_2 ) m
c     ( 0  1  1 )   ( α_3 )   ( β_3 ) s

          ^A          ^x        ^b
```
and so

`A·x = b`

is solved for `x`, which yields the coefficients `α_i` required to perform
the conversion of a measurement `X` with value `V` units `U` from system `A` to system `B`
by

`V_B = V_A Π_i A_i^α_i`

In addition a check for the units is performed by verifying that the input unit
`U` shows up with powers `α_i` of the remaining units in the system.
]#



##################################################
# Definitions of a simple Matrix type
##################################################

type
  Matrix = object
    N: int
    M: int
    data: seq[float]
  Vector = seq[float]

proc initVector(N: int): Vector = Vector(newSeq[float](N))
proc toVector(x: openArray[float]): Vector = Vector(@x)

proc initMatrix(N, M: int): Matrix =
  result = Matrix(N: N, M: M, data: newSeq[float](N*M))

proc `[]`(m: Matrix, i, j: int): float = m.data[i * m.M + j]
proc `[]`(m: var Matrix, i, j: int): var float = m.data[i * m.M + j]
proc `[]=`(m: var Matrix, i, j: int, val: float) = m.data[i * m.M + j] = val
proc toMatrix(ar: seq[seq[float]]): Matrix =
  let N = ar.len
  let M = ar[0].len
  result = initMatrix(N, M)
  for i, row in ar:
    for j, val in row:
      result[i,j] = val

proc transpose(m: Matrix): Matrix =
  let N = m.N
  let M = m.M
  result = initMatrix(M, N)
  for i in 0 ..< N:
    for j in 0 ..< M:
      result[j,i] = m[i,j]

proc `$`(m: Matrix): string =
  let N = m.N
  let M = m.M
  for i in 0 ..< N:
    result.add "["
    for j in 0 ..< M:
      if j > 0:
        result.add " " & $m[i,j]
      else:
        result.add $m[i,j]
    result.add "]\n"

proc swap(m: var Matrix, row1, row2: int) =
  ## swap rows 1 by 2
  let M = m.M
  for i in 0 ..< M:
    swap(m[row1, i], m[row2, i])

## Instead of computing inverse, just solve linear system of equations
func gaussPartialScaled(a: Matrix; b: Vector): Vector =
  ## Solve linear system of equations described by
  ##
  ## `a · b = c`
  ##
  ## for `b` using gaussian elimination.
  ## Taken from rosetta code:
  ## https://rosettacode.org/wiki/Gaussian_elimination#Nim
  ## just to have something standalone for this code.
  doAssert a.N == b.len, "matrix and vector have incompatible dimensions"
  let N = a.N

  var m = initMatrix(N, N + 1)
  for i in 0 ..< N:
    for j in 0 ..< N:
      m[i,j] = a[i,j]
    m[i,N] = b[i]

  for k in 0..<N:
    var imax = 0
    var vmax = -1.0

    for i in k..<N:
      # Compute scale factor s = max abs in row.
      var s = -1.0
      for j in k..N:
        let e = abs(m[i,j])
        if e > s: s = e
      # Scale the abs used to pick the pivot.
      let val = abs(m[i,k]) / s
      if val > vmax:
        imax = i
        vmax = val

    if m[imax,k] == 0:
      raise newException(ValueError, "matrix is singular")

    swap m, imax, k

    for i in (k + 1)..<N:
      for j in (k + 1)..N:
        m[i,j] -= m[k,j] * m[i,k] / m[k,k]
      m[i,k] = 0

  result = initVector(N)
  for i in countdown(N - 1, 0):
    result[i] = m[i,N]
    for j in (i + 1)..<N:
      result[i] -= m[i,j] * result[j]
    result[i] /= m[i,i]

##################################################
# Definitions of a simple unit system
##################################################

import std / [tables, math, sets]
from std/strutils import removeSuffix
# define unit system
type
  Units = enum
    Energy
    Action
    Speed
    #SmallMass
    #SmallLength

    Mass
    Length
    Time

  ## Ideally we'd use a CountTable, but they remove entries with values of 0.
  ## We need to differentiate between 0 and non existing however.
  Unit = OrderedTable[Units, int]

  Measure = object
    val: float
    unit: Unit

  UnitSystem = OrderedSet[Units]

  ## How given units map to other units
  SystemConvert = proc(u: Units): Unit

  UnitError = object of CatchableError

proc initUnit(): Unit = initOrderedTable[Units, int]()

proc `$`(u: Unit): string =
  for k, v in u:
    result.add $k
    if v != 1:
      result.add "^" & $v
    result.add "·"
  result.removeSuffix("·")
proc `$`(m: Measure): string =
  result = $m.val & " " & $m.unit

proc toUnit(units: varargs[(Units, int)]): Unit =
  result = initOrderedTable[Units, int]()
  for (u, p) in units:
    result[u] = p

proc add(u: var Unit, unit: Units, power: int) =
  if unit notin u:
    u[unit] = 0
  u[unit] += power

proc toBase(u: Units): Unit =
  case u
  of Energy: toUnit([(Mass, 1), (Length, 2), (Time, -2)])
  of Action: toUnit([(Mass, 1), (Length, 2), (Time, -1)])
  of Speed:  toUnit([(Mass, 0), (Length, 1), (Time, -1)])
  #of SmallMass:   toUnit([(Mass, 1), (Length, 0), (Time, 0)])
  #of SmallLength: toUnit([(Mass, 0), (Length, 1), (Time, 0)])
  of Mass:   toUnit([(Mass, 1), (Length, 0), (Time,  0)])
  of Length: toUnit([(Mass, 0), (Length, 1), (Time,  0)])
  of Time:   toUnit([(Mass, 0), (Length, 0), (Time,  1)])

proc toValue(u: Units): float =
  case u
  of Energy: 1.602e-19
  of Action: 6.626e-34 / (2*PI)
  of Speed:  299792458.0
  #of SmallMass: 1e-3
  #of SmallLength: 1e-2
  of Mass:   1.0
  of Length: 1.0
  of Time:   1.0

proc isBase(u: Unit): bool =
  if u.len == 1:
    for v in values(u):
      if v == 1: return true

proc toMeasure(v: float, u: Unit): Measure = Measure(val: v, unit: u)

proc toMatrix(source, target: UnitSystem, toTarget: SystemConvert): Matrix = # this would be RT based of course
  let N = source.card
  result = initMatrix(N, N)
  var i = 0
  for x in source:
    let bu = x.toTarget()
    for j, base in target:
      let power = bu[base]
      result[i,j] = power.float
    inc i

proc toTarget(s: UnitSystem, unit: Unit): Vector =
  result = initVector(s.card)
  for i, t in s:
    if t in unit:
      result[i] = unit[t].float

proc getUnits(v: Vector, s: UnitSystem): Unit =
  ## Convers the given vector to a Unit given the unit system.
  result = initUnit()
  for i, t in s:
    if v[i] != 0:
      result.add t, v[i].int

proc checkUnits(u: Unit, m: Measure, desired: Unit) =
  ## Checks if the units of the `m` are compatible with `u`. This is
  ## done by checking if the powers of units `u` appearing in `m`
  ## are equivalent
  for unit, power in u:
    if unit in m.unit and power != m.unit[unit]:
      raise newException(UnitError, "Incompatible units of `" & $m & "` for " &
        "desired target " & $desired)

proc toSystem(m: Measure, source, target: UnitSystem, asUnit: Unit): Measure =
  ## This performs the actual unit converson from system `source` to `target`
  ## given a measurement `m` to desired units `asUnit`.
  # Note: `toBase` could become a RT based thing
  # transpose necessary, because I think the algorithm expects the indexing to
  # be inverted
  let sourceMat = toMatrix(source, target, toBase).transpose()
  echo sourceMat
  # generate vector based on desired target unit in `target`
  let v = target.toTarget(asUnit)
  echo "\tInput vector: ", v
  # solve the linear system
  let p = gaussPartialScaled(sourceMat, v)
  # compute the result

  # get units of result
  let resUnits = p.getUnits(source)
  echo "\tResultUnits: ", resUnits

  checkUnits(resUnits, m, asUnit)

  var res = m.val
  for i, t in source:
    res *= pow(t.toValue, p[i])
  result = toMeasure(res, asUnit)

block NaturalToSI:
  let u = toUnit([(Energy, 2)])
  let meas = toMeasure(1.0, u)

  # may be confusing, but effectively hbar & c are "units"
  # in the sense that they are required knowledge about how to
  # get back to SI
  const natural = [Energy, Action, Speed].toOrderedSet()
  const SI = [Mass, Length, Time].toOrderedSet()
  # (to make this runtime base, just replace the enum by
  # e.g. strings & a UnitSystem by a HashSet[Units = alias string]
  # or similar

  #echo toSystem(meas, natural, SI, toUnit([(Mass, 1)]))
  echo toSystem(meas, natural, SI, toUnit([(Mass, 1), (Time, -1)]))
  #echo toSystem(meas, natural, SI, toUnit([(Mass, 1), (Time, 1)]))
  #echo toSystem(meas, natural, SI, toUnit([(Length, 1), (Time, -1)]))

when false:
#block CGStoSI:
  let u = toUnit([(SmallLength, 2), (SmallMass, 1), (Time, -2)])
  let meas = toMeasure(1.0, u)

  const CGS = [SmallMass, SmallLength, Time].toOrderedSet()
  const SI  = [Mass, Length, Time].toOrderedSet()

  #echo toSystem(meas, natural, SI, toUnit([(Mass, 1)]))
  echo toSystem(meas, CGS, SI, toUnit([(Mass, 1), (Length, 2), (Time, -2)]))
  #echo toSystem(meas, natural, SI, toUnit([(Mass, 1), (Time, 1)]))
  #echo toSystem(meas, natural, SI, toUnit([(Length, 1), (Time, -1)]))

when false:
  # to cross check 1.g•cm²•s⁻² conversion back to SI units
  import unchained
  defUnit(g•cm²•s⁻²)
  echo 1.g•cm²•s⁻².to(J)