nn_batch.cpp

// Copyright (c) 2013, Manuel Blum
// All rights reserved.

#include <iostream>
#include <fstream>
#include <assert.h>

#include "nn.h"

// default parameters for Rprop
const RpropParams NeuralNet::p = {0.1, 50, 1e-6, 0.5, 1.2};

NeuralNet::NeuralNet(Eigen::VectorXi &topology)
{
	assert(topology.size() > 1);
	init_layer(topology);
	init_weights(0.5);
	autoscale_reset();
}

NeuralNet::NeuralNet(const char *filename)
{
	std::ifstream fs(filename, std::ios::in | std::ios::binary);
	if (fs)
	{
		// number of layers
		int num_layers;
		fs.read((char *)&num_layers, sizeof(int));
		Eigen::VectorXi topology(num_layers);
		// topology
		fs.read((char *)topology.data(), topology.rows() * sizeof(int));
		init_layer(topology);
		autoscale_reset();
		// scaling parameters
		fs.read((char *)Xscale.data(), Xscale.size() * sizeof(F_TYPE));
		fs.read((char *)Xshift.data(), Xshift.size() * sizeof(F_TYPE));
		fs.read((char *)Yscale.data(), Yscale.size() * sizeof(F_TYPE));
		fs.read((char *)Yshift.data(), Yshift.size() * sizeof(F_TYPE));
		// weights
		for (int i = 1; i < layer.size(); ++i)
		{
			fs.read((char *)layer[i].W.data(), layer[i].W.rows() * layer[i].W.cols() * sizeof(F_TYPE));
			fs.read((char *)layer[i].b.data(), layer[i].b.size() * sizeof(F_TYPE));
		}
	}
}

void NeuralNet::init_layer(Eigen::VectorXi &topology)
{
	// init input layer
	Layer l;
	l.size = topology(0);
	layer.push_back(l);
	// init hidden and output layer
	for (int i = 1; i < topology.size(); ++i)
	{
		Layer l;
		l.size = topology(i);
		l.W.setZero(l.size, layer[i - 1].size);
		l.b.setZero(l.size);
		l.dEdW.setZero(l.size, layer[i - 1].size);
		l.dEdb.setZero(l.size);
		// set initial Delta
		l.DeltaW.resize(l.size, layer[i - 1].size);
		l.Deltab.resize(l.size);
		l.DeltaW.setConstant(p.Delta_0);
		l.Deltab.setConstant(p.Delta_0);
		layer.push_back(l);
	}
}

void NeuralNet::init_weights(F_TYPE sd)
{
	for (int i = 1; i < layer.size(); ++i)
	{
		layer[i].W.setRandom();
		layer[i].b.setRandom();
		layer[i].W *= sd;
		layer[i].b *= sd;
	}
}

NeuralNet::~NeuralNet()
{
}

F_TYPE NeuralNet::loss(const matrix_t &X, const matrix_t &Y, F_TYPE lambda)
{
	assert(layer.front().size == X.cols());
	assert(layer.back().size == Y.cols());
	assert(X.rows() == Y.rows());
	// number of samples
	size_t m = X.rows();
	// forward pass
	forward_pass(X);
	// compute error
	matrix_t error = layer.back().a - ((Y.rowwise() - Yshift.transpose()) * Yscale.asDiagonal());
	// compute cost
	F_TYPE J = 0.5 * error.rowwise().squaredNorm().mean();
	// compute delta
	layer.back().delta = (error.array() * sigmoid_gradient(layer.back().a).array()).matrix();
	for (size_t i = layer.size() - 2; i > 0; --i)
	{
		matrix_t g = sigmoid_gradient(layer[i].a);
		layer[i].delta = (layer[i + 1].delta * layer[i + 1].W).cwiseProduct(g);
	}
	// compute partial derivatives and RPROP parameters
	for (int i = 1; i < layer.size(); ++i)
	{
		// add regularization to weights, bias weights are not regularized
		J += 0.5 * lambda * layer[i].W.array().square().sum() / m;
		matrix_t dEdW = (layer[i].delta.transpose() * layer[i - 1].a + lambda * layer[i].W) / m;
		vector_t dEdb = layer[i].delta.colwise().sum().transpose() / m;
		layer[i].directionW = layer[i].dEdW.cwiseProduct(dEdW);
		layer[i].directionb = layer[i].dEdb.cwiseProduct(dEdb);
		layer[i].dEdW = dEdW;
		layer[i].dEdb = dEdb;
	}
	return J;
}
double NeuralNet::measure_accuracy(matrix_t Y, matrix_t Y_test)
{
	int correct = 0;
	for (int i = 0; i < Y.rows(); i++)
	{
		int j = 0;
		for (; j < Y.cols(); j++)
		{
			int out = Y_test(i, j) > 0.7 ? 1 : 0;
			//printf("%3i   %4.4f", (int)Y(i, j), Y_test(i, j));
			if (out != Y(i, j))
				break;
		}
		//cout << endl;
		if (j == Y.cols())
			correct++;
	}
	//cout << "Correct answer " << correct << " among " << Y.rows() << " test" << endl;
	double accuracy=(double)correct / (double)Y.rows() * 100;
	return accuracy;
	//cout << "Accuracy: " <<  accuracy << "%" << endl;
}
void NeuralNet::rprop()
{
	for (int i = 1; i < layer.size(); ++i)
	{
		//std::cout << "#######Layer " << i << "###########" << std::endl;
		for (int j = 0; j < layer[i].size; ++j)
		{
			for (int k = 0; k < layer[i - 1].size; ++k)
			{
				F_TYPE u = rprop_update(layer[i].directionW(j, k), layer[i].DeltaW(j, k), layer[i].dEdW(j, k));
				layer[i].W(j, k) += u;
			}
			layer[i].b(j) += rprop_update(layer[i].directionb(j), layer[i].Deltab(j), layer[i].dEdb(j));
		}
	}
}

F_TYPE NeuralNet::rprop_update(F_TYPE &direction, F_TYPE &Delta, F_TYPE grad)
{
	if (direction > 0)
	{
		Delta = std::min(Delta * p.eta_plus, p.Delta_max);
		direction = grad;
		if (grad > 0)
			return -Delta;
		else
			return Delta;
	}
	else if (direction < 0)
	{
		Delta = std::max(Delta * p.eta_minus, p.Delta_min);
		direction = 0;
		return 0;
	}
	else
	{
		direction = grad;
		if (grad > 0)
			return -Delta;
		else
			return Delta;
	}
}

void NeuralNet::rprop_reset()
{
	for (int i = 1; i < layer.size(); ++i)
	{
		layer[i].dEdW.setZero();
		layer[i].dEdb.setZero();
		layer[i].DeltaW.setConstant(p.Delta_0);
		layer[i].Deltab.setConstant(p.Delta_0);
	}
}

void NeuralNet::forward_pass(const matrix_t &X)
{
	assert(layer.front().size == X.cols());
	// copy and scale data matrix
	layer[0].a = (X.rowwise() - Xshift.transpose()) * Xscale.asDiagonal();
	for (int i = 1; i < layer.size(); ++i)
	{
		// compute input for current layer
		layer[i].z = layer[i - 1].a * layer[i].W.transpose();
		// add bias
		layer[i].z.rowwise() += layer[i].b.transpose();
		// apply activation function
		layer[i].a = sigmoid(layer[i].z);
	}
}

matrix_t NeuralNet::get_activation()
{
	return (layer.back().a * Yscale.asDiagonal().inverse()).rowwise() + Yshift.transpose();
}

matrix_t NeuralNet::sigmoid(const matrix_t &x)
{
	return ((-x).array().exp() + 1.0).inverse().matrix();
}

matrix_t NeuralNet::sigmoid_gradient(const matrix_t &x)
{
	return x.cwiseProduct((1.0 - x.array()).matrix());
}

void NeuralNet::gradient_descent(F_TYPE alpha)
{
	for (int i = 1; i < layer.size(); ++i)
	{
		layer[i].W -= alpha * layer[i].dEdW;
		layer[i].b -= alpha * layer[i].dEdb;
	}
}

bool NeuralNet::write(const char *filename)
{
	// open file
	std::ofstream fs(filename, std::ios::out | std::ios::binary);
	// write everything to disk
	if (fs)
	{
		// number of layers
		int num_layers = layer.size();
		fs.write((char *)&num_layers, sizeof(int));
		// topology
		for (int i = 0; i < layer.size(); ++i)
		{
			int num_units = layer[i].size;
			fs.write((char *)&num_units, sizeof(int));
		}
		// scaling parameters
		fs.write((char *)Xscale.data(), Xscale.size() * sizeof(F_TYPE));
		fs.write((char *)Xshift.data(), Xshift.size() * sizeof(F_TYPE));
		fs.write((char *)Yscale.data(), Yscale.size() * sizeof(F_TYPE));
		fs.write((char *)Yshift.data(), Yshift.size() * sizeof(F_TYPE));
		// weights
		for (int i = 1; i < layer.size(); ++i)
		{
			fs.write((char *)layer[i].W.data(), layer[i].W.rows() * layer[i].W.cols() * sizeof(F_TYPE));
			fs.write((char *)layer[i].b.data(), layer[i].b.size() * sizeof(F_TYPE));
		}
	}
	else
	{
		return false;
	}
	fs.close();
	return true;
}

void NeuralNet::autoscale(const matrix_t &X, const matrix_t &Y)
{
	assert(layer.front().size == X.cols());
	assert(layer.back().size == Y.cols());
	assert(X.rows() == Y.rows());
	// compute the mean of the input data
	Xshift = X.colwise().mean();
	// compute the standard deviation of the input data
	Xscale = (X.rowwise() - Xshift.transpose()).array().square().colwise().mean().array().sqrt().inverse();
	for (size_t i = 0; i < Xscale.size(); ++i)
		if (Xscale(i) > 10e9)
			Xscale(i) = 1;
	// compute the minimum target values
	Yshift = Y.colwise().minCoeff();
	// compute the maximum shifted target values
	Yscale = (Y.colwise().maxCoeff() - Yshift.transpose()).array().inverse();
	for (size_t i = 0; i < Yscale.size(); ++i)
		if (Yscale(i) > 10e9)
			Yscale(i) = 1;
}

void NeuralNet::autoscale_reset()
{
	Xscale = vector_t::Ones(layer.front().size);
	Xshift = vector_t::Zero(layer.front().size);
	Yscale = vector_t::Ones(layer.back().size);
	Yshift = vector_t::Zero(layer.back().size);
}