手写神经网络：C++版

最近为了对比神经网络在不同语言中的运行效率，于是手搓了一个简单模型在不同语言中的代码，在这儿记录一下。
基于MINIST手写数字数据集，采用3层全连接网络，由Pytorch、Python、C++、CUDA四个版本实现。

0. 实验结果
1. Pytorch版
2. Python版
   2.1.Numpy版
   2.2.纯Python版
3. C++版
4. CUDA版

1. 采用C++实现的全连接神经网络，采用一维数组的方式对数据进行处理
2. 为了简化代码，目前batchsize只能设置为1
3. 在C++中，高维数组通常采用一维的方式进行处理。因为二维数组做函数形参时，需要写死数组的形状，例如float a[][10]或float (*a)[10]；或者以二级指针的方式传入：float** a，但是最终还是要以一维的方式进行索引。另外，一维数组的处理方式速度会更快，可以参照这个案例：https://zhuanlan.zhihu.com/p/26896837
4. 二维数组的一维索引方式。对于一个形状为(m,n)的二维数组，在一维中的索引为a[i][j]=a[i*n+j]。
5. 如果要实现多batchsize的版本，可以将代码前向传播、反向传播部分的矩阵相乘进行修改，将三维数组以一维的方式进行处理，思路与上面二维的类似。

#pragma once
using namespace std;

void read_Mnist_Images(string filename, float* img_array);
void read_Mnist_Label(string filename, float* label_array);

//dataloader.cpp
#include 
#include 
#include 
#include 
using namespace std;

int ReverseInt(int i)
{
	unsigned char ch1, ch2, ch3, ch4;
	ch1 = i & 255;
	ch2 = (i >> 8) & 255;
	ch3 = (i >> 16) & 255;
	ch4 = (i >> 24) & 255;
	return((int)ch1 << 24) + ((int)ch2 << 16) + ((int)ch3 << 8) + ch4;
}


void read_Mnist_Label(string filename, float* label_array)
{
	ifstream file(filename, ios::binary);
	if (file.is_open())
	{
		int magic_number = 0;
		int number_of_images = 0;
		file.read((char*)&magic_number, sizeof(magic_number));
		file.read((char*)&number_of_images, sizeof(number_of_images));
		magic_number = ReverseInt(magic_number);
		number_of_images = ReverseInt(number_of_images);
		cout << "magic number = " << magic_number << endl;
		cout << "number of images = " << number_of_images << endl;
		for (int i = 0; i < number_of_images; i++)
		{
			unsigned char label = 0;
			file.read((char*)&label, sizeof(label));
			label_array[i * 10 + (int)label] = 1;//在这里将标签转换为onehot格式
		}
	}
	else {
		cout << "open file failed." << endl;
	}
	file.close();
}

void read_Mnist_Images(string filename, float* img_array)
{
	ifstream file(filename, ios::binary);
	if (file.is_open())
	{
		int magic_number = 0;
		int number_of_images = 0;
		int n_rows = 0;
		int n_cols = 0;
		file.read((char*)&magic_number, sizeof(magic_number));
		file.read((char*)&number_of_images, sizeof(number_of_images));
		file.read((char*)&n_rows, sizeof(n_rows));
		file.read((char*)&n_cols, sizeof(n_cols));
		magic_number = ReverseInt(magic_number);
		number_of_images = ReverseInt(number_of_images);
		n_rows = ReverseInt(n_rows);
		n_cols = ReverseInt(n_cols);
		cout << "magic number = " << magic_number << endl;
		cout << "number of images = " << number_of_images << endl;
		cout << "rows = " << n_rows << endl;
		cout << "cols = " << n_cols << endl;
		for (int i = 0; i < number_of_images; i++)
		{
			int img_num = 0;
			for (int r = 0; r < n_rows; r++)
			{
				for (int c = 0; c < n_cols; c++)
				{

					unsigned char image = 0;
					file.read((char*)&image, sizeof(image));
					img_array[i * 784 + img_num] = (float)image/255.0;
					img_num++;
				}
			}
		}
	}
	else {
		cout << "open file failed." << endl;
	}
	file.close();
}

#pragma once

void matrix_transpose(float* a, int shape0, int shape1);
void softmax(float* input,float* output);
float cross_entropy_error(float* y, float* t);
void FC_init(float* params, int size);
void SGD_update(float* params, float* grads, int size);

#include 
#include 
#include 
#include 
#include 
#include 
#include 

using namespace cv;
using namespace std;

void matrix_transpose(float* a,int shape0,int shape1)
{
	//对以1维数组形式表示的2维数组进行转置
	float* b = new float[shape0 * shape1];
	for (int i = 0; i < shape1; i++)
	{
		for (int j = 0; j < shape0; j++)
		{
			b[i * shape0 + j] = a[i + shape1 * j];
		}
	}
	memcpy(a, b, sizeof(float) * shape0 * shape1);
	delete[]b;
}

void softmax(float* input,float* output)
{
	const int input_size = 10;//在这里，默认锁定输入尺寸为10
	int max = input[0];
	float* input_exp = new float[input_size];
	//查找max
	for (int i = 1; i < input_size; i++)
	{
		if (input[i] > max)
		{
			max = input[i];
		}
	}
	float exp_sum = 0;
	for (int i = 0; i < input_size; i++)
	{
		input_exp[i] = exp(input[i] - max);
		exp_sum += input_exp[i];
	}
	for (int i = 0; i < input_size; i++)
	{
		output[i] = input_exp[i] / exp_sum;
	}
	delete[]input_exp;
}

float cross_entropy_error(float* y, float* t)
{
	const int input_size = 10;//同上，锁定为10
	float delta = 1e-7;
	float out;
	for (int i = 0; i < input_size; i++)
	{
		if (t[i] != 0)
		{
			out = -(t[i] * log(y[i] + delta));
		}
	}
	return out;
}

void FC_init(float* params, int size)
{
	default_random_engine e(time(0));
	uniform_real_distribution u(-2, 2);
	for (int i = 0; i < size; i++)
	{
		params[i] = 0.01*u(e);
	}
}

void SGD_update(float* params, float* grads, int size)
{
	for (int i = 0; i < size; i++)
	{
		params[i] -= 0.01 * grads[i];
	}
}

#pragma once



class FC
{
public:
	float* self_w;//2维数组
	float* self_b;
	float* self_x;
	float* self_dw;//2维数组
	float* self_db;
	float* self_output;
	float* self_dx;
private:
	int self_input_size;
	int self_output_size;
public:
	FC();
	FC(float* w, int input_size, int output_size, float* b);
	~FC();
	void fc_forward(float* x);
	void fc_backward(float* dout);
};

class Relu
{
public:
	Relu();
	Relu(int size);
	~Relu();
	void r_forward(float* x);
	void r_backward(float* dout);
private:
	int self_size;
public:
	float* self_mask;
	float* self_dx;
};

class Softmax_with_loss
{
public:
	Softmax_with_loss();
	~Softmax_with_loss();
	float swl_forward(float* x, float* t);
	void swl_backward();

private:
	int self_input_size;
public:
	float* self_y;
	float* self_t;
	float* self_dx;
	float self_loss;
};

#include 
#include "layers.h"
#include "functions.h"
using namespace std;

FC::FC() {}

FC::FC(float* w,int input_size,int output_size, float* b)//以一维数组的方式传递二维数组w
	:self_w(w),
	self_input_size(input_size),//w的维度0
	self_output_size(output_size),//w的维度1
	self_b(b)
{
	self_x = new float[self_input_size];
	self_dw = new float[self_input_size*self_output_size];
	self_db = new float[self_output_size];
	self_output = new float[self_output_size];
	self_dx = new float[self_input_size];
}

FC::~FC()
{
	//delete[] self_x;
	//delete[] self_dw;
	//delete[] self_db;
	//delete[] self_output;
	//delete[]self_dx;
}

void FC::fc_forward(float* x)
{
	memset(self_output, 0, sizeof(float) * self_output_size);
	memcpy(self_x, x, sizeof(float) * self_input_size);//即self_x = x，将输入x的值给self_x
	for (int i = 0; i < self_output_size; i++)
	{
		for (int j = 0; j < self_input_size; j++)
		{
			self_output[i] += x[j] * self_w[i + j*self_output_size];
		}
		self_output[i] += self_b[i];
	}
}

void FC::fc_backward(float* dout)
{

	float* self_w_T = new float[self_input_size * self_output_size];
	memcpy(self_w_T, self_w, sizeof(float) * self_input_size * self_output_size);
	memset(self_dx, 0, sizeof(float) * self_input_size);
	matrix_transpose(self_w_T, self_input_size, self_output_size);

	for (int i = 0; i < self_input_size; i++)
	{
		for (int j = 0; j < self_output_size; j++)
		{
			self_dx[i] += dout[j] * self_w_T[i  + j*self_input_size];
		}
	}
	delete[]self_w_T;

	for (int i = 0; i < self_input_size; i++)
	{
		for (int j = 0; j < self_output_size; j++)
		{
			self_dw[i * self_output_size + j] = self_x[i] * dout[j];
		}
	}
	memcpy(self_db, dout, sizeof(float) * self_output_size);
}

Relu::Relu(){}
Relu::Relu(int size):self_size(size)
{
	self_mask = new float[size];
	self_dx = new float[size];
}
Relu::~Relu()
{
	delete[]self_mask;
	delete[]self_dx;
}
void Relu::r_forward(float* x)
{
	for (int i = 0; i < self_size; i++)
	{
		x[i] = (x[i] > 0) ? x[i] : 0;
		self_mask[i] = (x[i] > 0) ? 1 : 0;
	}
}
void Relu::r_backward(float* dout)
{
	for (int i = 0; i < self_size; i++)
	{
		self_dx[i] = (self_mask[i] > 0) ? dout[i] : 0;
	}
}




Softmax_with_loss::Softmax_with_loss():self_input_size(10)
{
	self_y = new float[self_input_size];
	self_t = new float[self_input_size];
	self_dx = new float[self_input_size];
}

Softmax_with_loss::~Softmax_with_loss()
{
	//delete[]self_t;
	//delete[]self_y;
	//delete[]self_dx;
}

float Softmax_with_loss::swl_forward(float* x, float* t)
{
	memcpy(self_t, t, sizeof(float) * self_input_size);
	softmax(x, self_y);
	self_loss = cross_entropy_error(self_y, self_t);
	return self_loss;
}

void Softmax_with_loss::swl_backward()
{
	for (int i = 0; i < self_input_size; i++)
	{
		self_dx[i] = self_y[i] - self_t[i];
	}
}

#pragma once
#include "layers.h"

class Model
{
public:
	float* self_params_w1;
	float* self_params_b1;
	float* self_params_w2;
	float* self_params_b2;
	float* self_params_w3;
	float* self_params_b3;

	float* self_grads_w1;
	float* self_grads_b1;
	float* self_grads_w2;
	float* self_grads_b2;
	float* self_grads_w3;
	float* self_grads_b3;


	FC self_layer1;
	Relu self_layer1_act;
	FC self_layer2;
	Relu self_layer2_act;
	FC self_layer3;
	//Relu self_layer3_act;
	Softmax_with_loss self_layer4;


private:
	int* self_model_size;
public:
	Model(int* model_size);
	~Model();
	float* model_forward(float* x);
	float loss(float* x, float* t);
	float gradient(float* x, float* t);
	void SGD();
};

#include 
#include "main.h"
#include "layers.h"
#include "functions.h"
#include "dataloader.h"
#include 
#include 
using namespace cv;
using namespace std;


Model::Model(int* model_size)
	:self_model_size(model_size)
{
	self_params_w1 = new float[self_model_size[0] * self_model_size[1]];
	self_params_b1 = new float[self_model_size[1]];
	self_params_w2 = new float[self_model_size[1] * self_model_size[2]];
	self_params_b2 = new float[self_model_size[2]];
	self_params_w3 = new float[self_model_size[2] * self_model_size[3]];
	self_params_b3 = new float[self_model_size[3]];

	self_grads_w1 = new float[self_model_size[0] * self_model_size[1]];
	self_grads_b1 = new float[self_model_size[1]];
	self_grads_w2 = new float[self_model_size[1] * self_model_size[2]];
	self_grads_b2 = new float[self_model_size[2]];
	self_grads_w3 = new float[self_model_size[2] * self_model_size[3]];
	self_grads_b3 = new float[self_model_size[3]];

	//给params赋值
	FC_init(self_params_w1, self_model_size[0] * self_model_size[1]);
	FC_init(self_params_w2, self_model_size[1] * self_model_size[2]);
	FC_init(self_params_w3, self_model_size[1] * self_model_size[3]);
	memset(self_params_b1, 0, sizeof(float) * self_model_size[1]);
	memset(self_params_b2, 0, sizeof(float) * self_model_size[2]);
	memset(self_params_b3, 0, sizeof(float) * self_model_size[3]);

	//实例化layers
	self_layer1 = FC(self_params_w1, self_model_size[0], self_model_size[1], self_params_b1);
	self_layer1_act = Relu(self_model_size[1]);
	self_layer2 = FC(self_params_w2, self_model_size[1], self_model_size[2], self_params_b2);
	self_layer2_act = Relu(self_model_size[2]);
	self_layer3 = FC(self_params_w3, self_model_size[2], self_model_size[3], self_params_b3);
	//self_layer3_act = Relu(self_model_size[3]);
	self_layer4 = Softmax_with_loss();
}

Model::~Model()
{
	
}

float* Model::model_forward(float* x)
{
	self_layer1.fc_forward(x);
	self_layer1_act.r_forward(self_layer1.self_output);
	self_layer2.fc_forward(self_layer1.self_output);
	self_layer2_act.r_forward(self_layer2.self_output);
	self_layer3.fc_forward(self_layer2.self_output);
	//cout << self_layer1.self_output[1] << endl;

	return self_layer3.self_output;

}

float Model::loss(float* x, float* t)
{
	return self_layer4.swl_forward(x, t);//返回的是一个float型
}

float Model::gradient(float* x, float* t)
{
	//前向传播部分
	float* y=model_forward(x);
	float loss_value = loss(y, t);
	//cout << x[1] << endl;

	//反向传播部分
	self_layer4.swl_backward();
	//self_layer3_act.backward(self_layer4.self_dx);
	self_layer3.fc_backward(self_layer4.self_dx);
	self_layer2_act.r_backward(self_layer3.self_dx);
	self_layer2.fc_backward(self_layer2_act.self_dx);
	self_layer1_act.r_backward(self_layer2.self_dx);
	self_layer1.fc_backward(self_layer1_act.self_dx);

	//存储grads
	memcpy(self_grads_w1, self_layer1.self_dw, sizeof(float) * self_model_size[0] * self_model_size[1]);
	memcpy(self_grads_b1, self_layer1.self_db, sizeof(float) * self_model_size[1]);
	memcpy(self_grads_w2, self_layer2.self_dw, sizeof(float) * self_model_size[1] * self_model_size[2]);
	memcpy(self_grads_b2, self_layer2.self_db, sizeof(float) * self_model_size[2]);
	memcpy(self_grads_w3, self_layer3.self_dw, sizeof(float) * self_model_size[2] * self_model_size[3]);
	memcpy(self_grads_b3, self_layer3.self_db, sizeof(float) * self_model_size[3]);

	return loss_value;
}

void Model::SGD()
{
	SGD_update(self_params_w1, self_grads_w1, self_model_size[0] * self_model_size[1]);
	SGD_update(self_params_b1, self_grads_b1, self_model_size[1]);
	SGD_update(self_params_w2, self_grads_w2, self_model_size[1] * self_model_size[2]);
	SGD_update(self_params_b2, self_grads_b2, self_model_size[2]);
	SGD_update(self_params_w3, self_grads_w3, self_model_size[2] * self_model_size[3]);
	SGD_update(self_params_b3, self_grads_b3, self_model_size[3]);
}


int main()
{
	//实例化Model
	int model_size[4] = { 784,1024,1024,10 };
	Model DNN(model_size);

	//加载数据集
	const int train_size = 60000;
	const int test_size = 10000;
	const int label_size = 10;
	const int image_size = 784;
	float* train_label = new float[train_size * label_size];//用于存储label的数组
	memset(train_label, 0, train_size * label_size * sizeof(float));
	float* train_image = new float[train_size * image_size];//用于存储iamge的二维数组
	read_Mnist_Label("./dataset/train-labels.idx1-ubyte", train_label);
	read_Mnist_Images("./dataset/train-images.idx3-ubyte", train_image);

	//训练
	const int epoch = 1;
	for (int e = 0; e < epoch; e++)
	{
		float loss_sum = 0;
		for (int i = 0; i < 3000; i++)
		{
			float* train_x = new float[image_size];
			float* train_t = new float[label_size];
			copy(train_label + (i * label_size), train_label + (i + 1) * label_size, train_t);
			copy(train_image + (i * image_size), train_image + (i + 1) * image_size, train_x);

			float loss_ = DNN.gradient(train_x, train_t);
			DNN.SGD();
			loss_sum += loss_;

			delete[]train_x;
			delete[]train_t;

			if (i % 500 == 0)
			{
				float loss_avg = loss_sum / 500;
				cout << i << " " << loss_avg << endl;
				loss_sum = 0;
			}
		}
	}
	
	//读入一张图进行测试
	Mat img = imread("D:\1AA\BirenHomework\v3\dataset\testData\0.jpg",0);
	if (!img.data) { cout << "测试图片读取错误" << endl; }
	float test_array[784] = {};
	int a = 0;
	for (int i = 0; i < img.rows; i++)
	{
		for (int j = 0; j < img.cols; j++)
		{
			float value = (img.at(i, j)) / 255.0;
			test_array[a] = value;
			a += 1;
		}
	}
	float* test_y=DNN.model_forward(test_array);
	float test_y_softmax[10] = {};
	softmax(test_y, test_y_softmax);

	int y_out = max_element(test_y_softmax, test_y_softmax + 10) - test_y_softmax;
	cout << "预测值：" << y_out << endl;



	return 0;
}