Julia, Python, F#, C#, C++で単位付き数値を扱う

技術者なら「プログラミングで単位の扱いの間違いはコンパイルエラーにして欲しい」と思った事があるでしょう。
mm と inch でごっちゃになったとか、次元(ディメンジョン)が不合理な演算をしたとか、そういう頭の痛いバグに悩まされることありますよね。
言語によっては実行時またコンパイル時に単位チェックできるパッケージもあります。
今回は、Julia, Python, F#, C#, C++ について見てみましょう。
# Julia
## Unitful
Unitfulパッケージのインストール
`julia`でREPLに入り`]`でパッケージ管理に入る。
```julia
(@v1.8) pkg> add Unitful
```
使ってみる。バックスペースでREPLに戻る。
```julia
julia> using Unitful

julia> 1.0u"m" == 1000.0u"mm"
true

julia> 0.9u"m" == 1000.0u"mm"
false
```
μ は、`\mu`と打鍵してから TABキーを押すと μ に変換されます。
```julia
julia> 1.0u"mm" + 1.0u"μm"
0.001001 m

julia> 1.0u"mm" + 1.0u"nm"
0.001000001 m

julia> len = 2.0u"m"
2.0 m

julia> t = 1.0u"s"
1.0 s

julia> mass = 3.0u"kg"
3.0 kg

julia> mass * (len/t)^2
12.0 kg m² s⁻²
```

# Python
units というパッケージがある様ですが、この`units`はなんと言うかイケてない、書きづらい。
単位の変換の方法もドキュメントが見つからず、よく分からない。開発はやめているのかも。
なので、`astropy`の方の`units`を使う。
## astropy.units
`astropy`は天文学関連のパッケージらしいです。単位関連だけ使います。
```shell_session
pip install astropy
```
REPLで
```python
>>> from astropy import units as u
>>> 1.0*u.m == 1000.0*u.mm
True
>>> 0.9*u.m == 1000.0*u.mm
False
```
Pythonでは μ は u です。
```python
>>> 1*u.mm + 1*u.um
<Quantity 1.001 mm>
>>> 1*u.mm + 1*u.nm
<Quantity 1.000001 mm>
>>> len = 2.0 * u.m
>>> mass = 3.0 * u.kg
>>> time = 1.0 * u.s
>>> mass * (len/time)**2
<Quantity 12. kg m2 / s2>
```

# F#

F#では標準で単位がサポートされています。
SI単位系とかも用意されています。
`dotnet fsi`でREPLに入ります。
```fsharp
❯ dotnet fsi

Microsoft (R) F# インタラクティブ バージョン F# 7.0 のための 12.5.0.0
Copyright (C) Microsoft Corporation. All rights reserved.

ヘルプを表示するには次を入力してください: #help;;

> open Microsoft.FSharp.Data.UnitSystems.SI.UnitSymbols;; // SI単位系を使わないなら要らない
> 1.0<kg> + 2.0<kg>;;
val it: float<kg> = 3.0
```
が、グラム`g`とかは定義されていません。SI単位系に含まれませんからね。
```fsharp
> 2.0<g>;;
  2.0<g>;;
  ----^
stdin(4,11): error FS0039: 型 'g' が定義されていません。
```
無いものは自分で定義できます。グラムを定義
```fsharp
> [<Measure>] type g;;
[<Measure>]
type g
```
当然 g から kg への暗黙の変換などは無く
```fsharp
> 1.0<kg> + 2.0<g>;;
  1.0<kg> + 2.0<g>;;
  ----------^^^^^^
stdin(16,11): error FS0001: 測定単位 'g' は測定単位 'kg' と一致しません

> 1.0<kg> + 2.0<g> * 1000.0;;
  1.0<kg> + 2.0<g> * 1000.0;;
  -------------------^^^^^^
/Users/yamana/fsharp/work/stdin(4,20): error FS0001: 測定単位 'g' は測定単位 'kg' と一致しません
```
変換する関数を作るか、変換レートを定義してやらないといけません。
```fsharp
> let g2kg x = x / 1000.0<g/kg>;; // 変換関数
val g2kg: x: float<'u> -> float<'u kg/g>

> 1.0<kg> + g2kg 2.0<g>;;
val it: float<kg> = 1.002

> let g_per_kg = 1.0 / 1000.0<g/kg>;; // 変換レート
val g_per_kg: float<kg/g> = 0.001

> 1.0<kg> + g_per_kg * 2.0<g>;;
val it: float<kg> = 1.002
```
単位付き値の2乗のやり方が分からず`v * v`なんてしています。`v^2`とかは駄目でした。誰かやり方教えて下さい。
```fsharp
> let len = 2.0<m>;;
val len: float<m> = 2.0

> let time = 1.0
;;
val time: float
 = 1.0

> let mass = 3.0<kg>;;
val mass: float<kg> = 3.0

> let v = len / time;;   
val v: float<m/s> = 2.0

> let energy = mass * (v*v);;
val energy: float<kg m ^ 2/s ^ 2> = 12.0
```
定義済みの単位から単位を定義する事もできます。Galを定義するなら
$1 Gal = 1 cm/s^2$
```fsharp
> [<Measure>] type cm;; // cmはSI単位系に含まれない
[<Measure>]
type cm

> [<Measure>] type Gal = cm/s^2;;
[<Measure>]
type Gal = cm/s ^ 2

> let x:float<Gal> = 12.3<cm>/2.0
/3.0
;;
val x: float<Gal> = 2.05
```
F#は金融系でも使われるらしいです。
通貨単位の間違いが減るので嬉しいんですかね。F#が使われるのは関数型だからというのが主な理由だとは思いますが。
```fsharp
> [<Measure>] type USD;;
[<Measure>]
type USD

> [<Measure>] type JPY;;
[<Measure>]
type JPY

> let rateUsdToJpy = 150M<JPY/USD>;; // 円安!!
val rateUsdToJpy: decimal<JPY/USD> = 150M

> let rateJpyToUsd = 1.0M / rateUsdToJpy;;
val rateJpyToUsd: decimal<USD/JPY> = 0.0066666666666666666666666667M

> let dollars: decimal<USD> = 123.45M<USD>;;
val dollars: decimal<USD> = 123.45M

> let yen:decimal<JPY> = dollars * rateUsdToJpy;;
val yen: decimal<JPY> = 18517.50M

> let dollars2:decimal<USD> = yen * rateJpyToUsd;;
val dollars2: decimal<USD> = 123.45000000000000000000000062M

> let dollars2:decimal<USD> = yen * rateUsdToJpy;; //レートを間違えたら

  let dollars2:decimal<USD> = yen * rateUsdToJpy;;
  ----------------------------------^^^^^^^^^^^^

/Users/yamana/fsharp/work/stdin(34,35): error FS0001: 測定単位 'USD' は測定単位 'JPY ^ 2/USD' と一致しません
```
小数点以下が長すぎるので`Decimal.Round`を使いたいですよね。`Decimal.Round`が単位付きの`decimal`型を受け付けないので、一時的に単位を外して`Round`を呼ばないといけません。
fsiでは長い関数を書けなかったので、以下はファイルに書いて`dotnet run`で実行します。
```fsharp:Program.fs
open System

[<Measure>]
type USD

[<Measure>]
type JPY

let rateUsdToJpy = 150.0M<JPY / USD>
let rateJpyToUsd = 1.0M / rateUsdToJpy

// 単位を取り除く関数
let inline stripUnit (value: decimal<'u>) : decimal = value |> decimal

// 単位を再適用する関数
let inline applyUnit (value: decimal) : decimal<'u> =
    value |> LanguagePrimitives.DecimalWithMeasure

// 単位を取り除いて丸める
let roundAmount (x:decimal<'u>) (deci) : decimal<'u> =
    x
    |> stripUnit
    |> fun x -> Decimal.Round(x, deci, MidpointRounding.ToEven) // 銀行丸め
    |> applyUnit


[<EntryPoint>]
let main argv =
    let dollars: decimal<USD> = 123.45M<USD>
    let yen: decimal<JPY> = dollars * rateUsdToJpy // OK
    //let yen:decimal<JPY> = dollars * rateJpyToUsd // コンパイルエラー
    let dollars: decimal<USD> = yen * rateJpyToUsd // OK
    printfn "yen=%O, dollars=%O" yen dollars
    //printf "USD+JPY=%O" (dollars + yen) // コンパイルエラー
    printfn "rate USD/JPY = %O, JPY/USD = %O" (dollars / yen) (yen / dollars)
    printfn "rate USD/JPY = %O, JPY/USD = %O" rateJpyToUsd rateUsdToJpy
    let amountInDollars: decimal<USD> = 123.4567M<USD>
    let yen2:decimal<JPY> = roundAmount yen 0
    printfn "%O yen" yen2
    let dollars2:decimal<USD> = roundAmount dollars 2
    printfn "%O dollers" dollars2
    0
```

```shell-session
❯ dotnet run
yen=18517.500, dollars=123.45000000000000000000000062
rate USD/JPY = 0.0066666666666666666666666667, JPY/USD = 149.99999999999999999999999925
rate USD/JPY = 0.0066666666666666666666666667, JPY/USD = 150.0
18518 yen
123.45 dollers
```

# C#
## UnitsNet
適当なプロジェクトを作ってから、プロジェクトに`UnitsNet`を追加します。
shellで
```shell_session
mkdir units //名前は何でも良い
cd units
dotnet new console
dotnet add package UnitsNet
```

```csharp:Program.cs
using u = UnitsNet;

var len = u.Length.FromMeters(2.0);
var time = u.Duration.FromSeconds(1.0);
var mass = u.Mass.FromKilograms(3.0);
var mass1 = u.Mass.FromGrams(2.0);
Console.WriteLine("{0:E}", mass + mass1);
var v = len / time;
Console.WriteLine("{0:E}", v);
var energy = mass * (v * v);  // 2乗の仕方が分からないのでこうしています
Console.WriteLine("{0:E}", energy);
```
C#も単位付き値の2乗のやり方が分からないんです。まあ`(v * v)`でも良いんですが。
shellで
```shell_session
dotnet run
3.002000E+000 kg
2.000000E+000 m/s
1.200000E+001 J
```
エネルギーの単位が`J`と表示されるのは良いですね。

# C++
## boost::units
boostです。boostのインストール方法は省きます。REPLは基本無いのでファイルに書きます。
```cpp:units_work.cpp
#include <iostream>
#include <boost/units/systems/si.hpp>
#include <boost/units/systems/cgs.hpp>
#include <boost/units/pow.hpp>
#include <boost/units/io.hpp> // std::ostream への出力に必要

using namespace boost::units;

int main()
{
  // boost::units::quantity に単位型と値型を渡す事で単位付きの値型となる。
  // ここでは、単位型としてboost::units::si::length、
  // 値型としてdoubleを使用する。
  quantity<si::length, double> len = 2.0 * si::meter;
  auto len1 = 2.0 * si::meter; // 型は明示しなくても良い

  auto time = 3.0 * si::second;

  auto weight = 4.0 * si::kilogram;
  auto weight1 = 5.0 * cgs::gram;

  auto wkg = 1.0 * si::kilogram;
  auto wg = 2.0 * cgs::gram;
  auto weight2 = wkg + static_cast<quantity<si::mass>>(wg); // 明にキャストする必要あり

  std::cout << len << std::endl
            << len1 << std::endl
            << time << std::endl
            << weight << std::endl
            << weight1 << std::endl
            << weight2 << std::endl
            << len1 / time << std::endl
            << weight * pow<2>(len / time) << std::endl;

  return 0;
}
```
shellで
```shell_session
clang++ -std=c++20 -stdlib=libc++ -I/opt/homebrew/include/ -Wall -O2 units_work.cpp -o units_work.exe && ./units_work.exe
2 m
2 m
3 s
4 kg
5 g
1.002 kg
0.666667 m s^-1
1.77778 m^2 kg s^-2
```

# 言語ごとの比較

## 単位の乗除
メートルを秒で割って速度の単位になるか
### Julia
```julia
julia> u"m" / u"s"
m s⁻¹
```
### Python
```python
>>> u.m / u.s
Unit("m / s")
```
JuliaやPythonは単位がインスタンスになる様です。
### F#
F#では単位はインスタンスでなく、あくまで「型(属性?)」なので、単位だけの乗除算はできなさそうです。ただ、単位`m`と`s`が定義されているなら`<m/s>`と書くことができます。
```fsharp
> 1<m/s>;;
val it: int<m/s> = 1
```

### C#
```csharp
var len = u.Length.FromMeters(2.0);
var time = u.Duration.FromSeconds(1.0);
Console.WriteLine((len / time).GetType());  // UnitsNet.Speed
```
C#も「単位」というインスタンスは無いので単位だけの乗除算はできなそうです。型表示させると`Speed`と表示されます。
### C++
```cpp
std::cout << si::meter / si::second << std::endl;  // m s^-1
```
単位を`m s^-1`と表示しますね。
## 次元チェック
- 次元が指定された変数への代入で次元チェックできるかどうか
- 次元が不合理な計算をエラーにできるかどうか
### Julia
変数の次元を`::typeof(1.0u"kg*(m/s)^2")`の様に定義しています。直接的な単位型の記述が分からなかったので、適当な単位付き数リテラルから`typeof`で型を引き出す様にしました。
この方法で変数毎に毎回`typeof`で記述すると、型のパースとtypeofのオーバーヘッドが馬鹿にならないので、同じ次元の単位は上の方で`speed_t = typeof(1.0u"m/s")`の様に単位型の変数を用意しておいて、`v::speed_t = len / t`の様にすべきと思います。
```julia:unitful_work.jl
using Unitful

len = 2.0u"m"
t = 1.0u"s"
v::typeof(1.0u"m/s") = len / t
println(v)

mass = 3.0u"kg"
energy::typeof(1.0u"kg*(m/s)^2") = mass * (len/t)^2
println(energy)

x3::typeof(1.0u"g*(m/s)^3") = mass * (len/t)^3
println(x3)

# 下の3行はディメンジョンが違うのでエラー
#v2::typeof(1.0u"m") = len / t
#energy3::typeof(1.0u"g*(m/s)^2") = mass * (len/t)^3
#mass + x3

# これはディメンジョン違いの加算なのでエラー
#println(len + t)
```
次元違いの代入も加減算もしっかりチェックされてエラーになります。
JuliaのREPLで
```julia
julia> include("unitful_work.jl")
2.0 m s⁻¹
12.0 kg m² s⁻²
24000.0 g m³ s⁻³
```
良いですね。
### Python
変数の次元を型ヒント`:u.Unit("kg m2 / s2")`の様に指定します。
```python
from astropy import units as u

len = 2.0 * u.m
mass = 3.0 * u.kg
time = 1.0 * u.s

energy:u.Unit("kg m2 / s2") = mass * (len/time)**2
print(energy)

x2:u.Unit("kg m2 / s2") = mass * (len/time)**3
print(x2)

#energy + x2  # 次元違いの加算はエラー

@u.quantity_input(myarg=u.deg)
def myfunction(myarg):
    return myarg.unit

myfunction(100*u.arcsec)
#myfunction(2*u.m)
```
Pythonは型ヒントのチェックをしないので、mypy を使ってチェックするらしいですが、mypy の使い方がよく分かりません。astropyの stubs が無いとか、そんなエラーになります。評価は保留。
ディメンジョン違いの加算はエラーになります。これは良いですね。
関数の引数に関してはデコレータ`quantity_input`で指定ができて、間違えた単位の値を渡すとエラーになります。
shellで
```shell_session
python units_work.py
12.0 kg m2 / s2
24.0 kg m3 / s3
```
二つめが`kg m3 / s3`と表示はちゃんとしているので、mypy が使えるとうまくエラーに引っ掛かるのかもしれません。
### F#
```fsharp
> open Microsoft.FSharp.Data.UnitSystems.SI.UnitSymbols;;

> let len = 2.0<m>;;
val len: float<m> = 2.0

> let time = 1.0
;;
val time: float
 = 1.0

> let mass = 3.0<kg>;;
val mass: float<kg> = 3.0

> let v:float<m/s> = len / time;;
val v: float<m/s> = 2.0

> let v:float<m/s^2> = len / time;;

  let v:float<m/s^2> = len / time;;
  ---------------------------^^^^

error FS0001: 測定単位 'm/s ^ 2' は測定単位 'm/s' と一致しません

> let energy = mass * (v*v);;
val energy: float<kg m ^ 2/s ^ 2> = 12.0

> let energy:float<J> = mass * (v*v);;
val energy: float<J> = 12.0

> let energy:float<J> = mass * (v*v*v);;

  let energy:float<J> = mass * (v*v*v);;
  ----------------------------------^

error FS0001: 測定単位 'J' は測定単位 'kg m ^ 3/s ^ 3' と一致しません

> 2.0<J> + 3.0<kg m^2 / s^2>;;
val it: float<J> = 5.0

> 2.0<J> + 3.0<kg m^2 / s^3>;;

  2.0<J> + 3.0<kg m^2 / s^3>;;
  ---------^^^^^^^^^^^^^^^^^

error FS0001: 測定単位 'kg m ^ 2/s ^ 3' は測定単位 'J' と一致しません
```
良いですね。

### C#
```csharp
u.Length x = u.Duration.FromSeconds(1.0);  //エラーになる
var len = u.Length.FromMeters(2.0);
var time = u.Duration.FromSeconds(1.0);
var x2 = len + time;  //エラーになる
```
良いです。
### C++

```cpp:units_w2.cpp
#include <iostream>
#include <boost/units/pow.hpp>
#include <boost/units/io.hpp>   // std::ostream への出力に必要
#include <boost/units/systems/si.hpp>

using namespace boost::units;

int main() {
  auto len = 2.0*si::meters;
  auto mass = 3.0*si::kilograms;
  auto time = 1.0*si::seconds;
  quantity<si::energy> E = mass * pow<2>(len / time);
  std::cout << "len=" << len << ", E=" << E << std::endl;

  //これは次元違いでコンパイルエラー
  //quantity<si::energy> E = mass * pow<3>(len / time);

  //これは次元違いの加算なのでコンパイルエラー
  //当然、乗除算ならOK
  //std::cout << len + time << std::endl;

  return 0;
}
```
次元違いの代入も加減算もコンパイルエラーになります。
shellで
```shell_session
% clang++ -std=c++20 -stdlib=libc++ -I/opt/homebrew/include/ -Wall -O2 units_w2.cpp -o units_w2.exe && ./units_w2.exe
len=2 m, E=12 m^2 kg s^-2
```

# まとめ
単位付き数を扱う事に限って言えば F# が良いのですが、単位付き数を使いたいのは理系・工学系の人でしょうから、総合的には Julia が一番良いのかな。
実行時のコストが気になるなら F# がゼロコストらしいです。Juliaがどのくらいのオーバーヘッドがあるのか暇があれば調べたいと思います。