diff --git a/libs/ml/CMakeLists.txt b/libs/ml/CMakeLists.txt
new file mode 100644
index 00000000..46124ad4
--- /dev/null
+++ b/libs/ml/CMakeLists.txt
@@ -0,0 +1,6 @@
+file(GLOB_RECURSE PSEMEK_ML_HEADERS RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "include/*.hpp")
+file(GLOB_RECURSE PSEMEK_ML_SOURCES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "source/*.cpp")
+
+psemek_add_library(psemek-ml ${PSEMEK_ML_HEADERS} ${PSEMEK_ML_SOURCES})
+target_include_directories(psemek-ml PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
+target_link_libraries(psemek-ml PUBLIC psemek-util psemek-geom psemek-random)
diff --git a/libs/ml/include/psemek/ml/activation.hpp b/libs/ml/include/psemek/ml/activation.hpp
new file mode 100644
index 00000000..186d645b
--- /dev/null
+++ b/libs/ml/include/psemek/ml/activation.hpp
@@ -0,0 +1,55 @@
+#pragma once
+
+#include <cmath>
+#include <exception>
+
+namespace psemek::ml
+{
+
+	// All activation functions are chosen in a way so that the derivative
+	//   can be expressed as a function of the activation function's value, i.e.
+	//   f'(x) = G(f(x)) for some G: R -> R
+	enum class activation_type
+	{
+		sigmoid,
+		tanh,
+		relu,
+	};
+
+	struct unknown_activation_type
+		: std::exception
+	{
+		char const * what() noexcept;
+	};
+
+	template <typename T>
+	T activation(T x, activation_type type)
+	{
+		switch (type) {
+		case activation_type::sigmoid:
+			return 1.f / (1.f + std::exp(-x));
+		case activation_type::tanh:
+			return 2.f / (1.f + std::exp(- 2.f * x)) - 1.f;
+		case activation_type::relu:
+			return std::max(T{0}, x);
+		default:
+			throw unknown_activation_type{};
+		}
+	}
+
+	template <typename T>
+	T activation_derivative(T value, activation_type type)
+	{
+		switch (type) {
+		case activation_type::sigmoid:
+			return value * (T{1} - value);
+		case activation_type::tanh:
+			return T{1} - value * value;
+		case activation_type::relu:
+			return value == T{0} ? T{0} : T{1};
+		default:
+			throw unknown_activation_type{};
+		}
+	}
+
+}
diff --git a/libs/ml/include/psemek/ml/neural_net.hpp b/libs/ml/include/psemek/ml/neural_net.hpp
new file mode 100644
index 00000000..24c976d6
--- /dev/null
+++ b/libs/ml/include/psemek/ml/neural_net.hpp
@@ -0,0 +1,109 @@
+#pragma once
+
+#include <psemek/ml/activation.hpp>
+#include <psemek/random/uniform.hpp>
+
+#include <vector>
+#include <exception>
+
+namespace psemek::ml
+{
+
+	struct empty_neural_net_error
+		: std::exception
+	{
+		char const * what() noexcept;
+	};
+
+	namespace detail
+	{
+
+		inline std::pair<std::vector<std::size_t>, std::vector<activation_type>> make_nn_ctor_args(std::vector<std::size_t> layer_sizes, activation_type type)
+		{
+			if (layer_sizes.empty())
+				throw empty_neural_net_error{};
+			std::vector<activation_type> activation_types(layer_sizes.size() - 1, type);
+			return {std::move(layer_sizes), std::move(activation_types)};
+		}
+
+		inline std::size_t weight_count(std::vector<std::size_t> const & layer_sizes)
+		{
+			std::size_t result = 0;
+			for (std::size_t l = 0; l + 1 < layer_sizes.size(); ++l)
+				result += (layer_sizes[l] + 1) * layer_sizes[l + 1];
+			return result;
+		}
+
+	}
+
+	template <typename T>
+	struct neural_net
+	{
+		neural_net() = default;
+		neural_net(std::vector<std::size_t> layer_sizes);
+		neural_net(std::vector<std::size_t> layer_sizes, activation_type type);
+		neural_net(std::vector<std::size_t> layer_sizes, std::vector<activation_type> activation_types);
+
+		// A non-empty neural net is basically unusable
+		bool empty() const { return layer_sizes_.empty(); }
+
+		std::size_t layer_count() const { return layer_sizes_.size(); }
+		std::size_t const * layer_sizes() const { return layer_sizes_.data(); }
+
+		std::size_t activation_type_count() const { return activation_types_.size(); }
+		activation_type const * activation_types() const { return activation_types_.data(); }
+
+		// Weights are stored in a sequential manner, in the order in which they
+		//   appear in evalution of the neural net, i.e. first come the weights
+		//   between layers 0 and 1 (including the bias) in row-major order, then
+		//   1-2 in row-major order, etc.
+		std::size_t weight_count() const { return weights_.size(); }
+		T const * weights() const { return weights_.data(); }
+		T * weights() { return weights_.data(); }
+
+	private:
+		std::vector<std::size_t> layer_sizes_;
+		std::vector<activation_type> activation_types_;
+		std::vector<T> weights_;
+
+		// proxy constructor to overcome unspecified evaluation order in
+		//   neural_net(std::move(layer_sizes), std::vector{layer_sizes.size(), activation_type::tanh})
+		neural_net(std::pair<std::vector<std::size_t>, std::vector<activation_type>> args);
+
+		void assert_nonempty() const
+		{
+			if (empty())
+				throw empty_neural_net_error{};
+		}
+	};
+
+	extern template struct neural_net<float>;
+	extern template struct neural_net<double>;
+
+	template <typename T>
+	neural_net<T>::neural_net(std::vector<std::size_t> layer_sizes)
+		: neural_net(std::move(layer_sizes), activation_type::tanh)
+	{}
+
+	template <typename T>
+	neural_net<T>::neural_net(std::vector<std::size_t> layer_sizes, activation_type type)
+		: neural_net(detail::make_nn_ctor_args(std::move(layer_sizes), type))
+	{}
+
+	template <typename T>
+	neural_net<T>::neural_net(std::vector<std::size_t> layer_sizes, std::vector<activation_type> activation_types)
+		: layer_sizes_(layer_sizes)
+		, activation_types_(activation_types)
+		, weights_(detail::weight_count(layer_sizes_))
+	{
+		if (layer_sizes_.empty())
+			throw empty_neural_net_error{};
+	}
+
+	template <typename T>
+	neural_net<T>::neural_net(std::pair<std::vector<std::size_t>, std::vector<activation_type>> args)
+		: neural_net(std::move(args.first), std::move(args.second))
+	{}
+
+
+}
diff --git a/libs/ml/include/psemek/ml/neural_net_evaluator.hpp b/libs/ml/include/psemek/ml/neural_net_evaluator.hpp
new file mode 100644
index 00000000..9f74ecb3
--- /dev/null
+++ b/libs/ml/include/psemek/ml/neural_net_evaluator.hpp
@@ -0,0 +1,57 @@
+#pragma once
+
+#include <psemek/ml/neural_net.hpp>
+
+#include <stdexcept>
+
+namespace psemek::ml
+{
+
+	struct wrong_neural_net_input_size
+		: std::runtime_error
+	{
+		wrong_neural_net_input_size(std::size_t expected, std::size_t actual);
+	};
+
+	// A class that stores temporary data to facilitate
+	//   allocation-free multiple evaluations of a neural net
+	template <typename T>
+	struct neural_net_evaluator
+	{
+		std::vector<T> evaluate(neural_net<T> const & nn, std::vector<T> input) const;
+	private:
+		mutable std::vector<T> temp_;
+	};
+
+	extern template struct neural_net_evaluator<float>;
+	extern template struct neural_net_evaluator<double>;
+
+	template <typename T>
+	std::vector<T> neural_net_evaluator<T>::evaluate(neural_net<T> const & nn, std::vector<T> input) const
+	{
+		std::size_t const * layer_sizes = nn.layer_sizes();
+		T const * weights = nn.weights();
+
+		if (layer_sizes[0] != input.size())
+			throw wrong_neural_net_input_size{layer_sizes[0], input.size()};
+
+		for (std::size_t l = 0; l + 1 < nn.layer_count(); ++l)
+		{
+			temp_.resize(layer_sizes[l + 1]);
+
+			for (std::size_t i = 0; i < layer_sizes[l + 1]; ++i)
+			{
+				temp_[i] = *weights++;
+
+				for (std::size_t j = 0; j < layer_sizes[l]; ++j)
+					temp_[i] += (*weights++) * input[j];
+
+				temp_[i] = activation(temp_[i], nn.activation_types()[l]);
+			}
+
+			std::swap(temp_, input);
+		}
+
+		return input;
+	}
+}
diff --git a/libs/ml/source/activation.cpp b/libs/ml/source/activation.cpp
new file mode 100644
index 00000000..76078957
--- /dev/null
+++ b/libs/ml/source/activation.cpp
@@ -0,0 +1,11 @@
+#include <psemek/ml/activation.hpp>
+
+namespace psemek::ml
+{
+
+	char const * unknown_activation_type::what() noexcept
+	{
+		return "unknown activation type";
+	}
+
+}
diff --git a/libs/ml/source/neural_net.cpp b/libs/ml/source/neural_net.cpp
new file mode 100644
index 00000000..94ce67a9
--- /dev/null
+++ b/libs/ml/source/neural_net.cpp
@@ -0,0 +1,14 @@
+#include <psemek/ml/neural_net.hpp>
+
+namespace psemek::ml
+{
+
+	char const * empty_neural_net_error::what() noexcept
+	{
+		return "neural net must have at least a single layer";
+	}
+
+	template struct neural_net<float>;
+	template struct neural_net<double>;
+
+}
diff --git a/libs/ml/source/neural_net_evaluator.cpp b/libs/ml/source/neural_net_evaluator.cpp
new file mode 100644
index 00000000..7e38e788
--- /dev/null
+++ b/libs/ml/source/neural_net_evaluator.cpp
@@ -0,0 +1,15 @@
+#include <psemek/ml/neural_net_evaluator.hpp>
+
+#include <psemek/util/to_string.hpp>
+
+namespace psemek::ml
+{
+
+	wrong_neural_net_input_size::wrong_neural_net_input_size(std::size_t expected, std::size_t actual)
+		: std::runtime_error(util::to_string("wrong neural net input size: expected ", expected, ", got ", actual))
+	{}
+
+	template struct neural_net_evaluator<float>;
+	template struct neural_net_evaluator<double>;
+
+}