ublas.cpp

#include <iostream>
#include <random>
#include <algorithm>
#include <fstream>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace boost::numeric;

const unsigned input_layer_size = 784;
const unsigned hidden_layer_size = 100;
const unsigned output_layer_size = 10;

const double learning_rate = 0.1;

static ublas::matrix<double> foreach(ublas::matrix<double> matrix, std::function<double(double&)> activate) {
    ublas::matrix<double> modified(matrix.size1(), matrix.size2());
    for (unsigned i = 0; i < modified.size1(); ++i) {
        for (unsigned j = 0; j < modified.size2(); ++j) {
            modified(i, j) = activate(matrix(i, j));
        }
    }
    return modified;
}

class NeuralNetwork {
    ublas::matrix<double> input_to_hidden_weights;
    ublas::matrix<double> hidden_to_output_weights;

    std::function<double(double&)> activate_neuron =
            [] (double& x) -> double {
                return 1.0 / (1.0 + exp(-x));
            };

    void init_input_to_hidden_weights() {
        std::random_device device;
        std::mt19937 mt(device());

        for (unsigned i = 0; i < input_to_hidden_weights.size1(); ++i) {
            for (unsigned j = 0; j < input_to_hidden_weights.size2(); ++j) {
                input_to_hidden_weights(i, j) = std::generate_canonical<double, 10>(mt) * pow(hidden_layer_size, -0.5);
            }
        }
    }

    void init_hidden_to_output_weights() {
        std::random_device device;
        std::mt19937 mt(device());

        for (unsigned i = 0; i < hidden_to_output_weights.size1(); ++i) {
            for (unsigned j = 0; j < hidden_to_output_weights.size2(); ++j) {
                hidden_to_output_weights(i, j) = std::generate_canonical<double, 10>(mt) * pow(output_layer_size, -0.5);
            }
        }
    }

    void init_weights() {
       init_input_to_hidden_weights();
       init_hidden_to_output_weights();
    }

public:
    NeuralNetwork() {
        input_to_hidden_weights.resize(hidden_layer_size, input_layer_size);
        hidden_to_output_weights.resize(output_layer_size, hidden_layer_size);
        init_weights();
    }

    void train(ublas::matrix<double> inputs, ublas::matrix<double> targets) {
        ublas::matrix<double> hidden_output = ublas::prod(input_to_hidden_weights, inputs);
        std::cout << hidden_output << std::endl;
        hidden_output = foreach(hidden_output, activate_neuron);

        ublas::matrix<double> final_output = ublas::prod(hidden_to_output_weights, hidden_output);
        final_output = foreach(final_output, activate_neuron);

        ublas::matrix<double> final_errors = targets - final_output;
        ublas::matrix<double> hidden_errors = ublas::prod(ublas::trans(hidden_to_output_weights), final_errors);

        ublas::matrix<double> one_minus_final_output = foreach(final_output, [] (double x) -> double { return 1.0 - x; });
        ublas::matrix<double> errors_by_final_output(final_output.size1(), final_output.size2());
        for (unsigned i = 0; i < output_layer_size; i++) {
            errors_by_final_output(i, 0) = final_output(i, 0) * final_errors(i, 0) * one_minus_final_output(i, 0);
        }
        hidden_to_output_weights += ublas::prod(errors_by_final_output, ublas::trans(hidden_output)) * learning_rate;

        ublas::matrix<double> one_minus_hidden_output = foreach(hidden_output, [] (double x) -> double { return 1.0 - x; });
        ublas::matrix<double> errors_by_hidden_output(hidden_output.size1(), hidden_output.size2());
        for (unsigned i = 0; i < hidden_layer_size; i++) {
            errors_by_hidden_output(i, 0) = hidden_output(i, 0) * hidden_errors(i, 0) * one_minus_hidden_output(i, 0);
        }
        input_to_hidden_weights += ublas::prod(errors_by_hidden_output, ublas::trans(inputs)) * learning_rate;
    }

    std::vector<double> query(ublas::matrix<double> inputs) {
        ublas::matrix<double> hidden_output = ublas::prod(input_to_hidden_weights, inputs);
        double max = input_to_hidden_weights(0, 0);
        for (unsigned i = 0; i < input_to_hidden_weights.size1(); ++i) {
            for (unsigned j = 0; j < input_to_hidden_weights.size2(); ++j) {
                if (input_to_hidden_weights(i, j) > max)
                    max = input_to_hidden_weights(i, j);
            }
        }
        hidden_output = foreach(hidden_output, activate_neuron);

        ublas::matrix<double> final_output = ublas::prod(hidden_to_output_weights, hidden_output);
        final_output = foreach(final_output, activate_neuron);

        std::vector<double> result(output_layer_size);
        for (unsigned i = 0; i < final_output.size1(); ++i) {
            result[i] = final_output(i, 0);
        }
        return result;
    }

};

std::vector<std::vector<int>> read_file(std::string file_name) {
    using namespace std;
    std::ifstream reader(file_name);
    if (!reader.is_open()) {
        std::cerr << "ERRR!";
        throw 1;
    }

    string line;
    vector<vector<int>> data;
    while (getline(reader, line)) {
        stringstream ss(line);
        vector<int> result;

        while(ss.good()) {
            string substr;
            getline(ss, substr, ',');
            result.push_back(atoi(substr.c_str()));
        }

        data.emplace_back(result);
    }
    return data;
}

int main() {
    NeuralNetwork neuralNetwork;
    //cout.precision(20);

    // LEARN
    auto train_set = read_file("train_data/mnist_train_100.csv");
    for (auto &&ex : train_set) {
        ublas::matrix<double> targets = ublas::zero_matrix<double>(10, 1);
        targets(ex[0], 0) = 1.0;

        ublas::matrix<double> inputs(784, 1);
        for (unsigned j = 1; j < ex.size(); j++) {
            inputs(j-1, 0) = ex[j] / 255.0 * 0.99 + 0.01;
        }

        neuralNetwork.train(inputs, targets);
    }

    auto test_set = read_file("train_data/mnist_test_10.csv");
    for (auto &&ex : test_set) {
        ublas::matrix<double> targets = ublas::zero_matrix<double>(10, 1);
        targets(ex[0], 0) = 1.0;

        ublas::matrix<double> inputs(784, 1);
        for (unsigned j = 1; j < ex.size(); j++) {
            inputs(j-1, 0) = ex[j] / 255.0 * 0.99 + 0.01;
        }

        neuralNetwork.train(inputs, targets);

        auto result = neuralNetwork.query(inputs);
        for (auto &&a : result) {
            std::cout << a << " ";
        }
        std::cout << std::endl;
        //cout << ex[0] << " " << distance(result.begin(), max_element(result.begin(), result.end())) << endl;
    }

    return 0;
}