n個の要素からr個を選ぶ組み合わせを得るためのJuliaプログラム

# はじめに

n個の要素の中からr個の組み合わせを得るプログラム（C言語やFortranなど）はいくつも報告されています。また、Juliaでも`Combinatorics`パッケージで組み合わせの配列を得ることができます（[参考URL](https://www.it-swarm.dev/ja/julia/%E7%B5%84%E3%81%BF%E5%90%88%E3%82%8F%E3%81%9B%E3%82%92%E6%95%B0%E3%81%88%E3%82%8B%E3%82%B8%E3%83%A5%E3%83%AA%E3%82%A2%E9%96%A2%E6%95%B0%E3%81%AF%E4%BD%95%E3%81%A7%E3%81%99%E3%81%8B%EF%BC%88n%E3%81%AFk%E3%82%92%E9%81%B8%E6%8A%9E%EF%BC%89%EF%BC%9F/838938326/)）。

ただし、これらのプログラムでは得られた組み合わせを画面上に表示するだけだったり、一度にすべての組み合わせの結果を得るようになっています。したがって、組み合わせの数が多い場合にはメモリを消費してしまいます。

一方、[このページ](https://slpr.sakura.ne.jp/qp/combination/) にはFORTRANで書かれた関数があります。これは、引数としてインデックスの配列を設定して関数を呼び出すと次の組み合わせが返るようになっており、ループの中で使用できるようになっています（実際に私も使用したことがあります）。

そこで、この関数を参考にしてJuliaで関数を作成してみました（ただし、`ini`と`key`を求めるための2行を除いて原型をとどめていません）。


# 内容

## 関数定義

以下の内容を `MyCombination.jl` というファイル名で保存します。ファイル名とモジュール名は適宜変更して下さい。

~~~~{.julia}

module MyCombination

export get_combination!, get_residual

function get_combination!(n, idx)
    r = length(idx)
    ini = 1:r
    key = @. (idx - ini) == (n - r)
    last0 = findlast(x -> x == 0, key)
    if !isnothing(last0)
        if last0 == r
            idx[r] += 1
        else
            idx[last0:r] = idx[last0] .+ collect(1:(r - last0 + 1))
        end
    else
        idx[:] = ini[:]
    end
    return
end


function get_residual(n, idx)
    filter(i -> !(i in idx), 1:n)
end

end # module
~~~~~~~~~~~~~




## 関数の説明

### `get_combination!`関数
#### 引数
- `n`: 全要素数
- `r`: 選択する要素数
- `idx`: `r`の長さの配列（通常は`collect(1:r)`で開始）

#### 使用方法
始めに`idx`に1からrまでの連続した数値を代入します。その後に`get_combination!`関数を呼び出すと`idx`が更新されて次の組み合わせが設定されます。

最後の組み合わせを与えた場合には、最初の組み合わせが戻ります。

n = 5, r = 3 の場合

- [1, 2, 3] → [1, 2, 4]
- [1, 2, 5] → [1, 3, 4]
- [1, 4, 5] → [2, 3, 4]
- [3, 4, 5] → [1, 2, 3]（最後の組み合わせ → 最初の組み合わせ）

### `get_residual`関数

1からnまでの配列の中から`idx`に含まれない数値（補集合）の配列を返す関数

#### 引数
- `n`: 全要素数
- `idx`: `r`の長さの配列

#### 返り値
- 長さ `n - r` の配列

## サンプル

初めに `idx`変数に最初の組み合わせを代入しておき、forループ内で何らかの処理を行い、ループの最後で `get_combinatin!`関数を呼び出します。

実際には`idx`変数を使用して目的の配列から要素を抜き出す場面が多いと思います。

ループの回数は n個の中からr個を選ぶ組み合わせ数のため、`binomial(n, r)` で求めることができます。


~~~~{.julia}

include("MyCombination.jl")
using .MyCombination

n = 5; r = 3
idx = collect(1:r)

for i in 1:binomial(n, r)
    println(idx, get_residual(n, idx))
    get_combination!(n, idx)
end
~~~~~~~~~~~~~


~~~~
[1, 2, 3][4, 5]
[1, 2, 4][3, 5]
[1, 2, 5][3, 4]
[1, 3, 4][2, 5]
[1, 3, 5][2, 4]
[1, 4, 5][2, 3]
[2, 3, 4][1, 5]
[2, 3, 5][1, 4]
[2, 4, 5][1, 3]
[3, 4, 5][1, 2]
~~~~





## `get_combination!`関数の簡単な解説

自分自身がこのアルゴリズムを理解できていないのですが（おぃ）、簡単に説明します。また、説明における用語の使い方が間違っている可能性があります。

ここでは n = 5, r = 3 の場合を例とします。

### `key`変数について
初めの数字を1、2番目の数字を2とした場合、最後の組み合わせは [1, 2, 5] となります（これを`idx`変数とします）。このとき、初めの組み合わせ（これを`ini`変数とします）である [1, 2, 3] との差をとると [0, 0, 2] となります。この3番目の要素の数字である2は n - r と一致します。

したがって、"`idx`変数と`ini`変数の値の差"の各要素がそれぞれ `n - r` と一致している場合には、その要素が取りうる組み合わせの最後の数字になっていることを示します。
逆に、`n - r` と異なる場合には、その要素ではまだ取りうる組み合わせが残っていることを示します。

この差が `n - r` と一致していれば1、一致していなければ0が`key`変数に代入されます。


同様に考えると、`idx = [1, 4, 5]` の場合には `key = [0, 1, 1]`、`idx = [3, 4, 5]` の場合には `key = [1, 1, 1]` となります。このように `key`変数は最後の要素から1になっていきます（[0, 1, 0] という値は取りません）。

### `last0`変数について
`key`変数の要素のうち、最後の0のインデックスを取得しています。`key`変数の要素がすべて1（最後の組み合わせ）の場合には `nothing` が返ります。

### 条件分岐について
1. `last0`変数の値が `nothing` の場合

    `idx`変数が最後の組み合わせのため、最初の組み合わせを返します。

2. `last0`変数の値が `r` と等しい場合

    取りうる組み合わせの最後まで達していないため、最後の要素 `idx[r]` の値を +1 します。
    - [1, 2, 3] の場合には last0 = 3 となるため、最後の要素を更新して [1, 2, 4]
    - [1, 3, 4] の場合にも同様に [1, 3, 5]

3. `last0`変数の値が `r` と異なる（小さい）場合

    取りうる組み合わせの最後まで達している場合です。この場合には繰り上がりが生じて、次の組み合わせのセットに移行します。このとき、`last0`変数で示す要素を+1し、さらに最後の要素まで1つずつ数字を大きくします。
    - [1, 2, 5] の場合には last0 = 2 となるため、2番目の要素から最後の要素まで更新して [1, 3, 4]
    - [1, 4, 5] の場合には last0 = 1 となるため、最初の要素から最後の要素まで更新して [2, 3, 4]


# 最後に

今回、組み合わせを選ぶアルゴリズムについて理解できたような気がします。まだほとんどJuliaでプログラムを書いたことはありませんが、少しずつ勉強していきたいと思います（実際にはメインで使うのはRなのでメインになることはないのですが……）。