MLJ.jl に自作のモデルを追加する

# 書いてあること
自作の回帰モデルを，julia の機械学習フレームワークである`MLJ.jl`から"最低限"使えるようにする方法が書いてあります．
ただし，公式パッケージとして使えるようにするのではなく，あくまでも自分の計算機の中だけで使えるようにします．
`MLJ.jl`そのものについてはほとんど書いていないので，知りたい人は公式ドキュメント（https://alan-turing-institute.github.io/MLJ.jl/dev/ ）を眺めたり，ドキュメントの Getting Started で遊んでみてください．

# 環境
```
julia> versioninfo()
Julia Version 1.7.0
Commit 3bf9d17731 (2021-11-30 12:12 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin19.5.0)
  CPU: Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, icelake-client)
Environment:
  JULIA_NUM_THREADS = 4
```

# 準備
まず，追加する自作モデルを用意します．
適当な作業用のディレクトリで REPL に入って，

```julia
(@v1.7) pkg> generate MyMLJModel
  Generating  project MyMLJModel:
    MyMLJModel/Project.toml
    MyMLJModel/src/MyMLJModel.jl

shell> cd MyMLJModel
```
として説明用のモジュールを作って,そのディレクトリに移動しておきます．

`./src/MyMLJModel.jl`の内容を以下のように編集します．
これはテキトーに作った説明用の自作モデルなので，モデルとしては非常にアレです．
このモデルは一応線形回帰になってますが,最小二乗法的なことはせずに，`a`回乱数を振って一番良かったものを重みとします．
パラメータを指定してデータを突っ込むと何か出てくるブラックボックスとでも思ってください．

```julia
module MyMLJModel

using Random

export MyModel, myfit!, mypredict

mutable struct MyModel
    a::Int
    weights::Vector{Float64}
end

MyModel(a) = MyModel(a, Float64[])

function myfit!(model, X, y)
    n_vars = size(X)[1]
    best_error = Inf
    for i = 1:model.a
        weights = rand(n_vars)
        y_hat = [sum(weights[i] * X[i, j] for i = 1:n_vars) for j = 1:size(X)[2]]
        tmp_error = rmse(y_hat, y)

        if tmp_error < best_error
            model.weights = weights
            best_error = tmp_error
        end
    end
end

mypredict(model, X) = [sum(model.weights[i] * X[i, j] for i = 1:size(X)[1]) for j = 1:size(X)[2]]

rmse(a, b) = sqrt(sum(x -> x^2, a[i] - b[i] for i in eachindex(a)) / length(a))

end # module

```

とりあえずこのモデルを動かしてみます．
`script`というディレクトリを作って，`script/example1.jl`でモデルを動かしてみます．

```julia
using MyMLJModel

f(x) = 0.2x[1] + 0.4x[2] + 0.6x[3]
X = rand(3, 10)
y = [f(X[:, i]) for i = 1:10]

mymodel = MyModel(100)
myfit!(mymodel, X, y)
yhat = mypredict(mymodel, X)

@show y
@show yhat
```

テキトーにデータ`X`, `y`を作って`myfit!`で学習して，`mypredict`で予測します．
ここで，`X`のサイズが `(3, 10)`になっていることに注意してください．`X`の行が特徴量に対応しています．つまり，1列（`X[:, i]`)が1つのインスタンスに対応しています．

```julia
(@v1.7) pkg> activate .
  Activating project at `~/hoge/MyMLJModel`

(MyMLJModel) pkg> add Random
   Resolving package versions...
    Updating `~/hoge/MyMLJModel/Project.toml`
  [9a3f8284] + Random
    Updating `~/hoge/MyMJLModel/Manifest.toml`
  [9a3f8284] + Random
  [ea8e919c] + SHA
  [9e88b42a] + Serialization

julia> using Revise

julia> include("script/example1.jl")
[ Info: Precompiling MyMLJModel [68c9983a-15f9-41e9-9249-e301e2c6b7d7]
┌ Warning: Package MyMLJModel does not have Random in its dependencies:
│ - If you have MyMLJModel checked out for development and have
│   added Random as a dependency but haven't updated your primary
│   environment's manifest file, try `Pkg.resolve()`.
│ - Otherwise you may need to report an issue with MyMLJModel
└ Loading Random into MyMLJModel from project dependency, future warnings for MyMLJModel are suppressed.
y = [0.5601453907761953, 0.6879910740562954, 0.3530135116286392, 0.5567237922697389, 0.5370888314572486, 1.1237237581585515, 0.8041827332557248, 0.2703739272387356, 0.1411133336353244, 0.7799328076341898]
yhat = [0.49971973364491074, 0.7010051972588169, 0.3424426329476714, 0.42568976463547314, 0.5650044702311595, 1.025868358799277, 0.7260945619574919, 0.25485599443926044, 0.12670020159703554, 0.7674023144526892]
10-element Vector{Float64}:
 0.49971973364491074
 0.7010051972588169
 0.3424426329476714
 0.42568976463547314
 0.5650044702311595
 1.025868358799277
 0.7260945619574919
 0.25485599443926044
 0.12670020159703554
 0.7674023144526892
```

warning が出てますが，ここでは問題ないので一旦無視して，`@show`で表示された`y`と`yhat`を確認してみてください．`yhat`がそれっぽい値になってると思います．


ここまでで，説明用の自作のモデルの準備は終了です．

# `MLJModelInterface.jl`を使って自作モデルを MLJ から使えるようにする
MLJ から自作のモデルを利用できるようにする機能は`MLJModelInterface.jl`で提供されています．
以下の3つを実装すると，とりあえず自分のモデルが MLJ から使えるようになります．

- MLJ で使用するモデル型
- `fit`
- `predict`

上から順番に実装していきます．

## モデル型を定義する
MLJ におけるモデルとは，学習に使用するハイパーパラメータを格納するものです．
ここで使っている自作モデルのハイパーパラメータは繰り返し回数`a`のみなので，MLJのモデルは`a`のみを持ちます．学習することで得られるパラメータ（ここでは重み）は含みません．

`@mlj_model` というマクロを使ってモデルを定義します．
上の方では`MyModel`という名前の構造体が定義されていましたが，ここでは`MyRegressor`として新しくMLJ用のモデルを定義します．

```julia
import MLJBase
import MLJModelInterface

const MMI = MLJModelInterface

export MyRegressor

MMI.@mlj_model mutable struct MyRegressor <: MMI.Deterministic
    a::Int = 100::(_ > 1)
end
```

パラメータには初期‍が設定されている必要がります．`a::Int = 100`のようにして初期値を定義します．また，オプションとしてパラメータの範囲も定義することができ，`::(_ > 1)` と追記することで，`a > 1`という制約を設けられます．

`@mlj_model`で生成されるコンストラクタの引数はすべてキーワード引数になります．

あと，外部から使うために`export MyRegressor`も追記しておいてください．

`MLJModelInterface`は長いので，本稿のプログラム中では`MMI`と書けるようにしておきます．

## `fit` を実装する
`MMI.fit`を実装します． まず実装例を示します．

```julia
function MMI.fit(model::MyRegressor, verbosity, X, y, w = nothing)
    Xmat = MLJBase.matrix(X, transpose = true)
    mymodel = MyModel(model.a)
    myfit!(mymodel, Xmat, y)
    fitresult = mymodel
    report = nothing
    cache = nothing

    return (fitresult, cache, report)
end
```

最初の引数`model`には自分のモデルの型を書いておきます．

`X`は入力で，行列や`DataFrame`などが入ってきます．ここで注意点があるのですが，MLJ では`X`の列が特徴量に対応しています．つまり，1行が1つのインスタンスに対応しています．自作のモデルでそうなっていない場合，`X`を変換しておく必要があります．`MLJBase．jl`というパッケージの`matrix()`を使うと，`DataFrame`を行列に変換したり，転置をとったりできます．

`MMI.fit`の内部では，今まで使っていた自作モデルの`MyModel`を学習させます．自作モデルで使っていた関数等を使って学習させます．

そして，`fitresult`,`cache`, `report`を定義します．
`fitresult`には`predict`に必要な情報を入れます．`predict`では自作モデル側の`mypredict`が呼びたくて，その`mypredict`は引数に`mymodel`が必要なので，`fitresult = mymodel`としておきます．

`cache`はここでは使わないので`nothing`とします．不要でも`nothing`として返す必要があります．

`report`は何かユーザーが知りたい情報を返すのに使う`NamedTuple`です．ここでは自作モデルの重みを返すことにします．

## `predict`を実装する

`MMI.predict`を実装します．まず実装例を示します．

```julia
function MMI.predict(model::MyRegressor, fitresult, Xnew)
    Xmat = MLJBase.matrix(Xnew, transpose = true)
    return mypredict(fitresult, Xmat)
end
```

最初の引数`model`には自分のモデルの型を書いておきます．
ここでも`fit`と同様に,`X`を`mypredict`に渡すために`MLJBase.matrix`で変換しておきます．
あとは`mypredict`を呼んでその結果を返すだけです．

# 自作モデルをMLJから使ってみる
これで自作モデルをMLJから使うための最低限の用意が整いました．
試しに，以下に示す`script/example2.jl`で動作を確認してみます．

```julia
using MyMLJModel
using MLJ

f(x) = 0.2x[1] + 0.4x[2] + 0.6x[3]
X = rand(10, 3) # n × p
y = [f(X[i, :]) for i = 1:10]

model = MyRegressor(a = 100)
mach = machine(model, X, y)

evaluate!(mach, measure = rms)
```

実行すると，以下のような結果が出てくると思います．

```julia
julia> include("script/example2.jl")
Evaluating over 6 folds: 100%[=========================] Time: 0:00:00
PerformanceEvaluation object with these fields:
  measure, measurement, operation, per_fold,
  per_observation, fitted_params_per_fold,
  report_per_fold, train_test_pairs
Extract:
┌────────────────────────┬─────────────┬───────────┬─────────────────────────────────────────────────┐
│ measure                │ measurement │ operation │ per_fold                                        │
├────────────────────────┼─────────────┼───────────┼─────────────────────────────────────────────────┤
│ RootMeanSquaredError() │ 0.0691      │ predict   │ [0.14, 0.0514, 0.0208, 0.0332, 0.00209, 0.0701] │
└────────────────────────┴─────────────┴───────────┴─────────────────────────────────────────────────┘
```

チューニングもしてみます．

```julia
model = MyRegressor(a = 100)
r = range(model, :a, lower = 10, upper = 1000)
tuned_model = TunedModel(model = model, resampling = CV(nfolds = 5), tuning = Grid(resolution = 10), range = r, measure = rms)
fit!(mach)
```

```julia
julia> mach.report.plotting.measurements
10-element Vector{Float64}:
 0.02442430654777976
 0.06935442889813258
 0.01598692079500873
 0.04161092341230241
 0.0268491382201212
 0.056923740156152744
 0.03366051708222645
 0.039598241639119326
 0.03212408357200946
 0.03651355084042113

julia> fitted_params(mach).best_model
MyRegressor(
    a = 890)

julia> mach.report.best_history_entry
(model = MyRegressor,
 measure = [RootMeanSquaredError()],
 measurement = [0.01598692079500873],
 per_fold = [[0.016621197285246875, 0.008638543600675262, 0.007427216316183303, 0.025226911542651647, 0.01534467137495609]],)
```
`a`回試行できるので`a`が大きい方が良さそうですが，10から1000までしか動かしてませんし運ゲーなので，それっぽい結果は得られてると思います．

これで自作モデルをMLJから最低限使えるようになりました．