SwiftでN次元行列演算ライブラリMatftを作ってみた

## はじめに
初投稿です．
タイトルの通りなのですが，SwiftでN次元行列演算ライブラリ**Matft([https://github.com/jjjkkkjjj/Matft](https://github.com/jjjkkkjjj/Matft))**を作ってみました．(**Mat** rix演算をする・Swi **ft**で，の略でMatftです笑)

事の発端は，会社の先輩の「Swiftで３行３列の逆行列を求めるコードを書いてほしい」という一言でした．
SwiftにはPythonのNumpyのようなN次元行列演算ライブラリがあるだろうと思ったのですが，調べると意外にもないんですよね．．．
公式のAccelerateは使い勝手悪そうだし，有名らしい[surge](https://github.com/Jounce/Surge)も２次元まで？みたいでした．そんなこんなで，せっかくなので自作のN次元行列演算ライブラリを作ってみようと思いました．（３行３列の逆行列を求めるコードに対して，完全にオーバースペックですが笑）

さらにそんなこんなで，Matftができました．そしてせっかくなので共有してみようということで，現在に至ります．

## 概要
基本的にはPythonのNumpyにならって作成したので，関数名や使い方はNumpyとほぼ同じです．

###宣言
宣言は`ndarray`なる`MfArray`で多次元配列を生成します．

```swift
let a = MfArray([[[ -8,  -7,  -6,  -5],
                  [ -4,  -3,  -2,  -1]],
        
                 [[ 0,  1,  2,  3],
                  [ 4,  5,  6,  7]]])
print(a)
/*
mfarray = 
[[[	-8.0,		-7.0,		-6.0,		-5.0],
[	-4.0,		-3.0,		-2.0,		-1.0]],

[[	0.0,		1.0,		2.0,		3.0],
[	4.0,		5.0,		6.0,		7.0]]], type=Float, shape=[2, 2, 4]
*/
```
## 型
いろいろな型に対応させたかったので，それなりの型を用意しました．`dtype`ならぬ`MfType`です．

```swift
let a = MfArray([[[ -8,  -7,  -6,  -5],
                  [ -4,  -3,  -2,  -1]],
            
                 [[ 0,  1,  2,  3],
                  [ 4,  5,  6,  7]]], mftype: .Float)
print(a)
/*
mfarray = 
[[[	-8.0,		-7.0,		-6.0,		-5.0],
[	-4.0,		-3.0,		-2.0,		-1.0]],

[[	0.0,		1.0,		2.0,		3.0],
[	4.0,		5.0,		6.0,		7.0]]], type=Float, shape=[2, 2, 4]
*/
let aa = MfArray([[[ -8,  -7,  -6,  -5],
                  [ -4,  -3,  -2,  -1]],
            
                 [[ 0,  1,  2,  3],
                  [ 4,  5,  6,  7]]], mftype: .UInt)
print(aa)
/*
mfarray = 
[[[	4294967288,		4294967289,		4294967290,		4294967291],
[	4294967292,		4294967293,		4294967294,		4294967295]],

[[	0,		1,		2,		3],
[	4,		5,		6,		7]]], type=UInt, shape=[2, 2, 4]
*/
//Above output is same as numpy!
/*
>>> np.arange(-8, 8, dtype=np.uint32).reshape(2,2,4)
array([[[4294967288, 4294967289, 4294967290, 4294967291],
        [4294967292, 4294967293, 4294967294, 4294967295]],

       [[         0,          1,          2,          3],
        [         4,          5,          6,          7]]], dtype=uint32)

```

型一覧は以下のEnum型で定義しました．
※実を言うと，裏では`Float`か`Double`で保存しているので，`UInt`なんかは値が大きいとオーバーフローします．ただ，実用上は問題ないと思います．

```swift
 public enum MfType: Int{
    case None // Unsupportted
    case Bool
    case UInt8
    case UInt16
    case UInt32
    case UInt64
    case UInt
    case Int8
    case Int16
    case Int32
    case Int64
    case Int
    case Float
    case Double
    case Object // Unsupported
}
```

## Indexing
Numpyでいう`a[:, ::-1]`のようなスライスも`~<`で実装しました．（当初`~`で実装していたのですが，`~2`などはビットNotになりますね．．．アホでした）
-1のような負のインデックスも実装しました．（これが一番苦労したかもしれません．．．）

```swift
let a = Matft.mfarray.arange(start: 0, to: 27, by: 1, shape: [3,3,3])
print(a)
/*
mfarray = 
[[[	0,		1,		2],
[	3,		4,		5],
[	6,		7,		8]],

[[	9,		10,		11],
[	12,		13,		14],
[	15,		16,		17]],

[[	18,		19,		20],
[	21,		22,		23],
[	24,		25,		26]]], type=Int, shape=[3, 3, 3]
*/
print(a[2,1,0])
// 21
print(a[1~<3]) //same as a[1:3] for numpy
/*
mfarray = 
[[[	9,		10,		11],
[	12,		13,		14],
[	15,		16,		17]],

[[	18,		19,		20],
[	21,		22,		23],
[	24,		25,		26]]], type=Int, shape=[2, 3, 3]
*/
print(a[-1~<-3])
/*
mfarray = 
	[], type=Int, shape=[0, 3, 3]
*/
print(a[~<~<-1]) // print(a[~<<-1]) エイリアスで ~<<も等価です
/*
mfarray = 
[[[	18,		19,		20],
[	21,		22,		23],
[	24,		25,		26]],

[[	9,		10,		11],
[	12,		13,		14],
[	15,		16,		17]],

[[	0,		1,		2],
[	3,		4,		5],
[	6,		7,		8]]], type=Int, shape=[3, 3, 3]*/
```

Boolean Indexingも実装しました．いい感じ．

```swift
let img = MfArray([[1, 2, 3],
                   [4, 5, 6],
                   [7, 8, 9]], mftype: .UInt8)
img[img > 3] = MfArray([10], mftype: .UInt8)
print(img)
/*
mfarray = 
[[	1,		2,		3],
[	10,		10,		10],
[	10,		10,		10]], type=UInt8, shape=[3, 3]
*/
```

## その他関数一覧
ここからは，具体的な関数一覧です．主な計算はAccelerateに任せているので，計算時間はある程度担保されていると思います．

\* はmethodも存在することを意味します．つまり，`a`が`MfArray`であれば，`a.shallowcopy()`が使えます．
^ はmethodのみの関数です．

- 生成系

| Matft                      | Numpy             |
| -------------------------- | ---------------- |
| *Matft.shallowcopy | *numpy.copy       |
| *Matft.deepcopy    | copy.deepcopy     |
| Matft.nums         | numpy.ones * N    |
| Matft.arange       | numpy.arange      |
| Matft.eye          | numpy.eye         |
| Matft.diag         | numpy.diag        |
| Matft.vstack       | numpy.vstack      |
| Matft.hstack       | numpy.hstack      |
| Matft.concatenate  | numpy.concatenate |


- 変換系
 
| Matft                       | Numpy                    |
| --------------------------- | ----------------------- |
| *Matft.astype       | *numpy.astype            |
| *Matft.transpose    | *numpy.transpose         |
| *Matft.expand_dims  | *numpy.expand_dims       |
| *Matft.squeeze      | *numpy.squeeze           |
| *Matft.broadcast_to | *numpy.broadcast_to      |
| *Matft.conv_order   | *numpy.ascontiguousarray |
| *Matft.flatten      | *numpy.flatten           |
| *Matft.flip         | *numpy.flip              |
| *Matft.clip         | *numpy.clip              |
| *Matft.swapaxes     | *numpy.swapaxes          |
| *Matft.moveaxis     | *numpy.moveaxis          |
| *Matft.sort         | *numpy.sort              |
| *Matft.argsort      | *numpy.argsort           |
| ^MfArray.toArray | ^numpy.ndarray.tolist |


- ファイル関係

saveが未完成です．

| Matft                         | Numpy            |
| ----------------------------- | --------------- |
| Matft.file.loadtxt    | numpy.loadtxt    |
| Matft.file.genfromtxt | numpy.genfromtxt |

- 演算系

2行目が演算子です．

| Matft                          | Numpy                      |
| ------------------------------ | ------------------------- |
| Matft.add<br />+       | numpy.add<br />+           |
| Matft.sub<br />-       | numpy.sub<br />-           |
| Matft.div<br />/       | numpy.div<br />.           |
| Matft.mul<br />*       | numpy.multiply<br />*      |
| Matft.inner<br />*+    | numpy.inner<br />n/a       |
| Matft.cross<br />*^    | numpy.cross<br />n/a       |
|Matft.matmul<br />*&　　　|numpy.matmul<br />@　|
| Matft.equal<br />===   | numpy.equal<br />==        |
| Matft.not_equal<br />!==   | numpy.not_equal<br />!=        |
| Matft.less<br />< | numpy.less<br />< |
| Matft.less_equal<br /><= | numpy.less_equal<br /><= |
| Matft.greater<br />> | numpy.greater<br />> |
| Matft.greater_equal<br />>= | numpy.greater_equal<br />>= |
| Matft.allEqual<br />== | numpy.array_equal<br />n/a |
| Matft.neg<br />-       | numpy.negative<br />-      |


- 初等関数系

| Matft                    | Numpy       |
| ------------------------ | ---------- |
| Matft.math.sin   | numpy.sin   |
| Matft.math.asin  | numpy.asin  |
| Matft.math.sinh  | numpy.sinh  |
| Matft.math.asinh | numpy.asinh |
| Matft.math.sin   | numpy.cos   |
| Matft.math.acos  | numpy.acos  |
| Matft.math.cosh  | numpy.cosh  |
| Matft.math.acosh | numpy.acosh |
| Matft.math.tan   | numpy.tan   |
| Matft.math.atan  | numpy.atan  |
| Matft.math.tanh  | numpy.tanh  |
| Matft.math.atanh | numpy.atanh |
| Matft.math.sqrt | numpy.sqrt |
| Matft.math.rsqrt | numpy.rsqrt |
| Matft.math.exp | numpy.exp |
| Matft.math.log | numpy.log |
| Matft.math.log2 | numpy.log2 |
| Matft.math.log10 | numpy.log10 |
| *Matft.math.ceil | numpy.ceil |
| *Matft.math.floor | numpy.floor |
| *Matft.math.trunc | numpy.trunc |
| *Matft.math.nearest | numpy.nearest |
| *Matft.math.round | numpy.round |
| Matft.math.abs | numpy.abs |
| Matft.math.reciprocal | numpy.reciprocal |
| Matft.math.power | numpy.power |
| Matft.math.square | numpy.square |
| Matft.math.sign | numpy.sign |

- 統計系

| Matft                       | Numpy         |
| --------------------------- | ------------ |
| *Matft.stats.mean   | *numpy.mean   |
| *Matft.stats.max    | *numpy.max    |
| *Matft.stats.argmax | *numpy.argmax |
| *Matft.stats.min    | *numpy.min    |
| *Matft.stats.argmin | *numpy.argmin |
| *Matft.stats.sum    | *numpy.sum    |
| Matft.stats.maximum | numpy.maximum |
| Matft.stats.minimum | numpy.minimum |
| Matft.stats.sumsqrt | n/a |
| Matft.stats.squaresum | n/a |

- 線形代数系

| Matft                            | Numpy              |
| -------------------------------- | ----------------- |
| Matft.linalg.solve       | numpy.linalg.solve |
| Matft.linalg.inv         | numpy.linalg.inv   |
| Matft.linalg.det         | numpy.linalg.det   |
| Matft.linalg.eigen       | numpy.linalg.eig   |
| Matft.linalg.svd         | numpy.linalg.svd   |
| Matft.linalg.polar_left  | scipy.linalg.polar |
| Matft.linalg.polar_right | scipy.linalg.polar |
| Matft.linalg.normlp_vec | scipy.linalg.norm |
| Matft.linalg.normfro_mat | scipy.linalg.norm |
| Matft.linalg.normnuc_mat | scipy.linalg.norm |

## インストール
SwiftPMとCocoaPodとCarthageに対応しました．

### SwiftPM

- Import
  - Project > Build Setting > + ![Build Setting](https://user-images.githubusercontent.com/16914891/77144994-b0c72280-6aca-11ea-8633-1fb1a13ec74d.png)
  - 適宜選択
    ![select](https://user-images.githubusercontent.com/16914891/77144995-b1f84f80-6aca-11ea-8f4d-911bd96013cb.png) 
- アップデート
  - File >Swift Packages >Update to Latest Package versions
    ![update](https://user-images.githubusercontent.com/16914891/77145225-4367c180-6acb-11ea-98ea-8d7a5a2a669f.png)


### CocoaPods

- Podfile作成 (すでにある場合は無視)

  ```bash
  pod init
  ```

- `pod 'Matft'`をPodfileに追記

  ```bash
  target 'your project' do
    pod 'Matft'
  end
  ```

- インストール

  ```bash
  pod install
  ```

### Carthage（未確認）

-　Cartfileに作成し，読み込みます．

```bash
echo 'github "realm/realm-cocoa"' > Cartfile
carthage update ###or append '--platform ios'
```

- 出来上がった`Matft.framework`をプロジェクトに読み込めばOKです．

## Performance

Accelerateに任せているので，計算は担保されていると言いましたが，足し算だけ速度を計算してみました．
時間があれば他の関数も調べます．．．

|case| Matft                      |          Numpy          |
|:--:|----------------------------|----------------------- |
|1   | 1.14ms                     |962 µs   |
|2   | 4.20ms                     |5.68 ms            |
|3   | 4.17ms                     |3.92 ms            |

- Matft

```swift
func testPefAdd1() {
        do{
            let a = Matft.mfarray.arange(start: 0, to: 10*10*10*10*10*10, by: 1, shape: [10,10,10,10,10,10])
            let b = Matft.mfarray.arange(start: 0, to: -10*10*10*10*10*10, by: -1, shape: [10,10,10,10,10,10])
            
            self.measure {
                let _ = a+b
            }
            /*
             '-[MatftTests.ArithmeticPefTests testPefAdd1]' measured [Time, seconds] average: 0.001, relative standard deviation: 23.418%, values: [0.001707, 0.001141, 0.000999, 0.000969, 0.001029, 0.000979, 0.001031, 0.000986, 0.000963, 0.001631]
            1.14ms
             */
        }
    }
    
    func testPefAdd2(){
        do{
            let a = Matft.arange(start: 0, to: 10*10*10*10*10*10, by: 1, shape: [10,10,10,10,10,10])
            let b = a.transpose(axes: [0,3,4,2,1,5])
            let c = a.T
            
            self.measure {
                let _ = b+c
            }
            /*
             '-[MatftTests.ArithmeticPefTests testPefAdd2]' measured [Time, seconds] average: 0.004, relative standard deviation: 5.842%, values: [0.004680, 0.003993, 0.004159, 0.004564, 0.003955, 0.004200, 0.003998, 0.004317, 0.003919, 0.004248]
            4.20ms
             */
        }
    }

    func testPefAdd3(){
        do{
            let a = Matft.arange(start: 0, to: 10*10*10*10*10*10, by: 1, shape: [10,10,10,10,10,10])
            let b = a.transpose(axes: [1,2,3,4,5,0])
            let c = a.T
            
            self.measure {
                let _ = b+c
            }
            /*
             '-[MatftTests.ArithmeticPefTests testPefAdd3]' measured [Time, seconds] average: 0.004, relative standard deviation: 16.815%, values: [0.004906, 0.003785, 0.003702, 0.005981, 0.004261, 0.003665, 0.004083, 0.003654, 0.003836, 0.003874]
            4.17ms
             */
        }
}
```
- Numpy

```python
In [1]:
import numpy as np
#import timeit

a = np.arange(10**6).reshape((10,10,10,10,10,10))
b = np.arange(0, -10**6, -1).reshape((10,10,10,10,10,10))

#timeit.timeit("b+c", repeat=10, globals=globals())
%timeit -n 10 a+b
962 µs ± 273 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [2]:
a = np.arange(10**6).reshape((10,10,10,10,10,10))
b = a.transpose((0,3,4,2,1,5))
c = a.T
#timeit.timeit("b+c", repeat=10, globals=globals())
%timeit -n 10 b+c
5.68 ms ± 1.45 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [3]:
a = np.arange(10**6).reshape((10,10,10,10,10,10))
b = a.transpose((1,2,3,4,5,0))
c = a.T
#timeit.timeit("b+c", repeat=10, globals=globals())
%timeit -n 10 b+c
3.92 ms ± 897 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

## 個人的に気に入っているところ
- 名前
- 複数型対応
- Numpyぽさ
- Indexing
- Row Major,Column Majorに対応した点（[ここ](https://github.com/jjjkkkjjj/Matft/blob/master/Sources/Matft/util/common/order.swift)）
- Accelerateに渡すアルゴリズム（[ここ](https://github.com/jjjkkkjjj/Matft/blob/master/Sources/Matft/util/common/offset.swift#L110)）

## 微妙なところ
- まだまだオーバーヘッドが存在する
    - せっかくMfArrayのOrderを保持する`flag`を実装したのに，無駄なコピーをしている箇所がいくつかある
    - 生成系が`Array`→`MfArray`に変換している（vDSPとかを使いたい）
- Protocolがうまく使えていない


## 最後に
気まぐれで作成しましたが，思ったよりいい感じのものができました．ただただ僕の検索力不足なだけで，より良いライブラリがあるかもしれませんが，是非試していただけると嬉しいです．（環境依存のチェックができていませんので，是非お願いします．．．）
とにもかくにもいい勉強になりました．ありがとうございました．

## 参考
[numpy](https://numpy.org/)
[scipy](https://www.scipy.org/)
[Accelerate](https://developer.apple.com/documentation/accelerate)
[SwiftでNDArray書く](https://qiita.com/t-ae/items/44dd88dadbe5ad7b27f4)（テストケースを参考にさせていただきました．）