C++で一次元、二次元配列を取り扱うクラスを作る

# C++で一次元、二次元配列を取り扱うクラスを作る
C++でnumpyみたいな物を作りたかったといった趣旨です。
ライブラリが有用かどうかは KMeans を作って有用性を確認しました。
ソースはこちら https://github.com/ikuo0/cpp_number/blob/main/numxd.h

## KMeans実行部

使用している seeds_dataset.txtは https://archive.ics.uci.edu/ml/datasets/seeds からダウンロードします
\t の数がおかしいので \t\t を \t に置き換えてから使用しましょう

``` cpp
void TestKMeans() {
    MemoryManager mm;
    auto tsvData = TSV::Read("./seeds_dataset.txt");
    auto data = TSV::ToDouble(mm, tsvData);
    
    Num1D<int> n1d(mm);
    Num2D<double> n2d(mm);
    
    auto xIndexes = n1d.Arange(0, 7);
    auto xData = n2d.IndexingT(data, xIndexes);
    auto yData = n2d.ValT(data, 7);
    
    StandardScaler scaler(xData);
    scaler.Fit();
    auto scaledX = scaler.Transform(mm);
    
    KMeans km(3);
    km.Initialize(scaledX, KMeans::enumInitializeRandom);
    km.Training(scaledX, 100, 1e-5);
    
    printf("########## centroids\n");
    Dump2D(km.Centroids);
    
    auto initCentroids = km.GetInitCentroids(mm);
    auto centroids = km.GetCentroids(mm);
    auto predict = km.GetPredict(mm, scaledX);
    
    TSV::Write("./cp_init_means.txt", initCentroids);
    TSV::Write("./cp_means.txt", centroids);
    TSV::Write("./cp_predict.txt", predict);
}
```

だいぶスマートに書けているんじゃないかと思います。（少なくとも過去に自分が試した中では一番まともに見える）
＞MemoryManager mm;
これはメモリを管理するクラスで確保したメモリをデストラクタで全て破棄するように作ってあります。
外部からメモリ管理クラスを提供することで計算スコープの外に出ても参照を維持し、かつ MemoryManager の宣言スコープ内ではメモリを自動的に破棄します。

## KMeans本体のソース

``` cpp

#pragma once

#include "numxd.h"

class KMeans {
    public:
    enum {
        enumInitializeRandom,
    };
    MemoryManager mm;
    const int Clusters;
    Num2D<double> InitCentroids;
    Num2D<double> Centroids;
    
    KMeans(const int clusters):
        Clusters(clusters),
        InitCentroids(mm),
        Centroids(mm)
    {
        Num2D<double> n2d(this->mm);
        this->InitCentroids = n2d.Create(1, 1);
        this->Centroids = n2d.Create(1, 1);
    }
    
    void InitializeRandom(Num2D<double> x) {
        SpotNum1D<int> n1d;
        SpotNum2D<double> n2d;
        auto indexes = n1d.Arange(0, x.Row);
        auto shuffled = n1d.Shuffle(indexes);
        auto selected = n1d.Slice(shuffled, 0, this->Clusters);
        auto initCentroids = n2d.Indexing(x, selected);
        this->InitCentroids.Release();
        this->InitCentroids = n2d.Clone(this->mm, initCentroids);
    }
    void Initialize(Num2D<double> x, const int init) {
        if(init == this->enumInitializeRandom) {
            this->InitializeRandom(x);
        } else {
            throw Format("error in %s: %d, unknown initialize parameter", __FUNCTION__, __LINE__);
        }
    }
    Num1D<int> EStep(Num2D<double> means, Num2D<double> x) {
        SpotNum1D<double> n1d;
        SpotNum2D<double> n2d;
        auto distances = n2d.Create(x.Row, this->Clusters);
        for(int i = 0; i < x.Row; i += 1) {
            for(int cluster = 0; cluster < this->Clusters; cluster += 1) {
                distances[i][cluster] = n1d.CalcDistance(x.Ref(i), means.Ref(cluster));
            }
        }
        
        Num1D<int> myN1d(this->mm);
        auto predict = myN1d.Create(x.Row);
        for(int i = 0; i < x.Row; i += 1) {
            predict[i] = n1d.ArgMin(distances.Ref(i));
        }
        return predict;
    }
    
    Num2D<double> MStep(Num1D<int> predict, Num2D<double> x) {
        Num2D<double> myN2d(this->mm);
        auto means = myN2d.Create(this->Clusters, x.Col);
        for(int cluster = 0; cluster < this->Clusters; cluster += 1) {
            SpotNum1D<int> n1d;
            SpotNum2D<double> n2d;
            auto idxs = n1d.WhereEq(predict, cluster);
            auto ix = n2d.Indexing(x, idxs);
            auto mean = ix.Mean();
            means.Ref(cluster).Copy(mean);
        }
        return means;
    }
    
    double CalcMeansDistance(Num2D<double> a, Num2D<double> b) {
        SpotNum2D<double> n2d;
        auto ia = n2d.Clone(a);
        auto ib = n2d.Clone(b);
        auto subtract = ia - ib;
        auto power = subtract.Power(2);
        auto total = power.TotalX();
        return sqrt(total / (double)a.Row);
    }
    
    void Training(Num2D<double> x, int maxIter=100, double threshold=1e-5) {
        Num2D<double> myN2d(this->mm);
        this->Centroids.Release();
        this->Centroids = myN2d.Create(this->Clusters, x.Col);
        SpotNum2D<double> n2d;
        auto means = n2d.Clone(this->InitCentroids);
        for(int i = 0; i < maxIter; i += 1) {
            auto predict = this->EStep(means, x);
            auto newMeans = this->MStep(predict, x);
            
            predict.Release();
            
            double distance = this->CalcMeansDistance(means, newMeans);
            means.Copy(newMeans);
            newMeans.Release();
            if(distance < threshold) {
                break;
            }
        }
        this->Centroids.Copy(means);
    }
    
    Num2D<double> GetInitCentroids(MemoryManager& mm) {
        Num2D<double> n2d(mm);
        return n2d.Clone(this->InitCentroids);
    }
    
    Num2D<double> GetCentroids(MemoryManager& mm) {
        Num2D<double> n2d(mm);
        return n2d.Clone(this->Centroids);
    }
    
    Num1D<int> GetPredict(MemoryManager& mm, Num2D<double> x) {
        auto predict = this->EStep(this->Centroids, x);
        Num1D<int> n1d(mm);
        return n1d.Clone(predict);
    }
};

```

１２０行そこそこで可読性もそれなりに確保されてます、python＋numpy で頑張ると５０行以下で書ける内容だと思いますがＣ＋＋にしては頑張って縮められたかなと思います。

## メモリ管理のルール
### SpotNum2D<double> n2d; 宣言
SpotNum1D, SpotNum2D はそのスコープ内で計算が完結しなければいけません、スコープ外に計算結果をコピーする時は Clone関数にメンバ、グローバルに宣言したMemoryManagerを渡してコピーを作成します。
KMeans実行部ソースの
＞auto initCentroids = km.GetInitCentroids(mm);
この処理が計算結果をスコープ外にコピーする処理になります

### Num2D<double> myN2d(this->mm); 宣言
前述のMemoryManagerを外部から与えるコンストラクタです、この場合は this->mm の寿命に合わせて計算を完結させるよう設計します。

Num1D、Num2D のコンストラクタに与える MemoryManager クラスの寿命（スコープ）を意識しなければいけないのがこのライブラリのルールになります。

### predict.Release(); 等の Release 呼び出しについて
MemoryManager クラスのメモリ使用中フラグをクリアします、クリアされた領域は free されず再利用され都度mallocするよりも高速化されます。
ループ内では都度Releaseしないとメモリリークを起こしますので確保～解放の記述を忘れずに書きます。

全てデストラクタで自動的に破棄して確保～解放を全く意識する必要の無いライブラリにしたかったのですが無理でした。

## std::vectorでいいんじゃない？
正直、速度気にしないならそれでいいかなという感想です。
二次元配列を密に確保して計算速度を稼ぎたいという所まで追い詰められてる人には本記事のライブラリは有用かと思いますが、安全性の担保の無いライブラリ使うくらいなら自作したほうマシなのでやはり自分以外に使い所のある人は居ないかなと思いました。
以上です。