nondominated sortを高速に計算する

## はじめに
表題の通り，高速なnondominated sortの計算に必要な関数を作りました．
nondominated sortとは主に多目的最適化においてサンプルの優劣をつけるために利用されるSort方法です．

## 追記 (2022年4月21日)
Githubに[コードを公開](https://github.com/nabenabe0928/fast_pareto)し，`pip install`で利用できるようにしました．

```shell
pip install fast-pareto
```

以下の方法でimport可能です．
```python
from fast_pareto import is_pareto_front, nondominated_rank
```

## 実装
以下の `nondominated_rank` 関数では入力 `costs (shape := (n_points, n_costs))` に対して，`ranks (shape := (n_points, ))` を返します．

```python
import numpy as np

def is_pareto_front(costs):
    on_front_indices = np.arange(costs.shape[0])
    (n_points, _) = costs.shape
    next_index = 0

    while next_index < len(costs):
        nd_mask = np.any(costs < costs[next_index], axis=1)
        nd_mask[next_index] = True
        # Remove dominated points
        on_front_indices, costs = on_front_indices[nd_mask], costs[nd_mask]
        next_index = np.sum(nd_mask[:next_index]) + 1

    mask = np.zeros(n_points, dtype=np.bool8)
    mask[on_front_indices] = True
    return mask

    
def nondominated_rank(costs):
    ranks = np.zeros(len(costs), dtype=np.int32)
    rank = 0
    indices = np.arange(len(costs))
    while indices.size > 0:
        on_front = is_pareto_front(costs)
        ranks[indices[on_front]] = rank
        # Remove pareto front points
        indices, costs = indices[~on_front], costs[~on_front]
        rank += 1

    return ranks
```

## ベンチマーク

### Input

```python
for n_points in [100, 1000, 10000]:
    for n_costs in [1, 5, 10, 50]:
        print(f"n_points={n_points}, n_costs={n_costs}")
        %time nondominated_rank(np.random.normal(size=(n_points, n_costs)))
        print("\n")
```

### Output

```
n_points=100, n_costs=1
CPU times: user 10.7 ms, sys: 0 ns, total: 10.7 ms
Wall time: 10.2 ms


n_points=100, n_costs=5
CPU times: user 3.29 ms, sys: 0 ns, total: 3.29 ms
Wall time: 3.3 ms


n_points=100, n_costs=10
CPU times: user 3 ms, sys: 0 ns, total: 3 ms
Wall time: 3 ms


n_points=100, n_costs=50
CPU times: user 3.57 ms, sys: 0 ns, total: 3.57 ms
Wall time: 3.57 ms


n_points=1000, n_costs=1
CPU times: user 105 ms, sys: 0 ns, total: 105 ms
Wall time: 105 ms


n_points=1000, n_costs=5
CPU times: user 37.8 ms, sys: 0 ns, total: 37.8 ms
Wall time: 37.8 ms


n_points=1000, n_costs=10
CPU times: user 53.5 ms, sys: 0 ns, total: 53.5 ms
Wall time: 53.5 ms


n_points=1000, n_costs=50
CPU times: user 90.4 ms, sys: 0 ns, total: 90.4 ms
Wall time: 90.4 ms


n_points=10000, n_costs=1
CPU times: user 3.36 s, sys: 0 ns, total: 3.36 s
Wall time: 3.36 s


n_points=10000, n_costs=5
CPU times: user 1.22 s, sys: 0 ns, total: 1.22 s
Wall time: 1.22 s


n_points=10000, n_costs=10
CPU times: user 2.72 s, sys: 0 ns, total: 2.72 s
Wall time: 2.72 s


n_points=10000, n_costs=50
CPU times: user 14.4 s, sys: 0 ns, total: 14.4 s
Wall time: 14.4 s
```

## (non-dominant set or 非劣解に関する) 補足
```math
\newcommand{\xv}{\boldsymbol{x}}
```

多出力関数 $f: \mathbb{R}^D \rightarrow \mathbb{R}^M$ の各目的最小化問題を考える．
任意の $i \in [1, M]$ に関して，$f_i(\xv) \leq f_i(\xv^\prime)$ であり，ある $i \in [1, M]$ が存在し，$f_i(\xv) < f_i(\xv^\prime)$ であるとき $f(\xv)$ は$f(\xv^\prime)$を **dominate (優越)** するという．
仮に$f(\xv)$ を優越するような $f(\xv^\prime)$ が存在しない場合にそのような $f(\xv)$ はパレート最適と呼ばれて，non-domination rankが 1 であると定義される．
また，non-domination rankが $n$ もしくは $n$ 以下であるような $f(\xv^\prime)$ を全て除外した集合において$f(\xv)$がパレート最適を達成する場合にそのような $f(\xv)$ は non-domination rankが $n + 1$ であると定義される．