C++標準ライブラリでNeural Network

C++の標準ライブラリで、手書き数字データセットを分類できるニューラルネットワークを作成し、手書き数字文字認識及びアヤメの品種推定を行いました。
精度が90%程で、手書き数字認識の割には、精度が低いかなと思うので、どこかで間違えている気がします。間違いを見つけてくださった方は、コメントいただけると幸いです。

コードはVisual Studioのプロジェクトとして、[GitHub](https://github.com/takoyaki-3/c-NeuralNetwork)にアップしてあります。

以下、序盤は手書き数字分類に関してです。

## 実行の様子

![](https://gyazo.com/dbc69790775131bcfdb18bf2c05c53a9.png)

## 精度に関して

数字認識の場合、ノード数が160を超えるあたりからあまり精度に変化が見られなくなりました。また、epoch数に関しても、３を超えたあたりから、乱高下を繰り返すようになりました。

![](https://gyazo.com/9cb138ff7e196e7c05aae5627714e0b4.png)

機能ごとに、何をしているのか追ってみます。

## 行列演算用のクラス
ニューラルネットワークを作成するにあたり、行列演算を行うことは必須だと考え、行列演算を行うためのクラスを定義しました。
今回は、次のような機能を持つクラスを定義しました。

|メソッド|役割|
|---|---|
|matrix()|初期化関数。1×1行列を作成する|
|matrix(int nn,int mm)|初期化関数。n×m行列を作成する|
|matrix(matrix* a)|行列を複製する|
|void newmat(int nn, int mm)|nn×mm行列を作成する|
|void copy(matrix* a)|行列をコピーする|
|void put(int i, int j, matrixelement e)|i,jにeを代入する|
|matrixelement get(int i, int j)|i,jの要素を取得する|
|void random_m05to05()|すべての要素に-0.5～0.5のランダムな値を代入する|
|bool plus(matrix a)|引数として与えられた行列を加算する|
|bool minus(matrix a)|引数として与えられた行列をマイナスする|
|bool dot(matrix a, matrix b)|行列の掛け算をする|
|void output()|行列の値を標準出力へ出力する|
|void sigmoid()|すべての要素をシグモイド関数にかけた結果に置き換える|
|void t_set(matrix* a)|引数の行列を転置して代入|
|void scalar(matrixelement k)|すべての要素をkでスカラー倍する|
|void scalar_plus(matrixelement k)|すべての要素にすべてkを加算する|
|void scalar_division(matrixelement k)|すべての要素をスカラーで除算する|
|bool multiplication(matrix* a)|引数として与えられた行列を、対応する要素に掛け合わせる|
|void clear()|行列用メモリを解放する|

for文など、基本的な文法だけで実装することができました。

```cpp
class matrix {
private:
public:
	// Matrix body
	matrixelement* mat;

	// Matrix size
	int n, m;

	// Init matrix
	matrix() {
		n = m = 1;
		mat = new matrixelement[n * m];
		memset(mat, 0, sizeof(matrixelement) * n * m);
	}
	matrix(int nn, int mm) {
		mat = new matrixelement[nn * mm];
		memset(mat, 0, sizeof(matrixelement) * nn * mm);
		n = nn;
		m = mm;
	}
	matrix(matrix* a) {
		mat = new matrixelement[a->n * a->m];
		n = a->n;
		m = a->m;
		memcpy(mat, a->mat, m * n * sizeof(matrixelement));
	}
	void newmat(int nn, int mm) {
		n = nn;
		m = mm;
		mat = new matrixelement[n * m];
		memset(mat, 0, sizeof(matrixelement) * n * m);
	}
	void copy(matrix* a) {
		n = a->n;
		m = a->m;
		delete[] mat;
		mat = new matrixelement[a->n * a->m];
		memcpy(mat, a->mat, m * n * sizeof(matrixelement));
	}
	void put(int i, int j, matrixelement e) {
		if (i >= n || j >= m || i < 0 || j < 0) return;
		mat[i * m + j] = e;
	}
	matrixelement get(int i, int j) {
		if (i >= n || j >= m || i < 0 || j < 0) return -1;
		return mat[i * m + j];
	}
	void random_m05to05() {
		std::random_device rnd;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] = (matrixelement)(rnd() % RAND) / (matrixelement)RAND - 0.5f;
	}
	bool plus(matrix a) {
		if (a.n != n || a.m != m) return false;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] += a.mat[i * m + j];
		return true;
	}
	bool minus(matrix a) {
		if (a.n != n || a.m != m) return false;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] -= a.mat[i * m + j];
		return true;
	}
	bool dot(matrix a, matrix b) {
		if (a.m != b.n) return false;
		n = a.n;
		m = b.m;
		delete[] mat;
		mat = new matrixelement[n * m];
		memset(mat, 0, sizeof(matrixelement) * n * m);
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < m; j++) {
				for (int k = 0; k < a.m; k++) {
					put(i, j, get(i, j) + a.get(i, k) * b.get(k, j));
				}
			}
		}
		return true;
	}
	void output() {
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < m; j++) {
				printf("[%f]	", get(i, j));
			}
			cout << "\n";
		}
		cout << "\n";
	}
	void sigmoid() {
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] = sigmoidfunction(mat[i * m + j]);
	}
	void t_set(matrix* a) {
		n = a->m;
		m = a->n;
		delete[] mat;
		mat = new matrixelement[n * m];
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] = a->mat[j * a->m + i];
	}
	void scalar(matrixelement k) {
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] *= k;
	}
	void scalar_plus(matrixelement k) {
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] += k;
	}
	void scalar_division(matrixelement k) {
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] /= k;
	}
	bool multiplication(matrix* a) {
		if (n != a->n) return false;
		if (m != a->m) return false;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				mat[i * m + j] *= a->mat[i * m + j];
	}
	void clear() {
		delete[] mat;
	}
};
```

## ニューラルネットワークのクラス

行列演算のクラスに機能を追加するのと並行で、ニューラルネットワークのプログラムを書いていきました。

### メンバー変数
|変数名|役割|
|---|---|
|int input_nodes|入力層のノード数|
|int hidden_nodes|隠れ層のノード数|
|int output_nodes|出力層のノード数|
|int hidden_layer|隠れ層を複数構成にしようとしていた際の残骸(常に1)|
|vector<matrix> w|隠れ層を複数構成にしようとしていた際の残骸（隠れ層の行列が格納される）|
|matrix wo|出力層の重みを記憶する行列|

### メンバー関数

|メソッド名|役割|
|---|---|
|void save(string fname)|訓練済みのニューラルネットワークを保存する関数|
|void load(string fname)|訓練済みのニューラルネットワークを読み込む関数|
|void init()|初期化関数。設定されたパラメータでネットワークを作成する|
|void train(matrix* input, matrix* ans, float rate)|入力と出力の行列を指定し、学習する関数|
|void query(matrix* input, matrix* output)|入力から結果を得る関数|

順に見ていきます。

#### save 関数

save関数では、学習済みネットワークの各ノードでの重みをファイルに書き出します。
各層のサイズを出力した後に、各層の重みを表す行列を書き出します。

```cpp
void save(string fname) {
	FILE* fp = fopen(fname.c_str(), "wb");
	fwrite(&input_nodes, sizeof(int), 1, fp);
	fwrite(&hidden_nodes, sizeof(int), 1, fp);
	fwrite(&output_nodes, sizeof(int), 1, fp);
	fwrite(&hidden_layer, sizeof(int), 1, fp);
	for (int i = 0; i < w.size(); i++) {
		fwrite(&w[i].n, sizeof(int), 1, fp);
		fwrite(&w[i].m, sizeof(int), 1, fp);
		fwrite(w[i].mat, sizeof(matrixelement), w[i].n * w[i].m, fp);
	}
	fwrite(&wo.n, sizeof(int), 1, fp);
	fwrite(&wo.m, sizeof(int), 1, fp);
	fwrite(wo.mat, sizeof(matrixelement), wo.n * wo.m, fp);
	fclose(fp);
}
```

#### load 関数

load関数では、上記save関数で書き出した学習済みネットワークの各ノードでの重みをファイルから読み出します。
各層のサイズを読み込んだ後に、各層の重みを表す行列を読み込みます。
この関数により、前回実行した際の学習結果を再利用することができます。

なお、Load関数に関しては、手書き数字の認識とアヤメの品種推定で、別々の関数を定義しています。下記コードは手書き数字認識のコードであり、アヤメの品種推定は後述しています。

```cpp
void load(string fname) {
	FILE* fp = fopen(fname.c_str(), "rb");
	fread(&input_nodes, sizeof(int), 1, fp);
	fread(&hidden_nodes, sizeof(int), 1, fp);
	fread(&output_nodes, sizeof(int), 1, fp);
	fread(&hidden_layer, sizeof(int), 1, fp);
	for (int i = 0; i < hidden_layer; i++) {
		w.push_back(matrix());
		fread(&w[i].n, sizeof(int), 1, fp);
		fread(&w[i].m, sizeof(int), 1, fp);
		w[i].newmat(w[i].n, w[i].m);
		fread(w[i].mat, sizeof(matrixelement), w[i].n * w[i].m, fp);
	}
	fread(&wo.n, sizeof(int), 1, fp);
	fread(&wo.m, sizeof(int), 1, fp);
	wo.newmat(wo.n, wo.m);
	fread(wo.mat, sizeof(matrixelement), wo.n * wo.m, fp);
}
```

#### init 関数

予め設定されたパラメータを基に、隠れ層のランダムな値での初期化といった処理を行います。

```cpp
void init() {
	int bef = input_nodes;
	for (int k = 0; k < hidden_layer; k++) {
		w.push_back(matrix());
		w[k].newmat(hidden_nodes, bef);
		bef = hidden_nodes;
		w[k].random_m05to05();
	}
	wo.newmat(output_nodes, hidden_nodes);
	wo.random_m05to05();
}
```

#### train 関数

ニューラルネットワークの学習を行う関数です。
入力データと答えのデータ、学習率を指定して、学習を行う関数です。
出力層に対しては、ニューラルネットワークで答えを出した後、正解との差に学習率をかけ、加算していきます。
また、隠れ層に関しては、誤差逆伝搬を用いてエラーを修正していきます。

```cpp
void train(matrix* input, matrix* ans, float rate) {
	matrix final_outputs;
	vector<matrix> hidden_inputs;
	vector<matrix> hidden_outputs;
	hidden_inputs.push_back(matrix());
	hidden_inputs[0].copy(input);
	for (int i = 0; i < w.size(); i++) {
		hidden_outputs.push_back(matrix());
		hidden_outputs[i].dot(w[i], hidden_inputs[i]);
		hidden_outputs[i].sigmoid();
		hidden_inputs.push_back(matrix());
		hidden_inputs[i + 1].copy(&hidden_outputs[i]);
	}
	final_outputs.dot(wo, hidden_inputs[w.size()]);
	final_outputs.sigmoid();
    
	matrix output_errors(ans);
	output_errors.minus(final_outputs);

	output_errors.multiplication(&final_outputs);
	matrix s(ans->n, ans->m), t, u;
	s.scalar_plus(1.0f);
	s.minus(final_outputs);
	output_errors.multiplication(&s);
	t.t_set(&hidden_outputs[w.size() - 1]);
	u.dot(output_errors, t);

	final_outputs.clear();
	s.clear();

	wo.plus(u);
	t.clear();
	u.clear();

    // Backpropagation
	for (int i = w.size() - 1; i >= 0; i--) {
		matrix hidden_errors;
		matrix wt;
		wt.t_set(&w[i]);
		hidden_errors.dot(wt, output_errors);
		output_errors.copy(&hidden_errors);

		matrix as(hidden_errors.n, hidden_errors.m);
		as.scalar_plus(1.0f);
		as.minus(hidden_outputs[i]);
		matrix at(&hidden_errors);
		at.multiplication(&hidden_outputs[i]);
		at.multiplication(&as);
		matrix au;
		au.t_set(&hidden_inputs[i]);
		matrix av;
		av.dot(au, av);
		av.scalar(rate);

		w[i].plus(av);
		wt.clear();
		as.clear();
		at.clear();
		au.clear();
		av.clear();
	}

	for (int i = 0; i < w.size(); i++) {
		hidden_inputs[i].clear();
		hidden_outputs[i].clear();
	}
	hidden_inputs[w.size()].clear();
	output_errors.clear();
}
```

#### query 関数

学習済みニューラルネットワークで、入力に対する出力を得ます。

```cpp
	void query(matrix* input, matrix* output) {
		matrix x(input);
		for (int i = 0; i < w.size(); i++) {
			matrix o;
			o.dot(w[i], x);
			o.sigmoid();
			x.copy(&o);
			o.clear();
		}
		output->dot(wo, x);
		output->sigmoid();
		x.clear();
	}
};
```

ニューラルネットに関するプログラムは、ここまでです。

## データの読み込み学習や入力などを行う部分

続いて、入出力などを行う処理を書いていきました。
手順としては

1. 入力数や出力数など、ニューラルネットワークの設定
2. 学習済みネットワークの読み込み
3. テストデータによるネットワークの評価

という手順で行っています。

また、defineでモードを切り替えられるようにしています。

|define|内容|
|---|---|
|MODE_USE_TRAINED_DATA|学習済みネットワークを利用する|
|MODE_OUTPUT_ON|テストデータの分類時に、詳細を表示する|


```cpp
void loadfile(string fname, vector<matrix>* inputs, vector<matrix>* outputs) {
	// Load file
	FILE* fp = NULL;
	fp = fopen(fname.c_str(), "r");

	while (!feof(fp))
	{
		char buf[1024 * 5];
		fgets(buf, 1024 * 5, fp);
		CellsTy cells = cells_split(buf, ',');
		if (cells.size() < 785) continue;

		matrix input(784, 1), output(10, 1);
		for (int i = 0; i < 784; i++) {
			input.put(i, 0, atof(cells[i + 1].c_str()));
		}
		input.scalar_division(255.0f);
		input.scalar(0.999f);

		output.scalar_plus(ZERO);
		output.put(atoi(cells[0].c_str()), 0, ONE);
		inputs->push_back(input);
		outputs->push_back(output);
//		output.output();
	}
	fclose(fp);
}

//#define MODE_USE_TRAINED_DATA
//#define MODE_OUTPUT_ON

int main() {

	// Setting neural network
	NeuralNetwork nn;
	nn.input_nodes = 784;
	nn.hidden_nodes = 200;
	nn.hidden_layer = 1;
	nn.output_nodes = 10;
	nn.learning_rate = 0.1;
	int epochs = 10;

	vector<matrix> inputs;
	vector<matrix> outputs;

	nn.init();

    // Load or training network 
#ifdef MODE_USE_TRAINED_DATA
	nn.load("w");
#else
	loadfile("./sample_data/mnist_train_small.csv", &inputs, &outputs);
	cout << "Start training" << endl;
	for (int i = 0; i < epochs; i++) {
		cout << "epoch:	" << i << endl;
		for (int j = 0; j < inputs.size(); j++) {
			nn.train(&inputs[j], &outputs[j], 0.1);
		}
	}
	cout << "Training Done !" << endl;
	nn.save("w");
	cout << "Saved." << endl;
#endif

    // Test network
	inputs.clear();
	outputs.clear();

	loadfile("./sample_data/mnist_test.csv", &inputs, &outputs);

	int cou = 0;
	int sum = 0;
	for (int j = 0; j < inputs.size(); j++) {
		int ans = 0;
		for (int i = 1; i < 10; i++) {
			if (outputs[j].get(i, 0) > outputs[j].get(ans, 0)) ans = i;
		}

		nn.query(&inputs[j], &outputs[j]);
		int mindex = 0;
		for (int i = 1; i < 10; i++) {
			if (outputs[j].get(i, 0) > outputs[j].get(mindex, 0)) mindex = i;
		}
#ifdef MODE_OUTPUT_ON
		outputs[j].output();
		cout << "Answer:" << ans << endl;
		if (ans == mindex) cout << "success" << endl;
		else cout << "failed" << endl;
#endif
		if (ans == mindex) cou++;
		sum++;
	}
	cout << "Correct answer：" << cou << endl;
	cout << "Fialed answer：" << sum - cou << endl;
	cout << "Success rate：" << (float)cou / (float)(sum) << endl;
}
```

# アヤメの品種推定

ソースコードを若干変更することで、アヤメの品種推定もできるようになりました。
アヤメの品種推定は、精度91%程でした。

![](https://gyazo.com/6f4bc3643e272e444d55b4facfceb1ba.png)

### 前処理

[Pandasのテスト用データとしてGithubに存在するCSVファイル](https://raw.githubusercontent.com/pandas-dev/pandas/master/pandas/tests/data/iris.csv)を、ランダムに並び替えた後、学習用データ7割とテスト用データ3割に分割しました。

```py
import pandas as pd
import numpy as np

df = pd.read_csv('iris.csv')
df = df.sample(frac=1)

df = df.replace('Iris-setosa', '0')
df = df.replace('Iris-versicolor', '1')
df = df.replace('Iris-virginica', '2')

df.insert(0, 'IrisName', df['Name'])

ratio = 0.7
p = int(ratio * len(df))
df1 = df.iloc[:p, :]
df2 = df.iloc[p:, :]

df1.to_csv("iris_train.csv", index=False, header=None)
df2.to_csv("iris_test.csv", index=False, header=None)
```

### 変更部分のコード
    
```cpp
void loadfile(string fname, vector<matrix>* inputs, vector<matrix>* outputs,int input_nodesize,int output_nodesize) {
	// Load file
	FILE* fp = NULL;
	fp = fopen(fname.c_str(), "r");

	while (!feof(fp))
	{
		char buf[1024 * 5];
		fgets(buf, 1024 * 5, fp);
		CellsTy cells = cells_split(buf, ',');
		if (cells.size() < input_nodesize+1) continue;

		matrix input(input_nodesize, 1), output(output_nodesize, 1);
		for (int i = 0; i < input_nodesize; i++) {
			input.put(i, 0, atof(cells[i + 1].c_str()));
		}
//		input.scalar_division(255.0f);
//		input.scalar(0.999f);

		output.scalar_plus(ZERO);
		output.put(atoi(cells[0].c_str()), 0, ONE);
		inputs->push_back(input);
		outputs->push_back(output);
	}
	fclose(fp);
}

//#define MODE_USE_TRAINED_DATA
//#define MODE_OUTPUT_ON

int main() {

	// Setting neural network
	NeuralNetwork nn;
	nn.input_nodes = 4;
	nn.hidden_nodes = 84;
	nn.hidden_layer = 1;
	nn.output_nodes = 3;
	int epochs = 200;

	vector<matrix> inputs;
	vector<matrix> outputs;

	nn.init();

	// Load or training network 
#ifdef MODE_USE_TRAINED_DATA
	nn.load("w");
#else
	loadfile("./sample_data/ayame.csv", &inputs, &outputs, nn.input_nodes, nn.output_nodes);
	cout << "Start training" << endl;
	for (int i = 0; i < epochs; i++) {
		cout << "epoch:	" << i << endl;
		for (int j = 0; j < inputs.size(); j++) {
			nn.train(&inputs[j], &outputs[j], 0.01);
		}
	}
	cout << "Training Done !" << endl;
	nn.save("w");
	cout << "Saved." << endl;
#endif

	// Test network
	inputs.clear();
	outputs.clear();

	loadfile("./sample_data/ayame_test.csv", &inputs, &outputs, nn.input_nodes, nn.output_nodes);

	int cou = 0;
	int sum = 0;
	for (int j = 0; j < inputs.size(); j++) {
		int ans = 0;
		for (int i = 1; i < nn.output_nodes; i++) {
			if (outputs[j].get(i, 0) > outputs[j].get(ans, 0)) ans = i;
		}

		nn.query(&inputs[j], &outputs[j]);
		int mindex = 0;
		for (int i = 1; i < nn.output_nodes; i++) {
			if (outputs[j].get(i, 0) > outputs[j].get(mindex, 0)) mindex = i;
		}
#ifdef MODE_OUTPUT_ON
		outputs[j].output();
		cout << "Answer:" << ans << endl;
		if (ans == mindex) cout << "success" << endl;
		else cout << "failed" << endl;
#endif
		if (ans == mindex) cou++;
		sum++;
	}
	cout << "Correct answer：" << cou << endl;
	cout << "Fialed answer：" << sum - cou << endl;
	cout << "Success rate：" << (float)cou / (float)(sum) << endl;
}
```
## ソースコードの全体像

最初に記した通り、[GitHub](https://github.com/takoyaki-3/c-NeuralNetwork)にアップしていますので、併せて参照してください。

読み込み部分を変更することで、Neural Network部分に関しては変更を加えずとも

- 手書き数字認識
- スイセンの品種分類

が9割ほどの精度でできたので、他のデータに関しても試せればと思います。