Elixir練習帳: .npyファイルの中を覗く

これは Elixir Advent Calendar 2020の5日目の記事です。

4日目の昨日は @myasu さんの"PLCをElixirでコントロールする"お話しでした。@myasu さん、お疲れさまでした。

小生は組み込み系の人なので聞きなれた話しでしたが、基幹系/WEB系の方にとっては異星人に見えたかも知れませんね。ですが、抽象度の高い関数型言語の筈なのに低レベルにも手が届く、そんなところが Elixirの面白いところであり、また大きな可能性を秘めているところかなと思います。

さて、今日も低レベルな話になりますが、暫しお付き合い願います。

#1.練習課題 - numpyのデータ保存ファイル.npyの中を覗く
先日、ぼっちProject: "Tflite YOLO V3 in Nerves"のデバッグ支援のために、Tfliteモデルの入出力tensorをnumpyデータ(.npy)として保存する機能(C++)を用意しました。これで、開発のお手本としたPython版 Tflite YOLO V3と動作を比較し易くはなったのですが‥‥.npyの中を覗くためにいちいちPython&numpyを起動するのは億劫だなぁと(^^;)我儘

という訳で、頭の老化防止かつElixirプログラミングの練習として、.npyの中を覗くElixirモジュールを書いてみることにしました。お手軽に.npyの中を覗きたいだけなので、凝ったことはせずに次のような簡単な仕様としましょう[*1]。

[設計仕様]
>1. モジュール名は "Npy"とする
2. numpy dtypeが '<f4'または '<i4'の .npyファイルが扱えれば良い
3. .npyファイル形式を解釈して、Npyデータ構造(適当に決めた内部表現)として load出来ること。また逆に、Npyデータ構造を.npyファイルとして saveできること。
4. loadした numpyデータの dtypeと shapeを見れること
5. loadした numpyデータを、Elixirのリストに変換できること。また、逆にElixirのリストをnumpyデータに変換できること

[*1]設計方針は、所謂KISS の原則「Keep it simple stupid」、YANGI「You ain't gonna need it」に従うことにします(^_-)

#2. .npyファイル形式
.npyファイルのフォーマットの詳細は参考文献[1]に公開されています。フォーマットには、今日現在 3つのバージョンがあるようですが、最もベーシックな Version 1.0 のみを対象にします。

参考文献[1]から、設計に必要な情報だけを拾ってまとめると以下の通りです。

####[.npyフォーマット・バージョン 1.0]

|offset |BYTE長 |意味|
|:------|------:|:----|
|+00    |6      |"\93NUMPY" - マジック・ワード(固定)|
|+06    |1      |"\x01" - フォーマット・バージョンの major番号 [unsigned char]|
|+07    |1      |"\x00" - フォーマット・バージョンの minor番号 [unsigned char]|
|+08    |2      |LEN - 後続のヘッダー部のBYTE長 [little endian/unsigned short]|
|+10    |LEN ※1|ヘッダー部 ["\n"終端文字列] |
|+10+LEN|*      |データ部|

※1) ファイル先頭からヘッダー部末尾の"\n"までの BYTE長(10+LEN)は 64でアライメントされています。BYTE長がアライメントに合うように、ヘッダー部の実データと"\n"終端の間に、空白文字"\x20"がパディングされます。LENは、そのパディングと\n"終端を含む長さです。

####ⅰ. ヘッダー部の詳細
ヘッダー部には、ファイルに格納された numpy配列の諸元が ASCII文字列(Pythonのdictionay表記)で記述されています。

```python
"{'descr': '<f4', 'fortran_order': False, 'shape': (1, 10647, 4), }"
```
<dl>
<dt>descr</dt>
  <dd>配列要素のデータ型(numpy.dtype)　例)'&lt;f4' : little endian 4byte float </dd>
<dt>fortran_order</dt>
<dd>配列要素の格納順序<br>
例) 配列 M[2][3]の場合↓<br>
　　fortran_order = True  -> M[1][1], M[2][1], M[1][2], M[2][2], M[1][3], M[2][3]<br>
　　fortran_order = False -> M[1][1], M[1][2], M[1][3], M[2][1], M[2][2], M[2][3]
<dt>shape</dt>
  <dd>配列の次元数・形状　例) (1,10647,4) : 1 x 10647 x 4 のtensor</dd>
</dl>

####ⅱ. データ部の詳細
ヘッダ部の情報に基づいて、配列要素がバイト列として隙間なく格納されています。

####ⅲ. .npyファイルの例
![npy.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/14158/5c77c14e-0f56-41a7-3be2-16a054c630bb.jpeg)

#3.設計&実装
KISSの原則に則って、loadしたデータを保持する Npyデータ構造は、.npyファイルの中身をほぼそのまま持つことにします。ヘッダ部の文字列(配列の諸元)は何かと扱い難いのでパースして、下記の Elixirデータ型に変換しますが、データ部のBYTE列はそのままバイナリとして持ちます。load時に Elixirのリストに変換するといったようなことはしません。

- 'descr' → String型
- 'fortran_order → boolean(Atom型)
- 'shape' → List型

Npyデータ構造の定義は下記の通りです。

```elixir
defmodule Npy do
  alias __MODULE__

  # npy data structure
  defstruct descr: "", fortran_order: false, shape: [], data: <<>>
```
◆Tips1: defstructで構造体を定義したら、`alias __MODULE__`も宣言しておこう◆

---
それでは、このNpyデータ構造を中心に据えて、上の仕様で挙げた４つの機能を実装していきましょう。
>(1) load("xxx.npy") ----------------- .npyファイルの読み込み 
(2) save("yyy.npy", %Npy{}) ------- .npyファイルへの保存
(3) to_list(%Npy{}) ----------------- Npyデータ構造からネストしたリストを生成
(4) from_list([[1,2],[3,4]], "<i4") ---- ネストしたリストからNpyデータ構造を生成

####(1)Npy.load関数
load/1関数は、引数fnameで指定した .npyファイルをオープンし、その読み込みデータから Npyデータ構造を生成して返します。

ファイルのオープンには、第３引数に関数func/1を取る File.open/3を使用しました(2行目)。この関数を使用すると、Pythonの`with open(fname) as f:`に似たマナーでファイルを扱うことができます。func/1にはオープンしたファイルが渡され実行されます。そして正常終了、エラー終了に係わらず、func/1が終了すると自動的にファイルがcloseされます。

3行目では、バイナリのパターンマッチを用いて、ファイルの先頭に.npyのマジック・ワードが有るかどうかをチェックしています。また同時に、フォーマットのバージョンmajor,minorの読み取りを行っています。この行の様に with関数を併用すると、パターンマッチに失敗した場合のエラー・ハンドリング(27行目)が簡単に行えます。

10～23行では、ヘッダー部に置かれている文字列をパースしています。正規表現によるパターンマッチで 'descr','fortran_order','shape'それぞれの値(文字列)を読み取り、先に述べた Elixirデータ型に変換しています。Regex.run/3が返すパターンマッチの結果はリストとなっており、その第１要素は正規表現全体にマッチした文字列、第２要素は一つ目のキャプチャリング括弧にマッチした文字列、以下同様です。ここでは、一つ目のキャプチャリング括弧にマッチした文字列が欲しいので、caseとパターンマッチを組み合わせて結果リストの第２要素を切り出しています。

ちなみに、Elixirの正規表現は Erlang譲りで Perlの正規表現相当の表現力があります。詳しくは参考文献[2][3]を参照してください。

25行目では、まだファイルに残っているデータ部をバイナリとして全て読み出し、ここまでに得られた 'descr','fortran_order','shape'の値と合わせて %Npyデータ構造を生成しています。

話が前後しますが、2行目の File.open/3の返値は、ファイル・オープンに成功した場合は{:ok, 第３引数の関数func/1の返値}、オープンに失敗した場合は{:error, :enoent}となります。31～34行では、その返値を分類&加工し、本関数 Npy.load/1の返値としています。


```elixir:Npy.load
 1:  def load(fname) do
 2:    res = File.open(fname, [:read], fn file ->
 3:      with <<0x93, "NUMPY", major, _minor>> <- IO.binread(file, 8)
 4:      do
 5:        header_len = case major do
 6:          1 -> <<len::little-16>> = IO.binread(file, 2); len
 7:          _ -> <<len::little-32>> = IO.binread(file, 4); len
 8:        end
 9:
10:        header = IO.binread(file, header_len)
11:        descr = case Regex.run(~r/'descr': '([<=|>]?\w+)',/, header) do
12:          [_, descr] -> descr
13:          _ -> nil
14:        end
15:        fortran_order = case Regex.run(~r/'fortran_order': (True|False),/, header) do
16:          [_, "True" ] -> true
17:          [_, "False"] -> false
18:          _ -> nil
19:        end
20:        shape = case Regex.run(~r/'shape': \((\d+(,\s*\d+)*)\),/, header) do
21:          [_, shape|_] -> String.split(shape, ~r/,\s*/) |> Enum.map(&String.to_integer/1)
22:          _ -> nil
23:        end
24:
25:        {:ok, %Npy{descr: descr, fortran_order: fortran_order, shape: shape, data: IO.binread(file, :all)}}
26:      else
27:        _ -> {:error, "illegal npy file"}
28:      end
29:    end)
30:
31:    case res do
32:      {:ok, fun_res} -> fun_res
33:      err -> err
34:    end
35:  end
```
◆Tips2: 第３引数に関数を取る File.open/3を使用すれば、ファイルの closeが楽になる◆
◆Tips3: ファイルをBYTE列として読み出す場合は、IO.read/2ではなく IO.binread/2◆
◆Tips4: 関数内のコードでパターンマッチの成否を調べたい場合は withを検討しよう◆
◆Tips5: リストの第２要素を切り出すパターンは [_, second | _]◆

####(2)Npy.save関数
save/2関数は、引数fnameで指定したファイルに、Npyデータ構造を.npy形式で保存します。

2～6行では、Npyデータ構造が持つ情報から .npyファイル頭部となるバイナリを作成しています。3行目はヘッダー文字列の作成、4行目は64アライメントの為のパディングを生成してヘッダー文字列に連結、5行目はマジック・ワード、フォーマット・バージョン、ヘッダー部バイト長にヘッダー部を連結しています。このコードは、少々強引にパイプを使用して中間変数を使わないスタイルで書いてみました。とても可読性が悪いコードなので、反面教師として見るのが良いでしょうね。

7～12行では、ファイルへの書き出しを行っています。ファイル・オープンの成否を withでチェックするスタイルは小生の個人的な好みです。withの代わりに caseを用いるコードもよく見かけます。

```elixir:Npy.save
 1:  def save(fname, %Npy{}=npy) do
 2:    meta =
 3:      "{'descr': '#{npy.descr}', 'fortran_order': #{if npy.fortran_order,do: "True",else: "False"}, 'shape': (#{Enum.join(npy.shape, ", ")}), }"
 4:      |> (&(&1 <> String.duplicate(" ", 63-rem(byte_size(&1)+10, 64)) <> "\n")).()
 5:      |> (&(<<0x93,"NUMPY",1,0,byte_size(&1)::little-integer-16>> <> &1)).()
 6:
 7:    with {:ok, file} <- File.open(fname, [:write])
 8:    do
 9:      IO.binwrite(file, meta)
10:      IO.binwrite(file, npy.data)
11:      File.close(file)
12:    end
13:  end
```
◆Tips6: 関数の引数に構造体を期待する場合は、`%Npy{}=npy`の様にパターンマッチ◆
◆Tips7: 無名関数をパイプで利用する場合は、関数定義全体を丸括弧で括り".()"を付ける◆

####(3)Npy.to_list関数
to_list/1関数は、引数%Npy{descr: descr, ...}に渡された Npyデータ構造から、そのshapeと同じ次元数・形状を持つネストしたリストを生成して返します。

2～6行では、dataが持っているバイナリを descrに応じたデータ型に切り出し、フラットなリストを作っています。`for <<x::little-float-32 <- data>> do x end`は、バイナリをリストに変換する鉄板コード・パターンです。

8行目および13～14行では、フラットなリストを、shapeの次元数・形状を持つネストしたリストに織り上げています。主役はプライベート関数list_forming/2で、最も深い次元(shapeの末尾の次元)から始めて最も浅い次元(shapeの先頭の次元)に向かい、リストの要素を Enum.chunk_every/2用いて小ブロックに分割し、順にネストしたリストに織り上げていきます。この仕事は本質的にはループ処理なので、list_forming/2は末尾再帰の関数として記述しています。

```elixir:Npy.to_list
 1:  def to_list(%Npy{descr: descr, shape: shape, data: data}) do
 2:    flat_list = case descr do
 3:      "<f4" -> for <<x::little-float-32 <- data>> do x end
 4:      "<i4" -> for <<x::little-integer-32 <- data>> do x end
 5:      _ -> nil
 6:    end
 7:
 8:    if (flat_list), do: list_forming(Enum.reverse(shape), flat_list)
 9:  end
11:  def to_list(_), do: nil
12:
13:  defp list_forming([_],          formed), do: formed
14:  defp list_forming([size|shape], formed), do: list_forming(shape, Enum.chunk_every(formed, size))
```
◆Tips8: バイナリをリストに変換する鉄板コードパターン`for <<x::float <- binary>> do x end`

####(4)Npy.from_list関数
from_list/2関数は、引数xに渡されたネストしたリストを　descrで指定したデータ型の要素を持つ多次元配列とみなし、それを Npyデータ構造に変換して返します。

2～6行目では、descrに応じてリストx の要素をバイナリに変換する関数を選択しています。12行目が本来なら"処理の分岐"となるところを、関数渡しで置き換えるためです。

12行目では、上で選択した関数と Enum.reduce/3、List.flatten/1を組み合わせて、リストx からバイナリを生成しています。

11行目および18～19行では、リストx の最左端の要素を深さ順に訪れ、shapeの値を計算しています。こちらのプライベート関数calc_shape/1は、末尾再帰ではない再帰関数としてコンパクトに書いています。再帰の深さが小さいと見積もれる場合は、素直に表現するのも良いですね。


```elixir:Npy.from_list
 1:  def from_list(x, descr) when length(x) > 0 do
 2:    to_binary = case descr do
 3:      "<f4" -> fn x,acc -> acc <> <<x::little-float-32>> end
 4:      "<i4" -> fn x,acc -> acc <> <<x::little-integer-32>> end
 5:      _ -> nil
 6:    end
 7:
 8:    if to_binary do
 9:      %Npy{
10:        descr: descr,
11:        shape: calc_shape(x),
12:        data:  Enum.reduce(List.flatten(x), <<>>, to_binary)
13:      }
14:    end
15:  end
16:  def from_list(_, _), do: nil
17:
18:  defp calc_shape([item|_]=x), do: [Enum.count(x)|calc_shape(item)]
19:  defp calc_shape(_),          do: []
```
◆Tips9: 処理が分岐するコードを、関数渡しで置き換えられないか考えみよう◆

#4.動作確認
では、iexを起動して動作確認をしてみましょう。
テストデータとして Python&numpyで保存した "test/pred0.npy"を用意しました。
https://github.com/shoz-f/npy_ex

‥‥期待した動作をしているようですね。

```elixir
C:\home\shozo\Elixir\npy>iex -S mix
Interactive Elixir (1.10.4) - press Ctrl+C to exit (type h() ENTER for help)
iex(1)> {:ok, a} = Npy.load("test/pred0.npy")
{:ok,
 %Npy{
   data: <<174, 81, 140, 64, 246, 199, 28, 64, 192, 100, 255, 64, 222, 75, 156,
     64, 246, 101, 51, 65, 174, 225, 14, 64, 78, 219, 192, 65, 27, 251, 105, 64,
     18, 75, 161, 65, 86, 80, 206, 63, 212, 213, 181, 65, 178, 6, 250, ...>>,
   descr: "<f4",
   fortran_order: false,
   shape: [1, 10647, 4]
 }}

iex(2)> a.shape
[1, 10647, 4]

iex(3)> a.descr
"<f4"

iex(4)> b = Npy.to_list(a)
[
  [
    [4.384970664978027, 2.449704647064209, 7.981048583984375, 4.884261131286621],
    [11.212392807006836, 2.2325243949890137, 24.10708236694336,
     3.6559512615203857],
    [20.16165542602539, 1.6118266582489014, 22.729408264160156,
     1.95332932472229],
<中略>
    [356.31878662109375, 1.6345915794372559, 27.936599731445313,
     1.9644194841384888],
    [363.70611572265625, 1.5674588680267334, 27.086915969848633, ...],
    [372.02618408203125, 1.481452226638794, ...],
    [379.8704833984375, ...],
    [...],
    ...
  ]
]

iex(5)> c = Npy.from_list(b, "<f4")
%Npy{
  data: <<174, 81, 140, 64, 246, 199, 28, 64, 192, 100, 255, 64, 222, 75, 156,
    64, 246, 101, 51, 65, 174, 225, 14, 64, 78, 219, 192, 65, 27, 251, 105, 64,
    18, 75, 161, 65, 86, 80, 206, 63, 212, 213, 181, 65, 178, 6, 250, 63, 185,
    ...>>,
  descr: "<f4",
  fortran_order: false,
  shape: [1, 10647, 4]
}

iex(6)> c == a
true

iex(7)> Npy.save("test/c.npy", c)
:ok

iex(8)> {:ok, d} = Npy.load("test/c.npy")
{:ok,
 %Npy{
   data: <<174, 81, 140, 64, 246, 199, 28, 64, 192, 100, 255, 64, 222, 75, 156,
     64, 246, 101, 51, 65, 174, 225, 14, 64, 78, 219, 192, 65, 27, 251, 105, 64,
     18, 75, 161, 65, 86, 80, 206, 63, 212, 213, 181, 65, 178, 6, 250, ...>>,
   descr: "<f4",
   fortran_order: false,
   shape: [1, 10647, 4]
 }}

iex(9)> d == c
true
```


#5.TDD (Test Driven Development)？ - 裏話？？
少々裏話になりますが、上に紹介したコードは最初からあのような姿になっていた訳ではありません。何度かリファクタリングを行い、試行錯誤を行った結果の姿です。

Elixirをはじめ、インタラクティブな実行環境(REPL)が整った言語処理系では、小さなコード断片から実行と改良を繰り返しながらアイデアを盛り込み、プログラムを仕上げていく開発スタイルをとることがあります[*2]。

そんなとき、コードの改良によってプログラム全体が壊れていないかどうかを確かめる足場が欲しくなります。Elixirには素敵なユニット・テストの機能が付属しているので、コードがそこそこ動き始めたら、並行してテストを少しずつ書いて行くことが出来ます。

今回の練習課題においても、舞台裏では、下のようなテストを順次書き足すTDDなアクションを回していました。

[*2]Lispはこの開発スタイルをとる代表格

```elixir:test/npy_test.exs
defmodule NpyTest do
  use ExUnit.Case
  doctest Npy

  test "load dog_output0.npy" do
    assert {:ok, npy} = Npy.load("test/dog_output0.npy")
    assert %Npy{descr: "<f4", fortran_order: false, shape: [1,10647,4]} = npy
  end

  test "load pred0.npy" do
    assert {:ok, npy} = Npy.load("test/pred0.npy")
    assert %Npy{descr: "<f4", fortran_order: false, shape: [1,10647,4]} = npy
  end

  test "load illegal file" do
    assert {:error, "illegal npy file"} = Npy.load("test/test_helper.exs")
  end

  test "load absent file" do
    assert {:error, :enoent} = Npy.load("test/nofile.npy")
  end

  test "to_list" do
    assert {:ok, npy} = Npy.load("test/dog_output0.npy")
    assert [[[4.177847385406494, 2.4681198596954346,7.7223896980285645, 4.816765785217285]|_]] = Npy.to_list(npy)
  end
  
  test "from_list: vector" do
    assert %Npy{descr: "<i4", shape: [5], data: data} = Npy.from_list([1,2,3,4,5], "<i4")
    assert data == <<1,0,0,0, 2,0,0,0, 3,0,0,0, 4,0,0,0, 5,0,0,0>>
  end
  
  test "from_list: matrix" do
    assert %Npy{descr: "<f4", shape: [2,2], data: data} = Npy.from_list([[1.0, 0.0],[0.0, 1.0]], "<f4")
    assert <<1.0::little-32, 0.0::little-32, 0.0::little-32, 1.0::little-32>> == data
  end

  test "from_list: tensor" do
    assert {:ok, npy} = Npy.load("test/dog_output0.npy")
    assert npy == Npy.from_list(Npy.to_list(npy), "<f4")
  end

  test "from_list: []" do
    assert nil == Npy.from_list([], "<f4")
  end
  
  test "from_list: <p16" do
    assert nil == Npy.from_list([1,2,3,4,5], "<i16")
  end
end
```

#6.あとがき
さて、年末故にいつもとは違う口調の記事にしてみました。はい、気まぐれです:stuck_out_tongue_winking_eye:

Npyモジュールの発展形は‥‥まぁ～～ったくやる気がありませんが、Tensor演算が出来る様にしてみるのも楽しいかも知れませんね。その場合、どう考えても不利なElixirのリストで演算を行うのではなく、%Npyが持つバイナリをNIFsに渡し、BLASなりCUDAなりでガシガシと演算する戦術をとるのが吉かと‥‥

明日、6日目の Elixir Advent Calendar 2020 の担当は @zacky1972 さんです。よろしくお願いします。

それでは、また。
まだ早いですが、メリークリスマス:christmas_tree: 素敵な日々を:gift::champagne:

#参考文献
[1] [numpy.lib.format](https://numpy.org/devdocs/reference/generated/numpy.lib.format.html)
[2] [Elixir Documentation - Regex](https://hexdocs.pm/elixir/Regex.html#content)
[3] [Erlang Reference Manual - re](http://erlang.org/doc/man/re.html)
[4] [Elixir Documentation - bitstring](https://hexdocs.pm/elixir/Kernel.SpecialForms.html#%3C%3C%3E%3E/1)
[5] [Elixir Documentation - for](https://hexdocs.pm/elixir/Kernel.SpecialForms.html#for/1)